/*
setup.c

frequency independent
*/

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

static void setup_cell_size_node(void)
{
	DXn = (double *)malloc((Nx + 1) * sizeof(double));
	DYn = (double *)malloc((Ny + 1) * sizeof(double));
	DZn = (double *)malloc((Nz + 1) * sizeof(double));

	for (int i = 1; i < Nx; i++) {
		DXn[i] = (Xn[i + 1] - Xn[i - 1]) / 2;
	}
	for (int j = 1; j < Ny; j++) {
		DYn[j] = (Yn[j + 1] - Yn[j - 1]) / 2;
	}
	for (int k = 1; k < Nz; k++) {
		DZn[k] = (Zn[k + 1] - Zn[k - 1]) / 2;
	}
	DXn[0]  = Xn[1]  - Xn[0];
	DYn[0]  = Yn[1]  - Yn[0];
	DZn[0]  = Zn[1]  - Zn[0];
	DXn[Nx] = Xn[Nx] - Xn[Nx - 1];
	DYn[Ny] = Yn[Ny] - Yn[Ny - 1];
	DZn[Nz] = Zn[Nz] - Zn[Nz - 1];
/*
	for (int i = 0; i <= Nx; i++) printf("DXn[%d]=%f\n", i, DXn[i]);
	for (int j = 0; j <= Ny; j++) printf("DYn[%d]=%f\n", j, DYn[j]);
	for (int k = 0; k <= Nz; k++) printf("DZn[%d]=%f\n", k, DZn[k]);
*/
}


static void setup_cell_size_center(void)
{
	DXc = (double *)malloc(Nx * sizeof(double));
	DYc = (double *)malloc(Ny * sizeof(double));
	DZc = (double *)malloc(Nz * sizeof(double));

	for (int i = 0; i < Nx; i++) {
		DXc[i] = Xn[i + 1] - Xn[i];
	}
	for (int j = 0; j < Ny; j++) {
		DYc[j] = Yn[j + 1] - Yn[j];
	}
	for (int k = 0; k < Nz; k++) {
		DZc[k] = Zn[k + 1] - Zn[k];
	}
/*
	for (int i = 0; i < Nx; i++) printf("DXc[%d]=%f\n", i, DXc[i]);
	for (int j = 0; j < Ny; j++) printf("DYc[%d]=%f\n", j, DYc[j]);
	for (int k = 0; k < Nz; k++) printf("DZc[%d]=%f\n", k, DZc[k]);
*/
}


static double pml_factor(double f, int n)
{
	double kpml = (Mpml + 1.0) / (2.0 * Lpml) * log(1.0 / R0pml);

	return (f <     EPS) ? kpml * pow(fabs(f    ) / Lpml, Mpml)
	     : (f > n - EPS) ? kpml * pow(fabs(f - n) / Lpml, Mpml) : 0;
}


static void setup_pml_array(void)
{
	NFxpml = 2 * (Nx + 2 * Lpml) + 1;
	NFypml = 2 * (Ny + 2 * Lpml) + 1;
	NFzpml = 2 * (Nz + 2 * Lpml) + 1;

	Fxpml = (double *)malloc(NFxpml * sizeof(double));
	Fypml = (double *)malloc(NFypml * sizeof(double));
	Fzpml = (double *)malloc(NFzpml * sizeof(double));
}


static void setup_pml_factor(void)
{
	for (int i = 0; i < NFxpml; i++) {
		Fxpml[i] = pml_factor(i / 2.0 - Lpml, Nx);
	}
	for (int j = 0; j < NFypml; j++) {
		Fypml[j] = pml_factor(j / 2.0 - Lpml, Ny);
	}
	for (int k = 0; k < NFzpml; k++) {
		Fzpml[k] = pml_factor(k / 2.0 - Lpml, Nz);
	}
/*
	for (int i = 0; i < NFxpml; i++) printf("%d %f %f\n", i, i/2.0-Lpml, Fxpml[i]);
	for (int j = 0; j < NFypml; j++) printf("%d %f %f\n", j, j/2.0-Lpml, Fypml[j]);
	for (int k = 0; k < NFzpml; k++) printf("%d %f %f\n", k, k/2.0-Lpml, Fzpml[k]);
*/
}

/*
// higher material ?
static void higher(int *id0, int id1)
{
	// type=1 -> 1st-frequency only
	const material_t m0 = Material[*id0];
	const material_t m1 = Material[ id1];

	if      ((*id0 == id1) ||
	         (*id0 == PEC) ||
	         (m0.type > 1) ||
	         (m1.type > 2)) {
		;
	}
	else if (( id1 == PEC) ||
	         (m1.type == 2) ||
	         (m0.cepsr.r + m0.cmur.r < m1.cepsr.r + m1.cmur.r)) {
		*id0 = id1;
	}
}
*/

// highest material
// type = 1(normal), 2(file), 3(R), 4(C), 5(L)
static inline int highest(int id0, int id1, int id2, int id3, int id4, int *type)
{
	int ret = id0;

	if      ((id1 == id0) && (id2 == id0) && (id3 == id0) && (id4 == id0)) {
		;
	}
	else if ((id0 == PEC) || (id1 == PEC) || (id2 == PEC) || (id3 == PEC) || (id4 == PEC)) {
		ret = PEC;
	}
	else if (type[id1] == 2) {
		ret = id1;
	}
	else if (type[id2] == 2) {
		ret = id2;
	}
	else if (type[id3] == 2) {
		ret = id3;
	}
	else if (type[id4] == 2) {
		ret = id4;
	}

	return ret;
}


// correct curved surface (E <- H)
static void setup_surface(void)
{
	//const int m0 = (IPlanewave && Planewave.ifull) ? Planewave.mfull : 0;
	//printf("%d %d %d %d %d %d\n", i0Ex, i1Ex, j0Ex, j1Ex, k0Ex, k1Ex);

	// alloc work array
	int *wk = (int *)malloc(NE * sizeof(int));
	memset(wk, 0, NE * sizeof(int));

	int i;

	// Ex
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (    i = i0Ex; i <= i1Ex; i++) {
	for (int j = j0Ex; j <= j1Ex; j++) {
	for (int k = k0Ex; k <= k1Ex; k++) {
		//assert(NEX(i, j, k) < NE);
		wk[NEX(i, j, k)] = Material[MEX(i, j, k)].type;
	}
	}
	}

	// Ey
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (    i = i0Ey; i <= i1Ey; i++) {
	for (int j = j0Ey; j <= j1Ey; j++) {
	for (int k = k0Ey; k <= k1Ey; k++) {
		//assert(NEY(i, j, k) < NE);
		wk[NEY(i, j, k)] = Material[MEY(i, j, k)].type;
	}
	}
	}

	// Ez
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (    i = i0Ez; i <= i1Ez; i++) {
	for (int j = j0Ez; j <= j1Ez; j++) {
	for (int k = k0Ez; k <= k1Ez; k++) {
		//assert(NEZ(i, j, k) < NE);
		wk[NEZ(i, j, k)] = Material[MEZ(i, j, k)].type;
	}
	}
	}

	// Ex
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (    i = i0Ex; i <= i1Ex; i++) {
	for (int j = j0Ex; j <= j1Ex; j++) {
	for (int k = k0Ex; k <= k1Ex; k++) {
/*
		if ((MEX(i, j, k) == m0) && (j > 0) && (j < Ny) && (k > 0) && (k < Nz)) {
			if (MHY(i,     j,     k    ) != m0) MEX(i, j, k) = MHY(i,     j,     k    );
			if (MHY(i,     j,     k - 1) != m0) MEX(i, j, k) = MHY(i,     j,     k - 1);
			if (MHZ(i,     j,     k    ) != m0) MEX(i, j, k) = MHZ(i,     j,     k    );
			if (MHZ(i,     j - 1, k    ) != m0) MEX(i, j, k) = MHZ(i,     j - 1, k    );
		}
*/
		if ((j > 0) && (j < Ny) && (k > 0) && (k < Nz)) {
/*
			higher(&MEX(i, j, k), MHY(i,     j,     k    ));
			higher(&MEX(i, j, k), MHY(i,     j,     k - 1));
			higher(&MEX(i, j, k), MHZ(i,     j,     k    ));
			higher(&MEX(i, j, k), MHZ(i,     j - 1, k    ));
*/
			MEX(i, j, k) = highest(MEX(i, j, k), MHY(i,     j,     k    ),
			                                     MHY(i,     j,     k - 1),
			                                     MHZ(i,     j,     k    ),
			                                     MHZ(i,     j - 1, k    ), wk);
		}

	}
	}
	}

	// Ey
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (    i = i0Ey; i <= i1Ey; i++) {
	for (int j = j0Ey; j <= j1Ey; j++) {
	for (int k = k0Ey; k <= k1Ey; k++) {
/*
		if ((MEY(i, j, k) == m0) && (k > 0) && (k < Nz) && (i > 0) && (i < Nx)) {
			if (MHZ(i,     j,     k    ) != m0) MEY(i, j, k) = MHZ(i,     j,     k    );
			if (MHZ(i - 1, j,     k    ) != m0) MEY(i, j, k) = MHZ(i - 1, j,     k    );
			if (MHX(i,     j,     k    ) != m0) MEY(i, j, k) = MHX(i,     j,     k    );
			if (MHX(i,     j,     k - 1) != m0) MEY(i, j, k) = MHX(i,     j,     k - 1);
		}
*/
		if ((k > 0) && (k < Nz) && (i > 0) && (i < Nx)) {
/*
			higher(&MEY(i, j, k), MHZ(i,     j,     k    ));
			higher(&MEY(i, j, k), MHZ(i - 1, j,     k    ));
			higher(&MEY(i, j, k), MHX(i,     j,     k    ));
			higher(&MEY(i, j, k), MHX(i,     j,     k - 1));
*/
			MEY(i, j, k) = highest(MEY(i, j, k), MHZ(i,     j,     k    ),
			                                     MHZ(i - 1, j,     k    ),
			                                     MHX(i,     j,     k    ),
			                                     MHX(i,     j,     k - 1), wk);
		}

	}
	}
	}

	// Ez
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (    i = i0Ez; i <= i1Ez; i++) {
	for (int j = j0Ez; j <= j1Ez; j++) {
	for (int k = k0Ez; k <= k1Ez; k++) {
/*
		if ((MEZ(i, j, k) == m0) && (i > 0) && (i < Nx) && (j > 0) && (j < Ny)) {
			if (MHX(i,     j,     k    ) != m0) MEZ(i, j, k) = MHX(i,     j,     k    );
			if (MHX(i,     j - 1, k    ) != m0) MEZ(i, j, k) = MHX(i,     j - 1, k    );
			if (MHY(i,     j,     k    ) != m0) MEZ(i, j, k) = MHY(i,     j,     k    );
			if (MHY(i - 1, j,     k    ) != m0) MEZ(i, j, k) = MHY(i - 1, j,     k    );
		}
*/
		if ((i > 0) && (i < Nx) && (j > 0) && (j < Ny)) {
/*
			higher(&MEZ(i, j, k), MHX(i,     j,     k    ));
			higher(&MEZ(i, j, k), MHX(i,     j - 1, k    ));
			higher(&MEZ(i, j, k), MHY(i,     j,     k    ));
			higher(&MEZ(i, j, k), MHY(i - 1, j,     k    ));
*/
			MEZ(i, j, k) = highest(MEZ(i, j, k), MHX(i,     j,     k    ),
			                                     MHX(i,     j - 1, k    ),
			                                     MHY(i,     j,     k    ),
			                                     MHY(i - 1, j,     k    ), wk);
		}

	}
	}
	}

	// free work array
	free(wk);
}


void setup(void)
{
	// cell size (node)
	setup_cell_size_node();

	// cell size (center)
	setup_cell_size_center();

	// PML factor
	setup_pml_array();   // include Mur : SX
	if (Abc == PML) {
		setup_pml_factor();
	}

	// background media
	if (IPlanewave && (Planewave.isemi || Planewave.ifull)) {
		setupBg();
	}

	// correct curved surface
	setup_surface();
}
