/*
planewave.c
*/

//#define MAIN

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

// exp(jkr*r), k : complex
d_complex_t expjkr(d_complex_t k, double rr)
{
	return d_rmul(exp(-k.i * rr), d_exp(+k.r * rr));
}


// v3 = v1 X v2
static void outprod(const double v1[3], const double v2[3], double v3[3])
{
	v3[0] = (v1[1] * v2[2]) - (v1[2] * v2[1]);
	v3[1] = (v1[2] * v2[0]) - (v1[0] * v2[2]);
	v3[2] = (v1[0] * v2[1]) - (v1[1] * v2[0]);
}


d_complex_t planewave(int ic, double x, double y, double z, int ifreq)
{
	assert(IPlanewave);
	assert((ic >= 0) && (ic <= 5));

	const double *r1 = Planewave.r1;
	const double *t1 = Planewave.t1;
	const double *p1 = Planewave.p1;
	const int pol = Planewave.pol;
	const double kwave = (2 * PI * Frequency[ifreq]) / C;
	d_complex_t eh[6], ei[3], efi[3];

	for (int m = 0; m < 6; m++) {
		eh[m] = d_complex(0, 0);
	}

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

	if (!Planewave.isemi && !Planewave.ifull) {
		// normal case
		// Ei factor
		double rri = (x * r1[0]) + (y * r1[1]) + (z * r1[2]);
		const d_complex_t fi = d_exp(+kwave * rri);
		//printf("%f %f %f %f\n", z, fi.r, fi.i, d_abs(fi));

		// Ei
		for (int m = 0; m < 3; m++) {
			efi[m] = d_mul(ei[m], fi);
		}

		// E-total
		for (int m = 0; m < 3; m++) {
			eh[m] = efi[m];
		}

		// Z0 * H-total = Ei X r1
		eh[3] = d_sub(d_rmul(r1[2], efi[1]), d_rmul(r1[1], efi[2]));
		eh[4] = d_sub(d_rmul(r1[0], efi[2]), d_rmul(r1[2], efi[0]));
		eh[5] = d_sub(d_rmul(r1[1], efi[0]), d_rmul(r1[0], efi[1]));
	}
	else if (Planewave.ifull) {
		// full space
		d_complex_t epsr0 = Material[Planewave.mfull].cepsr;
		d_complex_t mur0  = Material[Planewave.mfull].cmur;
		d_complex_t k0 = d_rmul(kwave, d_sqrt(d_mul(epsr0, mur0)));

		// Ei factor
		double rri = (x * r1[0]) + (y * r1[1]) + (z * r1[2]);
		//const d_complex_t fi = d_rmul(exp(-k0.i * rri), d_exp(+k0.r * rri));
		const d_complex_t fi = expjkr(k0, rri);

		// Ei
		for (int m = 0; m < 3; m++) {
			efi[m] = d_mul(ei[m], fi);
		}

		// E-total
		for (int m = 0; m < 3; m++) {
			eh[m] = efi[m];
		}

		// Z * H-total = Ei X r1
		eh[3] = d_sub(d_rmul(r1[2], efi[1]), d_rmul(r1[1], efi[2]));
		eh[4] = d_sub(d_rmul(r1[0], efi[2]), d_rmul(r1[2], efi[0]));
		eh[5] = d_sub(d_rmul(r1[1], efi[0]), d_rmul(r1[0], efi[1]));
	}
	else if (Planewave.isemi) {
		// semi space
		assert((Planewave.msemi1 < NMaterial) && (Planewave.msemi2 < NMaterial));

		const int isemi = Planewave.isemi;
		d_complex_t er[3], et[3], efr[3], eft[3];
		double ri[3], si[3], pi[3];
		double rr[3], sr[3], pr[3];
		double rt[3], st[3], pt[3];

		// normal vector
		double n[3] = {0, 0, 0};
		if (isemi == 1) n[0] = (r1[0] >= 0) ? +1 : -1;
		if (isemi == 2) n[1] = (r1[1] >= 0) ? +1 : -1;
		if (isemi == 3) n[2] = (r1[2] >= 0) ? +1 : -1;
		//printf("%f %f %f\n", n[0], n[1], n[2]);

		// reflection plane
		const double x0 = (isemi == 1) ? Planewave.psemi : 0;
		const double y0 = (isemi == 2) ? Planewave.psemi : 0;
		const double z0 = (isemi == 3) ? Planewave.psemi : 0;
		//printf("%f %f %f\n", x0, y0, z0);

		// ri = -r1, pi = -t1, si = -p1
		// r1, t1, p1 : RHS (r1 X t1 = p1)
		// ri, si, pi : RHS (pi = ri X si = -r1 X -p1 = -(p1 X r1) = -t1 !!)
		for (int m = 0; m < 3; m++) {
			ri[m] = -r1[m];
			pi[m] = -t1[m];
			si[m] = -p1[m];
		}
		//printf("ri : %f %f %f\n", ri[0], ri[1], ri[2]);
		//printf("pi : %f %f %f\n", pi[0], pi[1], pi[2]);
		//printf("si : %f %f %f\n", si[0], si[1], si[2]);

		// sr = st = si
		for (int m = 0; m < 3; m++) {
			sr[m] = st[m] = si[m];
		}
		//printf("sr : %f %f %f\n", sr[0], sr[1], sr[2]);
		//printf("st : %f %f %f\n", st[0], st[1], st[2]);

		// reflection ?
		const int ref = ((isemi == 1) && (((x - x0) * ri[0]) <= 0)) ||
		                ((isemi == 2) && (((y - y0) * ri[1]) <= 0)) ||
		                ((isemi == 3) && (((z - z0) * ri[2]) <= 0));
		//printf("%f %d\n", z, ref);

		// 1 -> 2
		const int m0 = ifreq * NMaterial;
		const d_complex_t epsr1 = Material[Planewave.msemi1 + m0].cepsr;
		const d_complex_t mur1  = Material[Planewave.msemi1 + m0].cmur;
		const d_complex_t epsr2 = Material[Planewave.msemi2 + m0].cepsr;
		const d_complex_t mur2  = Material[Planewave.msemi2 + m0].cmur;
		const d_complex_t eps21 = d_div(epsr2, epsr1);
		const d_complex_t mu21 = d_div(mur2, mur1);
		const double n12 = d_sqrt(d_mul(epsr1, mur1)).r / d_sqrt(d_mul(epsr2, mur2)).r;
		const d_complex_t k1 = d_rmul(kwave, d_sqrt(d_mul(epsr1, mur1)));
		//const d_complex_t k2 = d_rmul(kwave, d_sqrt(d_mul(epsr2, mur2)));
		//printf("%f\n", n12);
		//printf("%f %f %f %f %f\n", kwave, k1.r, k1.i, k2.r, k2.i);
		//printf("%f %f %f %f\n", eps21.r, eps21.i, mu21.r, mu21.i);

		// Ei factor
		const double rri = ((x - x0) * ri[0]) + ((y - y0) * ri[1]) + ((z - z0) * ri[2]);
		//const d_complex_t fi = d_rmul(exp(+k1.i * rri), d_exp(-k1.r * rri));
		//const d_complex_t fi = expjkr(k1, rri);
		//const double dli = (isemi == 1) ? ri[0] * (x - x0)
		//                 : (isemi == 2) ? ri[1] * (y - y0)
		//                 : (isemi == 3) ? ri[2] * (z - z0) : 0;
		//const d_complex_t fi = d_rmul(exp(k1.i * fabs(dli)), d_exp(-k1.r * rri));
		const d_complex_t fi = expjkr(d_neg(k1), rri);

		// Ei
		for (int m = 0; m < 3; m++) {
			efi[m] = d_mul(ei[m], fi);
		}
		//printf("ei : %f %f %f %f %f %f %f\n", z, efi[0].r, efi[0].i, efi[1].r, efi[1].i, efi[2].r, efi[2].i);
		//printf("ei : %f %f %f %f\n", z, efh[0].r, efi[0].i, d_abs(efi[0]));

		// cos-i = -n * ri
		// sin-i = sqrt(1 - cos-i * cos-i)
		const double cosi = -((n[0] * ri[0]) + (n[1] * ri[1]) + (n[2] * ri[2]));
		const double sini = sqrt(MAX(1 - (cosi * cosi), 0));
		const int total_reflection = (n12 * sini > 1);
		//printf("%f %f\n", cosi, sini);

		// rr = ri + (2 * cos-i) * n
		for (int m = 0; m < 3; m++) {
			rr[m] = ri[m] + (2 * cosi) * n[m];
		}
		//printf("rr : %f %f %f\n", rr[0], rr[1], rr[2]);

		// pr = rr X sr
		outprod(rr, sr, pr);
		//printf("pr : %f %f %f\n", pr[0], pr[1], pr[2]);

		// ei * pi, ei * si
		const d_complex_t eipi = d_add3(d_rmul(pi[0], ei[0]), d_rmul(pi[1], ei[1]), d_rmul(pi[2], ei[2]));
		const d_complex_t eisi = d_add3(d_rmul(si[0], ei[0]), d_rmul(si[1], ei[1]), d_rmul(si[2], ei[2]));
		//printf("eipi : %f %f %f %f\n", z, eipi.r, eipi.i, d_abs(eipi));
		//printf("eisi : %f %f %f %f\n", z, eisi.r, eisi.i, d_abs(eisi));

		const d_complex_t c1p = d_rmul(cosi, eps21);
		const d_complex_t c1s = d_rmul(cosi, mu21);

		d_complex_t root = d_sqrt(d_sub(d_mul(eps21, mu21), d_complex(sini * sini, 0)));
		if (total_reflection) {
			root = d_rmul(-1, root);  // negative principal
		}

		if (ref) {
			// reflection
			d_complex_t rp, rs;
			if (Planewave.msemi2 == PEC) {
				// 2 : PEC
				rp = d_complex(+1, 0);
				rs = d_complex(-1, 0);
			}
/*
			else if (n12 * sini > 1) {
				// total reflection
				rp = d_complex(-1, 0);
				rs = d_complex(+1, 0);
			}
*/
			else {
				// Fresnel reflection (Rp, Rs)
				rp = d_div(d_sub(c1p, root), d_add(c1p, root));
				rs = d_div(d_sub(c1s, root), d_add(c1s, root));
				//printf("Rp : %f %f %f %f\n", z, rp.r, rp.i, d_abs(rp));
				//printf("Rs : %f %f %f %f\n", z, rs.r, rs.i, d_abs(rs));
			}

			// er = Rp * (ei * pi) * pr + Rs * (ei * si) * sr
			for (int m = 0; m < 3; m++) {
				er[m] = d_add(d_rmul(pr[m], d_mul(rp, eipi)),
				              d_rmul(sr[m], d_mul(rs, eisi)));
			}
			//printf("er : %f %f %f %f %f %f %f\n", z, er[0].r, er[0].i, er[1].r, er[1].i, er[2].r, er[2].i);
			//printf("%f %f\n", z, sqrt((er[0].r * er[0].r) + (er[1].r * er[1].r) + (er[2].r * er[2].r)));

			// Er factor
			const double rrr = ((x - x0) * rr[0]) + ((y - y0) * rr[1]) + ((z - z0) * rr[2]);
			const d_complex_t fr = expjkr(d_neg(k1), rrr);

			// Er
			for (int m = 0; m < 3; m++) {
				efr[m] = d_mul(er[m], fr);
			}

			// E-total = Ei + Er
			for (int m = 0; m < 3; m++) {
				eh[m] = d_add(efi[m], efr[m]);
			}
			//printf("%f %f %f %f %f\n", z, ef[0].r, ef[0].i, d_mul(er[0], fr).r, d_mul(er[0], fr).i);
			//printf("%f %f %f %f\n", z, d_abs(ef[0]), d_abs(ef[1]), d_abs(ef[2]));

			// Z * H-total = ri X Ei + rr X Er
			eh[3] = d_add(d_sub(d_rmul(ri[1], efi[2]), d_rmul(ri[2], efi[1])),
			              d_sub(d_rmul(rr[1], efr[2]), d_rmul(rr[2], efr[1])));
			eh[4] = d_add(d_sub(d_rmul(ri[2], efi[0]), d_rmul(ri[0], efi[2])),
			              d_sub(d_rmul(rr[2], efr[0]), d_rmul(rr[0], efr[2])));
			eh[5] = d_add(d_sub(d_rmul(ri[0], efi[1]), d_rmul(ri[1], efi[0])),
			              d_sub(d_rmul(rr[0], efr[1]), d_rmul(rr[1], efr[0])));
		}
		else {
			// transmision
			if (Planewave.msemi2 == PEC) {
				// PEC : Et = 0
				for (int m = 0; m < 6; m++) {
					eh[m] = d_complex(0, 0);
				}
			}
/*
			else if (n12 * sini > 1) {
				// total reflection : Et = 0
				for (int m = 0; m < 6; m++) {
					eh[m] = d_complex(0, 0);
				}
			}
*/
			else {
				// Snell's law
				// sin-t = (n1 / n2) * sin-i
				// cos-t = sqrt(1 - sin-i * sin-i)
				//const double sint = n12 * sini;
				//const double cost = sqrt(MAX(1 - (sint * sint), 0));
				//printf("%f %f %f\n", sint, cost, asin(sint) / DTOR);

				// Fresnel transmission (Tp, Ts)
				const d_complex_t c2p = d_rmul(cosi, d_sqrt(d_mul(eps21, mu21)));
				const d_complex_t tp = d_div(d_rmul(2, c2p), d_add(c1p, root));
				const d_complex_t ts = d_div(d_rmul(2, c1s), d_add(c1s, root));
				//printf("Tp : %f %f %f %f\n", z, tp.r, tp.i, d_abs(tp));
				//printf("Ts : %f %f %f %f\n", z, ts.r, ts.i, d_abs(ts));

				// rt
				if (!total_reflection) {
					// Snell's law : sin-t = (n1 / n2) * sin-i
					// rt = (n1 / n2) * ri + ((n1 / n2) * cos-i - cos-t) * n
					const double sint = n12 * sini;
					const double cost = sqrt(MAX(1 - (sint * sint), 0));
					for (int m = 0; m < 3; m++) {
						rt[m] = (n12 * ri[m]) + (n12 * cosi - cost) * n[m];
					}
					//printf("%f %f %f %f\n", rt[0], rt[1], rt[2], sqrt((rt[0] * rt[0]) + (rt[1] * rt[1]) + (rt[2] * rt[2])));
				}
				else {
					// total reflection : rt = tangential vector
					for (int m = 0; m < 3; m++) {
						rt[m] = (ri[m] + (n[m] * cosi)) / sini;
					}
				}

				// pt = rt X st
				outprod(rt, st, pt);
				//printf("pt : %f %f %f\n", pt[0], pt[1], pt[2]);

				// et = Tp * (ei * pi) * pt + Ts * (ei * si) * st
				for (int m = 0; m < 3; m++) {
					et[m] = d_add(d_rmul(pt[m], d_mul(tp, eipi)),
					              d_rmul(st[m], d_mul(ts, eisi)));
				}
				//printf("%f %f %f %f %f %f\n", et[0].r, et[0].i, et[1].r, et[1].i, et[2].r, et[2].i);
				//printf("%f %f\n", z, sqrt((et[0].r * et[0].r) + (et[1].r * et[1].r) + (et[2].r * et[2].r)));

				// Et factor : ft
				//const double rrt = ((x - x0) * rt[0]) + ((y - y0) * rt[1]) + ((z - z0) * rt[2]);
				//const d_complex_t ft = expjkr(d_neg(k2), rrt);
				const double rn = ((x - x0) * n[0]) + ((y - y0) * n[1]) + ((z - z0) * n[2]);
				d_complex_t rrt = d_add(d_complex(rri, 0), d_rmul(rn, d_sub(d_complex(cosi, 0), root)));
				rrt = d_mul(d_complex(0, -1), d_mul(k1, rrt));
				const d_complex_t ft = d_rmul(exp(rrt.r), d_exp(rrt.i));

				// Et
				for (int m = 0; m < 3; m++) {
					eft[m] = d_mul(et[m], ft);
				}
				//printf("%f %f %f %f %f %f %f\n", z, eft[0].r, eft[0].i, eft[1].r, eft[1].i, eft[2].r, eft[2].i);

				// E-total = Et
				for (int m = 0; m < 3; m++) {
					eh[m] = eft[m];
				}

				// Z * H-total = rt X Et
				eh[3] = d_sub(d_rmul(rt[1], eft[2]), d_rmul(rt[2], eft[1]));
				eh[4] = d_sub(d_rmul(rt[2], eft[0]), d_rmul(rt[0], eft[2]));
				eh[5] = d_sub(d_rmul(rt[0], eft[1]), d_rmul(rt[1], eft[0]));
			}
		}
	}

	//printf("%d %f %f %f\n", ic, z, eh[ic].r, eh[ic].i);
	return eh[ic];
}

/*
// cl.exe -I../include /nologo /W4 planewave.c /Fea.exe
int main(void)
{
	IPlanewave = 1;
	Frequency = (double *)malloc(1 * sizeof(double));
	Frequency[0] = 2e9;

	const double theta = 0;
	const double phi   = 0;

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

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

	Planewave.pol = 1;

	Planewave.isemi = 3;
	Planewave.psemi = 0;
	Planewave.msemi1 = 0;
	Planewave.msemi2 = 2;

	NMaterial = 3;
	Material = (material_t *)malloc(NMaterial * sizeof(material_t));
	Material[0].cepsr = d_complex(1, -0);
	Material[0].cmur  = d_complex(1, -0);
	Material[2].cepsr = d_complex(2, -0);
	Material[2].cmur  = d_complex(1, -0);

	double x, y, z;
	x = 0.0e-3;
	y = 0.0e-3;
	for (int i = -20; i <= 20; i++) {
		z = 5e-3 * i;
		d_complex_t ex = planewave(0, x, y, z, 0);
		d_complex_t ey = planewave(1, x, y, z, 0);
		d_complex_t ez = planewave(2, x, y, z, 0);
		d_complex_t hx = planewave(3, x, y, z, 0);
		d_complex_t hy = planewave(4, x, y, z, 0);
		d_complex_t hz = planewave(5, x, y, z, 0);
		//printf("%f %f %f %f %f %f %f\n", z, ex.r, ex.i, d_abs(ex), ez.r, ez.i, d_abs(ez));
		//printf("%f %f %f %f %f %f %f\n", z, ex.r, ex.i, ey.r, ey.i, ez.r, ez.i);
		//printf("%f %f %f %f %f\n", z, d_abs(ex), d_abs(ey), d_abs(ez), sqrt(d_norm(ex) + d_norm(ey) + d_norm(ez)));
		printf("%f %f %f\n", z, sqrt(d_norm(ex) + d_norm(ey) + d_norm(ez)), sqrt(d_norm(hx) + d_norm(hy) + d_norm(hz)));
	}

	return 0;
}
*/
