Compute fake conductivity (for PML).

4 files changed, 120 insertions, 79 deletions

diff --git a/src/physics.c b/src/physics.c
index a56ad97..bbe2499 100644
--- a/src/physics.c
+++ b/src/physics.c
@@ -1,4 +1,4 @@
-#include "common.h"
+#include "common.h"                                 
 #include "itc.h"
 #include "physics.h"
 
@@ -16,24 +16,46 @@ void phy_sim_destroy(phy_sim *sim)
 	}
 }
 
+static float phy_calc_fake_cond(float *fake_cond, size_t *x, size_t *widths,
+                                size_t *dims, float dt)
+{
+	size_t i;
+
+	for (i = 0; i < 3; i++) {
+		float f;
+
+		if (x[i] <= widths[i])
+			f = 1.0f - (float)(x[i] - 1) / widths[i];
+		else if (x[i] >= dims[i] - widths[i] - 1)
+			f = (float)(x[i] + widths[i] - dims[i] + 2) / widths[i];
+		else
+			continue;
+
+		fake_cond[i] = EPS_ZERO / (2 * dt) * pow(f, 3);
+	}
+}
+
 void phy_sim_reset_aux(phy_sim *sim)
 {
-	size_t x, y, z;
+	size_t x[3], i;
 	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;
+	for (x[2] = 1; x[2] < fi->dims[2] - 1; x[2]++)
+	for (x[1] = 1; x[1] < fi->dims[1] - 1; x[1]++)
+	for (x[0] = 1; x[0] < fi->dims[0] - 1; x[0]++) {
+		phy_field_aux_point *aux;
+
+		aux = sim->aux + x[2] * fi->zstr1 + x[1] * fi->ystr1 +
+			x[0] * fi->xstr1;
 
-		point = sim->aux + z * fi->zstr1 + y * fi->ystr1 +
-			x * fi->xstr1;
+		aux->eps = 1.0f / EPS_ZERO;
+		aux->mu = -1.0f / MU_ZERO;
 
-		point->eps = 1.0f / EPS_ZERO;
-		point->mu = -1.0f / MU_ZERO;
+		v3_set(aux->int_curl_E, 0, 0, 0);
+		v3_set(aux->int_curl_H, 0, 0, 0);
 
-		v3_set(point->int_curl_E, 0, 0, 0);
-		v3_set(point->int_curl_H, 0, 0, 0);
+		phy_calc_fake_cond(aux->fake_cond, x, sim->pml_widths,
+		                   fi->dims, sim->time_delta);
 	}
 }
 
@@ -65,24 +87,28 @@ int phy_sim_create(phy_sim *sim)
 
 	fi = &sim->field_info;
 
-	fi->width = 100;
-	fi->height = 100;
-	fi->depth = 100;
-	fi->spacing = 10e+9f / max3(fi->width, fi->height, fi->depth);
+	fi->dims[0] = 50;
+	fi->dims[1] = 50;
+	fi->dims[2] = 50;
+	fi->spacing = 10e+9f / max3(fi->dims[0], fi->dims[1], fi->dims[2]);
 	sim->time_delta = 0.1;
 
 	printf("phy_sim_create: Courant number is %f\n",
 	       3 * LIGHT_SPEED * sim->time_delta / fi->spacing);
 
+	sim->pml_widths[0] = fi->dims[0] / 5;
+	sim->pml_widths[1] = fi->dims[1] / 5;
+	sim->pml_widths[2] = fi->dims[2] / 5;
+
 	fi->xstr = 3;
-	fi->ystr = fi->xstr * fi->width;
-	fi->zstr = fi->ystr * fi->height;
-	fi->size = fi->zstr * fi->depth;
+	fi->ystr = fi->xstr * fi->dims[0];
+	fi->zstr = fi->ystr * fi->dims[1];
+	fi->size = fi->zstr * fi->dims[2];
 
 	fi->xstr1 = 1;
-	fi->ystr1 = fi->xstr1 * fi->width;
-	fi->zstr1 = fi->ystr1 * fi->height;
-	fi->size1 = fi->zstr1 * fi->depth;
+	fi->ystr1 = fi->xstr1 * fi->dims[0];
+	fi->zstr1 = fi->ystr1 * fi->dims[1];
+	fi->size1 = fi->zstr1 * fi->dims[2];
 
 	for (i = 0; i < 3; i++) { 
 		sim->fields[i].E = malloc(fi->size * sizeof(float));
@@ -118,9 +144,9 @@ void phy_sim_update_aux_H(phy_sim *sim)
 	phy_field_info *fi = &sim->field_info;
 	size_t x, y, z;
 
-	for (z = 1; z < fi->depth - 1; z++)
-	for (y = 1; y < fi->height - 1; y++)
-	for (x = 1; x < fi->width - 1; x++) {
+	for (z = 1; z < fi->dims[2] - 1; z++)
+	for (y = 1; y < fi->dims[1] - 1; y++)
+	for (x = 1; x < fi->dims[0] - 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;
@@ -150,9 +176,9 @@ void phy_sim_update_aux_E(phy_sim *sim)
 	phy_field_info *fi = &sim->field_info;
 	size_t x, y, z;
 
-	for (z = 1; z < fi->depth - 1; z++)
-	for (y = 1; y < fi->height - 1; y++)
-	for (x = 1; x < fi->width - 1; x++) {
+	for (z = 1; z < fi->dims[2] - 1; z++)
+	for (y = 1; y < fi->dims[1] - 1; y++)
+	for (x = 1; x < fi->dims[0] - 1; x++) {
 		size_t offs, offs_nx, offs_ny, offs_nz;
 		phy_field_aux_point *aux;
 		float *E = sim->fields[1].E, *H = sim->fields[1].H;
@@ -187,19 +213,19 @@ void phy_sim_update_em_H(phy_sim *sim)
 	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++) {
+	for (z = 1; z < fi->dims[2] - 1; z++)
+	for (y = 1; y < fi->dims[1] - 1; y++)
+	for (x = 1; x < fi->dims[0] - 1; x++) {
 		size_t offs;
 		phy_field_aux_point *aux;
 
 		offs = x * fi->xstr + y * fi->ystr + z * fi->zstr;
 		aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 +
-		         z * fi->zstr1;
+		      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];
+		Hn[offs+0] = Hp[offs+0] + dt * aux->mu * aux->curl_E[0];
+		Hn[offs+1] = Hp[offs+1] + dt * aux->mu * aux->curl_E[1];
+		Hn[offs+2] = Hp[offs+2] + dt * aux->mu * aux->curl_E[2];
 	}
 }
 
@@ -213,9 +239,9 @@ void phy_sim_update_em_E(phy_sim *sim)
 	Ep = sim->fields[0].E;
 	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++) {
+	for (z = 1; z < fi->dims[2] - 1; z++)
+	for (y = 1; y < fi->dims[1] - 1; y++)
+	for (x = 1; x < fi->dims[0] - 1; x++) {
 		size_t offs;
 		phy_field_aux_point *aux;
 
@@ -223,9 +249,9 @@ void phy_sim_update_em_E(phy_sim *sim)
 		aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 +
 		         z * fi->zstr1;
 
-		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];
+		En[offs+0] = Ep[offs+0] + dt * aux->eps * aux->curl_H[0];
+		En[offs+1] = Ep[offs+1] + dt * aux->eps * aux->curl_H[1];
+		En[offs+2] = Ep[offs+2] + dt * aux->eps * aux->curl_H[2];
 	}
 }
 
@@ -235,34 +261,44 @@ void phy_sim_add_sources_E(phy_sim *sim)
 	float *E;
 	size_t x, y, z;
 
-	x = 0;
-	for (y = 1; y < fi->height - 1; y++)
-	for (z = 1; z < fi->depth - 1; z++) {
+
+	x = fi->dims[0]/2;
+	y = fi->dims[1]/2;
+	z = fi->dims[2]/2;
+	E = sim->fields[2].E + z * fi->zstr + y * fi->ystr +
+		     x * fi->xstr;
+
+	E[0] += exp(-40.0f * pow(sim->time - 3, 2)) * 50;
+
+/*
+	x = 10;
+	for (y = 1; y < fi->dims[1] - 1; y++)
+	for (z = 1; z < fi->dims[2] - 1; z++) {
 		E = sim->fields[2].E + z * fi->zstr + y * fi->ystr +
 		     x * fi->xstr;
 
-		E[0] += sin(sim->time) * 100;
+		E[0] += sin(sim->time * 3) * 1;
 		E[1] += 0;
 		E[2] += 0;
-	}
+	}*/
 }
 
 void phy_sim_add_sources_H(phy_sim *sim)
 {
-	phy_field_info *fi = &sim->field_info;
+	/*phy_field_info *fi = &sim->field_info;
 	float *H;
 	size_t x, y, z;
 
 	x = 0;
-	for (y = 1; y < fi->height - 1; y++)
-	for (z = 1; z < fi->depth - 1; z++) {
+	for (y = 1; y < fi->dims[1] - 1; y++)
+	for (z = 1; z < fi->dims[2] - 1; z++) {
 		H = sim->fields[2].H + z * fi->zstr + y * fi->ystr +
 		     x * fi->xstr;
 
 		H[0] += 0;
 		H[1] += 0;
 		H[2] += sin(sim->time + sim->time_delta / 2)  / IMPEDANCE_ZERO * 100;
-	}
+	}*/
 }
 
 
@@ -310,12 +346,14 @@ void phy_sim_debug(phy_sim *sim, size_t *coords)
 	       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],
+	printf("CH = [%E %E %E]\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("FC = [%E %E %E]\n", aux->fake_cond[0], aux->fake_cond[1],
+	       aux->fake_cond[2]);
 	printf("=======================\n");
 }
 
diff --git a/src/physics.h b/src/physics.h
index 2557d3c..06c5d4f 100644
--- a/src/physics.h
+++ b/src/physics.h
@@ -10,7 +10,7 @@
 #define IMPEDANCE_ZERO (sqrt(MU_ZERO/EPS_ZERO))
 
 typedef struct {
-	size_t width, height, depth;
+	size_t dims[3];
 	size_t zstr, ystr, xstr, size; // strides in no. of points (not bytes)
 	size_t zstr1, ystr1, xstr1, size1; // same as above but for aux
 	float spacing; // in meters
@@ -24,6 +24,8 @@ typedef struct {
 typedef struct {
 	float eps, mu; // permittivity and permeability
 
+	float fake_cond[3];
+
 	vec3_t curl_E, curl_H;
 
 	vec3_t int_E, int_H;
@@ -43,6 +45,7 @@ typedef struct {
 	phy_field_info field_info;
 	phy_field_em fields[3];
 	phy_field_aux_point *aux;
+	size_t pml_widths[3];
 	float time, time_delta;
 
 	// UI stuff
diff --git a/src/renderer.c b/src/renderer.c
index 9cf8db3..94f0333 100644
--- a/src/renderer.c
+++ b/src/renderer.c
@@ -376,14 +376,14 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
 		return 0; // nothing's changed
 
 	if (flags & XSECTION_XY) {
-		width = fi->width;
-		height = fi->height;
+		width = fi->dims[0];
+		height = fi->dims[1];
 	} else if (flags & XSECTION_XZ) {
-		width = fi->width;
-		height = fi->depth;
+		width = fi->dims[0];
+		height = fi->dims[2];
 	} else {
-		width = fi->height;
-		height = fi->depth;
+		width = fi->dims[1];
+		height = fi->dims[2];
 	}
 
 	if (flags & XSECTION_E) {
@@ -408,12 +408,12 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
 	SDL_LockTexture(xsection->texture, NULL, (void**)&pixels, &pitch);
 
 	if (flags & XSECTION_XY) {
-		z = frac * fi->depth;       
-		if (z >= fi->depth)
-			z = fi->depth - 1;
+		z = frac * fi->dims[2];       
+		if (z >= fi->dims[2])
+			z = fi->dims[2] - 1;
 
-		for (y = 0; y < fi->height; y++)
-		for (x = 0; x < fi->width; x++) {
+		for (y = 0; y < fi->dims[1]; y++)
+		for (x = 0; x < fi->dims[0]; x++) {
 			uint8_t *pixel;
 			float *point;
 
@@ -426,12 +426,12 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
 			pixel[2] = r_float_to_u8(point[2] * scale);
 		}
 	} else if (flags & XSECTION_XZ) {
-		y = frac * fi->height;
-		if (y >= fi->height)
-			y = fi->height - 1;
+		y = frac * fi->dims[1];
+		if (y >= fi->dims[1])
+			y = fi->dims[1] - 1;
 
-		for (z = 0; z < fi->depth; z++)
-		for (x = 0; x < fi->width; x++) {
+		for (z = 0; z < fi->dims[2]; z++)
+		for (x = 0; x < fi->dims[0]; x++) {
 			uint8_t *pixel;
 			float *point;
 
@@ -444,12 +444,12 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
 			pixel[2] = r_float_to_u8(point[2] * scale);
 		}
 	} else {
-		x = frac * fi->width;
-		if (x >= fi->width)
-			x = fi->width - 1;
+		x = frac * fi->dims[0];
+		if (x >= fi->dims[0])
+			x = fi->dims[0] - 1;
 
-		for (z = 0; z < fi->depth; z++)
-		for (y = 0; y < fi->height; y++) {
+		for (z = 0; z < fi->dims[2]; z++)
+		for (y = 0; y < fi->dims[1]; y++) {
 			uint8_t *pixel;
 			float *point;
 
diff --git a/src/ui.c b/src/ui.c
index f485cfe..4245603 100644
--- a/src/ui.c
+++ b/src/ui.c
@@ -268,9 +268,9 @@ void ui_event_window(SDL_Event *event, ui_window *uiw)
 				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;
+				coords[0] = sv->select[0] * fi->dims[0];
+				coords[1] = sv->select[1] * fi->dims[1];
+				coords[2] = sv->select[2] * fi->dims[2];
 
 				itc_chan_push(&sv->sim->ctl, PHY_CMD_DEBUG,
 				              coords);
@@ -464,9 +464,9 @@ void ui_draw_simview_bars(ui_window *uiw, float dt)
 		r_draw_rect(uiw->rw, 0, uiw->h - 4 * theme.font_size,
 			    uiw->w, 3 * theme.font_size, theme.color_main);
 
-		x = sv->select[0] * fi->width;
-		y = sv->select[1] * fi->height;
-		z = sv->select[2] * fi->depth;
+		x = sv->select[0] * fi->dims[0];
+		y = sv->select[1] * fi->dims[1];
+		z = sv->select[2] * fi->dims[2];
 		offs = z * fi->zstr + y * fi->ystr + x * fi->xstr;
 
 		r_draw_text(uiw->rw, 0, uiw->h - 4 * theme.font_size,