%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cooke_triplet.m
%
% This program can be used to initilize a triplet lens design.
% Algoritm is extracted from the book Lens Design, Laikin, Page 54, 4th Ed.
%
% We assume that
%   - the system has PNP structure and
%   - material of 1st and 3rd lens are the same.
%
% Inputs:
%    f      System focal length
%    f/#    f-number
%    Na,Va  Refractive index and Abbe number of first lens.
%    Nb,Vb  Refractive index and Abbe number of second lens.
%
% Outputs:
%    Pi     Power of each lens (i=1,2,3)
%    D      Enterance pupil diameter
%    Rij    Radius of curvatures of the lenses (i,j=1,2,3)
%    T1     Thickness between 1st and 2nd lens
%    T2     Thickness between 2nd and 3nd lens
%    T3     Thickness between 3rd lens and image plane
%
% Apr 2023
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clc;

% refractive indices and abbe values
% for BK7 and sf5 glasses at d-line
Na = 1.5168; Va = 64.17;
Nb = 1.6727; Vb = 32.25;
Nc = Na;     Vc = Va;
% desired focal length and f/#
f  = 100;
fno= 3.5;

% calculations start here -------
D  = f/fno;
phi= 1/f;
P  = 1/f;
TR = 1;
T1 = 0.1/P;
X  = 0.8;
T2 = TR*T1;
E  = 0.0001;

for i=1:10
   Pc = (1/X - sqrt( (X/TR^2+1)*Vb*Na/(Va*X*Nb*(X+1)) ) ) / T1; Fc = 1/Pc;
   Pb = -Pc*Nb*(X+1)/Na; Fb=1/Pb;
   Pa = X*Pc; Fa = 1/Pa;
   T3 = ((Fa-T1)*Fb-(Fa+Fb-T1)*T2)*Fc / ((Fa-T1)*Fb+(Fa+Fb-T1)*(Fc-T2));
   P  = Pa+Pb+Pc - T1*Pa*(Pb+Pc)-T2*Pc*(Pa+Pb)+T1*T2*Pa*Pb*Pc;
   X1 = T3*P*TR;
   if abs(X1-X)-E < 0
       fprintf('Loop is terminated since |X1-X| = %f < %f\n\n',abs(X1-X),E)
       break;
   end
   X = X1;
end

S  = P/phi;
Pa = Pa/S;
Pb = Pb/S;
Pc = Pc/S;
T1 = T1*S;
T3 = T3*S;
T2 = T1*TR;

R11= (Na-1)/Pa;
R12= inf;
R21= -2*(Nb-1)/Pa;
R22= -R21;
R31= inf;
R32= (Na-1)/Pc;

% Outputs to the user screen -------
fprintf('*** Brief Prescription ***\n')
fprintf('Powers  : P1  = %8.3f P2 = %8.3f  P3 = %8.3f\n',Pa,Pb,Pc)
fprintf('EFFL    : F   = %8.3f\n',f)
fprintf('F/number: F/# = %8.3f\n',fno)
fprintf('ENPD    : D   = %8.3f\n\n',D)

fprintf('*** Lens Data ***\n')
fprintf('# Surface     Radius    Thickness    Index   Abbe\n');
fprintf('0 Object    %8.3f     %8.3f \n',inf,inf);
fprintf('1 Standard  %8.3f     %8.3f %8.3f %8.3f\n',R11,4,Na,Va);
fprintf('2 Standard  %8.3f     %8.3f\n',R12,T1);
fprintf('3 Standard  %8.3f     %8.3f %8.3f %8.3f\n',R21,4,Nb,Vb);
fprintf('4 Standard  %8.3f     %8.3f\n',R22,T2/2);
fprintf('4 Stop      %8.3f     %8.3f\n',inf,T2/2);
fprintf('5 Standard  %8.3f     %8.3f %8.3f %8.3f\n',R31,4,Nc,Vc);
fprintf('6 Standard  %8.3f     %8.3f\n',R32,T3);
fprintf('7 Image     %8.3f     %8.3f\n',inf,0);
