/*
outputSpara.c

output S-parameters
*/

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

void calcSpara(void)
{
	if (!NPoint || !NFrequency) return;

	d_complex_t *ce = (d_complex_t *)malloc((NPoint + 2) * sizeof(d_complex_t));
	Spara = (d_complex_t *)malloc(NPoint * NFrequency * sizeof(d_complex_t));

	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		for (int n = 0; n < NPoint + 2; n++) {
			const int i = Point[n].i;
			const int j = Point[n].j;
			const int k = Point[n].k;
			const char dir = Point[n].dir;
			if      (dir == 'X') {
				ce[n] = d_cast(Ex[ifreq][i][j][k]);
			}
			else if (dir == 'Y') {
				ce[n] = d_cast(Ey[ifreq][i][j][k]);
			}
			else if (dir == 'Z') {
				ce[n] = d_cast(Ez[ifreq][i][j][k]);
			}
			// plane wave : add Einc
			if (IPlanewave) {
				d_complex_t einc = d_complex(0, 0);
				if      (dir == 'X') {
					einc = planewave(0, Xc[i], Yn[j], Zn[k], ifreq);
				}
				else if (dir == 'Y') {
					einc = planewave(1, Xn[i], Yc[j], Zn[k], ifreq);
				}
				else if (dir == 'Z') {
					einc = planewave(2, Xn[i], Yn[j], Zc[k], ifreq);
				}
				ce[n] = d_add(ce[n], einc);
			}
		}
		// E1, E1+, E1-
		d_complex_t ce0 = ce[0];
		d_complex_t cep = ce[NPoint];
		d_complex_t cem = ce[NPoint + 1];
		//printf("%d %f %f %f %f %f %f\n", ifreq, ce0.r, ce0.i, cep.r, cep.i, cem.r, cem.i);
		// S11
		d_complex_t c1 = d_div(d_add(cep, cem), ce0);
		d_complex_t c2 = d_sqrt(d_sub(d_mul(c1, c1), d_complex(4, 0)));
		d_complex_t c3 = d_add(c1, c2);
		if (c3.i < 0) c3 = d_sub(c1, c2);  // Im > 0
		c3 = d_rmul(0.5, c3);              // exp(+gd)
		d_complex_t c4 = d_div(d_complex(1, 0), c3);   // exp(-gd)
		d_complex_t c5 = d_sub(d_mul(c4, c4), d_mul(c3, c3));  // exp(-2gd) - exp(2gd)
		d_complex_t c6 = d_div(d_sub(d_mul(cep, c4), d_mul(cem, c3)), c5);  // E+
		d_complex_t c7 = d_div(d_sub(d_mul(cem, c4), d_mul(cep, c3)), c5);  // E-
		// S11 = E- / E+
		Spara[ifreq] = d_div(c7, c6);
		// Sn1 (n > 1)
		for (int ipoint = 1; ipoint < NPoint; ipoint++) {
			Spara[(ipoint * NFrequency) + ifreq] = d_div(ce[ipoint], c6);  // Sn1 = En / E+
		}
	}

	free(ce);
}

static void _outputSpara(FILE *fp)
{
	fprintf(fp, "=== S-parameters ===\n");

	fprintf(fp, "  freq[Hz] ");
	for (int ipoint = 0; ipoint < NPoint; ipoint++) {
		fprintf(fp, "  S%d1[dB] S%d1[deg]", ipoint + 1, ipoint + 1);
	}
	fprintf(fp, "\n");

	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		fprintf(fp, "%11.3e", Frequency[ifreq]);
		for (int ipoint = 0; ipoint < NPoint; ipoint++) {
			int id = (ipoint * NFrequency) + ifreq;
			fprintf(fp, "%9.3f%9.3f", 10 * log10(MAX(d_norm(Spara[id]), EPS2)), d_deg(Spara[id]));
		}
		fprintf(fp, "\n");
	}

	fflush(fp);
}

void outputSpara(FILE *fp)
{
	_outputSpara(stdout);
	_outputSpara(fp);
}
