#Compare the FDTD Solution results with theory 

#From FDTD Solutions, get the value of Ez 
Ez = getdata("monitor1","Ez");

#From theory, get the value of Ez 

z = getdata("monitor1","z")+ getdata("monitor1","delta_z"); #second term adjusts z to Yee Cell Ez positions
dz = z(2)-z(1);
mu = getdata("source1","moment");
f = getdata("monitor1","f");
k = 2*pi*f/c;

R = abs(z-dz/2)+1e-10; # to regularize
theta = (z-dz/2<0)*pi;

Ezt = mu/(4*pi*eps0)*exp(1i*k*R)/R*
       (
          k^2*sin(theta)^2+
         (1/R^2-1i*k/R)*(3*cos(theta)^2-1)
       );
#At dipole (z = 0), correct Re(Ezt) to give the theoretical electric field averaged over one mesh cell
V = dz^3; # Mesh cell volume
Ezt(find(z,dz/2)) = -mu/(3*eps0)/V + 1i*imag(Ezt(find(z,dz/2))); 
          
#plot Ez 

plot(z*1e6,real(Ez),real(Ezt),"z (microns)","real(E)","real(E)");
legend("FDTD","theory");
plot(z*1e6,imag(Ez),imag(Ezt),"z (microns)","imag(E)","imag(E)");
legend("FDTD","theory");