clear;
switchtolayout;

# graphene material parameters:
graph_name = "C (graphene) - broadband";
graph_temp = 300.0;
# frequency range:
freq_min = 15e12;
freq_max = 45e12;
freq_samples = 400;
# Set frequency range:
setglobalsource("set frequency", true);
setglobalsource("frequency start", freq_min);
setglobalsource("frequency stop", freq_max);
freq_val = linspace(freq_min,freq_max,freq_samples);

# ------------------------
# Full Conductivity Model:
# ------------------------

#setmaterial(addmaterial("Graphene"),"name",graph_name);
setmaterial(graph_name,"temperature (K)",graph_temp);
setmaterial(graph_name,"conductivity scaling",1);

# case 1:
graph_gama = 0.00102;
graph_mu = 0.265;
setmaterial(graph_name,"scattering rate (eV)",graph_gama);
setmaterial(graph_name,"chemical potential (eV)",graph_mu);
setmaterial(graph_name,"frequency samples", freq_samples);
# case 1 conductivity data:
cond_val1=getfdtdsurfaceconductivity(graph_name,freq_val,freq_min,freq_max);

# case 2:
graph_gama = 0.00099;
graph_mu = 0.217;
setmaterial(graph_name,"scattering rate (eV)",graph_gama);
setmaterial(graph_name,"chemical potential (eV)",graph_mu);
# case 2 conductivity data:
cond_val2=getfdtdsurfaceconductivity(graph_name,freq_val,freq_min,freq_max);

# -------------------------
# Drude Conductivity Model:
# -------------------------

imag_unit = 0+1i;

# case 1:
chu_mu = 0.25*e; 
chu_tau = 0.3e-12;
cond_chu1 = 1i*(e^2)*chu_mu/(pi*(hbar^2)*(2*pi*freq_val+(imag_unit/chu_tau)));

# Case 2:
chu_mu = 0.2*e;
chu_tau = 0.3e-12;
cond_chu2 = 1i*(e^2)*chu_mu/(pi*(hbar^2)*(2*pi*freq_val+(imag_unit/chu_tau)));

# ------
# Plots:
# ------

plot(freq_val*1e-12,1e6*real(cond_chu1),1e6*real(cond_val1),1e6*real(cond_chu2),1e6*real(cond_val2),"Frequency (THz)", "Re(sigma) (uS)");
legend("Drude, Ef = 0.25eV","Full, muc = 0.265eV","Drude, Ef = 0.20eV","Full, muc = 0.217eV");

plot(freq_val*1e-12,1e6*imag(cond_chu1),1e6*imag(cond_val1),1e6*imag(cond_chu2),1e6*imag(cond_val2),"Frequency (THz)", "Im(sigma) (uS)");
legend("Drude, Ef = 0.25eV","Full, mcu = 0.265eV","Drude, Ef = 0.20eV","Full, muc = 0.217eV");