summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--TODO2
-rw-r--r--src/common.h1
-rw-r--r--src/physics.c128
-rw-r--r--src/physics.h7
-rw-r--r--src/ui.c58
5 files changed, 153 insertions, 43 deletions
diff --git a/TODO b/TODO
index d5bb629..261e76f 100644
--- a/TODO
+++ b/TODO
@@ -1 +1,3 @@
+- unstable simulation
+- wrong impedance
itc_chan_destroy: free queue
diff --git a/src/common.h b/src/common.h
index 743dea4..9dc85e9 100644
--- a/src/common.h
+++ b/src/common.h
@@ -18,6 +18,7 @@
#define max3(a, b, c) (max2(max2(a, b), c))
#define min2(a, b) ((a) < (b) ? (a) : (b))
#define min3(a, b, c) (min2(min2(a, b), c))
+#define sq(x) ((x)*(x))
char *va(const char *fmt, ...);
int64_t get_time(void);
diff --git a/src/physics.c b/src/physics.c
index bbe2499..952e53a 100644
--- a/src/physics.c
+++ b/src/physics.c
@@ -16,8 +16,8 @@ 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)
+static void phy_calc_fc(float *fc, size_t *x, size_t *widths,
+ size_t *dims, float dt)
{
size_t i;
@@ -31,31 +31,88 @@ static float phy_calc_fake_cond(float *fake_cond, size_t *x, size_t *widths,
else
continue;
- fake_cond[i] = EPS_ZERO / (2 * dt) * pow(f, 3);
+ fc[i] = EPS_ZERO / (2 * dt) * pow(f, 3);
}
}
+static void phy_calc_uec_H(float *uec, float dt, float fc0, float fc1,
+ float fc2, float mu)
+{
+ float tmp;
+
+ tmp = 1.0f / dt + (fc1 + fc2) / (2 * EPS_ZERO) +
+ fc1 * fc2 * dt / sq(2.0f * EPS_ZERO);
+ uec[0] = 1.0f / tmp * (1.0f / dt - (fc1 + fc2) / (2 * EPS_ZERO)
+ - fc1 * fc2 * dt / sq(2.0f * EPS_ZERO));
+ uec[1] = -1.0f / tmp * LIGHT_SPEED / mu;
+ uec[2] = -1.0f / tmp * LIGHT_SPEED * dt / EPS_ZERO * fc0 / mu;
+ uec[3] = -1.0f / tmp * dt / sq(EPS_ZERO) * fc1 * fc2;
+}
+
+static void phy_calc_uec_E(float *uec, float dt, float fc0, float fc1,
+ float fc2, float eps)
+{
+ float tmp;
+
+ tmp = 1.0f / dt + (fc1 + fc2) / (2 * EPS_ZERO) +
+ fc1 * fc2 * dt / sq(2.0f * EPS_ZERO);
+ uec[0] = 1.0f / tmp * (1.0f / dt - (fc1 + fc2) / (2 * EPS_ZERO)
+ - fc1 * fc2 * dt / sq(2.0f * EPS_ZERO));
+ uec[1] = LIGHT_SPEED / tmp;
+ uec[2] = 1.0f / tmp * LIGHT_SPEED * dt * fc0 / EPS_ZERO;
+ 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;*/
+}
+
void phy_sim_reset_aux(phy_sim *sim)
{
- size_t x[3], i;
+ size_t x[3];
phy_field_info *fi = &sim->field_info;
+ float dt = sim->time_delta;
+
+ memset(sim->aux, 0, fi->size1 * sizeof(phy_field_aux_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;
+ vec3_t fc = {0, 0, 0};
aux = sim->aux + x[2] * fi->zstr1 + x[1] * fi->ystr1 +
x[0] * fi->xstr1;
- aux->eps = 1.0f / EPS_ZERO;
- aux->mu = -1.0f / MU_ZERO;
+ aux->eps = EPS_ZERO;
+ aux->mu = MU_ZERO;
v3_set(aux->int_curl_E, 0, 0, 0);
v3_set(aux->int_curl_H, 0, 0, 0);
- phy_calc_fake_cond(aux->fake_cond, x, sim->pml_widths,
- fi->dims, sim->time_delta);
+ 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);
}
}
@@ -96,9 +153,15 @@ int phy_sim_create(phy_sim *sim)
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;
+ #else
+ sim->pml_widths[0] = 0;
+ sim->pml_widths[1] = 0;
+ sim->pml_widths[2] = 0;
+ #endif
fi->xstr = 3;
fi->ystr = fi->xstr * fi->dims[0];
@@ -207,25 +270,30 @@ void phy_sim_update_em_H(phy_sim *sim)
{
phy_field_info *fi = &sim->field_info;
size_t x, y, z;
- float *Hn, *Hp, dt;
+ float *Hn, *Hp;
Hn = sim->fields[2].H;
Hp = sim->fields[0].H;
- dt = sim->time_delta;
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;
+ size_t offs, i;
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;
- 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];
+ for (i = 0; i < 3; i++)
+ Hn[offs+i] = aux->uec_H[i][0] * Hp[offs+i] +
+ 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];
}
}
@@ -233,25 +301,30 @@ 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;
+ float *En, *Ep;
En = sim->fields[2].E;
Ep = sim->fields[0].E;
- dt = sim->time_delta;
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;
+ size_t offs, i;
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;
- 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];
+ for (i = 0; i < 3; i++)
+ En[offs+i] = aux->uec_E[i][0] * Ep[offs+i] +
+ 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];
}
}
@@ -342,8 +415,9 @@ void phy_sim_debug(phy_sim *sim, size_t *coords)
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("x = [%zu %zu %zu], offs1=%zu, aux=%p\n", coords[0], coords[1],
+ coords[2], offs1, 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]);
printf("CH = [%E %E %E]\n", aux->curl_H[0], aux->curl_H[1],
@@ -352,8 +426,14 @@ void phy_sim_debug(phy_sim *sim, size_t *coords)
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("UEC_E = [[%E %E %E %E] [%E %E %E %E] [%E %E %E %E]]\n",
+ aux->uec_E[0][0], aux->uec_E[0][1], aux->uec_E[0][2], aux->uec_E[0][3],
+ aux->uec_E[1][0], aux->uec_E[1][1], aux->uec_E[1][2], aux->uec_E[1][3],
+ aux->uec_E[2][0], aux->uec_E[2][1], aux->uec_E[2][2], aux->uec_E[2][3]);
+ printf("UEC_H = [[%E %E %E %E] [%E %E %E %E] [%E %E %E %E]]\n",
+ 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("=======================\n");
}
diff --git a/src/physics.h b/src/physics.h
index 06c5d4f..09ee23e 100644
--- a/src/physics.h
+++ b/src/physics.h
@@ -5,7 +5,7 @@
#include "itc.h"
#define LIGHT_SPEED (299792458.0f)
-#define MU_ZERO (4.0f * M_PI * 10e-7f)
+#define MU_ZERO (4.0f * M_PI * 1e-7f)
#define EPS_ZERO (1.0f / (MU_ZERO * LIGHT_SPEED * LIGHT_SPEED))
#define IMPEDANCE_ZERO (sqrt(MU_ZERO/EPS_ZERO))
@@ -22,9 +22,10 @@ typedef struct {
} phy_field_em;
typedef struct {
- float eps, mu; // permittivity and permeability
+ float eps, mu; // for visualisation, not needed for stepping
- float fake_cond[3];
+ float uec_H[3][4];
+ float uec_E[3][4];
vec3_t curl_E, curl_H;
diff --git a/src/ui.c b/src/ui.c
index 4245603..3f2214a 100644
--- a/src/ui.c
+++ b/src/ui.c
@@ -24,7 +24,7 @@ struct {
typedef struct {
phy_sim *sim;
- bool dragging, selecting;
+ bool dragging, dragging_frac, selecting;
int xsection_flags;
float xsection_frac;
bool select_valid;
@@ -167,31 +167,54 @@ void ui_event_window(SDL_Event *event, ui_window *uiw)
case SDL_MOUSEBUTTONDOWN:
if (event->button.y > sv->margin_top &&
event->button.y < uiw->h - sv->margin_bottom) {
- if (event->button.button == SDL_BUTTON_LEFT)
+ switch (event->button.button) {
+ case SDL_BUTTON_LEFT:
sv->dragging = true;
-
- if (event->button.button == SDL_BUTTON_RIGHT)
+ break;
+ case SDL_BUTTON_RIGHT:
+ sv->dragging_frac = true;
+ break;
+ case SDL_BUTTON_MIDDLE:
sv->selecting = true;
+ break;
+ }
}
break;
case SDL_MOUSEBUTTONUP:
- if (event->button.button == SDL_BUTTON_LEFT)
+ switch (event->button.button) {
+ case SDL_BUTTON_LEFT:
sv->dragging = false;
- else if (event->button.button == SDL_BUTTON_RIGHT)
+ break;
+ case SDL_BUTTON_RIGHT:
+ sv->dragging_frac = false;
+ break;
+ case SDL_BUTTON_MIDDLE:
sv->selecting = false;
- break;
+ break;
+ }
case SDL_MOUSEMOTION:
uiw->mouse[0] = event->motion.x;
uiw->mouse[1] = event->motion.y;
- if (sv->dragging) {
+ if (sv->dragging || sv->dragging_frac) {
vec2_t delta_v;
v2_set(delta_v, event->motion.xrel, event->motion.yrel);
v2_div_mst2_nt(delta_v, delta_v, sv->tf_v2s);
- v2_add(sv->origin, sv->origin, delta_v);
+
+ if (sv->dragging)
+ v2_add(sv->origin, sv->origin, delta_v);
+
+ if (sv->dragging_frac) {
+ sv->xsection_frac += delta_v[1];
+
+ if (sv->xsection_frac > 1.0f)
+ sv->xsection_frac = 1.0f;
+ else if (sv->xsection_frac < 0.0f)
+ sv->xsection_frac = 0.0f;
+ }
}
break;
@@ -407,14 +430,17 @@ void ui_draw_simview_selection(ui_window *uiw)
char *ui_draw_simview_top_bar_text(ui_simview *sv)
{
int flags = sv->xsection_flags;
+ char *rv;
+
+ rv = va("Przekrój %s płaszczyzną %s=%f",
+ (flags & XSECTION_E ? "E" :
+ "H"),
+ (flags & XSECTION_XY ? "Z" :
+ flags & XSECTION_XZ ? "Y" :
+ "X"),
+ sv->xsection_frac);
- va("Przekrój %s płaszczyzną %s=%f",
- (flags & XSECTION_E ? "E" :
- "H"),
- (flags & XSECTION_XY ? "Z" :
- flags & XSECTION_XZ ? "Y" :
- "X"),
- sv->xsection_frac);
+ return rv;
}
void ui_draw_simview_bars(ui_window *uiw, float dt)
generated by cgit (git 2.34.1) at 2024-12-26 05:39:10 +0000