% Define lens parameters (All distances in mm)

f = 20; % This is a definition from paraxial approx.

n = 1.5;

R1 = 40;

R2 = (1/(f*(n-1))-1/R1)^-1;



T = 1; % Thickness on axis

H = fsolve(@(x) sqrt(R1^2-x^2)+sqrt(R2^2-x^2)-(R1+R2-T),T); % Corresponding element size (radius)

l1 = R1-sqrt(R1^2-H^2);

l2 = T-l1;



% Initial [position, angle (deg)]

%   Angle input is in deg, but kept in radians throughout

x0 = -4/4*H;

a0 = 0; a0 = radian(a0);



% Propagate distance z1, go through first interface

z1 = 10;

x1 = x0 + z1*tan(a0);

a1 = asin(1/n*sin(asin(x1/R1)+a0))-asin(x1/R1);



% Go through the thin/thick lens, past the interface

x2 = x1; % Thin

% x2 = something else; % Thick

a2 = asin(n*sin(a1-asin(x2/R2)))+asin(x2/R2);



% Propagate distance z3

z3 = 2*f;

x3 = x2 + z3*tan(a2);

a3 = a2;



% combine into nice vector form

zPos = [0 z1 z1+T z1+T+z3];

ray = [x0 a0;

       x1 a1;

       x2 a2;

       x3 a3;]';

ray(2,:) = deg(ray(2,:));



% Compare to paraxial

paraxRay = zeros(2,4);

paraxRay(:,1) = [x0 a0]';

paraxRay(:,2) = [1 z1; 0 1]*paraxRay(:,1);

paraxRay(:,3) = [1 0;-1/f 1]*paraxRay(:,2);

paraxRay(:,4) = [1 z3; 0 1]*paraxRay(:,3);

paraxRay(2,:) = deg(paraxRay(2,:));



%------------------------------------------------------------

% Draw results

%------------------------------------------------------------



% Draw the lens at z1

x = linspace(z1,z1+l1,1000);

y = sqrt(R1^2-(x-(z1+R1)).^2);

plot(x,y); hold on; plot(x,-y);

x = linspace(z1+l1,z1+T,1000);

y = sqrt(R2^2-(x-(z1+T-R2)).^2);

plot(x,y); plot(x,-y);

xlim([zPos(1) zPos(4)]);



% Draw ray

plot(zPos,ray(1,:),'r');

plot(zPos,paraxRay(1,:),'g');