From f1e955e7e9719aa82d1b32dac2bce62f26f75df7 Mon Sep 17 00:00:00 2001 From: Paweł Redman Date: Sat, 2 Apr 2016 21:10:03 +0200 Subject: Auxiliary fields (WIP). --- src/common.h | 28 +++++++ src/physics.c | 227 ++++++++++++++++++++++++++++++++++++++++----------------- src/physics.h | 22 +++++- src/renderer.c | 27 ++++--- src/ui.c | 18 +++++ 5 files changed, 241 insertions(+), 81 deletions(-) diff --git a/src/common.h b/src/common.h index 8933749..743dea4 100644 --- a/src/common.h +++ b/src/common.h @@ -85,6 +85,34 @@ static inline void v2_mul(vec_t *R, vec_t *A, vec_t b) R[1] = A[1] * b; } +static inline void v3_add(vec_t *R, vec_t *A, vec_t *B) +{ + R[0] = A[0] + B[0]; + R[1] = A[1] + B[1]; + R[2] = A[2] + B[2]; +} + +static inline void v3_acc(vec_t *R, vec_t *A) +{ + R[0] += A[0]; + R[1] += A[1]; + R[2] += A[2]; +} + +static inline void v3_sub(vec_t *R, vec_t *A, vec_t *B) +{ + R[0] = A[0] - B[0]; + R[1] = A[1] - B[1]; + R[2] = A[2] - B[2]; +} + +static inline void v3_mul(vec_t *R, vec_t *A, vec_t b) +{ + R[0] = A[0] * b; + R[1] = A[1] * b; + R[2] = A[2] * b; +} + static inline void mst2_set_identity(mst2_t R) { R[0] = 0; diff --git a/src/physics.c b/src/physics.c index 7926478..a56ad97 100644 --- a/src/physics.c +++ b/src/physics.c @@ -16,6 +16,27 @@ void phy_sim_destroy(phy_sim *sim) } } +void phy_sim_reset_aux(phy_sim *sim) +{ + size_t x, y, z; + phy_field_info *fi = &sim->field_info; + + for (z = 0; z < fi->depth; z++) + for (y = 0; y < fi->height; y++) + for (x = 0; x < fi->width; x++) { + phy_field_aux_point *point; + + point = sim->aux + z * fi->zstr1 + y * fi->ystr1 + + x * fi->xstr1; + + point->eps = 1.0f / EPS_ZERO; + point->mu = -1.0f / MU_ZERO; + + v3_set(point->int_curl_E, 0, 0, 0); + v3_set(point->int_curl_H, 0, 0, 0); + } +} + void phy_sim_reset(phy_sim *sim) { size_t i; @@ -28,9 +49,7 @@ void phy_sim_reset(phy_sim *sim) memset(sim->fields[i].H, 0, fi->size * sizeof(float)); } - memset(sim->field_eps, 0, fi->size1 * sizeof(float)); - memset(sim->field_eps, 0, fi->size1 * sizeof(float)); - phy_sim_compute_const_fields(sim); + phy_sim_reset_aux(sim); sim->time = 0; sim->running = false; @@ -49,8 +68,11 @@ int phy_sim_create(phy_sim *sim) fi->width = 100; fi->height = 100; fi->depth = 100; - fi->spacing = 1.5e10f / max3(fi->width, fi->height, fi->depth); - sim->time_delta = 0.2; + fi->spacing = 10e+9f / max3(fi->width, fi->height, fi->depth); + sim->time_delta = 0.1; + + printf("phy_sim_create: Courant number is %f\n", + 3 * LIGHT_SPEED * sim->time_delta / fi->spacing); fi->xstr = 3; fi->ystr = fi->xstr * fi->width; @@ -71,9 +93,8 @@ int phy_sim_create(phy_sim *sim) } } - sim->field_eps = malloc(fi->size1 * sizeof(float)); - sim->field_mu = malloc(fi->size1 * sizeof(float)); - if (!sim->field_eps || !sim->field_mu) { + sim->aux = malloc(fi->size1 * sizeof(phy_field_aux_point)); + if (!sim->aux) { con_printf("phy_sim_create: out of memory\n"); goto error; } @@ -92,83 +113,119 @@ error: return 1; } -void phy_sim_compute_const_fields(phy_sim *sim) +void phy_sim_update_aux_H(phy_sim *sim) { - size_t x, y, z; phy_field_info *fi = &sim->field_info; + size_t x, y, z; - for (z = 0; z < fi->depth; z++) - for (y = 0; y < fi->height; y++) - for (x = 0; x < fi->width; x++) { - float *eps, *mu; + for (z = 1; z < fi->depth - 1; z++) + for (y = 1; y < fi->height - 1; y++) + for (x = 1; x < fi->width - 1; x++) { + size_t offs, offs_px, offs_py, offs_pz; + phy_field_aux_point *aux; + float *E = sim->fields[1].E, *H = sim->fields[1].H; + + offs = x * fi->xstr + y * fi->ystr + z * fi->zstr; + offs_px = offs + fi->xstr; + offs_py = offs + fi->ystr; + offs_pz = offs + fi->zstr; - eps = sim->field_eps + z * fi->zstr1 + y * fi->ystr1 + - x * fi->xstr1; - mu = sim->field_mu + z * fi->zstr1 + y * fi->ystr1 + - x * fi->xstr1; + aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 + + z * fi->zstr1; - *eps = 1.129409066837282e+11; - *mu = -7.957747154594767e+5; + aux->curl_E[0] = (E[offs_py+2] - E[offs+2] - + E[offs_pz+1] + E[offs+1]) / fi->spacing; + aux->curl_E[1] = (E[offs_pz+0] - E[offs+0] - + 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); } } -void phy_sim_step_E(phy_field_info *fi, float *mul_const, float dt, - float *En, float *Ep, float *H) +void phy_sim_update_aux_E(phy_sim *sim) { + phy_field_info *fi = &sim->field_info; size_t x, y, z; - float dx = fi->spacing; for (z = 1; z < fi->depth - 1; z++) for (y = 1; y < fi->height - 1; y++) for (x = 1; x < fi->width - 1; x++) { size_t offs, offs_nx, offs_ny, offs_nz; - size_t offs1; - float dExdt, dEydt, dEzdt; + phy_field_aux_point *aux; + float *E = sim->fields[1].E, *H = sim->fields[1].H; offs = x * fi->xstr + y * fi->ystr + z * fi->zstr; offs_nx = offs - fi->xstr; offs_ny = offs - fi->ystr; offs_nz = offs - fi->zstr; - offs1 = x * fi->xstr1 + y * fi->ystr1 + z * fi->zstr1; + aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 + + z * fi->zstr1; - dExdt = (H[offs+2] - H[offs_ny+2] - H[offs+1] + H[offs_nz+1]); - dEydt = (H[offs+0] - H[offs_nz+0] - H[offs+2] + H[offs_nx+2]); - dEzdt = (H[offs+1] - H[offs_nx+1] - H[offs+0] + H[offs_ny+0]); + aux->curl_H[0] = (H[offs+2] - H[offs_ny+2] - + H[offs+1] + H[offs_nz+1]) / fi->spacing; + aux->curl_H[1] = (H[offs+0] - H[offs_nz+0] - + 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; - En[offs+0] = Ep[offs+0] + dt * mul_const[offs1] * dExdt / dx; - En[offs+1] = Ep[offs+1] + dt * mul_const[offs1] * dEydt / dx; - En[offs+2] = Ep[offs+2] + dt * mul_const[offs1] * dEzdt / dx; + v3_acc(aux->int_E, E + offs); + v3_acc(aux->int_curl_H, aux->curl_H); } } -void phy_sim_step_H(phy_field_info *fi, float *mul_const, float dt, - float *Hn, float *Hp, float *E) +void phy_sim_update_em_H(phy_sim *sim) { + phy_field_info *fi = &sim->field_info; size_t x, y, z; - float dx = fi->spacing; + float *Hn, *Hp, dt; + + Hn = sim->fields[2].H; + Hp = sim->fields[0].H; + dt = sim->time_delta; for (z = 1; z < fi->depth - 1; z++) for (y = 1; y < fi->height - 1; y++) for (x = 1; x < fi->width - 1; x++) { - size_t offs, offs_px, offs_py, offs_pz; - size_t offs1; - float dHxdt, dHydt, dHzdt; + size_t offs; + phy_field_aux_point *aux; offs = x * fi->xstr + y * fi->ystr + z * fi->zstr; - offs_px = offs + fi->xstr; - offs_py = offs + fi->ystr; - offs_pz = offs + fi->zstr; + aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 + + z * fi->zstr1; + + Hn[offs+0] = Hp[offs+0] + dt * aux->eps * aux->curl_E[0]; + Hn[offs+1] = Hp[offs+1] + dt * aux->eps * aux->curl_E[1]; + Hn[offs+2] = Hp[offs+2] + dt * aux->eps * aux->curl_E[2]; + } +} + +void phy_sim_update_em_E(phy_sim *sim) +{ + phy_field_info *fi = &sim->field_info; + size_t x, y, z; + float *En, *Ep, dt; + + En = sim->fields[2].E; + Ep = sim->fields[0].E; + dt = sim->time_delta; - offs1 = x * fi->xstr1 + y * fi->ystr1 + z * fi->zstr1; + for (z = 1; z < fi->depth - 1; z++) + for (y = 1; y < fi->height - 1; y++) + for (x = 1; x < fi->width - 1; x++) { + size_t offs; + phy_field_aux_point *aux; - dHxdt = (E[offs_py+2] - E[offs+2] - E[offs_pz+1] + E[offs+1]); - dHydt = (E[offs_pz+0] - E[offs+0] - E[offs_px+2] + E[offs+2]); - dHzdt = (E[offs_px+1] - E[offs+1] - E[offs_py+0] + E[offs+0]); + offs = x * fi->xstr + y * fi->ystr + z * fi->zstr; + aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 + + z * fi->zstr1; - Hn[offs+0] = Hp[offs+0] + dt * mul_const[offs1] * dHxdt / dx; - Hn[offs+1] = Hp[offs+1] + dt * mul_const[offs1] * dHydt / dx; - Hn[offs+2] = Hp[offs+2] + dt * mul_const[offs1] * dHzdt / dx; + En[offs+0] = Ep[offs+0] + dt * aux->mu * aux->curl_H[0]; + En[offs+1] = Ep[offs+1] + dt * aux->mu * aux->curl_H[1]; + En[offs+2] = Ep[offs+2] + dt * aux->mu * aux->curl_H[2]; } } @@ -178,13 +235,16 @@ void phy_sim_add_sources_E(phy_sim *sim) float *E; size_t x, y, z; - x = fi->width / 2; - y = fi->height / 2; - z = fi->depth / 2; + x = 0; + for (y = 1; y < fi->height - 1; y++) + for (z = 1; z < fi->depth - 1; z++) { + E = sim->fields[2].E + z * fi->zstr + y * fi->ystr + + x * fi->xstr; - E = sim->fields[2].E + z * fi->zstr + y * fi->ystr + x * fi->xstr; - E[0] += sin(sim->time * 0.5) * - exp(-0.1f * pow(sim->time - 3.0f, 2)) * 50; + E[0] += sin(sim->time) * 100; + E[1] += 0; + E[2] += 0; + } } void phy_sim_add_sources_H(phy_sim *sim) @@ -193,27 +253,31 @@ void phy_sim_add_sources_H(phy_sim *sim) float *H; size_t x, y, z; - x = fi->width / 2; - y = fi->height / 2; - z = fi->depth / 2; + x = 0; + for (y = 1; y < fi->height - 1; y++) + for (z = 1; z < fi->depth - 1; z++) { + H = sim->fields[2].H + z * fi->zstr + y * fi->ystr + + x * fi->xstr; - H = sim->fields[2].H + z * fi->zstr + y * fi->ystr + x * fi->xstr; - H[1] -= exp(-0.1f * pow(sim->time - 5.0f - sim->time_delta / 2, 2)) * 5; + 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_info *fi = &sim->field_info; phy_field_em *fields = sim->fields, tmp; int64_t start_time; start_time = get_time(); - phy_sim_step_E(fi, sim->field_mu, sim->time_delta, - fields[2].E, fields[0].E, fields[1].H); - phy_sim_add_sources_E(sim); - phy_sim_step_H(fi, sim->field_eps, sim->time_delta, - fields[2].H, fields[0].H, fields[1].E); + 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); @@ -230,6 +294,31 @@ void phy_sim_step(phy_sim *sim) { SDL_mutexV(sim->rotate_lock); } +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; + + offs1 = coords[2] * fi->zstr1 + coords[1] * fi->ystr1 + + coords[0] * fi->xstr1; + + aux = sim->aux + offs1; + + printf("===[ phy_sim_debug ]===\n"); + printf("x = [%zu %zu %zu], offs1=%zu\n", coords[0], coords[1], + coords[2], offs1); + printf("CE = [%E %E %E]\n", aux->curl_E[0], aux->curl_E[1], + aux->curl_E[2]); + printf("CH = [%E %E Ef]\n", aux->curl_H[0], aux->curl_H[1], + aux->curl_H[2]); + printf("ICE = [%E %E %E]\n", aux->int_curl_E[0], aux->int_curl_E[1], + aux->int_curl_E[2]); + printf("ICH = [%E %E %E]\n", aux->int_curl_H[0], aux->int_curl_H[1], + aux->int_curl_H[2]); + printf("=======================\n"); +} + void phy_run_command(phy_sim *sim, int number, void *data) { switch (number) { @@ -249,12 +338,17 @@ void phy_run_command(phy_sim *sim, int number, void *data) case PHY_CMD_RESET: phy_sim_reset(sim); break; + + case PHY_CMD_DEBUG: + phy_sim_debug(sim, data); + free(data); + break; } } int phy_thread(phy_sim *sim) { - int rv, cmd_num; + int cmd_num; void *cmd_data; while (1) { @@ -286,5 +380,6 @@ cont_main_loop: out: printf("phy_thread: quitting\n"); + return 0; } diff --git a/src/physics.h b/src/physics.h index 38a6100..2557d3c 100644 --- a/src/physics.h +++ b/src/physics.h @@ -4,11 +4,15 @@ #include "common.h" #include "itc.h" +#define LIGHT_SPEED (299792458.0f) +#define MU_ZERO (4.0f * M_PI * 10e-7f) +#define EPS_ZERO (1.0f / (MU_ZERO * LIGHT_SPEED * LIGHT_SPEED)) +#define IMPEDANCE_ZERO (sqrt(MU_ZERO/EPS_ZERO)) + typedef struct { size_t width, height, depth; size_t zstr, ystr, xstr, size; // strides in no. of points (not bytes) - size_t zstr1, ystr1, xstr1, size1; // same as above but for associated - // scalar fields + size_t zstr1, ystr1, xstr1, size1; // same as above but for aux float spacing; // in meters } phy_field_info; @@ -17,18 +21,28 @@ typedef struct { float *H; } phy_field_em; +typedef struct { + float eps, mu; // permittivity and permeability + + vec3_t curl_E, curl_H; + + vec3_t int_E, int_H; + vec3_t int_curl_E, int_curl_H; +} phy_field_aux_point; + enum { PHY_CMD_QUIT, PHY_CMD_PAUSE, PHY_CMD_RESUME, PHY_CMD_STEP, - PHY_CMD_RESET + PHY_CMD_RESET, + PHY_CMD_DEBUG }; typedef struct { phy_field_info field_info; phy_field_em fields[3]; - float *field_eps, *field_mu; //permittivity and permeability + phy_field_aux_point *aux; float time, time_delta; // UI stuff diff --git a/src/renderer.c b/src/renderer.c index eee6e6d..9cf8db3 100644 --- a/src/renderer.c +++ b/src/renderer.c @@ -367,6 +367,7 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim, size_t width, height, x, y, z; uint8_t *pixels; int pitch; + float scale; if (xsection->frame_index && xsection->frame_index == sim->frame_index && @@ -385,10 +386,14 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim, height = fi->depth; } - if (flags & XSECTION_E) + if (flags & XSECTION_E) { field = sim->fields[1].E; - else + scale = 1; + } + else { field = sim->fields[1].H; + scale = IMPEDANCE_ZERO; + } if (frac < 0.0f) frac = 0.0f; @@ -416,9 +421,9 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim, point = field + z * fi->zstr + y * fi->ystr + x * fi->xstr; - pixel[0] = r_float_to_u8(point[0]); - pixel[1] = r_float_to_u8(point[1]); - pixel[2] = r_float_to_u8(point[2]); + pixel[0] = r_float_to_u8(point[0] * scale); + pixel[1] = r_float_to_u8(point[1] * scale); + pixel[2] = r_float_to_u8(point[2] * scale); } } else if (flags & XSECTION_XZ) { y = frac * fi->height; @@ -434,9 +439,9 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim, point = field + z * fi->zstr + y * fi->ystr + x * fi->xstr; - pixel[0] = r_float_to_u8(point[0]); - pixel[1] = r_float_to_u8(point[1]); - pixel[2] = r_float_to_u8(point[2]); + pixel[0] = r_float_to_u8(point[0] * scale); + pixel[1] = r_float_to_u8(point[1] * scale); + pixel[2] = r_float_to_u8(point[2] * scale); } } else { x = frac * fi->width; @@ -452,9 +457,9 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim, point = field + z * fi->zstr + y * fi->ystr + x * fi->xstr; - pixel[0] = r_float_to_u8(point[0]); - pixel[1] = r_float_to_u8(point[1]); - pixel[2] = r_float_to_u8(point[2]); + pixel[0] = r_float_to_u8(point[0] * scale); + pixel[1] = r_float_to_u8(point[1] * scale); + pixel[2] = r_float_to_u8(point[2] * scale); } } diff --git a/src/ui.c b/src/ui.c index fac56ce..f485cfe 100644 --- a/src/ui.c +++ b/src/ui.c @@ -258,6 +258,24 @@ void ui_event_window(SDL_Event *event, ui_window *uiw) sv->xsection_flags &= ~XSECTION_E; sv->xsection_flags |= XSECTION_H; break; + + case SDLK_z: + if (sv->select_valid) + { + phy_field_info *fi = &sv->sim->field_info; + size_t *coords; + + coords = calloc(3, sizeof(size_t)); + assert(coords); + + coords[0] = sv->select[0] * fi->width; + coords[1] = sv->select[1] * fi->height; + coords[2] = sv->select[2] * fi->depth; + + itc_chan_push(&sv->sim->ctl, PHY_CMD_DEBUG, + coords); + break; + } } } } -- cgit