/*
planewave

plane wave field
*/

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

static d_complex_t einc(double t1, double p1, int pol, double a, double r)
{
	d_complex_t ret;

	if      (pol == 1) {
		// V
		ret = d_complex(-t1, 0);
	}
	else if (pol == 2) {
		// H
		ret = d_complex(+p1, 0);
	}
	else if (pol == 3) {
		// RHCP
		ret = d_rmul(sqrt(0.5), d_complex(+t1, +p1));
	}
	else if (pol == 4) {
		// LHCP
		ret = d_rmul(sqrt(0.5), d_complex(+t1, -p1));
	}
	else if (pol == 5) {
		// elliptical
		const double dtor = atan(1) / 45;
		const double major = + t1 * cos(a * dtor) + p1 * sin(a * dtor);
		const double minor = - t1 * sin(a * dtor) + p1 * cos(a * dtor);
		ret = d_rmul(1 / sqrt(1 + (r * r)), d_complex(major, r * minor));
	}

	return ret;
}

void planewave(const double pos[], d_complex_t ei[], d_complex_t hi[],
               double theta, double phi, int pol, double a, double r, double wavn, int iground)
{
	double r1[3], t1[3], p1[3];
	d_complex_t e[3], h[3], e2[3], h2[3];
	const double dtor = atan(1) / 45;

	if ((iground < 0) || (iground > 1)) return;
	if ((pol < 1) || (pol > 5)) return;

	const double sint = sin(theta * dtor);
	const double cost = cos(theta * dtor);
	const double sinp = sin(phi   * dtor);
	const double cosp = cos(phi   * dtor);

	// unit vector in r, theta and phi

	r1[0] = + sint * cosp;
	r1[1] = + sint * sinp;
	r1[2] = + cost;
	t1[0] = + cost * cosp;
	t1[1] = + cost * sinp;
	t1[2] = - sint;
	p1[0] =        - sinp;
	p1[1] =        + cosp;
	p1[2] =          0;

	// field

	double phs = wavn * ((pos[0] * r1[0]) + (pos[1] * r1[1]) + (pos[2] * r1[2]));
	d_complex_t efac = d_exp(phs);

	for (int k = 0; k < 3; k++) {
		e[k] = d_mul(einc(t1[k], p1[k], pol, a, r), efac);
/*
		if      (pol == 1) {
			e[k] = d_rmul(-t1[k], efac);
		}
		else if (pol == 2) {
			e[k] = d_rmul(+p1[k], efac);
		}
		else if (pol == 3) {
			e[k] = d_mul(d_complex(+root2i * t1[k], +root2i * p1[k]), efac);
		}
		else if (pol == 4) {
			e[k] = d_mul(d_complex(+root2i * t1[k], -root2i * p1[k]), efac);
		}
*/
	}
	h[0] = d_sub(d_rmul(-r1[1], e[2]), d_rmul(-r1[2], e[1]));
	h[1] = d_sub(d_rmul(-r1[2], e[0]), d_rmul(-r1[0], e[2]));
	h[2] = d_sub(d_rmul(-r1[0], e[1]), d_rmul(-r1[1], e[0]));

	// add ground image

	if (iground == 1) {
		r1[2] *= -1;
		t1[2] *= -1;
		phs = wavn * ((pos[0] * r1[0]) + (pos[1] * r1[1]) + (pos[2] * r1[2]));
		efac = d_exp(phs);

		for (int k = 0; k < 3; k++) {
			e2[k] = d_rmul(-1, d_mul(einc(t1[k], p1[k], pol, a, r), efac));
/*
			if      (pol == 1) {
				e2[k] = d_rmul(+t1[k], efac);
			}
			else if (pol == 2) {
				e2[k] = d_rmul(-p1[k], efac);
			}
			else if (pol == 3) {
				e2[k] = d_mul(d_complex(-root2i * t1[k], -root2i * p1[k]), efac);
			}
			else if (pol == 4) {
				e2[k] = d_mul(d_complex(-root2i * t1[k], +root2i * p1[k]), efac);
			}
*/
		}
		h2[0] = d_sub(d_rmul(-r1[1], e2[2]), d_rmul(-r1[2], e2[1]));
		h2[1] = d_sub(d_rmul(-r1[2], e2[0]), d_rmul(-r1[0], e2[2]));
		h2[2] = d_sub(d_rmul(-r1[0], e2[1]), d_rmul(-r1[1], e2[0]));

		for (int k = 0; k < 3; k++) {
			e[k] = d_add(e[k], e2[k]);
			h[k] = d_add(h[k], h2[k]);
		}
	}

	// output

	for (int k = 0; k < 3; k++) {
		ei[k] = e[k];
		hi[k] = h[k];
	}

}
