summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorPaweł Redman <pawel.redman@gmail.com>2016-04-02 21:10:03 +0200
committerPaweł Redman <pawel.redman@gmail.com>2016-04-02 21:10:03 +0200
commitf1e955e7e9719aa82d1b32dac2bce62f26f75df7 (patch)
treeeee5d1cd4a68c9cb065384bbd3f5584bd4445bd5
parente40ab41ae226e257d57a73c69d9bdd8d19f64cb5 (diff)
Auxiliary fields (WIP).
-rw-r--r--src/common.h28
-rw-r--r--src/physics.c227
-rw-r--r--src/physics.h22
-rw-r--r--src/renderer.c27
-rw-r--r--src/ui.c18
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;
+ }
}
}
}