Fix phy_calc_uec_E. Plane wave source still broken.

1 files changed, 54 insertions, 25 deletions

diff --git a/src/physics.c b/src/physics.c
index af251c9..947eb60 100644
--- a/src/physics.c
+++ b/src/physics.c
@@ -64,12 +64,17 @@ 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[0] *= 1.0f / eps;
 	uec[1] *= 1.0f / eps;
 	uec[2] *= 1.0f / eps;
 	uec[3] *= 1.0f / eps;
 } 
 
+static float remap(float x, float a0, float b0, float a1, float b1)
+{
+	return (x - a0) * (b1 - a1) / (b0 - a0) + a1;
+}
+
 void phy_sim_reset_aux(phy_sim *sim)
 {
 	size_t x[3];
@@ -87,8 +92,23 @@ void phy_sim_reset_aux(phy_sim *sim)
 		aux = sim->aux + x[2] * fi->zstr1 + x[1] * fi->ystr1 +
 			x[0] * fi->xstr1;
 
-		aux->eps = EPS_ZERO;
-		aux->mu = MU_ZERO;
+		aux->eps = 1;
+		aux->mu = 1;
+
+		{
+			float cx, cy, cz, cr, d;
+
+			cx = 25, cy = 25, cz = 65;
+			cr = 10; 
+
+			d = sqrt(sq(x[0] - cx) + sq(x[1] - cy) + sq(x[2] - cz));
+
+			if (d < cr) {
+				aux->eps = remap(d, cr, cr * 0.8, 1, 2);
+				if (aux->eps > 2)
+					aux->eps = 2;
+			}
+		}
 
 		v3_set(aux->int_curl_E, 0, 0, 0);
 		v3_set(aux->int_curl_H, 0, 0, 0);
@@ -96,13 +116,13 @@ void phy_sim_reset_aux(phy_sim *sim)
 		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_H(aux->uec_H[0], dt, fc_H[0], fc_H[1], fc_H[2], aux->mu);
+		phy_calc_uec_H(aux->uec_H[1], dt, fc_H[1], fc_H[0], fc_H[2], aux->mu);
+		phy_calc_uec_H(aux->uec_H[2], dt, fc_H[2], fc_H[0], fc_H[1], aux->mu);
 
-		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);
+		phy_calc_uec_E(aux->uec_E[0], dt, fc_E[0], fc_E[1], fc_E[2], aux->eps);
+		phy_calc_uec_E(aux->uec_E[1], dt, fc_E[1], fc_E[0], fc_E[2], aux->eps);
+		phy_calc_uec_E(aux->uec_E[2], dt, fc_E[2], fc_E[0], fc_E[1], aux->eps);
 	}
 }
 
@@ -137,13 +157,13 @@ int phy_sim_create(phy_sim *sim)
 	fi->dims[0] = 50;
 	fi->dims[1] = 50;
 	fi->dims[2] = 100;
-	fi->spacing = 10e+9f / min3(fi->dims[0], fi->dims[1], fi->dims[2]);
-	sim->time_delta = 0.15f;
+	fi->spacing = 7e+9f / min3(fi->dims[0], fi->dims[1], fi->dims[2]);
+	sim->time_delta = 0.10f;
 
 	printf("phy_sim_create: Courant number is %f\n",
 	       3 * LIGHT_SPEED * sim->time_delta / fi->spacing);
 
-	#if 1
+	#if 0
 	sim->pml_widths[0] = 0;
 	sim->pml_widths[1] = 0;
 	sim->pml_widths[2] = fi->dims[2] / 5 * 1;
@@ -156,6 +176,8 @@ int phy_sim_create(phy_sim *sim)
 	sim->source_xyplane.z = fi->dims[2] / 2;//sim->pml_widths[2] + 5;
 	sim->source_xyplane.Px = 0;
 	sim->source_xyplane.Py = 1;
+	sim->source_xyplane.lambda = 2;
+	sim->source_xyplane.time = 5;
 
 	fi->xstr = 3;
 	fi->ystr = fi->xstr * fi->dims[0];
@@ -196,27 +218,33 @@ error:
 	return 1;
 }
 
+static float gauss(phy_source_xyplane *s, float t)
+{
+	return exp(-s->lambda * sq(t - s->time));
+}
+
 void phy_sim_add_sources_H(phy_sim *sim)
 {
 	phy_field_info *fi = &sim->field_info;
 	size_t x, y, z;
+	float src_A, src_x, src_y, delay;
 
-	z = sim->source_xyplane.z - 1;
+	z = sim->source_xyplane.z;
+
+	delay = -fi->spacing / (2.0f * LIGHT_SPEED)
+		+ sim->time_delta / 2.0f;
+
+	src_A = gauss(&sim->source_xyplane, sim->time + delay);
+	src_x = -sim->source_xyplane.Py * src_A;
+	src_y = sim->source_xyplane.Px * src_A;
 
 	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;
 	}
@@ -226,20 +254,21 @@ void phy_sim_add_sources_E(phy_sim *sim)
 {
 	phy_field_info *fi = &sim->field_info;
 	size_t x, y, z;
+	float src_x, src_y, src_A;
 
-	z = sim->source_xyplane.z;
+	z = sim->source_xyplane.z - 1;
+
+	src_A = gauss(&sim->source_xyplane, sim->time);
+	src_x = sim->source_xyplane.Px * src_A;
+	src_y = sim->source_xyplane.Py * src_A;
 
 	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;
 	}