From 5e519dca5f8fbf567f8d07f72327f62659eee14d Mon Sep 17 00:00:00 2001 From: Paweł Redman Date: Mon, 4 Apr 2016 13:44:11 +0200 Subject: Fix phy_calc_uec_E. Plane wave source still broken. --- src/physics.c | 79 ++++++++++++++++++++++++++++++++++++++++------------------- 1 file changed, 54 insertions(+), 25 deletions(-) (limited to 'src/physics.c') 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; } -- cgit