run_sims=1;


# Do the 0 (P) degree polarization
switchtolayout;
setnamed("source1","polarization angle",0);
if (run_sims) {
  save("temp_0"); 
  run(1);
} else {
  load("temp_0");
}
T_0=-transmission("field_XY");
m="field_XZ";
x=getdata(m,"x");
z=getdata(m,"z");
E2_0=getelectric(m);


# Do the 90 (S) degree polarization
switchtolayout;
setnamed("source1","polarization angle",90);
if (run_sims) {
  save("temp_90"); 
  run(1);
} else {
  load("temp_90");
}
T_90=-transmission("field_XY");
m="field_XZ";
x=getdata(m,"x");
z=getdata(m,"z");
E2_90=getelectric(m);


# calculate the unpolarized transmission and field profile
T_unpol=0.5*(T_0+T_90);
E2_unpol=0.5*(E2_0+E2_90);



# plot the final results
colorbar_max=.3;
image(x*1e6,z*1e6,E2_0,"x (um)","z (um)","|E|^2 (P polarization)");
setplot("colorbar min",0);  setplot("colorbar max",colorbar_max); 
image(x*1e6,z*1e6,E2_90,"x (um)","z (um)","|E|^2 (S polarization)");
setplot("colorbar min",0);  setplot("colorbar max",colorbar_max); 
image(x*1e6,z*1e6,E2_unpol,"x (um)","z (um)","|E|^2 (unpolarized)");
setplot("colorbar min",0);  setplot("colorbar max",colorbar_max); 

?" Transmission (P polarization):   "+num2str(round(T_0*100))+"%";
?" Transmission (S polarization):   "+num2str(round(T_90*100))+"%";
?" Transmission (unpolarized):      "+num2str(round(T_unpol*100))+"%";