From 4bf61f911194bb0b40d302df8d7dfff2fda19892 Mon Sep 17 00:00:00 2001 From: Paweł Redman Date: Wed, 30 Mar 2016 08:25:04 +0200 Subject: Initial commit. --- src/physics.c | 244 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 244 insertions(+) create mode 100644 src/physics.c (limited to 'src/physics.c') diff --git a/src/physics.c b/src/physics.c new file mode 100644 index 0000000..be6a426 --- /dev/null +++ b/src/physics.c @@ -0,0 +1,244 @@ +#include "common.h" + +void phy_sim_destroy(phy_sim *sim) { + size_t i; + + SDL_DestroyMutex(sim->rotate_lock); + + for (i = 0; i < 3; i++) { + free(sim->fields[i].E); + free(sim->fields[i].H); + } +} + +int phy_sim_create(phy_sim *sim) { + size_t i; + phy_field_info *fi; + + memset(sim, 0, sizeof(phy_sim)); + + fi = &sim->field_info; + + fi->width = 30; + fi->height = 30; + fi->depth = 30; + fi->spacing = 1.5e10 / max3(fi->width, fi->height, fi->depth); + + fi->xstr = 3; + fi->ystr = fi->xstr * fi->width; + fi->zstr = fi->ystr * fi->height; + fi->size = fi->zstr * fi->depth; + + fi->xstr1 = 1; + fi->ystr1 = fi->xstr1 * fi->width; + fi->zstr1 = fi->ystr1 * fi->height; + fi->size1 = fi->zstr1 * fi->depth; + + for (i = 0; i < 3; i++) { + sim->fields[i].E = calloc(fi->size, sizeof(float)); + sim->fields[i].H = calloc(fi->size, sizeof(float)); + if (!sim->fields[i].E || !sim->fields[i].H) { + con_printf("phy_sim_create: out of memory\n"); + goto error; + } + } + + sim->field_eps = calloc(fi->size1, sizeof(float)); + sim->field_mu = calloc(fi->size1, sizeof(float)); + + phy_sim_compute_const_fields(sim); + + sim->rotate_lock = SDL_CreateMutex(); + assert(sim->rotate_lock); + + return 0; + +error: + phy_sim_destroy(sim); + return 1; +} + +void phy_sim_compute_const_fields(phy_sim *sim) +{ + size_t x, y, z; + phy_field_info *fi = &sim->field_info; + size_t border_width; + + border_width = 3; + + for (z = 0; z < fi->depth; z++) + for (y = 0; y < fi->height; y++) + for (x = 0; x < fi->width; x++) { + float *eps, *mu; + float border_fx, border_fy, border_fz, border_f; + + 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; + + *eps = 1.129409066837282e+11; + *mu = -7.957747154594767e+5; + +/* + if (x < border_width) { + border_fx = (float)x / border_width; + } else if (x >= fi->width - border_width) { + border_fx = 1.0f - ((float)x - fi->width + border_width + 1) / border_width; + } else { + border_fx = 1.0f; + } + + if (y < border_width) { + border_fy = (float)y / border_width; + } else if (y >= fi->height - border_width) { + border_fy = 1.0f - ((float)y - fi->height + border_width + 1) / border_width; + } else { + border_fy = 1.0f; + } + + if (z < border_width) { + border_fz = (float)y / border_width; + } else if (z >= fi->depth - border_width) { + border_fz = 1.0f - ((float)z - fi->depth + border_width + 1) / border_width; + } else { + border_fz = 1.0f; + } + + border_f = min3(border_fx, border_fy, border_fz); + + *eps *= border_f; + *mu *= border_f;*/ + } +} + +void phy_sim_step_1(phy_field_info *fi, float dt, float *mul_const, + const float *L0, float *L2, const float *R) +{ + size_t x, y, z; + +#define AT(F, x, y, z) ((F) + (z) * fi->zstr + (y) * fi->ystr + (x) * fi->xstr) +#define AT1(F, x, y, z) ((F)[(z) * fi->zstr1 + (y) * fi->ystr1 + \ + (x) * fi->xstr1]) + + // X + for (z = 1; z < fi->depth - 1; z++) + for (y = 1; y < fi->height - 1; y++) + for (x = 0; x < fi->width; x++) { + float dRzdy, dRydz, dLxdt; + + dRzdy = AT(R, x, y + 1, z)[2] - AT(R, x, y - 1, z)[2]; + dRzdy /= 2 * fi->spacing; + + dRydz = AT(R, x, y, z + 1)[1] - AT(R, x, y, z - 1)[1]; + dRydz /= 2 * fi->spacing; + + dLxdt = AT1(mul_const, x, y, z) * (dRzdy - dRydz); + + AT(L2, x, y, z)[0] = AT(L0, x, y, z)[0] + dLxdt * dt; + } + + // Y + for (z = 1; z < fi->depth - 1; z++) + for (y = 0; y < fi->height; y++) + for (x = 1; x < fi->width - 1; x++) { + float dRxdz, dRzdx, dLydt; + + dRxdz = AT(R, x, y, z + 1)[0] - AT(R, x, y, z - 1)[0]; + dRxdz /= 2 * fi->spacing; + + dRzdx = AT(R, x + 1, y, z)[2] - AT(R, x - 1, y, z)[2]; + dRzdx /= 2 * fi->spacing; + + dLydt = AT1(mul_const, x, y, z) * (dRxdz - dRzdx); + + AT(L2, x, y, z)[1] = AT(L0, x, y, z)[1] + dLydt * dt; + } + + // Z + for (z = 0; z < fi->depth; z++) + for (y = 1; y < fi->height - 1; y++) + for (x = 1; x < fi->width - 1; x++) { + float dRydx, dRxdy, dLzdt; + + dRydx = AT(R, x + 1, y, z)[1] - AT(R, x - 1, y, z)[1]; + dRydx /= 2 * fi->spacing; + + dRxdy = AT(R, x, y + 1, z)[0] - AT(R, x, y - 1, z)[0]; + dRxdy /= 2 * fi->spacing; + + dLzdt = AT1(mul_const, x, y, z) * (dRydx - dRxdy); + + AT(L2, x, y, z)[2] = AT(L0, x, y, z)[2] + dLzdt * dt; + } +#undef AT +#undef AT1 +} + +void phy_sim_add_sources(phy_sim *sim) +{ + 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 *E, *H, fx, fy, fz, chuj; + + fx = (float)x / fi->width * 2 - 1; + fy = (float)y / fi->height * 2 - 1; + fz = (float)z / fi->depth * 2 - 1; + + E = sim->fields[2].E + z * fi->zstr + y * fi->ystr + x * fi->xstr; + H = sim->fields[2].H + z * fi->zstr + y * fi->ystr + x * fi->xstr; + + chuj = (0.3 - sqrt(fx * fx + fy * fy + fz * fz)) / 0.3; + if (chuj < 0) + chuj = 0; + + chuj *= sin(sim->time) * 10; + //chuj *= exp(-10.0f * pow(sim->time - 0.5, 2)) * 2; + E[0] += chuj; + } +} + +void phy_sim_step(phy_sim *sim, float dt) { + phy_field_info *fi = &sim->field_info; + phy_field_em *fields = sim->fields; + + size_t i, size_in_bytes; + + phy_sim_add_sources(sim); + + size_in_bytes = fi->size * sizeof(float); + + SDL_mutexP(sim->rotate_lock); + for (i = 0; i < 2; i++) { + // TODO: avoid copying, just switch around pointers + memcpy(fields[i].E, fields[i + 1].E, size_in_bytes); + memcpy(fields[i].H, fields[i + 1].H, size_in_bytes); + } + SDL_mutexV(sim->rotate_lock); + + phy_sim_step_1(fi, dt, sim->field_eps, + fields[0].E, fields[2].E, fields[1].H); + phy_sim_step_1(fi, dt, sim->field_mu, + fields[0].H, fields[2].H, fields[1].E); + + sim->time += dt; +} + +/* +void phy_control(phy_command cmd, void *data) + +} */ + +int phy_thread(phy_sim *sim) +{ + while (1) { + phy_sim_step(sim, 0.1); + } + + printf("phy_thread: quitting\n"); +} + -- cgit