From d44eacf9e805904867dd28524ff2b76540cb6ae1 Mon Sep 17 00:00:00 2001 From: Paweł Redman Date: Wed, 6 Apr 2016 20:17:11 +0200 Subject: PBC on sides, fixes plane wave source. --- src/physics.c | 84 ++++++++++++++++++++++++++++++++++++++--------------------- src/physics.h | 3 ++- 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; -- cgit