% Spherical thick lenses
clear all; close all; clc;

% Coordinates
o  = 100; % object distance, light source distance
oh = 50;   % object height
oo = -oh/2;   % object vertical offset
i = o;  % image distance
dx = (i+0)/1e3;  % x resolution
x = -1.1*o:dx:i; % x-axis

% Optical parameters
n1 = 1.0; % refraction indice of start/end medium(air = 1.0)
n2 = 2.0; % refraction indice of lens
R1 = 100; % incoming curvature radius
R2 = 200; % outgoing curvature radius
d = 20; % lens thickness across lens axis
D = 55; % diameter of lens

% focus length
n = n2/n1;
f = 1 / ((n-1) * (1/R1 - 1/R2 + d*(n-1)/(n*R1*R2)));

% Lens
x1  = -d/2:dx:0;
cR1 = -d/2+R1; % center of incoming radius
y1 = sqrt(R1^2-(x1-cR1).^2);
for j = 1:length(y1); if y1(j) >= D/2; y1(j) = D/2; end; end; clear j
x2  = 0:dx:d/2;
cR2 = -R2+d/2; % center of incoming radius
y2 = sqrt(R2^2-(x2-cR2).^2);
for j = 1:length(y2); if y2(j) >= D/2; y2(j) = D/2; end; end; clear j

% % Horizontal Ray 1
% y0_1 = (oo + oh)
% x0_1 = -o
%     th1 = asin(y0_1/R1); theta1=th1/pi*180
% y1_1 = y0_1
% x1_1 = cR1-R1*cos(th1)
%     th2 = asin(1/n*sin(th1))-th1; theta2=th2/pi*180
% x2_1 = 0
% y2_1 = y1_1-x1_1*(th2)
%     L = sqrt((R2-d/2)^2+y2_1^2)
%     th3 = asin(y2_1/L); theta3 = th3/pi*180
%     th4 = pi-th3+th2; theta4=th4/pi*180
%     L2 = (L*cos(th4)+sqrt(L^2*(cos(th4))^2+R2^2-L^2))
% x3_1 = L2*cos(th2)
% y3_1 = y2_1+L2*sin(th2)

%% Light rays
% start conditions:
x0_1 = -o;
y0_1 = 5;
t0_1 = -2/180*pi;

% incoming crossing point on lens:
x1_1 = 1/2*(-2*cR1+2*t0_1*(-t0_1*x0_1+y0_1)-sqrt((2*cR1-2*t0_1*(-t0_1*x0_1+y0_1))^2-4*(-1-t0_1)*(R1^2-cR1^2-(-t0_1*x0_1+y0_1)^2)))/(-1-t0_1)
y1_1 = sqrt(R1^2-(cR1-x1_1)^2)
t1_1 = t0_1 + atan2(y0_1,cR1-x0_1) + asin(y1_1/R1)
t2_1 = t1_1*n1/n2;

% outgoing crossing point on lens:



theta0=t0_1/pi*180
theta1=t1_1/pi*180
% theta0=t0_1/pi*180
% theta0=t0_1/pi*180
% theta0=t0_1/pi*180



%% plots
% plot object
plot([-o -o],[oo oo+oh],'k','LineWidth',3)
hold on;

% plot lens and lens middle - part 1
plot(x1,y1,'k:')

% plot ray 1
plot([x0_1 x1_1],[y0_1 y1_1],'r')

% plot focus point
plot(f,0,'or','LineWidth',2)

% plot initial ray
plot([x0_1 x0_1+o*cos(t0_1)],[y0_1 y0_1+o*sin(t0_1)],'b:')

% plot x-axis
plot(x,zeros(1,length(x)),'k--');

% plot lens and lens middle - part 2
plot(x1,-y1,'k:')
plot(x2,y2,'k:',x2,-y2,'k:')
plot([0 0],[-D/2 D/2],'k--')

axis equal;
legend('Object','Lens','Light ray','Focus')
