commit a57cf548c99f5ea0af24d4df231b45f3d14d2df9
parent 1a1eb57276d2d57455b34cde392dc6013935f95a
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 13 Oct 2015 09:39:26 +0200
Test the Schiff estimator
Test the estimation of radiative properties on spheres whose the radius is
distributed with respect to a lognormal distribution.
Diffstat:
3 files changed, 316 insertions(+), 109 deletions(-)
diff --git a/src/sschiff.h b/src/sschiff.h
@@ -56,12 +56,10 @@ struct s3d_device;
struct s3d_shape;
struct ssp_rng;
-struct sschiff_refractive_index { double real, imag; };
-
-enum sschiff_material_property {
- SSCHIFF_MATERIAL_MEDIUM_REFRACTIVE_INDEX,
- SSCHIFF_MATERIAL_RELATIVE_REAL_REFRACTIVE_INDEX,
- SSCHIFF_MATERIAL_RELATIVE_IMAGINARY_REFRACTIVE_INDEX
+struct sschiff_material_properties {
+ double medium_refractive_index;
+ double relative_real_refractive_index;
+ double relative_imaginary_refractive_index;
};
struct sschiff_material {
@@ -69,8 +67,7 @@ struct sschiff_material {
void (*get_property)
(void* material,
const double wavelength,
- const enum sschiff_material_property property,
- double* value);
+ struct sschiff_material_properties* properties);
void* material; /* Opaque user defined representation of the material */
};
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -46,20 +46,21 @@
#define NDIRS 1 /* # Sampled directions */
+struct accum {
+ double weight;
+ double sqr_weight;
+ size_t naccums;
+};
+
struct sschiff_estimator {
size_t nsteps;
double wavelength;
+ struct accum extinction_cross_section;
struct sschiff_device* dev;
ref_T ref;
};
-struct accum {
- double weight;
- double sqr_weight;
- size_t naccums;
-};
-
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -160,7 +161,7 @@ compute_path_length
} while(!S3D_HIT_NONE(&hit));
- *length = len * 1.e-6/* Micron to meter */;
+ *length = (len - first_hit->distance) * 1.e-6/* Micron to meter */;
return RES_OK;
}
@@ -190,7 +191,7 @@ compute_monte_carlo_weight
ASSERT(scn && rng && mtl && basis && pos && lower && upper && accum);
f2_sub(plane_size, upper, lower);
- rcp_pdf = plane_size[0] * plane_size[1];
+ rcp_pdf = plane_size[0] * 1.e-6f * plane_size[1] * 1.e-6f;
/* Define the projection axis */
f3_mulf(axis_x, basis + 0, plane_size[0] * 0.5f);
@@ -212,6 +213,7 @@ compute_monte_carlo_weight
if(S3D_HIT_NONE(&hit)) {
weight = 0.0;
} else {
+ struct sschiff_material_properties props;
double n_r, k_r;
double n_e, k_e, lambda_e;
double path_length;
@@ -219,13 +221,11 @@ compute_monte_carlo_weight
res = compute_path_length(scn, &hit, ray_org, axis_z, &path_length);
if(res != RES_OK) goto error;
- /* Retrieve material properties */
- mtl->get_property(mtl->material, wavelength,
- SSCHIFF_MATERIAL_MEDIUM_REFRACTIVE_INDEX, &n_e);
- mtl->get_property(mtl->material, wavelength,
- SSCHIFF_MATERIAL_RELATIVE_REAL_REFRACTIVE_INDEX, &n_r);
- mtl->get_property(mtl->material, wavelength,
- SSCHIFF_MATERIAL_RELATIVE_IMAGINARY_REFRACTIVE_INDEX, &k_r);
+ mtl->get_property(mtl->material, wavelength, &props);
+
+ n_e = props.medium_refractive_index;
+ n_r = props.relative_real_refractive_index;
+ k_r = props.relative_imaginary_refractive_index;
lambda_e = wavelength / n_e;
k_e = 2.0 * PI / lambda_e;
@@ -325,10 +325,10 @@ radiative_properties
struct sschiff_integrator_desc* desc,
const double wavelength,
const int istep,
+ struct s3d_scene* scene,
+ struct s3d_shape* shape,
struct accum* accum)
{
- struct s3d_shape* shape = NULL;
- struct s3d_scene* scene = NULL;
struct sschiff_material material = SSCHIFF_NULL_MATERIAL;
float lower[3], upper[3];
float aabb_pt[8][3];
@@ -337,23 +337,6 @@ radiative_properties
res_T res = RES_OK;
ASSERT(dev && rng && desc && accum);
- res = s3d_shape_create_mesh(dev->s3d, &shape);
- if(res != RES_OK) {
- log_error(dev, "Couldn't create the Star-3D Schiff shape.\n");
- goto error;
- }
- res = s3d_scene_create(dev->s3d, &scene);
- if(res != RES_OK) {
- log_error(dev, "Couldn't create the Star-3D Schiff scene.\n");
- goto error;
- }
- res = s3d_scene_attach_shape(scene, shape);
- if(res != RES_OK) {
- log_error(dev,
- "Couldn't attach the Star-3D Schiff shape to the Star-3D Schiff scene.\n");
- goto error;
- }
-
res = desc->sample_geometry(rng, &material, shape, desc->sampler_context);
if(res != RES_OK) {
log_error(dev, "Error in sampling a Schiff geometry.\n");
@@ -407,7 +390,7 @@ radiative_properties
f3_max(upper, upper, pt[i]);
}
-#if 0
+#if 1
(void)istep;
res = compute_monte_carlo_weight
(scene, rng, &material, wavelength, basis, centroid, lower, upper, accum);
@@ -425,8 +408,6 @@ radiative_properties
S3D(scene_end_session(scene));
exit:
- if(scene) S3D(scene_ref_put(scene));
- if(shape) S3D(shape_ref_put(shape));
return res;
error:
goto exit;
@@ -495,7 +476,8 @@ sschiff_integrate
struct sschiff_estimator** out_estimator)
{
struct sschiff_estimator* estimator = NULL;
- struct accum accum;
+ struct s3d_shape* shape = NULL;
+ struct s3d_scene* scene = NULL;
int i;
res_T res = RES_OK;
@@ -510,14 +492,34 @@ sschiff_integrate
estimator->nsteps = steps_count;
estimator->wavelength = wavelength;
- memset(&accum, 0, sizeof(struct accum));
+ res = s3d_shape_create_mesh(dev->s3d, &shape);
+ if(res != RES_OK) {
+ log_error(dev, "Couldn't create the Star-3D Schiff shape.\n");
+ goto error;
+ }
+ res = s3d_scene_create(dev->s3d, &scene);
+ if(res != RES_OK) {
+ log_error(dev, "Couldn't create the Star-3D Schiff scene.\n");
+ goto error;
+ }
+ res = s3d_scene_attach_shape(scene, shape);
+ if(res != RES_OK) {
+ log_error(dev,
+ "Couldn't attach the Star-3D Schiff shape to the Star-3D Schiff scene.\n");
+ goto error;
+ }
+
+ memset(&estimator->extinction_cross_section, 0, sizeof(struct accum));
FOR_EACH(i, 0, (int)steps_count) {
- res = radiative_properties(dev, rng, desc, wavelength, i, &accum);
+ res = radiative_properties(dev, rng, desc, wavelength, i, scene, shape,
+ &estimator->extinction_cross_section);
if(res != RES_OK) goto error;
}
exit:
if(out_estimator) *out_estimator = estimator;
+ if(shape) S3D(shape_ref_put(shape));
+ if(scene) S3D(scene_ref_put(scene));
return res;
error:
if(estimator) {
@@ -543,3 +545,21 @@ sschiff_estimator_ref_put(struct sschiff_estimator* estimator)
return RES_OK;
}
+res_T
+sschiff_estimator_get_extinction_cross_section
+ (struct sschiff_estimator* estimator,
+ struct sschiff_estimator_status* status)
+{
+ const struct accum* acc;
+ if(!estimator || !status) return RES_BAD_ARG;
+
+ acc = &estimator->extinction_cross_section;
+
+ status->E = acc->weight / (double)acc->naccums;
+ status->V = acc->sqr_weight / (double)acc->naccums
+ - (acc->weight * acc->weight) / (double)(acc->naccums * acc->naccums);
+ status->SE = sqrt(status->V / (double)acc->naccums);
+ status->nsteps = acc->naccums;
+ return RES_OK;
+}
+
diff --git a/src/test_sschiff_estimator.c b/src/test_sschiff_estimator.c
@@ -29,32 +29,65 @@
#include "sschiff.h"
#include "test_sschiff_utils.h"
+#include <rsys/clock_time.h>
#include <rsys/float3.h>
+#include <rsys/stretchy_array.h>
#include <star/s3d.h>
#include <star/ssp.h>
-#define STEP 0.05f
+#define STEP 0.1f
+#define NREALISATIONS 1000
-enum { A, B, M, N0, N1, N2 };
+struct optical_properties { double m, n_r, kappa_r, n_medium; };
-struct super_shape {
- float super_formulas[2][6];
+struct sphere_context {
+ struct optical_properties props;
+ double mean_radius;
+};
+
+struct sphere {
+ double radius;
size_t ntheta, nphi;
};
+struct cylinder {
+ float* vertices;
+ unsigned* indices;
+};
+
static void
-super_shape_get_indices(const unsigned itri, unsigned ids[3], void* ctx)
+get_material_property
+ (void* mtl,
+ const double wavelength,
+ struct sschiff_material_properties* properties)
{
- struct super_shape* sshape = ctx;
+ struct optical_properties* props = mtl;
+
+ NCHECK(mtl, NULL);
+ NCHECK(properties, NULL);
+ CHECK(eq_eps(props->m, wavelength, 1.e-8), 1);
+
+ properties->medium_refractive_index = props->n_medium;
+ properties->relative_real_refractive_index = props->n_r;
+ properties->relative_imaginary_refractive_index = props->kappa_r;
+}
+
+/*******************************************************************************
+ * Sphere geometry
+ ******************************************************************************/
+static void
+sphere_get_indices(const unsigned itri, unsigned ids[3], void* ctx)
+{
+ struct sphere* sphere = ctx;
const size_t iquad = itri / 2;
const size_t iquad_tri = itri % 2;
- const size_t i = iquad / (sshape->nphi - 1);
- const size_t j = iquad % (sshape->nphi - 1);
+ const size_t i = iquad / (sphere->nphi - 1);
+ const size_t j = iquad % (sphere->nphi - 1);
size_t A, B;
- B = i * sshape->nphi + j + 1;
- A = ((i + 1) % sshape->ntheta) * sshape->nphi + j;
+ B = i * sphere->nphi + j + 1;
+ A = ((i + 1) % sphere->ntheta) * sphere->nphi + j;
if(iquad_tri == 0) {
ids[0] = (unsigned)B - 1;
ids[1] = (unsigned)A;
@@ -66,20 +99,17 @@ super_shape_get_indices(const unsigned itri, unsigned ids[3], void* ctx)
}
}
-static void
-super_shape_get_vertex_position(const unsigned ivert, float vertex[3], void* ctx)
+static INLINE void
+sphere_get_vertex_position(const unsigned ivert, float vertex[3], void* ctx)
{
- struct super_shape* sshape = ctx;
- const size_t itheta = ivert / sshape->nphi;
- const size_t iphi = ivert % sshape->nphi;
+ struct sphere* sphere = ctx;
+ const size_t itheta = ivert / sphere->nphi;
+ const size_t iphi = ivert % sphere->nphi;
double angle[2];
double sin_angle[2];
double cos_angle[2];
- double uv[2];
- int iform;
- ASSERT(itheta < sshape->ntheta);
- if(iphi == sshape->nphi - 1) { /* Polar vertices: Force polar coordinates */
+ if(iphi == sphere->nphi - 1) { /* Polar vertices: Force polar coordinates */
angle[0] = PI;
angle[1] = PI/2;
} else {
@@ -87,40 +117,35 @@ super_shape_get_vertex_position(const unsigned ivert, float vertex[3], void* ctx
angle[1] = -PI/2 + (double)iphi * STEP;
}
- FOR_EACH(iform, 0, 2) {
- double m, k, g;
- float* form = sshape->super_formulas[iform];
- m = cos(form[M] * angle[iform] / 4.0);
- m = fabs(m) / form[A];
- k = sin(form[M] * angle[iform] / 4.0);
- k = fabs(k) / form[B];
- g = pow(m, form[N1]) + pow(k, form[N2]);
- uv[iform] = pow(g, (-1.0/form[N0]));
-
- cos_angle[iform] = cos(angle[iform]);
- sin_angle[iform] = sin(angle[iform]);
- }
+ cos_angle[0] = cos(angle[0]);
+ cos_angle[1] = cos(angle[1]);
+ sin_angle[0] = sin(angle[0]);
+ sin_angle[1] = sin(angle[1]);
- vertex[0] = (float)(uv[0] * cos_angle[0] * uv[1] * cos_angle[1]);
- vertex[1] = (float)(uv[0] * sin_angle[0] * uv[1] * cos_angle[1]);
- vertex[2] = (float)(uv[1] * sin_angle[1]);
+ vertex[0] = (float)(cos_angle[0] * cos_angle[1] * sphere->radius);
+ vertex[1] = (float)(sin_angle[0] * cos_angle[1] * sphere->radius);
+ vertex[2] = (float)(sin_angle[1] * sphere->radius);
}
static res_T
-sample_geometry
+sample_sphere
(struct ssp_rng* rng,
struct sschiff_material* mtl,
struct s3d_shape* shape,
void* sampler_context)
{
struct s3d_vertex_data attrib;
- struct super_shape sshape;
+ struct sphere sphere;
+ struct sphere_context* ctx = sampler_context;
+
const size_t ntheta = (int)((2*PI)/STEP) + 1;
const size_t nphi = (int)((PI)/STEP) + 1;
const size_t nverts = ntheta * nphi;
const size_t nprims = ntheta * (nphi - 1)/*#quads*/ * 2/*#triangles*/;
- (void)sampler_context;
+ sphere.radius = ssp_ran_lognormal(rng, log(ctx->mean_radius), log(1.18));
+ sphere.ntheta = ntheta;
+ sphere.nphi = nphi;
NCHECK(rng, NULL);
NCHECK(mtl, NULL);
@@ -128,38 +153,158 @@ sample_geometry
attrib.usage = S3D_POSITION;
attrib.type = S3D_FLOAT3;
- attrib.get = super_shape_get_vertex_position;
-
- sshape.super_formulas[0][A] = 1;
- sshape.super_formulas[0][B] = 1;
- sshape.super_formulas[0][M] = (float)ssp_rng_uniform_double(rng, 0.1, 6);
- sshape.super_formulas[0][N0] = 1;
- sshape.super_formulas[0][N1] = 1;
- sshape.super_formulas[0][N2] = (float)ssp_rng_uniform_double(rng, 0.1, 6);
-
- sshape.super_formulas[1][A] = 1;
- sshape.super_formulas[1][B] = 1;
- sshape.super_formulas[1][M] = (float)ssp_rng_uniform_double(rng, 0.1, 6);
- sshape.super_formulas[1][N0] = 1;
- sshape.super_formulas[1][N1] = 1;
- sshape.super_formulas[1][N2] = (float)ssp_rng_uniform_double(rng, 0.1, 6);
-
- sshape.ntheta = ntheta;
- sshape.nphi = nphi;
-
- *mtl = SSCHIFF_NULL_MATERIAL;
+ attrib.get = sphere_get_vertex_position;
+
+ mtl->get_property = get_material_property;
+ mtl->material = ctx;
+
return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims,
- super_shape_get_indices, (unsigned)nverts, &attrib, 1, &sshape);
+ sphere_get_indices, (unsigned)nverts, &attrib, 1, &sphere);
+}
+
+/*******************************************************************************
+ * Cylinder geometry
+ ******************************************************************************/
+static void
+cylinder_get_indices(const unsigned itri, unsigned ids[3], void* ctx)
+{
+ struct cylinder* cylinder = ctx;
+ const size_t id = itri * 3;
+
+ CHECK(itri < sa_size(cylinder->indices) / 3, 1);
+ ids[0] = cylinder->indices[id+0];
+ ids[1] = cylinder->indices[id+1];
+ ids[2] = cylinder->indices[id+2];
+}
+
+static void
+cylinder_get_vertex_position(const unsigned ivert, float vertex[3], void* ctx)
+{
+ struct cylinder* cylinder = ctx;
+ const size_t i = ivert * 3;
+
+ CHECK(ivert < sa_size(cylinder->vertices) / 3, 1);
+ f3_set(vertex, cylinder->vertices + i);
+}
+
+static res_T
+sample_cylinder
+ (struct ssp_rng* rng,
+ struct sschiff_material* mtl,
+ struct s3d_shape* shape,
+ void* sampler_context)
+{
+ struct s3d_vertex_data attrib;
+ struct cylinder cylinder = { NULL, NULL };
+ const unsigned nsteps = (unsigned)((2*PI)/STEP);
+ unsigned istep;
+ double sample;
+ const float aspect_ratio = 0.263f;
+ float radius;
+ float height;
+ res_T res = RES_OK;
+ (void)rng, (void)sampler_context;
+
+ sample = ssp_ran_lognormal(rng, log(0.983), log(1.1374));
+ radius = (float)(sample * pow(aspect_ratio, 1.0/3.0));
+ height = 2 * radius * aspect_ratio;
+
+ /* Generate the vertex coordinates */
+ FOR_EACH(istep, 0, nsteps) {
+ const float theta = (float)istep * STEP;
+ const float x = (float)(radius * cos(theta));
+ const float y = (float)(radius * sin(theta));
+ f3(sa_add(cylinder.vertices, 3), x, 0.f, y);
+ f3(sa_add(cylinder.vertices, 3), x, height, y);
+ }
+ /* "Polar" vertices */
+ f3_splat(sa_add(cylinder.vertices, 3), 0.f);
+ f3(sa_add(cylinder.vertices, 3), 0.f, height, 0.f);
+
+ /* Cylinder primitives */
+ FOR_EACH(istep, 0, nsteps) {
+ const unsigned id = istep * 2;
+ unsigned* iprim;
+
+ iprim = sa_add(cylinder.indices, 3);
+ iprim[0] = (id + 0);
+ iprim[1] = (id + 1);
+ iprim[2] = (id + 2) % (nsteps*2);
+
+ iprim = sa_add(cylinder.indices, 3);
+ iprim[0] = (id + 2) % (nsteps*2);
+ iprim[1] = (id + 1);
+ iprim[2] = (id + 3) % (nsteps*2);
+ }
+ /* Cap primitives */
+ FOR_EACH(istep, 0, nsteps) {
+ const unsigned id = istep * 2;
+ unsigned* iprim;
+
+ iprim = sa_add(cylinder.indices, 3);
+ iprim[0] = (nsteps * 2);
+ iprim[1] = (id + 0);
+ iprim[2] = (id + 2) % (nsteps*2);
+
+ iprim = sa_add(cylinder.indices, 3);
+ iprim[0] = (nsteps * 2) + 1;
+ iprim[1] = (id + 3) % (nsteps*2);
+ iprim[2] = (id + 1);
+ }
+
+ attrib.usage = S3D_POSITION;
+ attrib.type = S3D_FLOAT3;
+ attrib.get = cylinder_get_vertex_position;
+
+ mtl->get_property = get_material_property;
+ mtl->material = NULL;
+
+ res = s3d_mesh_setup_indexed_vertices(shape, (unsigned)nsteps*4,
+ cylinder_get_indices, (unsigned)nsteps*2 + 2, &attrib, 1, &cylinder);
+
+ sa_release(cylinder.vertices);
+ sa_release(cylinder.indices);
+
+ return res;
}
int
main(int argc, char** argv)
{
+ char buf[64];
struct mem_allocator allocator;
+ struct sphere_context sphere_ctx;
struct sschiff_device* dev;
struct sschiff_estimator* estimator;
struct sschiff_integrator_desc desc = SSCHIFF_NULL_INTEGRATOR_DESC;
+ struct sschiff_estimator_status status;
struct ssp_rng* rng;
+ struct time t0, t1;
+ size_t i;
+ const double sphere_data_1_01[][2] = { /* mean_radius, result */
+ { 1.00, 6.43e-13 },
+ { 3.22, 2.66e-11 },
+ { 5.44, 1.34e-10 },
+ { 7.67, 3.56e-10 },
+ { 9.89, 6.83e-10 },
+ { 12.1, 1.08e-9 },
+ { 14.3, 1.53e-9 },
+ { 16.6, 2.00e-9 },
+ { 18.8, 2.51e-9 },
+ { 21.0, 3.07e-9 }
+ };
+ const double sphere_data_1_1[][2] = { /* mean_radius, result */
+ { 1.00, 8.39e-12 },
+ { 3.22, 6.93e-11 },
+ { 5.44, 1.97e-10 },
+ { 7.67, 3.92e-10 },
+ { 9.89, 6.51e-10 },
+ { 12.1, 9.75e-10 },
+ { 14.3, 1.37e-9 },
+ { 16.6, 1.82e-9 },
+ { 18.8, 2.34e-9 },
+ { 21.0, 2.93e-9 }
+ };
(void)argc, (void)argv;
mem_init_proxy_allocator(&allocator, &mem_default_allocator);
@@ -167,7 +312,7 @@ main(int argc, char** argv)
CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng), RES_OK);
CHECK(sschiff_device_create(NULL, &allocator, 0, NULL, &dev), RES_OK);
- desc.sample_geometry = sample_geometry;
+ desc.sample_geometry = sample_sphere;
desc.scattering_angles_count = 1;
CHECK(sschiff_integrate(NULL, NULL, NULL, 0, 0, NULL), RES_BAD_ARG);
CHECK(sschiff_integrate(dev, NULL, NULL, 0, 0, NULL), RES_BAD_ARG);
@@ -184,7 +329,29 @@ main(int argc, char** argv)
CHECK(sschiff_integrate(NULL, NULL, &desc, 0, 1, &estimator), RES_BAD_ARG);
CHECK(sschiff_integrate(dev, NULL, &desc, 0, 1, &estimator), RES_BAD_ARG);
CHECK(sschiff_integrate(NULL, rng, &desc, 0, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, &desc, 0, 100, &estimator), RES_OK);
+
+ desc.sampler_context = &sphere_ctx;
+
+ sphere_ctx.props.m = 0.6e-6;
+ sphere_ctx.props.n_r = 1.1;
+ sphere_ctx.props.kappa_r = 0.004;
+ sphere_ctx.props.n_medium = 4.0/3.0;
+ sphere_ctx.mean_radius = 1.0;
+
+ time_current(&t0);
+ CHECK(sschiff_integrate(dev, rng, &desc, 6.e-7, NREALISATIONS, &estimator), RES_OK);
+ time_current(&t1);
+
+ CHECK(sschiff_estimator_get_extinction_cross_section(NULL, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_extinction_cross_section(estimator, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_extinction_cross_section(NULL, &status), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_extinction_cross_section(estimator, &status), RES_OK);
+
+ printf("Extinction cross section: %g +/- %g\n", status.E, status.SE);
+ CHECK(eq_eps(status.E, 8.3e-12, 2*status.SE), 1);
+ time_sub(&t0, &t1, &t0);
+ time_dump(&t0, TIME_MIN|TIME_SEC|TIME_MSEC, NULL, buf, sizeof(buf));
+ printf("Elapsed time: %s\n", buf);
CHECK(sschiff_estimator_ref_get(NULL), RES_BAD_ARG);
CHECK(sschiff_estimator_ref_get(estimator), RES_OK);
@@ -192,6 +359,29 @@ main(int argc, char** argv)
CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
+ sphere_ctx.props.n_r = 1.1;
+ printf("\nMean radius: 1.1\n");
+ FOR_EACH(i, 0, sizeof(sphere_data_1_1)/sizeof(double[2])) {
+ sphere_ctx.mean_radius = sphere_data_1_1[i][0];
+ CHECK(sschiff_integrate(dev, rng, &desc, 6.e-7, NREALISATIONS, &estimator), RES_OK);
+ CHECK(sschiff_estimator_get_extinction_cross_section(estimator, &status), RES_OK);
+ printf("%d - Extinction cross section = %g ~ %g +/- %g\n",
+ (int)i, sphere_data_1_1[i][1], status.E, status.SE);
+ CHECK(eq_eps(status.E, sphere_data_1_1[i][1], 2*status.SE), 1);
+ CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
+ }
+ sphere_ctx.props.n_r = 1.01;
+ printf("\nMean radius: 1.01\n");
+ FOR_EACH(i, 0, sizeof(sphere_data_1_01)/sizeof(double[2])) {
+ sphere_ctx.mean_radius = sphere_data_1_01[i][0];
+ CHECK(sschiff_integrate(dev, rng, &desc, 6.e-7, NREALISATIONS, &estimator), RES_OK);
+ CHECK(sschiff_estimator_get_extinction_cross_section(estimator, &status), RES_OK);
+ printf("%d - Extinction cross section = %g ~ %g +/- %g\n",
+ (int)i, sphere_data_1_01[i][1], status.E, status.SE);
+ CHECK(eq_eps(status.E, sphere_data_1_01[i][1], 2*status.SE), 1);
+ CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
+ }
+
CHECK(sschiff_device_ref_put(dev), RES_OK);
CHECK(ssp_rng_ref_put(rng), RES_OK);
