commit 6f0c04f4efc298cb78ffccd492849d61e6cd8e4c
parent 381594965e7475474acb11b1d66615853228bbfa
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 6 Jul 2022 16:26:51 +0200
Refactor the realisation of the sensitivity computation
Diffstat:
1 file changed, 90 insertions(+), 38 deletions(-)
diff --git a/src/sgs_compute_sensitivity_translation.c b/src/sgs_compute_sensitivity_translation.c
@@ -28,28 +28,82 @@
/* Sugar syntax */
enum {X, Y, Z};
-enum {WEIGHT, SENSIT, WEIGHTS_COUNT__};
+enum {WEIGHT, SENSIB, WEIGHTS_COUNT__};
/* FIXME this should be variables */
#if 0
+/* Il semblerait que le cube soit de taille {4, 4, 3} */
static const double RECV_MIN[2] = {0.5, 1.5};
static const double RECV_MAX[2] = {1.5, 2.5};
static const double EMIT_E_THRESHOLD = 2;
static const double EMIT_E_SZ[3] = {0, 4, 2};
static const double EMIT_S_SZ[3] = {4, 4, 0};
#else
-static const double RECV_MIN[2] = {0.125, 0.375};
-static const double RECV_MAX[2] = {0.375, 0.625};
+static const double RECV_MIN[2] = {0.125/*<=>1/8*/, 0.375/*<=>3/8*/};
+static const double RECV_MAX[2] = {0.375/*<=>3/8*/, 0.625/*<=>5/8*/};
static const double EMIT_E_THRESHOLD = 2.0/3.0;
static const double EMIT_E_SZ[3] = {0, 1, 2.0/3.0};
static const double EMIT_S_SZ[3] = {1, 1, 0};
#endif
+/* Vecteur de déformation (à renommer en gamma ?) */
static const double V[3] = {0, 0, 1};
/*******************************************************************************
* Helper functions
******************************************************************************/
+static double
+compute_rho
+ (const double position[2], /* Position sur la surface réfléchissante */
+ const double size[2]) /* Taille de la surface réfléchissante */
+{
+ ASSERT(position && size);
+ return 0.25
+ * (1 - cos(2*PI*position[X]/size[X]))
+ * (1 - cos(2*PI*position[Y]/size[Y]));
+}
+
+static void
+compute_grad_rho
+ (const double position[2], /* Position sur la surface réfléchissante */
+ const double size[2], /* Taille de la surface réfléchissante */
+ double out_gradient[2])
+{
+ ASSERT(position && size && out_gradient);
+ out_gradient[X] = 0.25
+ * (((2*PI)/size[X])*sin(2*PI*position[X]/size[X]))
+ * (1 - cos(2*PI*position[Y]/size[Y]));
+ out_gradient[Y] = 0.25
+ * (((2*PI)/size[Y])*sin(2*PI*position[Y]/size[Y]))
+ * (1 - cos(2*PI*position[X]/size[X]));
+}
+
+static double
+compute_I
+ (const double position[2], /* Position sur la surface émettrice */
+ const double size[2]) /* Taille de la surface émettrice */
+{
+ ASSERT(position && size);
+ return
+ (1 - cos(2*PI*position[X]/size[X]))
+ * (1 - cos(2*PI*position[Y]/size[Y]));
+}
+
+static void
+compute_grad_I
+ (const double position[2], /* Position sur la surface réfléchissante */
+ const double size[2], /* Taille de la surface réfléchissante */
+ double out_gradient[2])
+{
+ ASSERT(position && size && out_gradient);
+ out_gradient[X] =
+ (((2*PI)/size[X])*sin(2*PI*position[X]/size[X]))
+ * (1 - cos(2*PI*position[Y]/size[Y]));
+ out_gradient[Y] =
+ (((2*PI)/size[Y])*sin(2*PI*position[Y]/size[Y]))
+ * (1 - cos(2*PI*position[X]/size[X]));
+}
+
static res_T
realisation
(void* weights,
@@ -65,8 +119,8 @@ realisation
double recv_min[2];
double recv_max[2];
double emit_e_threshold;
- double emit_e_sz[3];
- double emit_s_sz[3];
+ double emit_e_sz[3], emit_e_sz_2d[2];
+ double emit_s_sz[3], emit_s_sz_2d[2];
double box_low[3];
double box_upp[3];
@@ -78,8 +132,8 @@ realisation
double dir_emit_s[3];
double dir_spec_e[3];
double dir_spec_s[3];
- double pos_emit_e[3];
- double pos_emit_s[3];
+ double pos_emit_e[3], pos_emit_e_2d[2];
+ double pos_emit_s[3], pos_emit_s_2d[2];
double pos_recv[3];
double range[2];
@@ -97,11 +151,11 @@ realisation
double beta_spec_e;
double beta_emit_s;
double beta_spec_s;
-
+
double rho;
- double grad_rho[3];
+ double grad_rho[3], grad_rho_2d[2];
double I;
- double grad_I[3];
+ double grad_I[3], grad_I_2d[2];
double S;
@@ -143,12 +197,12 @@ realisation
/* Sample the cosine weighted sampling of the emissive direction */
ssp_ran_hemisphere_cos(rng, frag.normal, dir_emit_s, NULL);
-
+
/* Trace the ray from the sensitivity emissive surface */
sgs_geometry_trace_ray
(sgs->geom, pos_emit_s, dir_emit_s, range, surf_emit_s, &hit);
if(SGS_HIT_NONE(&hit)) {
- res = RES_OK;
+ res = RES_BAD_ARG;
goto error;
}
pos_recv[X] = pos_emit_s[X] + hit.distance*dir_emit_s[X];
@@ -167,7 +221,7 @@ realisation
sgs_geometry_trace_ray
(sgs->geom, pos_emit_s, dir_spec_s, range, surf_emit_s, &hit);
if(SGS_HIT_NONE(&hit)) {
- res = RES_OK;
+ res = RES_BAD_ARG;
goto error;
}
pos_emit_e[X] = pos_emit_s[X] + hit.distance*dir_spec_s[X];
@@ -209,30 +263,28 @@ realisation
beta_emit_e = d3_normalize(u_emit_e, u_emit_e);
beta_spec_e = d3_normalize(u_spec_e, u_spec_e);
- rho = 0.25
- * (1 - cos(2*PI*pos_emit_s[X]/emit_s_sz[X]))
- * (1 - cos(2*PI*pos_emit_s[Y]/emit_s_sz[Y]));
- grad_rho[X] = 0.25
- * (((2*PI)/emit_s_sz[X])*sin(2*PI*pos_emit_s[X]/emit_s_sz[X]))
- * (1 - cos(2*PI*pos_emit_s[Y]/emit_s_sz[Y]));
- grad_rho[Y] = 0.25
- * (((2*PI)/emit_s_sz[Y])*sin(2*PI*pos_emit_s[Y]/emit_s_sz[Y]))
- * (1 - cos(2*PI*pos_emit_s[X]/emit_s_sz[X]));
+ pos_emit_s_2d[X] = pos_emit_s[X];
+ pos_emit_s_2d[Y] = pos_emit_s[Y];
+ emit_s_sz_2d[X] = emit_s_sz[X];
+ emit_s_sz_2d[Y] = emit_s_sz[Y];
+ rho = compute_rho(pos_emit_s_2d, emit_s_sz_2d);
+ compute_grad_rho(pos_emit_s_2d, emit_s_sz_2d, grad_rho_2d);
+ grad_rho[X] = grad_rho_2d[X];
+ grad_rho[Y] = grad_rho_2d[Y];
grad_rho[Z] = 0;
- I =
- (1 - cos(2*PI*pos_emit_e[Y]/emit_e_sz[Y]))
- * (1 - cos(2*PI*pos_emit_e[Z]/emit_e_sz[Z]));
+ pos_emit_e_2d[X] = pos_emit_e[Y];
+ pos_emit_e_2d[Y] = pos_emit_e[Z];
+ emit_e_sz_2d[X] = emit_e_sz[Y];
+ emit_e_sz_2d[Y] = emit_e_sz[Z];
+ I = compute_I(pos_emit_e_2d, emit_e_sz_2d);
+ compute_grad_I(pos_emit_e_2d, emit_e_sz_2d, grad_I_2d);
grad_I[X] = 0;
- grad_I[Y] =
- (((2*PI)/emit_e_sz[Y])*sin(2*PI*pos_emit_e[Y]/emit_e_sz[Y]))
- * (1 - cos(2*PI*pos_emit_e[Z]/emit_e_sz[Z]));
- grad_I[Z] =
- (((2*PI)/emit_e_sz[Z])*sin(2*PI*pos_emit_e[Z]/emit_e_sz[Z]))
- * (1 - cos(2*PI*pos_emit_e[Y]/emit_e_sz[Y]));
-
- S =
- - beta_emit_s * d3_dot(grad_rho, u_emit_s) * I
+ grad_I[Y] = grad_I_2d[X];
+ grad_I[Z] = grad_I_2d[Y];
+
+ S =
+ - beta_emit_s * d3_dot(grad_rho, u_emit_s) * I
- rho * beta_emit_s * beta_emit_e * d3_dot(grad_I, u_emit_e)
+ rho * beta_spec_s * beta_spec_e * d3_dot(grad_I, u_spec_e);
@@ -241,7 +293,7 @@ realisation
exit:
SMC_DOUBLEN(weights)[WEIGHT] = w;
- SMC_DOUBLEN(weights)[SENSIT] = s;
+ SMC_DOUBLEN(weights)[SENSIB] = s;
return res;
error:
goto exit;
@@ -314,11 +366,11 @@ compute_sensitivity_translation(struct sgs* sgs)
}
/* Print the result */
- sgs_log(sgs, "estim ~ %g +/- %g; sensib ~ %g +/- %g\n",
+ sgs_log(sgs, "estim ~ %g +/- %g; sensib ~ %g +/- %g\n",
SMC_DOUBLEN(status.E)[WEIGHT],
SMC_DOUBLEN(status.SE)[WEIGHT],
- SMC_DOUBLEN(status.E)[SENSIT],
- SMC_DOUBLEN(status.SE)[SENSIT]);
+ SMC_DOUBLEN(status.E)[SENSIB],
+ SMC_DOUBLEN(status.SE)[SENSIB]);
sgs_log(sgs, "#failures = %lu/%lu\n",
(unsigned long)status.NF,
(unsigned long)sgs->nrealisations);
