#include "common.h"

void phy_sim_destroy(phy_sim *sim) {
	size_t i;

	SDL_DestroyMutex(sim->rotate_lock);

	for (i = 0; i < 3; i++) {
		free(sim->fields[i].E);
		free(sim->fields[i].H);
	}
}

int phy_sim_create(phy_sim *sim) {
	size_t i;
	phy_field_info *fi;

	memset(sim, 0, sizeof(phy_sim));

	fi = &sim->field_info;

	fi->width = 30;
	fi->height = 30;
	fi->depth = 30;
	fi->spacing = 1.5e10 / max3(fi->width, fi->height, fi->depth);

	fi->xstr = 3;
	fi->ystr = fi->xstr * fi->width;
	fi->zstr = fi->ystr * fi->height;
	fi->size = fi->zstr * fi->depth;

	fi->xstr1 = 1;
	fi->ystr1 = fi->xstr1 * fi->width;
	fi->zstr1 = fi->ystr1 * fi->height;
	fi->size1 = fi->zstr1 * fi->depth;

	for (i = 0; i < 3; i++) { 
		sim->fields[i].E = calloc(fi->size, sizeof(float));
		sim->fields[i].H = calloc(fi->size, sizeof(float));
		if (!sim->fields[i].E || !sim->fields[i].H) {
			con_printf("phy_sim_create: out of memory\n");
			goto error;
		}
	}

	sim->field_eps = calloc(fi->size1, sizeof(float));
	sim->field_mu = calloc(fi->size1, sizeof(float));

	phy_sim_compute_const_fields(sim);

	sim->rotate_lock = SDL_CreateMutex();
	assert(sim->rotate_lock);

	return 0;

error:
	phy_sim_destroy(sim);
	return 1;
}

void phy_sim_compute_const_fields(phy_sim *sim)
{
	size_t x, y, z;
	phy_field_info *fi = &sim->field_info;
	size_t border_width;

	border_width = 3;

	for (z = 0; z < fi->depth; z++)
	for (y = 0; y < fi->height; y++)
	for (x = 0; x < fi->width; x++) {
		float *eps, *mu;
		float border_fx, border_fy, border_fz, border_f;

		eps = sim->field_eps + z * fi->zstr1 + y * fi->ystr1 +
		      x * fi->xstr1;
		mu = sim->field_mu + z * fi->zstr1 + y * fi->ystr1 +
		     x * fi->xstr1;

		*eps = 1.129409066837282e+11;
		*mu = -7.957747154594767e+5;

/*
		if (x < border_width) {
			border_fx = (float)x / border_width;
		} else if (x >= fi->width - border_width) {
			border_fx = 1.0f - ((float)x - fi->width + border_width + 1) / border_width;
		} else {
			border_fx = 1.0f;
		}

		if (y < border_width) {
			border_fy = (float)y / border_width;
		} else if (y >= fi->height - border_width) {
			border_fy = 1.0f - ((float)y - fi->height + border_width + 1) / border_width;
		} else {
			border_fy = 1.0f;
		}

		if (z < border_width) {
			border_fz = (float)y / border_width;
		} else if (z >= fi->depth - border_width) {
			border_fz = 1.0f - ((float)z - fi->depth + border_width + 1) / border_width;
		} else {
			border_fz = 1.0f;
		}

		border_f = min3(border_fx, border_fy, border_fz);

		*eps *= border_f;
		*mu *= border_f;*/
	}
}

void phy_sim_step_1(phy_field_info *fi, float dt, float *mul_const,
                    const float *L0, float *L2, const float *R)
{
	size_t x, y, z;

#define AT(F, x, y, z) ((F) + (z) * fi->zstr + (y) * fi->ystr + (x) * fi->xstr)
#define AT1(F, x, y, z) ((F)[(z) * fi->zstr1 + (y) * fi->ystr1 + \
                         (x) * fi->xstr1])

	// X
	for (z = 1; z < fi->depth - 1; z++)
	for (y = 1; y < fi->height - 1; y++)
	for (x = 0; x < fi->width; x++) {
		float dRzdy, dRydz, dLxdt;

		dRzdy = AT(R, x, y + 1, z)[2] - AT(R, x, y - 1, z)[2];
		dRzdy /= 2 * fi->spacing;

		dRydz = AT(R, x, y, z + 1)[1] - AT(R, x, y, z - 1)[1];
		dRydz /= 2 * fi->spacing;

		dLxdt = AT1(mul_const, x, y, z) * (dRzdy - dRydz);

		AT(L2, x, y, z)[0] = AT(L0, x, y, z)[0] + dLxdt * dt;
	}

	// Y
	for (z = 1; z < fi->depth - 1; z++)
	for (y = 0; y < fi->height; y++)
	for (x = 1; x < fi->width - 1; x++) {
		float dRxdz, dRzdx, dLydt;

		dRxdz = AT(R, x, y, z + 1)[0] - AT(R, x, y, z - 1)[0];
		dRxdz /= 2 * fi->spacing;

		dRzdx = AT(R, x + 1, y, z)[2] - AT(R, x - 1, y, z)[2];
		dRzdx /= 2 * fi->spacing;

		dLydt = AT1(mul_const, x, y, z) * (dRxdz - dRzdx);

		AT(L2, x, y, z)[1] = AT(L0, x, y, z)[1] + dLydt * dt;
	}

	// Z
	for (z = 0; z < fi->depth; z++)
	for (y = 1; y < fi->height - 1; y++)
	for (x = 1; x < fi->width - 1; x++) {
		float dRydx, dRxdy, dLzdt;

		dRydx = AT(R, x + 1, y, z)[1] - AT(R, x - 1, y, z)[1];
		dRydx /= 2 * fi->spacing;

		dRxdy = AT(R, x, y + 1, z)[0] - AT(R, x, y - 1, z)[0];
		dRxdy /= 2 * fi->spacing;

		dLzdt = AT1(mul_const, x, y, z) * (dRydx - dRxdy);

		AT(L2, x, y, z)[2] = AT(L0, x, y, z)[2] + dLzdt * dt;
	}
#undef AT
#undef AT1
}

void phy_sim_add_sources(phy_sim *sim)
{
	phy_field_info *fi = &sim->field_info;
	size_t x, y, z;

	for (z = 0; z < fi->depth; z++)
	for (y = 0; y < fi->height; y++)
	for (x = 0; x < fi->width; x++) {
		float *E, *H, fx, fy, fz, chuj;

		fx = (float)x / fi->width * 2 - 1;
		fy = (float)y / fi->height * 2 - 1;
		fz = (float)z / fi->depth * 2 - 1;

		E = sim->fields[2].E + z * fi->zstr + y * fi->ystr + x * fi->xstr;
		H = sim->fields[2].H + z * fi->zstr + y * fi->ystr + x * fi->xstr;

		chuj = (0.3 - sqrt(fx * fx + fy * fy + fz * fz)) / 0.3;
		if (chuj < 0)
			chuj = 0;

		chuj *= sin(sim->time) * 10;
		//chuj *= exp(-10.0f * pow(sim->time - 0.5, 2)) * 2;
		E[0] += chuj;
	}
}

void phy_sim_step(phy_sim *sim, float dt) {
	phy_field_info *fi = &sim->field_info;
	phy_field_em *fields = sim->fields;

	size_t i, size_in_bytes;

	phy_sim_add_sources(sim);

	size_in_bytes = fi->size * sizeof(float);

	SDL_mutexP(sim->rotate_lock);
	for (i = 0; i < 2; i++) {
		// TODO: avoid copying, just switch around pointers
		memcpy(fields[i].E, fields[i + 1].E, size_in_bytes);
		memcpy(fields[i].H, fields[i + 1].H, size_in_bytes);
	}
	SDL_mutexV(sim->rotate_lock);

	phy_sim_step_1(fi, dt, sim->field_eps,
	               fields[0].E, fields[2].E, fields[1].H);
	phy_sim_step_1(fi, dt, sim->field_mu,
	               fields[0].H, fields[2].H, fields[1].E);

	sim->time += dt;
}

/*
void phy_control(phy_command cmd, void *data)

}     */

int phy_thread(phy_sim *sim)
{
	while (1) {
		phy_sim_step(sim, 0.1);
	}

	printf("phy_thread: quitting\n");
}