% Task1: Use ray transfer matrices in MATLAB to model the raysthrough the system. Demonstrate multiple rays (e.g. raysparallel to optical axis, oblique incidence, point source onopical axis, point source off-axis, or any other interestingcase). You should also model rays that are focused to asingle point in the Fourier plane (both on-axis and off-axiscase). 
%% Case1: Parallel rays that are parallel with optic axis.
f1 = 10; %focal length of the first lens 
f2 = 10; %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=linspace(-10, 10, 20);
angle=0*rays;

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);
subplot(2,2,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;

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 100]);
ylim([-80 30]);

title('Parallel rays that are parallel with optic axis');
%% Case2: Parallel rays with 30° input angle.
f1 = 10; %focal length of the first lens 
f2 = 10; %focal length of the second lens
o = 40; %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=linspace(-10, 10, 20);
angle=0*rays+(pi/6);

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);
subplot(2,2,2);

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;

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 100]);
ylim([-80 30]);

title('Parallel rays with -30° input angle.');
%% Case3: Parallel rays with -60° input angle.
f1 = 10; %focal length of the first lens 
f2 = 10; %focal length of the second lens
o = 40; %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=linspace(-10, 10, 20);
angle=0*rays-(pi/3);

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);
subplot(2,2,3);

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;

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 100]);
ylim([-80 30]);

title('Parallel rays with -60° input angle.');
%% Case4: Rays from a point source.
f1 = 10; %focal length of the first lens 
f2 = 10; %focal length of the second lens
o = 15; %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
  
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);
subplot(2,2,4);

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;

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 100]);
ylim([-80 30]);

title('Rays from a point source.');