# Layer thickness:
d_a = 226.0e-9;
d_b = 220.0e-9;
d_PMMA = 100.0e-9;

N_per = 5;

d_PC = N_per*(d_a+d_b);

# source angles:

src_theta_TIR = [43.0:0.1:43.9];
src_theta_rest = [44.0:1:48.0];
src_theta_dip = [49:0.1:51];

TIR_len = length(src_theta_TIR);
rest_len = length(src_theta_rest);
dip_len = length(src_theta_dip);

src_theta = matrix(TIR_len+rest_len+dip_len,1);
src_theta(1:TIR_len) = src_theta_TIR;
src_theta(TIR_len+1:TIR_len+rest_len) = src_theta_rest;
src_theta(TIR_len+rest_len+1:TIR_len+rest_len+dip_len) = src_theta_dip;

# T/R results to collect:
R_num = matrix(length(src_theta));
T_num = matrix(length(src_theta));

n_E2 = 400;
xgrid_E2 = linspace(-d_PC,300e-9,n_E2);
E2 = matrix(n_E2,length(src_theta));

# ------------
# Simulations:
# ------------

# source angle sweep:
for(I = 1:length(src_theta))
{
    switchtolayout;
    ?msg = "Source angle: " + num2str(src_theta(I));
    setnamed("PlaneWaveSrc","angle theta",src_theta(I));
    run;
    
    R_num(I) = -transmission("BackPwrMonitor");
    T_num(I) = transmission("FrontPwrMonitor");

    E2(1:n_E2,I) = interp(getelectric("Field_Monitor"),getdata("Field_Monitor","x"),xgrid_E2);
}

# ==================================
# THEORETICAL CALCULATION (STACKRT):
# ==================================

fsim = getdata("PlaneWaveSrc","f"); #frequency used in simulation
?msg = "Used source f = " + num2str(fsim);

graph_name = "C (graphene) - broadband";
graph_thick = 0.1e-9; #Tickness of graphene (for stackrt)

add_graphene = true;
add_PMMA = true;

# Refractive indices:
index_a = 1.4468; #SiO2
index_b = 2.7204; #TiO2
index_cladding = 1; #air
index_substrate = 1.4468; #SiO2
index_PMMA = 1.481;

# Calculate refractive index of graphene:
graph_cond =getsurfaceconductivity(graph_name,fsim);
?msg = "Graphene: cond = " + num2str(graph_cond);

imag_unit = 0+1i;
omega_val = 2*pi*fsim;
perm_val = 1 + imag_unit*graph_cond/(eps0*omega_val*graph_thick);
index_val = sqrt(perm_val);

# Use stackrt:
nd = 2*N_per + 2; 

if (add_PMMA){
 ?msg = "Added PMMA layer";
 nd = nd +1;
}

if (add_graphene){
 ?msg = "Added graphene layer";
 nd = nd +1;
}

d_layers = matrix(nd,1);

nf = length(fsim);
index_layers = matrix(nd,nf);

d_layers(1) = 0; #Substrate
d_layers(nd) = 0; #Cladding

index_layers(1) = index_substrate; #Substrate
index_layers(nd) = index_cladding; #Cladding

if (add_PMMA){
 index_layers(nd-1) = index_PMMA;
 d_layers(nd-1) = 100.0e-9;

if (add_graphene){
 index_layers(nd-2) = index_val;
 d_layers(nd-2) = graph_thick;
}

}
else{

if (add_graphene){
 index_layers(nd-1) = index_val;
 d_layers(nd-1) = graph_thick;
}

}

if (N_per>0){

for (I = 1:N_per){ 
 layer_b = 2*I;
 layer_a = layer_b+1;

 d_layers(layer_a) = d_a;
 d_layers(layer_b) = d_b;

 index_layers(layer_a) = index_a;
 index_layers(layer_b) = index_b;
}

}

RT = stackrt(index_layers,d_layers,fsim,src_theta);

# ======
# PLOTS:
# ======

plot(src_theta,R_num,RT.Rs,T_num,RT.Ts,"Source angle (deg)", "Power");
legend("Reflected - Numerical","Reflected - Analytical","Transmitted - Numerical","Transmitted - Analytical");

leg=cell(length(src_theta));
for (I=1:length(src_theta))
{
 leg{I}=num2str(src_theta(I))+"deg";
}

image(xgrid_E2*1e6,src_theta,E2,"x (um)","Source angle (deg)");