Contents

ASTIG demonstrates field curvature and astigmatism

requires RAYSPHERE

See Also:  ABDEMO, CHROMATIC, COMA, DISTORTION, SPHERICAL

by Chuck DiMarzio Northeastern University 2011

!! 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.

Generate pupil space in polar coordinates in rho phi

philist=(0:360)*pi/180; % radians
rholist=(0.05:0.05:1);
%

Set up parameters for the problem

zobject=-10;
xobject=2.5;
daperture=0.5; % Aperture diameter in length units
endplanes=(8:0.3:10.1); % endplanes to plot
p0=[xobject,0,zobject]'; % object location on the z axis
r=2; % radius of sphere
pc=[0,0,r]'; % sphere center
n=1.5; % index of refraction of glass
zimage=n/((n-1)/r-1/(-zobject));
%
zplot=[];xplot=[]; % arrays for plot variables
%
hmax=1;  % Maximum height of ray above lens vertex

Trace rays exactly

count=0;
phiplot=zeros(1,length(rholist)*length(philist));
rhoplot=zeros(1,length(rholist)*length(philist));
xshot=zeros(length(endplanes),length(rholist)*length(philist));
yshot=zeros(length(endplanes),length(rholist)*length(philist));
for rho=rholist; % Loop on rays: h is height above vertex
  for phi=philist;
    count=count+1;
    % point in vertex plane
    pt=[rho*cos(phi)*daperture/2,rho*sin(phi)*daperture/2,0]';
    v=pt-p0; % vector from object to lens
    v=v/sqrt(v'*v); % normalize it
    [s1,s2]=raysphere(p0,pc,v,r);
    p1=p0+s1*v;
    nor=-(p1-pc)/sqrt((p1-pc)'*(p1-pc)); % Unit normal
    vnew=v/n+(sqrt(1-(1/n)^2*(1-(v'*nor)^2))-v'*nor/n)*nor;
    p2=p1+vnew/vnew(3)*(zimage-p1(3));
    phiplot(count)=phi;
    rhoplot(count)=rho;
    for nend=1:length(endplanes);
      endplane=endplanes(nend);
      p3=p1+vnew/vnew(3)*(endplane-p1(3));
      xshot(nend,count)=p3(1);
      yshot(nend,count)=p3(2);
    end;
  end;
end;

% Don't modify figcount without changing Ch5_printfig.m
figcount=0;
for nend=1:length(endplanes);
  % Save the size for plotting vs z and the center for shot patterns
  ysize(nend)=max(yshot(nend,:))-min(yshot(nend,:));
  ycenter(nend)=(max(yshot(nend,:))+min(yshot(nend,:)))/2;
  xsize(nend)=max(xshot(nend,:))-min(xshot(nend,:));
  xcenter(nend)=(max(xshot(nend,:))+min(xshot(nend,:)))/2;
  ymin=ycenter(nend)-max(ysize)/2;
  ymax=ycenter(nend)+max(ysize)/2;
  xmin=xcenter(nend)-max(xsize)/2;
  xmax=xcenter(nend)+max(xsize)/2;
  % Don't modify figcount without changing Ch5_printfig.m
  figcount=figcount+1;
  fig(figcount)=figure;plot(yshot(nend,:),xshot(nend,:),'.');
  axis([ymin,ymax,xmin,xmax]);axis equal;
  xlabel('y'),ylabel('x');
%  title(['z = ',sprintf('%0.2f',endplanes(nend))]);
end;

% Don't modify figcount without changing Ch5_printfig.m
figcount=figcount+1;
fig(figcount)=figure;plot(endplanes,xsize,'-x',endplanes,ysize,'-o');