/*
setup2.c

frequency dependent
*/

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

// k * cellsize
static void setup_cell_factor(int ifreq)
{
	const double kwave = (2 * PI * Frequency[ifreq]) / C;

	for (int i = 0; i <= Nx; i++) {
		KRXn[i] = 1 / (kwave * DXn[i]);
	}
	for (int j = 0; j <= Ny; j++) {
		KRYn[j] = 1 / (kwave * DYn[j]);
	}
	for (int k = 0; k <= Nz; k++) {
		KRZn[k] = 1 / (kwave * DZn[k]);
	}

	for (int i = 0; i < Nx; i++) {
		KRXc[i] = 1 / (kwave * DXc[i]);
	}
	for (int j = 0; j < Ny; j++) {
		KRYc[j] = 1 / (kwave * DYc[j]);
	}
	for (int k = 0; k < Nz; k++) {
		KRZc[k] = 1 / (kwave * DZc[k]);
	}

	for (int i = 0; i <= Nx; i++) {
		KDXn[i] = kwave * DXn[i];
	}
	for (int j = 0; j <= Ny; j++) {
		KDYn[j] = kwave * DYn[j];
	}
	for (int k = 0; k <= Nz; k++) {
		KDZn[k] = kwave * DZn[k];
	}

	for (int i = 0; i < Nx; i++) {
		KDXc[i] = kwave * DXc[i];
	}
	for (int j = 0; j < Ny; j++) {
		KDYc[j] = kwave * DYc[j];
	}
	for (int k = 0; k < Nz; k++) {
		KDZc[k] = kwave * DZc[k];
	}

	for (int i = 0; i <= Nx; i++) {
		KRXp[i] = 1 / (kwave * kwave * DXn[i] * DXc[(i < Nx) ? i : i - 1]);
		KRXm[i] = 1 / (kwave * kwave * DXn[i] * DXc[(i >  0) ? i - 1 : i]);
	}
	for (int j = 0; j <= Ny; j++) {
		KRYp[j] = 1 / (kwave * kwave * DYn[j] * DYc[(j < Ny) ? j : j - 1]);
		KRYm[j] = 1 / (kwave * kwave * DYn[j] * DYc[(j >  0) ? j - 1 : j]);
	}
	for (int k = 0; k <= Nz; k++) {
		KRZp[k] = 1 / (kwave * kwave * DZn[k] * DZc[(k < Nz) ? k : k - 1]);
		KRZm[k] = 1 / (kwave * kwave * DZn[k] * DZc[(k >  0) ? k - 1 : k]);
	}
}

// PBC
static void setup_pbc(int ifreq)
{
	const double kwave = (2 * PI * Frequency[ifreq]) / C;

	if (NFeed) {
		// feed
		if (NFeed > 1) {
			double x[2], y[2], z[2];
			for (int m = 0; m < 2; m++) {
				x[m] = (Feed[m].dir == 'X') ? Xc[Feed[m].i] : Xn[Feed[m].i];
				y[m] = (Feed[m].dir == 'Y') ? Yc[Feed[m].j] : Yn[Feed[m].j];
				z[m] = (Feed[m].dir == 'Z') ? Zc[Feed[m].k] : Zn[Feed[m].k];
			}
			double dphase = (d_deg(Feed[1].volt) - d_deg(Feed[0].volt)) * DTOR;
			if (dphase < -PI) dphase += 2 * PI;
			if (dphase > +PI) dphase -= 2 * PI;
			const double px = (fabs(x[1] - x[0]) > EPS) ? dphase * (Xn[Nx] - Xn[0]) / (x[1] - x[0]) : 0;
			const double py = (fabs(y[1] - y[0]) > EPS) ? dphase * (Yn[Ny] - Yn[0]) / (y[1] - y[0]) : 0;
			const double pz = (fabs(z[1] - z[0]) > EPS) ? dphase * (Zn[Nz] - Zn[0]) / (z[1] - z[0]) : 0;
			Pbc.fx = d_exp(px);
			Pbc.fy = d_exp(py);
			Pbc.fz = d_exp(pz);
		}
		else {
			Pbc.fx = d_complex(1, 0);
			Pbc.fy = d_complex(1, 0);
			Pbc.fz = d_complex(1, 0);
		}
	}
	else if (IPlanewave) {
		// planewave
		d_complex_t k = d_complex(0, 0);
		if (!Planewave.isemi && !Planewave.ifull) {
			// normal case
			k = d_complex(kwave, 0);
		}
		else if (Planewave.ifull) {
			// full space
			const int m0 = ifreq * NMaterial;
			const d_complex_t epsr0 = Material[Planewave.mfull + m0].cepsr;
			const d_complex_t mur0  = Material[Planewave.mfull + m0].cmur;
			k = d_rmul(kwave, d_sqrt(d_mul(epsr0, mur0)));
		}
		else if (Planewave.isemi) {
			// semi space
			// 1 = 2 by Snell's law
			const int m0 = ifreq * NMaterial;
			const d_complex_t epsr1 = Material[Planewave.msemi1 + m0].cepsr;
			const d_complex_t mur1  = Material[Planewave.msemi1 + m0].cmur;
			k = d_rmul(kwave, d_sqrt(d_mul(epsr1, mur1)));
		}
		const double dlx = -Planewave.r1[0] * (Xn[Nx] - Xn[0]);
		const double dly = -Planewave.r1[1] * (Yn[Ny] - Yn[0]);
		const double dlz = -Planewave.r1[2] * (Zn[Nz] - Zn[0]);
		//printf("%f %f %f\n", dlx, dly, dlz);
		Pbc.fx = expjkr(d_neg(k), dlx);
		Pbc.fy = expjkr(d_neg(k), dly);
		Pbc.fz = expjkr(d_neg(k), dlz);
		//printf("%f %f %f %f %f %f\n", Pbc.fx.r, Pbc.fx.i, Pbc.fy.r, Pbc.fy.i, Pbc.fz.r, Pbc.fz.i);
	}
}

void setup2(int ifreq)
{
	setup_cell_factor(ifreq);

	if (Pbc.xpbc || Pbc.ypbc || Pbc.zpbc) {
		setup_pbc(ifreq);
	}
}
