/*
nearfield.c

near field
*/

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


static void fcurrent(
	const double pos[3], d_complex_t e[3], d_complex_t h[3],
	const d_complex_t cm, const double pc[3], const double ut[3], const double lng, const double kwave)
{
	double r1c[3], r1p[3], r1m[3], d[3];
	double rc, rp, rm, krc, krp, krm, kl, rtmp1, rtmp2, rtmp3;
	d_complex_t gc, gp, gm, fp, fm, ctmp1, ctmp2, ctmp3;

	const double r0 = 0.25 * lng;

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

	for (int k = 0; k < 3; k++) {
		d[k] = pos[k] - pc[k];
	}

	rc = sqrt((d[0] * d[0]) + (d[1] * d[1]) + (d[2] * d[2]));
	if (rc < r0) rc = r0;
	for (int k = 0; k < 3; k++) {
		r1c[k] = d[k] / rc;
	}

	krc = kwave * rc;
	kl = kwave * lng;
/*
	ctmp1 = d_exp(-krc);
	gc = d_rmul(1 / krc, ctmp1);
*/
	ctmp1.r = +cos(krc);
	ctmp1.i = -sin(krc);
	gc.r = ctmp1.r / krc;
	gc.i = ctmp1.i / krc;

	// E field

	// +node
	for (int k = 0; k < 3; k++) {
		d[k] = pos[k] - (pc[k] + (0.5 * lng * ut[k]));
	}
	rp = sqrt((d[0] * d[0]) + (d[1] * d[1]) + (d[2] * d[2]));
	if (rp < r0) rp = r0;
	for (int k = 0; k < 3; k++) {
		r1p[k] = d[k] / rp;
	}

	// -node
	for (int k = 0; k < 3; k++) {
		d[k] = pos[k] - (pc[k] - (0.5 * lng * ut[k]));
	}
	rm = sqrt((d[0] * d[0]) + (d[1] * d[1]) + (d[2] * d[2]));
	if (rm < r0) rm = r0;
	for (int k = 0; k < 3; k++) {
		r1m[k] = d[k] / rm;
	}

	krp = kwave * rp;
	krm = kwave * rm;
/*
	ctmp1 = d_exp(-krp);
	ctmp2 = d_exp(-krm);
	gp = d_rmul(1 / krp, ctmp1);
	gm = d_rmul(1 / krm, ctmp2);
	fp = d_complex(1 / krp, 1);
	fm = d_complex(1 / krm, 1);
*/
	ctmp1.r = +cos(krp);
	ctmp1.i = -sin(krp);
	ctmp2.r = +cos(krm);
	ctmp2.i = -sin(krm);
	gp.r = ctmp1.r / krp;
	gp.i = ctmp1.i / krp;
	gm.r = ctmp2.r / krm;
	gm.i = ctmp2.i / krm;
	fp.r = 1 / krp;
	fp.i = 1;
	fm.r = 1 / krm;
	fm.i = 1;

	for (int k = 0; k < 3; k++) {
/*
		ctmp1 = d_rmul(-ut[k] * kl, gc);
		ctmp2 = d_rmul(-r1p[k], d_mul(fp, gp));
		ctmp3 = d_rmul(+r1m[k], d_mul(fm, gm));
		e[k] = d_add3(ctmp1, ctmp2, ctmp3);
*/
		ctmp1.r = -ut[k] * kl * gc.r;
		ctmp1.i = -ut[k] * kl * gc.i;
		ctmp2.r = -r1p[k] * ((fp.r * gp.r) - (fp.i * gp.i));
		ctmp2.i = -r1p[k] * ((fp.r * gp.i) + (fp.i * gp.r));
		ctmp3.r = +r1m[k] * ((fm.r * gm.r) - (fm.i * gm.i));
		ctmp3.i = +r1m[k] * ((fm.r * gm.i) + (fm.i * gm.r));
		e[k].r = ctmp1.r + ctmp2.r + ctmp3.r;
		e[k].i = ctmp1.i + ctmp2.i + ctmp3.i;
	}

	// E[V/m]
	ctmp1 = d_complex(0, (ETA0 * kwave) / (4 * PI));
	for (int k = 0; k < 3; k++) {
/*
		ctmp2 = d_mul(d_cast(cm), e[k]);
		e[k] = d_mul(ctmp1, ctmp2);
*/
		ctmp2.r = (cm.r * e[k].r) - (cm.i * e[k].i);
		ctmp2.i = (cm.r * e[k].i) + (cm.i * e[k].r);
		e[k].r = (ctmp1.r * ctmp2.r) - (ctmp1.i * ctmp2.i);
		e[k].i = (ctmp1.r * ctmp2.i) + (ctmp1.i * ctmp2.r);
	}

	// H field

	rtmp1 = (r1c[1] * ut[2]) - (r1c[2] * ut[1]);
	rtmp2 = (r1c[2] * ut[0]) - (r1c[0] * ut[2]);
	rtmp3 = (r1c[0] * ut[1]) - (r1c[1] * ut[0]);
/*
	ctmp1 = d_complex(1 / krc, 1);
	ctmp2 = d_rmul(kl, gc);
	ctmp3 = d_mul(ctmp1, ctmp2);
	h[0] = d_rmul(rtmp1, ctmp3);
	h[1] = d_rmul(rtmp2, ctmp3);
	h[2] = d_rmul(rtmp3, ctmp3);
*/
	ctmp1.r = 1 / krc;
	ctmp1.i = 1;
	ctmp2.r = kl * gc.r;
	ctmp2.i = kl * gc.i;
	ctmp3.r = (ctmp1.r * ctmp2.r) - (ctmp1.i * ctmp2.i);
	ctmp3.i = (ctmp1.r * ctmp2.i) + (ctmp1.i * ctmp2.r);
	h[0].r = rtmp1 * ctmp3.r;
	h[0].i = rtmp1 * ctmp3.i;
	h[1].r = rtmp2 * ctmp3.r;
	h[1].i = rtmp2 * ctmp3.i;
	h[2].r = rtmp3 * ctmp3.r;
	h[2].i = rtmp3 * ctmp3.i;

	// H[A/m]
	ctmp1 = d_complex(-kwave / (4 * PI), 0);
	for (int k = 0; k < 3; k++) {
/*
		ctmp2 = d_mul(d_cast(cm), h[k]);
		h[k] = d_mul(ctmp1, ctmp2);
*/
		ctmp2.r = (cm.r * h[k].r) - (cm.i * h[k].i);
		ctmp2.i = (cm.r * h[k].i) + (cm.i * h[k].r);
		h[k].r = (ctmp1.r * ctmp2.r) - (ctmp1.i * ctmp2.i);
		h[k].i = (ctmp1.r * ctmp2.i) + (ctmp1.i * ctmp2.r);
	}
}


void nearfield(const double p[], int ifreq, double kwave, int noinc, d_complex_t ef[], d_complex_t hf[])
{
	d_complex_t e[3], h[3];
	double pos[3], tan[3];

	for (int k = 0; k < 3; k++) {
		ef[k] = hf[k] = d_complex(0, 0);
	}

	// under ground
	if (IGround && ((kwave * p[2]) < -EPS)) {
		return;
	}

	for (int ie = 0; ie < NElement; ie++) {
		d_complex_t c = Cv[(ifreq * NElement) + ie];
		double lng = Element[ie].lng[0];

		for (int k = 0; k < 3; k++) {
			pos[k] = Element[ie].pos[0][k];
			tan[k] = Element[ie].tan[k];
		}

		fcurrent(p, e, h, c, pos, tan, lng, kwave);
		for (int k = 0; k < 3; k++) {
			ef[k] = d_add(ef[k], e[k]);
			hf[k] = d_add(hf[k], h[k]);
		}

		// ground image
		if (IGround) {
			pos[2] *= -1;
			tan[2] *= -1;

			fcurrent(p, e, h, c, pos, tan, lng, kwave);
			for (int k = 0; k < 3; k++) {
				ef[k] = d_sub(ef[k], e[k]);
				hf[k] = d_sub(hf[k], h[k]);
			}
		}
	}

	// plane wave incidence
	if (IPlanewave && !noinc) {
		d_complex_t ei[3], hi[3];
		planewave(p, ei, hi, Planewave.theta, Planewave.phi, Planewave.pol, Planewave.a, Planewave.r, kwave, IGround);
		for (int k = 0; k < 3; k++) {
			ef[k] = d_add(ef[k], ei[k]);
			hf[k] = d_add(hf[k], d_rmul(1 / ETA0, hi[k]));
		}
	}
}
