summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/physics.c209
-rw-r--r--src/physics.h7
-rw-r--r--src/renderer.c2
-rw-r--r--src/ui.c14
4 files changed, 112 insertions, 120 deletions
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;
generated by cgit (git 2.34.1) at 2024-12-26 05:37:41 +0000