# scriptfile: plotcs.lsf
# This script calculates the scattering, absorption 
# and extinction cross-sections in 2D.

# get cross-sections from analysis objects
sigmascat = getresult("scat","sigma");
sigmaabs = getresult("total","sigma");
sigmaext = -(sigmaabs.sigma) + (sigmascat.sigma);
lam = c/(sigmascat.f)*1e9;#in nanometer

#load cross section theory from text file
tmp=readdata("nanowire_theory.csv");
lam_theory=tmp(1:100,1)/1e9;
f_theory = c/lam_theory;
r23=tmp(1:100,2:4)*2*23*1e-9;
r24=tmp(1:100,5:7)*2*24*1e-9;
r25=tmp(1:100,8:10)*2*25*1e-9;
r26=tmp(1:100,11:13)*2*26*1e-9;
r27=tmp(1:100,14:16)*2*27*1e-9;

for (i=1:3) {  #interpolate onto FDTD's frequency points
	r23(1:100,i)=interp(r23(1:100,i),f_theory,sigmascat.f);
	r24(1:100,i)=interp(r24(1:100,i),f_theory,sigmascat.f);
	r25(1:100,i)=interp(r25(1:100,i),f_theory,sigmascat.f);
	r26(1:100,i)=interp(r26(1:100,i),f_theory,sigmascat.f);
	r27(1:100,i)=interp(r27(1:100,i),f_theory,sigmascat.f);
}

# Compare FDTD to theory
plot(lam, sigmaext*1e9, -(sigmaabs.sigma)*1e9, sigmascat.sigma*1e9, 
          r25(1:100,1)*1e9, r25(1:100,2)*1e9, r25(1:100,3)*1e9, 
          "wavelength (nm)","cross section (nm)","Cross sections");
legend("Extinction", "Absorption", "Scattering", "r25 theory", "r25 theory", "r25 theory");
# compare FDTD to theory of 24 and 26nm particle
plot(lam, sigmaext*1e9, -(sigmaabs.sigma)*1e9, sigmascat.sigma*1e9, 
          r24(1:100,1)*1e9, r24(1:100,2)*1e9, r24(1:100,3)*1e9, 
          r26(1:100,1)*1e9, r26(1:100,2)*1e9, r26(1:100,3)*1e9, 
          "wavelength (nm)", "cross section (nm)", "Cross sections");
legend("Extinction", "Absorption", "Scattering", "r24 theory", "r24 theory", 
       "r24 theory", "r26 theory", "r26 theory", "r26 theory");