PBC on sides, fixes plane wave source.

2 files changed, 56 insertions, 31 deletions

diff --git a/src/physics.c b/src/physics.c
index a50ba1f..0946b06 100644
--- a/src/physics.c
+++ b/src/physics.c
@@ -186,9 +186,9 @@ int phy_sim_create(phy_sim *sim)
 	assert(sim->rotate_lock);
 	itc_chan_create(&sim->ctl);
 
-	con_var_set("width", "50");
-	con_var_set("height", "50");
-	con_var_set("depth", "50");
+	con_var_set("width", "40");
+	con_var_set("height", "40");
+	con_var_set("depth", "70");
 	con_var_set("scale", "7e+9"); 
 	con_var_set("time_delta", "0.1");
 	con_var_set("pml_x", "0");
@@ -197,6 +197,8 @@ int phy_sim_create(phy_sim *sim)
 
 	con_var_set("src_z", "10"); 
 	con_var_set("src_angle", "0");
+	con_var_set("src_f", "0.1");  
+	con_var_set("src_wtf", "0");
                               
 	con_var_set("csg_shape", "sphere");
 	con_var_set("csg_v", NULL);
@@ -231,11 +233,11 @@ static void phy_sim_configure(phy_sim *sim)
 	con_var_set("courant", va("%f", 3 * LIGHT_SPEED * sim->time_delta /
 	                          fi->spacing));
 
-	sim->source_xyplane.z = con_var_get_size("src_z");
-	sim->source_xyplane.Px = cos(DEG2RAD(con_var_get_float("src_angle")));
-	sim->source_xyplane.Py = sin(DEG2RAD(con_var_get_float("src_angle")));
-	sim->source_xyplane.enabled = true;
-
+	sim->source_z = con_var_get_size("src_z");
+	sim->source_f = con_var_get_float("src_f");
+	sim->source_px = cos(DEG2RAD(con_var_get_float("src_angle")));
+	sim->source_py = sin(DEG2RAD(con_var_get_float("src_angle")));
+	sim->source_wtf = con_var_get_float("src_wtf");
 }
 
 int phy_sim_create_fields(phy_sim *sim)
@@ -307,9 +309,9 @@ void phy_sim_destroy_fields(phy_sim *sim)
 	sim->valid = false;
 }
 
-static float gauss(phy_source_xyplane *s, float t)
+static float phy_sim_eval_source(phy_sim *sim, float t)
 {
-	return sin(t); //exp(-s->lambda * sq(t - s->time));
+	return sin(t * 2.0f * M_PI * sim->source_f);
 }
 
 void phy_sim_add_sources_H(phy_sim *sim)
@@ -318,20 +320,13 @@ void phy_sim_add_sources_H(phy_sim *sim)
 	size_t x, y, z;
 	float src_A, src_x, src_y, delay;
 
-	if (!sim->source_xyplane.enabled)
-		return;
-
-	z = sim->source_xyplane.z;
-
-	delay = fi->spacing / (2.0f * LIGHT_SPEED)
-		- sim->time_delta / 2.0f + 0.35;
+	z = sim->source_z;
 
-	/*delay = - fi->spacing / (2.0f * LIGHT_SPEED)
-	        + sim->time_delta / 2.0f;*/
+	delay = sim->source_wtf;
 
-	src_A = gauss(&sim->source_xyplane, sim->time + delay) * 1;
-	src_x = -sim->source_xyplane.Py * src_A;
-	src_y = sim->source_xyplane.Px * src_A;
+	src_A = phy_sim_eval_source(sim, sim->time + delay);
+	src_x = -sim->source_py * src_A;
+	src_y = sim->source_px * src_A;
 
 	for (y = 1; y < fi->dims[1] - 1; y++)
 	for (x = 1; x < fi->dims[0] - 1; x++) {
@@ -351,14 +346,11 @@ void phy_sim_add_sources_E(phy_sim *sim)
 	size_t x, y, z;
 	float src_x, src_y, src_A;
 
-	if (!sim->source_xyplane.enabled)
-		return;
+	z = sim->source_z - 1;
 
-	z = sim->source_xyplane.z - 1;
-
-	src_A = gauss(&sim->source_xyplane, sim->time) * 1;
-	src_x = sim->source_xyplane.Px * src_A;
-	src_y = sim->source_xyplane.Py * src_A;
+	src_A = phy_sim_eval_source(sim, sim->time);
+	src_x = sim->source_px * src_A;
+	src_y = sim->source_py * src_A;
 
 	for (y = 1; y < fi->dims[1] - 1; y++)
 	for (x = 1; x < fi->dims[0] - 1; x++) {
@@ -494,6 +486,36 @@ void phy_sim_update_em_E(phy_sim *sim)
 	}
 }
 
+void phy_sim_wrap_field(phy_sim *sim, float *field)
+{
+	phy_field_info *fi = &sim->field_info;
+	size_t x[3];
+
+#define AT(x,y,z) (field + (x) * fi->xstr + (y) * fi->ystr + \
+                   (z) * fi->zstr)
+
+	x[0] = 0;
+	for (x[2] = 0; x[2] < fi->dims[2]; x[2]++)
+	for (x[1] = 0; x[1] < fi->dims[1]; x[1]++)
+		v3_copy(AT(x[0], x[1], x[2]), AT(fi->dims[0] - 2, x[1], x[2]));
+
+	x[0] = fi->dims[0] - 1;                   
+	for (x[2] = 0; x[2] < fi->dims[2]; x[2]++)
+	for (x[1] = 0; x[1] < fi->dims[1]; x[1]++)
+		v3_copy(AT(fi->dims[0] - 1, x[1], x[2]), AT(1, x[1], x[2]));
+
+	x[1] = 0;                                 
+	for (x[2] = 0; x[2] < fi->dims[2]; x[2]++)
+	for (x[0] = 0; x[0] < fi->dims[0]; x[0]++)
+		v3_copy(AT(x[0], 0, x[2]), AT(x[0], fi->dims[1] - 2, x[2]));
+
+	x[1] = fi->dims[1] - 1;
+	for (x[2] = 0; x[2] < fi->dims[2]; x[2]++)
+	for (x[0] = 0; x[0] < fi->dims[0]; x[0]++)
+		v3_copy(AT(x[0], fi->dims[1] - 1, x[2]), AT(x[0], 1, x[2]));
+#undef AT
+}
+
 void phy_sim_step(phy_sim *sim) {
 	phy_field_em *fields = sim->fields, tmp;
 	int64_t start_time;
@@ -502,8 +524,10 @@ void phy_sim_step(phy_sim *sim) {
 
 	phy_sim_update_aux_H(sim);
 	phy_sim_update_em_H(sim);
+	phy_sim_wrap_field(sim, fields[2].H);
 	phy_sim_update_aux_E(sim);
-	phy_sim_update_em_E(sim);
+	phy_sim_update_em_E(sim);            
+	phy_sim_wrap_field(sim, fields[2].E);
 
 	SDL_mutexP(sim->rotate_lock);
 
diff --git a/src/physics.h b/src/physics.h
index 97c44d4..ac9f511 100644
--- a/src/physics.h
+++ b/src/physics.h
@@ -75,7 +75,8 @@ typedef struct {
 	size_t pml_widths[3];
 	float time, time_delta;
 
-	phy_source_xyplane source_xyplane;
+	size_t source_z;
+	float source_f, source_px, source_py, source_wtf;
 
 	phy_csg_step *csg_stack;
 	float eps_min, eps_max;