/*
outputNear2d.c

plot and log of near2d (2D+3D)
*/

#include "oth.h"
#include "complex.h"
#include "oth_prototype.h"
#include "ev.h"

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

void outputNear2d(void)
{
	if (!NFrequency || !NNear2d || (!Near2dDim[0] && !Near2dDim[1])) return;

	const double eps = 1e-10;
	const int ngline = Near2dObj ? NGline : 0;
	const char fmt[] = "%5d%5d%12.3e%12.3e%12.3e%12.4e%12.4e%9.3f%12.4e%9.3f%12.4e%9.3f\n";

	// log
	FILE *fp;
	if ((fp = fopen(FN_near2d, "w")) == NULL) {
		printf("*** %s open error.\n", FN_near2d);
		return;
	}

	for (int n = 0; n < NNear2d; n++) {
		int    start1 = 0, stop1 = 0;
		int    start2 = 0, stop2 = 0;
		double pos0 = 0, *pos1 = NULL, *pos2 = NULL;
		if      (Near2d[n].dir == 'X') {
			// Y-Z
			start1 = !Near2dIzoom ?  0 : nearest(Near2dHzoom[0], 0, Ny, Yn);
			stop1  = !Near2dIzoom ? Ny : nearest(Near2dHzoom[1], 0, Ny, Yn);
			start2 = !Near2dIzoom ?  0 : nearest(Near2dVzoom[0], 0, Nz, Zn);
			stop2  = !Near2dIzoom ? Nz : nearest(Near2dVzoom[1], 0, Nz, Zn);
			pos0 = Xn[Near2d[n].id0];
			pos1 = &Yn[start1];
			pos2 = &Zn[start2];
		}
		else if (Near2d[n].dir == 'Y') {
			// X-Z
			start1 = !Near2dIzoom ?  0 : nearest(Near2dHzoom[0], 0, Nx, Xn);
			stop1  = !Near2dIzoom ? Nx : nearest(Near2dHzoom[1], 0, Nx, Xn);
			start2 = !Near2dIzoom ?  0 : nearest(Near2dVzoom[0], 0, Nz, Zn);
			stop2  = !Near2dIzoom ? Nz : nearest(Near2dVzoom[1], 0, Nz, Zn);
			pos0 = Yn[Near2d[n].id0];
			pos1 = &Xn[start1];
			pos2 = &Zn[start2];
		}
		else if (Near2d[n].dir == 'Z') {
			// X-Y
			start1 = !Near2dIzoom ?  0 : nearest(Near2dHzoom[0], 0, Nx, Xn);
			stop1  = !Near2dIzoom ? Nx : nearest(Near2dHzoom[1], 0, Nx, Xn);
			start2 = !Near2dIzoom ?  0 : nearest(Near2dVzoom[0], 0, Ny, Yn);
			stop2  = !Near2dIzoom ? Ny : nearest(Near2dVzoom[1], 0, Ny, Yn);
			pos0 = Zn[Near2d[n].id0];
			pos1 = &Xn[start1];
			pos2 = &Yn[start2];
		}
		const int div1 = stop1 - start1;
		const int div2 = stop2 - start2;

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

		for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
			// log
			const int cmp = toupper(Near2d[n].cmp[0]);
			fprintf(fp, "#%d : frequency[Hz] = %.3e\n", (n + 1), Frequency[ifreq]);
			fprintf(fp, "%s", "  No.  No.     X[m]        Y[m]        Z[m]       ");
			fprintf(fp, "%s\n",
				(cmp == 'E') ? "E[V/m]      Ex[V/m]   Ex[deg]    Ey[V/m]   Ey[deg]    Ez[V/m]   Ez[deg]" :
				(cmp == 'H') ? "H[A/m]      Hx[A/m]   Hx[deg]    Hy[A/m]   Hy[deg]    Hz[A/m]   Hz[deg]" :
				(cmp == 'I') ? "I[A]        Ix[A]     Ix[deg]    Iy[A]     Iy[deg]    Iz[A]     Iz[deg]" : "");

			// data
			for (int n1 = 0; n1 <= div1; n1++) {
			for (int n2 = 0; n2 <= div2; n2++) {
				// field
				d_complex_t f[3];
				double pos[3];
				field(n, ifreq, n1, n2, start1, start2, f, pos);
				fcomp(Near2d[n].cmp, f, &amp[n1][n2], &deg[n1][n2]);

				// log
				fprintf(fp, fmt,
					start1 + n1, start2 + n2, pos[0], pos[1], pos[2],
					sqrt(d_norm(f[0]) + d_norm(f[1]) + d_norm(f[2])),
					d_abs(f[0]), d_deg(f[0]),
					d_abs(f[1]), d_deg(f[1]),
					d_abs(f[2]), d_deg(f[2]));
			}
			}

			// 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], EPS2));
				}
				}
			}

			// 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, pos0, Near2d[n].cmp,
						ngline, Gline,
						Width2d, Height2d, Fontsize2d);
					// phase
					if (strcmp(Near2d[n].cmp, "E") && strcmp(Near2d[n].cmp, "H") && strcmp(Near2d[n].cmp, "I")) {
						plot2dNear2d0(
							div1, div2, deg, pos1, pos2,
							1, 0, 1, -180, +180, Near2dContour,
							Title, Frequency[ifreq],
							Near2d[n].dir, 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++) {
								d_complex_t f[3];
								double pos[3], wf;
								field(n, ifreq, n1, n2, start1, start2, f, pos);
								const double wt = 2 * PI * frame / Near2dFrame;
								fcomp2(wt, Near2d[n].cmp, f, &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++) {
							d_complex_t f[3];
							double pos[3];
							field(n, ifreq, n1, n2, start1, start2, f, pos);
							const double wt = 2 * PI * frame / Near2dFrame;
							fcomp2(wt, Near2d[n].cmp, f, &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, pos0, Near2d[n].cmp,
					ngline, Gline);
				ev3d_setColor(160, 160, 160);
				ev3d_drawBox(Xn[0], Yn[0], Zn[0], Xn[Nx], Yn[Ny], Zn[Nz]);
				// phase
				if (strcmp(Near2d[n].cmp, "E") && strcmp(Near2d[n].cmp, "H") && strcmp(Near2d[n].cmp, "I")) {
					plot3dNear2d0(
						div1, div2, deg, pos1, pos2,
						1, 0, 1, -180, +180, Near2dContour,
						Title, Frequency[ifreq],
						Near2d[n].dir, pos0, Near2d[n].cmp,
						ngline, Gline);
					ev3d_setColor(160, 160, 160);
					ev3d_drawBox(Xn[0], Yn[0], Zn[0], Xn[Nx], Yn[Ny], Zn[Nz]);
				}
			}
		}

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

	fclose(fp);
}

// field (complex)
static void field(int n, int ifreq, int n1, int n2, int start1, int start2, d_complex_t *f, double *pos)
{
	// i, j, k
	int i = 0, j = 0, k = 0;
	if      (Near2d[n].dir == 'X') {
		i = Near2d[n].id0;
		j = start1 + n1;
		k = start2 + n2;
	}
	else if (Near2d[n].dir == 'Y') {
		i = start1 + n1;
		j = Near2d[n].id0;
		k = start2 + n2;
	}
	else if (Near2d[n].dir == 'Z') {
		i = start1 + n1;
		j = start2 + n2;
		k = Near2d[n].id0;
	}

	// field
	if      (Near2d[n].cmp[0] == 'E') {
		enode(ifreq, i, j, k, Near2dExclude, f);
	}
	else if (Near2d[n].cmp[0] == 'H') {
		hnode(ifreq, i, j, k, Near2dExclude, f);
		for (int m = 0; m < 3; m++) {
			f[m] = d_rmul(1 / ETA0, f[m]);  // to A/m
		}
	}
	else if (Near2d[n].cmp[0] == 'I') {
		cnode(ifreq, i, j, k, f);
	}

	// position
	pos[0] = Xn[i];
	pos[1] = Yn[j];
	pos[2] = Zn[k];
}

// amplitude and phase
static void fcomp(const char comp[], const d_complex_t f[], double *amp, double *deg)
{
	const int c = toupper(comp[1]);

	if      (c == 'X') {
		*amp = d_abs(f[0]);
		*deg = d_deg(f[0]);
	}
	else if (c == 'Y') {
		*amp = d_abs(f[1]);
		*deg = d_deg(f[1]);
	}
	else if (c == 'Z') {
		*amp = d_abs(f[2]);
		*deg = d_deg(f[2]);
	}
	else {
		*amp = sqrt(d_norm(f[0])
		          + d_norm(f[1])
		          + d_norm(f[2]));
		*deg = 0;
	}
}

// waveform
static void fcomp2(double wt, const char comp[], const d_complex_t f[], double *mag)
{
	const int c = toupper(comp[1]);

	if      (c == 'X') {
		*mag = d_abs(f[0]) * cos(d_rad(f[0]) + wt);
	}
	else if (c == 'Y') {
		*mag = d_abs(f[1]) * cos(d_rad(f[1]) + wt);
	}
	else if (c == 'Z') {
		*mag = d_abs(f[2]) * cos(d_rad(f[2]) + wt);
	}
	else {
		*mag = sqrt(d_norm(f[0]) * cos(d_rad(f[0]) + wt) * cos(d_rad(f[0]) + wt)
		          + d_norm(f[1]) * cos(d_rad(f[1]) + wt) * cos(d_rad(f[1]) + wt)
		          + d_norm(f[2]) * cos(d_rad(f[2]) + wt) * cos(d_rad(f[2]) + wt));
	}
}
