/*
outputZin.c

output input impedance
*/

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


void calcZin(void)
{
	if (!NFeed || !NFrequency) return;

	Zin = (d_complex_t *)malloc(NFeed * NFrequency * sizeof(d_complex_t));
	Ref =      (double *)malloc(NFeed * NFrequency * sizeof(double));

	for (int ifeed = 0; ifeed < NFeed; ifeed++) {
		int i = Feed[ifeed].i;
		int j = Feed[ifeed].j;
		int k = Feed[ifeed].k;
		d_complex_t v = Feed[ifeed].volt;
		d_complex_t z0 = d_complex(Feed[ifeed].z0, 0);
		for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
			d_complex_t c = (Feed[ifeed].dir == 'X') ? cxyee(ifreq, i, j, k)
			              : (Feed[ifeed].dir == 'Y') ? cyyee(ifreq, i, j, k)
			                                         : czyee(ifreq, i, j, k);

			int id = (ifeed * NFrequency) + ifreq;

			// Zin
			Zin[id] = d_div(v, c);
			//printf("%d %d %f %f %f %f %f %f\n", ifeed, ifreq, v.r, v.i, c.r, c.i, Zin[id].r, Zin[id].i);

			// Reflection = (Zin - Z0) / (Zin + Z0)
			d_complex_t ref = d_div(d_sub(Zin[id], z0), d_add(Zin[id], z0));
			Ref[id] = 20 * log10(MAX(d_abs(ref), EPS));
		}
	}
}


static void _outputZin(FILE *fp)
{
	fprintf(fp, "%s\n", "=== input impedance ===");

	for (int ifeed = 0; ifeed < NFeed; ifeed++) {
		fprintf(fp, "feed #%d (Z0[ohm] = %.2f)\n", ifeed + 1, Feed[ifeed].z0);
		fprintf(fp, "  %s\n", "frequency[Hz] Rin[ohm]   Xin[ohm]    Gin[mS]    Bin[mS]    Ref[dB]       VSWR");
		for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
			const int id = (ifeed * NFrequency) + ifreq;
			const d_complex_t yin = d_inv(Zin[id]);
			const double gamma = pow(10, Ref[id] / 20);
			const double vswr = (fabs(1 - gamma) > EPS) ? (1 + gamma) / (1 - gamma) : 1000;
			fprintf(fp, "%13.5e%11.3f%11.3f%11.3f%11.3f%11.3f%11.3f\n",
				Frequency[ifreq], Zin[id].r, Zin[id].i, yin.r * 1e3, yin.i * 1e3, Ref[id], vswr);
		}
	}

	fflush(fp);
}


void outputZin(FILE *fp)
{
	_outputZin(stdout);
	_outputZin(fp);
}
