/*
input2.c
*/

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

// index and length
static void getindex(char dir, double x, double y, double z,
	int *i, int *j, int *k, double *dx, double *dy, double *dz)
{
	if      (dir == 'X') {
		*i = nearest(x, 0, Nx - 1, Xc);
		*j = nearest(y, 0, Ny,     Yn);
		*k = nearest(z, 0, Nz,     Zn);
		if ((*j <= 0) || (*j >= Ny) || (*k <= 0) || (*k >= Nz)) {
			printf("%s\n", "*** getindex : invalid position data");
			exit(1);
		}
		*dx = Xn[*i + 1] - Xn[*i    ];
		*dy = Yc[*j    ] - Yc[*j - 1];
		*dz = Zc[*k    ] - Zc[*k - 1];
	}
	else if (dir == 'Y') {
		*i = nearest(x, 0, Nx,     Xn);
		*j = nearest(y, 0, Ny - 1, Yc);
		*k = nearest(z, 0, Nz,     Zn);
		if ((*k <= 0) || (*k >= Nz) || (*i <= 0) || (*i >= Nx)) {
			printf("%s\n", "*** getindex : invalid position data");
			exit(1);
		}
		*dx = Xc[*i    ] - Xc[*i - 1];
		*dy = Yn[*j + 1] - Yn[*j    ];
		*dz = Zc[*k    ] - Zc[*k - 1];
	}
	else if (dir == 'Z') {
		*i = nearest(x, 0, Nx,     Xn);
		*j = nearest(y, 0, Ny,     Yn);
		*k = nearest(z, 0, Nz - 1, Zc);
		if ((*i <= 0) || (*i >= Nx) || (*j <= 0) || (*j >= Ny)) {
			printf("%s\n", "*** getindex : invalid position data");
			exit(1);
		}
		*dx = Xc[*i    ] - Xc[*i - 1];
		*dy = Yc[*j    ] - Yc[*j - 1];
		*dz = Zn[*k + 1] - Zn[*k    ];
	}
}


// number of cells
void setup_cells(int nxr, int nyr, int nzr, int *dxr, int *dyr, int *dzr)
{
	int xsum = 0;
	for (int i = 0; i < nxr; i++) {
		xsum += dxr[i];
	}
	Nx = xsum;

	int ysum = 0;
	for (int j = 0; j < nyr; j++) {
		ysum += dyr[j];
	}
	Ny = ysum;

	int zsum = 0;
	for (int k = 0; k < nzr; k++) {
		zsum += dzr[k];
	}
	Nz = zsum;
}


// node
void setup_node(int nxr, int nyr, int nzr, double *xr, double *yr, double *zr, int *dxr, int *dyr, int *dzr)
{
	if ((nxr <= 0) || (nyr <= 0) || (nzr <= 0)) return;

	Xn = (double *)malloc((Nx + 1) * sizeof(double));
	Yn = (double *)malloc((Ny + 1) * sizeof(double));
	Zn = (double *)malloc((Nz + 1) * sizeof(double));

	int xid = 0;
	for (int mx = 0; mx < nxr; mx++) {
		double dx = (xr[mx + 1] - xr[mx]) / dxr[mx];
		for (int i = 0; i < dxr[mx]; i++) {
			Xn[xid++] = xr[mx] + (i * dx);
		}
	}
	Xn[xid] = xr[nxr];

	int yid = 0;
	for (int my = 0; my < nyr; my++) {
		double dy = (yr[my + 1] - yr[my]) / dyr[my];
		for (int j = 0; j < dyr[my]; j++) {
			Yn[yid++] = yr[my] + (j * dy);
		}
	}
	Yn[yid] = yr[nyr];

	int zid = 0;
	for (int mz = 0; mz < nzr; mz++) {
		double dz = (zr[mz + 1] - zr[mz]) / dzr[mz];
		for (int k = 0; k < dzr[mz]; k++) {
			Zn[zid++] = zr[mz] + (k * dz);
		}
	}
	Zn[zid] = zr[nzr];
/*
	// debug
	for (int i = 0; i <= Nx; i++) {
		printf("Xn[%d]=%.5f\n", i, Xn[i] * 1e3);
	}
	for (int j = 0; j <= Ny; j++) {
		printf("Yn[%d]=%.5f\n", j, Yn[j] * 1e3);
	}
	for (int k = 0; k <= Nz; k++) {
		printf("Zn[%d]=%.5f\n", k, Zn[k] * 1e3);
	}
*/
}


// cell center
void setup_center(void)
{
	Xc = (double *)malloc(Nx * sizeof(double));
	Yc = (double *)malloc(Ny * sizeof(double));
	Zc = (double *)malloc(Nz * sizeof(double));

	for (int i = 0; i < Nx; i++) {
		Xc[i] = (Xn[i] + Xn[i + 1]) / 2;
	}
	for (int j = 0; j < Ny; j++) {
		Yc[j] = (Yn[j] + Yn[j + 1]) / 2;
	}
	for (int k = 0; k < Nz; k++) {
		Zc[k] = (Zn[k] + Zn[k + 1]) / 2;
	}
/*
	// debug
	for (int i = 0; i < Nx; i++) {
		printf("Xc[%d]=%.5f\n", i, Xc[i] * 1e3);
	}
	for (int j = 0; j < Ny; j++) {
		printf("Yc[%d]=%.5f\n", j, Yc[j] * 1e3);
	}
	for (int k = 0; k < Nz; k++) {
		printf("Zc[%d]=%.5f\n", k, Zc[k] * 1e3);
	}
*/
}


// fit geometry without thickness to the nearest node
void fitgeometry(void)
{
	if ((Nx <= 0) || (Ny <= 0) || (Nz <= 0)) return;

	double d0 = EPS * (
		fabs(Xn[Nx] - Xn[0]) +
		fabs(Yn[Ny] - Yn[0]) +
		fabs(Zn[Nz] - Zn[0]));

	for (int n = 0; n < NGeometry; n++) {
		int    shape = Geometry[n].shape;
		double *g = Geometry[n].g;

		if      ((shape == 1) || (shape == 2)) {
			if (fabs(g[0] - g[1]) < d0) {
				int i = nearest(g[0], 0, Nx, Xn);
				g[0] = g[1] = Xn[i];
			}
			if (fabs(g[2] - g[3]) < d0) {
				int j = nearest(g[2], 0, Ny, Yn);
				g[2] = g[3] = Yn[j];
			}
			if (fabs(g[4] - g[5]) < d0) {
				int k = nearest(g[4], 0, Nz, Zn);
				g[4] = g[5] = Zn[k];
			}
		}
		else if ((shape == 11) || (shape == 64) || (shape == 65)) {
			if (fabs(g[0] - g[1]) < d0) {
				int i = nearest(g[0], 0, Nx, Xn);
				g[0] = g[1] = Xn[i];
			}
		}
		else if ((shape == 12) || (shape == 61) || (shape == 66)) {
			if (fabs(g[2] - g[3]) < d0) {
				int j = nearest(g[2], 0, Ny, Yn);
				g[2] = g[3] = Yn[j];
			}
		}
		else if ((shape == 13) || (shape == 62) || (shape == 63)) {
			if (fabs(g[4] - g[5]) < d0) {
				int k = nearest(g[4], 0, Nz, Zn);
				g[4] = g[5] = Zn[k];
			}
		}
	}
}


// plane wave
void setup_planewave(void)
{
	const double cost = cos(Planewave.theta * DTOR);
	const double sint = sin(Planewave.theta * DTOR);
	const double cosp = cos(Planewave.phi   * DTOR);
	const double sinp = sin(Planewave.phi   * DTOR);

	Planewave.r1[0] = + sint * cosp;
	Planewave.r1[1] = + sint * sinp;
	Planewave.r1[2] = + cost;
	Planewave.t1[0] = + cost * cosp;
	Planewave.t1[1] = + cost * sinp;
	Planewave.t1[2] = - sint;
	Planewave.p1[0] = - sinp;
	Planewave.p1[1] = + cosp;
	Planewave.p1[2] = 0;
}


// feed
void setup_feed(const double *x, const double *y, const double *z)
{
	for (int n = 0; n < NFeed; n++) {
		getindex(Feed[n].dir, x[n], y[n], z[n],
			&Feed[n].i, &Feed[n].j, &Feed[n].k, &Feed[n].dx, &Feed[n].dy, &Feed[n].dz);
	}
}


// point
void setup_point(const double *x, const double *y, const double *z, const char strprop[])
{
	for (int n = 0; n < NPoint; n++) {
		double dummy;
		getindex(Point[n].dir, x[n], y[n], z[n],
			&Point[n].i, &Point[n].j, &Point[n].k, &dummy, &dummy, &dummy);
	}

	// add 1+, 1- points
	Point = (point_t *)realloc(Point, (NPoint + 2) * sizeof(point_t));
	Point[NPoint].dir = Point[NPoint + 1].dir = Point[0].dir;
	if      (!strcmp(strprop, "+X") || !strcmp(strprop, "-X")) {
		Point[NPoint].i     = Point[0].i + (!strcmp(strprop, "+X") ? +1 : -1);
		Point[NPoint + 1].i = Point[0].i + (!strcmp(strprop, "+X") ? -1 : +1);
		Point[NPoint].j = Point[NPoint + 1].j = Point[0].j;
		Point[NPoint].k = Point[NPoint + 1].k = Point[0].k;
	}
	else if (!strcmp(strprop, "+Y") || !strcmp(strprop, "-Y")) {
		Point[NPoint].j     = Point[0].j + (!strcmp(strprop, "+Y") ? +1 : -1);
		Point[NPoint + 1].j = Point[0].j + (!strcmp(strprop, "+Y") ? -1 : +1);
		Point[NPoint].k = Point[NPoint + 1].k = Point[0].k;
		Point[NPoint].i = Point[NPoint + 1].i = Point[0].i;
	}
	else if (!strcmp(strprop, "+Z") || !strcmp(strprop, "-Z")) {
		Point[NPoint].k     = Point[0].k + (!strcmp(strprop, "+Z") ? +1 : -1);
		Point[NPoint + 1].k = Point[0].k + (!strcmp(strprop, "+Z") ? -1 : +1);
		Point[NPoint].i = Point[NPoint + 1].i = Point[0].i;
		Point[NPoint].j = Point[NPoint + 1].j = Point[0].j;
	}
/*
	for (int n = 0; n < NPoint + 2; n++) {
		printf("%d %c %d %d %d\n", n, Point[n].dir, Point[n].i, Point[n].j, Point[n].k);
	}
*/
}


// load
void setup_load(material_data_t *material_data, int array_inc)
{
	for (int n = 0; n < NMaterial; n++) {
		const int type = material_data[n].type;
		const char dir = material_data[n].dirload;
		if ((type == LOAD_R) || (type == LOAD_C) || (type == LOAD_L)) {
			if (NGeometry % array_inc == 0) {
				Geometry = (geometry_t *)realloc(Geometry, (NGeometry + array_inc) * sizeof(geometry_t));
			}
			int i = 0, j = 0, k = 0;
			double dummy;
			getindex(dir, material_data[n].posload[0], material_data[n].posload[1], material_data[n].posload[2], &i, &j, &k, &dummy, &dummy, &dummy);
			geometry_t *ptr = &Geometry[NGeometry];
			ptr->m = n;
			ptr->shape = 1;
			if      (dir == 'X') {
				ptr->g[0] = Xn[i    ];
				ptr->g[1] = Xn[i + 1];
				ptr->g[2] =
				ptr->g[3] = Yn[j];
				ptr->g[4] =
				ptr->g[5] = Zn[k];
				material_data[n].dsdl = (Yc[j] - Yc[j - 1]) * (Zc[k] - Zc[k - 1]) / (Xn[i + 1] - Xn[i]);
			}
			else if (dir == 'Y') {
				ptr->g[2] = Yn[j    ];
				ptr->g[3] = Yn[j + 1];
				ptr->g[4] =
				ptr->g[5] = Zn[k];
				ptr->g[0] =
				ptr->g[1] = Xn[i];
				material_data[n].dsdl = (Zc[k] - Zc[k - 1]) * (Xc[i] - Xc[i - 1]) / (Yn[j + 1] - Yn[j]);
			}
			else if (dir == 'Z') {
				ptr->g[4] = Zn[k    ];
				ptr->g[5] = Zn[k + 1];
				ptr->g[0] =
				ptr->g[1] = Xn[i];
				ptr->g[2] =
				ptr->g[3] = Yn[j];
				material_data[n].dsdl = (Xc[i] - Xc[i - 1]) * (Yc[j] - Yc[j - 1]) / (Zn[k + 1] - Zn[k]);
			}
			NGeometry++;
		}
	}
/*
	// debug
	for (int n = 0; n < NGeometry; n++) {
		printf("geometry : %d %d %d %.3e %.3e %.3e %.3e %.3e %.3e\n", n, Geometry[n].m, Geometry[n].shape, Geometry[n].g[0], Geometry[n].g[1], Geometry[n].g[2], Geometry[n].g[3], Geometry[n].g[4], Geometry[n].g[5]);
	}
*/
}


// material
void setup_material(material_data_t *material_data)
{
	Material = (material_t *)malloc(NMaterial * NFrequency * sizeof(material_t));

	for (int m = 0; m < NMaterial; m++) {
		const int type = material_data[m].type;
		for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
			const double freq = Frequency[ifreq];
			const double omega = 2 * PI * freq;
			d_complex_t ce = d_complex(1, 0);
			d_complex_t cm = d_complex(1, 0);
			if      (type == 1) {
				ce = material_data[m].cepsr;
				cm = material_data[m].cmur;
			}
			else if (type == 2) {
				// dispersion file
				double p[4] = {0, 0, 0, 0};
				interpolation(freq, material_data[m].ntable, 5, material_data[m].table, p);
				ce = d_complex(p[0], -p[1]);
				cm = d_complex(p[2], -p[3]);
			}
			else if (type == LOAD_R) {
				// R
				double r = material_data[m].rcl;
				ce = d_complex(1, -1 / (r * omega * EPS0 * material_data[m].dsdl));
				cm = d_complex(1, 0);
			}
			else if (type == LOAD_C) {
				// C
				double c = material_data[m].rcl;
				ce = d_complex(c / (EPS0 * material_data[m].dsdl), 0);
				cm = d_complex(1, 0);
			}
			else if (type == LOAD_L) {
				// L
				double l = material_data[m].rcl;
				ce = d_complex(-1 / (l * omega * omega * EPS0 * material_data[m].dsdl), 0);
				cm = d_complex(1, 0);
			}
			int id = m + (ifreq * NMaterial);
			Material[id].type  = type;
			Material[id].cepsr = ce;
			Material[id].cmur  = cm;
			Material[id].cnur  = d_inv(cm);
			Material[id].cv    = d_sqrt(d_mul(ce, cm));
			//printf("%d %d %f %f %f %f %f %f %f %f\n", m, ifreq, Material[id].cepsr.r, Material[id].cepsr.i, Material[id].cmur.r, Material[id].cmur.i, Material[id].cnur.r, Material[id].cnur.i, Material[id].cv.r, Material[id].cv.i);
		}
	}
}


// setup geometry lines
void setup_geomline3d(void)
{
	if (!NGeometry) return;

	// alloc
	int *shape     =         (int *)malloc(NGeometry * sizeof(int));
	int *mid       =         (int *)malloc(NGeometry * sizeof(int));
	double (*g)[8] = (double (*)[8])malloc(NGeometry * 8 * sizeof(double));
	for (int n = 0; n < NGeometry; n++) {
		shape[n] = Geometry[n].shape;
		mid[n] = Geometry[n].m;
		memcpy(g[n], Geometry[n].g, 8 * sizeof(double));
	}

	// number of lines
	NGline = geomlines(0, NGeometry, shape, NULL, NULL, NULL, NULL, 0);

	// alloc
	Gline  = (double (*)[2][3])malloc(NGline * 2 * 3 * sizeof(double));
	MGline =            (int *)malloc(NGline * sizeof(int));
	memset(Gline, 0, NGline * 2 * 3 * sizeof(double));
	memset(MGline, 0, NGline * sizeof(int));

	// set geometry
	const double eps = EPS * (fabs(Xn[Nx] - Xn[0]) + fabs(Yn[Ny] - Yn[0]) + fabs(Zn[Nz] - Zn[0]));
	geomlines(1, NGeometry, shape, mid, g, Gline, MGline, eps);

	// free
	free(shape);
	free(mid);
	free(g);
}
