%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ynu.m
%
% Simple y-nu ray tracing program for N surfaces
%
% As an example, we trace first patented Cooke Triplet with
% f/#  = 4 and EFFL = 100 mm. Note that the seecond surface
% of the second lens (in the middle) is STOP.
%
% Inputs:
%           R  Radius of curvature of each surface
%           n  Index of refraction between surfaces
%           t  Thickness (distance) between surfaces
%    z0,y0,u0  Initial z-position, height and angle
%
% Outputs:
%        EFFL  Effective focal length of the system
%    zk,yk,uk  Paraxial ray trace data at kth iteration
%
% Nov 2021
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clc;

% Radius of curvature of each surface
% For OBJ and IMA planes R = inf.
R = [inf 26.4 -150.7 -29.8 24.2 150.7 -26.4 inf];

% index of refraction between surfaces (for air n=1)
n = [1 1.5108 1 1.6042 1 1.5108 1 1];

% Thickness (distance) between surfaces
t = [1 5.9 10.9 0.2 12.5 5.9 81];
N = length(t);

% Initial values (z-position, height and angle)
z(1) = 0.0;
y(1) = 12.5; % semi-diameter of entrance pupil
u(1) = 0.0;  % parallel ray

% y-nu trace
for  k = 1:N
    z(k+1) = z(k) + t(k);
    y(k+1) = y(k) + u(k)*t(k);
    u(k+1) = (n(k)*u(k) - (n(k+1)-n(k))*y(k+1)/R(k+1))/n(k+1);
end
EFFL = -y(1)/u(end);

% print the results
fprintf('EFFL = %f\n',EFFL)
result = [z' y' u']
