/*
farfield.c

far field
*/

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

typedef struct {
	double nx, ny, nz;
	double x, y, z;
	double ds;
	int    m;
} surface_t;
surface_t *surface;
d_complex_t **surface_ex, **surface_ey, **surface_ez;
d_complex_t **surface_hx, **surface_hy, **surface_hz;
int64_t nsurface;

void farfield(int ifreq, double theta, double phi, double ffctr, d_complex_t *etheta, d_complex_t *ephi)
{
	// wave number

	const double kwave = (2 * PI * Frequency[ifreq]) / C;

	// unit vector in r, theta, phi

	double r1[3], t1[3], p1[3];

	const double cos_t = cos(theta * DTOR);
	const double sin_t = sin(theta * DTOR);
	const double cos_p = cos(phi   * DTOR);
	const double sin_p = sin(phi   * DTOR);

	r1[0] = +sin_t * cos_p;
	r1[1] = +sin_t * sin_p;
	r1[2] = +cos_t;

	t1[0] = +cos_t * cos_p;
	t1[1] = +cos_t * sin_p;
	t1[2] = -sin_t;

	p1[0] = -sin_p;
	p1[1] = +cos_p;
	p1[2] = 0;

	d_complex_t plx = d_complex(0, 0);
	d_complex_t ply = d_complex(0, 0);
	d_complex_t plz = d_complex(0, 0);
	d_complex_t pnx = d_complex(0, 0);
	d_complex_t pny = d_complex(0, 0);
	d_complex_t pnz = d_complex(0, 0);

	d_complex_t *ex = surface_ex[ifreq];
	d_complex_t *ey = surface_ey[ifreq];
	d_complex_t *ez = surface_ez[ifreq];
	d_complex_t *hx = surface_hx[ifreq];
	d_complex_t *hy = surface_hy[ifreq];
	d_complex_t *hz = surface_hz[ifreq];

	for (int n = 0; n < nsurface; n++) {
		surface_t *p = &surface[n];
		const material_t m = Material[p->m + (ifreq * NMaterial)];

		// Z0 * J = n X (Z0 * H)
		const d_complex_t cjx = d_sub(d_rmul(+p->ny, hz[n]), d_rmul(+p->nz, hy[n]));
		const d_complex_t cjy = d_sub(d_rmul(+p->nz, hx[n]), d_rmul(+p->nx, hz[n]));
		const d_complex_t cjz = d_sub(d_rmul(+p->nx, hy[n]), d_rmul(+p->ny, hx[n]));

		// M = -n X E
		const d_complex_t cmx = d_sub(d_rmul(-p->ny, ez[n]), d_rmul(-p->nz, ey[n]));
		const d_complex_t cmy = d_sub(d_rmul(-p->nz, ex[n]), d_rmul(-p->nx, ez[n]));
		const d_complex_t cmz = d_sub(d_rmul(-p->nx, ey[n]), d_rmul(-p->ny, ex[n]));

		// exp(jkr * r) * dS
		const d_complex_t k = d_rmul(kwave, d_sqrt(d_mul(m.cepsr, m.cmur)));
		const double rr = (r1[0] * p->x) + (r1[1] * p->y) + (r1[2] * p->z);
		const d_complex_t expds = d_rmul(p->ds * exp(-k.i * rr), d_exp(k.r * rr));

		// L += M * exp(jkr * r) * dS
		plx = d_add(plx, d_mul(cmx, expds));
		ply = d_add(ply, d_mul(cmy, expds));
		plz = d_add(plz, d_mul(cmz, expds));

		// Z * N += (Z * J) * exp(jkr * r) * dS
		d_complex_t z = d_sqrt(d_div(m.cmur, m.cepsr));
		pnx = d_add(pnx, d_mul(z, d_mul(cjx, expds)));
		pny = d_add(pny, d_mul(z, d_mul(cjy, expds)));
		pnz = d_add(pnz, d_mul(z, d_mul(cjz, expds)));
	}

	// Z * N-theta, Z * N-phi
	const d_complex_t pnt = d_add3(d_rmul(t1[0], pnx), d_rmul(t1[1], pny), d_rmul(t1[2], pnz));
	const d_complex_t pnp = d_add3(d_rmul(p1[0], pnx), d_rmul(p1[1], pny), d_rmul(p1[2], pnz));

	// L-theta, L-phi
	const d_complex_t plt = d_add3(d_rmul(t1[0], plx), d_rmul(t1[1], ply), d_rmul(t1[2], plz));
	const d_complex_t plp = d_add3(d_rmul(p1[0], plx), d_rmul(p1[1], ply), d_rmul(p1[2], plz));

	// F-theta, F-phi
	*etheta = d_rmul(ffctr, d_add(pnt, plp));
	*ephi   = d_rmul(ffctr, d_sub(pnp, plt));
}


// alloc far field array
void alloc_farfield(void)
{
	const int64_t nsx = !Pbc.xpbc ? (Ny - 1) * (Nz - 1) : 0;
	const int64_t nsy = !Pbc.ypbc ? (Nz - 1) * (Nx - 1) : 0;
	const int64_t nsz = !Pbc.zpbc ? (Nx - 1) * (Ny - 1) : 0;
	nsurface = 2 * (nsx + nsy + nsz);

	surface = (surface_t *)malloc(nsurface * sizeof(surface_t));

	surface_ex = (d_complex_t **)malloc(NFrequency * sizeof(d_complex_t *));
	surface_ey = (d_complex_t **)malloc(NFrequency * sizeof(d_complex_t *));
	surface_ez = (d_complex_t **)malloc(NFrequency * sizeof(d_complex_t *));
	surface_hx = (d_complex_t **)malloc(NFrequency * sizeof(d_complex_t *));
	surface_hy = (d_complex_t **)malloc(NFrequency * sizeof(d_complex_t *));
	surface_hz = (d_complex_t **)malloc(NFrequency * sizeof(d_complex_t *));
	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		surface_ex[ifreq] = (d_complex_t *)malloc(nsurface * sizeof(d_complex_t));
		surface_ey[ifreq] = (d_complex_t *)malloc(nsurface * sizeof(d_complex_t));
		surface_ez[ifreq] = (d_complex_t *)malloc(nsurface * sizeof(d_complex_t));
		surface_hx[ifreq] = (d_complex_t *)malloc(nsurface * sizeof(d_complex_t));
		surface_hy[ifreq] = (d_complex_t *)malloc(nsurface * sizeof(d_complex_t));
		surface_hz[ifreq] = (d_complex_t *)malloc(nsurface * sizeof(d_complex_t));
	}
}


// setup surface E and H
void setup_farfield(void)
{
	const int exclude = 1;  // 1:Es

	// counter
	int64_t n = 0;

	// X surface
	for (int side = 0; side < 2; side++) {
		if (Pbc.xpbc) continue;
		const int i = (side == 0) ? 0 : Nx - 1;
		for (int j = 1; j < Ny; j++) {
		for (int k = 1; k < Nz; k++) {
			for (int f = 0; f < NFrequency; f++) {
				surface_ex[f][n] = d_complex(0, 0);
				surface_ey[f][n] = d_rmul(0.25, d_add4(eyyee(f, i,     j,     k,     exclude),
				                                       eyyee(f, i + 1, j,     k,     exclude),
				                                       eyyee(f, i,     j - 1, k,     exclude),
				                                       eyyee(f, i + 1, j - 1, k,     exclude)));
				surface_ez[f][n] = d_rmul(0.25, d_add4(ezyee(f, i,     j,     k,     exclude),
				                                       ezyee(f, i + 1, j,     k,     exclude),
				                                       ezyee(f, i,     j,     k - 1, exclude),
				                                       ezyee(f, i + 1, j,     k - 1, exclude)));
				surface_hx[f][n] = d_complex(0, 0);
				surface_hy[f][n] = d_rmul(0.5,   d_add(hyyee(f, i,     j,     k,     exclude),
				                                       hyyee(f, i,     j,     k - 1, exclude)));
				surface_hz[f][n] = d_rmul(0.5,   d_add(hzyee(f, i,     j,     k,     exclude),
				                                       hzyee(f, i,     j - 1, k,     exclude)));
			}
			surface[n].nx = (side == 0) ? -1 : +1;
			surface[n].ny = 0;
			surface[n].nz = 0;
			surface[n].x = Xc[i];
			surface[n].y = Yn[j];
			surface[n].z = Zn[k];
			surface[n].ds = DYn[j] * DZn[k];
			surface[n].m = MEX(i, j, k);
			n++;
		}
		}
	}

	// Y surface
	for (int side = 0; side < 2; side++) {
		if (Pbc.ypbc) continue;
		const int j = (side == 0) ? 0 : Ny - 1;
		for (int k = 1; k < Nz; k++) {
		for (int i = 1; i < Nx; i++) {
			for (int f = 0; f < NFrequency; f++) {
				surface_ey[f][n] = d_complex(0, 0);
				surface_ez[f][n] = d_rmul(0.25, d_add4(ezyee(f, i,     j,     k,     exclude),
				                                       ezyee(f, i,     j + 1, k,     exclude),
				                                       ezyee(f, i,     j,     k - 1, exclude),
				                                       ezyee(f, i,     j + 1, k - 1, exclude)));
				surface_ex[f][n] = d_rmul(0.25, d_add4(exyee(f, i,     j,     k,     exclude),
				                                       exyee(f, i,     j + 1, k,     exclude),
				                                       exyee(f, i - 1, j,     k,     exclude),
				                                       exyee(f, i - 1, j + 1, k,     exclude)));
				surface_hy[f][n] = d_complex(0, 0);
				surface_hz[f][n] = d_rmul(0.5,   d_add(hzyee(f, i,     j,     k,     exclude),
				                                       hzyee(f, i - 1, j,     k,     exclude)));
				surface_hx[f][n] = d_rmul(0.5,   d_add(hxyee(f, i,     j,     k,     exclude),
				                                       hxyee(f, i,     j,     k - 1, exclude)));
			}
			surface[n].nx = 0;
			surface[n].ny = (side == 0) ? -1 : +1;
			surface[n].nz = 0;
			surface[n].x = Xn[i];
			surface[n].y = Yc[j];
			surface[n].z = Zn[k];
			surface[n].ds = DZn[k] * DXn[i];
			surface[n].m = MEY(i, j, k);
			n++;
		}
		}
	}

	// Z surface
	for (int side = 0; side < 2; side++) {
		if (Pbc.zpbc) continue;
		const int k = (side == 0) ? 0 : Nz - 1;
		for (int i = 1; i < Nx; i++) {
		for (int j = 1; j < Ny; j++) {
			for (int f = 0; f < NFrequency; f++) {
				surface_ez[f][n] = d_complex(0, 0);
				surface_ex[f][n] = d_rmul(0.25, d_add4(exyee(f, i,     j,     k,     exclude),
				                                       exyee(f, i,     j,     k + 1, exclude),
				                                       exyee(f, i - 1, j,     k,     exclude),
				                                       exyee(f, i - 1, j,     k + 1, exclude)));
				surface_ey[f][n] = d_rmul(0.25, d_add4(eyyee(f, i,     j,     k,     exclude),
				                                       eyyee(f, i,     j,     k + 1, exclude),
				                                       eyyee(f, i,     j - 1, k,     exclude),
				                                       eyyee(f, i,     j - 1, k + 1, exclude)));
				surface_hz[f][n] = d_complex(0, 0);
				surface_hx[f][n] = d_rmul(0.5,   d_add(hxyee(f, i,     j,     k,     exclude),
				                                       hxyee(f, i,     j - 1, k,     exclude)));
				surface_hy[f][n] = d_rmul(0.5,   d_add(hyyee(f, i,     j,     k,     exclude),
				                                       hyyee(f, i - 1, j,     k,     exclude)));
			}
			surface[n].nx = 0;
			surface[n].ny = 0;
			surface[n].nz = (side == 0) ? -1 : +1;
			surface[n].x = Xn[i];
			surface[n].y = Yn[j];
			surface[n].z = Zc[k];
			surface[n].ds = DXn[i] * DYn[j];
			surface[n].m = MEZ(i, j, k);
			n++;
		}
		}
	}
}


// far field components
void farComponent(d_complex_t etheta, d_complex_t ephi, double e[])
{
	// abs
	e[0] = sqrt(d_norm(etheta) + d_norm(ephi));

	// theat/phi
	e[1] = d_abs(etheta);
	e[2] = d_abs(ephi);

	// major/minor
	double tmp = d_abs(d_add(d_mul(etheta, etheta), d_mul(ephi, ephi)));
	e[3] = sqrt((d_norm(etheta) + d_norm(ephi) + tmp) / 2);
	e[4] = sqrt((d_norm(etheta) + d_norm(ephi) - tmp) / 2);

	// RHCP/LHCP
	e[5] = d_abs(d_add(etheta, d_mul(d_complex(0, 1), ephi))) / sqrt(2);
	e[6] = d_abs(d_sub(etheta, d_mul(d_complex(0, 1), ephi))) / sqrt(2);
}


// far field factor
void farfactor(double *ffctr)
{
	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		const double kwave = (2 * PI * Frequency[ifreq]) / C;
		if (NFeed) {
			double sum = 0;
			for (int ifeed = 0; ifeed < NFeed; ifeed++) {
				const d_complex_t zin = Zin[(ifeed * NFrequency) + ifreq];
				double denom = 1;
				if (MatchingLoss) {
					const double z0 = Feed[ifeed].z0;
					denom = 1
						 - ((zin.r - z0) * (zin.r - z0) + (zin.i * zin.i))
						 / ((zin.r + z0) * (zin.r + z0) + (zin.i * zin.i));
					denom = MAX(denom, EPS);
				}
				const double gin = d_inv(zin).r;
				sum += 0.5 * gin * d_norm(Feed[ifeed].volt) / denom;
			}
			ffctr[ifreq] = kwave / sqrt(8 * PI * ETA0 * sum);
		}
		else if (IPlanewave) {
			const double einc = 1;
			ffctr[ifreq] = kwave / (einc * sqrt(4 * PI));
		}
	}
}
