% This example calculates the through-focus intensity point-spread function for a
% lambda=0.20 micron, NA=0.50, optical system. In addition the finite pinhole size
% of 0.25 micron is taken into account. It is assumed that the optical
% system suffers from astigmatism aberration (m=2; n=2).
%
% The PSF for three focus values are shown.
%
% All setting are easily changed by modifing the parameters in this routine.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m=2;
n=2;

Z=0.1;

na=0.50;        % the numerical aperture
diam=0.25;      % diameter in microns
lambda=0.200;   % wavelength in microns
eps=0.00001;    % a small offset epsilon, to avoid deviding by zero

radius=diam/2;
ap=2*pi*(na/lambda)*radius;
d=1/8*ap^2+1/384*ap^4+1/10240*ap^6; % optimal d.

scale_z=lambda/(2*pi)*1/(1-sqrt(1-na^2)); % scaling factor for defocus

r_list=[0:0.01:0.5];f_list=[-1:1:1];
[r,f]=meshgrid((r_list+eps)*2*pi*na/lambda,f_list/scale_z);

f=f+i*d;        % by using a complex defocus parameter, the finite diameter is taken into account.

[x,y]=meshgrid([-r_list(end:-1:2) r_list],[-r_list(end:-1:2) r_list]);
phi=angle(x+i*y);

u_0=zeros(size(x));
u_aber=zeros(size(x));
u_0=Vnm(r,f,0,0);
u_aber=Vnm(r,f,n,m);

for foc=1:length(f_list)
    u_0xy=interp1(r_list,u_0(foc,:),sqrt(x(:).^2+y(:).^2));
    u_0xy=reshape(u_0xy,size(x));
    u_aberxy=interp1(r_list,u_aber(foc,:),sqrt(x(:).^2+y(:).^2));
    u_aberxy=i^(m+1)*2*pi*Z*reshape(u_aberxy,size(x)).*cos(m*phi);
    field=u_0xy+u_aberxy;
    intensity=real(field.*conj(field));
    subplot(1,3,foc)
    pcolor(x,y,intensity);shading interp;axis square;axis off;
    title(['Focus = ' num2str(f_list(foc)) ' [um]'])
end