/*
rhs.c
*/

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

void rhs(d_complex_t *b, int ifreq)
{
	if (NFeed) {
		// feed
		for (int n = 0; n < NFeed; n++) {
			int i = Feed[n].i;
			int j = Feed[n].j;
			int k = Feed[n].k;
			d_complex_t v = Feed[n].volt;
			if      ((Feed[n].dir == 'X') &&
				(i >= i0Ex) && (i <= i1Ex) &&
				(j >= j0Ex) && (j <= j1Ex) &&
				(k >= k0Ex) && (k <= k1Ex)) {
				b[NEX(i, j, k)] = d_rmul(-1 / DXc[i], v);
			}
			else if ((Feed[n].dir == 'Y') &&
				(i >= i0Ey) && (i <= i1Ey) &&
				(j >= j0Ey) && (j <= j1Ey) &&
				(k >= k0Ey) && (k <= k1Ey)) {
				b[NEY(i, j, k)] = d_rmul(-1 / DYc[j], v);
			}
			else if ((Feed[n].dir == 'Z') &&
				(i >= i0Ez) && (i <= i1Ez) &&
				(j >= j0Ez) && (j <= j1Ez) &&
				(k >= k0Ez) && (k <= k1Ez)) {
				b[NEZ(i, j, k)] = d_rmul(-1 / DZc[k], v);
			}
			//printf("%d %c %d %d %d %f %f\n", n, Feed[n].dir, i, j, k, v.r, v.i);
		}
	}
	else if (IPlanewave) {
		// plane wave
		planewaveEx(b, ifreq);
		planewaveEy(b, ifreq);
		planewaveEz(b, ifreq);
	}
}
