
% vliv vlnovych aberaci na zobrazeni

function sfa5_1


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=1*(x.^2+y.^2)/(D/2)^2;
W=1*(x.^2)/(D/2)^2;    % astigmatismus (x)
%W=1*(x.^2+y.^2).*x/(D/2)^3;  %koma

%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);
U=lam*z/(N*del);

% obrazove souradnice

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

% 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)

[fx,fy]=meshgrid((-N/2:N/2-1)*del/lam/z);
H=ft2(abs(h).^2,del);



% prenos v ose x a y

figure(4)
% H(fy=0)
idx=(fy==0);
Hx=H(idx);
osafx=fx(idx);
Hx=Hx/Hx(osafx==0);
% H(fx=0)
idx=(fx==0);
Hy=H(idx);
osafy=fy(idx);
Hy=Hy/Hy(osafy==0);
Hy
plot(osafx,real(Hx),"-","LineWidth",2,osafy,real(Hy),"--","LineWidth",2);
xlim([0,100000]);ylim([-0.2,1.1]);
grid("on");








max(max(u))
1/lam/z
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


