function sfa4

graphics_toolkit("gnuplot")
close all


N=256;
L=0.1;
D=0.02;
del=L/N;
lam=1e-6; % vlnova delka
z=0.25; % obrazova vzdalenost

% pupilove souradnice
[x,y]=meshgrid((-N/2:N/2-1)*del);
[theta,t]=cart2pol(x,y);

% aberace vlnoplochy
W=0.5*(x.^2+y.^2)/(D/2)^2;

%pupilova fce
P=circ(x,y,D);
%P=P.*exp(I*2*pi*W); % zakomentovat pro vliv pouze pupily
%amplitudova rozptylova fce

h=ft2(P,del);
delU=lam*z/(N*del);

% obrazove souradnice

[u,v]=meshgrid((-N/2:N/2-1)*delU);

% geometricky obraz 

obj=(rect((u-1e-4)/5e-5)+rect(u/5e-5)+rect((u+1e-4)/5e-5)).*rect(v/2e-4);

figure(1)
pcolor(u,v,abs(obj).^2);
colormap(gray());shading("flat");axis("square");



% PSF 

figure(2)
pcolor(u,v,abs(h));
colormap(gray());shading("flat");axis("square");


% obraz 

img=ift2(ft2(abs(obj).^2,1).*ft2(abs(h).^2,1),1);
figure(3)
pcolor(u,v,real(img));
colormap(gray());shading("flat");axis("square");


% prenosova fce (pupila)

[fu,fv]=meshgrid((-N/2:N/2-1)*del/lam/z);
figure(4)
H=ft2(abs(h).^2,delU);
pcolor(fu,fv,real(H))
colormap(gray());shading("flat");axis("square");
xlim([-250000,250000]);ylim([-250000,250000]);


% prenos v ose x

figure(5)
idx=(fv==0);
Hu=H(idx);
fu=fu(idx);
Hu=Hu/Hu(fu==0);
plot(fu,real(Hu),"LineWidth",2);
xlim([0,100000]);ylim([-0.1,1.1]);
grid("on");




end





function G=ft(g,del)
 G=fftshift(fft(fftshift(g)))*del;
end


function g=ift(G,delf)
 g=ifftshift(ifft(ifftshift(G)))*delf*length(G);
end

function G=ft2(g,del)
 G=fftshift(fft2(fftshift(g)))*del*del;
end


function g=ift2(G,delf)
 g=ifftshift(ifft2(ifftshift(G)))*(delf*size(G,1))^2;
end


function z=circ(x,y,D)
 r=sqrt(x.^2+y.^2);
 z=double(r<D/2);
 z(r==D/2)=0.5;
end

function y=rect(x,D)
  if (nargin==1) D=1; end
  x=abs(x);
  y=double(x<D/2);
  y(x==D/2)=0.5;
end


