function [p,q,r,s,v,w,f1,f2] = Thick_Lens_gen(nn1,nnL,nn2,RR1,RR2,dthick)

% -------------------------------------------------------------------------
%  Example:
%   [p,q,r,s,v,w,f1,f2] = Thick_Lens_gen(1,1.50,1,50,-50,5);
% -------------------------------------------------------------------------
%   [p,q,r,s,v,w,f1,f2] = Thick_Lens(n1,nL,n2,R1,R2,dthick)
%
%    Calculate cardinal points from A, B, C, D
%
%    It calculates symbolically the A, B, C, D (it needs symbolic Toolbox)
%    and then it determines the numerical values
%
%    Input:
%
%    nn1    = refractive index of the left surrounding the lens material
%    nnL    = refractive index of the lens material
%    nn2    = refractive index of the right surrounding the lens material
%    RR1    = Radius of curvature of the left spherical lens surface (mm)
%    RR2    = Radius of curvature of the right spherical lens surface (mm)
%    dthick = Lens thickness (mm)
%    
%    Output:  p,q,r,s,v,w,f1,f2 (all in mm)
% -------------------------------------------------------------------------

%
%    Determine symbolically the A, B, C, D elements 
%
syms  n1 nL n2 R1 R2 d

M1 = [1 0; ((n1/nL)-1)*(1/R1) (n1/nL)];
M2 = [1 d; 0 1];
M3 = [1 0; ((nL/n2)-1)*(1/R2) (nL/n2)];

MM = M3*M2*M1;
ABCD = simplify(MM)

AA = MM(1,1);
AA = simplify(AA);
BB = MM(1,2);
BB = simplify(BB);
CC = MM(2,1);
CC = simplify(CC);
DD = MM(2,2);
DD = simplify(DD);

A  = subs(AA,{n1,nL,n2,R1,R2,d},{nn1,nnL,nn2,RR1,RR2,dthick});
B  = subs(BB,{n1,nL,n2,R1,R2,d},{nn1,nnL,nn2,RR1,RR2,dthick});
C  = subs(CC,{n1,nL,n2,R1,R2,d},{nn1,nnL,nn2,RR1,RR2,dthick});
D  = subs(DD,{n1,nL,n2,R1,R2,d},{nn1,nnL,nn2,RR1,RR2,dthick});

ABCD_num = [A B; C D]

%   A  =  1 + ((n1/nL)-1)*(dthick/R1);
%   B  =  dthick * (n1/nL);
%   C  =  -((nL/n1)-1) * ((1/R1)-(1/R2)) - (nL-n1)^2*dthick/(n1*nL*R1*R2);
%   D  =  (1-(n1/nL))*dthick/R2 +1;


p  = D/C;
q  = -A/C;
n1n2 = nn1/nn2;
r  = (D-n1n2)/C;
s  = (1-A)/C;
v  = (D-1)/C;
w  = (n1n2-A)/C;
f1 = n1n2/C;
f2 = -1/C;

p = double(p);
q = double(q);
r = double(r);
s = double(s);
v = double(v);
w = double(w);
f1 = double(f1);
f2 = double(f2);