

% vyuziti Fraunhoferovy difrakce
% zobrazeni optikou iPhone4
% zkoumani vlivu apertury na opticke zobrazeni
% originálni clonove cislo je 2.8, vetsim zaclonenim lze zkoumat vliv difrakce na kvalitu zobrazeni 
% nekonecne vzdaleneho bodoveho predmetu

% autor Bohumil Stoklasa (mirne upraveno)

function difrakce_iPhone

f=3.85; %ohniskova vzdalenost
A=f/2.8; %apertura
px=0.00175; %velikost pixelu na kamere
lambda=0.000555; %vlnova delka
N=513; %pocet deleni vstupni apertury

%vypocet velikosti vstupniho pole pro zajisteni 
%spravneho vzorkovani ve vystupni rovine=velikost pixelu kamery
L=lambda*f/(px); 
delta=L/N; %vzorkovani ve vstupni rovine

[x1 y1]=meshgrid((-N/2:N/2-1)*delta);
%definice vstupniho pole jako kruhova apertura
Uin=circ(x1,y1,A);
%vypocet pole v rovine kamery (ohniskova rovina optiky)
%Fraunhoferova difrakce
[Uout x2 y2]=lens_against_ft(Uin,lambda,delta,f);

%intenzita ve vystupni rovine +-5 pixelu kolem stredu kamery
Iout=abs(Uout(round(N/2)-5:round(N/2)+5,round(N/2)-5:round(N/2)+5)).^2;
figure,imshow(Iout/(max(max(Iout))))


endfunction


%vypocet pole v ohniskove rovine cocky pro vstupni pole tesne pred cockou
function [Uout x2 y2]=lens_against_ft(Uin,lam,delta,f)
	N=size(Uin,1);
	k=2*pi/lam; %vlnovy vektor
	fx=(-N/2:1:N/2-1)/(N*delta); %samplovani ve frekvencni domene ohniskove roviny
	[x2 y2]=meshgrid(lam*f*fx); %samplovani v prostorove domene ohniskove roviny
	clear('fx');
	%Franhoferova difrakce
	Uout=exp(i*k/(2*f)*(x2.^2+y2.^2))/(i*lam*f).*ft2(Uin,delta);
endfunction


%kruhova apertura o velikosti D na matici o souradnicích x,y
function z=circ(x,y,D)
r=sqrt(x.^2+y.^2);
z=double(r<D/2);
z(r==D/2)=0.5;
endfunction

%2D FT
function G=ft2(g,delta)
G=fftshift(fft2(ifftshift(g)))*delta^2;
endfunction

