%Example 4.1: reflection off a mirror
%pol_locus_plot.m and calc_ellipticity.m and
%find_global_reflection_from_surface.m
%must be on the Matlab path for 
%this script to function properly

%coordinates are    z in horiz w pos direction on right
%                   y vertical w pos up
%                   x into plane of drawing

%Parameters
N_i         = 1;
N_t         = 1.02+6.85j;
J1          = [1;-j];       %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     = -125*pi/180;  %define surface normal direction
                            %assume in yz plane - this is rotation from y axis

%plot E locus for j1
state = [J1(1) J1(2)];
figure(1);pol_locus_plot
title('Locus plot for J1')
disp('FOR INCIDENT WAVE:')
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, J2, a2] ...
    = find_global_reflection_from_surface(J1,N_i,N_t,khat_i,nhat);

%plot E locus for J2
state = [J2(1) J2(2)];
figure(2);pol_locus_plot
title('Locus plot for J2')
disp('FOR REFLECTED WAVE:')
calc_ellipticity

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

%find visibility
I1          = abs(norm(J1))^2;
I2          = abs(norm(J2))^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;%angle(a1*conj(a2));
yvec        = I1 + I2 + 2*sqrt(I1*I2)*p12*cos(2*pi*xvec + phi_delta +phi_p);
figure(3);plot(xvec,yvec);axis([0 xlim 0 max(yvec)]);grid
title('Interferogram Profile')
xlabel('Distance Along k_\Delta in Waves')
ylabel('Relative Irradiance')