summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorPaweł Redman <pawel.redman@gmail.com>2016-04-04 13:44:11 +0200
committerPaweł Redman <pawel.redman@gmail.com>2016-04-04 13:44:11 +0200
commit5e519dca5f8fbf567f8d07f72327f62659eee14d (patch)
treed11b2afc0b9aca1c75eebe03d30d38b530e7fbff
parenta681d0a7fe91293728052cad405214dcbeb46ad5 (diff)
Fix phy_calc_uec_E. Plane wave source still broken.
-rw-r--r--src/physics.c79
-rw-r--r--src/physics.h2
-rw-r--r--src/renderer.c2
3 files changed, 57 insertions, 26 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;
}
diff --git a/src/physics.h b/src/physics.h
index 798fa8c..73ea4d7 100644
--- a/src/physics.h
+++ b/src/physics.h
@@ -36,6 +36,8 @@ typedef struct {
typedef struct {
size_t z;
float Px, Py;
+
+ float lambda, time;
} phy_source_xyplane;
enum {
diff --git a/src/renderer.c b/src/renderer.c
index 52f4208..3f01e88 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 = 1; //scale = IMPEDANCE_ZERO;
+ scale = 1;
}
if (frac < 0.0f)