summaryrefslogtreecommitdiff
path: root/src/physics.c
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 /src/physics.c
parente40ab41ae226e257d57a73c69d9bdd8d19f64cb5 (diff)
Auxiliary fields (WIP).
Diffstat (limited to 'src/physics.c')
-rw-r--r--src/physics.c227
1 files changed, 161 insertions, 66 deletions
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;
}