Initial commit.

1 files changed, 244 insertions, 0 deletions

diff --git a/src/physics.c b/src/physics.c
new file mode 100644
index 0000000..be6a426
--- /dev/null
+++ b/src/physics.c
@@ -0,0 +1,244 @@
+#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");
+}
+