Finish CSG, permittivity visualisation.

1 files changed, 142 insertions, 48 deletions

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;
 	}
 }