summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/physics.c140
-rw-r--r--src/physics.h5
-rw-r--r--src/renderer.c42
-rw-r--r--src/ui.c12
4 files changed, 120 insertions, 79 deletions
diff --git a/src/physics.c b/src/physics.c
index a56ad97..bbe2499 100644
--- a/src/physics.c
+++ b/src/physics.c
@@ -1,4 +1,4 @@
-#include "common.h"
+#include "common.h"
#include "itc.h"
#include "physics.h"
@@ -16,24 +16,46 @@ void phy_sim_destroy(phy_sim *sim)
}
}
+static float phy_calc_fake_cond(float *fake_cond, size_t *x, size_t *widths,
+ size_t *dims, float dt)
+{
+ 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];
+ else
+ continue;
+
+ fake_cond[i] = EPS_ZERO / (2 * dt) * pow(f, 3);
+ }
+}
+
void phy_sim_reset_aux(phy_sim *sim)
{
- size_t x, y, z;
+ size_t x[3], i;
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;
+ for (x[2] = 1; x[2] < fi->dims[2] - 1; x[2]++)
+ 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;
+
+ aux = sim->aux + x[2] * fi->zstr1 + x[1] * fi->ystr1 +
+ x[0] * fi->xstr1;
- point = sim->aux + z * fi->zstr1 + y * fi->ystr1 +
- x * fi->xstr1;
+ aux->eps = 1.0f / EPS_ZERO;
+ aux->mu = -1.0f / MU_ZERO;
- point->eps = 1.0f / EPS_ZERO;
- point->mu = -1.0f / MU_ZERO;
+ v3_set(aux->int_curl_E, 0, 0, 0);
+ v3_set(aux->int_curl_H, 0, 0, 0);
- v3_set(point->int_curl_E, 0, 0, 0);
- v3_set(point->int_curl_H, 0, 0, 0);
+ phy_calc_fake_cond(aux->fake_cond, x, sim->pml_widths,
+ fi->dims, sim->time_delta);
}
}
@@ -65,24 +87,28 @@ int phy_sim_create(phy_sim *sim)
fi = &sim->field_info;
- fi->width = 100;
- fi->height = 100;
- fi->depth = 100;
- fi->spacing = 10e+9f / max3(fi->width, fi->height, fi->depth);
+ 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;
printf("phy_sim_create: Courant number is %f\n",
3 * LIGHT_SPEED * sim->time_delta / fi->spacing);
+ sim->pml_widths[0] = fi->dims[0] / 5;
+ sim->pml_widths[1] = fi->dims[1] / 5;
+ sim->pml_widths[2] = fi->dims[2] / 5;
+
fi->xstr = 3;
- fi->ystr = fi->xstr * fi->width;
- fi->zstr = fi->ystr * fi->height;
- fi->size = fi->zstr * fi->depth;
+ fi->ystr = fi->xstr * fi->dims[0];
+ fi->zstr = fi->ystr * fi->dims[1];
+ fi->size = fi->zstr * fi->dims[2];
fi->xstr1 = 1;
- fi->ystr1 = fi->xstr1 * fi->width;
- fi->zstr1 = fi->ystr1 * fi->height;
- fi->size1 = fi->zstr1 * fi->depth;
+ fi->ystr1 = fi->xstr1 * fi->dims[0];
+ fi->zstr1 = fi->ystr1 * fi->dims[1];
+ fi->size1 = fi->zstr1 * fi->dims[2];
for (i = 0; i < 3; i++) {
sim->fields[i].E = malloc(fi->size * sizeof(float));
@@ -118,9 +144,9 @@ void phy_sim_update_aux_H(phy_sim *sim)
phy_field_info *fi = &sim->field_info;
size_t x, y, z;
- for (z = 1; z < fi->depth - 1; z++)
- for (y = 1; y < fi->height - 1; y++)
- for (x = 1; x < fi->width - 1; x++) {
+ for (z = 1; z < fi->dims[2] - 1; z++)
+ for (y = 1; y < fi->dims[1] - 1; y++)
+ for (x = 1; x < fi->dims[0] - 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;
@@ -150,9 +176,9 @@ void phy_sim_update_aux_E(phy_sim *sim)
phy_field_info *fi = &sim->field_info;
size_t x, y, z;
- for (z = 1; z < fi->depth - 1; z++)
- for (y = 1; y < fi->height - 1; y++)
- for (x = 1; x < fi->width - 1; x++) {
+ for (z = 1; z < fi->dims[2] - 1; z++)
+ for (y = 1; y < fi->dims[1] - 1; y++)
+ for (x = 1; x < fi->dims[0] - 1; x++) {
size_t offs, offs_nx, offs_ny, offs_nz;
phy_field_aux_point *aux;
float *E = sim->fields[1].E, *H = sim->fields[1].H;
@@ -187,19 +213,19 @@ void phy_sim_update_em_H(phy_sim *sim)
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++) {
+ for (z = 1; z < fi->dims[2] - 1; z++)
+ for (y = 1; y < fi->dims[1] - 1; y++)
+ for (x = 1; x < fi->dims[0] - 1; x++) {
size_t offs;
phy_field_aux_point *aux;
offs = x * fi->xstr + y * fi->ystr + z * fi->zstr;
aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 +
- z * fi->zstr1;
+ 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];
+ 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];
}
}
@@ -213,9 +239,9 @@ void phy_sim_update_em_E(phy_sim *sim)
Ep = sim->fields[0].E;
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++) {
+ for (z = 1; z < fi->dims[2] - 1; z++)
+ for (y = 1; y < fi->dims[1] - 1; y++)
+ for (x = 1; x < fi->dims[0] - 1; x++) {
size_t offs;
phy_field_aux_point *aux;
@@ -223,9 +249,9 @@ void phy_sim_update_em_E(phy_sim *sim)
aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 +
z * fi->zstr1;
- 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];
+ 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];
}
}
@@ -235,34 +261,44 @@ void phy_sim_add_sources_E(phy_sim *sim)
float *E;
size_t x, y, z;
- x = 0;
- for (y = 1; y < fi->height - 1; y++)
- for (z = 1; z < fi->depth - 1; 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) * 100;
+ 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;
+ /*phy_field_info *fi = &sim->field_info;
float *H;
size_t x, y, z;
x = 0;
- for (y = 1; y < fi->height - 1; y++)
- for (z = 1; z < fi->depth - 1; z++) {
+ 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;
- }
+ }*/
}
@@ -310,12 +346,14 @@ void phy_sim_debug(phy_sim *sim, size_t *coords)
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],
+ printf("CH = [%E %E %E]\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("FC = [%E %E %E]\n", aux->fake_cond[0], aux->fake_cond[1],
+ aux->fake_cond[2]);
printf("=======================\n");
}
diff --git a/src/physics.h b/src/physics.h
index 2557d3c..06c5d4f 100644
--- a/src/physics.h
+++ b/src/physics.h
@@ -10,7 +10,7 @@
#define IMPEDANCE_ZERO (sqrt(MU_ZERO/EPS_ZERO))
typedef struct {
- size_t width, height, depth;
+ size_t dims[3];
size_t zstr, ystr, xstr, size; // strides in no. of points (not bytes)
size_t zstr1, ystr1, xstr1, size1; // same as above but for aux
float spacing; // in meters
@@ -24,6 +24,8 @@ typedef struct {
typedef struct {
float eps, mu; // permittivity and permeability
+ float fake_cond[3];
+
vec3_t curl_E, curl_H;
vec3_t int_E, int_H;
@@ -43,6 +45,7 @@ typedef struct {
phy_field_info field_info;
phy_field_em fields[3];
phy_field_aux_point *aux;
+ size_t pml_widths[3];
float time, time_delta;
// UI stuff
diff --git a/src/renderer.c b/src/renderer.c
index 9cf8db3..94f0333 100644
--- a/src/renderer.c
+++ b/src/renderer.c
@@ -376,14 +376,14 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
return 0; // nothing's changed
if (flags & XSECTION_XY) {
- width = fi->width;
- height = fi->height;
+ width = fi->dims[0];
+ height = fi->dims[1];
} else if (flags & XSECTION_XZ) {
- width = fi->width;
- height = fi->depth;
+ width = fi->dims[0];
+ height = fi->dims[2];
} else {
- width = fi->height;
- height = fi->depth;
+ width = fi->dims[1];
+ height = fi->dims[2];
}
if (flags & XSECTION_E) {
@@ -408,12 +408,12 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
SDL_LockTexture(xsection->texture, NULL, (void**)&pixels, &pitch);
if (flags & XSECTION_XY) {
- z = frac * fi->depth;
- if (z >= fi->depth)
- z = fi->depth - 1;
+ z = frac * fi->dims[2];
+ if (z >= fi->dims[2])
+ z = fi->dims[2] - 1;
- for (y = 0; y < fi->height; y++)
- for (x = 0; x < fi->width; x++) {
+ for (y = 0; y < fi->dims[1]; y++)
+ for (x = 0; x < fi->dims[0]; x++) {
uint8_t *pixel;
float *point;
@@ -426,12 +426,12 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
pixel[2] = r_float_to_u8(point[2] * scale);
}
} else if (flags & XSECTION_XZ) {
- y = frac * fi->height;
- if (y >= fi->height)
- y = fi->height - 1;
+ y = frac * fi->dims[1];
+ if (y >= fi->dims[1])
+ y = fi->dims[1] - 1;
- for (z = 0; z < fi->depth; z++)
- for (x = 0; x < fi->width; x++) {
+ for (z = 0; z < fi->dims[2]; z++)
+ for (x = 0; x < fi->dims[0]; x++) {
uint8_t *pixel;
float *point;
@@ -444,12 +444,12 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
pixel[2] = r_float_to_u8(point[2] * scale);
}
} else {
- x = frac * fi->width;
- if (x >= fi->width)
- x = fi->width - 1;
+ x = frac * fi->dims[0];
+ if (x >= fi->dims[0])
+ x = fi->dims[0] - 1;
- for (z = 0; z < fi->depth; z++)
- for (y = 0; y < fi->height; y++) {
+ for (z = 0; z < fi->dims[2]; z++)
+ for (y = 0; y < fi->dims[1]; y++) {
uint8_t *pixel;
float *point;
diff --git a/src/ui.c b/src/ui.c
index f485cfe..4245603 100644
--- a/src/ui.c
+++ b/src/ui.c
@@ -268,9 +268,9 @@ void ui_event_window(SDL_Event *event, ui_window *uiw)
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;
+ coords[0] = sv->select[0] * fi->dims[0];
+ coords[1] = sv->select[1] * fi->dims[1];
+ coords[2] = sv->select[2] * fi->dims[2];
itc_chan_push(&sv->sim->ctl, PHY_CMD_DEBUG,
coords);
@@ -464,9 +464,9 @@ void ui_draw_simview_bars(ui_window *uiw, float dt)
r_draw_rect(uiw->rw, 0, uiw->h - 4 * theme.font_size,
uiw->w, 3 * theme.font_size, theme.color_main);
- x = sv->select[0] * fi->width;
- y = sv->select[1] * fi->height;
- z = sv->select[2] * fi->depth;
+ x = sv->select[0] * fi->dims[0];
+ y = sv->select[1] * fi->dims[1];
+ z = sv->select[2] * fi->dims[2];
offs = z * fi->zstr + y * fi->ystr + x * fi->xstr;
r_draw_text(uiw->rw, 0, uiw->h - 4 * theme.font_size,
generated by cgit (git 2.34.1) at 2025-11-04 06:41:50 +0000