Plane wave source (broken).

4 files changed, 112 insertions, 120 deletions

diff --git a/src/physics.c b/src/physics.c
index 952e53a..af251c9 100644
--- a/src/physics.c
+++ b/src/physics.c
@@ -17,21 +17,22 @@ void phy_sim_destroy(phy_sim *sim)
 }
 
 static void phy_calc_fc(float *fc, size_t *x, size_t *widths,
-                        size_t *dims, float dt)
+                        size_t *dims, float dt, float offset)
 {
 	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];
+		if (x[i] + offset <= widths[i])
+			f = 1.0f - (float)(x[i] + offset - 1) / widths[i];
+		else if (x[i] + offset >= dims[i] - widths[i] - 1)
+			f = (float)(x[i] + offset + widths[i] - dims[i] + 2)
+			    / widths[i];
 		else
 			continue;
 
-		fc[i] = EPS_ZERO / (2 * dt) * pow(f, 3);
+		fc[i] = EPS_ZERO / (2 * dt) * sin(0.5f * M_PI * f);
 	}
 }
 
@@ -63,11 +64,10 @@ static void phy_calc_uec_E(float *uec, float dt, float fc0, float fc1,
 	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;*/
+	uec[3] *= 1.0f / eps;
 } 
 
 void phy_sim_reset_aux(phy_sim *sim)
@@ -82,7 +82,7 @@ void phy_sim_reset_aux(phy_sim *sim)
 	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};
+		vec3_t fc_H = {0, 0, 0}, fc_E = {0, 0, 0};
 
 		aux = sim->aux + x[2] * fi->zstr1 + x[1] * fi->ystr1 +
 			x[0] * fi->xstr1;
@@ -93,26 +93,16 @@ void phy_sim_reset_aux(phy_sim *sim)
 		v3_set(aux->int_curl_E, 0, 0, 0);
 		v3_set(aux->int_curl_H, 0, 0, 0);
 
-		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);
+		phy_calc_fc(fc_H, x, sim->pml_widths, fi->dims, dt, 0.0f);
+		phy_calc_fc(fc_E, x, sim->pml_widths, fi->dims, dt, 0.5f);
+
+		phy_calc_uec_H(aux->uec_H[0], dt, fc_H[0], fc_H[1], fc_H[2], 1);
+		phy_calc_uec_H(aux->uec_H[1], dt, fc_H[1], fc_H[0], fc_H[2], 1);
+		phy_calc_uec_H(aux->uec_H[2], dt, fc_H[2], fc_H[0], fc_H[1], 1);
+
+		phy_calc_uec_E(aux->uec_E[0], dt, fc_E[0], fc_E[1], fc_E[2], 1);
+		phy_calc_uec_E(aux->uec_E[1], dt, fc_E[1], fc_E[0], fc_E[2], 1);
+		phy_calc_uec_E(aux->uec_E[2], dt, fc_E[2], fc_E[0], fc_E[1], 1);
 	}
 }
 
@@ -146,23 +136,27 @@ int phy_sim_create(phy_sim *sim)
 
 	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;
+	fi->dims[2] = 100;
+	fi->spacing = 10e+9f / min3(fi->dims[0], fi->dims[1], fi->dims[2]);
+	sim->time_delta = 0.15f;
 
 	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;
+	sim->pml_widths[0] = 0;
+	sim->pml_widths[1] = 0;
+	sim->pml_widths[2] = fi->dims[2] / 5 * 1;
 	#else
 	sim->pml_widths[0] = 0;
 	sim->pml_widths[1] = 0;
 	sim->pml_widths[2] = 0;
 	#endif
 
+	sim->source_xyplane.z = fi->dims[2] / 2;//sim->pml_widths[2] + 5;
+	sim->source_xyplane.Px = 0;
+	sim->source_xyplane.Py = 1;
+
 	fi->xstr = 3;
 	fi->ystr = fi->xstr * fi->dims[0];
 	fi->zstr = fi->ystr * fi->dims[1];
@@ -202,6 +196,55 @@ error:
 	return 1;
 }
 
+void phy_sim_add_sources_H(phy_sim *sim)
+{
+	phy_field_info *fi = &sim->field_info;
+	size_t x, y, z;
+
+	z = sim->source_xyplane.z - 1;
+
+	for (y = 1; y < fi->dims[1] - 1; y++)
+	for (x = 1; x < fi->dims[0] - 1; x++) {
+		phy_field_aux_point *aux;
+		float src_x, src_y, delay;
+
+		aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 +
+		         z * fi->zstr1;
+
+		delay = -fi->spacing / (2.0f * LIGHT_SPEED)
+		        + sim->time_delta / 2.0f;
+
+		src_x = -sim->source_xyplane.Py * cos(sim->time + delay);
+		src_y = sim->source_xyplane.Px * cos(sim->time + delay);
+
+		aux->curl_H[0] += 1.0f / fi->spacing * src_y;
+		aux->curl_H[1] += -1.0f / fi->spacing * src_x;
+	}
+}
+
+void phy_sim_add_sources_E(phy_sim *sim)
+{
+	phy_field_info *fi = &sim->field_info;
+	size_t x, y, z;
+
+	z = sim->source_xyplane.z;
+
+	for (y = 1; y < fi->dims[1] - 1; y++)
+	for (x = 1; x < fi->dims[0] - 1; x++) {
+		phy_field_aux_point *aux;
+		float src_x, src_y;
+
+		aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 +
+		      z * fi->zstr1;
+
+		src_x = sim->source_xyplane.Px * cos(sim->time);
+		src_y = sim->source_xyplane.Py * cos(sim->time);
+
+		aux->curl_E[0] += 1.0f / fi->spacing * src_y;
+		aux->curl_E[1] += -1.0f / fi->spacing * src_x;
+	}
+}
+
 void phy_sim_update_aux_H(phy_sim *sim)
 {
 	phy_field_info *fi = &sim->field_info;
@@ -220,7 +263,10 @@ void phy_sim_update_aux_H(phy_sim *sim)
 		offs_pz = offs + fi->zstr;
 
 		aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 +
-		         z * fi->zstr1;
+		      z * fi->zstr1;
+
+		v3_acc(aux->int_H, H + offs);
+		v3_acc(aux->int_curl_E, aux->curl_E);
 
 		aux->curl_E[0] = (E[offs_py+2] - E[offs+2] -
 		                  E[offs_pz+1] + E[offs+1]) / fi->spacing;
@@ -228,10 +274,9 @@ void phy_sim_update_aux_H(phy_sim *sim)
 		                  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);
 	}
+
+	phy_sim_add_sources_E(sim);
 }
 
 void phy_sim_update_aux_E(phy_sim *sim)
@@ -252,7 +297,10 @@ void phy_sim_update_aux_E(phy_sim *sim)
 		offs_nz = offs - fi->zstr;
 
 		aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 +
-		         z * fi->zstr1;
+		      z * fi->zstr1;
+
+		v3_acc(aux->int_E, E + offs);
+		v3_acc(aux->int_curl_H, aux->curl_H);
 
 		aux->curl_H[0] = (H[offs+2] - H[offs_ny+2] -
 		                  H[offs+1] + H[offs_nz+1]) / fi->spacing;
@@ -260,10 +308,9 @@ void phy_sim_update_aux_E(phy_sim *sim)
 		                  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;
-
-		v3_acc(aux->int_E, E + offs);
-		v3_acc(aux->int_curl_H, aux->curl_H);
 	}
+
+	phy_sim_add_sources_H(sim);
 }
 
 void phy_sim_update_em_H(phy_sim *sim)
@@ -290,10 +337,6 @@ void phy_sim_update_em_H(phy_sim *sim)
 			             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];
 	}
 }
 
@@ -321,60 +364,9 @@ void phy_sim_update_em_E(phy_sim *sim)
 			             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];
 	}
 }
 
-void phy_sim_add_sources_E(phy_sim *sim)
-{
-	phy_field_info *fi = &sim->field_info;
-	float *E;
-	size_t x, y, 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 * 3) * 1;
-		E[1] += 0;
-		E[2] += 0;
-	}*/
-}
-
-void phy_sim_add_sources_H(phy_sim *sim)
-{
-	/*phy_field_info *fi = &sim->field_info;
-	float *H;
-	size_t x, y, z;
-
-	x = 0;
-	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;
-	}*/
-}
-
-
 void phy_sim_step(phy_sim *sim) {
 	phy_field_em *fields = sim->fields, tmp;
 	int64_t start_time;
@@ -383,10 +375,8 @@ void phy_sim_step(phy_sim *sim) {
 
 	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);
 
@@ -406,17 +396,23 @@ void phy_sim_step(phy_sim *sim) {
 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;
+	float *E, *H, Emag, Hmag;
+
+	aux = sim->aux + coords[2] * fi->zstr1 + coords[1] * fi->ystr1 +
+	      coords[0] * fi->xstr1;
 
-	offs1 = coords[2] * fi->zstr1 + coords[1] * fi->ystr1 +
-	        coords[0] * fi->xstr1;
+	E = sim->fields[1].E + coords[2] * fi->zstr + coords[1] * fi->ystr +
+	      coords[0] * fi->xstr;
+	H = sim->fields[1].H + coords[2] * fi->zstr + coords[1] * fi->ystr +
+	      coords[0] * fi->xstr;
 
-	aux = sim->aux + offs1;
+	Emag = sqrt(sq(E[0]) + sq(E[1]) + sq(E[2]));
+	Hmag = sqrt(sq(H[0]) + sq(H[1]) + sq(H[2]));
 
 	printf("===[ phy_sim_debug ]===\n");
-	printf("x = [%zu %zu %zu], offs1=%zu, aux=%p\n", coords[0], coords[1],
-	       coords[2], offs1, aux);
+	printf("x = [%zu %zu %zu], E=%p, H=%p, aux=%p\n", coords[0], coords[1],
+	       coords[2], E, H, 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]);
@@ -434,6 +430,7 @@ void phy_sim_debug(phy_sim *sim, size_t *coords)
 	       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("Emag = %E, Hmag = %E, Z = %E\n", Emag, Hmag, Emag/Hmag);
 	printf("=======================\n");
 }
 
diff --git a/src/physics.h b/src/physics.h
index 09ee23e..798fa8c 100644
--- a/src/physics.h
+++ b/src/physics.h
@@ -33,6 +33,11 @@ typedef struct {
 	vec3_t int_curl_E, int_curl_H;
 } phy_field_aux_point;
 
+typedef struct {
+	size_t z;
+	float Px, Py;
+} phy_source_xyplane;
+
 enum {
 	PHY_CMD_QUIT,
 	PHY_CMD_PAUSE,
@@ -49,6 +54,8 @@ typedef struct {
 	size_t pml_widths[3];
 	float time, time_delta;
 
+	phy_source_xyplane source_xyplane;
+
 	// UI stuff
 	bool running;
 	itc_chan ctl;
diff --git a/src/renderer.c b/src/renderer.c
index 94f0333..52f4208 100644
--- a/src/renderer.c
+++ b/src/renderer.c
@@ -392,7 +392,7 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
 	}
 	else {
 		field = sim->fields[1].H;
-		scale = IMPEDANCE_ZERO;
+		scale = 1; //scale = IMPEDANCE_ZERO;
 	}
 
 	if (frac < 0.0f)
diff --git a/src/ui.c b/src/ui.c
index 2674756..4c02e53 100644
--- a/src/ui.c
+++ b/src/ui.c
@@ -240,18 +240,6 @@ void ui_event_window(SDL_Event *event, ui_window *uiw)
 			sv->xsection_flags |= XSECTION_YZ;
 			break;
 
-		case SDLK_w:
-			sv->xsection_frac += 0.1f;
-			if (sv->xsection_frac > 1.0f)
-				sv->xsection_frac = 1.0f;
-			break;
-
-		case SDLK_s:
-			sv->xsection_frac -= 0.1f;
-			if (sv->xsection_frac < 0.0f)
-				sv->xsection_frac = 0.0f;
-			break;
-
 		case SDLK_o:
 			ui_infof(uiw, "Symulacja wstrzymana");
 			itc_chan_push(&sv->sim->ctl, PHY_CMD_PAUSE, NULL);
@@ -575,7 +563,7 @@ void ui_draw_window(ui_window *uiw)
 
 		v2_set(sv->origin, 0, 0);
 		v2_set(sv->scale, 1, 1);
-		sv->xsection_flags = XSECTION_XY | XSECTION_E;
+		sv->xsection_flags = XSECTION_XZ | XSECTION_E;
 		sv->xsection_frac = 0.5f;
 
 		sv->select_valid = false;