/*
getvalue.c
*/

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

// (static) far field
// e[0-6] = Eabs/Etheta/Ephi/Emajor/Eminor/RHCP/LHCP
static void efar(int ifreq, double theta, double phi, double e[])
{
	assert(ifreq >= 0);
	d_complex_t etheta, ephi;
	double ffctr = farfactor(ifreq);
	farfield(ifreq, theta, phi, ffctr, &etheta, &ephi);
	farComponent(etheta, ephi, e);
}


// input impedance
d_complex_t zin(int ifeed, int ifreq)
{
	assert(ifeed >= 0);
	assert(ifreq >= 0);
	return Zin[(ifeed * NFrequency) + ifreq];
}


// input impedance difference
double zindiff(int ifeed, int ifreq, double z0)
{
	assert(ifeed >= 0);
	assert(ifreq >= 0);
	assert(z0 > 0);
	const d_complex_t z = zin(ifeed, ifreq);
	return d_abs(d_sub(z, d_complex(z0, 0)));
}


// reflection
double reflection(int ifeed, int ifreq, double z0, int db)
{
	assert(ifeed >= 0);
	assert(ifreq >= 0);
	assert(z0 > 0);
	const d_complex_t z = zin(ifeed, ifreq);
	double ret = d_abs(d_div(d_sub(z, d_complex(Z0, 0)), d_add(z, d_complex(Z0, 0))));
	if (db) {
		ret = 20 * log10(MAX(ret, 1e-2));
	}
	return ret;
}


// VSWR
double vswr(int ifeed, int ifreq, double z0)
{
	assert(ifeed >= 0);
	assert(ifreq >= 0);
	assert(z0 > 0);
	const d_complex_t z = zin(ifeed, ifreq);
	const double gamma = d_abs(d_div(d_sub(z, d_complex(Z0, 0)), d_add(z, d_complex(Z0, 0))));
	double ret = (fabs(1 - gamma) > EPS) ? (1 + gamma) / (1 - gamma) : 1000;
	return ret;
}


// gain, component = 0-6
double gain(int ifreq, double theta, double phi, int component, int db)
{
	assert(ifreq >= 0);
	assert((component >= 0) && (component <= 6));
	double e[7];
	efar(ifreq, theta, phi, e);
	double ret = e[component];
	if (db) {
		ret = 20 * log10(MAX(ret, EPS));
	}
	return ret;
}


double getvalue(int problem)
{
	double f = 0;

	const double z0 = Z0;

	if (!problem) {
		const int ifeed = 0;
		const int ifreq = 0;
		const double theta = 0;
		const double phi = 0;
		const int pol = 1;
		const int db = 1;

		// gain
		f -= gain(ifreq, theta, phi, pol, db);
		//f -= gain(ifreq, 45,   0, 1, 0);
		//f += 0.02 * gain(ifreq, 30, 180, 0, 0);
		//f += 0.02 * gain(ifreq, 50, 180, 0, 0);

/*
		// gain : frequency average
		for (ifreq = 0; ifreq < NFrequency; ifreq++) {
			f -= gain(ifreq, theta, phi, pol, db) / NFrequency;
		}
*/

		// input impedanece
		f += 0.03 * zindiff(ifeed, ifreq, z0);
		//f += 0.01 * zindiff(ifeed, ifreq, z0);

/*
		// input impedanece : frequency average
		for (ifreq = 0; ifreq < NFrequency; ifreq++) {
			f += 0.03 * zindiff(ifeed, ifreq, z0) / NFrequency;
		}
*/
	}


	// sample problem : No. 1...9
	switch (problem) {
		case 1:
			f -= gain(0, 0, 0, 1, 1);
			f += 0.03 * zindiff(0, 0, z0);
			break;
		case 2:
			f -= gain(0, 0, 0, 5, 1);
			f += 0.03 * zindiff(0, 0, z0);
			break;
		case 3:
			f -= gain(0, 45,   0, 1, 1) / 2;
			f -= gain(0, 45, 180, 2, 1) / 2;
			f += 0.02 * zindiff(0, 0, z0);
			break;
		case 4:
			f -= gain(0, 45,   0, 5, 1) / 2;
			f -= gain(0, 45, 180, 6, 1) / 2;
			f += 0.02 * zindiff(0, 0, z0);
			break;
		case 5:
			f -= gain(0, 0, 0, 1, 1) / 2;
			f -= gain(1, 0, 0, 2, 1) / 2;
			for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
				f += 0.03 * zindiff(0, ifreq, z0) / NFrequency;
			}
			break;
		case 6:
			f -= gain(0, 0, 0, 5, 1) / 2;
			f -= gain(1, 0, 0, 6, 1) / 2;
			for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
				f += 0.03 * zindiff(0, ifreq, z0) / NFrequency;
			}
			break;
		case 7:
			for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
				f -= gain(ifreq, 0, 0, 0, 1) / NFrequency;
				f += 0.04 * zindiff(0, ifreq, z0) / NFrequency;
			}
			break;
		case 8:
			f -= gain(0, 45,   0, 1, 1) / 4;
			f -= gain(0, 45,  90, 1, 1) / 4;
			f -= gain(0, 45, 180, 1, 1) / 4;
			f -= gain(0, 45, 270, 1, 1) / 4;
			f += 0.02 * zindiff(0, 0, z0);
			break;
		case 9:
			{
				const int div = 36;
				double gmax = -1000;
				double gmin = +1000;
				for (int i = 0; i < div; i++) {
					const double theta = 45;
					const double phi = 360.0 * i / div;
					const double g = gain(0, theta, phi, 2, 1);
					f -= g / div;
					gmax = MAX(gmax, g);
					gmin = MIN(gmin, g);
				}
				f += 0.3 * fabs(gmax - gmin);
			}
			f += 0.02 * zindiff(0, 0, z0);
			break;
		default:
			break;
	}

	return f;
}
