#Compare the FDTD Solution results with theory 

#From FDTD Solutions, get the value of Ez ( and eventually green's function)
Ez = getdata("monitor1","Ez");

#From theory, get the value of Ez ( and eventually green's function)

z = getdata("monitor1","z");
mu = getdata("source1","moment");
f = getdata("monitor1","f");
k = 2*pi*f/c;

R = abs(z)+1e-10; # to regularize
theta = (z<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
dz = z(2)-z(1);
V = dz^3; # Mesh cell volume
Ezt(find(z,0)) = -mu/(3*eps0)/V + 1i*imag(Ezt(find(z,0))); 
          
#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");