/*
cyee.c

current on Yee E point
I = dS * (rot H)
*/

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

d_complex_t cxyee(int f, int i, int j, int k)
{
	assert((f >= 0) && (f < NFrequency));
	assert((i >= 0) && (i < Nx) && (j >= 0) && (j <= Ny) && (k >= 0) && (k <= Nz));

	d_complex_t hzyp = (j < Ny) ? hzyee(f, i,     j,     k,     0) : hzyee(f, i,     j - 1, k,     0);
	d_complex_t hzym = (j >  0) ? hzyee(f, i,     j - 1, k,     0) : hzyee(f, i,     j,     k,     0);
	d_complex_t hyzp = (k < Nz) ? hyyee(f, i,     j,     k,     0) : hyyee(f, i,     j,     k - 1, 0);
	d_complex_t hyzm = (k >  0) ? hyyee(f, i,     j,     k - 1, 0) : hyyee(f, i,     j,     k,     0);

	return d_sub(d_rmul(DZn[k] / ETA0, d_sub(hzyp, hzym)),
	             d_rmul(DYn[j] / ETA0, d_sub(hyzp, hyzm)));
}

d_complex_t cyyee(int f, int i, int j, int k)
{
	assert((f >= 0) && (f < NFrequency));
	assert((j >= 0) && (j < Ny) && (k >= 0) && (k <= Nz) && (i >= 0) && (i <= Nx));

	d_complex_t hxzp = (k < Nz) ? hxyee(f, i,     j,     k,     0) : hxyee(f, i,     j,     k - 1, 0);
	d_complex_t hxzm = (k >  0) ? hxyee(f, i,     j,     k - 1, 0) : hxyee(f, i,     j,     k,     0);
	d_complex_t hzxp = (i < Nx) ? hzyee(f, i,     j,     k,     0) : hzyee(f, i - 1, j,     k,     0);
	d_complex_t hzxm = (i >  0) ? hzyee(f, i - 1, j,     k,     0) : hzyee(f, i,     j,     k,     0);

	return d_sub(d_rmul(DXn[i] / ETA0, d_sub(hxzp, hxzm)),
	             d_rmul(DZn[k] / ETA0, d_sub(hzxp, hzxm)));
}

d_complex_t czyee(int f, int i, int j, int k)
{
	assert((f >= 0) && (f < NFrequency));
	assert((k >= 0) && (k < Nz) && (i >= 0) && (i <= Nx) && (j >= 0) && (j <= Ny));

	d_complex_t hyxp = (i < Nx) ? hyyee(f, i,     j,     k,     0) : hyyee(f, i - 1, j,     k,     0);
	d_complex_t hyxm = (i >  0) ? hyyee(f, i - 1, j,     k,     0) : hyyee(f, i,     j,     k,     0);
	d_complex_t hxyp = (j < Ny) ? hxyee(f, i,     j,     k,     0) : hxyee(f, i,     j - 1, k,     0);
	d_complex_t hxym = (j >  0) ? hxyee(f, i,     j - 1, k,     0) : hxyee(f, i,     j,     k,     0);

	return d_sub(d_rmul(DYn[j] / ETA0, d_sub(hyxp, hyxm)),
	             d_rmul(DXn[i] / ETA0, d_sub(hxyp, hxym)));
}
