/*
plot3dFar2d.c

plot far2d field (3D)
*/

#include "omm.h"
#include "complex.h"
#include "omm_prototype.h"

static void statistics(int, d_complex_t **, d_complex_t **, char [], char []);
static double opterror(int);

void plot3dFar2d(d_complex_t ***etheta, d_complex_t ***ephi)
{
	const char scomp[][BUFSIZ] = {"E-abs", "E-theta", "E-phi", "E-major", "E-minor", "E-RHCP", "E-LHCP"};
	char unit[BUFSIZ], stat1[BUFSIZ], stat2[BUFSIZ];
	double e[7];

	// unit
	if (Far2dScale.db) {
		sprintf(unit, "[dB%s]", (IPlanewave ? "sm" : ""));
	}
	else {
		strcpy(unit, (IPlanewave ? "[m*m]" : ""));
	}

	// alloc
	double **ef = (double **)malloc((Far2d.divtheta + 1) * sizeof(double *));
	for (int itheta = 0; itheta <= Far2d.divtheta; itheta++) {
		ef[itheta] = (double *)malloc((Far2d.divphi + 1) * sizeof(double));
	}
	char *comment[5];
	for (int n = 0; n < 5; n++) {
		comment[n] = (char *)malloc(BUFSIZ * sizeof(char));
	}

	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		// calculate statistics
		statistics(ifreq, etheta[ifreq], ephi[ifreq], stat1, stat2);

		for (int icomp = 0; icomp < 7; icomp++) {
			if (Far2dComp[icomp]) {
				double emax = 0;
				double tmax = 0;
				double pmax = 0;
				for (int itheta = 0; itheta <= Far2d.divtheta; itheta++) {
				for (int iphi   = 0; iphi   <= Far2d.divphi;   iphi++  ) {
					farComponent(etheta[ifreq][itheta][iphi], ephi[ifreq][itheta][iphi], e);
					ef[itheta][iphi] = e[icomp];
					// max
					if (e[icomp] > emax) {
						emax = e[icomp];
						tmax = 180.0 * itheta / Far2d.divtheta;
						pmax = 360.0 * iphi   / Far2d.divphi;
					}
				}
				}

				// comment
				sprintf(comment[0], "%s%s f=%.3e[Hz]", scomp[icomp], unit, Frequency[ifreq]);
				strcpy(comment[1], stat1);
				strcpy(comment[2], stat2);
				sprintf(comment[3], "max[deg] @ theta=%.1f, phi=%.1f", tmax, pmax);
				strcpy(comment[4], Title);

				// plot
				plot3dFar2d0(
					Far2d.divtheta, Far2d.divphi, ef,
					Far2dScale.db, Far2dScale.user, Far2dScale.min, Far2dScale.max,
					NGline, Gline, Far2dObj,
					5, comment);
			}
		}
	}

	// free
	for (int itheta = 0; itheta <= Far2d.divtheta; itheta++) {
		free(ef[itheta]);
	}
	free(ef);
	for (int n = 0; n < 5; n++) {
		free(comment[n]);
	}
}

static void statistics(int ifreq, d_complex_t **etheta, d_complex_t **ephi, char stat1[], char stat2[])
{
	const double eps = 1e-20;

	double pmax = 0;
	double psum = 0;
	const double sumfactor = (PI / Far2d.divtheta) * (2 * PI / Far2d.divphi) / (4 * PI);
	for (int itheta = 0; itheta <= Far2d.divtheta; itheta++) {
	for (int iphi   = 0; iphi   <  Far2d.divphi;   iphi++  ) {
		const double pow = d_norm(etheta[itheta][iphi])
		                 + d_norm(  ephi[itheta][iphi]);
		psum += sumfactor * sin(PI * itheta / Far2d.divtheta) * pow;
		pmax = MAX(pow, pmax);
	}
	}

	if (NFeed) {
		pmax /= psum;
		if (Far2dScale.db) {
			const double pmaxdb = 10 * log10(MAX(pmax, eps));
			sprintf(stat1, "directive gain = %.3f[dBi]", pmaxdb);
		}
		else {
			sprintf(stat1, "directive gain = %.3f", pmax);
		}
		sprintf(stat2, "efficiency = %.3f[%%]", psum * 100);
	}
	else if (IPlanewave) {
		if (Far2dScale.db) {
			const double psumdb = 10 * log10(MAX(psum, eps));
			sprintf(stat1, "total cross section = %.3f[dBsm]", psumdb);
		}
		else {
			sprintf(stat1, "total cross section = %.3e[m*m]", psum);
		}
		double loss = 0;
		for (int ie = 0; ie < NElement; ie++) {
			if (Element[ie].iload == 1) {
				loss += 0.5 * Element[ie].load * d_norm(Cv[(ifreq * NElement) + ie]);
			}
		}
		const double ei = 1;
		const double rhs = opterror(ifreq);
		const double lhs = psum + (2 * ETA0 * loss / (ei * ei));
		//printf("%e %e\n", rhs, lhs);
		const double oerr = fabs(1 - (rhs / lhs));
		sprintf(stat2, "optical theorem error = %.3f[%%]", oerr * 100);
	}
}

static double opterror(int ifreq)
{
	assert(IPlanewave == 1);
	assert((ifreq >= 0) && (ifreq < NFrequency));

	const double kwave = (2 * PI * Frequency[ifreq]) / C;
	const double fctr = farfactor(ifreq);

	d_complex_t etheta, ephi;
	farfield(ifreq, 180 - Planewave.theta, Planewave.phi + 180, fctr, &etheta, &ephi);

	double polfctr = 1;
	if      (Planewave.pol == 1) {
		polfctr = etheta.i;
	}
	else if (Planewave.pol == 2) {
		polfctr = ephi.i;
	}
	else if (Planewave.pol == 3) {
		polfctr = sqrt(0.5) * (-etheta.i - ephi.r);
	}
	else if (Planewave.pol == 4) {
		polfctr = sqrt(0.5) * (-etheta.i + ephi.r);
	}
	else if (Planewave.pol == 5) {
		const d_complex_t j = d_complex(0, 1);
		const double dtor = atan(1) / 45;
		const double cosa = cos(Planewave.a * dtor);
		const double sina = sin(Planewave.a * dtor);
		const double r = Planewave.r;
		const d_complex_t emajor = d_mul(j, d_add(d_rmul(-cosa, etheta), d_rmul(+sina, ephi)));
		const d_complex_t eminor = d_mul(j, d_add(d_rmul(-sina, etheta), d_rmul(-cosa, ephi)));
		polfctr = -d_add(emajor, d_mul(d_complex(0, r), eminor)).r / sqrt(1 + (r * r));
	}

	return (sqrt(4 * PI) / kwave) * polfctr;
}
