# load project file
load("graphene_THz_metamaterial.fsp");

# graphene related parameters:
graph_name = "C (graphene) - broadband";
graph_mu = [0.265,0.217];
graph_gama = [0.00102,0.00099];
graph_layers = [1,2,4];

# cell to collect transmission results:
T_power = cell(6);
T_norm_fac = 0.0;

# run simulations:
for(J=1:3)
{
  for (I=1:2)
  {
    switchtolayout;
    setmaterial(graph_name,"scattering rate (eV)",graph_gama(I));
    setmaterial(graph_name,"chemical potential (eV)",graph_mu(I));
    setmaterial(graph_name,"conductivity scaling",graph_layers(J));
    run;
    T_power{I+length(graph_mu)*(J-1)}=getresult("Transmission-Monitor","T");
    T_norm_fac_tmp = max(T_power{I+length(graph_mu)*(J-1)}.T);
    if(T_norm_fac_tmp>T_norm_fac)
    {
      T_norm_fac = T_norm_fac_tmp;
    }
  }
}

# transmitted power plot:
plot(T_power{1}.lambda*1e6,T_power{1}.T/T_norm_fac);
holdon;
plot(T_power{2}.lambda*1e6,T_power{2}.T/T_norm_fac);
plot(T_power{3}.lambda*1e6,T_power{3}.T/T_norm_fac);
plot(T_power{4}.lambda*1e6,T_power{4}.T/T_norm_fac);
plot(T_power{5}.lambda*1e6,T_power{5}.T/T_norm_fac);
plot(T_power{6}.lambda*1e6,T_power{6}.T/T_norm_fac);
setplot("x label","Wavelength (microns)");
setplot("y label","Transmission");
legend("muc = 0.265, N = 1","muc = 0.217, N = 1","muc = 0.265, N = 2","muc = 0.217, N = 2","muc = 0.265, N = 4","muc = 0.217, N = 4");
holdoff;