# This script calculates a compares the
# FDTD far field radar cross section
# with theory.
##########################################


runanalysis;
XY = getresult("scat_ff","XY");
XZ = getresult("scat_ff","XZ");
YZ = getresult("scat_ff","YZ");

phi = XY.phi;
x = cos(phi);
z = sin(phi);
E2_XZ = XZ.E2;
E2_XY = XY.E2;
E2_YZ = YZ.E2;
E2_XZ = E2_XZ/max(E2_XZ);    # set max value to 1


##########################################
# calculate theory.  We only have analytic solution
# for half and full wavelength antenna.
half_wavelength=1;  # 1 for half wavelength, 0 for full wavelength

phi2 = (phi+90)*pi/180; # measured from the z-axis
phi2 = phi2+1e-10; # avoid divide by zero problems

if (half_wavelength) {  # half wavelength
  E2_xz_theory= ( (cos (pi/2*cos(phi2)) )^2 ) /
                ( (sin(phi2)            )^2 );
} else {                # full wavelength
  E2_xz_theory= ( 4*(cos (pi/2*cos(phi2)) )^4 ) /
                ( (sin(phi2)            )^2 );
}

E2_xz_theory=E2_xz_theory/max(E2_xz_theory);   # set max value to 1


##########################################
# plot comparison between FDTD and theory
plot(phi,E2_XZ,E2_xz_theory,
     "phi (deg)","Intensity","XZ radiation cross section");
legend("FDTD","Theory");

polar(phi/180*pi,E2_XZ,E2_xz_theory, 
      "phi (deg)","Intensity","XZ radiation pattern");
legend("FDTD","Theory");