clear;
closeall;

##################################################################################
N_layers = 3;
Nfreqs = 1;

d = [0; 50e-7; 0]; # air/ birefringent slab / air
f = c/451e-9;
theta = 0:0.1:89.9;


## Using the "diagonal" anisotropy
n1  = 1; # Incidence Medium, must be isotropic and loss-less
n3  = 1;


## permittivity matrix
eps_x = 2.25 + 0.01*1i; 
eps_y = 6.25+ 2.4e-4*1i; 
eps_z = 10.24;
#?sqrt(eps_x);
#?sqrt(eps_y);
#?sqrt(eps_z);
eps2_diag = zeros(3,3);
eps2_diag(1,1)=eps_x;
eps2_diag(2,2)=eps_y;
eps2_diag(3,3)=eps_z;


## define Euler angle based on Z1X2Z3 and generate a permittivity matrix
alpha = 28;  # degree
beta  = -35;
gamma = 90;

alpha = alpha*pi/180; # radian
beta = beta*pi/180;
gamma = gamma*pi/180;

s1 = sin(alpha); s2 = sin(beta); s3 = sin(gamma);
c1 = cos(alpha); c2 = cos(beta); c3 = cos(gamma);

rotationMatrix = [ c1*c3-c2*s1*s3,  -c1*s3-c2*c3*s1, s1*s2;
                   c3*s1+c1*c2*s3, c1*c2*c3-s1*s3, -c1*s2;
                   s2*s3, c3*s2, c2];
eps2_diag = mult(inv(rotationMatrix), mult(eps2_diag,rotationMatrix));


## Using the full anisotropy
n = matrix(N_layers, Nfreqs, 9);
n(1, :, 1) = n1; # Incidence Medium, must be isotropic and loss-less
n(1, :, 5) = n1; # Incidence Medium, must be isotropic and loss-less
n(1, :, 9) = n1; # Incidence Medium, must be isotropic and loss-less



## Convert the permittivity to refractive index 
for( i=1:9) {  
    n(2, :, i) = sqrt(eps2_diag(i)); 
    # n11 -> n(2, :, 1)
    # n21 -> n(2, :, 2)
    # n31 -> n(2, :, 3)
    # n12 -> n(2, :, 4)
    # n22 -> n(2, :, 5)
    # n32 -> n(2, :, 6)
}

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


RT = stackrt(n,d,f,theta);


plot(theta, RT.Rpp, RT.Rps, RT.Rss, RT.Rsp, RT.Tp, RT.Ts,'Theta [deg]','R/T','Fully anisotropic stackrt',"linewidth=2");
legend('R_pp','R_ps','R_ss','R_sp','T_p','T_s');
setplot('y min',0);
setplot('y max',1);
holdoff;