%
% twolens.m  
% Various analyses of a two-element lens
% by Chuck DiMarzio
%    Northeastern University
%    September 2008    
%
%    The front vertex of the first lens defines the zero of z
%
echo on
f1=10;  % Distance in centimeters
f2=5;
z12=7;
s=8;
zS=-s; % Object to the left of the first vertex
x=1;    % Object height
%
% variable with first letter z and second letter capitalized
% is the z location of the corresponding point or plane
% eg:  zF1 is the location of the front focal point of lens 1
%      zF1prime is the location of the back focal point of same
%      V denotes vertex
%      H principal plane
%      S object
%
zF1=-f1  
zF1prime=f1
zF2=z12-f2
zF2prime=z12+f2
zV1=0
zV1prime=zV1
zV2=z12
zV2prime=zV1
%
% First, calculate the image position using thin lenses
%
zS1=zS                    % Location of object on z axis
s1=-zS
s1prime=1/(1/f1-1/s1);s1prime
m1=-s1prime/s1
x1=x
x1prime=x1*m1
x2=x1prime
zS1prime=s1prime
s2=z12-s1prime
s2prime=1/(1/f2-1/s2);s2prime
m2=-s2prime/s2
zS2prime=z12+s2prime      %  Location of image on z axis
x2prime=x2*m2
m=m1*m2                   %  Total magnification
%
%  Plot object and image for two thin lenses in sequence.
%
echo off;
allz=[zS,zS1,zS1prime,zS2prime,zF1,zF1prime,zF2,zF2prime];
extender=0.05;
zlo=min(allz)*(1+extender)-max(allz)*extender
zhi=max(allz)*(1+extender)-min(allz)*extender
figure;plot([zlo,zhi],[0,0],'k-',...
            zV1*[1,1],0.5*x*[-1,1],'k-',...
            zV1prime*[1,1],0.5*x1*[-1,1],'k-',...
            zV2*[1,1],0.5*x*[-1,1],'k-',...
            zV2prime*[1,1],0.5*x*[-1,1],'k-',...
            zS1*[1,1],[0,x1],'m-',...
            zS1prime*[1,1],[0,x1prime],'b--',...
            zS2prime*[1,1],[0,x2prime],'b-',...
            zF1,0,'ko',...
            zF1prime,0,'ko',...
            zF2,0,'ko',...
            zF2prime,0,'ko');
text(zF1,-x*0.4,'F1');
text(zF1prime,-x*0.4,'F1''');
text(zF2,x*0.2,'F2');
text(zF2prime,x*0.2,'F2''');
moose=text(zS1,x*1.2,'S1');set(moose,'color','magenta');
moose=text(zS1prime,x1prime+x*0.2,'S1''');set(moose,'color','blue');
moose=text(zS2prime,x2prime*1.2,'S2''');set(moose,'color','blue');
disp(['Plotting Figure ',num2str(gcf)]);
print -depsc 2lens1.eps
echo on;
%
%  Now solve with matrices
%
L_1=[1,0;-1/f1,1]         % Matrix from H1 to H1'
T_12=[1,z12;0,1]          % Matrix from H1' to H2
L_2=[1,0;-1/f2,1]          % Matrix from H2 to H2'
%
%  The statement that lens matrices Lx are from Hx to Hx' don't matter yet,
%  because Hx, Vx, Hx' and Vx' are all in the same place for a thin lens. 
%
M_H1H2prime=L_2*T_12*L_1  % Computed Matrix from H1 to H2'
h=(1-M_H1H2prime(2,2))/M_H1H2prime(2,1) % Distance by which front principal
                                        % plane is in front of H1
hprime=...
  (1-M_H1H2prime(1,1))/M_H1H2prime(2,1) % Distance by which back principal
                                        % plane is in front of H1
f=-1/M_H1H2prime(2,1)                   % Focal length
zH=-h                    % Front principal plane location on z axis
zHprime=z12+hprime       %  Back principal plane location on z axis
zF=zH-f
zFprime=zHprime+f
%
%  Now use thick lenses.  Lenses are assumed symmetric biconvex
%  with thickness of 6 millimeters.  Thus the distance between principal
%  planes is 2 mm, and the principal planes are 2 mm in from their
%  respective vertices.
%
%  Change the vertex spacing so that the focal length of the combination
%  stays the same.
%
t=0.6;
zV1=0;   % Origin of coordinates at front vertex of first lens
%
% Locate all vertices, principal planes, and focal points
%
zH1=t/3
zF1=zH1-f1
zH1prime=2*t/3
zV1prime=t   %  back vertex of first lens
zF1prime=zH1prime+f1
%  Keep s the same as before, so move the object.
zS=zH1-s
zS1=zS;
zH2=zH1prime+z12
zF2=zH2-f2
zH2prime=zH2+t/3
zV2=zH2-t/3  %  front vertex of second lens
zV2prime=zV2+t
zF2prime=zH2prime+f2
%
% Solve the thick lens problem
s1=zH1-zS1
s1prime=1/(1/f1-1/s1)
s2=(zH2-zH1prime)-s1prime
s2prime=1/(1/f2-1/s2)
zS2prime=zH2prime+s2prime
%
%  Remember the matrix stays the same as before, but now we interpret it as
%  going from H1 to H2' so we need to recompute the physical locations.
%
zH=zH1-h   %  Remember h was computed from the matrix M_H1H2prime
zHprime=zH2prime+hprime
zF=zH-f
zFprime=zHprime+f
%
%  Plot object and image for two thick lenses in sequence.
%
echo off;
allz=[zS,zS1,zS2prime,zF1,zF1prime,zF2,zF2prime];
extender=0.05;
zlo=min(allz)*(1+extender)-max(allz)*extender
zhi=max(allz)*(1+extender)-min(allz)*extender
figure;plot([zlo,zhi],[0,0],'k-',...
            zV1*[1,1],0.5*x*[-1,1],'k-',...
            zV1prime*[1,1],0.5*x1*[-1,1],'k-',...
            zV2*[1,1],0.5*x*[-1,1],'k-',...
            zV2prime*[1,1],0.5*x*[-1,1],'k-',...
            zH1*[1,1],1.2*x*[-1,1],'g--',...
            zH1prime*[1,1],1.2*x*[-1,1],'g-.',...
            zH2*[1,1],1.2*x*[-1,1],'g--',...
            zH2prime*[1,1],1.2*x*[-1,1],'g-.',...
            zS1*[1,1],[0,x1],'m-',...
            zS2prime*[1,1],[0,x2prime],'b-',...
            [zS1,zH1],x*0.8*[1,1],'m-',...
            [zS2prime,zH2prime],x2prime*0.9*[1,1],'b-',...
            [zH1prime,zH2],x*0.9*[1,1],'k-',...
            zF1,0,'ko',...
            zF1prime,0,'ko',...
            zF2,0,'ko',...
            zF2prime,0,'ko');
text(zF1,-x*0.2,'F1');
text(zF1prime,-x*0.2,'F1''');
text(zF2,x*0.2,'F2');
text(zF2prime,x*0.2,'F2''');
text(zS1,x*1.2,'S1');set(moose,'color','magenta');
moose=text(zS1prime,x1prime+x*0.2,'S1''');set(moose,'color','blue');
moose=text(zS2prime,x2prime*1.2,'S2''');set(moose,'color','blue');
moose=text((zS1+zH1)/2,0.85*x,['s1=',num2str(s)]);set(moose,'color','magenta');
moose=text((zS2prime+zH2prime)/2,1.1*x2prime,'s2''');set(moose,'color','blue');
text((zH2+zH1prime)/2,0.95*x,['z12=',num2str(z12)]);
disp(['Plotting Figure ',num2str(gcf)]);
print -depsc 2lens2.eps
echo on;
%
% Plot the locations of principal planes and focal points for the whole
% compound lens.  The original matrix went from H1 to H2prime, so H is
% measured from H1 and Hprime from H2prime.  
%
echo off;
allz=[zH,zHprime,zH1,zH1prime,...
      zH2,zH2prime,zV1,zV1prime,zV2,zV2prime,zF,zFprime];
extender=0.05;
zlo=min(allz)*(1+extender)-max(allz)*extender;
zhi=max(allz)*(1+extender)-min(allz)*extender;
figure;plot([zlo,zhi],[0,0],'k-',...
            [zH1,zH],x*[1,1],'g-',...
            [zF,zH],0.2*x*[1,1],'k-',...
            [zFprime,zHprime],-0.2*x*[1,1],'k-',...
            [zH2prime,zHprime],-x*[1,1],'g-',...
            zV1*[1,1],0.5*x*[-1,1],'k-',...
            zV1prime*[1,1],0.5*x*[-1,1],'k-',...
            zV2*[1,1],0.5*x*[-1,1],'k-',...
            zV2prime*[1,1],0.5*x*[-1,1],'k-',...
            zH1*[1,1],1.2*x*[-1,1],'g--',...
            zH1prime*[1,1],1.2*x*[-1,1],'g-.',...
            zH2*[1,1],1.2*x*[-1,1],'g--',...
            zH2prime*[1,1],1.2*x*[-1,1],'g-.',...
            zF,0,'ko',...
            zFprime,0,'ko',...
            zH*[1,1],1.2*x*[-1,1],'g--',...
            zHprime*[1,1],1.2*x*[-1,1],'g-.');
moose=text((zH1+zH)/2,x*1.1,'-h');set(moose,'color','green');
moose=text((zH2prime+zHprime)/2,-x*1.1,'-h''');set(moose,'color','green');
text((zF+zH)/2,x*0.3,'f');
text((zFprime+zHprime)/2,-x*0.3,'f');
moose=text(zH,-x*1.3,'H');set(moose,'color','green');
moose=text(zHprime,-x*1.3,'H''');set(moose,'color','green');
text(zF,-x*0.2,'F');
text(zFprime,-x*0.2,'F''');
disp(['Plotting Figure ',num2str(gcf)]);
print -depsc 2lens3.eps
echo on;
%
%  Now use the lens equation, measuring s from the front principal plane
%
s=zH-zS
sprime=1/(1/f-1/s);sprime
m=-sprime/s
xprime=x*m
zSprime=zHprime+sprime
%
% Plot of object and image for the compound lens, showing object and image
% distances.
%
echo off;
allz=[zS,zSprime,zF,zFprime,zH,zHprime];
extender=0.05;
zlo=min(allz)*(1+extender)-max(allz)*extender;
zhi=max(allz)*(1+extender)-min(allz)*extender;
figure;plot([zlo,zhi],[0,0],'k-',...
            [zS,zH],0.9*x*[1,1],'m-',...
            [zSprime,zHprime],0.9*xprime*[1,1],'b-',...
            zV1*[1,1],0.5*x*[-1,1],'k-',...
            zV1prime*[1,1],0.5*x*[-1,1],'k-',...
            zV2*[1,1],0.5*x*[-1,1],'k-',...
            zV2prime*[1,1],0.5*x*[-1,1],'k-',...
            zS*[1,1],[0,x],'m-',...
            zSprime*[1,1],[0,xprime],'b-',...
            zF,0,'ko',...
            zFprime,0,'ko',...
            zH*[1,1],1.2*x*[-1,1],'g--',...
            zHprime*[1,1],1.2*x*[-1,1],'g-.');
moose=text((zS+zH)/2,x,'s');set(moose,'color','magenta');
moose=text((zSprime+zHprime)/2,xprime*(1.1),'s''');...
            set(moose,'color','blue');
moose=text(zH,-x*1.3,'H');set(moose,'color','green');
moose=text(zHprime,-x*1.3,'H''');set(moose,'color','green');
text(zF,-x*0.2,'F');
text(zFprime,-x*0.2,'F''');
moose=text(zS,x*1.2,'S');set(moose,'color','magenta');
moose=text(zSprime,xprime-x*0.2,'S''');set(moose,'color','blue');
disp(['Plotting Figure ',num2str(gcf)]);
print -depsc 2lens4.eps
echo on;
%
% Now let's look at apertures in object space
% Each lens is 1 cm in diameter.  Find the image of lens 2 as seen through
% lens 1. Treat lens 2 as an object and find its image.  Choose the first
% vertex to define the location of the lens.  Remember that we are
% reversing the normal direction of our diagrams so the light now goes from
% right to left.
%
dia1=1;
dia2=1;
s_2_1=zV2-zH1prime % Object distance
x_2_1=dia1/2
sprime_2_1=1/(1/f1-1/s_2_1)
zSprime_2_1=zH1-sprime_2_1 
xprime_2_1=-x_2_1*sprime_2_1/s_2_1
%
%  Plot the object and image location for lens 2 as an object
%
echo off;
allz=[zS,zV1,zV1prime,zV2,zV2prime,zH1,zH1prime,zSprime_2_1];
extender=0.05;
zlo=min(allz)*(1+extender)-max(allz)*extender;
zhi=max(allz)*(1+extender)-min(allz)*extender;

figure;plot([zlo,zhi],[0,0],'k-',...
            zV1*[1,1],0.5*x*[-1,1],'k-',...
            zV1prime*[1,1],0.5*x1*[-1,1],'k-',...
            zV2*[1,1],0.5*x*[-1,1],'k-',...
            zV2prime*[1,1],0.5*x*[-1,1],'k-',...
            zF1,0,'ko',...
            zF1prime,0,'ko',...
            zV2*[1,1],[0,x_2_1],'m-',...  % object
            zSprime_2_1*[1,1],[0,xprime_2_1],'b--')
disp(['Plotting Figure ',num2str(gcf)]);
print -depsc 2lens5.eps
echo on;
fov_edge=dia1/2*abs((zS-zSprime_2_1)/(zV1-zSprime_2_1))
%  The equation above is to extend the cyan line to the object.
%  This determines the largest object that fits the field of view.
%
%  Compute the angles  subtended by the apertures
%
tantheta=abs(xprime_2_1/(zSprime_2_1-zS))
na=sin(atan(tantheta))   % Numerical Aperture
tanfov=abs(dia1/2/(zV1-zSprime_2_1))
full_fov=2*atan(tanfov)  % Full field of view in radians
%
% Plot all apetures and angles
%
echo off;
allz=[zS,zV1,zV1prime,zV2,zV2prime,zH1,zH1prime,zSprime_2_1];
extender=0.05;
zlo=min(allz)*(1+extender)-max(allz)*extender;
zhi=max(allz)*(1+extender)-min(allz)*extender;
figure;plot([zlo,zhi],[0,0],'k-',...
            zV1*[1,1],dia1/2*[1,1.5],'k-',...  
            zV1*[1,1],-dia1/2*[1,1.5],'k-',...  
            zSprime_2_1*[1,1],xprime_2_1*[1,1.5],'k-',...  
            zSprime_2_1*[1,1],-xprime_2_1*[1,1.5],'k-',... 
            [zS,zV1],[0,dia1/2],'r:',...
            [zSprime_2_1,zV1],[0,dia1/2],'c-',...     % Cyan for window
            [zSprime_2_1,zS],...  
               [0,fov_edge],'c-',...  
            [zS,zSprime_2_1],[0,xprime_2_1],'r-',...  % Red for pupil
            zS*[1,1],[0,x],'m-');
%
moose=text(zS,x*1.2,'S');set(moose,'color','magenta');
text(zV1,dia1,'Lens1');
text(zSprime_2_1,xprime_2_1*1.5,'Lens2');
disp(['Plotting Figure ',num2str(gcf)]);
print -depsc 2lens6.eps
echo on;
