/*
setupBg_mpi.c (MPI)

setup background media
*/

#include "oth.h"
//#include "oth_mpi.h"
#include "oth_prototype.h"

void setupBg(void)
{
	assert(IPlanewave);
	assert(Planewave.isemi || Planewave.ifull);

	int64_t *nsemi1 = (int64_t *)calloc(NMaterial, sizeof(int64_t));
	int64_t *nsemi2 = (int64_t *)calloc(NMaterial, sizeof(int64_t));
	int64_t *nfull  = (int64_t *)calloc(NMaterial, sizeof(int64_t));

	const int isemi = Planewave.isemi;
	const double x0 = (isemi == 1) ? Planewave.psemi : 0;
	const double y0 = (isemi == 2) ? Planewave.psemi : 0;
	const double z0 = (isemi == 3) ? Planewave.psemi : 0;

	double ri[3];
	for (int m = 0; m < 3; m++) {
		ri[m] = -Planewave.r1[m];
	}

	// Ex
	for (int i = i0Ex; i <= i1Ex; i++) {
	for (int j = j0Ex; j <= j1Ex; j++) {
	for (int k = k0Ex; k <= k1Ex; k++) {
		const int m = MEX(i, j, k);
		const double x = Xc[i];
		const double y = Yn[j];
		const double z = Zn[k];
		if (Planewave.isemi) {
			if (((isemi == 1) && (((x - x0) * ri[0]) <= 0)) ||
			    ((isemi == 2) && (((y - y0) * ri[1]) <= 0)) ||
			    ((isemi == 3) && (((z - z0) * ri[2]) <= 0))) {
				nsemi1[m]++;
			}
			else {
				nsemi2[m]++;
			}
		}
		else if (Planewave.ifull) {
			nfull[m]++;
		}
	}
	}
	}

	// Ey
	for (int i = i0Ey; i <= i1Ey; i++) {
	for (int j = j0Ey; j <= j1Ey; j++) {
	for (int k = k0Ey; k <= k1Ey; k++) {
		const int m = MEY(i, j, k);
		const double x = Xn[i];
		const double y = Yc[j];
		const double z = Zn[k];
		if (Planewave.isemi) {
			if (((isemi == 1) && (((x - x0) * ri[0]) <= 0)) ||
			    ((isemi == 2) && (((y - y0) * ri[1]) <= 0)) ||
			    ((isemi == 3) && (((z - z0) * ri[2]) <= 0))) {
				nsemi1[m]++;
			}
			else {
				nsemi2[m]++;
			}
		}
		else if (Planewave.ifull) {
			nfull[m]++;
		}
	}
	}
	}

	// Ez
	for (int i = i0Ez; i <= i1Ez; i++) {
	for (int j = j0Ez; j <= j1Ez; j++) {
	for (int k = k0Ez; k <= k1Ez; k++) {
		const int m = MEZ(i, j, k);
		const double x = Xn[i];
		const double y = Yn[j];
		const double z = Zc[k];
		if (Planewave.isemi) {
			if (((isemi == 1) && (((x - x0) * ri[0]) <= 0)) ||
			    ((isemi == 2) && (((y - y0) * ri[1]) <= 0)) ||
			    ((isemi == 3) && (((z - z0) * ri[2]) <= 0))) {
				nsemi1[m]++;
			}
			else {
				nsemi2[m]++;
			}
		}
		else if (Planewave.ifull) {
			nfull[m]++;
		}
	}
	}
	}

	// sum (MPI)
	for (int m = 0; m < NMaterial; m++) {
		//printf("%d %d %d %zd %zd %zd\n", commSize, commRank, m, nsemi1[m], nsemi2[m], nfull[m]); fflush(stdout);
		nsemi1[m] = comm_isum((int)nsemi1[m]);
		nsemi2[m] = comm_isum((int)nsemi2[m]);
		nfull[m]  = comm_isum((int)nfull[m]);
		//printf("%d %d %d %zd %zd %zd\n", commSize, commRank, m, nsemi1[m], nsemi2[m], nfull[m]); fflush(stdout);
	}

	// find max
	int msemi1 = 0;
	int msemi2 = 0;
	int mfull  = 0;
	int64_t maxsemi1 = 0;
	int64_t maxsemi2 = 0;
	int64_t maxfull  = 0;
	for (int m = 0; m < NMaterial; m++) {
		//printf("%d %zd %zd %zd\n", m, nsemi1[m], nsemi2[m], nfull[m]);
		if (nsemi1[m] > maxsemi1) {
			msemi1 = m;
			maxsemi1 = nsemi1[m];
		}
		if (nsemi2[m] > maxsemi2) {
			msemi2 = m;
			maxsemi2 = nsemi2[m];
		}
		if (nfull[m] > maxfull) {
			mfull = m;
			maxfull = nfull[m];
		}
	}
	Planewave.msemi1 = msemi1;
	Planewave.msemi2 = msemi2;
	Planewave.mfull  = mfull;
	//printf("%d %d %d %d %d\n", commSize, commRank, Planewave.msemi1, Planewave.msemi2, Planewave.mfull); fflush(stdout);

	free(nsemi1);
	free(nsemi2);
	free(nfull);
}
