############################################
# file: PSL_Cu_scattering_helper.lsf
# Reproduces the results in Kim et al. 
# figures number 6 to 9.
# 
# Copyright 2008 Lumerical Solutions Inc
############################################


#############################################
# run simulations 
run_sim=fileexists(file+"ref.fsp"); 
if (run_sim==0) {  # run sim if it doesn't already exist
  switchtolayout;
  setnamed("source1","center wavelength",lambda);
  select("sphere");          
  delete;                    
  save(file+"ref.fsp");
  run;
}

run_sim=fileexists(file+"def.fsp"); 
if (run_sim==0) {  # run sim if it doesn't already exist
  switchtolayout;
  addsphere;
  set("x",0);
  set("y",0);
  set("z",rad +0* 1.5e-9);
  set("radius",rad);
  if (mat=="Cu") {
    set("material","Cu (Copper) - Palik");
  } else {
    set("material",mat);
  }
  save(file+"def.fsp");
  run;
}


#############################################
# Load files and get raw data

load(file+"ref"); # Reference sim
# near field data            
Ex_ref=getdata(mname,"Ex");
Ey_ref=getdata(mname,"Ey");
Ez_ref=getdata(mname,"Ez");
#far field data
E_far_ref=farfieldpolar3d(mname,1,resX,resY);

load(file+"def"); # Defect sim
# near field data
Ex_def=getdata(mname,"Ex");
Ey_def=getdata(mname,"Ey");
Ez_def=getdata(mname,"Ez");
x=getdata(mname,"x");
y=getdata(mname,"y");
# far field data
E_far_def = farfieldpolar3d(mname,1,resX,resY);
ux = farfieldux(mname,1,resX,resY);
uy = farfielduy(mname,1,resX,resY);


#############################################
# Near field 
Ex_scat=Ex_def-Ex_ref;
Ey_scat=Ey_def-Ey_ref;
Ez_scat=Ez_def-Ez_ref;
E2_scat=abs(Ex_scat)^2+abs(Ey_scat)^2+abs(Ez_scat)^2;


#############################################
# far field 
E_far_scat = E_far_def - E_far_ref;
Er_far_scat= pinch(E_far_scat,3,1);
Et_far_scat= pinch(E_far_scat,3,2);
Ep_far_scat= pinch(E_far_scat,3,3);
E2_far_scat = abs(Er_far_scat)^2 + abs(Et_far_scat)^2 +abs(Ep_far_scat)^2;
	
costheta=sqrt(1-ux^2-0^2);
theta=acos(costheta);
theta=theta* (2*(ux>=0)-1);
theta=theta*180/pi;

###############################################
# Differential scattering cross section dSigma/dOmega

# E2far is proportional to dSigma/dOmega
dS_dO=E2_far_scat;


###############################################
# Degree of circular polarization

# decompose Etheta, Ephi to E+,E- circular polarization
Ecpos=sqrt(2)/2*(Et_far_scat-1i*Ep_far_scat);
Ecneg=sqrt(2)/2*(Et_far_scat+1i*Ep_far_scat);

# calculate Pc
Icpos=abs(Ecpos)^2;
Icneg=abs(Ecneg)^2;
Pc=(Icpos-Icneg)/(Icpos+Icneg);


###############################################
# Polarization angle (eta)

temp=matrix(resX);  # polarization angle
kappa = linspace(0,359,360)*pi/180;

for (i=1:resX) {
  ES = real(exp(1i*kappa)*Ep_far_scat(i));
  EP = real(exp(1i*kappa)*Et_far_scat(i));  

  diameter = sqrt( (ES)^2+(EP)^2 );
  major_axis = max( diameter);
  pos_major = find(diameter,major_axis);

  major_angle = atan2(ES(pos_major),EP(pos_major))*180/pi;
  if(major_angle < -90) { major_angle = major_angle+180; }
  if(major_angle > 90) { major_angle = major_angle-180; }

  temp(i)=major_angle;
}

# unwrap and shift by 90 degrees to match
# definition in reference
eta=180/pi/2*unwrap(temp*2*pi/180)+90;