commit a3ac5a3a7fe8378facbb497a53d9918889ddb0ba
parent 60c23cce2e3d1a847cd65fec2f08e1ec3fbf9e7a
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 26 Feb 2016 16:34:21 +0100
Finalise the computation of the phase function
Compute the [cumulative] phase function from the computed [cumulative]
differential cross section and the scattering cross section. Note that
the inverse cumulative phase function is still not computed.
Furthermore, this computation is still disabled since the phase function
model for wide angles are still not correctly solved.
Diffstat:
2 files changed, 233 insertions(+), 35 deletions(-)
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -79,6 +79,14 @@ struct mc_accum {
struct mc_data* differential_cross_section_cumulative;
};
+/* Estimated phase function */
+struct phase_function {
+ /* Per angle values */
+ double* values;
+ double* cumulative;
+ double* inverse_cumulative;
+};
+
/* The geometry context stores data that are constants for a given geometry
* distribution. It stores limit scattering angles as well as geometry material
* since the particle of a distribution shared the same optical properties.
@@ -115,7 +123,8 @@ struct integrator_context {
struct sschiff_estimator {
double* wavelengths; /* List of wavelengths in vacuum */
double* angles; /* List of scattering angles */
- struct mc_accum* accums; /* Monte Carlo estimation */
+ struct phase_function* phase_functions; /* Per wavelength phase function */
+ struct mc_accum* accums; /* Per wavelength Monte Carlo estimation */
struct sschiff_device* dev;
size_t nrealisations;
ref_T ref;
@@ -192,14 +201,20 @@ build_transform
transform[11] = -f3_dot(basis+6, pos);
}
-static void
+static FINLINE double
+expected_value(const struct mc_data* data, const size_t nrealisations)
+{
+ return data->weight / (double)nrealisations;
+}
+
+static INLINE void
get_mc_value
(const struct mc_data* data,
const size_t nrealisations,
struct sschiff_state* value)
{
ASSERT(data && value);
- value->E = data->weight / (double)nrealisations;
+ value->E = expected_value(data, nrealisations);
value->V = data->sqr_weight / (double)nrealisations
- (data->weight * data->weight) / (double)(nrealisations * nrealisations);
value->SE = sqrt(value->V / (double)nrealisations);
@@ -326,6 +341,7 @@ compute_path_length
double len = 0;
float dst;
ASSERT(scn && first_hit && ray_org && ray_dir && !S3D_HIT_NONE(first_hit));
+ (void)dev;
dst = range[0] = first_hit->distance;
hit_from = *first_hit;
@@ -516,26 +532,35 @@ function_fdf(const gsl_vector* vec, void* params, gsl_vector* res, gsl_matrix* J
return GSL_SUCCESS;
}
+/* Solve the parameters used to connect the Monte-Carlo estimation of the
+ * [cumulative] differential cross sections for small angles with their
+ * analytical estimation for wide angles */
static res_T
-phase_function_post_process
+solve_wide_angles_phase_function_model_parameters
(struct sschiff_device* dev,
- const size_t iwavelength,
- const struct geometry_context* ctx,
- struct sschiff_estimator* estimator)
+ const double angle,
+ const double sin_half_angle,
+ const double cos_angle,
+ const double cos_2_angle,
+ const double scattering_cross_section,
+ const double differential_cross_section,
+ const double differential_cross_section_cumulative,
+ double* n,
+ double* r)
{
gsl_multiroot_fdfsolver* solver = NULL;
gsl_multiroot_function_fdf func;
gsl_vector* init_val = NULL;
gsl_vector* root = NULL;
- struct mc_accum* accums;
struct function_arg arg;
- double n, r; /* Estimated value */
- size_t ilimit_angle; /* Index of the limit angle of the current wavelength */
int status;
int iter;
const int max_iter = 1000; /* Maximum number of Newton iterations */
res_T res = RES_OK;
- ASSERT(dev);
+ ASSERT(dev && n && r);
+ ASSERT(eq_eps(sin(angle*0.5), sin_half_angle, 1.e-6));
+ ASSERT(eq_eps(cos(angle), cos_angle, 1.e-6));
+ ASSERT(eq_eps(cos(2.0*angle), cos_2_angle, 1.e-6));
/* Create a Newton-Raphson nonlinear solver without derivatives */
solver = gsl_multiroot_fdfsolver_alloc(gsl_multiroot_fdfsolver_newton, 1);
@@ -555,27 +580,17 @@ phase_function_post_process
goto error;
}
- /* Fetch attribute */
- accums = estimator->accums + iwavelength;
- ilimit_angle = ctx->limit_angles[iwavelength]-1;
-
/* Fill system arguments */
- arg.angle = estimator->angles[ilimit_angle];
- arg.scattering_cross_section =
- accums->scattering_cross_section.weight
- / (double)estimator->nrealisations;
- arg.differential_cross_section =
- accums->differential_cross_section[ilimit_angle].weight
- / (double)estimator->nrealisations;
- arg.differential_cross_section_cumulative =
- accums->differential_cross_section_cumulative[ilimit_angle].weight
- / (double)estimator->nrealisations;
+ arg.angle = angle;
+ arg.scattering_cross_section = scattering_cross_section;
+ arg.differential_cross_section = differential_cross_section;
+ arg.differential_cross_section_cumulative = differential_cross_section_cumulative;
/* Precompute some values to speed up Newton iterations */
- arg.cos_2_angle = cos(2.0 * arg.angle);
- arg.cos_angle = cos(arg.angle);
- arg.cos_half_angle = cos(arg.angle * 0.5);
- arg.sin_half_angle = sin(arg.angle * 0.5);
+ arg.cos_2_angle = cos_2_angle;
+ arg.cos_angle = cos_angle;
+ arg.cos_half_angle = cos(angle * 0.5);
+ arg.sin_half_angle = sin_half_angle;
arg.sqr_cos_angle = arg.cos_angle * arg.cos_angle;
arg.sqr_sin_half_angle = arg.sin_half_angle * arg.sin_half_angle;
@@ -595,6 +610,7 @@ phase_function_post_process
/* Launch the Newton estimation. Iterate up to `max_iter'or until the
* estimation is "good enough". */
iter = 0;
+
do {
status = gsl_multiroot_fdfsolver_iterate(solver);
if(status == GSL_SUCCESS) {
@@ -612,11 +628,11 @@ phase_function_post_process
/* Retrieve the estimated "n" value */
root = gsl_multiroot_fdfsolver_root(solver);
- n = gsl_vector_get(root, 0);
- printf("%g\n", n);
+ *n = gsl_vector_get(root, 0);
+ printf("theta_limit = %g, n = %g\n", arg.angle, *n); /* TODO remove this */
/* Compute r from the estimated n */
- r = 2 * arg.differential_cross_section * pow(arg.sin_half_angle, n)
+ *r = 2 * arg.differential_cross_section * pow(arg.sin_half_angle, *n)
/ (1 + arg.sqr_cos_angle);
exit:
@@ -627,6 +643,144 @@ error:
goto exit;
}
+/* Helper function used to compute one term of the analytic evaluation of the
+ * cumulative differential cross section */
+static FINLINE double
+compute_differential_cross_section_cumulative_term
+ (const double sin_half_angle,
+ const double cos_angle,
+ const double cos_2_angle,
+ const double n)
+{
+ return
+ -pow(sin_half_angle, 2-n)
+ * ((n*n - 6*n + 8)*cos_2_angle + 3*n*n - 8*(n-2)*cos_angle - 26*n + 72)
+ / ((n-2)*(n-4)*(n-6));
+}
+
+static res_T
+compute_phase_function
+ (struct sschiff_device* dev,
+ const size_t iwlen,
+ const size_t nangles,
+ const struct geometry_context* ctx,
+ struct sschiff_estimator* estimator)
+{
+ double limit_angle;
+ double sin_half_limit_angle;
+ double cos_limit_angle;
+ double cos_2_limit_angle;
+ double coef_limit;
+ double n, r; /* Connector values */
+ double scattering_cross_section;
+ double rcp_scattering_cross_section;
+ double limit_differential_cross_section;
+ double limit_differential_cross_section_cumulative;
+ size_t ilimit_angle; /* Index of the limit angle of the current wavelength */
+ size_t iangle;
+ res_T res = RES_OK;
+ ASSERT(dev && ctx && estimator);
+
+ /* Fetch the limit angle and precompute some values */
+ ilimit_angle = ctx->limit_angles[iwlen]-1;
+ limit_angle = estimator->angles[ilimit_angle];
+ sin_half_limit_angle = sin(0.5*limit_angle);
+ cos_limit_angle = cos(limit_angle);
+ cos_2_limit_angle = cos(2.0*limit_angle);
+
+ /* Find the connector values between small and whide angles */
+ scattering_cross_section = expected_value
+ (&estimator->accums[iwlen].scattering_cross_section,
+ estimator->nrealisations);
+ limit_differential_cross_section = expected_value
+ (&estimator->accums[iwlen].differential_cross_section[ilimit_angle],
+ estimator->nrealisations);
+ limit_differential_cross_section_cumulative = expected_value
+ (&estimator->accums[iwlen].differential_cross_section_cumulative[ilimit_angle],
+ estimator->nrealisations);
+ res = solve_wide_angles_phase_function_model_parameters
+ (dev,
+ limit_angle,
+ sin_half_limit_angle,
+ cos_limit_angle,
+ cos_2_limit_angle,
+ scattering_cross_section,
+ limit_differential_cross_section,
+ limit_differential_cross_section_cumulative,
+ &n, &r);
+ if(res != RES_OK) goto error;
+
+ /* Precomputed vlaue */
+ coef_limit = compute_differential_cross_section_cumulative_term
+ (sin_half_limit_angle, cos_limit_angle, cos_2_limit_angle, n);
+
+ /* Fill the [cumulative] differential cross sections of wide angles with the
+ * analytic model */
+ FOR_EACH(iangle, ilimit_angle+1, nangles) {
+ struct mc_accum* accum = estimator->accums + iwlen;
+ double angle, sin_half_angle, cos_angle, cos_2_angle;
+ double coef;
+
+ /* Precompute some values */
+ angle = estimator->angles[iangle];
+ cos_angle = cos(angle);
+ cos_2_angle = cos(2.0*angle);
+ sin_half_angle = sin(angle*0.5);
+
+ /* Analytically compute the differential cross section */
+ accum->differential_cross_section[iangle].sqr_weight = -1.f;
+ accum->differential_cross_section[iangle].weight =
+ r * (1+cos_angle*cos_angle) / (2.0*pow(sin_half_angle, n));
+
+ /* Analytically compute the differential cross section cumulative */
+ coef = compute_differential_cross_section_cumulative_term
+ (sin_half_angle, cos_angle, cos_2_angle, n);
+ accum->differential_cross_section_cumulative[iangle].sqr_weight = -1.f;
+ accum->differential_cross_section_cumulative[iangle].weight =
+ limit_differential_cross_section_cumulative + 2*PI*r*(coef-coef_limit);
+ }
+
+ /* Check the post condition of the cumulative differential cross section */
+ /*ASSERT(eq_eps
+ (estimator->accums[iwlen].differential_cross_section_cumulative[nangles-1].weight,
+ scattering_cross_section, 1.e-3));*/
+
+ /* Compute the [cumulative] phase function from the [cumulative] differential
+ * cross section */
+ rcp_scattering_cross_section = 1.0 / scattering_cross_section;
+ FOR_EACH(iangle, 0, ilimit_angle) {
+ estimator->phase_functions[iwlen].values[iangle] = expected_value
+ (&estimator->accums[iwlen].differential_cross_section[iangle],
+ estimator->nrealisations);
+ estimator->phase_functions[iwlen].values[iangle] *=
+ rcp_scattering_cross_section;
+
+ estimator->phase_functions[iwlen].cumulative[iangle] = expected_value
+ (&estimator->accums[iwlen].differential_cross_section_cumulative[iangle],
+ estimator->nrealisations);
+ estimator->phase_functions[iwlen].cumulative[iangle] *=
+ rcp_scattering_cross_section;
+ }
+ FOR_EACH(iangle, ilimit_angle + 1, nangles) {
+ estimator->phase_functions[iwlen].values[iangle] =
+ estimator->accums[iwlen].differential_cross_section[iangle].weight
+ * rcp_scattering_cross_section;
+
+ estimator->phase_functions[iwlen].values[iangle] =
+ estimator->accums[iwlen].differential_cross_section_cumulative[iangle].weight
+ * rcp_scattering_cross_section;
+ }
+ estimator->phase_functions[iwlen].values[ilimit_angle] =
+ limit_differential_cross_section * rcp_scattering_cross_section;
+ estimator->phase_functions[iwlen].cumulative[ilimit_angle]=
+ limit_differential_cross_section_cumulative * rcp_scattering_cross_section;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
static res_T
accum_monte_carlo_weight
(struct sschiff_device* dev,
@@ -1080,6 +1234,19 @@ estimator_release(ref_T* ref)
}
sa_release(estimator->accums);
}
+ if(estimator->phase_functions) {
+ const size_t nwavelengths = sa_size(estimator->phase_functions);
+ size_t i;
+ FOR_EACH(i, 0, nwavelengths) {
+ if(estimator->phase_functions[i].values)
+ sa_release(estimator->phase_functions[i].values);
+ if(estimator->phase_functions[i].cumulative)
+ sa_release(estimator->phase_functions[i].cumulative);
+ if(estimator->phase_functions[i].inverse_cumulative)
+ sa_release(estimator->phase_functions[i].inverse_cumulative);
+ }
+ sa_release(estimator->phase_functions);
+ }
MEM_RM(dev->allocator, estimator);
SSCHIFF(device_ref_put(dev));
}
@@ -1153,7 +1320,6 @@ estimator_create
}
memset(estimator->accums, 0, sizeof(estimator->accums[0])*nwavelengths);
FOR_EACH(i, 0, nwavelengths) {
- /* FIXME rename the following estimation */
if(!sa_add(estimator->accums[i].differential_cross_section, nangles)) {
log_error(dev,
"Couldn't allocate the differential cross sections to estimate.\n");
@@ -1167,6 +1333,37 @@ estimator_create
goto error;
}
}
+
+ /* Allocate the phase function data */
+ if(!sa_add(estimator->phase_functions, nwavelengths)) {
+ log_error(dev,
+ "Couldn't allocate the per wavelength phase functions data.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ memset(estimator->phase_functions, 0,
+ sizeof(estimator->phase_functions[0])*nwavelengths);
+ FOR_EACH(i, 0, nwavelengths) {
+ if(!sa_add(estimator->phase_functions[i].values, nangles)) {
+ log_error(dev,
+ "Couldn't allocate the per angle phase function values.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ if(!sa_add(estimator->phase_functions[i].cumulative, nangles)) {
+ log_error(dev,
+ "Couldn't allocate the per angle cumulative phase function.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ if(!sa_add(estimator->phase_functions[i].inverse_cumulative, nangles)) {
+ log_error(dev,
+ "Couldn't allocate the per angle inverse cumulative phase function.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ }
+
exit:
*out_estimator = estimator;
return res;
@@ -1307,9 +1504,10 @@ sschiff_integrate
}
}
+ /* TODO use parallelism */
#if 0
FOR_EACH(i, 0, nwavelengths) {
- res = phase_function_post_process(dev, i, &geom_ctx, estimator);
+ res = compute_phase_function(dev, i, nangles, &geom_ctx, estimator);
if(res != RES_OK) goto error;
}
#endif
diff --git a/src/test_sschiff_estimator_sphere.c b/src/test_sschiff_estimator_sphere.c
@@ -252,7 +252,6 @@ check_schiff_estimation
distrib.material.get_property = get_material_property;
distrib.material.material = sampler_ctx;
- distrib.caracteristic_length = 1.0;
distrib.sample = sample_sphere;
distrib.context = sampler_ctx;
@@ -267,6 +266,7 @@ check_schiff_estimation
double dst;
sampler_ctx->mean_radius = results[i].mean_radius;
+ distrib.caracteristic_length = sampler_ctx->mean_radius * 2.0;
time_current(&t0);
CHECK(sschiff_integrate(dev, rng, &distrib, &wavelength, 1,
