commit 1da84d8ec665aeff9b43f4e20f6d4d2d0708fa2d
parent 525cbe51e6ffca609d823e51aeaabfa9061962d8
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 1 Mar 2016 16:28:08 +0100
Fix issues in solving the small/wide angle parameters
Use the cylindrical rather than the spherical Bessel functions and fix
up the derivative used in the Newton Solver.
Add a first draft of the cumulative phase function inversion.
Diffstat:
4 files changed, 257 insertions(+), 70 deletions(-)
diff --git a/src/sschiff.h b/src/sschiff.h
@@ -192,6 +192,13 @@ sschiff_estimator_get_cross_sections
(const struct sschiff_estimator* estimator,
struct sschiff_cross_section cross_sections[]);
+SSCHIFF_API res_T
+sschiff_estimator_inverse_cumulative_phase_function
+ (const struct sschiff_estimator* estimator,
+ const size_t iwlen,
+ double* thetas,
+ const size_t nthetas); /* Must be >= 2 */
+
END_DECLS
#endif /* SSCHIFF_H */
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -81,10 +81,11 @@ struct mc_accum {
/* Estimated phase function */
struct phase_function {
+ size_t ilimit_angle; /* Index of the firs wide angle */
+
/* Per angle values */
- double* values;
- double* cumulative;
- double* inverse_cumulative;
+ struct sschiff_state* values;
+ struct sschiff_state* cumulative;
};
/* The geometry context stores data that are constants for a given geometry
@@ -99,10 +100,12 @@ struct geometry_context {
* of the phase function is no more valid. Must be < nangles. */
size_t* limit_angles;
- double* medium_refractive_indices;
+ double* wavenumbers;
struct sschiff_material_properties* properties; /* Material properties */
};
+static struct geometry_context GEOMETRY_CONTEXT_NULL = {NULL, NULL, NULL, NULL};
+
/* The integrator context stores per thread data */
struct integrator_context {
/* Geometry of the micro particle */
@@ -153,6 +156,15 @@ cmp_double(const void* a, const void* b)
return i;
}
+static FINLINE int
+cmp_schiff_state(const void* a, const void* b)
+{
+ const struct sschiff_state* op0 = a;
+ const struct sschiff_state* op1 = b;
+ const int i = op0->E < op1->E ? -1 : (op0->E > op1->E ? 1 : 0);
+ return i;
+}
+
/* Compute an orthonormal basis where `dir' is the Z axis. */
static void
build_basis(const float dir[3], float basis[9])
@@ -391,7 +403,7 @@ accum_cross_section
double k_e;
double w_extinction, w_absorption, w_scattering;
- k_e = ctx->geom_ctx->medium_refractive_indices[iwlen];
+ k_e = ctx->geom_ctx->wavenumbers[iwlen];
n_r = ctx->geom_ctx->properties[iwlen].relative_real_refractive_index;
k_r = ctx->geom_ctx->properties[iwlen].relative_imaginary_refractive_index;
@@ -437,7 +449,7 @@ accum_differential_cross_section
size_t iangle;
/* Fetch optical properties */
- k_e = ctx->geom_ctx->medium_refractive_indices[iwlen];
+ k_e = ctx->geom_ctx->wavenumbers[iwlen];
n_r = ctx->geom_ctx->properties[iwlen].relative_real_refractive_index;
k_r = ctx->geom_ctx->properties[iwlen].relative_imaginary_refractive_index;
@@ -450,10 +462,10 @@ accum_differential_cross_section
beta_i[1] = k_e * k_r * path_length[1];
tmp = (k_e * plane_area) / (2*PI);
tmp *= tmp;
- tmp *= 1
+ tmp *= (1
+ exp(-(beta_i[0] + beta_i[1])) * cos(beta_r[0] - beta_r[1])
- exp(-beta_i[0]) * cos(beta_r[0])
- - exp(-beta_i[1]) * cos(beta_r[1]);
+ - exp(-beta_i[1]) * cos(beta_r[1]));
/* Compute and accumulate the MC weights of the differential cross section
* and its cumulative. Note that ctx->limit_angles store the index of the
@@ -466,11 +478,11 @@ accum_differential_cross_section
k_e_angle_delta_r = ctx->geom_ctx->angles[iangle] * k_e_delta_r;
- bessel_j0 = gsl_sf_bessel_j0(k_e_angle_delta_r);
+ bessel_j0 = gsl_sf_bessel_J0(k_e_angle_delta_r);
weight = bessel_j0 * tmp;
ctx->accums[iwlen].differential_cross_section[iangle] += weight;
- bessel_j1 = gsl_sf_bessel_j1(k_e_angle_delta_r);
+ bessel_j1 = gsl_sf_bessel_J1(k_e_angle_delta_r);
weight = 2*PI * ctx->geom_ctx->angles[iangle] * bessel_j1 / k_e_delta_r * tmp;
ctx->accums[iwlen].differential_cross_section_cumulative[iangle] += weight;
}
@@ -517,8 +529,8 @@ function_fdf(const gsl_vector* vec, void* params, gsl_vector* res, gsl_matrix* J
u_prime = arg->sqr_sin_half_angle
* ((2*n - 6) * arg->cos_2_angle + 6*n - 8*arg->cos_angle - 26)
- - n/2.0 * arg->cos_half_angle * pow(arg->sin_half_angle, n-1)
- * (4*n*n - 24*n + 40) - pow(arg->sin_half_angle, n) * (8*n - 24);
+ - pow(arg->sin_half_angle, n)
+ * (log(arg->sin_half_angle) * (4*n*n - 24*n + 40) + 8*n-24);
v_prime = (n-4)*(n-6) + (n-2)*(n-6) + (n-2)*(n-4);
f = arg->differential_cross_section_cumulative
@@ -603,7 +615,7 @@ solve_wide_angles_phase_function_model_parameters
/* Initialise the solver. TODO perform this in an initialisation step of a
* phase_function context */
- gsl_vector_set(init_val, 0, 3.0); /* n */
+ gsl_vector_set(init_val, 0, 3);
gsl_multiroot_fdfsolver_set(solver, &func, init_val);
gsl_set_error_handler_off();
@@ -614,10 +626,15 @@ solve_wide_angles_phase_function_model_parameters
do {
status = gsl_multiroot_fdfsolver_iterate(solver);
if(status == GSL_SUCCESS) {
- status = gsl_multiroot_test_residual(solver->f, 1.e-3);
+ status = gsl_multiroot_test_residual(solver->f, 1.e-6);
}
} while(status == GSL_CONTINUE && ++iter < max_iter);
+ /* Retrieve the estimated "n" value */
+ root = gsl_multiroot_fdfsolver_root(solver);
+ *n = gsl_vector_get(root, 0);
+ printf("theta_limit = %g, n = %g\n", arg.angle, *n); /* TODO remove this */
+
if(status != GSL_SUCCESS) {
/* TODO enhance the error message. */
log_error(dev, "Error in Newton estimation - GSL error: %s - %d.\n",
@@ -626,11 +643,6 @@ solve_wide_angles_phase_function_model_parameters
goto error;
}
- /* Retrieve the estimated "n" value */
- root = gsl_multiroot_fdfsolver_root(solver);
- *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)
/ (1 + arg.sqr_cos_angle);
@@ -674,8 +686,8 @@ compute_phase_function
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;
+ struct sschiff_state limit_differential_cross_section;
+ struct sschiff_state 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;
@@ -692,12 +704,21 @@ compute_phase_function
scattering_cross_section = expected_value
(&estimator->accums[iwlen].scattering_cross_section,
estimator->nrealisations);
- limit_differential_cross_section = expected_value
+ get_mc_value
(&estimator->accums[iwlen].differential_cross_section[ilimit_angle],
- estimator->nrealisations);
- limit_differential_cross_section_cumulative = expected_value
+ estimator->nrealisations,
+ &limit_differential_cross_section);
+ get_mc_value
(&estimator->accums[iwlen].differential_cross_section_cumulative[ilimit_angle],
- estimator->nrealisations);
+ estimator->nrealisations,
+ &limit_differential_cross_section_cumulative);
+ printf("sigma_s = %g, Ws = %g +/- %g, Wc = %g +/- %g\n",
+ scattering_cross_section,
+ limit_differential_cross_section.E,
+ limit_differential_cross_section.SE,
+ limit_differential_cross_section_cumulative.E,
+ limit_differential_cross_section_cumulative.SE);
+
res = solve_wide_angles_phase_function_model_parameters
(dev,
limit_angle,
@@ -705,8 +726,8 @@ compute_phase_function
cos_limit_angle,
cos_2_limit_angle,
scattering_cross_section,
- limit_differential_cross_section,
- limit_differential_cross_section_cumulative,
+ limit_differential_cross_section.E,
+ limit_differential_cross_section_cumulative.E,
&n, &r);
if(res != RES_OK) goto error;
@@ -737,7 +758,7 @@ compute_phase_function
(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);
+ limit_differential_cross_section_cumulative.E + 2*PI*r*(coef-coef_limit);
}
/* Check the post condition of the cumulative differential cross section */
@@ -745,35 +766,50 @@ compute_phase_function
(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;
+ /* Compute the [cumulative] phase function for short angles */
+ rcp_scattering_cross_section = 1.0/scattering_cross_section;
FOR_EACH(iangle, 0, ilimit_angle) {
- estimator->phase_functions[iwlen].values[iangle] = expected_value
+ get_mc_value
(&estimator->accums[iwlen].differential_cross_section[iangle],
- estimator->nrealisations);
- estimator->phase_functions[iwlen].values[iangle] *=
+ estimator->nrealisations,
+ &estimator->phase_functions[iwlen].values[iangle]);
+ estimator->phase_functions[iwlen].values[iangle].E *=
rcp_scattering_cross_section;
- estimator->phase_functions[iwlen].cumulative[iangle] = expected_value
+ get_mc_value
(&estimator->accums[iwlen].differential_cross_section_cumulative[iangle],
- estimator->nrealisations);
- estimator->phase_functions[iwlen].cumulative[iangle] *=
+ estimator->nrealisations,
+ &estimator->phase_functions[iwlen].cumulative[iangle]);
+ estimator->phase_functions[iwlen].cumulative[iangle].E *=
rcp_scattering_cross_section;
}
+ /* Compute the [cumulative] phase function for wide angles */
FOR_EACH(iangle, ilimit_angle + 1, nangles) {
- estimator->phase_functions[iwlen].values[iangle] =
+ estimator->phase_functions[iwlen].values[iangle].E =
estimator->accums[iwlen].differential_cross_section[iangle].weight
* rcp_scattering_cross_section;
- estimator->phase_functions[iwlen].values[iangle] =
+ estimator->phase_functions[iwlen].cumulative[iangle].E =
estimator->accums[iwlen].differential_cross_section_cumulative[iangle].weight
* rcp_scattering_cross_section;
+
+ /* THe phase function for wide angles is analitically computed, i.e. there
+ * is no variance or standard error */
+ estimator->phase_functions[iwlen].values[iangle].V = -1;
+ estimator->phase_functions[iwlen].values[iangle].SE = -1;
+ estimator->phase_functions[iwlen].cumulative[iangle].V = -1;
+ estimator->phase_functions[iwlen].cumulative[iangle].SE = -1;
}
+
+ /* Setup the [cumulative] phase function for the limit angle */
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;
+ limit_differential_cross_section;
+ estimator->phase_functions[iwlen].values[ilimit_angle].E *=
+ rcp_scattering_cross_section;
+ estimator->phase_functions[iwlen].cumulative[ilimit_angle] =
+ limit_differential_cross_section_cumulative;
+ estimator->phase_functions[iwlen].cumulative[ilimit_angle].E *=
+ rcp_scattering_cross_section;
exit:
return res;
@@ -1000,6 +1036,64 @@ error:
goto exit;
}
+static res_T
+inverse_cumulative_phase_function
+ (const struct sschiff_estimator* estimator,
+ const double* cumulative_small_angles,
+ const size_t iwlen,
+ const double cumulative,
+ double* theta)
+{
+ double u, lower, upper;
+ size_t iangle, nsmall_angles;
+ ASSERT(estimator && cumulative_small_angles && theta);
+ ASSERT(cumulative >= 0.0 && cumulative <= 1.0);
+ ASSERT(iwlen < sa_size(estimator->phase_functions));
+
+ nsmall_angles = estimator->phase_functions[iwlen].ilimit_angle;
+ if(cumulative < cumulative_small_angles[nsmall_angles-1]) {
+ /* Look for the the cumulative in the filtered small angles array */
+ double* found = search_lower_bound(&cumulative, cumulative_small_angles,
+ nsmall_angles, sizeof(double), cmp_double);
+ if(!found) {
+ FATAL("Error in inverting the phase function cumulative for small angles.\n");
+ }
+ iangle = (size_t)(found - cumulative_small_angles);
+ upper = cumulative_small_angles[iangle];
+ lower = cumulative_small_angles[iangle-1];
+ } else {
+ /* Look for the cumulative in the wide angles cumulative */
+ struct sschiff_state cumul;
+ struct sschiff_state* found = NULL;
+ size_t nangles;
+ cumul.E = cumulative;
+
+ nangles = sa_size(estimator->angles);
+ found = search_lower_bound
+ (&cumul,
+ estimator->phase_functions[iwlen].cumulative + nsmall_angles,
+ nangles - nsmall_angles,
+ sizeof(struct sschiff_state), cmp_schiff_state);
+ if(!found) {
+ FATAL("Error in inverting the phase function cumulative for wide angles.\n");
+ }
+ iangle = (size_t)(found - estimator->phase_functions[iwlen].cumulative);
+ ASSERT(iangle >= nsmall_angles);
+ upper = estimator->phase_functions[iwlen].cumulative[iangle].E;
+ if(iangle == nsmall_angles) {
+ lower = cumulative_small_angles[nsmall_angles-1];
+ } else {
+ lower = estimator->phase_functions[iwlen].cumulative[iangle-1].E;
+ }
+ }
+
+ /* Use the cumulative bounds to linearly interpolate the angles */
+ u = (cumulative - lower) / (upper - lower);
+ *theta = u*estimator->angles[iangle] + (1.0 - u)*estimator->angles[iangle-1];
+
+ return RES_OK;
+}
+
static char
check_distribution(struct sschiff_geometry_distribution* distrib)
{
@@ -1012,10 +1106,11 @@ geometry_context_release(struct geometry_context* ctx)
{
ASSERT(ctx);
if(ctx->limit_angles) sa_release(ctx->limit_angles);
- if(ctx->medium_refractive_indices) sa_release(ctx->medium_refractive_indices);
+ if(ctx->wavenumbers) sa_release(ctx->wavenumbers);
if(ctx->properties) sa_release(ctx->properties);
}
+/* Return RES_BAD_OP if the limit angle of the phase function is invalid */
static res_T
geometry_context_init
(struct sschiff_device* dev,
@@ -1028,6 +1123,7 @@ geometry_context_init
{
double rcp_PI_Lc;
size_t iwlen;
+ int skip_differential_cross_section = 0;
res_T res = RES_OK;
ASSERT(dev && distrib && wavelengths && nwavelengths && angles && nangles && ctx);
memset(ctx, 0, sizeof(struct geometry_context));
@@ -1041,9 +1137,8 @@ geometry_context_init
res = RES_MEM_ERR;
goto error;
}
- if(!sa_add(ctx->medium_refractive_indices, nwavelengths)) {
- log_error(dev,
-"Couldn't allocate the per wavelength refractive indices of the medium.\n");
+ if(!sa_add(ctx->wavenumbers, nwavelengths)) {
+ log_error(dev, "Couldn't allocate the list of wavenumbers.\n");
res = RES_MEM_ERR;
goto error;
}
@@ -1064,19 +1159,32 @@ geometry_context_init
/* Precompute the refractive index in the medium */
lambda_e = wavelengths[iwlen] / ctx->properties[iwlen].medium_refractive_index;
- ctx->medium_refractive_indices[iwlen] = 2.0*PI / lambda_e;
+ ctx->wavenumbers[iwlen] = 2.0*PI / lambda_e;
+
+ if(skip_differential_cross_section)
+ continue;
/* Search the limit angle */
theta_l = sqrt(lambda_e * rcp_PI_Lc);
+ if(theta_l <= 0 || theta_l >= PI) {
+ log_error(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, wavelengths[iwlen]);
+ res = RES_BAD_ARG;
+ goto error;
+ }
angle = search_lower_bound
(&theta_l, angles, nangles, sizeof(double), cmp_double);
if(!angle) {
- log_error(dev, "Can't find a limit angle for theta = %g.\n", theta_l);
+ 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 angle \"%g\".\n", *angle);
+ log_error(dev, "Invalid limit scattering angle `%g'.\n", *angle);
res = RES_BAD_ARG;
goto error;
}
@@ -1086,12 +1194,17 @@ geometry_context_init
ASSERT(ctx->limit_angles[iwlen] < nangles);
}
+ if(skip_differential_cross_section) {
+
+ }
+
ctx->angles = angles;
exit:
return res;
error:
geometry_context_release(ctx);
+ *ctx = GEOMETRY_CONTEXT_NULL;
goto exit;
}
@@ -1242,8 +1355,6 @@ estimator_release(ref_T* ref)
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);
}
@@ -1356,12 +1467,6 @@ estimator_create
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:
@@ -1391,7 +1496,7 @@ sschiff_integrate
const size_t ndirs,
struct sschiff_estimator** out_estimator)
{
- struct geometry_context geom_ctx;
+ struct geometry_context geom_ctx = GEOMETRY_CONTEXT_NULL;
struct integrator_context* ctxs = NULL;
struct sschiff_estimator* estimator = NULL;
struct ssp_rng** rngs = NULL;
@@ -1504,8 +1609,8 @@ sschiff_integrate
}
}
+#if 1
/* TODO use parallelism */
-#if 0
FOR_EACH(i, 0, nwavelengths) {
res = compute_phase_function(dev, i, nangles, &geom_ctx, estimator);
if(res != RES_OK) goto error;
@@ -1624,3 +1729,78 @@ sschiff_estimator_get_cross_sections
return RES_OK;
}
+res_T
+sschiff_estimator_inverse_cumulative_phase_function
+ (const struct sschiff_estimator* estimator,
+ const size_t iwlen,
+ double* thetas,
+ const size_t nthetas)
+{
+ struct sschiff_state* cumul[2];
+ double* cumulative_small_angles = NULL;
+ double cumulative;
+ double step;
+ size_t itheta;
+ size_t iangle, nsmall_angles;
+ res_T res = RES_OK;
+
+ if(!estimator || !nthetas < 2 || !thetas) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Allocate the "filtered" phase function cumulative of small angles, i.e.
+ * the increasing cumulative defined from the Monte-Carlo estimation. */
+ nsmall_angles = estimator->phase_functions[iwlen].ilimit_angle;
+ if(!sa_add(cumulative_small_angles, nsmall_angles)) {
+ log_error(estimator->dev,
+"Couldn't allocate the filtered phase function cumulative of small angles.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ /* Fill the phase function cumulative for small angles. Ensure that the
+ * resulting array is sorted in increasing order, i.e. increasing function */
+ cumul[0] = &estimator->phase_functions[iwlen].cumulative[0];
+ cumulative_small_angles[0] = cumul[0]->E;
+ FOR_EACH(iangle, 1, nsmall_angles) {
+ cumul[1] = &estimator->phase_functions[iwlen].cumulative[1];
+
+ if(cumul[0]->E <= cumul[1]->E) {
+ cumulative_small_angles[iangle] = cumul[1]->E;
+ } else {
+ /* Numerical imprecisions may lead to a cumulative that is not an
+ * increasing function. In such case, check that the standard-error
+ * ensure at least a threshold between the decreasing entry an the
+ * previous one. If not, return an error */
+ cumulative = MMIN(cumul[1]->E + cumul[1]->SE, cumul[0]->E);
+ if(cumulative < cumul[0]->E) {
+ log_error(estimator->dev,
+"The phase function cumulative of small angles is not an increasing function.\n"
+"This may be due to numerical imprecisions.\n");
+ res = RES_BAD_OP;
+ goto error;
+ }
+ cumulative_small_angles[iangle] = cumulative;
+ }
+ cumul[0] = cumul[1];
+ }
+
+ /* Inverse the phase function cumulative */
+ thetas[0] = 0.0;
+ thetas[nthetas-1] = PI;
+ cumulative = step = 1.0/(double)nthetas;
+ FOR_EACH(itheta, 1, nthetas-1) {
+ /* TODO since the submitted cumulative are strictly increasing, one can
+ * speed up the inversion by using the previous search result to reduce the
+ * search domain in the cumulative arrays */
+ inverse_cumulative_phase_function
+ (estimator, cumulative_small_angles, iwlen, cumulative, &thetas[itheta]);
+ }
+
+exit:
+ if(cumulative_small_angles) sa_release(cumulative_small_angles);
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/test_sschiff_estimator_cylinder.c b/src/test_sschiff_estimator_cylinder.c
@@ -441,7 +441,7 @@ main(int argc, char** argv)
sampler_ctx.mean_radius = (x[i] * 0.450) / (2*PI);
sampler_ctx.sigma = 1.18;
- distrib.caracteristic_length = 1;
+ distrib.caracteristic_length = sampler_ctx.mean_radius;
distrib.material.get_property = get_material_property;
distrib.material.material = &sampler_ctx;
distrib.sample = sample_cylinder;
diff --git a/src/test_sschiff_estimator_rhodo.c b/src/test_sschiff_estimator_rhodo.c
@@ -48,14 +48,14 @@ get_material_property
struct sschiff_material_properties* props)
{
(void)mtl;
-#if 0
- props->medium_refractive_index = 1.34177984360538;
- props->relative_imaginary_refractive_index = 2.20382600910679e-03;
- props->relative_real_refractive_index = 1.38912721283454;
+#if 1
+ props->medium_refractive_index = 1.341779;
+ props->relative_imaginary_refractive_index = 2.20382600910679e-03/1.341779;
+ props->relative_real_refractive_index = 1.38912721283454/1.341779;
#else
- props->medium_refractive_index = 1.331073;
- props->relative_imaginary_refractive_index = 8.53534081490802e-04;
- props->relative_real_refractive_index = 1.39071033668243;
+ props->medium_refractive_index = 1.33;
+ props->relative_imaginary_refractive_index = 8.5e-4/1.33;
+ props->relative_real_refractive_index = 1.39/1.33;
#endif
}
@@ -86,16 +86,16 @@ main(int argc, char** argv)
struct sschiff_estimator* estimator;
struct sschiff_cross_section cross_section;
struct ssp_rng* rng;
- const double wavelength = 0.6; /* In micron */
+ const double wavelength = 0.4; /* In micron */
const size_t nscatt_angles = 1000;
const size_t ngeoms = 10000;
- const size_t ndirs = 10;
+ const size_t ndirs = 100;
(void)argc, (void)argv;
mem_init_proxy_allocator(&allocator, &mem_default_allocator);
CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng), RES_OK);
CHECK(sschiff_device_create
- (NULL, &allocator, 1, 1, NULL, &dev), RES_OK);
+ (NULL, &allocator, SSCHIFF_NTHREADS_DEFAULT, 1, NULL, &dev), RES_OK);
geometry_init_cylinder(&sampler_ctx.geometry, 64);
