Auxiliary fields (WIP).

5 files changed, 241 insertions, 81 deletions

diff --git a/src/common.h b/src/common.h
index 8933749..743dea4 100644
--- a/src/common.h
+++ b/src/common.h
@@ -85,6 +85,34 @@ static inline void v2_mul(vec_t *R, vec_t *A, vec_t b)
 	R[1] = A[1] * b;
 }
 
+static inline void v3_add(vec_t *R, vec_t *A, vec_t *B)
+{
+	R[0] = A[0] + B[0];
+	R[1] = A[1] + B[1];
+	R[2] = A[2] + B[2];
+}
+
+static inline void v3_acc(vec_t *R, vec_t *A)
+{
+	R[0] += A[0];
+	R[1] += A[1];
+	R[2] += A[2];
+}
+
+static inline void v3_sub(vec_t *R, vec_t *A, vec_t *B)
+{
+	R[0] = A[0] - B[0];
+	R[1] = A[1] - B[1];
+	R[2] = A[2] - B[2];
+}
+
+static inline void v3_mul(vec_t *R, vec_t *A, vec_t b)
+{
+	R[0] = A[0] * b;
+	R[1] = A[1] * b;
+	R[2] = A[2] * b;
+}
+
 static inline void mst2_set_identity(mst2_t R)
 {
 	R[0] = 0;
diff --git a/src/physics.c b/src/physics.c
index 7926478..a56ad97 100644
--- a/src/physics.c
+++ b/src/physics.c
@@ -16,6 +16,27 @@ void phy_sim_destroy(phy_sim *sim)
 	}
 }
 
+void phy_sim_reset_aux(phy_sim *sim)
+{
+	size_t x, y, z;
+	phy_field_info *fi = &sim->field_info;
+
+	for (z = 0; z < fi->depth; z++)
+	for (y = 0; y < fi->height; y++)
+	for (x = 0; x < fi->width; x++) {
+		phy_field_aux_point *point;
+
+		point = sim->aux + z * fi->zstr1 + y * fi->ystr1 +
+			x * fi->xstr1;
+
+		point->eps = 1.0f / EPS_ZERO;
+		point->mu = -1.0f / MU_ZERO;
+
+		v3_set(point->int_curl_E, 0, 0, 0);
+		v3_set(point->int_curl_H, 0, 0, 0);
+	}
+}
+
 void phy_sim_reset(phy_sim *sim)
 {
 	size_t i;
@@ -28,9 +49,7 @@ void phy_sim_reset(phy_sim *sim)
 		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);
+	phy_sim_reset_aux(sim);
 
 	sim->time = 0;
 	sim->running = false;
@@ -49,8 +68,11 @@ int phy_sim_create(phy_sim *sim)
 	fi->width = 100;
 	fi->height = 100;
 	fi->depth = 100;
-	fi->spacing = 1.5e10f / max3(fi->width, fi->height, fi->depth);
-	sim->time_delta = 0.2;
+	fi->spacing = 10e+9f / max3(fi->width, fi->height, fi->depth);
+	sim->time_delta = 0.1;
+
+	printf("phy_sim_create: Courant number is %f\n",
+	       3 * LIGHT_SPEED * sim->time_delta / fi->spacing);
 
 	fi->xstr = 3;
 	fi->ystr = fi->xstr * fi->width;
@@ -71,9 +93,8 @@ int phy_sim_create(phy_sim *sim)
 		}
 	}
 
-	sim->field_eps = malloc(fi->size1 * sizeof(float));
-	sim->field_mu = malloc(fi->size1 * sizeof(float));
-	if (!sim->field_eps || !sim->field_mu) {
+	sim->aux = malloc(fi->size1 * sizeof(phy_field_aux_point));
+	if (!sim->aux) {
 		con_printf("phy_sim_create: out of memory\n");
 		goto error;
 	}
@@ -92,83 +113,119 @@ error:
 	return 1;
 }
 
-void phy_sim_compute_const_fields(phy_sim *sim)
+void phy_sim_update_aux_H(phy_sim *sim)
 {
-	size_t x, y, z;
 	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 *eps, *mu;
+	for (z = 1; z < fi->depth - 1; z++)
+	for (y = 1; y < fi->height - 1; y++)
+	for (x = 1; x < fi->width - 1; x++) {
+		size_t offs, offs_px, offs_py, offs_pz;
+		phy_field_aux_point *aux;
+		float *E = sim->fields[1].E, *H = sim->fields[1].H;
+
+		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;
 
-		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;
+		aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 +
+		         z * fi->zstr1;
 
-		*eps = 1.129409066837282e+11;
-		*mu = -7.957747154594767e+5;
+		aux->curl_E[0] = (E[offs_py+2] - E[offs+2] -
+		                  E[offs_pz+1] + E[offs+1]) / fi->spacing;
+		aux->curl_E[1] = (E[offs_pz+0] - E[offs+0] -
+		                  E[offs_px+2] + E[offs+2]) / fi->spacing;
+		aux->curl_E[2] = (E[offs_px+1] - E[offs+1] -
+		                  E[offs_py+0] + E[offs+0]) / fi->spacing;
+
+		v3_acc(aux->int_H, H + offs);
+		v3_acc(aux->int_curl_E, aux->curl_E);
 	}
 }
 
-void phy_sim_step_E(phy_field_info *fi, float *mul_const, float dt,
-                    float *En, float *Ep, float *H)
+void phy_sim_update_aux_E(phy_sim *sim)
 {
+	phy_field_info *fi = &sim->field_info;
 	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++) {
 		size_t offs, offs_nx, offs_ny, offs_nz;
-		size_t offs1;
-		float dExdt, dEydt, dEzdt;
+		phy_field_aux_point *aux;
+		float *E = sim->fields[1].E, *H = sim->fields[1].H;
 
 		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;
 
-		offs1 = x * fi->xstr1 + y * fi->ystr1 + z * fi->zstr1;
+		aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 +
+		         z * fi->zstr1;
 
-		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]);
+		aux->curl_H[0] = (H[offs+2] - H[offs_ny+2] -
+		                  H[offs+1] + H[offs_nz+1]) / fi->spacing;
+		aux->curl_H[1] = (H[offs+0] - H[offs_nz+0] -
+		                  H[offs+2] + H[offs_nx+2]) / fi->spacing;
+		aux->curl_H[2] = (H[offs+1] - H[offs_nx+1] -
+		                  H[offs+0] + H[offs_ny+0]) / fi->spacing;
 
-		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;
+		v3_acc(aux->int_E, E + offs);
+		v3_acc(aux->int_curl_H, aux->curl_H);
 	}
 }
 
-void phy_sim_step_H(phy_field_info *fi, float *mul_const, float dt,
-                    float *Hn, float *Hp, float *E)
+void phy_sim_update_em_H(phy_sim *sim)
 {
+	phy_field_info *fi = &sim->field_info;
 	size_t x, y, z;
-	float dx = fi->spacing;
+	float *Hn, *Hp, dt;
+
+	Hn = sim->fields[2].H;
+	Hp = sim->fields[0].H;
+	dt = sim->time_delta;
 
 	for (z = 1; z < fi->depth - 1; z++)
 	for (y = 1; y < fi->height - 1; y++)
 	for (x = 1; x < fi->width - 1; x++) {
-		size_t offs, offs_px, offs_py, offs_pz;
-		size_t offs1;
-		float dHxdt, dHydt, dHzdt;
+		size_t offs;
+		phy_field_aux_point *aux;
 
 		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;
+		aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 +
+		         z * fi->zstr1;
+
+		Hn[offs+0] = Hp[offs+0] + dt * aux->eps * aux->curl_E[0];
+		Hn[offs+1] = Hp[offs+1] + dt * aux->eps * aux->curl_E[1];
+		Hn[offs+2] = Hp[offs+2] + dt * aux->eps * aux->curl_E[2];
+	}
+}
+
+void phy_sim_update_em_E(phy_sim *sim)
+{
+	phy_field_info *fi = &sim->field_info;
+	size_t x, y, z;
+	float *En, *Ep, dt;
+
+	En = sim->fields[2].E;
+	Ep = sim->fields[0].E;
+	dt = sim->time_delta;
 
-		offs1 = x * fi->xstr1 + y * fi->ystr1 + z * fi->zstr1;
+	for (z = 1; z < fi->depth - 1; z++)
+	for (y = 1; y < fi->height - 1; y++)
+	for (x = 1; x < fi->width - 1; x++) {
+		size_t offs;
+		phy_field_aux_point *aux;
 
-		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]);
+		offs = x * fi->xstr + y * fi->ystr + z * fi->zstr;
+		aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 +
+		         z * fi->zstr1;
 
-		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;
+		En[offs+0] = Ep[offs+0] + dt * aux->mu * aux->curl_H[0];
+		En[offs+1] = Ep[offs+1] + dt * aux->mu * aux->curl_H[1];
+		En[offs+2] = Ep[offs+2] + dt * aux->mu * aux->curl_H[2];
 	}
 }
 
@@ -178,13 +235,16 @@ void phy_sim_add_sources_E(phy_sim *sim)
 	float *E;
 	size_t x, y, z;
 
-	x = fi->width / 2;
-	y = fi->height / 2;
-	z = fi->depth / 2;
+	x = 0;
+	for (y = 1; y < fi->height - 1; y++)
+	for (z = 1; z < fi->depth - 1; z++) {
+		E = sim->fields[2].E + z * fi->zstr + y * fi->ystr +
+		     x * fi->xstr;
 
-	E = sim->fields[2].E + z * fi->zstr + y * fi->ystr + x * fi->xstr;
-	E[0] += sin(sim->time * 0.5) *
-	        exp(-0.1f * pow(sim->time - 3.0f, 2)) * 50;
+		E[0] += sin(sim->time) * 100;
+		E[1] += 0;
+		E[2] += 0;
+	}
 }
 
 void phy_sim_add_sources_H(phy_sim *sim)
@@ -193,27 +253,31 @@ void phy_sim_add_sources_H(phy_sim *sim)
 	float *H;
 	size_t x, y, z;
 
-	x = fi->width / 2;
-	y = fi->height / 2;
-	z = fi->depth / 2;
+	x = 0;
+	for (y = 1; y < fi->height - 1; y++)
+	for (z = 1; z < fi->depth - 1; z++) {
+		H = sim->fields[2].H + z * fi->zstr + y * fi->ystr +
+		     x * fi->xstr;
 
-	H = sim->fields[2].H + z * fi->zstr + y * fi->ystr + x * fi->xstr;
-	H[1] -= exp(-0.1f * pow(sim->time - 5.0f - sim->time_delta / 2, 2)) * 5;
+		H[0] += 0;
+		H[1] += 0;
+		H[2] += sin(sim->time + sim->time_delta / 2)  / IMPEDANCE_ZERO * 100;
+	}
 }
 
+
 void phy_sim_step(phy_sim *sim) {
-	phy_field_info *fi = &sim->field_info;
 	phy_field_em *fields = sim->fields, tmp;
 	int64_t start_time;
 
 	start_time = get_time();
 
-	phy_sim_step_E(fi, sim->field_mu, sim->time_delta,
-	               fields[2].E, fields[0].E, fields[1].H);
-	phy_sim_add_sources_E(sim);
-	phy_sim_step_H(fi, sim->field_eps, sim->time_delta,
-	               fields[2].H, fields[0].H, fields[1].E);
+	phy_sim_update_aux_H(sim);
+	phy_sim_update_em_H(sim);
 	phy_sim_add_sources_H(sim);
+	phy_sim_update_aux_E(sim);
+	phy_sim_update_em_E(sim);
+	phy_sim_add_sources_E(sim);
 
 	SDL_mutexP(sim->rotate_lock);
 
@@ -230,6 +294,31 @@ void phy_sim_step(phy_sim *sim) {
 	SDL_mutexV(sim->rotate_lock);
 }
 
+void phy_sim_debug(phy_sim *sim, size_t *coords)
+{
+	phy_field_info *fi = &sim->field_info;
+	size_t offs1;
+	phy_field_aux_point *aux;
+
+	offs1 = coords[2] * fi->zstr1 + coords[1] * fi->ystr1 +
+	        coords[0] * fi->xstr1;
+
+	aux = sim->aux + offs1;
+
+	printf("===[ phy_sim_debug ]===\n");
+	printf("x = [%zu %zu %zu], offs1=%zu\n", coords[0], coords[1],
+	       coords[2], offs1);
+	printf("CE = [%E %E %E]\n", aux->curl_E[0], aux->curl_E[1],
+	       aux->curl_E[2]);
+	printf("CH = [%E %E Ef]\n", aux->curl_H[0], aux->curl_H[1],
+	       aux->curl_H[2]);
+	printf("ICE = [%E %E %E]\n", aux->int_curl_E[0], aux->int_curl_E[1],
+	       aux->int_curl_E[2]);
+	printf("ICH = [%E %E %E]\n", aux->int_curl_H[0], aux->int_curl_H[1],
+	       aux->int_curl_H[2]);
+	printf("=======================\n");
+}
+
 void phy_run_command(phy_sim *sim, int number, void *data)
 {
 	switch (number) {
@@ -249,12 +338,17 @@ void phy_run_command(phy_sim *sim, int number, void *data)
 	case PHY_CMD_RESET:
 		phy_sim_reset(sim);
 		break;
+
+	case PHY_CMD_DEBUG:
+		phy_sim_debug(sim, data);
+		free(data);
+		break;
 	}
 }
 
 int phy_thread(phy_sim *sim)
 {
-	int rv, cmd_num;
+	int cmd_num;
 	void *cmd_data;
 
 	while (1) {
@@ -286,5 +380,6 @@ cont_main_loop:
 out:
 
 	printf("phy_thread: quitting\n");
+	return 0;
 }
 
diff --git a/src/physics.h b/src/physics.h
index 38a6100..2557d3c 100644
--- a/src/physics.h
+++ b/src/physics.h
@@ -4,11 +4,15 @@
 #include "common.h"
 #include "itc.h"
 
+#define LIGHT_SPEED (299792458.0f)
+#define MU_ZERO (4.0f * M_PI * 10e-7f)
+#define EPS_ZERO (1.0f / (MU_ZERO * LIGHT_SPEED * LIGHT_SPEED))
+#define IMPEDANCE_ZERO (sqrt(MU_ZERO/EPS_ZERO))
+
 typedef struct {
 	size_t width, height, depth;
 	size_t zstr, ystr, xstr, size; // strides in no. of points (not bytes)
-	size_t zstr1, ystr1, xstr1, size1; // same as above but for associated
-	                                   // scalar fields
+	size_t zstr1, ystr1, xstr1, size1; // same as above but for aux
 	float spacing; // in meters
 } phy_field_info;
 
@@ -17,18 +21,28 @@ typedef struct {
 	float *H;
 } phy_field_em;
 
+typedef struct {
+	float eps, mu; // permittivity and permeability
+
+	vec3_t curl_E, curl_H;
+
+	vec3_t int_E, int_H;
+	vec3_t int_curl_E, int_curl_H;
+} phy_field_aux_point;
+
 enum {
 	PHY_CMD_QUIT,
 	PHY_CMD_PAUSE,
 	PHY_CMD_RESUME,
 	PHY_CMD_STEP,
-	PHY_CMD_RESET
+	PHY_CMD_RESET,
+	PHY_CMD_DEBUG
 };
 
 typedef struct {
 	phy_field_info field_info;
 	phy_field_em fields[3];
-	float *field_eps, *field_mu; //permittivity and permeability
+	phy_field_aux_point *aux;
 	float time, time_delta;
 
 	// UI stuff
diff --git a/src/renderer.c b/src/renderer.c
index eee6e6d..9cf8db3 100644
--- a/src/renderer.c
+++ b/src/renderer.c
@@ -367,6 +367,7 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
 	size_t width, height, x, y, z;
 	uint8_t *pixels;
 	int pitch;
+	float scale;
 
 	if (xsection->frame_index &&
 	    xsection->frame_index == sim->frame_index &&
@@ -385,10 +386,14 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
 		height = fi->depth;
 	}
 
-	if (flags & XSECTION_E)
+	if (flags & XSECTION_E) {
 		field = sim->fields[1].E;
-	else
+		scale = 1;
+	}
+	else {
 		field = sim->fields[1].H;
+		scale = IMPEDANCE_ZERO;
+	}
 
 	if (frac < 0.0f)
 		frac = 0.0f;
@@ -416,9 +421,9 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
 			point = field + z * fi->zstr + y * fi->ystr +
 			        x * fi->xstr;
 
-			pixel[0] = r_float_to_u8(point[0]);
-			pixel[1] = r_float_to_u8(point[1]);
-			pixel[2] = r_float_to_u8(point[2]);
+			pixel[0] = r_float_to_u8(point[0] * scale);
+			pixel[1] = r_float_to_u8(point[1] * scale);
+			pixel[2] = r_float_to_u8(point[2] * scale);
 		}
 	} else if (flags & XSECTION_XZ) {
 		y = frac * fi->height;
@@ -434,9 +439,9 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
 			point = field + z * fi->zstr + y * fi->ystr +
 			        x * fi->xstr;
 
-			pixel[0] = r_float_to_u8(point[0]);
-			pixel[1] = r_float_to_u8(point[1]);
-			pixel[2] = r_float_to_u8(point[2]);
+			pixel[0] = r_float_to_u8(point[0] * scale);
+			pixel[1] = r_float_to_u8(point[1] * scale);
+			pixel[2] = r_float_to_u8(point[2] * scale);
 		}
 	} else {
 		x = frac * fi->width;
@@ -452,9 +457,9 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
 			point = field + z * fi->zstr + y * fi->ystr +
 			        x * fi->xstr;
 
-			pixel[0] = r_float_to_u8(point[0]);
-			pixel[1] = r_float_to_u8(point[1]);
-			pixel[2] = r_float_to_u8(point[2]);
+			pixel[0] = r_float_to_u8(point[0] * scale);
+			pixel[1] = r_float_to_u8(point[1] * scale);
+			pixel[2] = r_float_to_u8(point[2] * scale);
 		}
 	}
 
diff --git a/src/ui.c b/src/ui.c
index fac56ce..f485cfe 100644
--- a/src/ui.c
+++ b/src/ui.c
@@ -258,6 +258,24 @@ void ui_event_window(SDL_Event *event, ui_window *uiw)
 			sv->xsection_flags &= ~XSECTION_E;
 			sv->xsection_flags |= XSECTION_H;
 			break;
+
+		case SDLK_z:
+			if (sv->select_valid)
+			{
+				phy_field_info *fi = &sv->sim->field_info;
+				size_t *coords;
+
+				coords = calloc(3, sizeof(size_t));
+				assert(coords);
+
+				coords[0] = sv->select[0] * fi->width;
+				coords[1] = sv->select[1] * fi->height;
+				coords[2] = sv->select[2] * fi->depth;
+
+				itc_chan_push(&sv->sim->ctl, PHY_CMD_DEBUG,
+				              coords);
+				break;
+			}
 		}
 	}
 }