
% vliv apodizace na zobrazeni

function sfa5_2


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,r]=cart2pol(x,y);

%pupilove fce

P=circ(x,y,D);
Pinv=P.*r.^4; Pinv=Pinv/max(max(Pinv)); % inverzni apodizace
Papod=P.*cos(pi/2*r/(D/2)); %apodizace


%amplitudova rozptylova fce

h=ft2(P,del);
hinv=ft2(Pinv,del);
hapod=ft2(Papod,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);
obj(100:120,50:51)=1;obj(100:120,56:57)=1;
figure(1)
pcolor(u,v,abs(obj).^2);
colormap(gray());shading("flat");axis("square");



% PSF 


figure(2)
subplot(2,3,1)
pcolor(u,v,abs(h));
colormap(gray());shading("flat");axis("square");
xlim([-0.0001,0.0001]);ylim([-0.0001,0.0001]);
title("plna apertura")
subplot(2,3,2)
pcolor(u,v,abs(hapod));
colormap(gray());shading("flat");axis("square");
xlim([-0.0001,0.0001]);ylim([-0.0001,0.0001]);
title("apodizace")
subplot(2,3,3)
pcolor(u,v,abs(hinv));
colormap(gray());shading("flat");axis("square");
xlim([-0.0001,0.0001]);ylim([-0.0001,0.0001]);
title("inverzni apodizace")




% obraz 

img=ift2(ft2(abs(obj).^2,1).*ft2(abs(h).^2,1),1);
imgapod=ift2(ft2(abs(obj).^2,1).*ft2(abs(hapod).^2,1),1);
imginv=ift2(ft2(abs(obj).^2,1).*ft2(abs(hinv).^2,1),1);

subplot(2,3,4)
pcolor(u,v,real(img));
colormap(gray());shading("flat");axis("square");
subplot(2,3,5)
pcolor(u,v,real(imgapod));
colormap(gray());shading("flat");axis("square");
subplot(2,3,6)
pcolor(u,v,real(imginv));
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);
Hapod=ft2(abs(hapod).^2,del);
Hinv=ft2(abs(hinv).^2,del);

% prenos v ose x

figure(3)
idx=(fy==0);
osafx=fx(idx);
Hx=H(idx); Hx=Hx/Hx(osafx==0);
Hxapod=Hapod(idx); Hxapod=Hxapod/Hxapod(osafx==0);
Hxinv=Hinv(idx); Hxinv=Hxinv/Hxinv(osafx==0);

plot(osafx,real(Hx),"-","LineWidth",2,osafx,real(Hxapod),"--","LineWidth",2,
osafx,real(Hxinv),"-.","LineWidth",2);
xlim([0,90000]);ylim([-0.1,1.1]);
grid("on");
legend("plna","apodizace","inverzni");








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


