/*
hyee.c

H on Yee H point
Z0 * H = j / k0 * nu-r * (rot E)
*/

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

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

	double kwave = (2 * PI * Frequency[f]) / C;
	d_complex_t dey = d_rmul(1 / (kwave * DZc[k]), d_sub(eyyee(f, i,    j,     k + 1, exclude),
	                                                     eyyee(f, i,    j,     k,     exclude)));
	d_complex_t dez = d_rmul(1 / (kwave * DYc[j]), d_sub(ezyee(f, i,    j + 1, k,     exclude),
	                                                     ezyee(f, i,    j,     k,     exclude)));
	d_complex_t jnu = d_mul(d_complex(0, 1), Material[MHX(i, j, k) + (f * NMaterial)].cnur);

	return d_mul(jnu, d_sub(dez, dey));
}

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

	double kwave = (2 * PI * Frequency[f]) / C;
	d_complex_t dez = d_rmul(1 / (kwave * DXc[i]), d_sub(ezyee(f, i + 1, j,    k,     exclude),
	                                                     ezyee(f, i,     j,    k,     exclude)));
	d_complex_t dex = d_rmul(1 / (kwave * DZc[k]), d_sub(exyee(f, i,     j,    k + 1, exclude),
	                                                     exyee(f, i,     j,    k,     exclude)));
	d_complex_t jnu = d_mul(d_complex(0, 1), Material[MHY(i, j, k) + (f * NMaterial)].cnur);

	return d_mul(jnu, d_sub(dex, dez));
}

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

	double kwave = (2 * PI * Frequency[f]) / C;
	d_complex_t dex = d_rmul(1 / (kwave * DYc[j]), d_sub(exyee(f, i,     j + 1, k,    exclude),
	                                                     exyee(f, i,     j,     k,    exclude)));
	d_complex_t dey = d_rmul(1 / (kwave * DXc[i]), d_sub(eyyee(f, i + 1, j,     k,    exclude),
	                                                     eyyee(f, i,     j,     k,    exclude)));
	d_complex_t jnu = d_mul(d_complex(0, 1), Material[MHZ(i, j, k) + (f * NMaterial)].cnur);

	return d_mul(jnu, d_sub(dey, dex));
}
