%Example 4.2: Laser Isolator
%pol_locus_plot.m must be on the Matlab path for this script to function
%properly

%Parameters
N_i         = 1;
N_t         = 1.02+6.85j;
J1a         = [1;0];       %incident Jones vector for wave 1
theta_ki    = 0*pi/180;   %define incident k vector direction
                            %assume in yz plane  -this is rotation from z
                            %axis
theta_n     = -90*pi/180;  %define surface normal direction
                            %assume in yz plane - this is rotation from y axis

%plot E locus for j1a
state = [J1a(1) J1a(2)];
figure(1);subplot(3,2,1);pol_locus_plot
title('Locus of j1a, on LHS of QWP')
disp(' ')
disp('FOR j1a, on LHS of QWP:')
calc_ellipticity

%Define Rotated QWP matrix
QWP_45      = pol_rotate(QWP,45);

%Define state J1b with Jones Calculus
J1b         = QWP_45*J1a;

%plot E locus for j1b
state = [J1b(1) J1b(2)];
subplot(3,2,2);pol_locus_plot
title('Locus of j1b, on RHS of QWP')
disp(' ')
disp('FOR j1b, on RHS of QWP:')
calc_ellipticity

%define incident k vector direction
% theta_ki    = -30*pi/180;  %assume in yz plane  -this is rotation from z axis
alpha_i     = 0;
beta_i      = sin(theta_ki);
gamma_i     = cos(theta_ki);
khat_i      = [alpha_i beta_i gamma_i];

% %define surface normal direction
% theta_n     = -135*pi/180;  %assume in yz plane - this is rotation from y axis
alpha_n     = 0;
beta_n      = cos(theta_n);
gamma_n     = sin(theta_n);
nhat        = [alpha_n beta_n gamma_n];

[khat_r, xhat_p, yhat_pi, yhat_pr, MP1, MP2, a1, J2a, a2] ...
    = find_global_reflection_from_surface(J1b,N_i,N_t,khat_i,nhat);

%plot E locus for J2a
state       = [J2a(1) J2a(2)];
subplot(3,2,3);pol_locus_plot
title('Locus of j_2a, after relfection')
disp(' ')
disp('FOR j2a, after reflection:')
calc_ellipticity

%find p12 and phi_p
p12         = abs(dot(a1,a2));
phi_p       = angle(dot(a1,a2));

%find visibility
I1          = abs(norm(J1b))^2;
I2          = abs(norm(J2a))^2;
V           = 2*sqrt(I1*I2)/(I1+I2)*p12;
disp(['I1 = ' num2str(I1)])
disp(['I2 = ' num2str(I2)])
disp(['p12 = ' num2str(p12)])
disp(['V = ' num2str(V)])

%plot interferogram profile 
xlim        = 10;
xvec        = linspace(0,xlim,500);
phi_delta   = 0;
yvec        = I1 + I2 + 2*sqrt(I1*I2)*p12*cos(2*pi*xvec + phi_delta +phi_p);
subplot(3,2,4);plot(xvec,yvec);axis([0 xlim 0 max(yvec)]);grid
title('Interferogram Between QWP and Mirror')
xlabel('Distance Along k_\Delta in Waves')
ylabel('Relative Irradiance')

%Define Rotated QWP matrix (second pass)
QWP_m45     = pol_rotate(QWP,-45);

%Define state J2b with Jones Calculus
J2b         = QWP_m45*J2a;

%plot E locus for J2b
state       = [J2b(1) J2b(2)];
subplot(3,2,5);pol_locus_plot
title('Locus of j_2b, after 2nd pass through QWP')
disp(' ')
disp('FOR j2b after 2nd pass through QWP:')
calc_ellipticity

%find p12 and phi_p
a1          = J1a/norm(J1a);
a2          = J2b/norm(J2b);
p12         = abs(dot(a1,a2));
phi_p       = angle(dot(a1,a2));

%find visibility
I1          = abs(norm(J1b))^2;
I2          = abs(norm(J2a))^2;
V           = 2*sqrt(I1*I2)/(I1+I2)*p12;
disp(['I1 = ' num2str(I1)])
disp(['I2 = ' num2str(I2)])
disp(['On LHS of QWP, p12 = ' num2str(p12)])
disp(['On LHS of QWP, V   = ' num2str(V)])

%plot interferogram profile for waves on LHS of QWP
xlim        = 10;
xvec        = linspace(0,xlim,500);
phi_delta   = 0;
yvec        = I1 + I2 + 2*sqrt(I1*I2)*p12*cos(2*pi*xvec + phi_delta + phi_p);
subplot(3,2,6);plot(xvec,yvec);axis([0 xlim 0 max(yvec)]);grid
title('Interferogram Profile on LHS of QWP')
xlabel('Distance Along k_\Delta in Waves')
ylabel('Relative Irradiance')