/*
hnode.c

Z0 * H on node
*/

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

void hnode(int f, int i, int j, int k, int exclude, d_complex_t h[])
{
	assert((f >= 0) && (f < NFrequency));
	assert((i >= 0) && (i <= Nx) && (j >= 0) && (j <= Ny) && (k >= 0) && (k <= Nz));

	d_complex_t hpp, hpm, hmp, hmm;

	// Z0 * Hx
	hpp = (j < Ny) && (k < Nz) ? hxyee(f, i,     j,     k,     exclude)
	                : (j < Ny) ? hxyee(f, i,     j,     k - 1, exclude)
	                : (k < Nz) ? hxyee(f, i,     j - 1, k,     exclude)
	                           : hxyee(f, i,     j - 1, k - 1, exclude);
	hpm = (j < Ny) && (k >  0) ? hxyee(f, i,     j,     k - 1, exclude)
	                : (j < Ny) ? hxyee(f, i,     j,     k,     exclude)
	                : (k >  0) ? hxyee(f, i,     j - 1, k - 1, exclude)
	                           : hxyee(f, i,     j - 1, k,     exclude);
	hmp = (j >  0) && (k < Nz) ? hxyee(f, i,     j - 1, k,     exclude)
	                : (j >  0) ? hxyee(f, i,     j - 1, k - 1, exclude)
	                : (k < Nz) ? hxyee(f, i,     j,     k,     exclude)
	                           : hxyee(f, i,     j,     k - 1, exclude);
	hmm = (j >  0) && (k >  0) ? hxyee(f, i,     j - 1, k - 1, exclude)
	                : (j >  0) ? hxyee(f, i,     j - 1, k,     exclude)
	                : (k >  0) ? hxyee(f, i,     j,     k - 1, exclude)
	                           : hxyee(f, i,     j,     k,     exclude);
	h[0] = d_rmul(0.25, d_add4(hpp, hpm, hmp, hmm));

	// Z0 * Hy
	hpp = (k < Nz) && (i < Nx) ? hyyee(f, i,     j,     k,     exclude)
	                : (k < Nz) ? hyyee(f, i - 1, j,     k,     exclude)
	                : (i < Nx) ? hyyee(f, i,     j,     k - 1, exclude)
	                           : hyyee(f, i - 1, j,     k - 1, exclude);
	hpm = (k < Nz) && (i >  0) ? hyyee(f, i - 1, j,     k,     exclude)
	                : (k < Nz) ? hyyee(f, i,     j,     k,     exclude)
	                : (i >  0) ? hyyee(f, i - 1, j,     k - 1, exclude)
	                           : hyyee(f, i,     j,     k - 1, exclude);
	hmp = (k >  0) && (i < Nx) ? hyyee(f, i,     j,     k - 1, exclude)
	                : (k >  0) ? hyyee(f, i - 1, j,     k - 1, exclude)
	                : (i < Nx) ? hyyee(f, i,     j,     k,     exclude)
	                           : hyyee(f, i - 1, j,     k,     exclude);
	hmm = (k >  0) && (i >  0) ? hyyee(f, i - 1, j,     k - 1, exclude)
	                : (k >  0) ? hyyee(f, i,     j,     k - 1, exclude)
	                : (i >  0) ? hyyee(f, i - 1, j,     k,     exclude)
	                           : hyyee(f, i,     j,     k,     exclude);
	h[1] = d_rmul(0.25, d_add4(hpp, hpm, hmp, hmm));

	// Z0 * Hz
	hpp = (i < Nx) && (j < Ny) ? hzyee(f, i,     j,     k,     exclude)
	                : (i < Nx) ? hzyee(f, i,     j - 1, k,     exclude)
	                : (j < Ny) ? hzyee(f, i - 1, j,     k,     exclude)
	                           : hzyee(f, i - 1, j - 1, k,     exclude);
	hpm = (i < Nx) && (j >  0) ? hzyee(f, i,     j - 1, k,     exclude)
	                : (i < Nx) ? hzyee(f, i,     j,     k,     exclude)
	                : (j >  0) ? hzyee(f, i - 1, j - 1, k,     exclude)
	                           : hzyee(f, i - 1, j,     k,     exclude);
	hmp = (i >  0) && (j < Ny) ? hzyee(f, i - 1, j,     k,     exclude)
	                : (i >  0) ? hzyee(f, i - 1, j - 1, k,     exclude)
	                : (j < Ny) ? hzyee(f, i,     j,     k,     exclude)
	                           : hzyee(f, i,     j - 1, k,     exclude);
	hmm = (i >  0) && (j >  0) ? hzyee(f, i - 1, j - 1, k,     exclude)
	                : (i >  0) ? hzyee(f, i - 1, j,     k,     exclude)
	                : (j >  0) ? hzyee(f, i,     j - 1, k,     exclude)
	                           : hzyee(f, i,     j,     k,     exclude);
	h[2] = d_rmul(0.25, d_add4(hpp, hpm, hmp, hmm));
}
