%
% Optical Imaging and Spectroscopy
%
% David J. Brady
% Duke University
% www.opticalimaging.org
%
%
%
% PSF of the Log Asphere lens
%
%
%
%
close all;
figure(1);set(gcf,'color','white');
X=500; %microns
lambda=1; %micron
N=2048; % samples per axis
xrange=linspace(-X,X,N);
urange=(-N/(4*X)):(1/(2*X)):(N/2-1)/(2*X);
[x,y]=meshgrid(xrange);
[u,v]=meshgrid(urange);
%
%input pattern is the focal modulation of a lens
%
F=4*X;
f=exp(pi*i*(x.^2+y.^2)/(F*lambda)).*P(sqrt(x.^2+y.^2.5)/X);
d=F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ftm=max(max(abs(ft)));
ft=ft/ftm;
subplot(2,1,1);plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+2.5,'-k'); 
axis([ -15 15 0 7]);axis square; hold on;title('(a)');xlabel('\lambda');
d=.9875*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2,'-k');
d=.99*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+.5,'-k');
d=.9925*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+1,'-k');
d=.995*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+1.5,'-k');
d=.9975*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+2,'-k');
d=1.0025*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+3,'-k');
d=1.005*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+3.5,'-k');
d=1.0075*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+4,'-k');
d=1.01*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+4.5,'-k');
d=1.0125*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+5,'-k');
d=1.015*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+5.5,'-k');
d=1.0175*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+6,'-k');
%
%
% now a logarithmic lens
%
%
alpha=0.215/(F*lambda);gamma=1;dF=100*F;
f=exp(i*alpha*(x.^2+y.^2+F^2).*log(gamma*((x.^2+y.^2)/X^2+1))).*exp(pi*i*(x.^2+y.^2)/((F+dF)*lambda)).*P(sqrt(x.^2+y.^2)/X);
d=F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ftm=max(max(abs(ft)));
ft=ft/ftm;
subplot(2,1,2);plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+2.5,'-k'); 
axis([ -15 15 0 7]);axis square; hold on; title('(b)');xlabel('\lambda');
d=.9875*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2,'-k');
d=.99*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+.5,'-k');
d=.9925*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+1,'-k');
d=.995*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+1.5,'-k');
d=.9975*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+2,'-k');
d=1.0025*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+3,'-k');
d=1.005*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+3.5,'-k');
d=1.0075*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+4,'-k');
d=1.01*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+4.5,'-k');
d=1.0125*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+5,'-k');
d=1.015*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+5.5,'-k');
d=1.0175*F; 
fresnelT=exp(i*pi*d*lambda*(u.^2+v.^2));
ft=ifft2(fftshift(fresnelT).*fft2(f));
ft=ft/ftm;
plot(xrange((N/2-30):(N/2+30)),abs(ft(N/2+1,(N/2-30):(N/2+30))).^2+6,'-k');