CODDEMO demonstrates spherical aberrations using coddington.m

For a given object and image distance, compute the lateral aberration
and compare to the diffraction limit.  Examples are given with different
wavelengths and different indices of refraction.

requires CODDINGTON

by Chuck DiMarzio
   Northeastern University
   April 2009
!! This file may be copied, used, or modified for educational and
!! research purposes provided that this header information is not
!! removed or altered, and provided that the book is cited in
!! publications, as DiMarzio, Charles A., Optics for Engineers,
!! CRC Press, Boca Raton, FL, 2011.
!! http://www.crcpress.com
!! Other distribution is prohibited without permission.

Contents

Set up the parameters of the problem

s=1;sp=0.04;  % Use meters for length units
f=1/(1/s+1/sp); % focal length
h=0.004; % ray height at edge of lens
p=(sp-s)/(sp+s);
[q,n]=meshgrid(*-2:0.1:7),[1.5,2.4,4]);
[ls,qmin]=coddington(p,q,f,h,n);
ds=sp^2*ls; % displacement of focus of edge ray
dx=ds*h/sp;  % height of edge ray at paraxial focus

Compute the diffraction limit (see chapter 8)

lambda=[500e-9,1.064e-6,10.59e-6]';  % Three wavelengths
dl=0.61*lambda*sp/h;   % Respective diffraction limits

fig1=figure;semilogy(q',dx'*1e6,'-',...
                     q(:,[1,end])',repmat(dl,1,2)'*1e6,'--');
axis([-2,10,1,10000]);
grid on;
xlabel('q, Shape Factor');
ylabel('\Delta x, Transverse Aberration, \mu m');
% Label the plots for different indices of refraction
for ii=1:size(n,1);
text(7.1,dx(ii,end)*1e6,sprintf('n = %5.1f',n(ii,end)));
end;
% label the diffraction limits for different wavlengths
for ii=1:length(lambda);
text(6.75,dl(ii)*1e6,[sprintf('%0.2f',lambda(ii)*1e6),' \mu m']);
end;
text(-1.5,5000,['s=',sprintf('%0.1f',s*100),' cm, ',...
             's''=',sprintf('%0.1f',sp*100),' cm, ',...
             'x_1=',sprintf('%0.1f',h*100),' cm.']);

Find the best q for given p

[ptest,ntest]=meshgrid([-1.1,1.1],n(1:3,1));
ftest=f;htest=h;
[dummy,qtest]=coddington(ptest,[],ftest,htest,ntest);
fig2=figure;plot(ptest',qtest','',[1,0,-1],[-1,0,1],'k+');grid on;
xlabel('p, Position Factor');
ylabel('q, Shape Factor');
legend(['n=',sprintf('%0.1f',n(1,1))],...
        ['    ',sprintf('%0.1f',n(2,1))],...
        ['    ',sprintf('%0.1f',n(3,1))]);
%
% Draw lens shapes along q axis
%
hold on;
pplot=-1.15;lensrad=0.5;qplotlist=[-4:2:4];
for qplot=qplotlist;
  r1=10/(qplot+1);r2=10/(qplot-1);
  [x,y,edge]=lensshape(r1,r2,lensrad,0.05);
  plot(pplot+x,qplot+y,'k-');
end;
pplot=-1.35;lensrad=0.5;qplotlist=(-3:2:5);
for qplot=qplotlist;
  r1=10/(qplot+1);r2=10/(qplot-1);
  [x,y,edge]=lensshape(r1,r2,lensrad,0.05);
  plot(pplot+x,qplot+y,'k-');
end;
%
% Draw imaging along p axis
%
w=0.2;
qplot=-5.25;
plot(-1.3+[0,0],qplot+1.4*[-w,w],'k-',...
     -1.3+[-w,w]/2,qplot+[0,0],'k-',...
     -1.3+[-w,0,w]/2,qplot+[1.1*w,w,0],'k-');
plot(-1+[0,0],qplot+1.4*[-w,w],'k-',...
     -1+[-w,w]/2,qplot+[0,0],'k-',...
     -1+[-w,0,w]/2,qplot+[w,w,0],'k-');
plot(-0.5+[0,0],qplot+1.4*[-w,w],'k-',...
     -0.5+[-3*w,w,0,-3*w]/2,qplot+[0,0,w,0],'k-');
plot(0+[0,0],qplot+1.4*[-w,w],'k-',...
     0+[-w,w,0,-w],qplot+[0,0,w,0],'k-');
plot(0.5+[0,0],qplot+1.4*[-w,w],'k-',...
     0.5+[-w,3*w,0,-w]/2,qplot+[0,0,w,0],'k-');
plot(1+[0,0],qplot+1.4*[-w,w],'k-',...
     1+[-w,w]/2,qplot+[0,0],'k-',...
     1+[-w,0,w]/2,qplot+[0,w,w],'k-');
plot(1.3+[0,0],qplot+1.4*[-w,w],'k-',...
     1.3+[-w,w]/2,qplot+[0,0],'k-',...
     1.3+[-w,0,w]/2,qplot+[0,w,1.1*w],'k-');
hold off;

Show how aberrations increase with lens diameter (2h)

p2=1;
n2=2.4;
f2=f;
lambda2=632.8e-9;
h2=(0.1:.01:1)*f/2; % Go up to f/1 lens (2h=f);
[ls2,q2]=coddington(p2,[],f2,h2,n2);
ds2=f2^2*ls2; % displacement of focus of edge ray
dx2=ds2.*h2/f2;  % height of edge ray at paraxial focus
fig3=figure;loglog(h2*1e3,dx2*1e6,'-',h2*1e3,0.61*f2*lambda2./h2*1e6,'--');
grid on;
xlabel('x_1,  mm');
ylabel('\Delta x, \mu m');
%legend('\Delta x','D.L.');
moose=axis;
% $$$ text(1.2,4000,...
% $$$      ['n=',sprintf('%0.1f',n2),', ',...
% $$$       'p=',sprintf('%0.1f',p2),', ',...
% $$$       'f=',sprintf('%0.1f',f2*100),' cm, ',...
% $$$       '\lambda=',sprintf('%0.1f',lambda2*1e9),' nm']);
% $$$ text(1.2,2000,['Optimal q = ',sprintf('%0.2f',q2)]);