/*
planewaveEz.c
*/

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

void planewaveEz(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 = i0Ez; i <= i1Ez; i++) {
	for (int j = j0Ez; j <= j1Ez; j++) {
	for (int k = k0Ez; k <= k1Ez; k++) {
		const int64_t n = NEZ(i, j, k);
		if (MEZ(i, j, k) == PEC) {
			b[n] = d_neg(planewave(2, Xn[i], Yn[j], Zc[k], ifreq));
		}
		else {
			// exclude Mur and PML points, cf. matrixEz.c
			if ((Abc == MUR) && (
				(!Pbc.xpbc && ((i == 0) || (i == Nx))) ||
				(!Pbc.ypbc && ((j == 0) || (j == Ny))))) {
			}
			else if ((Abc == PML) &&
				((k < 0) || (k >= Nz) || (i <= 0) || (i >= Nx) || (j <= 0) || (j >= Ny))) {
			}
			else if (!Planewave.isemi) {
				// normal case
				const d_complex_t c0 = Material[MEZ(i,     j,     k    ) + m0].cepsr;
				const d_complex_t c1 = Material[MHY(i,     j,     k    ) + m0].cnur;
				const d_complex_t c2 = Material[MHY(i - 1, j,     k    ) + m0].cnur;
				const d_complex_t c3 = Material[MHX(i,     j,     k    ) + m0].cnur;
				const d_complex_t c4 = Material[MHX(i,     j - 1, 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 ez = planewave(2, Xn[i], Yn[j], Zc[k], ifreq);
				d_complex_t hx = planewave(3, Xn[i], Yn[j], Zc[k], ifreq);
				d_complex_t hy = planewave(4, Xn[i], Yn[j], Zc[k], ifreq);
				ez = d_mul(d_sub(c0, nur), ez);
				hy = d_mul(hy, d_mul(kf, d_mul(d_complex(0, KRXn[i]), d_sub(c1, c2))));
				hx = d_mul(hx, d_mul(kf, d_mul(d_complex(0, KRYn[j]), d_sub(c3, c4))));
				b[n] = d_add(ez, d_sub(hy, hx));
			}
			else {
				// semi-space
				const d_complex_t c0 = Material[MEZ(i,     j,     k    ) + m0].cepsr;
				const d_complex_t c1 = Material[MHY(i,     j,     k    ) + m0].cnur;
				const d_complex_t c2 = Material[MHY(i - 1, j,     k    ) + m0].cnur;
				const d_complex_t c3 = Material[MHX(i,     j,     k    ) + m0].cnur;
				const d_complex_t c4 = Material[MHX(i,     j - 1, k    ) + m0].cnur;

				d_complex_t c[NNZ];
				c[ 0] = d_neg(c0);
				c[ 1] = d_rmul(- KRXp[i],           c1);
				c[ 2] = d_rmul(- KRXm[i],           c2);
				c[ 3] = d_rmul(- KRYp[j],           c3);
				c[ 4] = d_rmul(- KRYm[j],           c4);
				c[ 5] = d_rmul(+ KRZc[k] * KRXn[i], c1);
				c[ 6] = d_neg(c[ 5]);
				c[ 7] = d_rmul(- KRZc[k] * KRXn[i], c2);
				c[ 8] = d_neg(c[ 7]);
				c[ 9] = d_rmul(+ KRZc[k] * KRYn[j], c3);
				c[10] = d_neg(c[ 9]);
				c[11] = d_rmul(- KRZc[k] * KRYn[j], 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 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]);
				const double ync = Yn[j];
				const double ynm = (j >  0) ? Yn[j - 1] : (2 * Yn[ 0] - Yn[     1]);
				const double ynp = (j < Ny) ? Yn[j + 1] : (2 * Yn[Ny] - Yn[Ny - 1]);
				const double ycm = (j >  0) ? Yc[j - 1] : (2 * Yn[ 0] - Yc[     0]);
				const double ycp = (j < Ny) ? Yc[j    ] : (2 * Yn[Ny] - Yc[Ny - 1]);

				d_complex_t e[NNZ];
				e[ 0] = planewave(2, xnc, ync, Zc[k    ], ifreq);
				e[ 1] = planewave(2, xnp, ync, Zc[k    ], ifreq);
				e[ 2] = planewave(2, xnm, ync, Zc[k    ], ifreq);
				e[ 3] = planewave(2, xnc, ynp, Zc[k    ], ifreq);
				e[ 4] = planewave(2, xnc, ynm, Zc[k    ], ifreq);
				e[ 5] = planewave(0, xcp, ync, Zn[k + 1], ifreq);
				e[ 6] = planewave(0, xcp, ync, Zn[k    ], ifreq);
				e[ 7] = planewave(0, xcm, ync, Zn[k + 1], ifreq);
				e[ 8] = planewave(0, xcm, ync, Zn[k    ], ifreq);
				e[ 9] = planewave(1, xnc, ycp, Zn[k + 1], ifreq);
				e[10] = planewave(1, xnc, ycp, Zn[k    ], ifreq);
				e[11] = planewave(1, xnc, ycm, Zn[k + 1], ifreq);
				e[12] = planewave(1, xnc, ycm, Zn[k    ], 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);
			}
		}
	}
	}
	}
}
