/*
planewaveEy.c
*/

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

void planewaveEy(d_complex_t *b, int ifreq)
{
	assert(IPlanewave);
	const int64_t m0 = ifreq * NMaterial;

	int i;
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (    i = i0Ey; i <= i1Ey; i++) {
	for (int j = j0Ey; j <= j1Ey; j++) {
	for (int k = k0Ey; k <= k1Ey; k++) {
		const int64_t n = NEY(i, j, k);
		if (MEY(i, j, k) == PEC) {
			b[n] = d_neg(planewave(1, Xn[i], Yc[j], Zn[k], ifreq));
		}
		else {
			// exclude Mur and PML points, cf. matrixEy.c
			if ((Abc == MUR) && (
				(!Pbc.zpbc && ((k == 0) || (k == Nz))) ||
				(!Pbc.xpbc && ((i == 0) || (i == Nx))))) {
			}
			else if ((Abc == PML) &&
				((j < 0) || (j >= Ny) || (k <= 0) || (k >= Nz) || (i <= 0) || (i >= Nx))) {
			}
			else if (!Planewave.isemi) {
				// normal case
				const d_complex_t c0 = Material[MEY(i,     j,     k    ) + m0].cepsr;
				const d_complex_t c1 = Material[MHX(i,     j,     k    ) + m0].cnur;
				const d_complex_t c2 = Material[MHX(i,     j,     k - 1) + m0].cnur;
				const d_complex_t c3 = Material[MHZ(i,     j,     k    ) + m0].cnur;
				const d_complex_t c4 = Material[MHZ(i - 1, j,     k    ) + m0].cnur;

				// background
				const d_complex_t kf = Planewave.ifull
					? d_sqrt(d_mul(Material[Planewave.mfull].cepsr, Material[Planewave.mfull].cmur))
					: d_complex(1, 0);

				const d_complex_t nur = d_mul(d_rmul(0.25, d_add4(c1, c2, c3, c4)), d_mul(kf, kf));
				d_complex_t ey = planewave(1, Xn[i], Yc[j], Zn[k], ifreq);
				d_complex_t hz = planewave(5, Xn[i], Yc[j], Zn[k], ifreq);
				d_complex_t hx = planewave(3, Xn[i], Yc[j], Zn[k], ifreq);
				ey = d_mul(d_sub(c0, nur), ey);
				hx = d_mul(hx, d_mul(kf, d_mul(d_complex(0, KRZn[k]), d_sub(c1, c2))));
				hz = d_mul(hz, d_mul(kf, d_mul(d_complex(0, KRXn[i]), d_sub(c3, c4))));
				b[n] = d_add(ey, d_sub(hx, hz));
			}
			else {
				// semi-space
				const d_complex_t c0 = Material[MEY(i,     j,     k    ) + m0].cepsr;
				const d_complex_t c1 = Material[MHX(i,     j,     k    ) + m0].cnur;
				const d_complex_t c2 = Material[MHX(i,     j,     k - 1) + m0].cnur;
				const d_complex_t c3 = Material[MHZ(i,     j,     k    ) + m0].cnur;
				const d_complex_t c4 = Material[MHZ(i - 1, j,     k    ) + m0].cnur;

				d_complex_t c[NNZ];
				c[ 0] = d_neg(c0);
				c[ 1] = d_rmul(- KRZp[k],           c1);
				c[ 2] = d_rmul(- KRZm[k],           c2);
				c[ 3] = d_rmul(- KRXp[i],           c3);
				c[ 4] = d_rmul(- KRXm[i],           c4);
				c[ 5] = d_rmul(+ KRYc[j] * KRZn[k], c1);
				c[ 6] = d_neg(c[ 5]);
				c[ 7] = d_rmul(- KRYc[j] * KRZn[k], c2);
				c[ 8] = d_neg(c[ 7]);
				c[ 9] = d_rmul(+ KRYc[j] * KRXn[i], c3);
				c[10] = d_neg(c[ 9]);
				c[11] = d_rmul(- KRYc[j] * KRXn[i], c4);
				c[12] = d_neg(c[11]);
				c[ 0] = d_sub(c[ 0], d_add4(c[ 1], c[ 2], c[ 3], c[ 4]));

				const double znc = Zn[k];
				const double znm = (k >  0) ? Zn[k - 1] : (2 * Zn[ 0] - Zn[     1]);
				const double znp = (k < Nz) ? Zn[k + 1] : (2 * Zn[Nz] - Zn[Nz - 1]);
				const double zcm = (k >  0) ? Zc[k - 1] : (2 * Zn[ 0] - Zc[     0]);
				const double zcp = (k < Nz) ? Zc[k    ] : (2 * Zn[Nz] - Zc[Nz - 1]);
				const double xnc = Xn[i];
				const double xnm = (i >  0) ? Xn[i - 1] : (2 * Xn[ 0] - Xn[     1]);
				const double xnp = (i < Nx) ? Xn[i + 1] : (2 * Xn[Nx] - Xn[Nx - 1]);
				const double xcm = (i >  0) ? Xc[i - 1] : (2 * Xn[ 0] - Xc[     0]);
				const double xcp = (i < Nx) ? Xc[i    ] : (2 * Xn[Nx] - Xc[Nx - 1]);

				d_complex_t e[NNZ];
				e[ 0] = planewave(1, xnc, Yc[j    ], znc, ifreq);
				e[ 1] = planewave(1, xnc, Yc[j    ], znp, ifreq);
				e[ 2] = planewave(1, xnc, Yc[j    ], znm, ifreq);
				e[ 3] = planewave(1, xnp, Yc[j    ], znc, ifreq);
				e[ 4] = planewave(1, xnm, Yc[j    ], znc, ifreq);
				e[ 5] = planewave(2, xnc, Yn[j + 1], zcp, ifreq);
				e[ 6] = planewave(2, xnc, Yn[j    ], zcp, ifreq);
				e[ 7] = planewave(2, xnc, Yn[j + 1], zcm, ifreq);
				e[ 8] = planewave(2, xnc, Yn[j    ], zcm, ifreq);
				e[ 9] = planewave(0, xcp, Yn[j + 1], znc, ifreq);
				e[10] = planewave(0, xcp, Yn[j    ], znc, ifreq);
				e[11] = planewave(0, xcm, Yn[j + 1], znc, ifreq);
				e[12] = planewave(0, xcm, Yn[j    ], znc, ifreq);

				d_complex_t sum = d_complex(0, 0);
				for (int m = 0; m < NNZ; m++) {
					sum = d_add(sum, d_mul(c[m], e[m]));
				}
				b[n] = d_neg(sum);
			}
		}
	}
	}
	}
}
