# choose to rerun the simulations (must be 1 the first time)
rerun_simulations = 1;

# choose the weighting factor for the 0 and 90 degree incident polarizations
weight_0 = 1;
weight_90 = exp( 1i*45*pi/180); #45 degrees, elliptical polarization

# run the simulations if requested
if(rerun_simulations) {
    load("usr_polarization2.fsp");
    switchtolayout;
    select("source1");
    set("polarization angle",0);
    save("usr_polarization2_0.fsp");
    run;
    switchtolayout;
    select("source1");
    set("polarization angle",90);
    save("usr_polarization2_90.fsp");
    run;
}

# collect the field data from each simulation
load("usr_polarization2_0.fsp");
    Ex_0 = getdata("XY","Ex");
    Ey_0 = getdata("XY","Ey");
load("usr_polarization2_90.fsp");
    Ex_90 = getdata("XY","Ex");
    Ey_90 = getdata("XY","Ey");
    x = getdata("XY","x");
    y = getdata("XY","y");

# create the total field based on the weighting factors
Ex = weight_0*Ex_0 + weight_90*Ex_90;
Ey = weight_0*Ey_0 + weight_90*Ey_90;

# find the middle points to take cross sections
nx = length(x);
pos_y = round(length(y)/2);

# plot the field amplitudes
plot(x*1e6,abs(Ex(1:nx,pos_y)),abs(Ey(1:nx,pos_y)),"x (microns)","Amplitude of E-field (V/m)","","linewidth=2");
legend("|Ex|","|Ey|");

# plot the phases
plot(x*1e6,180/pi*unwrap(angle(Ey(1:nx,pos_y))-angle(Ex(1:nx,pos_y))),
                   "x (microns)","angle(Ey)-angle(Ex) (degrees)",
                   "","linewidth=2");
legend("");

# reload usr_polarization2.fsp
load("usr_polarization2.fsp");