Yee grid.

3 files changed, 179 insertions, 122 deletions

diff --git a/src/physics.c b/src/physics.c
index 49c6552..3981258 100644
--- a/src/physics.c
+++ b/src/physics.c
@@ -2,7 +2,8 @@
 #include "itc.h"
 #include "physics.h"
 
-void phy_sim_destroy(phy_sim *sim) {
+void phy_sim_destroy(phy_sim *sim)
+{
 	size_t i;
 
 	SDL_DestroyMutex(sim->rotate_lock);
@@ -15,7 +16,28 @@ void phy_sim_destroy(phy_sim *sim) {
 	}
 }
 
-int phy_sim_create(phy_sim *sim) {
+void phy_sim_reset(phy_sim *sim)
+{
+	size_t i;
+	phy_field_info *fi;
+
+	fi = &sim->field_info;
+
+	for (i = 0; i < 3; i++) {
+		memset(sim->fields[i].E, 0, fi->size * sizeof(float));
+		memset(sim->fields[i].H, 0, fi->size * sizeof(float));
+	}
+
+	memset(sim->field_eps, 0, fi->size1 * sizeof(float));
+	memset(sim->field_eps, 0, fi->size1 * sizeof(float));
+	phy_sim_compute_const_fields(sim);
+
+	sim->time = 0;
+	sim->running = false;
+}
+
+int phy_sim_create(phy_sim *sim)
+{
 	size_t i;
 	phy_field_info *fi;
 
@@ -23,10 +45,11 @@ int phy_sim_create(phy_sim *sim) {
 
 	fi = &sim->field_info;
 
-	fi->width = 50;
-	fi->height = 50;
+	fi->width = 100;
+	fi->height = 100;
 	fi->depth = 50;
-	fi->spacing = 1.5e9 / max3(fi->width, fi->height, fi->depth);
+	fi->spacing = 1.5e10f / max3(fi->width, fi->height, fi->depth);
+	sim->time_delta = 0.1;
 
 	fi->xstr = 3;
 	fi->ystr = fi->xstr * fi->width;
@@ -39,25 +62,28 @@ int phy_sim_create(phy_sim *sim) {
 	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));
+		sim->fields[i].E = malloc(fi->size * sizeof(float));
+		sim->fields[i].H = malloc(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));
+	sim->field_eps = malloc(fi->size1 * sizeof(float));
+	sim->field_mu = malloc(fi->size1 * sizeof(float));
+	if (!sim->field_eps || !sim->field_mu) {
+		con_printf("phy_sim_create: out of memory\n");
+		goto error;
+	}
 
-	phy_sim_compute_const_fields(sim);
+	phy_sim_reset(sim);
 
 	itc_chan_create(&sim->ctl);
 	sim->running = false;
 
 	sim->rotate_lock = SDL_CreateMutex();
 	assert(sim->rotate_lock);
-
 	return 0;
 
 error:
@@ -85,104 +111,89 @@ void phy_sim_compute_const_fields(phy_sim *sim)
 	}
 }
 
-void phy_sim_step_1(phy_field_info *fi, float dt, float *mul_const,
-                    const float *L0, float *L2, const float *R)
+void phy_sim_step_E(phy_field_info *fi, float *mul_const, float dt,
+                    float *En, float *Ep, float *H)
 {
 	size_t x, y, z;
+	float dx = fi->spacing;
 
-#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;
+		size_t offs, offs_nx, offs_ny, offs_nz;
+		size_t offs1;
+		float dExdt, dEydt, dEzdt;
 
-		dRxdz = AT(R, x, y, z + 1)[0] - AT(R, x, y, z - 1)[0];
-		dRxdz /= 2 * fi->spacing;
+		offs = x * fi->xstr + y * fi->ystr + z * fi->zstr;
+		offs_nx = offs - fi->xstr;
+		offs_ny = offs - fi->ystr;
+		offs_nz = offs - fi->zstr;
 
-		dRzdx = AT(R, x + 1, y, z)[2] - AT(R, x - 1, y, z)[2];
-		dRzdx /= 2 * fi->spacing;
+		offs1 = x * fi->xstr1 + y * fi->ystr1 + z * fi->zstr1;
 
-		dLydt = AT1(mul_const, x, y, z) * (dRxdz - dRzdx);
+		dExdt = (H[offs+2] - H[offs_ny+2] - H[offs+1] + H[offs_nz+1]);
+		dEydt = (H[offs+0] - H[offs_nz+0] - H[offs+2] + H[offs_nx+2]);
+		dEzdt = (H[offs+1] - H[offs_nx+1] - H[offs+0] + H[offs_ny+0]);
 
-		AT(L2, x, y, z)[1] = AT(L0, x, y, z)[1] + dLydt * dt;
+		En[offs+0] = Ep[offs+0] + dt * mul_const[offs1] * dExdt / dx;
+		En[offs+1] = Ep[offs+1] + dt * mul_const[offs1] * dEydt / dx;
+		En[offs+2] = Ep[offs+2] + dt * mul_const[offs1] * dEzdt / dx;
 	}
+}
 
-	// Z
-	for (z = 0; z < fi->depth; z++)
+void phy_sim_step_H(phy_field_info *fi, float *mul_const, float dt,
+                    float *Hn, float *Hp, float *E)
+{
+	size_t x, y, z;
+	float dx = fi->spacing;
+
+	for (z = 1; z < fi->depth - 1; z++)
 	for (y = 1; y < fi->height - 1; y++)
 	for (x = 1; x < fi->width - 1; x++) {
-		float dRydx, dRxdy, dLzdt;
+		size_t offs, offs_px, offs_py, offs_pz;
+		size_t offs1;
+		float dHxdt, dHydt, dHzdt;
 
-		dRydx = AT(R, x + 1, y, z)[1] - AT(R, x - 1, y, z)[1];
-		dRydx /= 2 * fi->spacing;
+		offs = x * fi->xstr + y * fi->ystr + z * fi->zstr;
+		offs_px = offs + fi->xstr;
+		offs_py = offs + fi->ystr;
+		offs_pz = offs + fi->zstr;
 
-		dRxdy = AT(R, x, y + 1, z)[0] - AT(R, x, y - 1, z)[0];
-		dRxdy /= 2 * fi->spacing;
+		offs1 = x * fi->xstr1 + y * fi->ystr1 + z * fi->zstr1;
 
-		dLzdt = AT1(mul_const, x, y, z) * (dRydx - dRxdy);
+		dHxdt = (E[offs_py+2] - E[offs+2] - E[offs_pz+1] + E[offs+1]);
+		dHydt = (E[offs_pz+0] - E[offs+0] - E[offs_px+2] + E[offs+2]);
+		dHzdt = (E[offs_px+1] - E[offs+1] - E[offs_py+0] + E[offs+0]);
 
-		AT(L2, x, y, z)[2] = AT(L0, x, y, z)[2] + dLzdt * dt;
+		Hn[offs+0] = Hp[offs+0] + dt * mul_const[offs1] * dHxdt / dx;
+		Hn[offs+1] = Hp[offs+1] + dt * mul_const[offs1] * dHydt / dx;
+		Hn[offs+2] = Hp[offs+2] + dt * mul_const[offs1] * dHzdt / dx;
 	}
-#undef AT
-#undef AT1
 }
 
 void phy_sim_add_sources(phy_sim *sim)
 {
 	phy_field_info *fi = &sim->field_info;
+	float *E;
 	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;
+	for (z = 1; z < fi->depth - 1; z++)
+	for (y = 15; y < fi->height - 15; y++)
+	for (x = 1; x < fi->width - 1; x++) {
+		float s;
 
 		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.5 - sqrt(fx * fx + fy * fy + fz * fz)) / 0.5;
-		if (chuj < 0)
-			chuj = 0;
-		else
-			chuj = 1;
-
-		chuj *= sin(sim->time * 3) * 1;
-		//chuj *= exp(-10.0f * pow(sim->time - 0.5, 2)) * 2;
-		E[0] += chuj;
+		s = exp(-600.0f * pow((float)x / fi->width - 0.15f, 2));
+		s *= sin(sim->time);
+		//s *= exp(-2.0f * pow(sim->time - 2.0f, 2)) * 1000;
+		E[0] += s;
 	}
 }
 
-void phy_sim_step(phy_sim *sim, float dt) {
+void phy_sim_step(phy_sim *sim) {
 	phy_field_info *fi = &sim->field_info;
 	phy_field_em *fields = sim->fields, tmp;
-	size_t i;
-
-	phy_sim_add_sources(sim);
 
 	// note: only pointers are rotated, not the entire fields
 	SDL_mutexP(sim->rotate_lock);
@@ -192,12 +203,13 @@ void phy_sim_step(phy_sim *sim, float dt) {
 	memcpy(fields + 2, &tmp, sizeof(phy_field_em));
 	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);
+	phy_sim_step_H(fi, sim->field_eps, sim->time_delta,
+	               fields[2].H, fields[0].H, fields[1].E);
+	phy_sim_step_E(fi, sim->field_mu, sim->time_delta,
+	               fields[2].E, fields[0].E, fields[1].H);
+	phy_sim_add_sources(sim);
 
-	sim->time += dt;
+	sim->time += sim->time_delta;
 }
 
 void phy_run_command(phy_sim *sim, int number, void *data)
@@ -212,7 +224,12 @@ void phy_run_command(phy_sim *sim, int number, void *data)
 		break;
 
 	case PHY_CMD_STEP:
-		phy_sim_step(sim, 0.1);
+		sim->running = false;
+		phy_sim_step(sim);
+		break;
+
+	case PHY_CMD_RESET:
+		phy_sim_reset(sim);
 		break;
 	}
 }
@@ -223,15 +240,19 @@ int phy_thread(phy_sim *sim)
 	void *cmd_data;
 
 	while (1) {
+cont_main_loop:
 		if (sim->running) {
 			while (!itc_chan_pop(&sim->ctl, &cmd_num, &cmd_data)) {
 				if (cmd_num == PHY_CMD_QUIT)
 					goto out;
 
 				phy_run_command(sim, cmd_num, cmd_data);
+
+				if (!sim->running)
+					goto cont_main_loop;
 			}
 
-			phy_sim_step(sim, 0.1);
+			phy_sim_step(sim);
 		} else {
 			// do it only once per iteration
 			if (!itc_chan_pop_block(&sim->ctl, &cmd_num,
diff --git a/src/physics.h b/src/physics.h
index 533d8a4..92ec6fd 100644
--- a/src/physics.h
+++ b/src/physics.h
@@ -21,14 +21,15 @@ enum {
 	PHY_CMD_QUIT,
 	PHY_CMD_PAUSE,
 	PHY_CMD_RESUME,
-	PHY_CMD_STEP
+	PHY_CMD_STEP,
+	PHY_CMD_RESET
 };
 
 typedef struct {
 	phy_field_info field_info;
 	phy_field_em fields[3];
 	float *field_eps, *field_mu; //permittivity and permeability
-	float time;
+	float time, time_delta;
 
 	// UI stuff
 	bool running;
@@ -39,7 +40,7 @@ typedef struct {
 void phy_sim_destroy(phy_sim *sim);
 int phy_sim_create(phy_sim *sim);
 void phy_sim_compute_const_fields(phy_sim *sim);
-void phy_sim_step(phy_sim *sim, float dt);
+void phy_sim_step(phy_sim *sim);
 
 int phy_thread(phy_sim *sim);
 
diff --git a/src/ui.c b/src/ui.c
index 4e2f791..72a7c7b 100644
--- a/src/ui.c
+++ b/src/ui.c
@@ -24,7 +24,7 @@ struct {
 typedef struct {
 	phy_sim *sim;
 
-	bool dragging;
+	bool dragging, selecting;
 	r_xsection_type xsection_type;
 	float xsection_frac;
 	bool select_valid;
@@ -143,27 +143,7 @@ void ui_infof(ui_window *uiw, const char *fmt, ...)
 
 static void ui_simview_set_select(ui_simview *sv, vec2_t sel_2d)
 {
-	if (sel_2d[0] < 0.0f || sel_2d[0] > 1.0f ||
-	    sel_2d[1] < 0.0f || sel_2d[1] > 1.0f) {
-		sv->select_valid = false;
-		return;
-	}
-
-	switch (sv->xsection_type) {
-	case XSECTION_XY:
-		v3_set(sv->select, sel_2d[0], sel_2d[1], sv->xsection_frac);
-		break;
-
-	case XSECTION_XZ:
-		v3_set(sv->select, sel_2d[0], sv->xsection_frac, sel_2d[1]);
-		break;
 
-	case XSECTION_YZ:
-		v3_set(sv->select, sv->xsection_frac, sel_2d[0], sel_2d[1]);
-		break;
-	}
-
-	sv->select_valid = true;
 }
 
 void ui_event_window(SDL_Event *event, ui_window *uiw)
@@ -192,22 +172,21 @@ void ui_event_window(SDL_Event *event, ui_window *uiw)
 		break;
 
 	case SDL_MOUSEBUTTONDOWN:
-		if (event->button.button == SDL_BUTTON_LEFT &&
-		    event->button.y > sv->margin_top &&
-		    event->button.y < uiw->h - sv->margin_bottom)
-			sv->dragging = true;
-
-		if (event->button.button == SDL_BUTTON_RIGHT) {
-			vec2_t sel_2d;
+		if (event->button.y > sv->margin_top &&
+		    event->button.y < uiw->h - sv->margin_bottom) {
+			if (event->button.button == SDL_BUTTON_LEFT)
+				sv->dragging = true;
 
-			v2_div_mst2(sel_2d, uiw->mouse, sv->tf_x2s);
-			ui_simview_set_select(sv, sel_2d);
+			if (event->button.button == SDL_BUTTON_RIGHT)
+				sv->selecting = true;
 		}
 		break;
 
 	case SDL_MOUSEBUTTONUP:
 		if (event->button.button == SDL_BUTTON_LEFT)
 			sv->dragging = false;
+		else if (event->button.button == SDL_BUTTON_RIGHT)
+			sv->selecting = false;
 		break;
 
 	case SDL_MOUSEMOTION:
@@ -256,6 +235,14 @@ void ui_event_window(SDL_Event *event, ui_window *uiw)
 		case SDLK_p:
 			itc_chan_push(&sv->sim->ctl, PHY_CMD_RESUME, NULL);
 			break;
+
+		case SDLK_l:
+			itc_chan_push(&sv->sim->ctl, PHY_CMD_STEP, NULL);
+			break;
+
+		case SDLK_r:
+			itc_chan_push(&sv->sim->ctl, PHY_CMD_RESET, NULL);
+			break;
 		}
 	}
 }
@@ -342,6 +329,38 @@ void ui_draw_simview_selection(ui_window *uiw)
 	phy_sim *sim = sv->sim;
 	vec2_t select_s;
 
+	if (sv->selecting)
+	{
+		vec2_t select_2d;
+
+		v2_div_mst2(select_2d, uiw->mouse, sv->tf_x2s);
+
+		if (select_2d[0] < 0.0f || select_2d[0] > 1.0f ||
+		    select_2d[1] < 0.0f || select_2d[1] > 1.0f) {
+			sv->select_valid = false;
+			return;
+		}
+
+		switch (sv->xsection_type) {
+		case XSECTION_XY:
+			v3_set(sv->select,
+			      select_2d[0], select_2d[1], sv->xsection_frac);
+			break;
+
+		case XSECTION_XZ:
+			v3_set(sv->select,
+			       select_2d[0], sv->xsection_frac, select_2d[1]);
+			break;
+
+		case XSECTION_YZ:
+			v3_set(sv->select,
+			       sv->xsection_frac, select_2d[0], select_2d[1]);
+			break;
+		}
+
+		sv->select_valid = true;
+	}
+
 	if (!sv->select_valid)
 		return;
 
@@ -377,6 +396,8 @@ void ui_draw_simview_bars(ui_window *uiw, float dt)
 
 	ui_simview *sv = &uiw->simview;
 	phy_sim *sim = sv->sim;
+	phy_field_em *field_em = sim->fields + 1;
+	phy_field_info *fi = &sim->field_info;
 
 	// upper
 
@@ -406,16 +427,30 @@ void ui_draw_simview_bars(ui_window *uiw, float dt)
 	if (sv->select_valid) {
 		size_t x, y, z, offs;
 
-		sv->margin_bottom += 2 * theme.font_size;
+		sv->margin_bottom += 3 * theme.font_size;
 
-		r_draw_rect(uiw->rw, 0, uiw->h - 3 * theme.font_size,
-			    uiw->w, 2 * theme.font_size, theme.color_main);
+		r_draw_rect(uiw->rw, 0, uiw->h - 4 * theme.font_size,
+			    uiw->w, 3 * theme.font_size, theme.color_main);
 
-		r_draw_text(uiw->rw, 0, uiw->h - theme.font_size,
-			    theme.font_size,
-			    va("x = [%f %f %f]", sv->select[0],
-			    sv->select[1], sv->select[2]),
+		x = sv->select[0] * fi->width;
+		y = sv->select[1] * fi->height;
+		z = sv->select[2] * fi->depth;
+		offs = z * fi->zstr + y * fi->ystr + x * fi->xstr;
+
+		r_draw_text(uiw->rw, 0, uiw->h - 4 * theme.font_size,
+			    theme.font_size, va("x = [%zu %zu %zu]", x, y, z),
 			    theme.color_text, 0);
+
+
+		r_draw_text(uiw->rw, 0, uiw->h - 3 * theme.font_size,
+			    theme.font_size, va("E(x) = [%+E %+E %+E]",
+			    field_em->E[offs], field_em->E[offs + 1],
+			    field_em->E[offs + 2]), theme.color_text, 0);
+
+		r_draw_text(uiw->rw, 0, uiw->h - 2 * theme.font_size,
+			    theme.font_size, va("H(x) = [%+E %+E %+E]",
+			    field_em->H[offs], field_em->H[offs + 1],
+			    field_em->H[offs + 2]), theme.color_text, 0);
 	}
 }