From a681d0a7fe91293728052cad405214dcbeb46ad5 Mon Sep 17 00:00:00 2001 From: Paweł Redman Date: Mon, 4 Apr 2016 02:25:26 +0200 Subject: Plane wave source (broken). --- src/physics.c | 209 ++++++++++++++++++++++++++++----------------------------- src/physics.h | 7 ++ src/renderer.c | 2 +- src/ui.c | 14 +--- 4 files changed, 112 insertions(+), 120 deletions(-) (limited to 'src') diff --git a/src/physics.c b/src/physics.c index 952e53a..af251c9 100644 --- a/src/physics.c +++ b/src/physics.c @@ -17,21 +17,22 @@ void phy_sim_destroy(phy_sim *sim) } static void phy_calc_fc(float *fc, size_t *x, size_t *widths, - size_t *dims, float dt) + size_t *dims, float dt, float offset) { size_t i; for (i = 0; i < 3; i++) { float f; - if (x[i] <= widths[i]) - f = 1.0f - (float)(x[i] - 1) / widths[i]; - else if (x[i] >= dims[i] - widths[i] - 1) - f = (float)(x[i] + widths[i] - dims[i] + 2) / widths[i]; + if (x[i] + offset <= widths[i]) + f = 1.0f - (float)(x[i] + offset - 1) / widths[i]; + else if (x[i] + offset >= dims[i] - widths[i] - 1) + f = (float)(x[i] + offset + widths[i] - dims[i] + 2) + / widths[i]; else continue; - fc[i] = EPS_ZERO / (2 * dt) * pow(f, 3); + fc[i] = EPS_ZERO / (2 * dt) * sin(0.5f * M_PI * f); } } @@ -63,11 +64,10 @@ 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[1] *= 1.0f / eps; uec[2] *= 1.0f / eps; - uec[3] *= 1.0f / eps;*/ + uec[3] *= 1.0f / eps; } void phy_sim_reset_aux(phy_sim *sim) @@ -82,7 +82,7 @@ void phy_sim_reset_aux(phy_sim *sim) for (x[1] = 1; x[1] < fi->dims[1] - 1; x[1]++) for (x[0] = 1; x[0] < fi->dims[0] - 1; x[0]++) { phy_field_aux_point *aux; - vec3_t fc = {0, 0, 0}; + vec3_t fc_H = {0, 0, 0}, fc_E = {0, 0, 0}; aux = sim->aux + x[2] * fi->zstr1 + x[1] * fi->ystr1 + x[0] * fi->xstr1; @@ -93,26 +93,16 @@ void phy_sim_reset_aux(phy_sim *sim) v3_set(aux->int_curl_E, 0, 0, 0); v3_set(aux->int_curl_H, 0, 0, 0); - phy_calc_fc(fc, x, sim->pml_widths, fi->dims, dt); - /* - for (size_t i = 0; i < 3; i++) { - aux->uec_H[i][0] = 1; - aux->uec_H[i][1] = -dt / MU_ZERO; - aux->uec_H[i][2] = 0; - aux->uec_H[i][3] = 0; - aux->uec_E[i][0] = 1; - aux->uec_E[i][1] = +dt / EPS_ZERO; - aux->uec_E[i][2] = 0; - aux->uec_E[i][3] = 0; - }*/ - - phy_calc_uec_H(aux->uec_H[0], dt, fc[0], fc[1], fc[2], 1); - phy_calc_uec_H(aux->uec_H[1], dt, fc[1], fc[0], fc[2], 1); - phy_calc_uec_H(aux->uec_H[2], dt, fc[2], fc[0], fc[1], 1); - - phy_calc_uec_E(aux->uec_E[0], dt, fc[0], fc[1], fc[2], 1); - phy_calc_uec_E(aux->uec_E[1], dt, fc[1], fc[0], fc[2], 1); - phy_calc_uec_E(aux->uec_E[2], dt, fc[2], fc[0], fc[1], 1); + 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_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); } } @@ -146,23 +136,27 @@ int phy_sim_create(phy_sim *sim) fi->dims[0] = 50; fi->dims[1] = 50; - fi->dims[2] = 50; - fi->spacing = 10e+9f / max3(fi->dims[0], fi->dims[1], fi->dims[2]); - sim->time_delta = 0.1; + fi->dims[2] = 100; + fi->spacing = 10e+9f / min3(fi->dims[0], fi->dims[1], fi->dims[2]); + sim->time_delta = 0.15f; printf("phy_sim_create: Courant number is %f\n", 3 * LIGHT_SPEED * sim->time_delta / fi->spacing); #if 1 - sim->pml_widths[0] = fi->dims[0] / 5; - sim->pml_widths[1] = fi->dims[1] / 5; - sim->pml_widths[2] = fi->dims[2] / 5; + sim->pml_widths[0] = 0; + sim->pml_widths[1] = 0; + sim->pml_widths[2] = fi->dims[2] / 5 * 1; #else sim->pml_widths[0] = 0; sim->pml_widths[1] = 0; sim->pml_widths[2] = 0; #endif + sim->source_xyplane.z = fi->dims[2] / 2;//sim->pml_widths[2] + 5; + sim->source_xyplane.Px = 0; + sim->source_xyplane.Py = 1; + fi->xstr = 3; fi->ystr = fi->xstr * fi->dims[0]; fi->zstr = fi->ystr * fi->dims[1]; @@ -202,6 +196,55 @@ error: return 1; } +void phy_sim_add_sources_H(phy_sim *sim) +{ + phy_field_info *fi = &sim->field_info; + size_t x, y, z; + + z = sim->source_xyplane.z - 1; + + 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; + } +} + +void phy_sim_add_sources_E(phy_sim *sim) +{ + phy_field_info *fi = &sim->field_info; + size_t x, y, z; + + z = sim->source_xyplane.z; + + 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; + } +} + void phy_sim_update_aux_H(phy_sim *sim) { phy_field_info *fi = &sim->field_info; @@ -220,7 +263,10 @@ void phy_sim_update_aux_H(phy_sim *sim) offs_pz = offs + fi->zstr; aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 + - z * fi->zstr1; + z * fi->zstr1; + + v3_acc(aux->int_H, H + offs); + v3_acc(aux->int_curl_E, aux->curl_E); aux->curl_E[0] = (E[offs_py+2] - E[offs+2] - E[offs_pz+1] + E[offs+1]) / fi->spacing; @@ -228,10 +274,9 @@ void phy_sim_update_aux_H(phy_sim *sim) E[offs_px+2] + E[offs+2]) / fi->spacing; aux->curl_E[2] = (E[offs_px+1] - E[offs+1] - E[offs_py+0] + E[offs+0]) / fi->spacing; - - v3_acc(aux->int_H, H + offs); - v3_acc(aux->int_curl_E, aux->curl_E); } + + phy_sim_add_sources_E(sim); } void phy_sim_update_aux_E(phy_sim *sim) @@ -252,7 +297,10 @@ void phy_sim_update_aux_E(phy_sim *sim) offs_nz = offs - fi->zstr; aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 + - z * fi->zstr1; + z * fi->zstr1; + + v3_acc(aux->int_E, E + offs); + v3_acc(aux->int_curl_H, aux->curl_H); aux->curl_H[0] = (H[offs+2] - H[offs_ny+2] - H[offs+1] + H[offs_nz+1]) / fi->spacing; @@ -260,10 +308,9 @@ void phy_sim_update_aux_E(phy_sim *sim) H[offs+2] + H[offs_nx+2]) / fi->spacing; aux->curl_H[2] = (H[offs+1] - H[offs_nx+1] - H[offs+0] + H[offs_ny+0]) / fi->spacing; - - v3_acc(aux->int_E, E + offs); - v3_acc(aux->int_curl_H, aux->curl_H); } + + phy_sim_add_sources_H(sim); } void phy_sim_update_em_H(phy_sim *sim) @@ -290,10 +337,6 @@ void phy_sim_update_em_H(phy_sim *sim) aux->uec_H[i][1] * aux->curl_E[i] + aux->uec_H[i][2] * aux->int_curl_E[i] + aux->uec_H[i][3] * aux->int_H[i]; - - //Hn[offs+0] = Hp[offs+0] + dt * aux->mu * aux->curl_E[0]; - //Hn[offs+1] = Hp[offs+1] + dt * aux->mu * aux->curl_E[1]; - //Hn[offs+2] = Hp[offs+2] + dt * aux->mu * aux->curl_E[2]; } } @@ -321,60 +364,9 @@ void phy_sim_update_em_E(phy_sim *sim) aux->uec_E[i][1] * aux->curl_H[i] + aux->uec_E[i][2] * aux->int_curl_H[i] + aux->uec_E[i][3] * aux->int_E[i]; - - //En[offs+0] = Ep[offs+0] + dt * aux->eps * aux->curl_H[0]; - //En[offs+1] = Ep[offs+1] + dt * aux->eps * aux->curl_H[1]; - //En[offs+2] = Ep[offs+2] + dt * aux->eps * aux->curl_H[2]; } } -void phy_sim_add_sources_E(phy_sim *sim) -{ - phy_field_info *fi = &sim->field_info; - float *E; - size_t x, y, z; - - - x = fi->dims[0]/2; - y = fi->dims[1]/2; - z = fi->dims[2]/2; - E = sim->fields[2].E + z * fi->zstr + y * fi->ystr + - x * fi->xstr; - - E[0] += exp(-40.0f * pow(sim->time - 3, 2)) * 50; - -/* - x = 10; - for (y = 1; y < fi->dims[1] - 1; y++) - for (z = 1; z < fi->dims[2] - 1; z++) { - E = sim->fields[2].E + z * fi->zstr + y * fi->ystr + - x * fi->xstr; - - E[0] += sin(sim->time * 3) * 1; - E[1] += 0; - E[2] += 0; - }*/ -} - -void phy_sim_add_sources_H(phy_sim *sim) -{ - /*phy_field_info *fi = &sim->field_info; - float *H; - size_t x, y, z; - - x = 0; - for (y = 1; y < fi->dims[1] - 1; y++) - for (z = 1; z < fi->dims[2] - 1; z++) { - H = sim->fields[2].H + z * fi->zstr + y * fi->ystr + - x * fi->xstr; - - H[0] += 0; - H[1] += 0; - H[2] += sin(sim->time + sim->time_delta / 2) / IMPEDANCE_ZERO * 100; - }*/ -} - - void phy_sim_step(phy_sim *sim) { phy_field_em *fields = sim->fields, tmp; int64_t start_time; @@ -383,10 +375,8 @@ void phy_sim_step(phy_sim *sim) { phy_sim_update_aux_H(sim); phy_sim_update_em_H(sim); - phy_sim_add_sources_H(sim); phy_sim_update_aux_E(sim); phy_sim_update_em_E(sim); - phy_sim_add_sources_E(sim); SDL_mutexP(sim->rotate_lock); @@ -406,17 +396,23 @@ void phy_sim_step(phy_sim *sim) { void phy_sim_debug(phy_sim *sim, size_t *coords) { phy_field_info *fi = &sim->field_info; - size_t offs1; phy_field_aux_point *aux; + float *E, *H, Emag, Hmag; + + aux = sim->aux + coords[2] * fi->zstr1 + coords[1] * fi->ystr1 + + coords[0] * fi->xstr1; - offs1 = coords[2] * fi->zstr1 + coords[1] * fi->ystr1 + - coords[0] * fi->xstr1; + E = sim->fields[1].E + coords[2] * fi->zstr + coords[1] * fi->ystr + + coords[0] * fi->xstr; + H = sim->fields[1].H + coords[2] * fi->zstr + coords[1] * fi->ystr + + coords[0] * fi->xstr; - aux = sim->aux + offs1; + Emag = sqrt(sq(E[0]) + sq(E[1]) + sq(E[2])); + Hmag = sqrt(sq(H[0]) + sq(H[1]) + sq(H[2])); printf("===[ phy_sim_debug ]===\n"); - printf("x = [%zu %zu %zu], offs1=%zu, aux=%p\n", coords[0], coords[1], - coords[2], offs1, aux); + printf("x = [%zu %zu %zu], E=%p, H=%p, aux=%p\n", coords[0], coords[1], + coords[2], E, H, aux); printf("eps = %f, mu = %f\n", aux->eps, aux->mu); printf("CE = [%E %E %E]\n", aux->curl_E[0], aux->curl_E[1], aux->curl_E[2]); @@ -434,6 +430,7 @@ void phy_sim_debug(phy_sim *sim, size_t *coords) aux->uec_H[0][0], aux->uec_H[0][1], aux->uec_H[0][2], aux->uec_H[0][3], aux->uec_H[1][0], aux->uec_H[1][1], aux->uec_H[1][2], aux->uec_H[1][3], aux->uec_H[2][0], aux->uec_H[2][1], aux->uec_H[2][2], aux->uec_H[2][3]); + printf("Emag = %E, Hmag = %E, Z = %E\n", Emag, Hmag, Emag/Hmag); printf("=======================\n"); } diff --git a/src/physics.h b/src/physics.h index 09ee23e..798fa8c 100644 --- a/src/physics.h +++ b/src/physics.h @@ -33,6 +33,11 @@ typedef struct { vec3_t int_curl_E, int_curl_H; } phy_field_aux_point; +typedef struct { + size_t z; + float Px, Py; +} phy_source_xyplane; + enum { PHY_CMD_QUIT, PHY_CMD_PAUSE, @@ -49,6 +54,8 @@ typedef struct { size_t pml_widths[3]; float time, time_delta; + phy_source_xyplane source_xyplane; + // UI stuff bool running; itc_chan ctl; diff --git a/src/renderer.c b/src/renderer.c index 94f0333..52f4208 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 = IMPEDANCE_ZERO; + scale = 1; //scale = IMPEDANCE_ZERO; } if (frac < 0.0f) diff --git a/src/ui.c b/src/ui.c index 2674756..4c02e53 100644 --- a/src/ui.c +++ b/src/ui.c @@ -240,18 +240,6 @@ void ui_event_window(SDL_Event *event, ui_window *uiw) sv->xsection_flags |= XSECTION_YZ; break; - case SDLK_w: - sv->xsection_frac += 0.1f; - if (sv->xsection_frac > 1.0f) - sv->xsection_frac = 1.0f; - break; - - case SDLK_s: - sv->xsection_frac -= 0.1f; - if (sv->xsection_frac < 0.0f) - sv->xsection_frac = 0.0f; - break; - case SDLK_o: ui_infof(uiw, "Symulacja wstrzymana"); itc_chan_push(&sv->sim->ctl, PHY_CMD_PAUSE, NULL); @@ -575,7 +563,7 @@ void ui_draw_window(ui_window *uiw) v2_set(sv->origin, 0, 0); v2_set(sv->scale, 1, 1); - sv->xsection_flags = XSECTION_XY | XSECTION_E; + sv->xsection_flags = XSECTION_XZ | XSECTION_E; sv->xsection_frac = 0.5f; sv->select_valid = false; -- cgit