/*
planewaveEx.c
*/

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

void planewaveEx(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 = i0Ex; i <= i1Ex; i++) {
	for (int j = j0Ex; j <= j1Ex; j++) {
	for (int k = k0Ex; k <= k1Ex; k++) {
		const int64_t n = NEX(i, j, k);
		if (MEX(i, j, k) == PEC) {
			b[n] = d_neg(planewave(0, Xc[i], Yn[j], Zn[k], ifreq));
		}
		else {
			// exclude Mur and PML points, cf. matrixEx.c
			if ((Abc == MUR) && (
				(!Pbc.ypbc && ((j == 0) || (j == Ny))) ||
				(!Pbc.zpbc && ((k == 0) || (k == Nz))))) {
			}
			else if ((Abc == PML) &&
				((i < 0) || (i >= Nx) || (j <= 0) || (j >= Ny) || (k <= 0) || (k >= Nz))) {
			}
			else if (!Planewave.isemi) {
				// normal case
				const d_complex_t c0 = Material[MEX(i,     j,     k    ) + m0].cepsr;
				const d_complex_t c1 = Material[MHZ(i,     j,     k    ) + m0].cnur;
				const d_complex_t c2 = Material[MHZ(i,     j - 1, k    ) + m0].cnur;
				const d_complex_t c3 = Material[MHY(i,     j,     k    ) + m0].cnur;
				const d_complex_t c4 = Material[MHY(i,     j,     k - 1) + 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 ex = planewave(0, Xc[i], Yn[j], Zn[k], ifreq);
				d_complex_t hy = planewave(4, Xc[i], Yn[j], Zn[k], ifreq);
				d_complex_t hz = planewave(5, Xc[i], Yn[j], Zn[k], ifreq);
				ex = d_mul(d_sub(c0, nur), ex);
				hz = d_mul(hz, d_mul(kf, d_mul(d_complex(0, KRYn[j]), d_sub(c1, c2))));
				hy = d_mul(hy, d_mul(kf, d_mul(d_complex(0, KRZn[k]), d_sub(c3, c4))));
				b[n] = d_add(ex, d_sub(hz, hy));
			}
			else {
				// semi-space
				const d_complex_t c0 = Material[MEX(i,     j,     k    ) + m0].cepsr;
				const d_complex_t c1 = Material[MHZ(i,     j,     k    ) + m0].cnur;
				const d_complex_t c2 = Material[MHZ(i,     j - 1, k    ) + m0].cnur;
				const d_complex_t c3 = Material[MHY(i,     j,     k    ) + m0].cnur;
				const d_complex_t c4 = Material[MHY(i,     j,     k - 1) + m0].cnur;

				d_complex_t c[NNZ];
				c[ 0] = d_neg(c0);
				c[ 1] = d_rmul(- KRYp[j],           c1);
				c[ 2] = d_rmul(- KRYm[j],           c2);
				c[ 3] = d_rmul(- KRZp[k],           c3);
				c[ 4] = d_rmul(- KRZm[k],           c4);
				c[ 5] = d_rmul(+ KRXc[i] * KRYn[j], c1);
				c[ 6] = d_neg(c[ 5]);
				c[ 7] = d_rmul(- KRXc[i] * KRYn[j], c2);
				c[ 8] = d_neg(c[ 7]);
				c[ 9] = d_rmul(+ KRXc[i] * KRZn[k], c3);
				c[10] = d_neg(c[ 9]);
				c[11] = d_rmul(- KRXc[i] * KRZn[k], 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 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]);
				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]);

				d_complex_t e[NNZ];
				e[ 0] = planewave(0, Xc[i    ], ync, znc, ifreq);
				e[ 1] = planewave(0, Xc[i    ], ynp, znc, ifreq);
				e[ 2] = planewave(0, Xc[i    ], ynm, znc, ifreq);
				e[ 3] = planewave(0, Xc[i    ], ync, znp, ifreq);
				e[ 4] = planewave(0, Xc[i    ], ync, znm, ifreq);
				e[ 5] = planewave(1, Xn[i + 1], ycp, znc, ifreq);
				e[ 6] = planewave(1, Xn[i    ], ycp, znc, ifreq);
				e[ 7] = planewave(1, Xn[i + 1], ycm, znc, ifreq);
				e[ 8] = planewave(1, Xn[i    ], ycm, znc, ifreq);
				e[ 9] = planewave(2, Xn[i + 1], ync, zcp, ifreq);
				e[10] = planewave(2, Xn[i    ], ync, zcp, ifreq);
				e[11] = planewave(2, Xn[i + 1], ync, zcm, ifreq);
				e[12] = planewave(2, Xn[i    ], ync, zcm, 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);
			}
		}
	}
	}
	}
}
