clear;
closeall;

####################################################################
## material settings

N_layers = 5;
Nfreqs = 11;

n1  = 1.5; # Incidence Medium, must be isotropic and loss-less
n5=1;
nx_a = 1.511;
ny_a = 1.5095;
nz_a = nx_a-(nx_a-ny_a)*1;
ne_p = 1.5+1i*0.019027;
no_p = 1.5+1i*0.000306;

f = c/linspace(400e-9,700e-9,Nfreqs);
d = [0; 24e-6;0.55e-6/2/(nx_a-ny_a);0.55e-6/4/(nx_a-ny_a); 0];
# glass(1.5), polarizer, 1/2 plate, 1/4 plate, air

theta_out = 0:1:85;# angles in air
theta = asin(sin(theta_out*pi/180)/n1)*180/pi;  # angles in glass


phi2 = 90*pi/180; # slow axis angle of polarizer
phi3 = 15*pi/180; # slow axis angle of 1/2 wave plate
phi4 = 75*pi/180; # slow axis angle of 1/4 wave plate
####################################################################
## Derive refractive indices
n = matrix(N_layers, Nfreqs, 9);

eps2_diag = [ne_p^2,0,0;0,no_p^2,0;0,0,no_p^2];
eps3_diag = [nx_a^2,0,0;0,ny_a^2,0;0,0,nz_a^2];
eps4_diag = eps3_diag;


rotation_matrix =[cos(phi2), sin(phi2),0;-sin(phi2),cos(phi2),0;0,0,1]; #rotate around z axis;
eps2_diag = mult(transpose(rotation_matrix), mult(eps2_diag,rotation_matrix));

rotation_matrix =[cos(phi3), sin(phi3),0;-sin(phi3),cos(phi3),0;0,0,1]; #rotate around z axis;
eps3_diag = mult(transpose(rotation_matrix), mult(eps3_diag,rotation_matrix));

rotation_matrix =[cos(phi4), sin(phi4),0;-sin(phi4),cos(phi4),0;0,0,1]; #rotate around z axis;
eps4_diag = mult(transpose(rotation_matrix), mult(eps4_diag,rotation_matrix));

n(1, :, 1) = n1; # Incidence Medium, must be isotropic and loss-less
n(1, :, 5) = n1;
n(1, :, 9) = n1;

n(5, :, 1) = n5; # Outgoing material can be lossy
n(5, :, 5) = n5; 
n(5, :, 9) = n5;

## Convert the permittivity to refractive index 
for (i=1:9) {  
    n(2, :, i) = sqrt(eps2_diag(i)); 
    n(3, :, i) = sqrt(eps3_diag(i));
    n(4, :, i) = sqrt(eps4_diag(i));  
}
####################################################################
## stack solver
RT = stackrt(n,d,f,theta);