/*
plot2d3dNear2d.c

plot near2d field (2D/3D)
*/

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

static void fcomp(const char [], const d_complex_t [], const d_complex_t [], double *, double *);
static void fcomp2(double, const char [], const d_complex_t [], const d_complex_t [], double *);


void plot2d3dNear2d(d_complex_t *****e, d_complex_t *****h)
{
	if (!NFrequency || !NNear2d || (!Near2dDim[0] && !Near2dDim[1])) return;

	const double eps = 1e-10;
	const int ngline = Near2dObj ? NGline : 0;

	for (int n = 0; n < NNear2d; n++) {

		int div1 = Near2d[n].div1;
		int div2 = Near2d[n].div2;

		// alloc
		double **amp = (double **)malloc((div1 + 1) * sizeof(double *));
		double **deg = (double **)malloc((div1 + 1) * sizeof(double *));
		for (int i = 0; i <= div1; i++) {
			amp[i] = (double *)malloc((div2 + 1) * sizeof(double));
			deg[i] = (double *)malloc((div2 + 1) * sizeof(double));
		}
		double *pos1 = (double *)malloc((div1 + 1) * sizeof(double));
		double *pos2 = (double *)malloc((div2 + 1) * sizeof(double));

		// position
		for (int n1 = 0; n1 <= div1; n1++) {
			pos1[n1] = Near2d[n].pos1[0] + ((double)n1 / div1) * (Near2d[n].pos1[1] - Near2d[n].pos1[0]);
		}
		for (int n2 = 0; n2 <= div2; n2++) {
			pos2[n2] = Near2d[n].pos2[0] + ((double)n2 / div2) * (Near2d[n].pos2[1] - Near2d[n].pos2[0]);
		}

		for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
			// E/H
			for (int n1 = 0; n1 <= div1; n1++) {
			for (int n2 = 0; n2 <= div2; n2++) {
				fcomp(Near2d[n].cmp, e[n][ifreq][n1][n2], h[n][ifreq][n1][n2], &amp[n1][n2], &deg[n1][n2]);
			}
			}

			// to dB
			if (Near2dScale.db) {
				for (int n1 = 0; n1 <= div1; n1++) {
				for (int n2 = 0; n2 <= div2; n2++) {
					amp[n1][n2] = 20 * log10(MAX(amp[n1][n2], eps));
				}
				}
			}

			// plot (2D)
			if (Near2dDim[0]) {
				// harmonic
				if (!Near2dFrame) {
					plot2dNear2d0(
						div1, div2, amp, pos1, pos2,
						0, Near2dScale.db, Near2dScale.user, Near2dScale.min, Near2dScale.max, Near2dContour,
						Title, Frequency[ifreq],
						Near2d[n].dir, Near2d[n].pos0, Near2d[n].cmp,
						ngline, Gline,
						Width2d, Height2d, Fontsize2d);
					// phase
					if (strcmp(Near2d[n].cmp, "E") && strcmp(Near2d[n].cmp, "H")) {
						plot2dNear2d0(
							div1, div2, deg, pos1, pos2,
							1, 0, 1, -180, +180, Near2dContour,
							Title, Frequency[ifreq],
							Near2d[n].dir, Near2d[n].pos0, Near2d[n].cmp,
							ngline, Gline,
							Width2d, Height2d, Fontsize2d);
					}
				}
				// animation
				else if ((n == 0) && (ifreq == 0)) {
					// alloc
					double **mag = (double **)malloc((div1 + 1) * sizeof(double *));
					for (int n1 = 0; n1 <= div1; n1++) {
						mag[n1] = (double *)malloc((div2 + 1) * sizeof(double));
					}
					// auto scale -> user scale
					scale_t scale = Near2dScale;
					if (scale.user == 0) {
						// data : max
						double dmax = 0;
						for (int frame = 0; frame < Near2dFrame; frame++) {
							for (int n1 = 0; n1 <= div1; n1++) {
							for (int n2 = 0; n2 <= div2; n2++) {
								double wf;
								const double wt = 2 * PI * frame / Near2dFrame;
								fcomp2(wt, Near2d[n].cmp, e[n][ifreq][n1][n2], h[n][ifreq][n1][n2], &wf);
								dmax = MAX(dmax, wf);
							}
							}
						}
						scale.user = 1;
						if (scale.db) {
							// dB
							scale.max = 20 * log10(MAX(dmax, eps));
							scale.min = scale.max - 40;
						}
						else {
							// linear
							scale.max = dmax;
							scale.min = 0;
						}
					}
					// frame
					for (int frame = 0; frame < Near2dFrame; frame++) {
						// data
						for (int n1 = 0; n1 <= div1; n1++) {
						for (int n2 = 0; n2 <= div2; n2++) {
							const double wt = 2 * PI * frame / Near2dFrame;
							fcomp2(wt, Near2d[n].cmp, e[n][ifreq][n1][n2], h[n][ifreq][n1][n2], &mag[n1][n2]);
						}
						}
						// to dB
						if (scale.db) {
							for (int n1 = 0; n1 <= div1; n1++) {
							for (int n2 = 0; n2 <= div2; n2++) {
								mag[n1][n2] = 20 * log10(MAX(mag[n1][n2], eps));
							}
							}
						}
						// plot 2D
						plot2dNear2d0(
							div1, div2, mag, pos1, pos2,
							0, scale.db, scale.user, scale.min, scale.max, Near2dContour,
							Title, Frequency[ifreq],
							Near2d[n].dir, Near2d[n].pos0, Near2d[n].cmp,
							ngline, Gline,
							Width2d, Height2d, Fontsize2d);
					}
					// free
					for (int n1 = 0; n1 <= div1; n1++) {
						free(mag[n1]);
					}
					free(mag);
				}
			}

			// plot (3D)
			if (Near2dDim[1]) {
				plot3dNear2d0(
					div1, div2, amp, pos1, pos2,
					0, Near2dScale.db, Near2dScale.user, Near2dScale.min, Near2dScale.max, Near2dContour,
					Title, Frequency[ifreq],
					Near2d[n].dir, Near2d[n].pos0, Near2d[n].cmp,
					ngline, Gline);
				// phase
				if (strcmp(Near2d[n].cmp, "E") && strcmp(Near2d[n].cmp, "H")) {
					plot3dNear2d0(
						div1, div2, deg, pos1, pos2,
						1, 0, 1, -180, +180, Near2dContour,
						Title, Frequency[ifreq],
						Near2d[n].dir, Near2d[n].pos0, Near2d[n].cmp,
						ngline, Gline);
				}
			}
		}

		// free
		for (int n1 = 0; n1 <= div1; n1++) {
			free(amp[n1]);
			free(deg[n1]);
		}
		free(amp);
		free(deg);
		free(pos1);
		free(pos2);
	}
}


// amplitude and phase
static void fcomp(const char comp[], const d_complex_t e[], const d_complex_t h[], double *amp, double *deg)
{
	if      (!strcmp(comp, "E")) {
		*amp = sqrt(d_norm(e[0])
		          + d_norm(e[1])
		          + d_norm(e[2]));
		*deg = 0;
	}
	else if (!strcmp(comp, "Ex")) {
		*amp = d_abs(e[0]);
		*deg = d_deg(e[0]);
	}
	else if (!strcmp(comp, "Ey")) {
		*amp = d_abs(e[1]);
		*deg = d_deg(e[1]);
	}
	else if (!strcmp(comp, "Ez")) {
		*amp = d_abs(e[2]);
		*deg = d_deg(e[2]);
	}
	else if (!strcmp(comp, "H")) {
		*amp = sqrt(d_norm(h[0])
		          + d_norm(h[1])
		          + d_norm(h[2]));
		*deg = 0;
	}
	else if (!strcmp(comp, "Hx")) {
		*amp = d_abs(h[0]);
		*deg = d_deg(h[0]);
	}
	else if (!strcmp(comp, "Hy")) {
		*amp = d_abs(h[1]);
		*deg = d_deg(h[1]);
	}
	else if (!strcmp(comp, "Hz")) {
		*amp = d_abs(h[2]);
		*deg = d_deg(h[2]);
	}
}


// waveform
static void fcomp2(double wt, const char comp[], const d_complex_t e[], const d_complex_t h[], double *mag)
{
	if      (!strcmp(comp, "E")) {
		*mag = sqrt(d_norm(e[0]) * cos(d_rad(e[0]) + wt) * cos(d_rad(e[0]) + wt)
		          + d_norm(e[1]) * cos(d_rad(e[1]) + wt) * cos(d_rad(e[1]) + wt)
		          + d_norm(e[2]) * cos(d_rad(e[2]) + wt) * cos(d_rad(e[2]) + wt));
	}
	else if (!strcmp(comp, "Ex")) {
		*mag = d_abs(e[0]) * cos(d_rad(e[0]) + wt);
	}
	else if (!strcmp(comp, "Ey")) {
		*mag = d_abs(e[1]) * cos(d_rad(e[1]) + wt);
	}
	else if (!strcmp(comp, "Ez")) {
		*mag = d_abs(e[2]) * cos(d_rad(e[2]) + wt);
	}
	else if (!strcmp(comp, "H")) {
		*mag = sqrt(d_norm(h[0]) * cos(d_rad(h[0]) + wt) * cos(d_rad(h[0]) + wt)
		          + d_norm(h[1]) * cos(d_rad(h[1]) + wt) * cos(d_rad(h[1]) + wt)
		          + d_norm(h[2]) * cos(d_rad(h[2]) + wt) * cos(d_rad(h[2]) + wt));
	}
	else if (!strcmp(comp, "Hx")) {
		*mag = d_abs(h[0]) * cos(d_rad(h[0]) + wt);
	}
	else if (!strcmp(comp, "Hy")) {
		*mag = d_abs(h[1]) * cos(d_rad(h[1]) + wt);
	}
	else if (!strcmp(comp, "Hz")) {
		*mag = d_abs(h[2]) * cos(d_rad(h[2]) + wt);
	}
}
