/*
outputPloss.c

output total Power Loss
*/

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


void calcPloss(void)
{
	const int exclude = 0;

	// alloc
	PLoss = (double *)malloc(NFrequency * sizeof(double));

	// calculation
	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		int64_t m0 = ifreq * NMaterial;
		int i;

		// check loss ?
		double sum = 0;
#ifdef _OPENMP
#pragma omp parallel for reduction (+: sum)
#endif
		for (    i = 0; i <  Nx; i++) {
		for (int j = 0; j <  Ny; j++) {
		for (int k = 0; k <  Nz; k++) {
			sum += Material[MEX(i, j, k) + m0].cepsr.i
			     + Material[MEY(i, j, k) + m0].cepsr.i
			     + Material[MEZ(i, j, k) + m0].cepsr.i
			     + Material[MHX(i, j, k) + m0].cmur.i
			     + Material[MHY(i, j, k) + m0].cmur.i
			     + Material[MHZ(i, j, k) + m0].cmur.i;
		}
		}
		}
		if (fabs(sum) < EPS2) {
			PLoss[ifreq] = 0;
			continue;
		}

		// Ex
		double exsum = 0;
#ifdef _OPENMP
#pragma omp parallel for reduction (+: exsum)
#endif
		for (    i = 0; i <  Nx; i++) {
		for (int j = 0; j <= Ny; j++) {
		for (int k = 0; k <= Nz; k++) {
			exsum += d_norm(exyee(ifreq, i, j, k, exclude))
			       * Material[MEX(i, j, k) + m0].cepsr.i
			       * (DXc[i] * DYn[j] * DZn[k]);
		}
		}
		}
		// Ey
		double eysum = 0;
#ifdef _OPENMP
#pragma omp parallel for reduction (+: eysum)
#endif
		for (    i = 0; i <= Nx; i++) {
		for (int j = 0; j <  Ny; j++) {
		for (int k = 0; k <= Nz; k++) {
			eysum += d_norm(eyyee(ifreq, i, j, k, exclude))
			       * Material[MEY(i, j, k) + m0].cepsr.i
			       * (DXn[i] * DYc[j] * DZn[k]);
		}
		}
		}
		// Ez
		double ezsum = 0;
#ifdef _OPENMP
#pragma omp parallel for reduction (+: ezsum)
#endif
		for (    i = 0; i <= Nx; i++) {
		for (int j = 0; j <= Ny; j++) {
		for (int k = 0; k <  Nz; k++) {
			ezsum += d_norm(ezyee(ifreq, i, j, k, exclude))
			       * Material[MEZ(i, j, k) + m0].cepsr.i
			       * (DXn[i] * DYn[j] * DZc[k]);
		}
		}
		}
		// Hx
		double hxsum = 0;
#ifdef _OPENMP
#pragma omp parallel for reduction (+: hxsum)
#endif
		for (    i = 0; i <= Nx; i++) {
		for (int j = 0; j <  Ny; j++) {
		for (int k = 0; k <  Nz; k++) {
			hxsum += d_norm(hxyee(ifreq, i, j, k, exclude))
			       * Material[MHX(i, j, k) + m0].cmur.i
			       * (DXn[i] * DYc[j] * DZc[k]);
		}
		}
		}
		// Hy
		double hysum = 0;
#ifdef _OPENMP
#pragma omp parallel for reduction (+: hysum)
#endif
		for (    i = 0; i <  Nx; i++) {
		for (int j = 0; j <= Ny; j++) {
		for (int k = 0; k <  Nz; k++) {
			hysum += d_norm(hyyee(ifreq, i, j, k, exclude))
			       * Material[MHY(i, j, k) + m0].cmur.i
			       * (DXc[i] * DYn[j] * DZc[k]);
		}
		}
		}
		// Hz
		double hzsum = 0;
#ifdef _OPENMP
#pragma omp parallel for reduction (+: hzsum)
#endif
		for (    i = 0; i <  Nx; i++) {
		for (int j = 0; j <  Ny; j++) {
		for (int k = 0; k <= Nz; k++) {
			hzsum += d_norm(hzyee(ifreq, i, j, k, exclude))
			       * Material[MHZ(i, j, k) + m0].cmur.i
			       * (DXc[i] * DYc[j] * DZn[k]);
		}
		}
		}

		// power loss
		//printf("%d %e %e %e %e %e %e\n", ifreq, exsum, eysum, ezsum, hxsum, hysum, hzsum);
		double omega = 2 * PI * Frequency[ifreq];
		PLoss[ifreq] = 0.5 * omega * EPS0 * fabs(exsum + eysum + ezsum + hxsum + hysum + hzsum);
	}
}

static void _outputPloss(FILE *fp)
{
	fprintf(fp, "=== power loss ===\n");
	fprintf(fp, "%s\n", "  freq[Hz]     Loss[dBW]");
	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		const double p = 10 * log10(MAX(PLoss[ifreq], 1e-30));
		fprintf(fp, "%11.3e%12.3f\n", Frequency[ifreq], p);
	}
	fflush(fp);
}


void outputPloss(FILE *fp)
{
	// sum
	double sum = 0;
	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		//printf("%d %e\n", ifreq, PLoss[ifreq]);
		sum += PLoss[ifreq];
	}

	// output (if exists)
	if (sum > 1e-30) {
		_outputPloss(stdout);
		_outputPloss(fp);
	}
}
