% Problem 2: Write a simple Matlab code to calculate the imaging of one of the designed lenses using geometrical optics and Fermat’s principle. 
   
f1 = 10; %focal length of the first lens 
f2 = 1.5; %focal length of the second lens
o = 5; %the object's distance from the first lens
d = f1+f2;  % distance of the lenses
f = 1/(1/f1+1/f2-d/f1/f2);
s=abs((o*f)/(o-f))*2;  % distance of the screen

if(f==Inf) 
    s=f1;
end

if(s==Inf) 
    s=f;
end

T_o=[1, o; 0, 1];       % object translation matrix
L1=[1, 0; -1/f1, 1];   % lens transfer matrix
TD=[1, d; 0, 1];       % separation translation matrix
L2=[1, 0; -1/f2, 1];   % lens transfer matrix
T_s=[1, s; 0, 1];       % screen translation matrix

rays = 20;
angle=[-(f1+rays)/o:f1/o/10:(f1-rays)/o]; % angle
rays=rays(1)*(angle*0+1);     % lateral position


v_0=[rays;angle]; 
v_1=T_o*v_0; % object to lens1
v_2=L1*v_1; % lens1
v_3=TD*v_2; % lens1 to lens2
v_4=L2*v_3; % lens2
v_5=T_s*v_4; % lens2 to screen

figure(1);

for i = 1:length(angle)
    plot([-o,0,0,d,d,s+d], [v_0(1,i),v_1(1,i),v_2(1,i),v_3(1,i),v_4(1,i),v_5(1,i)], 'r-');
    hold on;
end

grid on;
%axis equal
%axis tight
xlabel('x axis');
ylabel('y axis');

plot([0,0],[-f1,f1],'g.-','linewidth',2);
plot([d,d],[-f2,f2],'g.-','linewidth',2);
plot([-f1,f1,d-f2,d+f2],[0,0,0,0],'k*','linewidth',2); 
plot([-o-f1,s+d],[0,0],'k-','linewidth',1);       % plot focal point
hold off;

xlim manual;
ylim manual;

xlim([-20 20]);
ylim([-80 30]);

%if((xl(2)-xl(1))/1>(yl(2)-yl(1)))
%    xlim([-2*f1,d+2*f2])
%end
%axis equal
%axis tight

