########################################################################
# Script file: sp_ebeam.lsf
# This script file calculates and plots the scattering spectrum of 
# light from an electron beam. It will run 2 simulations, one of which is 
# a reference simulation.
########################################################################

run_simulations = 1;
save; 	# set to 1 to run simulations, 0 to load existing simulations.

fname = "sp_ebeam.fsp";
fname_ref = "sp_ebeam_ref.fsp";

# run simulations
if(run_simulations){
	load(fname);
	switchtolayout;
	save(fname);
	run;
	switchtolayout;
	save(fname_ref);
	select("structures");
	set("enabled",0);
	run;
}

# load the original simulation
load(fname);

# switch to no normalization state
nonorm; 

# collect the data
mname = "fieldXY";
E = getresult(mname,"E");
H = getresult(mname,"H");

# load the reference simulation
load(fname_ref);

# switch to no normalization state
nonorm;

# collect the data
E0 = getresult(mname,"E");
H0 = getresult(mname,"H");


# take difference of fields
Ex = E.Ex - E0.Ex;
Ey = E.Ey - E0.Ey;
Ez = E.Ez - E0.Ez;
Hx = H.Hx - H0.Hx;
Hy = H.Hy - H0.Hy;
Hz = H.Hz - H0.Hz;


# get position and frequency vectors
x = E.x;
y = E.y;
z = E.z;
f = E.f;
w = 2*pi*f;

# take the spectrum of 1 source to remove the effects of the finite source pulse
spectral_norm = abs(sourcenorm(f,"ebeam::s1"))^2;

# calculate the Poynting vector
Pz = (Ex*conj(Hy) - Ey*conj(Hx));

# Calculate the transmission, as it will be observed experimentally
T = integrate(0.5*real(Pz),1:2,x,y)/spectral_norm/w^2;

# units are arbitrary
T = T/max(T);

# plot results
plot(c/f*1e9,T,"wavelength (nm)","T_upper (a.u.)");

# image the Poynting vector at the peak transmission
p0 = findpeaks(T);
image(x*1e6,y*1e6,pinch(real(Pz),4,p0),"x (microns)","y (microns)","real(Pz) (a.u.) at " + num2str(c/f(p0)*1e9)+ "nm");


# calculate far field now at the peak wavelength
res = 201;
load(fname);
Efar = farfieldpolar3d(mname,p0,res,res,1,1,1,1,1); # the last "+l" is to make sure forward projection
load(fname_ref);
Efar_ref = farfieldpolar3d(mname,p0,res,res,1,1,1,1,1); # the last "+l" is to make sure forward projection
Efar = Efar - Efar_ref;
E2_far = sum(abs(Efar)^2,3);
ux = farfieldux(mname,p0,res,res);
uy = farfielduy(mname,p0,res,res);
image(ux,uy,E2_far,"","","|E|^2 in far field (a.u.)","polar");

load(fname);