commit ae3f5f8178fb7c601a061671d03a1a7ada6f092f
parent 3ff881a2accdc907ad9b5d7de412e8ab88d366a0
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 2 Mar 2016 15:24:32 +0100
Conditionally discard phase function evaluation
Do not estimate the differential cross section [cumulative] if no valid
limit scattering angle can be found.
Diffstat:
2 files changed, 96 insertions(+), 42 deletions(-)
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -50,6 +50,7 @@
#include <gsl/gsl_multiroots.h>
#define MAX_FAILURES 10
+#define INVALID_LIMIT_ANGLE SIZE_MAX
/* Accumulator of double precision values */
struct accum {
@@ -120,6 +121,7 @@ struct sschiff_estimator {
struct mc_accum* accums; /* Monte Carlo estimation */
size_t nrealisations; /* # realisation used by MC estimations */
+ int no_phase_function; /* Set to 1 when no phase function is computed */
struct sschiff_device* dev;
ref_T ref;
@@ -366,6 +368,10 @@ accum_differential_cross_section
double tmp;
size_t iangle;
+ /* Avoid evaluating the differential cross sections for wavelengths with no
+ * valid limit scattering angle */
+ if(ctx->estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE) continue;
+
/* Fetch optical properties */
k_e = ctx->estimator->wavenumbers[iwlen];
n_r = ctx->estimator->properties[iwlen].relative_real_refractive_index;
@@ -400,7 +406,7 @@ accum_differential_cross_section
ctx->accums[iwlen].differential_cross_section[iangle] += weight;
bessel_j1 = gsl_sf_bessel_J1(k_e_angle_delta_r);
- weight = 2*PI * ctx->estimator->angles[iangle] * bessel_j1 / k_e_delta_r * tmp;
+ weight = 2*PI*ctx->estimator->angles[iangle]*bessel_j1/k_e_delta_r*tmp;
ctx->accums[iwlen].differential_cross_section_cumulative[iangle] += weight;
}
}
@@ -464,6 +470,9 @@ accum_monte_carlo_weight
/* Compute and accumulate the cross section weight */
accum_cross_section(ctx, plane_area, path_length[0]);
+ /* Avoid the estimation of the phase function */
+ if(ctx->estimator->no_phase_function) break;
+
/* Uniformly sample a position onto the projection plane and use it as the
* origin of the 2nd ray to trace */
sample[0] = ssp_rng_uniform_double(ctx->rng, -1.0, 1.0);
@@ -542,6 +551,9 @@ radiative_properties
ctx->accums[iwlen].scattering_cross_section = 0;
ctx->accums[iwlen].avg_projected_area = 0;
+ /* Do not clean up accumulators of invalid differential cross section */
+ if(ctx->estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE) continue;
+
/* Clean up to "limit_angle" since great angles are analytically computed */
memset(ctx->accums[iwlen].differential_cross_section, 0,
sizeof(double)*ctx->estimator->limit_angles[iwlen]);
@@ -592,6 +604,9 @@ radiative_properties
MC_ACCUM(scattering_cross_section);
MC_ACCUM(avg_projected_area);
+ if(ctx->estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE)
+ continue;
+
/* Accum up to "limit angle" since great angles are analitically computed */
FOR_EACH(iangle, 0, ctx->estimator->limit_angles[iwlen]) {
MC_ACCUM(differential_cross_section[iangle]);
@@ -608,10 +623,14 @@ error:
ctx->mc_accums[iwlen].absorption_cross_section = MC_DATA_NULL;
ctx->mc_accums[iwlen].scattering_cross_section = MC_DATA_NULL;
ctx->mc_accums[iwlen].avg_projected_area = MC_DATA_NULL;
+
+ if(ctx->estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE)
+ continue;
+
memset(ctx->mc_accums[iwlen].differential_cross_section, 0,
- sizeof(double)*ctx->nangles);
+ sizeof(double)*ctx->estimator->limit_angles[iwlen]);
memset(ctx->mc_accums[iwlen].differential_cross_section_cumulative, 0,
- sizeof(double)*ctx->nangles);
+ sizeof(double)*ctx->estimator->limit_angles[iwlen]);
}
goto exit;
}
@@ -872,7 +891,7 @@ compute_phase_function
if(n < 2 || n > 4) {
log_warning(estimator->dev,
"The wide scattering angles phase function model parameter `n' is not in\n"
-"[0, 4]. One might increase the number of realisations.\n"
+"[2, 4]. One might increase the number of realisations.\n"
" n = %g\n"
" scattering cross section = %g +/- %g\n"
" differential cross section = %g +/- %g\n"
@@ -919,7 +938,7 @@ compute_phase_function
/* 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));*/
+ scattering_cross_section.E, 1.e-3));*/
/* Compute the [cumulative] phase function for small angles */
rcp_scattering_cross_section = 1.0/scattering_cross_section.E;
@@ -961,6 +980,8 @@ compute_phase_function
exit:
return res;
error:
+ memset(estimator->phase_functions[iwlen].values, 0, nangles);
+ memset(estimator->phase_functions[iwlen].cumulative, 0, nangles);
goto exit;
}
@@ -1087,11 +1108,16 @@ integrator_context_init
} (void)0
RESIZE(ctx->accums, nwavelengths);
FOR_EACH(iwlen, 0, nwavelengths) {
+ /* Do not allocate the accumulators if no limit angle was found */
+ if(estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE) continue;
RESIZE(ctx->accums[iwlen].differential_cross_section, nangles);
RESIZE(ctx->accums[iwlen].differential_cross_section_cumulative, nangles);
}
+
RESIZE(ctx->mc_accums, nwavelengths);
FOR_EACH(iwlen, 0, nwavelengths) {
+ /* Do not allocate the accumulators if no limit angle was found */
+ if(estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE) continue;
RESIZE(ctx->mc_accums[iwlen].differential_cross_section, nangles);
RESIZE(ctx->mc_accums[iwlen].differential_cross_section_cumulative, nangles);
}
@@ -1202,6 +1228,7 @@ estimator_create
struct sschiff_estimator* estimator = NULL;
double rcp_PI_Lc;
size_t i;
+ int hold_valid_limit_angle = 0;
res_T res = RES_OK;
ASSERT(dev && distrib && wavelengths && nwavelengths);
ASSERT(angles_distrib && nangles && out_estimator);
@@ -1288,54 +1315,69 @@ estimator_create
/* Search for limit scattering angle */
theta_l = sqrt(lambda_e * rcp_PI_Lc);
if(theta_l <= 0 || theta_l >= PI) {
- log_error(dev,
+ log_warning(dev,
"Invalid theta limit, i.e. angle between small and wide scattering angles\n"
"`%g'. The phase function for the wavelengths `%g' and its [inverse]\n"
"cumulative will be not computed.\n",
theta_l, estimator->wavelengths[i]);
- res = RES_BAD_ARG;
- goto error;
- }
- angle = search_lower_bound(&theta_l, estimator->angles, nangles,
- sizeof(double), cmp_double);
- if(!angle) {
- log_error(dev,
+ estimator->limit_angles[i] = INVALID_LIMIT_ANGLE;
+ } else {
+ angle = search_lower_bound(&theta_l, estimator->angles, nangles,
+ sizeof(double), cmp_double);
+ if(!angle) {
+ log_error(dev,
"Can't find a limit scattering angle for theta limit `%g'.\n", theta_l);
- res = RES_BAD_ARG;
- goto error;
- }
- if(eq_eps(*angle, PI, 1.e-6)) {
- log_error(dev, "Invalid limit scattering angle `%g'.\n", *angle);
- res = RES_BAD_ARG;
- goto error;
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ if(eq_eps(*angle, PI, 1.e-6)) {
+ log_error(dev, "Invalid limit scattering angle `%g'.\n", *angle);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ /* Define the index of the first "wide scattering angle" */
+ estimator->limit_angles[i] = (size_t)(angle - estimator->angles);
+ ASSERT(estimator->limit_angles[i] < nangles);
+ hold_valid_limit_angle = 1;
}
- /* Define the index of the first "wide scattering angle" */
- estimator->limit_angles[i] = (size_t)(angle - estimator->angles);
- ASSERT(estimator->limit_angles[i] < nangles);
+ }
+ estimator->no_phase_function = !hold_valid_limit_angle;
+ if(estimator->no_phase_function) {
+ /* No phase function => The scattering angles are useless */
+ sa_release(estimator->angles);
+ estimator->angles = NULL;
}
/* Allocate the estimator accumulators */
RESIZE(estimator->accums, nwavelengths,
"Couldn't allocate the Monte Carlo accumulator of the estimator.\n");
- memset(estimator->accums, 0, sizeof(estimator->accums[0])*nwavelengths);
- FOR_EACH(i, 0, nwavelengths) {
- RESIZE(estimator->accums[i].differential_cross_section, nangles,
- "Couldn't allocate the differential cross sections to estimate.\n");
- RESIZE(estimator->accums[i].differential_cross_section_cumulative, nangles,
- "Couldn't allocate the cumulative differential cross sections to estimate.\n");
+ if(!estimator->no_phase_function) {
+ FOR_EACH(i, 0, nwavelengths) {
+ /* Do not allocate the MC accumulator if no limit angle was found */
+ if(estimator->limit_angles[i] == INVALID_LIMIT_ANGLE) continue;
+ RESIZE(estimator->accums[i].differential_cross_section, nangles,
+ "Couldn't allocate the differential cross sections to estimate.\n");
+ RESIZE(estimator->accums[i].differential_cross_section_cumulative, nangles,
+ "Couldn't allocate the cumulative differential cross sections to estimate.\n");
+ }
}
- /* Allocate the phase function data */
- RESIZE(estimator->phase_functions, nwavelengths,
- "Couldn't allocate the per wavelength phase functions data.\n");
- FOR_EACH(i, 0, nwavelengths) {
- RESIZE(estimator->phase_functions[i].values, nangles,
- "Couldn't allocate the per angle phase function values.\n");
- RESIZE(estimator->phase_functions[i].cumulative, nangles,
- "Couldn't allocate the per angle cumulative phase function.\n");
+ if(!estimator->no_phase_function) {
+ /* Allocate the phase function data */
+ RESIZE(estimator->phase_functions, nwavelengths,
+ "Couldn't allocate the per wavelength phase functions data.\n");
+ FOR_EACH(i, 0, nwavelengths) {
+ /* Do not allocate the phase functions data if no limit angle was found */
+ if(estimator->limit_angles[i] == INVALID_LIMIT_ANGLE) continue;
+ RESIZE(estimator->phase_functions[i].values, nangles,
+ "Couldn't allocate the per angle phase function values.\n");
+ RESIZE(estimator->phase_functions[i].cumulative, nangles,
+ "Couldn't allocate the per angle cumulative phase function.\n");
+ }
}
#undef RESIZE
+
exit:
*out_estimator = estimator;
return res;
@@ -1461,6 +1503,9 @@ sschiff_integrate
MC_ACCUM(scattering_cross_section);
MC_ACCUM(avg_projected_area);
+ /* Discard differential cross section accumulation for invalid angle */
+ if(estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE) continue;
+
/* Accum up to "limit angle"; Great angles are analitically computed */
FOR_EACH(iangle, 0, estimator->limit_angles[iwlen]) {
MC_ACCUM(differential_cross_section[iangle]);
@@ -1469,13 +1514,12 @@ sschiff_integrate
}
}
-#if 1
/* TODO use parallelism */
FOR_EACH(i, 0, nwavelengths) {
- res = compute_phase_function(dev, i, nangles, estimator);
- if(res != RES_OK) goto error;
+ if(estimator->limit_angles[i] == INVALID_LIMIT_ANGLE) continue;
+ /* Do not handle phase function errors */
+ compute_phase_function(dev, i, nangles, estimator);
}
-#endif
exit:
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
@@ -1599,6 +1643,9 @@ sschiff_estimator_get_phase_function
{
if(!estimator || !states || iwlen >= sa_size(estimator->wavelengths))
return RES_BAD_ARG;
+ /* No phase function was computed */
+ if(estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE)
+ return RES_OK;
*states = estimator->phase_functions[iwlen].values;
return RES_OK;
}
@@ -1612,6 +1659,9 @@ sschiff_estimator_get_phase_function_cumulative
{
if(!estimator || !states || iwlen >= sa_size(estimator->wavelengths))
return RES_BAD_ARG;
+ /* No phase function cumulative was computed */
+ if(estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE)
+ return RES_OK;
*states = estimator->phase_functions[iwlen].cumulative;
return RES_OK;
}
@@ -1636,6 +1686,10 @@ sschiff_estimator_inverse_cumulative_phase_function
goto error;
}
+ /* No phase function cumulative is computed => nothing to inverse */
+ if(estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE)
+ goto exit;
+
/* Allocate the "filtered" phase function cumulative of small angles, i.e.
* the increasing cumulative defined from the Monte-Carlo estimation. */
nsmall_angles = estimator->limit_angles[iwlen];
diff --git a/src/test_sschiff_estimator_sphere.c b/src/test_sschiff_estimator_sphere.c
@@ -161,7 +161,7 @@ check_schiff_estimation
CHECK(eq_eps(val->E, results[i].scattering_cross_section, 4*val->SE), 1);
val = &cross_section.average_projected_area;
- printf(" Averavege projected area = %12.6g +/- %12.7g\n", val->E, val->SE);
+ printf(" Averavege projected area = %12.7g +/- %12.7g\n", val->E, val->SE);
CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
}
