# run the simulation
run;

# get Ex and Ey from monitor Z
Ex = getdata("Z","Ex");
Ey = getdata("Z","Ey");
z = getdata("Z","z");
f = getdata("Z","f");

# calculate theory
k = 2*pi*f/c;
phase_x = 0*pi/180;
Ex_theory = exp(1i*(k*z+phase_x));
phase_y = 90*pi/180;
Ey_theory = exp(1i*(k*z+phase_y));

# compare phases between FDTD and theory
plot(z*1e6,unwrap(angle(Ex)),unwrap(angle(Ex_theory)),"z (microns)","phase Ex (radians)","FDTD vs theory");
legend("FDTD","+kz+phi");
plot(z*1e6,unwrap(angle(Ey)),unwrap(angle(Ey_theory)),"z (microns)","phase Ey (radians)","FDTD vs theory");
legend("FDTD","+kz+phi");

# plot phases of Ex and Ey
plot(z*1e6,unwrap(angle(Ex)),unwrap(angle(Ey)),"z (microns)","phase (radians)","Phase of Ex and Ey along z");
legend("phase Ex = +kz","phase Ey = +kz + pi/2");

# verify amplitude
plot(z*1e6,abs(Ex),abs(Ey),"z (microns)","abs(E)","Amplitudes of Ex and Ey");
legend("|Ex|","|Ey|");