% example of ray tracing for a two lens microscope or inverted telescope;
% inspired by a script by A. Cheville posted on the LASER CULT web-site

% define lenses  (use mm units)
f1 = 40;
f2 = 200;

% define distances between object and first lens, distance between first
% lens and second lens, and distance between second lens and image
d1 = f1;
d2 = f1+f2;
d3 = f2;

% these are the matrices that describe how a ray propogates
Md1 = [1 d1;0 1];
Md2 = [1 d2;0 1];
Md3 = [1 d3;0 1];
Mlens1 = [1 0;-1/f1 1];
Mlens2 = [1 0;-1/f2 1];

% rays are described by a two component vector made from the displacement from the
% optical axis and the angle with respect to the optical axis
% set this up for five rays that start out on the optical axis but have
% diffferent angles with respect to the optical axis
Rin = [zeros(1,5);-0.1:0.05:0.1];

% propogate the ray by matrix multiplication
R1 = Md1*Rin;   % after displacement of the distance to the first lens
R2 = Mlens1*R1; % after the first lens
R3 = Md2*R2;    % after displacement of the distance to the second lens
R4 = Mlens2*R3; % after  the second lens
Rout = Md3*R4;  % after displacement from the second lens to the detector

% want to make a plot of this so write down the relative position of each
% element in the problem
ds = [0 d1 0 d2 0 d3];
% then calaculte the absolute position along the optical axis
zpos = cumsum(ds);
% do the plot as a loop (probably a more elegant way to do this but I'm
% just learning)
hold on;
for k=1:5;
  % create a vector that describes the position of the ray at each point
  rpos = [Rin(1,k) R1(1,k) R2(1,k) R3(1,k) R4(1,k) Rout(1,k)];
  % and plot the ray going through the system
  plot(zpos,rpos);
end
xlabel('Z (mm)');
ylabel('r (mm)');
axis([ 0 2*(f1+f2) -10 10])

%  do this again but with rays that leave at some point 1 mm off the optical
%  axis

Rin = [ones(1,5);-0.1:0.05:0.1];

% propogate the ray by matrix multiplication
R1 = Md1*Rin;   % after displacement of the distance to the first lens
R2 = Mlens1*R1; % after the first lens
R3 = Md2*R2;    % after displacement of the distance to the second lens
R4 = Mlens2*R3; % after  the second lens
Rout = Md3*R4;  % after displacement from the second lens to the detector

% do the plot as a loop
hold on;
for k=1:5;
  % create a vector that describes the position of the ray at each point
  rpos = [Rin(1,k) R1(1,k) R2(1,k) R3(1,k) R4(1,k) Rout(1,k)];
  % and plot the ray going through the system, this time use a red line
  plot(zpos,rpos,'r');
end