PML.

5 files changed, 153 insertions, 43 deletions

diff --git a/TODO b/TODO
index d5bb629..261e76f 100644
--- a/TODO
+++ b/TODO
@@ -1 +1,3 @@
+- unstable simulation
+- wrong impedance
 itc_chan_destroy: free queue
diff --git a/src/common.h b/src/common.h
index 743dea4..9dc85e9 100644
--- a/src/common.h
+++ b/src/common.h
@@ -18,6 +18,7 @@
 #define max3(a, b, c) (max2(max2(a, b), c))
 #define min2(a, b) ((a) < (b) ? (a) : (b))
 #define min3(a, b, c) (min2(min2(a, b), c))
+#define sq(x) ((x)*(x))
 
 char *va(const char *fmt, ...);
 int64_t get_time(void);
diff --git a/src/physics.c b/src/physics.c
index bbe2499..952e53a 100644
--- a/src/physics.c
+++ b/src/physics.c
@@ -16,8 +16,8 @@ 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)
+static void phy_calc_fc(float *fc, size_t *x, size_t *widths,
+                        size_t *dims, float dt)
 {
 	size_t i;
 
@@ -31,31 +31,88 @@ static float phy_calc_fake_cond(float *fake_cond, size_t *x, size_t *widths,
 		else
 			continue;
 
-		fake_cond[i] = EPS_ZERO / (2 * dt) * pow(f, 3);
+		fc[i] = EPS_ZERO / (2 * dt) * pow(f, 3);
 	}
 }
 
+static void phy_calc_uec_H(float *uec, float dt, float fc0, float fc1,
+                           float fc2, float mu)
+{
+	float tmp;
+
+	tmp = 1.0f / dt + (fc1 + fc2) / (2 * EPS_ZERO) +
+	      fc1 * fc2 * dt / sq(2.0f * EPS_ZERO);
+	uec[0] = 1.0f / tmp * (1.0f / dt - (fc1 + fc2) / (2 * EPS_ZERO)
+	         - fc1 * fc2 * dt / sq(2.0f * EPS_ZERO));
+	uec[1] = -1.0f / tmp * LIGHT_SPEED / mu;
+	uec[2] = -1.0f / tmp * LIGHT_SPEED * dt / EPS_ZERO * fc0 / mu;
+	uec[3] = -1.0f / tmp * dt / sq(EPS_ZERO) * fc1 * fc2;
+} 
+
+static void phy_calc_uec_E(float *uec, float dt, float fc0, float fc1,
+                           float fc2, float eps)
+{
+	float tmp;
+
+	tmp = 1.0f / dt + (fc1 + fc2) / (2 * EPS_ZERO) +
+	      fc1 * fc2 * dt / sq(2.0f * EPS_ZERO);
+	uec[0] = 1.0f / tmp * (1.0f / dt - (fc1 + fc2) / (2 * EPS_ZERO)
+	         - fc1 * fc2 * dt / sq(2.0f * EPS_ZERO));
+	uec[1] = LIGHT_SPEED / tmp;
+	uec[2] = 1.0f / tmp * LIGHT_SPEED * dt * fc0 / EPS_ZERO;
+	uec[3] = 1.0f / tmp * dt / sq(EPS_ZERO) * fc1 * fc2;
+
+	// equations above are for the D field and E = 1/eps * D
+                               /*
+	uec[0] *= 1.0f / eps;
+	uec[1] *= 1.0f / eps;
+	uec[2] *= 1.0f / eps;
+	uec[3] *= 1.0f / eps;*/
+} 
+
 void phy_sim_reset_aux(phy_sim *sim)
 {
-	size_t x[3], i;
+	size_t x[3];
 	phy_field_info *fi = &sim->field_info;
+	float dt = sim->time_delta;
+
+	memset(sim->aux, 0, fi->size1 * sizeof(phy_field_aux_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;
+		vec3_t fc = {0, 0, 0};
 
 		aux = sim->aux + x[2] * fi->zstr1 + x[1] * fi->ystr1 +
 			x[0] * fi->xstr1;
 
-		aux->eps = 1.0f / EPS_ZERO;
-		aux->mu = -1.0f / MU_ZERO;
+		aux->eps = EPS_ZERO;
+		aux->mu = MU_ZERO;
 
 		v3_set(aux->int_curl_E, 0, 0, 0);
 		v3_set(aux->int_curl_H, 0, 0, 0);
 
-		phy_calc_fake_cond(aux->fake_cond, x, sim->pml_widths,
-		                   fi->dims, sim->time_delta);
+		phy_calc_fc(fc, x, sim->pml_widths, fi->dims, dt);
+                   /*
+		for (size_t i = 0; i < 3; i++) {
+			aux->uec_H[i][0] = 1;
+			aux->uec_H[i][1] = -dt / MU_ZERO;
+			aux->uec_H[i][2] = 0;
+			aux->uec_H[i][3] = 0;
+			aux->uec_E[i][0] = 1;
+			aux->uec_E[i][1] = +dt / EPS_ZERO;
+			aux->uec_E[i][2] = 0;
+			aux->uec_E[i][3] = 0;
+		}*/
+
+		phy_calc_uec_H(aux->uec_H[0], dt, fc[0], fc[1], fc[2], 1);
+		phy_calc_uec_H(aux->uec_H[1], dt, fc[1], fc[0], fc[2], 1);
+		phy_calc_uec_H(aux->uec_H[2], dt, fc[2], fc[0], fc[1], 1);
+
+		phy_calc_uec_E(aux->uec_E[0], dt, fc[0], fc[1], fc[2], 1);
+		phy_calc_uec_E(aux->uec_E[1], dt, fc[1], fc[0], fc[2], 1);
+		phy_calc_uec_E(aux->uec_E[2], dt, fc[2], fc[0], fc[1], 1);
 	}
 }
 
@@ -96,9 +153,15 @@ int phy_sim_create(phy_sim *sim)
 	printf("phy_sim_create: Courant number is %f\n",
 	       3 * LIGHT_SPEED * sim->time_delta / fi->spacing);
 
+	#if 1
 	sim->pml_widths[0] = fi->dims[0] / 5;
 	sim->pml_widths[1] = fi->dims[1] / 5;
 	sim->pml_widths[2] = fi->dims[2] / 5;
+	#else
+	sim->pml_widths[0] = 0;
+	sim->pml_widths[1] = 0;
+	sim->pml_widths[2] = 0;
+	#endif
 
 	fi->xstr = 3;
 	fi->ystr = fi->xstr * fi->dims[0];
@@ -207,25 +270,30 @@ void phy_sim_update_em_H(phy_sim *sim)
 {
 	phy_field_info *fi = &sim->field_info;
 	size_t x, y, z;
-	float *Hn, *Hp, dt;
+	float *Hn, *Hp;
 
 	Hn = sim->fields[2].H;
 	Hp = sim->fields[0].H;
-	dt = sim->time_delta;
 
 	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;
+		size_t offs, i;
 		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;
 
-		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];
+		for (i = 0; i < 3; i++)
+			Hn[offs+i] = aux->uec_H[i][0] * Hp[offs+i] +
+			             aux->uec_H[i][1] * aux->curl_E[i] +
+			             aux->uec_H[i][2] * aux->int_curl_E[i] +
+			             aux->uec_H[i][3] * aux->int_H[i];
+
+		//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];
 	}
 }
 
@@ -233,25 +301,30 @@ 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;
+	float *En, *Ep;
 
 	En = sim->fields[2].E;
 	Ep = sim->fields[0].E;
-	dt = sim->time_delta;
 
 	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;
+		size_t offs, i;
 		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;
 
-		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];
+		for (i = 0; i < 3; i++)
+			En[offs+i] = aux->uec_E[i][0] * Ep[offs+i] +
+			             aux->uec_E[i][1] * aux->curl_H[i] +
+			             aux->uec_E[i][2] * aux->int_curl_H[i] +
+			             aux->uec_E[i][3] * aux->int_E[i];
+
+		//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];
 	}
 }
 
@@ -342,8 +415,9 @@ void phy_sim_debug(phy_sim *sim, size_t *coords)
 	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("x = [%zu %zu %zu], offs1=%zu, aux=%p\n", coords[0], coords[1],
+	       coords[2], offs1, aux);
+	printf("eps = %f, mu = %f\n", aux->eps, aux->mu);
 	printf("CE = [%E %E %E]\n", aux->curl_E[0], aux->curl_E[1],
 	       aux->curl_E[2]);
 	printf("CH = [%E %E %E]\n", aux->curl_H[0], aux->curl_H[1],
@@ -352,8 +426,14 @@ void phy_sim_debug(phy_sim *sim, size_t *coords)
 	       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("UEC_E = [[%E %E %E %E] [%E %E %E %E] [%E %E %E %E]]\n",
+	       aux->uec_E[0][0], aux->uec_E[0][1], aux->uec_E[0][2], aux->uec_E[0][3],
+	       aux->uec_E[1][0], aux->uec_E[1][1], aux->uec_E[1][2], aux->uec_E[1][3],
+	       aux->uec_E[2][0], aux->uec_E[2][1], aux->uec_E[2][2], aux->uec_E[2][3]);
+	printf("UEC_H = [[%E %E %E %E] [%E %E %E %E] [%E %E %E %E]]\n",
+	       aux->uec_H[0][0], aux->uec_H[0][1], aux->uec_H[0][2], aux->uec_H[0][3],
+	       aux->uec_H[1][0], aux->uec_H[1][1], aux->uec_H[1][2], aux->uec_H[1][3],
+	       aux->uec_H[2][0], aux->uec_H[2][1], aux->uec_H[2][2], aux->uec_H[2][3]);
 	printf("=======================\n");
 }
 
diff --git a/src/physics.h b/src/physics.h
index 06c5d4f..09ee23e 100644
--- a/src/physics.h
+++ b/src/physics.h
@@ -5,7 +5,7 @@
 #include "itc.h"
 
 #define LIGHT_SPEED (299792458.0f)
-#define MU_ZERO (4.0f * M_PI * 10e-7f)
+#define MU_ZERO (4.0f * M_PI * 1e-7f)
 #define EPS_ZERO (1.0f / (MU_ZERO * LIGHT_SPEED * LIGHT_SPEED))
 #define IMPEDANCE_ZERO (sqrt(MU_ZERO/EPS_ZERO))
 
@@ -22,9 +22,10 @@ typedef struct {
 } phy_field_em;
 
 typedef struct {
-	float eps, mu; // permittivity and permeability
+	float eps, mu; // for visualisation, not needed for stepping
 
-	float fake_cond[3];
+	float uec_H[3][4];
+	float uec_E[3][4];
 
 	vec3_t curl_E, curl_H;
 
diff --git a/src/ui.c b/src/ui.c
index 4245603..3f2214a 100644
--- a/src/ui.c
+++ b/src/ui.c
@@ -24,7 +24,7 @@ struct {
 typedef struct {
 	phy_sim *sim;
 
-	bool dragging, selecting;
+	bool dragging, dragging_frac, selecting;
 	int xsection_flags;
 	float xsection_frac;
 	bool select_valid;
@@ -167,31 +167,54 @@ void ui_event_window(SDL_Event *event, ui_window *uiw)
 	case SDL_MOUSEBUTTONDOWN:
 		if (event->button.y > sv->margin_top &&
 		    event->button.y < uiw->h - sv->margin_bottom) {
-			if (event->button.button == SDL_BUTTON_LEFT)
+			switch (event->button.button) {
+			case SDL_BUTTON_LEFT:
 				sv->dragging = true;
-
-			if (event->button.button == SDL_BUTTON_RIGHT)
+				break;
+			case SDL_BUTTON_RIGHT:
+				sv->dragging_frac = true;
+				break;
+			case SDL_BUTTON_MIDDLE:
 				sv->selecting = true;
+				break;
+			}
 		}
 		break;
 
 	case SDL_MOUSEBUTTONUP:
-		if (event->button.button == SDL_BUTTON_LEFT)
+		switch (event->button.button) {
+		case SDL_BUTTON_LEFT:
 			sv->dragging = false;
-		else if (event->button.button == SDL_BUTTON_RIGHT)
+			break;
+		case SDL_BUTTON_RIGHT:
+			sv->dragging_frac = false;
+			break;
+		case SDL_BUTTON_MIDDLE:
 			sv->selecting = false;
-		break;
+			break;
+		}
 
 	case SDL_MOUSEMOTION:
 		uiw->mouse[0] = event->motion.x;
 		uiw->mouse[1] = event->motion.y;
 
-		if (sv->dragging) {
+		if (sv->dragging || sv->dragging_frac) {
 			vec2_t delta_v;
 
 			v2_set(delta_v, event->motion.xrel, event->motion.yrel);
 			v2_div_mst2_nt(delta_v, delta_v, sv->tf_v2s);
-			v2_add(sv->origin, sv->origin, delta_v);
+
+			if (sv->dragging)
+				v2_add(sv->origin, sv->origin, delta_v);
+
+			if (sv->dragging_frac) {
+				sv->xsection_frac += delta_v[1];
+
+				if (sv->xsection_frac > 1.0f)
+					sv->xsection_frac = 1.0f;
+				else if (sv->xsection_frac < 0.0f)
+					sv->xsection_frac = 0.0f;
+			}
 		}
 		break;
 
@@ -407,14 +430,17 @@ void ui_draw_simview_selection(ui_window *uiw)
 char *ui_draw_simview_top_bar_text(ui_simview *sv)
 {
 	int flags = sv->xsection_flags;
+	char *rv;
+
+	rv = va("Przekrój %s płaszczyzną %s=%f",
+	        (flags & XSECTION_E ? "E" :
+	                              "H"),
+	        (flags & XSECTION_XY ? "Z" :
+	         flags & XSECTION_XZ ? "Y" :
+	                               "X"),
+	        sv->xsection_frac);
 
-	va("Przekrój %s płaszczyzną %s=%f",
-	   (flags & XSECTION_E ? "E" :
-	                         "H"),
-	   (flags & XSECTION_XY ? "Z" :
-	    flags & XSECTION_XZ ? "Y" :
-	                          "X"),
-	   sv->xsection_frac);
+	return rv;
 }
 
 void ui_draw_simview_bars(ui_window *uiw, float dt)