Finish CSG, permittivity visualisation.

7 files changed, 223 insertions, 96 deletions

diff --git a/src/common.h b/src/common.h
index 86d98ad..b6e90d2 100644
--- a/src/common.h
+++ b/src/common.h
@@ -178,4 +178,11 @@ static inline void mst2_dump(mst2_t M)
 	printf("[%f\t%f\t%f]\n", 0.0, 0.0, 1.0);
 }
 
+static inline float remap(float x, float a0, float b0, float a1, float b1)
+{
+	return (x - a0) * (b1 - a1) / (b0 - a0) + a1;
+}
+
+#define DEG2RAD(x) ((x) / (180.0f / M_PI))
+
 #endif // _COMMON_H
diff --git a/src/console.c b/src/console.c
index d415500..9c7cebb 100644
--- a/src/console.c
+++ b/src/console.c
@@ -155,6 +155,18 @@ char *cmd_load(char *arg)
 	return va("Załadowano plik %s", arg);
 }
 
+char *cmd_csg_clear(char *arg)
+{
+	itc_chan_push(&con_sim->ctl, PHY_CMD_CSG_CLEAR, NULL);
+	return "Wyzerowano przenikalność elektryczną";
+}
+
+char *cmd_csg_push(char *arg)
+{
+	itc_chan_push(&con_sim->ctl, PHY_CMD_CSG_PUSH, NULL);
+	return "Wykonano operację na przenikalności elektrycznej";
+}
+
 char *cmd_get(char *arg)
 {
 	if (!arg)
@@ -205,6 +217,8 @@ char *cmd_zero(char *arg)
 
 
 static const con_cmd con_cmds[] = {
+	{"csg-clear", cmd_csg_clear},
+	{"csg-push", cmd_csg_push},   
 	{"get", cmd_get},
 	{"getf", cmd_getf},
 	{"getz", cmd_getz},
diff --git a/src/physics.c b/src/physics.c
index fedc4c4..c7cc477 100644
--- a/src/physics.c
+++ b/src/physics.c
@@ -55,9 +55,48 @@ static void phy_calc_uec_E(float *uec, float dt, float fc0, float fc1,
 	uec[3] *= 1.0f / eps;
 } 
 
-static float remap(float x, float a0, float b0, float a1, float b1)
+
+static float phy_sim_csg_subsample(phy_csg_step *stack, float *x)
+{
+	phy_csg_step *step;
+	float rv = 1.0f;
+
+	for (step = stack; step; step = step->next) {
+		switch (step->shape) {
+		case CSG_SPHERE:
+			if (sq(x[0] - step->x[0]) + sq(x[1] - step->x[1]) +
+			    sq(x[2] - step->x[2]) > sq(step->d[0]))
+				goto next_step;
+			break;
+		}
+
+		rv = step->value;
+
+	next_step:
+		;
+	}
+
+	return rv;
+}
+
+static float phy_sim_csg_sample(phy_csg_step *stack, size_t *x)
 {
-	return (x - a0) * (b1 - a1) / (b0 - a0) + a1;
+	vec3_t xf = {x[0], x[1], x[2]};
+	size_t dx, dy, dz;
+	float total;
+
+	for (dz = 0; dz < 2; dz++)
+	for (dy = 0; dy < 2; dy++)
+	for (dx = 0; dx < 2; dx++) {
+		vec3_t xf;
+
+		v3_set(xf, x[0] + 0.5f * dx, x[1] + 0.5f * dy,
+		       x[2] + 0.5f * dz);
+
+		total += phy_sim_csg_subsample(stack, xf);
+	}
+	
+	return total / 8.0f;
 }
 
 static void phy_sim_zero_aux(phy_sim *sim)
@@ -68,6 +107,8 @@ static void phy_sim_zero_aux(phy_sim *sim)
 
 	memset(sim->aux, 0, fi->size1 * sizeof(phy_field_aux_point));
 
+	sim->eps_min = sim->eps_max = 1.0f;
+
 	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]++) {
@@ -77,23 +118,14 @@ static void phy_sim_zero_aux(phy_sim *sim)
 		aux = sim->aux + x[2] * fi->zstr1 + x[1] * fi->ystr1 +
 			x[0] * fi->xstr1;
 
-		aux->eps = 1;
-		aux->mu = 1;
-
-		if (0){
-			float cx, cy, cz, cr, d;
+		aux->eps = phy_sim_csg_sample(sim->csg_stack, x);
 
-			cx = 25, cy = 25, cz = 65;
-			cr = 10; 
+		if (aux->eps > sim->eps_max)
+			sim->eps_max = aux->eps;
+		if (aux->eps < sim->eps_min)
+			sim->eps_min = aux->eps;
 
-			d = sqrt(sq(x[0] - cx) + sq(x[1] - cy) + sq(x[2] - cz));
-
-			if (d < cr) {
-				aux->eps = remap(d, cr, cr * 0.8, 1, 2);
-				if (aux->eps > 2)
-					aux->eps = 2;
-			}
-		}
+		aux->mu = 1;
 
 		v3_set(aux->int_curl_E, 0, 0, 0);
 		v3_set(aux->int_curl_H, 0, 0, 0);
@@ -122,6 +154,10 @@ static void phy_sim_zero(phy_sim *sim)
 	size_t i;
 	phy_field_info *fi;
 
+	sim->pml_widths[0] = con_var_get_size("pml_x");
+	sim->pml_widths[1] = con_var_get_size("pml_y");
+	sim->pml_widths[2] = con_var_get_size("pml_z");
+
 	fi = &sim->field_info;
 
 	for (i = 0; i < 3; i++) {
@@ -149,6 +185,18 @@ int phy_sim_create(phy_sim *sim)
 	con_var_set("depth", "50");
 	con_var_set("scale", "7e+9"); 
 	con_var_set("time_delta", "0.1");
+	con_var_set("pml_x", "0");
+	con_var_set("pml_y", "0");
+	con_var_set("pml_z", "0");
+
+	con_var_set("src_z", "10"); 
+	con_var_set("src_angle", "0");
+
+	con_var_set("csg_x", NULL);
+	con_var_set("csg_y", NULL);
+	con_var_set("csg_z", NULL);
+	con_var_set("csg_r", NULL);
+	con_var_set("csg_v", NULL);
 
 	phy_sim_create_fields(sim);
 
@@ -157,16 +205,28 @@ int phy_sim_create(phy_sim *sim)
 
 void phy_sim_destroy(phy_sim *sim)
 {
-	size_t i;
-
+        phy_sim_csg_clear(sim);
+	phy_sim_destroy_fields(sim);
 	SDL_DestroyMutex(sim->rotate_lock);
- 
 	itc_chan_destroy(&sim->ctl);
+}
+
+static void phy_sim_configure(phy_sim *sim)
+{
+	phy_field_info *fi = &sim->field_info;
+
+	fi->spacing = con_var_get_float("scale") /
+	              min3(fi->dims[0], fi->dims[1], fi->dims[2]);
+	sim->time_delta = con_var_get_float("time_delta");
+
+	con_var_set("courant", va("%f", 3 * LIGHT_SPEED * sim->time_delta /
+	                          fi->spacing));
+
+	sim->source_xyplane.z = con_var_get_size("src_z");
+	sim->source_xyplane.Px = cos(DEG2RAD(con_var_get_float("src_angle")));
+	sim->source_xyplane.Py = sin(DEG2RAD(con_var_get_float("src_angle")));
+	sim->source_xyplane.enabled = true;
 
-	for (i = 0; i < 3; i++) {
-		free(sim->fields[i].E);
-		free(sim->fields[i].H);
-	}
 }
 
 int phy_sim_create_fields(phy_sim *sim)
@@ -179,29 +239,9 @@ int phy_sim_create_fields(phy_sim *sim)
 	fi->dims[0] = con_var_get_size("width");
 	fi->dims[1] = con_var_get_size("height");
 	fi->dims[2] = con_var_get_size("depth");
-	fi->spacing = con_var_get_float("scale") /
-	              min3(fi->dims[0], fi->dims[1], fi->dims[2]);
-	sim->time_delta = con_var_get_float("time_delta");
 
-	con_var_set("courant", va("%f", 3 * LIGHT_SPEED * sim->time_delta /
-	                          fi->spacing));
-
-	#if 0
-	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;
-	sim->source_xyplane.lambda = 20;
-	sim->source_xyplane.time = 2;
-	sim->source_xyplane.enabled = true;
+	printf("phy_sim_create_fields: %zux%zux%zu\n",
+	       fi->dims[0], fi->dims[1], fi->dims[2]);
 
 	fi->xstr = 3;
 	fi->ystr = fi->xstr * fi->dims[0];
@@ -228,6 +268,7 @@ int phy_sim_create_fields(phy_sim *sim)
 		goto error;
 	}
 
+	phy_sim_configure(sim);
 	phy_sim_zero(sim);
 
 	sim->running = false;
@@ -470,6 +511,44 @@ void phy_sim_step(phy_sim *sim) {
 	SDL_mutexV(sim->rotate_lock);
 }
 
+void phy_sim_csg_push(phy_sim *sim)
+{
+	char *shape_str;
+	phy_csg_shape shape;
+	phy_csg_step *step;
+
+	// TODO shape_str = con_var_get("csg_shape");
+
+	shape = CSG_SPHERE;
+
+	step = malloc(sizeof(phy_csg_step));
+	assert(step);
+
+	step->next = sim->csg_stack;
+	sim->csg_stack = step;
+
+	step->value = con_var_get_float("csg_v");
+
+	switch (shape) {
+	case CSG_SPHERE:
+		step->x[0] = con_var_get_float("csg_x");
+		step->x[1] = con_var_get_float("csg_y");
+		step->x[2] = con_var_get_float("csg_z");
+		step->d[0] = con_var_get_float("csg_r");
+		break;
+	}
+}
+
+void phy_sim_csg_clear(phy_sim *sim)
+{
+	phy_csg_step *step, *next;
+
+	for (step = sim->csg_stack; step; step = next) {
+		next = step->next;
+		free(step);
+	}
+}
+
 void phy_sim_debug(phy_sim *sim, size_t *coords)
 {
 	phy_field_info *fi = &sim->field_info;
@@ -538,10 +617,25 @@ void phy_run_command(phy_sim *sim, int number, void *data)
 
 	case PHY_CMD_RESET:
 		SDL_LockMutex(sim->rotate_lock);
-		phy_sim_destroy_fields(sim);
-		phy_sim_create_fields(sim);
+		if (sim->field_info.dims[0] == con_var_get_size("width") &&
+		    sim->field_info.dims[1] == con_var_get_size("height") &&
+		    sim->field_info.dims[2] == con_var_get_size("depth")) {
+			phy_sim_configure(sim);
+			phy_sim_zero(sim);
+		} else {
+			phy_sim_destroy_fields(sim);
+			phy_sim_create_fields(sim);
+		}
 		SDL_UnlockMutex(sim->rotate_lock);
 		break;
+
+	case PHY_CMD_CSG_PUSH:
+		phy_sim_csg_push(sim);
+		break;
+
+	case PHY_CMD_CSG_CLEAR:
+		phy_sim_csg_clear(sim);
+		break;
 	}
 }
 
diff --git a/src/physics.h b/src/physics.h
index 60db8bf..80f3fdf 100644
--- a/src/physics.h
+++ b/src/physics.h
@@ -38,8 +38,6 @@ typedef struct {
 
 	size_t z;
 	float Px, Py;
-
-	float lambda, time;
 } phy_source_xyplane;
 
 enum {
@@ -49,7 +47,22 @@ enum {
 	PHY_CMD_STEP,
 	PHY_CMD_ZERO,
 	PHY_CMD_DEBUG,
-	PHY_CMD_RESET
+	PHY_CMD_RESET,
+	PHY_CMD_CSG_PUSH,
+	PHY_CMD_CSG_CLEAR
+};
+
+typedef enum {
+	CSG_SPHERE
+} phy_csg_shape;
+
+typedef struct phy_csg_step_s phy_csg_step;
+struct phy_csg_step_s {
+	phy_csg_shape shape;
+	vec3_t x;
+	vec3_t d;
+	float value;
+	phy_csg_step *next;
 };
 
 typedef struct {
@@ -63,6 +76,9 @@ typedef struct {
 
 	phy_source_xyplane source_xyplane;
 
+	phy_csg_step *csg_stack;
+	float eps_min, eps_max;
+
 	// UI stuff
 	bool running;
 	itc_chan ctl;
@@ -79,6 +95,7 @@ int phy_sim_create_fields(phy_sim *sim);
 void phy_sim_destroy_fields(phy_sim *sim);
 void phy_sim_compute_const_fields(phy_sim *sim);
 void phy_sim_step(phy_sim *sim);
+void phy_sim_csg_clear(phy_sim *sim);
 
 int phy_thread(phy_sim *sim);
 
diff --git a/src/renderer.c b/src/renderer.c
index d68287c..6d159e7 100644
--- a/src/renderer.c
+++ b/src/renderer.c
@@ -373,7 +373,7 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
 {
 	phy_field_info *fi = &sim->field_info;
 	float *field;
-	size_t width, height, x, y, z;
+	size_t width, height, x[3], dx, dy, d0;
 	uint8_t *pixels;
 	int pitch;
 	float scale;
@@ -387,12 +387,21 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim,
 	if (flags & XSECTION_XY) {
 		width = fi->dims[0];
 		height = fi->dims[1];
+		d0 = 2;
+		dx = 0;
+		dy = 1;
 	} else if (flags & XSECTION_XZ) {
 		width = fi->dims[0];
 		height = fi->dims[2];
+		d0 = 1;
+		dx = 0;
+		dy = 2;
 	} else {
 		width = fi->dims[1];
 		height = fi->dims[2];
+		d0 = 0;
+		dx = 1;
+		dy = 2;
 	}
 
 	if (flags & XSECTION_E) {
@@ -419,61 +428,39 @@ 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->dims[2];       
-		if (z >= fi->dims[2])
-			z = fi->dims[2] - 1;
+	x[d0] = frac * fi->dims[d0];
+	if (x[d0] >= fi->dims[d0])
+		x[d0] = fi->dims[d0] - 1;
 
-		for (y = 0; y < fi->dims[1]; y++)
-		for (x = 0; x < fi->dims[0]; x++) {
+	if (flags & XSECTION_EPS)
+		for (x[dy] = 0; x[dy] < fi->dims[dy]; x[dy]++)
+		for (x[dx] = 0; x[dx] < fi->dims[dx]; x[dx]++) {
 			uint8_t *pixel;
-			float *point;
+			float gray;
+			phy_field_aux_point *aux;
 
-			pixel = pixels + (y * pitch + x * 4);
-			point = field + z * fi->zstr + y * fi->ystr +
-			        x * fi->xstr;
+			pixel = pixels + (x[dy] * pitch + x[dx] * 4);
+			aux = sim->aux + x[2] * fi->zstr1 + x[1] * fi->ystr1 +
+			      x[0] * fi->xstr1;
 
-			pixel[0] = r_float_to_u8(point[0] * scale);
-			pixel[1] = r_float_to_u8(point[1] * scale);
-			pixel[2] = r_float_to_u8(point[2] * scale);
+			gray = remap(aux->eps, sim->eps_min, sim->eps_max,
+			             0.0f, 1.0f);
+			pixel[0] = pixel[1] = pixel[2] = r_float_to_u8(gray);
 		}
-	} else if (flags & XSECTION_XZ) {
-		y = frac * fi->dims[1];
-		if (y >= fi->dims[1])
-			y = fi->dims[1] - 1;
-
-		for (z = 0; z < fi->dims[2]; z++)
-		for (x = 0; x < fi->dims[0]; x++) {
-			uint8_t *pixel;
-			float *point;
-
-			pixel = pixels + (z * pitch + x * 4);
-			point = field + z * fi->zstr + y * fi->ystr +
-			        x * fi->xstr;
-
-			pixel[0] = r_float_to_u8(point[0] * scale);
-			pixel[1] = r_float_to_u8(point[1] * scale);
-			pixel[2] = r_float_to_u8(point[2] * scale);
-		}
-	} else {
-		x = frac * fi->dims[0];
-		if (x >= fi->dims[0])
-			x = fi->dims[0] - 1;
-
-		for (z = 0; z < fi->dims[2]; z++)
-		for (y = 0; y < fi->dims[1]; y++) {
+	else
+		for (x[dy] = 0; x[dy] < fi->dims[dy]; x[dy]++)
+		for (x[dx] = 0; x[dx] < fi->dims[dx]; x[dx]++) {
 			uint8_t *pixel;
 			float *point;
 
-			pixel = pixels + (z * pitch + y * 4);
-			point = field + z * fi->zstr + y * fi->ystr +
-			        x * fi->xstr;
+			pixel = pixels + (x[dy] * pitch + x[dx] * 4);
+			point = field + x[2] * fi->zstr + x[1] * fi->ystr +
+				x[0] * fi->xstr;
 
 			pixel[0] = r_float_to_u8(point[0] * scale);
 			pixel[1] = r_float_to_u8(point[1] * scale);
 			pixel[2] = r_float_to_u8(point[2] * scale);
 		}
-	}
 
 	SDL_UnlockTexture(xsection->texture);
 
diff --git a/src/renderer.h b/src/renderer.h
index fadcd68..3ac10b8 100644
--- a/src/renderer.h
+++ b/src/renderer.h
@@ -47,6 +47,8 @@ void r_draw_text(r_window *rw, float x, float y, float h, char *text,
 #define XSECTION_PLANES  (XSECTION_XY|XSECTION_XZ|XSECTION_YZ)
 #define XSECTION_E       0x0008
 #define XSECTION_H       0x0010
+#define XSECTION_EPS     0x0020
+#define XSECTION_FIELDS  (XSECTION_E|XSECTION_H|XSECTION_EPS)
 
 
 typedef struct {
diff --git a/src/ui.c b/src/ui.c
index 5b58be7..868b6a1 100644
--- a/src/ui.c
+++ b/src/ui.c
@@ -338,16 +338,22 @@ void ui_event_window(SDL_Event *event, ui_window *uiw)
 
 		case SDLK_4:
 			ui_infof(uiw, "Przekrój: pola elektrycznego");
-			sv->xsection_flags &= ~XSECTION_H;
+			sv->xsection_flags &= ~XSECTION_FIELDS;
 			sv->xsection_flags |= XSECTION_E;
 			break;
 
 		case SDLK_5:
 			ui_infof(uiw, "Przekrój: pola magnetycznego");
-			sv->xsection_flags &= ~XSECTION_E;
+			sv->xsection_flags &= ~XSECTION_FIELDS;
 			sv->xsection_flags |= XSECTION_H;
 			break;
 
+		case SDLK_6:
+			ui_infof(uiw, "Przekrój: pola przenikalności");
+			sv->xsection_flags &= ~XSECTION_FIELDS;
+			sv->xsection_flags |= XSECTION_EPS;
+			break;
+
 		case SDLK_z:
 			if (sv->select_valid)
 			{