star-schiff

sschiff_estimator.c (65003B)

---

 /* Copyright (C) 2015, 2016, 2026 Centre National de la Recherche Scientifique
  * Copyright (C) 2026 Clermont Auvergne INP
  * Copyright (C) 2026 Institut Mines Télécom Albi-Carmaux
  * Copyright (C) 2020, 2021, 2023, 2026 |Méso|Star> (contact@meso-star.com)
  * Copyright (C) 2026 Université de Lorraine
  * Copyright (C) 2026 Université de Toulouse
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation, either version 3 of the License, or
  * (at your option) any later version.
  *
  * This program is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  * GNU General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
 
 #define _POSIX_C_SOURCE 200112L /* nextafterf support */
 
 #include "sschiff.h"
 #include "sschiff_device.h"
 
 #include <rsys/algorithm.h>
 #include <rsys/float2.h>
 #include <rsys/float3.h>
 #include <rsys/image.h>
 #include <rsys/logger.h>
 #include <rsys/math.h>
 #include <rsys/mem_allocator.h>
 #include <rsys/stretchy_array.h>
 
 #include <star/s3d.h>
 #include <star/ssp.h>
 
 #include <math.h>
 #include <omp.h>
 
 #include <gsl/gsl_sf.h>
 #include <gsl/gsl_multiroots.h>
 
 #define MAX_FAILURES 10
 #define INVALID_LIMIT_ANGLE SIZE_MAX
 
 /* Accumulator of double precision values */
 struct accum {
   double extinction_cross_section;
   double absorption_cross_section;
   double scattering_cross_section;
   double avg_projected_area;
 
   /* Per scattering angles weights */
   double* differential_cross_section;
   double* differential_cross_section_cumulative;
 };
 
 /* Monte carlo data */
 struct mc_data { double weight; double sqr_weight; };
 static const struct mc_data MC_DATA_NULL = { 0.0, 0.0 };
 
 /* Accumulator of Monte Carlo data */
 struct mc_accum {
   struct mc_data extinction_cross_section;
   struct mc_data absorption_cross_section;
   struct mc_data scattering_cross_section;
   struct mc_data avg_projected_area;
 
   /* Per scattering angles weight */
   struct mc_data* differential_cross_section;
   struct mc_data* differential_cross_section_cumulative;
 };
 
 /* Estimated phase function */
 struct phase_function {
   /* Per angle values */
   struct sschiff_state* values;
   struct sschiff_state* cumulative;
 };
 
 /* The integrator context stores per thread data */
 struct integrator_context {
   /* Geometry of the micro particle */
   struct s3d_shape* shape;
   struct s3d_scene* scene;
   struct s3d_scene_view* view;
   void* shape_data; /* Client side data of the shape */
 
   struct ssp_rng* rng; /* Random Number Generator */
   size_t nwavelengths; /* #wavelengths */
   size_t nangles; /* #scattering angles */
 
    /* Per wavelengths data */
   struct accum* accums; /* Temporary accumulators */
   struct mc_accum* mc_accums; /* Monte Carlo accumulators */
 
   /* Pointer toward the estimator */
   struct sschiff_estimator* estimator;
 };
 
 /* Store per thread data of the post integration process */
 struct solver_context {
   gsl_multiroot_fdfsolver* solver;
   gsl_vector* init_val;
 };
 
 struct sschiff_estimator {
   double* wavelengths; /* List of wavelengths in vacuum */
   double* angles; /* List of scattering angles */
 
   /* Per wavelength index of the scattering angle from which the MC estimation
    * of the phase function is no more valid. Must be < nangles. */
   size_t* limit_angles;
   double* wide_angles_parameter;
 
   double* wavenumbers;
   struct sschiff_material_properties* properties; /* Material properties */
 
   /* Per wavelength results */
   struct phase_function* phase_functions;
   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 */
   int discard_large_angles; /* Avoid analytic model for wide angles */
 
   struct sschiff_device* dev;
   ref_T ref;
 };
 
 /*******************************************************************************
  * Helper functions
  ******************************************************************************/
 static FINLINE int
 cmp_double(const void* a, const void* b)
 {
   const double* op0 = a;
   const double* op1 = b;
   const int i = *op0 < *op1 ? -1 : (*op0 > *op1 ? 1 : 0);
   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])
 {
   float dir_abs[3];
   float axis_min[3];
   int i;
   ASSERT(dir && basis && f3_is_normalized(dir));
 
   /* Define which axis direction is minimal and use its unit vector to compute
    * a vector ortogonal to `dir'. This ensures that the unit vector and `dir'
    * are not colinear and thus that their cross product is not a zero vector. */
   FOR_EACH(i, 0, 3) dir_abs[i] = absf(dir[i]);
   if(dir_abs[0] < dir_abs[1]) {
     if(dir_abs[0] < dir_abs[2]) f3(axis_min, 1.f, 0.f, 0.f);
     else f3(axis_min, 0.f, 0.f, 1.f);
   } else {
     if(dir_abs[1] < dir_abs[2]) f3(axis_min, 0.f, 1.f, 0.f);
     else f3(axis_min, 0.f, 0.f, 1.f);
   }
 
   f3_normalize(basis + 0, f3_cross(basis + 0, dir, axis_min));
   f3_cross(basis + 3, dir, basis + 0);
   f3_set(basis + 6, dir);
 }
 
 /* Compute a 3x4 affine transformation from a orthonormal basis and the
  * position of the origin */
 static FINLINE void
 build_transform
   (const float pos[3], /* Position of the origin */
    const float basis[9], /* Column major */
    float transform[12]) /* Column major */
 {
   ASSERT(pos && basis && transform);
 
    /* The linear part of the transform matrix is the inverse of the
     * orthonormal basis (i.e. its transpose) */
   f3(transform + 0, basis[0], basis[3], basis[6]);
   f3(transform + 3, basis[1], basis[4], basis[7]);
   f3(transform + 6, basis[2], basis[5], basis[8]);
 
   /* Affine part of the transform matrix */
   transform[9 ] = -f3_dot(basis+0, pos);
   transform[10] = -f3_dot(basis+3, pos);
   transform[11] = -f3_dot(basis+6, pos);
 }
 
 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 = 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);
 }
 
 static FINLINE int
 hit_on_edge(const struct s3d_hit* hit)
 {
   const float on_edge_eps = 1.e-4f;
   float w;
   ASSERT(hit && !S3D_HIT_NONE(hit));
   w = 1.f - hit->uv[0] - hit->uv[1];
   return eq_epsf(hit->uv[0], 0.f, on_edge_eps)
       || eq_epsf(hit->uv[0], 1.f, on_edge_eps)
       || eq_epsf(hit->uv[1], 0.f, on_edge_eps)
       || eq_epsf(hit->uv[1], 1.f, on_edge_eps)
       || eq_epsf(w, 0.f, on_edge_eps)
       || eq_epsf(w, 1.f, on_edge_eps);
 }
 
 static FINLINE int
 self_hit
   (const struct s3d_hit* hit_from,
    const struct s3d_hit* hit_to,
    const float epsilon)
 {
   ASSERT(hit_from && hit_to);
 
   if(S3D_HIT_NONE(hit_from) || S3D_HIT_NONE(hit_to))
     return 0;
   if(S3D_PRIMITIVE_EQ(&hit_from->prim, &hit_to->prim))
     return 1;
   if(eq_epsf(hit_to->distance - hit_from->distance, 0.f, epsilon))
     return hit_on_edge(hit_from) && hit_on_edge(hit_to);
   return 0;
 }
 
 /* Compute the length in micron of the ray part that traverses the scene */
 static res_T
 compute_path_length
   (struct sschiff_device* dev,
    struct s3d_scene_view* view,
    const struct s3d_hit* first_hit,
    const float ray_org[3],
    const float ray_dir[3],
    const double scale,
    double* length)
 {
   float range[2] = { 0, FLT_MAX };
   struct s3d_hit hit_from, hit;
   double len = 0;
   float dst;
   ASSERT(view && first_hit && ray_org && ray_dir && !S3D_HIT_NONE(first_hit));
   ASSERT(scale > 0);
   (void)dev;
 
   dst = range[0] = first_hit->distance;
   hit_from = *first_hit;
 
   do {
     do {
       range[0] = nextafterf(range[0], range[1]);
       S3D(scene_view_trace_ray(view, ray_org, ray_dir, range, NULL, &hit));
     } while(hit.distance < range[0] || self_hit(&hit_from, &hit, 1.e-6f));
 
     if(S3D_HIT_NONE(&hit))
       return RES_BAD_OP;
 
     len += hit.distance - dst;
     dst = range[0] = hit.distance;
     hit_from = hit;
 
     do {
       range[0] = nextafterf(range[0], range[1]);
       S3D(scene_view_trace_ray(view, ray_org, ray_dir, range, NULL, &hit));
     } while(hit.distance < range[0] || self_hit(&hit_from, &hit, 1.e-6f));
 
     range[0] = hit.distance;
     hit_from = hit;
 
   } while(!S3D_HIT_NONE(&hit));
 
   *length = len * scale;
   return RES_OK;
 }
 
 static INLINE void
 accum_cross_section
   (struct integrator_context* ctx,
    const double plane_area,
    const double path_length)
 {
   size_t iwlen;
   ASSERT(ctx && plane_area > 0 && path_length > 0);
 
   FOR_EACH(iwlen, 0, ctx->nwavelengths) {
     double n_r, k_r;
     double k_e;
     double w_extinction, w_absorption, w_scattering;
 
     k_e = ctx->estimator->wavenumbers[iwlen];
     n_r = ctx->estimator->properties[iwlen].relative_real_refractive_index;
     k_r = ctx->estimator->properties[iwlen].relative_imaginary_refractive_index;
 
     /* Extinction cross section */
     w_extinction = 2.0 * plane_area *
       (1.0 - exp(-k_e*k_r*path_length) * cos(k_e*(n_r - 1.0)*path_length));
     ctx->accums[iwlen].extinction_cross_section += w_extinction;
 
     /* Absorption cross section */
     w_absorption = plane_area * (1.0 - exp(-2.0*k_e*k_r*path_length));
     ctx->accums[iwlen].absorption_cross_section += w_absorption;
 
     /* Scattering cross section */
     w_scattering = w_extinction - w_absorption;
     ctx->accums[iwlen].scattering_cross_section += w_scattering;
 
     /* Average projected area */
     ctx->accums[iwlen].avg_projected_area += plane_area;
   }
 }
 
 static INLINE void
 accum_differential_cross_section
   (struct integrator_context* ctx,
    const double plane_area,
    const double path_length[2],
    const double delta_r)/* Length between the 2 ray origins in footprint space */
 {
   size_t iwlen;
 
   FOR_EACH(iwlen, 0, ctx->nwavelengths) {
     double k_e_delta_r;
     double n_r, k_r;
     double k_e;
     double beta_r[2];
     double beta_i[2];
     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;
     k_r = ctx->estimator->properties[iwlen].relative_imaginary_refractive_index;
 
     /* Precompute some values. TODO can be stored on the context */
     k_e_delta_r = k_e * delta_r;
     beta_r[0] = k_e * (n_r - 1) * path_length[0];
     beta_r[1] = k_e * (n_r - 1) * path_length[1];
     beta_i[0] = k_e * k_r * path_length[0];
     beta_i[1] = k_e * k_r * path_length[1];
     tmp  = (k_e * plane_area) / (2*PI);
     tmp *= tmp;
     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]));
 
     /* Compute and accumulate the MC weights of the differential cross section
      * and its cumulative. Note that ctx->limit_angles store the index of the
      * first scattering angle whose phase function is not estimated by
      * Monte-Carlo */
     FOR_EACH(iangle, 0, ctx->estimator->limit_angles[iwlen]) {
       double bessel_j0, bessel_j1;
       double weight;
       double k_e_angle_delta_r;
 
       k_e_angle_delta_r = ctx->estimator->angles[iangle] * k_e_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);
       weight = 2*PI*ctx->estimator->angles[iangle]*bessel_j1/k_e_delta_r*tmp;
       ctx->accums[iwlen].differential_cross_section_cumulative[iangle] += weight;
     }
   }
 }
 
 static res_T
 accum_monte_carlo_weight
   (struct sschiff_device* dev,
    struct integrator_context* ctx,
    const float basis[9], /* World space basis of the RT volume */
    const float pos[3], /* World space centroid of the RT volume */
    const float lower[3], /* Lower boundary of the RT volume */
    const float upper[3], /* Upper boundary of the RT volume */
    const double area_scaling, /* Scale factor to apply to the areas */
    const double dist_scaling) /* Scale factor to apply to the distances */
 {
   struct s3d_hit hit;
   double sample[2];
   const float range[2] = { 0, FLT_MAX };
   size_t nfailures = 0;
   float axis_x[3], axis_y[3], axis_z[3];
   float plane_size[2];
   float ray_org[2][3];
   float org[3];
   float x[3], y[3];
   float plane_area;
   res_T res = RES_OK;
   ASSERT(dev && ctx &&  basis && pos && lower && upper);
   ASSERT(area_scaling > 0 && dist_scaling > 0);
 
   f2_sub(plane_size, upper, lower); /* In micron */
   plane_area = plane_size[0] * plane_size[1] * (float)area_scaling;
 
   /* Define the projection axis */
   f3_mulf(axis_x, basis + 0, plane_size[0] * 0.5f);
   f3_mulf(axis_y, basis + 3, plane_size[1] * 0.5f);
   f3_set(axis_z, basis + 6);
 
   /* Compute the origin of the projection plane */
   f3_add(org, pos, f3_mulf(org, axis_z, lower[2]));
 
   do {
     float vec[3];
     double path_length[2];
     double delta_r;
 
     /* Uniformly sample a position onto the projection plane and use it as the
      * origin of the 1st ray to trace */
     sample[0] = ssp_rng_uniform_double(ctx->rng, -1.0, 1.0);
     sample[1] = ssp_rng_uniform_double(ctx->rng, -1.0, 1.0);
     f3_mulf(x, axis_x, (float)sample[0]);
     f3_mulf(y, axis_y, (float)sample[1]);
     f3_add(ray_org[0], f3_add(ray_org[0], x, y), org);
 
     S3D(scene_view_trace_ray(ctx->view, ray_org[0], axis_z, range, NULL, &hit));
     /* NULL cross section and differential cross section weight */
     if(S3D_HIT_NONE(&hit))
       break;
 
     res = compute_path_length(dev, ctx->view, &hit, ray_org[0], axis_z,
       dist_scaling, &path_length[0]);
     if(res != RES_OK) { /* Handle numerical/geometry issues */
       ++nfailures;
       continue;
     }
     /* 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);
     sample[1] = ssp_rng_uniform_double(ctx->rng, -1.0, 1.0);
     f3_mulf(x, axis_x, (float)sample[0]);
     f3_mulf(y, axis_y, (float)sample[1]);
     f3_add(ray_org[1], f3_add(ray_org[1], x, y), org);
 
     S3D(scene_view_trace_ray(ctx->view, ray_org[1], axis_z, range, NULL, &hit));
     if(S3D_HIT_NONE(&hit)) /* NULL differential cross section weight */
       break;
 
     res = compute_path_length(dev, ctx->view, &hit, ray_org[1], axis_z,
       dist_scaling, &path_length[1]);
     if(res != RES_OK) { /* Handle numerical/geometry issues */
       ++nfailures;
       continue;
     }
     /* Compute and accumulate the per scattering angle differential cross
      * section weight */
     delta_r = f3_len(f3_sub(vec, ray_org[0], ray_org[1]));
     delta_r *= dist_scaling;
     accum_differential_cross_section(ctx, plane_area, path_length, delta_r);
 
   } while(res != RES_OK && nfailures < MAX_FAILURES);
 
   if(nfailures < MAX_FAILURES) {
     res = RES_OK;
   } else {
     log_error(dev,
 "Too many failures in computing the radiative path length. The sampled geometry\n"
 "might not defined a closed volume.\n");
   }
   return res;
 }
 
 static res_T
 radiative_properties
   (struct sschiff_device* dev,
    const int istep,
    const size_t ndirs,
    struct sschiff_geometry_distribution* distrib,
    struct integrator_context* ctx)
 {
   float lower[3], upper[3];
   float aabb_pt[8][3];
   float centroid[3];
   size_t idir;
   size_t iwlen;
   int i;
   res_T res = RES_OK;
   ASSERT(dev && ndirs && distrib && ctx);
   (void)istep;
 
   S3D(scene_view_get_aabb(ctx->view, lower, upper));
 
   /* AABB vertex layout
    *    6-------7
    *   /'      /|
    *  2-------3 |
    *  | 4.....|.5
    *  |,      |/
    *  0-------1 */
   f3(aabb_pt[0], lower[0], lower[1], lower[2]);
   f3(aabb_pt[1], upper[0], lower[1], lower[2]);
   f3(aabb_pt[2], lower[0], upper[1], lower[2]);
   f3(aabb_pt[3], upper[0], upper[1], lower[2]);
   f3(aabb_pt[4], lower[0], lower[1], upper[2]);
   f3(aabb_pt[5], upper[0], lower[1], upper[2]);
   f3(aabb_pt[6], lower[0], upper[1], upper[2]);
   f3(aabb_pt[7], upper[0], upper[1], upper[2]);
 
   f3_mulf(centroid, f3_add(centroid, lower, upper), 0.5f);
 
   /* Clear direction accumulators */
   FOR_EACH(iwlen, 0, ctx->nwavelengths) {
     ctx->accums[iwlen].extinction_cross_section = 0;
     ctx->accums[iwlen].absorption_cross_section = 0;
     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]);
     memset(ctx->accums[iwlen].differential_cross_section_cumulative, 0,
       sizeof(double)*ctx->estimator->limit_angles[iwlen]);
   }
 
   FOR_EACH(idir, 0, ndirs) {
     float dir[3];
     float basis[9];
     float transform[12];
     float pt[8][3];
     double volume_scaling = 1.0;
     double area_scaling = 1.0;
     double dist_scaling = 1.0;
 
     /* Sample a volume scaling factor if necessary */
     if(distrib->sample_volume_scaling) {
       res = distrib->sample_volume_scaling
         (ctx->rng, ctx->shape_data, &volume_scaling, distrib->context);
       if(res != RES_OK) {
         log_error(dev, "Error in sampling a volume scaling.\n");
         goto error;
       }
       if(volume_scaling <= 0) {
         log_error(dev, "Invalid volume scale factor `%g'.\n", volume_scaling);
         res = RES_BAD_ARG;
         goto error;
       }
 
       /* Define the scale factor to apply the the area and the distance from the
        * scale factor of the volume */
       area_scaling = pow(volume_scaling, 2.0/3.0);
       dist_scaling = pow(volume_scaling, 1.0/3.0);
     }
 
     /* Sample an orientation */
     ssp_ran_sphere_uniform_float(ctx->rng, dir, NULL);
 
     /* Build the transformation matrix from world to footprint space. Use the
      * AABB centroid as the origin of the footprint space. */
     build_basis(dir, basis);
     build_transform(centroid, basis, transform);
 
     /* Transform the AABB from world to footprint space */
     FOR_EACH(i, 0, 8) {
       f3_add(pt[i], f33_mulf3(pt[i], transform, aabb_pt[i]), transform + 9);
     }
 
     /* Compute the AABB of the footprint */
     f3_splat(lower, FLT_MAX); f3_splat(upper,-FLT_MAX);
     FOR_EACH(i, 0, 8) {
       f3_min(lower, lower, pt[i]);
       f3_max(upper, upper, pt[i]);
     }
 
     res = accum_monte_carlo_weight
       (dev, ctx, basis, centroid, lower, upper, area_scaling, dist_scaling);
     if(res != RES_OK) goto error;
   }
   /* Compute the Monte Carlo weight of the temporary accumulator */
   FOR_EACH(iwlen, 0, ctx->nwavelengths) {
     size_t iangle;
 
     #define MC_ACCUM(Data) {                                                   \
       const double w = ctx->accums[iwlen].Data / (double)ndirs;                \
       ctx->mc_accums[iwlen].Data.weight += w;                                  \
       ctx->mc_accums[iwlen].Data.sqr_weight += w*w;                            \
     } (void)0
 
     MC_ACCUM(extinction_cross_section);
     MC_ACCUM(absorption_cross_section);
     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]);
       MC_ACCUM(differential_cross_section_cumulative[iangle]);
     }
 
     #undef MC_ACCUM
   }
 exit:
   return res;
 error:
   FOR_EACH(iwlen, 0, ctx->nwavelengths) {
     ctx->mc_accums[iwlen].extinction_cross_section = MC_DATA_NULL;
     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->estimator->limit_angles[iwlen]);
     memset(ctx->mc_accums[iwlen].differential_cross_section_cumulative, 0,
       sizeof(double)*ctx->estimator->limit_angles[iwlen]);
   }
   goto exit;
 }
 
 struct function_arg {
   double scattering_cross_section;
   double angle;
   double differential_cross_section;
   double differential_cross_section_cumulative;
 
   /* Precomputed values */
   double cos_2_angle;
   double cos_angle;
   double cos_half_angle;
   double sin_half_angle;
   double sqr_cos_angle;
   double sqr_sin_half_angle;
 };
 
 static int
 function_fdf(const gsl_vector* vec, void* params, gsl_vector* res, gsl_matrix* J)
 {
   struct function_arg* arg = params;
   const double n = gsl_vector_get(vec, 0);
   double u, u_prime;
   double v, v_prime;
   double coef;
   double f, df;
 
   if(n==2 || n==4 || n==6)
     return GSL_EDOM;
 
   coef = (4*PI * arg->differential_cross_section) / (1+arg->sqr_cos_angle);
 
   u = ( arg->sqr_sin_half_angle
        * ( (n*n - 6*n + 8) * arg->cos_2_angle
          + 3*n*n
          - 8*(n - 2) * arg->cos_angle
          - 26*n + 72 )
        - pow(arg->sin_half_angle, n) * (4*n*n - 24*n + 64) );
   v = ((n - 2) * (n - 4) * (n - 6));
 
   u_prime = arg->sqr_sin_half_angle
           * ((2*n - 6) * arg->cos_2_angle + 6*n - 8*arg->cos_angle - 26)
           - pow(arg->sin_half_angle, n)
           * (log(arg->sin_half_angle) * (4*n*n - 24*n + 64) + 8*n-24);
   v_prime = (n-4)*(n-6) + (n-2)*(n-6) + (n-2)*(n-4);
 
   f  = arg->differential_cross_section_cumulative
      - arg->scattering_cross_section
      + coef * u / v;
   df = (v*u_prime - u*v_prime) / (v*v);
   df = df * coef;
 
   gsl_vector_set(res, 0, f);
   gsl_matrix_set(J, 0, 0, df);
   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
 solve_wide_angles_phase_function_model_parameters
   (struct sschiff_device* dev,
    struct solver_context* ctx,
    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_function_fdf func;
   gsl_vector* root = NULL;
   struct function_arg arg;
   int status;
   int iter;
   const int max_iter = 1000; /* Maximum number of Newton iterations */
   res_T res = RES_OK;
   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));
 
     /* Fill system arguments */
   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_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;
 
   /* Setup system functions */
   func.f = NULL;
   func.df = NULL;
   func.fdf = function_fdf; /* Compute the function *and* its derivative */
   func.n = 1; /* Number of equations */
   func.params = &arg; /* System parameters */
 
   /* Initialise the solver. TODO perform this in an initialisation step of a
    * phase_function context */
   gsl_multiroot_fdfsolver_set(ctx->solver, &func, ctx->init_val);
   gsl_set_error_handler_off();
 
   /* Launch the Newton estimation. Iterate up to `max_iter'or until the
    * estimation is "good enough". */
   iter = 0;
   do {
     status = gsl_multiroot_fdfsolver_iterate(ctx->solver);
     if(status == GSL_SUCCESS) {
       status = gsl_multiroot_test_residual(ctx->solver->f, 1.e-6);
     }
   } while(status == GSL_CONTINUE && ++iter < max_iter);
 
   /* Retrieve the estimated "n" value */
   root = gsl_multiroot_fdfsolver_root(ctx->solver);
   *n = gsl_vector_get(root, 0);
 
   if(status != GSL_SUCCESS) {
     log_error(dev,
 "Cannot estimate the parameters of the wide scattering phase function model.\n"
 "GSL error: %s - %d.\n",
       gsl_strerror(status), iter);
     res = RES_BAD_OP;
     goto error;
   }
 
   if(*n < 0) {
     log_error(dev,
 "The estimated parameter `n' of the wide scattering angles phase function\n"
 "model cannot be negative.\n"
 "n = %g\n", *n);
     res = RES_BAD_OP;
     goto error;
   }
 
   /* Compute r from the estimated n */
   *r = 2 * arg.differential_cross_section * pow(arg.sin_half_angle, *n)
     / (1 + arg.sqr_cos_angle);
 
 exit:
   return res;
 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_small_scattering_angles_properties
   (struct sschiff_device* dev,
    struct solver_context* ctx,
    const size_t iwlen,
    const size_t nangles,
    struct sschiff_estimator* estimator)
 {
   double rcp_scattering_cross_section;
   double rcp_sqr_scattering_cross_section;
   struct sschiff_state scattering_cross_section;
   size_t ilimit_angle;
   size_t iangle;
   ASSERT(dev && estimator && iwlen < sa_size(estimator->wavelengths) && ctx);
   (void)ctx, (void)dev, (void)nangles;
 
   /* Fetch the limit angle and precompute some values */
   ilimit_angle = estimator->limit_angles[iwlen]-1;
 
   get_mc_value
     (&estimator->accums[iwlen].scattering_cross_section,
      estimator->nrealisations,
      &scattering_cross_section);
 
   /* Compute the [cumulative] phase function for small angles */
   rcp_scattering_cross_section = 1.0/scattering_cross_section.E;
   rcp_sqr_scattering_cross_section =
     rcp_scattering_cross_section * rcp_scattering_cross_section;
 
   FOR_EACH(iangle, 0, ilimit_angle+1) {
     struct mc_data mc_data;
 
     mc_data = estimator->accums[iwlen].differential_cross_section[iangle];
     mc_data.weight *= rcp_scattering_cross_section;
     mc_data.sqr_weight *= rcp_sqr_scattering_cross_section;
     get_mc_value(&mc_data, estimator->nrealisations,
        &estimator->phase_functions[iwlen].values[iangle]);
 
     mc_data = estimator->accums[iwlen].differential_cross_section_cumulative[iangle];
     mc_data.weight *= rcp_scattering_cross_section;
     mc_data.sqr_weight *= rcp_sqr_scattering_cross_section;
     get_mc_value(&mc_data, estimator->nrealisations,
        &estimator->phase_functions[iwlen].cumulative[iangle]);
   }
   return RES_OK;
 }
 
 static res_T
 compute_large_scattering_angles_properties
   (struct sschiff_device* dev,
    struct solver_context* ctx,
    const size_t iwlen,
    const size_t nangles,
    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 rcp_scattering_cross_section;
   struct sschiff_state scattering_cross_section;
   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;
   ASSERT(dev && estimator && iwlen < sa_size(estimator->wavelengths) && ctx);
 
   /* Fetch the limit angle and precompute some values */
   ilimit_angle = estimator->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);
 
   get_mc_value
     (&estimator->accums[iwlen].scattering_cross_section,
      estimator->nrealisations,
      &scattering_cross_section);
   get_mc_value
     (&estimator->accums[iwlen].differential_cross_section[ilimit_angle],
      estimator->nrealisations,
      &limit_differential_cross_section);
   get_mc_value
     (&estimator->accums[iwlen].differential_cross_section_cumulative[ilimit_angle],
      estimator->nrealisations,
      &limit_differential_cross_section_cumulative);
 
   /* Find the connector values between small and wide angles */
   res = solve_wide_angles_phase_function_model_parameters
     (dev,
      ctx,
      limit_angle,
      sin_half_limit_angle,
      cos_limit_angle,
      cos_2_limit_angle,
      scattering_cross_section.E,
      limit_differential_cross_section.E,
      limit_differential_cross_section_cumulative.E,
      &n, &r);
   if(res != RES_OK) {
     log_error(estimator->dev,
 "Couldn't estimate the parameters of the phase function for wide scattering\n"
 "angles. The phase function is thus not computed.\n"
 "Wavelength = %g micron\n",
       estimator->wavelengths[iwlen]);
     goto error;
   }
 
   if(n < 2 || n > 4) {
     log_warning(estimator->dev,
 "The wide scattering angles phase function model parameter `n' is not in\n"
 "[2, 4]. One might increase the number of realisations and/or adjust the\n"
 "characteristic length of the particles (that sets the limit between small\n"
 "and large scattering angles) and/or increase the number of scattering\n"
 "angles computed.\n"
 "  wavelength = %g micron\n"
 "  n = %g\n"
 "  scattering cross section = %g +/- %g\n"
 "  differential cross section = %g +/- %g\n"
 "  differential cross section cumulative = %g +/- %g\n",
       estimator->wavelengths[iwlen],
       n,
       scattering_cross_section.E,
       scattering_cross_section.SE,
       limit_differential_cross_section.E,
       limit_differential_cross_section.SE,
       limit_differential_cross_section_cumulative.E,
       limit_differential_cross_section_cumulative.SE);
   }
 
   /* Save the `n' wide angle parameter */
   estimator->wide_angles_parameter[iwlen] = n;
 
   /* 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.E + 2*PI*r*(coef-coef_limit);
   }
 
   /* Compute the [cumulative] phase function for wide angles */
   rcp_scattering_cross_section = 1.0/scattering_cross_section.E;
   FOR_EACH(iangle, ilimit_angle + 1, nangles) {
     estimator->phase_functions[iwlen].values[iangle].E =
       estimator->accums[iwlen].differential_cross_section[iangle].weight
     * rcp_scattering_cross_section;
 
     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 = 0;
     estimator->phase_functions[iwlen].values[iangle].SE = 0;
     estimator->phase_functions[iwlen].cumulative[iangle].V = 0;
     estimator->phase_functions[iwlen].cumulative[iangle].SE = 0;
   }
 
 exit:
   return res;
 error:
   FOR_EACH(iangle, ilimit_angle+1, nangles) {
     estimator->phase_functions[iwlen].values[iangle].E = -1;
     estimator->phase_functions[iwlen].values[iangle].V = -1;
     estimator->phase_functions[iwlen].values[iangle].SE = -1;
     estimator->phase_functions[iwlen].cumulative[iangle].E = -1;
     estimator->phase_functions[iwlen].cumulative[iangle].V = -1;
     estimator->phase_functions[iwlen].cumulative[iangle].SE = -1;
   }
   goto exit;
 }
 
 static res_T
 compute_scattering_angles_properties
   (struct sschiff_device* dev,
    struct solver_context* ctx,
    const size_t iwlen,
    const size_t nangles,
    struct sschiff_estimator* estimator)
 {
   res_T res = RES_OK;
 
   res = compute_small_scattering_angles_properties
     (dev, ctx, iwlen, nangles, estimator);
   if(res != RES_OK) goto error;
 
   if(!estimator->discard_large_angles) {
     res = compute_large_scattering_angles_properties
       (dev, ctx, iwlen, nangles, estimator);
     if(res != RES_OK) goto error;
   }
 
 exit:
   return res;
 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->limit_angles[iwlen];
   if(cumulative < cumulative_small_angles[nsmall_angles-1]) {
     /* Look for 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) {
       log_error(estimator->dev,
         "Error in inverting the phase function cumulative for small angles.\n");
       return RES_BAD_OP;
     }
     iangle = (size_t)(found - cumulative_small_angles);
     upper = cumulative_small_angles[iangle];
     lower = cumulative_small_angles[iangle-1];
 
     /* 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];
   } else if(estimator->discard_large_angles) {
     *theta = -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) {
       log_error(estimator->dev,
         "Error in inverting the phase function cumulative for wide angles.\n");
       return RES_BAD_OP;
     }
     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)
 {
   ASSERT(distrib);
   return distrib->sample != NULL
       && distrib->characteristic_length > 0;
 }
 
 static void
 integrator_context_release(struct integrator_context* ctx)
 {
   size_t iwlen;
   ASSERT(ctx);
 
   if(ctx->rng) SSP(rng_ref_put(ctx->rng));
   if(ctx->estimator) SSCHIFF(estimator_ref_put(ctx->estimator));
   if(ctx->shape) S3D(shape_ref_put(ctx->shape));
   if(ctx->scene) S3D(scene_ref_put(ctx->scene));
   if(ctx->view) S3D(scene_view_ref_put(ctx->view));
 
   #define RELEASE(Data) if(Data) sa_release(Data)
   if(ctx->accums) {
     FOR_EACH(iwlen, 0, ctx->nwavelengths) {
       RELEASE(ctx->accums[iwlen].differential_cross_section);
       RELEASE(ctx->accums[iwlen].differential_cross_section_cumulative);
     }
     sa_release(ctx->accums);
   }
   if(ctx->mc_accums) {
     FOR_EACH(iwlen, 0, ctx->nwavelengths) {
       RELEASE(ctx->mc_accums[iwlen].differential_cross_section);
       RELEASE(ctx->mc_accums[iwlen].differential_cross_section_cumulative);
     }
     sa_release(ctx->mc_accums);
   }
   #undef RELEASE
 }
 
 static res_T
 integrator_context_init
   (struct integrator_context* ctx,
    struct s3d_device* s3d,
    struct ssp_rng* rng,
    struct sschiff_estimator* estimator,
    const size_t nwavelengths,
    const size_t nangles)
 {
   size_t iwlen;
   res_T res = RES_OK;
   ASSERT(ctx && s3d && rng && nwavelengths && nangles >= 2);
 
   memset(ctx, 0, sizeof(struct integrator_context));
 
   if(RES_OK!=(res = s3d_shape_create_mesh(s3d, &ctx->shape))) goto error;
   if(RES_OK!=(res = s3d_scene_create(s3d, &ctx->scene))) goto error;
   if(RES_OK!=(res = s3d_scene_attach_shape(ctx->scene, ctx->shape))) goto error;
 
   #define RESIZE(Data, N) {                                                    \
     if(!sa_add((Data), (N))) {                                                 \
       res = RES_MEM_ERR;                                                       \
       goto error;                                                              \
     }                                                                          \
     memset((Data), 0, sizeof((Data)[0])*(N));                                  \
   } (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);
   }
   #undef RESIZE
 
   SSP(rng_ref_get(rng));
   ctx->rng = rng;
   SSCHIFF(estimator_ref_get(estimator));
   ctx->estimator = estimator;
   ctx->nwavelengths = nwavelengths;
   ctx->nangles = nangles;
 
 exit:
   return res;
 error:
   integrator_context_release(ctx);
   goto exit;
 }
 
 static void
 solver_context_release(struct solver_context* ctx)
 {
   ASSERT(ctx);
   if(ctx->solver) gsl_multiroot_fdfsolver_free(ctx->solver);
   if(ctx->init_val) gsl_vector_free(ctx->init_val);
 }
 
 static res_T
 solver_context_init(struct sschiff_device* dev, struct solver_context* ctx)
 {
   res_T res = RES_OK;
   ASSERT(ctx);
   memset(ctx, 0, sizeof(*ctx));
 
   /* Create a Newton-Raphson nonlinear solver with derivatives */
   ctx->solver = gsl_multiroot_fdfsolver_alloc(gsl_multiroot_fdfsolver_newton, 1);
   if(!ctx->solver) {
     log_error(dev,
       "Not enough memory. Couldn't allocate the GSL Newton solver.\n");
     res = RES_MEM_ERR;
     goto error;
   }
 
   /* Allocate the initial value vector */
   ctx->init_val = gsl_vector_alloc(1);
   if(!ctx->init_val) {
     log_error(dev,
       "Not enough memory. Couldn't allocate the GSL init vector.\n");
     res = RES_MEM_ERR;
     goto error;
   }
   gsl_vector_set(ctx->init_val, 0, 3);
 
 exit:
   return res;
 error:
   solver_context_release(ctx);
   goto exit;
 }
 
 static res_T
 begin_realisation
   (struct sschiff_device* dev,
    struct sschiff_geometry_distribution* distrib,
    struct integrator_context* ctx)
 {
   struct s3d_accel_struct_conf accel_struct_conf = S3D_ACCEL_STRUCT_CONF_DEFAULT;
   res_T res = RES_OK;
   ASSERT(dev && distrib && ctx && !ctx->view);
 
   /* Sample a particle */
   res = distrib->sample
     (ctx->rng, &ctx->shape_data, ctx->shape, distrib->context);
   if(res != RES_OK) {
     log_error(dev, "Error in sampling a particle.\n");
     goto error;
   }
 
   accel_struct_conf.quality = S3D_ACCEL_STRUCT_QUALITY_LOW;
   accel_struct_conf.mask =
     S3D_ACCEL_STRUCT_FLAG_ROBUST
   | S3D_ACCEL_STRUCT_FLAG_DYNAMIC;
 
   /* Build the Star-3D representation of the sampled shape */
   S3D(scene_view_create2(ctx->scene, S3D_TRACE, &accel_struct_conf, &ctx->view));
 
 exit:
   return res;
 error:
   if(ctx->view) {
     S3D(scene_view_ref_put(ctx->view));
     ctx->view = NULL;
   }
   goto exit;
 }
 
 static FINLINE void
 end_realisation(struct integrator_context* ctx)
 {
   ASSERT(ctx && ctx->view);
   S3D(scene_view_ref_put(ctx->view));
   ctx->view = NULL;
 }
 
 static void
 estimator_release(ref_T* ref)
 {
   struct sschiff_estimator* estimator;
   struct sschiff_device* dev;
   size_t nwavelengths;
   size_t i;
   ASSERT(ref);
 
   estimator = CONTAINER_OF(ref, struct sschiff_estimator, ref);
   dev = estimator->dev;
   nwavelengths = sa_size(estimator->accums);
 
   #define RELEASE(Data) if(Data) sa_release(Data)
   RELEASE(estimator->wavelengths);
   RELEASE(estimator->angles);
   RELEASE(estimator->limit_angles);
   RELEASE(estimator->wide_angles_parameter);
   RELEASE(estimator->wavenumbers);
   RELEASE(estimator->properties);
   if(estimator->accums) {
     FOR_EACH(i, 0, nwavelengths) {
       RELEASE(estimator->accums[i].differential_cross_section);
       RELEASE(estimator->accums[i].differential_cross_section_cumulative);
     }
     sa_release(estimator->accums);
   }
   if(estimator->phase_functions) {
     FOR_EACH(i, 0, nwavelengths) {
       RELEASE(estimator->phase_functions[i].values);
       RELEASE(estimator->phase_functions[i].cumulative);
     }
     sa_release(estimator->phase_functions);
   }
   #undef RELEASE
 
   MEM_RM(dev->allocator, estimator);
   SSCHIFF(device_ref_put(dev));
 }
 
 static res_T
 estimator_create
   (struct sschiff_device* dev,
    struct sschiff_geometry_distribution* distrib,
    const double* wavelengths,
    const size_t nwavelengths,
    sschiff_scattering_angles_distribution_T angles_distrib,
    const size_t nangles,
    const int discard_large_angles,
    struct sschiff_estimator** out_estimator)
 {
   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);
 
   estimator = MEM_CALLOC(dev->allocator, 1, sizeof(struct sschiff_estimator));
   if(!estimator) {
     log_error
       (dev, "Couldn't allocate the Star-Schiff estimator.\n");
     res = RES_MEM_ERR;
     goto error;
   }
   ref_init(&estimator->ref);
   SSCHIFF(device_ref_get(dev));
   estimator->dev = dev;
   estimator->discard_large_angles = discard_large_angles;
 
   #define RESIZE(Array, Count, ErrMsg) {                                       \
     if(!sa_add(Array, Count)) {                                                \
       log_error(dev, ErrMsg);                                                  \
       res = RES_MEM_ERR;                                                       \
       goto error;                                                              \
     }                                                                          \
     memset(Array, 0, sizeof(Array[0])*Count);                                  \
   } (void)0
 
   /* Check the wavelengths */
   FOR_EACH(i, 1, nwavelengths) {
     if(wavelengths[i] <= wavelengths[i-1]) {
       log_error(dev,
         "The submitted wavelengths are not sorted in ascending order.\n");
       res = RES_BAD_ARG;
       goto error;
     }
   }
   /* Copy the wavelengths */
   RESIZE(estimator->wavelengths, nwavelengths,
     "Couldn't allocate the list of estimated wavelengths.\n");
   memcpy(estimator->wavelengths, wavelengths, nwavelengths*sizeof(double));
 
   /* Generate the scattering angles */
   RESIZE(estimator->angles, nangles,
     "Couldn't allocate the list of scattering angles.\n");
   angles_distrib(estimator->angles, nangles);
   if(estimator->angles[0] != 0.f
   || !eq_eps(estimator->angles[nangles-1], PI, 1.e-6)) {
     log_error(dev, "Invalid scattering angle distribution.\n");
     log_error(dev, "The first and last angles must be 0 and PI, respectively.\n");
     res = RES_BAD_ARG;
     goto error;
   }
   FOR_EACH(i, 1, nangles) {
     if(estimator->angles[i] <= estimator->angles[i-1]) {
       log_error(dev, "Invalid scattering angle distribution.\n");
       log_error(dev, "Angles must be sorted in ascending order in [0, PI]\n");
       res = RES_BAD_ARG;
       goto error;
     }
   }
 
   /* Allocate miscellaneous array of temporary/precomputed data */
   RESIZE(estimator->limit_angles, nwavelengths,
     "Couldn't allocate the per wavelength indices of the limit scattering angle.\n");
   RESIZE(estimator->wide_angles_parameter, nwavelengths,
     "Couldn't allocate the per wavelength `n' parameter of the wide scattering\n"
     "angles model.\n");
   RESIZE(estimator->wavenumbers, nwavelengths,
     "Couldn't allocate the list of wavenumbers.\n");
   RESIZE(estimator->properties, nwavelengths,
     "Couldn't allocate the per wavelength optical properties.\n");
 
   rcp_PI_Lc = 1.0 / (PI*distrib->characteristic_length);
   FOR_EACH(i, 0, nwavelengths) {
     double lambda_e, theta_l;
     double* angle;
 
     /* Fetch the particle optical properties */
     distrib->material.get_property
       (distrib->material.material,
        estimator->wavelengths[i],
        &estimator->properties[i]);
 
     /* Precompute the wavenumbers */
     lambda_e =
       estimator->wavelengths[i]
     / estimator->properties[i].medium_refractive_index;
     estimator->wavenumbers[i] = 2.0*PI/lambda_e;
 
     /* Search for limit scattering angle */
     theta_l = sqrt(lambda_e * rcp_PI_Lc);
     if(theta_l <= 0 || theta_l >= PI) {
       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]);
       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;
       }*/
       /* 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;
     }
   }
   estimator->no_phase_function = !hold_valid_limit_angle;
   if(estimator->no_phase_function) {
     /* No phase function => scattering angles and wide angle parameter are
      * useless */
     sa_release(estimator->angles);
     sa_release(estimator->wide_angles_parameter);
     estimator->angles = NULL;
     estimator->wide_angles_parameter = NULL;
   }
 
   /* Allocate the estimator accumulators */
   RESIZE(estimator->accums, nwavelengths,
     "Couldn't allocate the Monte Carlo accumulator of the estimator.\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");
     }
   }
 
   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;
 error:
   if(estimator) {
     SSCHIFF(estimator_ref_put(estimator));
     estimator = NULL;
   }
   goto exit;
 }
 
 /*******************************************************************************
  * Exported functions
  ******************************************************************************/
 res_T
 sschiff_integrate
   (struct sschiff_device* dev,
    struct ssp_rng* rng,
    struct sschiff_geometry_distribution* distrib,
    const double* wavelengths,
    const size_t nwavelengths,
    const sschiff_scattering_angles_distribution_T angles_distrib,
    const size_t nangles,
    const size_t ngeoms,
    const size_t ndirs,
    const int discard_large_angles,
    struct sschiff_estimator** out_estimator)
 {
   struct integrator_context* ctxs = NULL;
   struct solver_context* solver_ctxs = NULL;
   struct sschiff_estimator* estimator = NULL;
   struct ssp_rng** rngs = NULL;
   struct ssp_rng_proxy* rng_proxy = NULL;
   size_t i;
   int iwlen;
   int igeom;
   ATOMIC res = (ATOMIC)RES_OK;
 
   if(!dev || !rng || !distrib || !wavelengths || !nwavelengths
   || !angles_distrib || nangles < 2 || !ngeoms || !ndirs || !out_estimator
   || !check_distribution(distrib)) {
     return RES_BAD_ARG;
   }
 
   /* Create the Schiff estimator */
   res = estimator_create(dev, distrib, wavelengths, nwavelengths,
     angles_distrib, nangles, discard_large_angles, &estimator);
   if(res != RES_OK) goto error;
   estimator->nrealisations = ngeoms;
 
   /* Create a RNG proxy from the submitted RNG state */
   res = ssp_rng_proxy_create_from_rng
     (dev->allocator, rng, dev->nthreads, &rng_proxy);
   if(res != RES_OK) {
     log_error(dev, "Couldn't create the proxy allocator\n");
     goto error;
   }
 
   /* Create per thread data structures */
   if(!sa_add(ctxs, dev->nthreads)) {
     log_error(dev, "Couldn't allocate the integrator contexts.\n");
     res = RES_MEM_ERR;
     goto error;
   }
   if(!sa_add(solver_ctxs, dev->nthreads)) {
     log_error(dev, "Couldn'nt allocate the solver contexts.\n");
     res = RES_MEM_ERR;
     goto error;
   }
   if(!sa_add(rngs, dev->nthreads)) {
     log_error(dev, "Couldn't allocate the list of RNGs.\n");
     res = RES_MEM_ERR;
     goto error;
   }
   memset(ctxs, 0, sizeof(ctxs[0])*dev->nthreads);
   memset(solver_ctxs, 0, sizeof(solver_ctxs[0])*dev->nthreads);
   memset(rngs, 0, sizeof(rngs[0])*dev->nthreads);
 
   /* Init the per thread data structures */
   FOR_EACH(i, 0, dev->nthreads) {
     res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs + i);
     if(res != RES_OK) {
       log_error(dev,
         "Couldn't create the RNG of the thread \"%lu\".\n", (unsigned long)i);
       goto error;
     }
     res = integrator_context_init(ctxs + i, dev->s3d[i], rngs[i], estimator,
       nwavelengths, nangles);
     if(res != RES_OK) {
       log_error(dev,
         "Couldn't initialise the integrator context of the thread \"%lu\".\n",
         (unsigned long)i);
       goto error;
     }
     res = solver_context_init(dev, solver_ctxs + i);
     if(res != RES_OK) {
       log_error(dev,
         "Couldn't initialise the solver context of the thread \"%lu\".\n",
         (unsigned long)i);
       goto error;
     }
   }
 
   /* Paralell Schiff integration */
   #pragma omp parallel for schedule(static)
   for(igeom=0; igeom < (int)ngeoms; ++igeom) {
     const int ithread = omp_get_thread_num();
     ATOMIC res_local = RES_OK;
 
     if(ATOMIC_GET(&res) != RES_OK) continue;
 
     /* Setup the data for the current realisation */
     res_local = begin_realisation(dev, distrib, ctxs+ithread);
     if(res_local != RES_OK) {
       ATOMIC_CAS(&res, res_local, RES_OK);
       continue;
     }
 
     /* Schiff Estimation */
     res_local = radiative_properties(dev, igeom, ndirs, distrib, ctxs+ithread);
     if(res_local != RES_OK) ATOMIC_CAS(&res, res_local, RES_OK);
 
     end_realisation(ctxs + ithread);
   }
   if(res != RES_OK) goto error; /* Handle integration error */
 
   /* Merge the per thread integration results */
   FOR_EACH(i, 0, dev->nthreads) {
     FOR_EACH(iwlen, 0, (int)nwavelengths) {
       size_t iangle;
       #define MC_ACCUM(Data) {                                                 \
         const struct mc_accum* mc_accum = ctxs[i].mc_accums + iwlen;           \
         estimator->accums[iwlen].Data.weight += mc_accum->Data.weight;         \
         estimator->accums[iwlen].Data.sqr_weight += mc_accum->Data.sqr_weight; \
       } (void)0
       MC_ACCUM(extinction_cross_section);
       MC_ACCUM(absorption_cross_section);
       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]);
         MC_ACCUM(differential_cross_section_cumulative[iangle]);
       }
     }
   }
 
   /* Static scheduling is not necessary to ensure the reproductability;
    * computations are purely analytics */
   #pragma omp parallel for
   for(iwlen=0; iwlen < (int)nwavelengths; ++iwlen) {
     const int ithread = omp_get_thread_num();
     if(estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE) continue;
     /* Do not handle phase function errors */
     compute_scattering_angles_properties
       (dev, &solver_ctxs[ithread], (size_t)iwlen, nangles, estimator);
   }
 
 exit:
   if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
   if(ctxs) {
     FOR_EACH(i, 0, dev->nthreads) integrator_context_release(ctxs+i);
     sa_release(ctxs);
   }
   if(solver_ctxs) {
     FOR_EACH(i, 0, dev->nthreads) solver_context_release(solver_ctxs+i);
     sa_release(solver_ctxs);
   }
   if(rngs) {
     FOR_EACH(i, 0, dev->nthreads) if(rngs[i]) SSP(rng_ref_put(rngs[i]));
     sa_release(rngs);
   }
   if(out_estimator) *out_estimator = estimator;
   return (res_T)res;
 error:
   if(estimator) {
     SSCHIFF(estimator_ref_put(estimator));
     estimator = NULL;
   }
   goto exit;
 }
 
 res_T
 sschiff_estimator_ref_get(struct sschiff_estimator* estimator)
 {
   if(!estimator) return RES_BAD_ARG;
   ref_get(&estimator->ref);
   return RES_OK;
 }
 
 res_T
 sschiff_estimator_ref_put(struct sschiff_estimator* estimator)
 {
   if(!estimator) return RES_BAD_ARG;
   ref_put(&estimator->ref, estimator_release);
   return RES_OK;
 }
 
 res_T
 sschiff_estimator_get_wavelengths
   (const struct sschiff_estimator* estimator,
    const double** wavelengths,
    size_t* count)
 {
   if(!estimator) return RES_BAD_ARG;
   if(wavelengths) *wavelengths = estimator->wavelengths;
   if(count) *count = sa_size(estimator->wavelengths);
   return RES_OK;
 }
 
 res_T
 sschiff_estimator_get_scattering_angles
   (const struct sschiff_estimator* estimator,
    const double** angles,
    size_t* count)
 {
   if(!estimator) return RES_BAD_ARG;
   if(angles) *angles = estimator->angles;
   if(count) *count = sa_size(estimator->angles);
   return RES_OK;
 }
 
 res_T
 sschiff_estimator_get_realisations_count
   (const struct sschiff_estimator* estimator, size_t* count)
 {
   if(!estimator || !count) return RES_BAD_ARG;
   *count = estimator->nrealisations;
   return RES_OK;
 }
 
 res_T
 sschiff_estimator_get_cross_section
   (const struct sschiff_estimator* estimator,
    const size_t iwlen,
    struct sschiff_cross_section* cross_section)
 {
   const struct mc_accum* acc;
   size_t N;
 
   if(!estimator || !cross_section || iwlen >= sa_size(estimator->wavelengths))
     return RES_BAD_ARG;
 
   acc = estimator->accums + iwlen;
   N = estimator->nrealisations;
   get_mc_value(&acc->extinction_cross_section, N, &cross_section->extinction);
   get_mc_value(&acc->absorption_cross_section, N, &cross_section->absorption);
   get_mc_value(&acc->scattering_cross_section, N, &cross_section->scattering);
   get_mc_value(&acc->avg_projected_area, N,
     &cross_section->average_projected_area);
   return RES_OK;
 }
 
 res_T
 sschiff_estimator_get_cross_sections
   (const struct sschiff_estimator* estimator,
    struct sschiff_cross_section* cross_sections)
 {
   size_t nwavelengths;
   size_t iwlen;
   size_t N;
   if(!estimator || !cross_sections) return RES_BAD_ARG;
 
   nwavelengths = sa_size(estimator->wavelengths);
   N = estimator->nrealisations;
   FOR_EACH(iwlen, 0, nwavelengths) {
     const struct mc_accum* acc = estimator->accums + iwlen;
     get_mc_value(&acc->extinction_cross_section, N, &cross_sections->extinction);
     get_mc_value(&acc->absorption_cross_section, N, &cross_sections->absorption);
     get_mc_value(&acc->scattering_cross_section, N, &cross_sections->scattering);
     get_mc_value(&acc->avg_projected_area, N,
       &cross_sections->average_projected_area);
   }
   return RES_OK;
 }
 
 res_T
 sschiff_estimator_get_phase_function
   (const struct sschiff_estimator* estimator,
    const size_t iwlen,
    const struct sschiff_state** states)
 {
   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_BAD_OP;
   /* Invalid phase function */
   if(estimator->phase_functions[iwlen].values[0].V < 0)
     return RES_BAD_OP;
   *states = estimator->phase_functions[iwlen].values;
   return RES_OK;
 }
 
 /* Retrieve a pointer onto the estimated phase function cumulative. */
 SSCHIFF_API res_T
 sschiff_estimator_get_phase_function_cumulative
   (const struct sschiff_estimator* estimator,
    const size_t iwlen,
    const struct sschiff_state** states)
 {
   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_BAD_OP;
   /* Invalid cumulative phase function */
   if(estimator->phase_functions[iwlen].cumulative[0].V < 0)
     return RES_BAD_OP;
   *states = estimator->phase_functions[iwlen].cumulative;
   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
   || iwlen >= sa_size(estimator->wavelengths)) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   /* No phase function cumulative is computed => nothing to inverse */
   if(estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE) {
     res = RES_BAD_OP;
     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->limit_angles[iwlen];
   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[iangle];
 
     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-1);
   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 */
     res = inverse_cumulative_phase_function
       (estimator, cumulative_small_angles, iwlen, cumulative, &thetas[itheta]);
     if(res != RES_OK) goto error;
     cumulative += step;
   }
 
 exit:
   if(cumulative_small_angles) sa_release(cumulative_small_angles);
   return res;
 error:
   goto exit;
 }
 
 res_T
 sschiff_estimator_get_limit_scattering_angle_index
   (const struct sschiff_estimator* estimator,
    const size_t iwlen,
    size_t* iangle)
 {
   size_t limit_angle;
   if(!estimator || iwlen >= sa_size(estimator->wavelengths) || !iangle)
      return RES_BAD_ARG;
   limit_angle = estimator->limit_angles[iwlen];
   if(limit_angle == INVALID_LIMIT_ANGLE)
     return RES_BAD_OP;
   *iangle = limit_angle - 1/*<=>index of the last small angle*/;
   return RES_OK;
 }
 
 res_T
 sschiff_estimator_get_wide_scattering_angle_model_parameter
   (const struct sschiff_estimator* estimator,
    const size_t iwlen,
    double* out_n)
 {
   if(!estimator || iwlen >= sa_size(estimator->wavelengths) || !out_n)
     return RES_BAD_ARG;
   if(estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE)
     return RES_BAD_OP;
   *out_n = estimator->wide_angles_parameter[iwlen];
   return RES_OK;
 }
 
 res_T
 sschiff_estimator_get_differential_cross_section
   (const struct sschiff_estimator* estimator,
    const size_t iwlen,
    const size_t iangle,
    struct sschiff_state* state)
 {
   if(!estimator
   || iwlen >= sa_size(estimator->wavelengths)
   || iangle >= sa_size(estimator->angles)
   || !state)
     return RES_BAD_ARG;
 
   if(estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE)
     return RES_BAD_OP;
 
   get_mc_value
     (&estimator->accums[iwlen].differential_cross_section[iangle],
      estimator->nrealisations,
      state);
   return RES_OK;
 }
 
 res_T
 sschiff_estimator_get_differential_cross_section_cumulative
   (const struct sschiff_estimator* estimator,
    const size_t iwlen,
    const size_t iangle,
    struct sschiff_state* state)
 {
   if(!estimator
   || iwlen >= sa_size(estimator->wavelengths)
   || iangle >= sa_size(estimator->angles)
   || !state)
     return RES_BAD_ARG;
 
   if(estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE)
     return RES_BAD_OP;
 
   get_mc_value
     (&estimator->accums[iwlen].differential_cross_section_cumulative[iangle],
      estimator->nrealisations,
      state);
   return RES_OK;
 }