commit 9899bb031bd2ec30d32691d48f6d4e2715b534fe
parent 831d21d76495d54c465f5c3fcdcc47c679b93af1
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 22 Oct 2015 09:16:51 +0200
New wavelengths integration API
The Schiff integration is performed on a list of wavelengths. All
wavelenghts share the same geometry and direction sampling as well as
the computation of the radiative path length.
Diffstat:
4 files changed, 502 insertions(+), 288 deletions(-)
diff --git a/src/sschiff.h b/src/sschiff.h
@@ -56,6 +56,15 @@ struct s3d_device;
struct s3d_shape;
struct ssp_rng;
+/* Type of estimated data */
+enum sschiff_data {
+ SSCHIFF_EXTINCTION_CROSS_SECTION,
+ SSCHIFF_ABSORPTION_CROSS_SECTION,
+ SSCHIFF_SCATTERING_CROSS_SECTION,
+ SSCHIFF_AVERAGE_PROJECTIVE_AREA,
+ SSCHIFF_DATA_COUNT__
+};
+
/* Optical properties of a material handled by the Schiff integrator */
struct sschiff_material_properties {
double medium_refractive_index;
@@ -88,12 +97,14 @@ struct sschiff_geometry_distribution {
static const struct sschiff_geometry_distribution
SSCHIFF_NULL_GEOMETRY_DISTRIBUTION = { NULL, NULL };
-/* Result of the Monte Carlo estimation */
-struct sschiff_estimator_status {
- double E; /* Expected value */
- double V; /* Variance */
- double SE; /* Standard error */
- size_t nsteps; /* # Monte Carlo weights */
+/* State of the Monte Carlo estimation */
+struct sschiff_estimator_state {
+ double wavelength; /* Estimated wavelength */
+ struct sschiff_estimator_value { /* Values of the estimation */
+ double E; /* Expected value */
+ double V; /* Variance */
+ double SE; /* Standard error */
+ } values[SSCHIFF_DATA_COUNT__];
};
/* Opaque types */
@@ -126,7 +137,8 @@ sschiff_integrate
(struct sschiff_device* dev,
struct ssp_rng* rng,
struct sschiff_geometry_distribution* distribution,
- const double wavelength, /* Wavelenght to estimate */
+ const double* wavelengths, /* list of wavelengths to estimate */
+ const size_t wavelengths_count, /* # wavelengths to estimate */
const size_t scattering_angles_count, /* # of scattering angles to handle */
const size_t sampled_geometries_count, /* # geometries to sample */
const size_t sampled_directions_count, /* # directions to sample per geometry */
@@ -140,32 +152,33 @@ SSCHIFF_API res_T
sschiff_estimator_ref_put
(struct sschiff_estimator* estimator);
+/* Return the number of estimated wavelength */
SSCHIFF_API res_T
-sschiff_estimator_get_extinction_cross_section
- (struct sschiff_estimator* estimator,
- struct sschiff_estimator_status* status);
+sschiff_estimator_get_wavelengths_count
+ (const struct sschiff_estimator* estimator,
+ size_t* count);
+/* Return the number of Monte Carlo realisations */
SSCHIFF_API res_T
-sschiff_estimator_get_absorption_cross_section
- (struct sschiff_estimator* estimator,
- struct sschiff_estimator_status* status);
+sschiff_estimator_get_realisations_count
+ (const struct sschiff_estimator* estimator,
+ size_t* count);
+/* Retrieve the estimation state of a given wavelength */
SSCHIFF_API res_T
-sschiff_estimator_get_scattering_cross_section
- (struct sschiff_estimator* estimator,
- struct sschiff_estimator_status* status);
-
-SSCHIFF_API res_T
-sschiff_estimator_get_average_projected_area
- (struct sschiff_estimator* estimator,
- struct sschiff_estimator_status* status);
-
-/* TODO Not implemented yet */
+sschiff_estimator_get_wavelength_state
+ (const struct sschiff_estimator* estimator,
+ const double wavelength,
+ struct sschiff_estimator_state* state);
+
+/* Retrieve the estimation state of all wavelengths. The length of `states'
+ * must be at least equal to the count of integrated wavelengths. One can use
+ * the sschiff_estimator_get_wavelengths_count function to retrieve this
+ * information. */
SSCHIFF_API res_T
-sschiff_estimator_get_phase_function
- (struct sschiff_estimator* estimator,
- const size_t iscattering_angle,
- struct sschiff_estimator_status* status);
+sschiff_estimator_get_states
+ (const struct sschiff_estimator* estimator,
+ struct sschiff_estimator_state states[]);
END_DECLS
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -31,6 +31,7 @@
#include "sschiff.h"
#include "sschiff_device.h"
+#include <rsys/dynamic_array_double.h>
#include <rsys/float2.h>
#include <rsys/float3.h>
#include <rsys/image.h>
@@ -46,27 +47,37 @@
#define MAX_FAILURES 10
-/* Type of computed Monte Carlo weights */
-enum weight_type {
- WEIGHT_EXTINCTION_CROSS_SECTION,
- WEIGHT_ABSORPTION_CROSS_SECTION,
- WEIGHT_SCATTERING_CROSS_SECTION,
- WEIGHT_AVERAGE_PROJECTED_AREA,
- WEIGHTS_COUNT
-};
+/* Accumulator of double precision values */
+struct accum { double weights[SSCHIFF_DATA_COUNT__]; };
-/* Accumulator of Monte Carlo weights */
-struct accum {
- double weight;
- double sqr_weight;
- size_t naccums;
-};
+/* Monte carlo data */
+struct mc_data { double weight; double sqr_weight; };
-struct sschiff_estimator {
- double wavelength;
+static const struct mc_data MC_DATA_NULL = { 0.0, 0.0 };
+
+/* Accumulator of Monte Carlo data */
+struct mc_accum { struct mc_data mc_data[SSCHIFF_DATA_COUNT__]; };
+
+static const struct mc_accum MC_ACCUM_NULL =
+{{{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}}};
+
+#define DARRAY_NAME mcaccum
+#define DARRAY_DATA struct mc_accum
+#include <rsys/dynamic_array.h>
+
+#define DARRAY_NAME accum
+#define DARRAY_DATA struct accum
+#include <rsys/dynamic_array.h>
- struct accum accums[WEIGHTS_COUNT];
+#define DARRAY_NAME mtl
+#define DARRAY_DATA struct sschiff_material_properties
+#include <rsys/dynamic_array.h>
+
+struct sschiff_estimator {
+ struct darray_double wavelengths;
+ struct darray_mcaccum accums;
struct sschiff_device* dev;
+ size_t nrealisations;
ref_T ref;
};
@@ -84,6 +95,15 @@ morton2D_decode(const uint32_t u32)
return (uint16_t)x;
}
+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;
+}
+
/* Compute an orthonormal basis where `dir' is the Z axis. */
static void
build_basis(const float dir[3], float basis[9])
@@ -132,6 +152,80 @@ build_transform
transform[11] = -f3_dot(basis+6, pos);
}
+/* Dump a thumbnail of the sampled geometry in given sampled direction */
+static INLINE void
+draw_thumbnail
+ (struct s3d_scene* scn,
+ 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 char* name) /* Filename of the thumbnail */
+{
+#define DEFINITION 128
+ STATIC_ASSERT(IS_POW2(DEFINITION), Invalid_thumbnail_definition);
+ const float range[2] = { 0, FLT_MAX };
+ const float pixel_size = 1.f/(float)DEFINITION;
+ float axis_x[3], axis_y[3], axis_z[3];
+ float plane_size[2];
+ float ray_org[3];
+ float org[3];
+ uint32_t npixels = (uint32_t)(DEFINITION * DEFINITION);
+ uint32_t mcode;
+ unsigned char* img = NULL;
+ ASSERT(basis && pos && lower && upper);
+
+ /* Allocate the thumbnail buffer */
+ img = sa_add(img, DEFINITION*DEFINITION);
+ ASSERT(img != NULL);
+
+ /* Define the projection axis */
+ f2_sub(plane_size, upper, lower);
+ 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]));
+
+ FOR_EACH(mcode, 0, npixels) {
+ struct s3d_hit hit;
+ size_t ipixel;
+ float sample[2], x[3], y[3];
+ uint16_t ipixel_x, ipixel_y;
+
+ /* Compute the sample position in [0, 1)^2 onto the projection plane */
+ ipixel_x = morton2D_decode(mcode);
+ ipixel_y = morton2D_decode(mcode>>1);
+ sample[0] = ((float)ipixel_x + 0.5f) * pixel_size;
+ sample[1] = ((float)ipixel_y + 0.5f) * pixel_size;
+
+ /* Compute the ray origin */
+ f3_mulf(x, axis_x, sample[0]*2.f - 1.f);
+ f3_mulf(y, axis_y, sample[1]*2.f - 1.f);
+ f3_add(ray_org, f3_add(ray_org, x, y), org);
+
+ S3D(scene_trace_ray(scn, ray_org, axis_z, range, &hit));
+
+ /* Simple shading */
+ ipixel = (size_t)(ipixel_y * DEFINITION + ipixel_x);
+ if(S3D_HIT_NONE(&hit)) {
+ img[ipixel] = 0;
+ } else {
+ float N[3] = {0.f,0.f,0.f};
+ float cosine;
+ f3_normalize(N, hit.normal);
+ if(0 > (cosine = f3_dot(N, axis_z))) {
+ cosine = f3_dot(f3_minus(N, N), axis_z);
+ }
+ img[ipixel] = (unsigned char)(cosine*255.f + 0.5f/*Round*/);
+ }
+ }
+ image_ppm_write(name, DEFINITION, DEFINITION, 1, img);
+ sa_release(img);
+#undef DEFINITION
+}
+
/* Compute the length of radiative path that traverses the scene */
static res_T
compute_path_length
@@ -183,17 +277,18 @@ compute_path_length
}
static res_T
-compute_monte_carlo_weight
+accum_monte_carlo_weight
(struct sschiff_device* dev,
struct s3d_scene* scn,
struct ssp_rng* rng,
struct sschiff_material_properties* props,
- const double wavelength,
+ const double* wavelengths,
+ const size_t nwavelengths,
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 */
- double weights[WEIGHTS_COUNT])
+ struct accum* accums)
{
struct s3d_hit hit;
double sample[2];
@@ -205,10 +300,9 @@ compute_monte_carlo_weight
float org[3];
float x[3], y[3];
float rcp_pdf;
- int i;
res_T res = RES_OK;
- ASSERT(scn && rng && props && basis && pos && lower && upper);
- ASSERT(weights);
+ ASSERT(scn && rng && props && wavelengths && accums && basis && pos);
+ ASSERT(lower && upper);
f2_sub(plane_size, upper, lower);
rcp_pdf = plane_size[0] * plane_size[1] * 1.e-12f/*Meter^2 to micron^2*/;
@@ -221,9 +315,6 @@ compute_monte_carlo_weight
/* Compute the origin of the projection plane */
f3_add(org, pos, f3_mulf(org, axis_z, lower[2]));
- FOR_EACH(i, 0, WEIGHTS_COUNT)
- weights[i] = 0;
-
do {
/* Uniformly sample a position onto the projection plane and use it as the
* origin of the ray to trace */
@@ -236,9 +327,8 @@ compute_monte_carlo_weight
S3D(scene_trace_ray(scn, ray_org, axis_z, range, &hit));
if(!S3D_HIT_NONE(&hit)) {
- double n_r, k_r;
- double n_e, k_e, lambda_e;
double path_length;
+ size_t iwlen;
res = compute_path_length(dev, scn, &hit, ray_org, axis_z, &path_length);
if(res != RES_OK) { /* Handle numerical/geometry issues */
@@ -246,24 +336,34 @@ compute_monte_carlo_weight
continue;
}
- n_e = props->medium_refractive_index;
- n_r = props->relative_real_refractive_index;
- k_r = props->relative_imaginary_refractive_index;
+ FOR_EACH(iwlen, 0, nwavelengths) {
+ double n_r, k_r;
+ double n_e, k_e, lambda_e;
+ double w_extinction, w_absorption, w_scattering;
+
+ n_e = props[iwlen].medium_refractive_index;
+ n_r = props[iwlen].relative_real_refractive_index;
+ k_r = props[iwlen].relative_imaginary_refractive_index;
- lambda_e = wavelength / n_e;
- k_e = 2.0 * PI / lambda_e;
+ lambda_e = wavelengths[iwlen] / n_e;
+ k_e = 2.0 * PI / lambda_e;
- weights[WEIGHT_EXTINCTION_CROSS_SECTION] = 2.0 * rcp_pdf
- * (1.0 - exp(-k_e*k_r*path_length) * cos(k_e*(n_r - 1.0)*path_length));
+ /* Extinction cross section */
+ w_extinction = 2.0 * rcp_pdf *
+ (1.0 - exp(-k_e*k_r*path_length) * cos(k_e*(n_r - 1.0)*path_length));
+ accums[iwlen].weights[SSCHIFF_EXTINCTION_CROSS_SECTION] += w_extinction;
- weights[WEIGHT_ABSORPTION_CROSS_SECTION] = rcp_pdf
- * (1.0 - exp(-2.0*k_e*k_r*path_length));
+ /* Absorption cross section */
+ w_absorption = rcp_pdf * (1.0 - exp(-2.0*k_e*k_r*path_length));
+ accums[iwlen].weights[SSCHIFF_ABSORPTION_CROSS_SECTION] += w_absorption;
- weights[WEIGHT_SCATTERING_CROSS_SECTION] =
- weights[WEIGHT_EXTINCTION_CROSS_SECTION]
- - weights[WEIGHT_ABSORPTION_CROSS_SECTION];
+ /* Scattering cross section */
+ w_scattering = w_extinction - w_absorption;
+ accums[iwlen].weights[SSCHIFF_SCATTERING_CROSS_SECTION] += w_scattering;
- weights[WEIGHT_AVERAGE_PROJECTED_AREA] = rcp_pdf;
+ /* Average projected area */
+ accums[iwlen].weights[SSCHIFF_AVERAGE_PROJECTIVE_AREA] += rcp_pdf;
+ }
}
} while(res != RES_OK && nfailures < MAX_FAILURES);
@@ -275,99 +375,28 @@ compute_monte_carlo_weight
return res;
}
-/* Dump a thumbnail of the sampled geometry in given sampled direction */
-static INLINE void
-draw_thumbnail
- (struct s3d_scene* scn,
- 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 char* name) /* Filename of the thumbnail */
-{
-#define DEFINITION 128
- STATIC_ASSERT(IS_POW2(DEFINITION), Invalid_thumbnail_definition);
- const float range[2] = { 0, FLT_MAX };
- const float pixel_size = 1.f/(float)DEFINITION;
- float axis_x[3], axis_y[3], axis_z[3];
- float plane_size[2];
- float ray_org[3];
- float org[3];
- uint32_t npixels = (uint32_t)(DEFINITION * DEFINITION);
- uint32_t mcode;
- unsigned char* img = NULL;
- ASSERT(basis && pos && lower && upper);
-
- /* Allocate the thumbnail buffer */
- img = sa_add(img, DEFINITION*DEFINITION);
- ASSERT(img != NULL);
-
- /* Define the projection axis */
- f2_sub(plane_size, upper, lower);
- 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]));
-
- FOR_EACH(mcode, 0, npixels) {
- struct s3d_hit hit;
- size_t ipixel;
- float sample[2], x[3], y[3];
- uint16_t ipixel_x, ipixel_y;
-
- /* Compute the sample position in [0, 1)^2 onto the projection plane */
- ipixel_x = morton2D_decode(mcode);
- ipixel_y = morton2D_decode(mcode>>1);
- sample[0] = ((float)ipixel_x + 0.5f) * pixel_size;
- sample[1] = ((float)ipixel_y + 0.5f) * pixel_size;
-
- /* Compute the ray origin */
- f3_mulf(x, axis_x, sample[0]*2.f - 1.f);
- f3_mulf(y, axis_y, sample[1]*2.f - 1.f);
- f3_add(ray_org, f3_add(ray_org, x, y), org);
-
- S3D(scene_trace_ray(scn, ray_org, axis_z, range, &hit));
-
- /* Simple shading */
- ipixel = (size_t)(ipixel_y * DEFINITION + ipixel_x);
- if(S3D_HIT_NONE(&hit)) {
- img[ipixel] = 0;
- } else {
- float N[3] = {0.f,0.f,0.f};
- float cosine;
- f3_normalize(N, hit.normal);
- if(0 > (cosine = f3_dot(N, axis_z))) {
- cosine = f3_dot(f3_minus(N, N), axis_z);
- }
- img[ipixel] = (unsigned char)(cosine*255.f + 0.5f/*Round*/);
- }
- }
- image_ppm_write(name, DEFINITION, DEFINITION, 1, img);
- sa_release(img);
-#undef DEFINITION
-}
-
static res_T
radiative_properties
(struct sschiff_device* dev,
struct ssp_rng* rng,
const int istep,
- const double wavelength,
+ const double* wavelengths, /* Sorted in ascending order */
+ const size_t nwavelengths,
const size_t ndirs,
struct s3d_scene* scene,
struct sschiff_material_properties* props,
- struct accum accums[WEIGHTS_COUNT])
+ struct mc_accum* mc_accums,
+ struct accum* accums)
{
- double weights_dirs[WEIGHTS_COUNT];
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 && rng && ndirs && props && accums);
+ ASSERT(dev && rng && wavelengths && nwavelengths && ndirs && props);
+ ASSERT(mc_accums && accums);
(void)istep;
S3D(scene_get_aabb(scene, lower, upper));
@@ -390,9 +419,8 @@ radiative_properties
f3_mulf(centroid, f3_add(centroid, lower, upper), 0.5f);
- memset(weights_dirs, 0, sizeof(weights_dirs));
+ memset(accums, 0, sizeof(struct accum)*nwavelengths);
FOR_EACH(idir, 0, ndirs) {
- double weights[WEIGHTS_COUNT];
float dir[4];
float basis[9];
float transform[12];
@@ -417,11 +445,10 @@ radiative_properties
f3_max(upper, upper, pt[i]);
}
- res = compute_monte_carlo_weight(dev, scene, rng, props, wavelength, basis,
- centroid, lower, upper, weights);
+ res = accum_monte_carlo_weight(dev, scene, rng, props, wavelengths,
+ nwavelengths, basis, centroid, lower, upper, accums);
if(res != RES_OK) goto error;
- FOR_EACH(i, 0, WEIGHTS_COUNT) weights_dirs[i] += weights[i];
#if 0
{
/* Compute an image of the shape transformed in the footprint space */
@@ -431,28 +458,32 @@ radiative_properties
}
#endif
}
- FOR_EACH(i, 0, WEIGHTS_COUNT) {
- weights_dirs[i] /= (double)ndirs;
- accums[i].weight += weights_dirs[i];
- accums[i].sqr_weight += weights_dirs[i] * weights_dirs[i];
- accums[i].naccums += 1;
+ FOR_EACH(iwlen, 0, nwavelengths) {
+ int iweight;
+ FOR_EACH(iweight, 0, SSCHIFF_DATA_COUNT__) {
+ const double w = accums[iwlen].weights[iweight] / (double)ndirs;
+ mc_accums[iwlen].mc_data[iweight].weight += w;
+ mc_accums[iwlen].mc_data[iweight].sqr_weight += w * w;
+ }
}
exit:
return res;
error:
- memset(accums, 0, sizeof(struct accum)*WEIGHTS_COUNT);
+ memset(mc_accums, 0, sizeof(struct mc_accum)*nwavelengths);
goto exit;
}
static void
-get_status(const struct accum* acc, struct sschiff_estimator_status* status)
+get_mc_value
+ (const struct mc_data* data,
+ const size_t nrealisations,
+ struct sschiff_estimator_value* value)
{
- ASSERT(acc && status);
- 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;
+ ASSERT(data && value);
+ value->E = data->weight / (double)nrealisations;
+ value->V = data->sqr_weight / (double)nrealisations
+ - (data->weight * data->weight) / (double)(nrealisations * nrealisations);
+ value->SE = sqrt(value->V / (double)nrealisations);
}
static char
@@ -470,13 +501,17 @@ estimator_release(ref_T* ref)
ASSERT(ref);
estimator = CONTAINER_OF(ref, struct sschiff_estimator, ref);
dev = estimator->dev;
+ darray_double_release(&estimator->wavelengths);
+ darray_mcaccum_release(&estimator->accums);
MEM_RM(dev->allocator, estimator);
SSCHIFF(device_ref_put(dev));
}
static res_T
estimator_create
- (struct sschiff_device* dev, struct sschiff_estimator** out_estimator)
+ (struct sschiff_device* dev,
+ const size_t nwavelengths,
+ struct sschiff_estimator** out_estimator)
{
struct sschiff_estimator* estimator = NULL;
res_T res = RES_OK;
@@ -492,6 +527,22 @@ estimator_create
ref_init(&estimator->ref);
SSCHIFF(device_ref_get(dev));
estimator->dev = dev;
+ darray_double_init(dev->allocator, &estimator->wavelengths);
+ darray_mcaccum_init(dev->allocator, &estimator->accums);
+
+ res = darray_double_resize(&estimator->wavelengths, nwavelengths);
+ if(res != RES_OK) {
+ log_error(dev,
+ "Not enough memory: couldn't allocate the list of wavelengths.\n");
+ goto error;
+ }
+ res = darray_mcaccum_resize(&estimator->accums, nwavelengths);
+ if(res != RES_OK) {
+ log_error(dev,
+ "Not enough memory: couldn't allocate the per wavelength "
+ "Monte Carlo accumulators.\n");
+ goto error;
+ }
exit:
*out_estimator = estimator;
@@ -512,27 +563,49 @@ sschiff_integrate
(struct sschiff_device* dev,
struct ssp_rng* rng,
struct sschiff_geometry_distribution* distrib,
- const double wavelength,
- const size_t scattering_angles_count,
- const size_t sampled_geometries_count,
- const size_t sampled_directions_count,
+ const double* wavelengths,
+ const size_t nwavelengths,
+ const size_t nangles,
+ const size_t ngeoms,
+ const size_t ndirs,
struct sschiff_estimator** out_estimator)
{
+ struct darray_mtl mtls;
+ struct darray_accum accums;
struct sschiff_estimator* estimator = NULL;
struct s3d_shape* shape = NULL;
struct s3d_scene* scene = NULL;
size_t i;
res_T res = RES_OK;
- (void)scattering_angles_count;
+ (void)nangles;
+
+ if(!dev || !rng || !distrib || !wavelengths || !nwavelengths || !nangles
+ || !ngeoms || !ndirs || !out_estimator || !check_distribution(distrib)) {
+ return RES_BAD_ARG;
+ }
- if(!dev || !rng || !distrib || !check_distribution(distrib) || !out_estimator) {
- res = RES_BAD_ARG;
+ darray_mtl_init(dev->allocator, &mtls);
+ darray_accum_init(dev->allocator, &accums);
+
+ res = darray_mtl_resize(&mtls, nwavelengths);
+ if(res != RES_OK) {
+ log_error(dev,
+ "No enough memory: couldn't allocate the per wavelength"
+ "optical properties.\n");
goto error;
}
- res = estimator_create(dev, &estimator);
+ res = darray_accum_resize(&accums, nwavelengths);
+ if(res != RES_OK) {
+ log_error(dev,
+ "No enough memory: couldn't allocate the per wavelength "
+ "temporary accumulators.\n");
+ goto error;
+ }
+
+ res = estimator_create(dev, nwavelengths, &estimator);
if(res != RES_OK) goto error;
- estimator->wavelength = wavelength;
+ estimator->nrealisations = ngeoms;
res = s3d_shape_create_mesh(dev->s3d, &shape);
if(res != RES_OK) {
@@ -551,10 +624,24 @@ sschiff_integrate
goto error;
}
- memset(estimator->accums, 0, sizeof(struct accum)*WEIGHTS_COUNT);
- FOR_EACH(i, 0, sampled_geometries_count) {
+ /* Clear the estimator accumulators */
+ memset(darray_mcaccum_data_get(&estimator->accums), 0,
+ sizeof(struct mc_accum)*nwavelengths);
+
+ /* Setup and sort the wavelengths */
+ FOR_EACH(i, 0, nwavelengths) {
+ darray_double_data_get(&estimator->wavelengths)[i] = wavelengths[i];
+ }
+ qsort(darray_double_data_get(&estimator->wavelengths), nwavelengths,
+ sizeof(double), cmp_double);
+
+ FOR_EACH(i, 0, ngeoms) {
+ struct sschiff_material_properties* props = darray_mtl_data_get(&mtls);
+ const double* wlengths = darray_double_cdata_get(&estimator->wavelengths);
+ struct mc_accum* mc_accums = darray_mcaccum_data_get(&estimator->accums);
+ struct accum* accs = darray_accum_data_get(&accums);
struct sschiff_material material = SSCHIFF_NULL_MATERIAL;
- struct sschiff_material_properties props;
+ size_t iwavelength;
/* Sample a geometry, i.e. a shape and its associated material */
res = distrib->sample(rng, &material, shape, distrib->context);
@@ -562,20 +649,26 @@ sschiff_integrate
log_error(dev, "Error in sampling a Schiff geometry.\n");
goto error;
}
- material.get_property(material.material, wavelength, &props);
+ /* Fetch the per wavelength material properties */
+ FOR_EACH(iwavelength, 0, nwavelengths) {
+ material.get_property
+ (material.material, wavelengths[iwavelength], props + iwavelength);
+ }
/* Build the Star-3D representation of the sampled shape */
S3D(scene_begin_session(scene, S3D_TRACE));
/* Schiff Estimation */
- res = radiative_properties(dev, rng, (int)i, wavelength,
- sampled_directions_count, scene, &props, estimator->accums);
+ res = radiative_properties(dev, rng, (int)i, wlengths, nwavelengths, ndirs,
+ scene, props, mc_accums, accs);
if(res != RES_OK) goto error;
S3D(scene_end_session(scene));
}
exit:
+ darray_mtl_release(&mtls);
+ darray_accum_release(&accums);
if(out_estimator) *out_estimator = estimator;
if(shape) S3D(shape_ref_put(shape));
if(scene) S3D(scene_ref_put(scene));
@@ -605,42 +698,80 @@ sschiff_estimator_ref_put(struct sschiff_estimator* estimator)
}
res_T
-sschiff_estimator_get_extinction_cross_section
- (struct sschiff_estimator* estimator,
- struct sschiff_estimator_status* status)
+sschiff_estimator_get_wavelengths_count
+ (const struct sschiff_estimator* estimator, size_t* count)
{
- if(!estimator || !status) return RES_BAD_ARG;
- get_status(estimator->accums + WEIGHT_EXTINCTION_CROSS_SECTION, status);
+ if(!estimator || !count) return RES_BAD_ARG;
+ *count = darray_double_size_get(&estimator->wavelengths);
return RES_OK;
}
res_T
-sschiff_estimator_get_absorption_cross_section
- (struct sschiff_estimator* estimator,
- struct sschiff_estimator_status* status)
+sschiff_estimator_get_realisations_count
+ (const struct sschiff_estimator* estimator, size_t* count)
{
- if(!estimator || !status) return RES_BAD_ARG;
- get_status(estimator->accums + WEIGHT_ABSORPTION_CROSS_SECTION, status);
+ if(!estimator || !count) return RES_BAD_ARG;
+ *count = estimator->nrealisations;
return RES_OK;
}
res_T
-sschiff_estimator_get_scattering_cross_section
- (struct sschiff_estimator* estimator,
- struct sschiff_estimator_status* status)
+sschiff_estimator_get_wavelength_state
+ (const struct sschiff_estimator* estimator,
+ const double wavelength,
+ struct sschiff_estimator_state* state)
{
- if(!estimator || !status) return RES_BAD_ARG;
- get_status(estimator->accums + WEIGHT_SCATTERING_CROSS_SECTION, status);
+ const struct mc_accum* accum;
+ const double* wavelengths;
+ const double* find;
+ size_t iwavelength;
+ size_t nwavelengths;
+ int i;
+
+ if(!estimator || !state)
+ return RES_BAD_ARG;
+
+ wavelengths = darray_double_cdata_get(&estimator->wavelengths);
+ nwavelengths = darray_double_size_get(&estimator->wavelengths);
+ find = bsearch
+ (&wavelength, wavelengths, nwavelengths, sizeof(double), cmp_double);
+ if(!find) return RES_BAD_ARG;
+
+ iwavelength = (size_t)(find - wavelengths);
+ accum = darray_mcaccum_cdata_get(&estimator->accums) + iwavelength;
+
+ FOR_EACH(i, 0, SSCHIFF_DATA_COUNT__) {
+ get_mc_value(accum->mc_data+i, estimator->nrealisations, state->values+i);
+ }
+ state->wavelength = wavelength;
return RES_OK;
}
res_T
-sschiff_estimator_get_average_projected_area
- (struct sschiff_estimator* estimator,
- struct sschiff_estimator_status* status)
+sschiff_estimator_get_states
+ (const struct sschiff_estimator* estimator,
+ struct sschiff_estimator_state* states)
{
- if(!estimator || !status) return RES_BAD_ARG;
- get_status(estimator->accums + WEIGHT_AVERAGE_PROJECTED_AREA, status);
+ const double* wavelengths;
+ const struct mc_accum* accums;
+ size_t nwavelengths;
+ size_t iwlen;
+ int i;
+ if(!estimator || !states) return RES_BAD_ARG;
+
+ nwavelengths = darray_double_size_get(&estimator->wavelengths);
+ wavelengths = darray_double_cdata_get(&estimator->wavelengths);
+ accums = darray_mcaccum_cdata_get(&estimator->accums);
+
+ FOR_EACH(iwlen, 0, nwavelengths) {
+ states[iwlen].wavelength = wavelengths[iwlen];
+ FOR_EACH(i, 0, SSCHIFF_DATA_COUNT__) {
+ get_mc_value
+ (accums[iwlen].mc_data + i,
+ estimator->nrealisations,
+ states[iwlen].values + i);
+ }
+ }
return RES_OK;
}
diff --git a/src/test_sschiff_estimator_cylinder.c b/src/test_sschiff_estimator_cylinder.c
@@ -37,7 +37,6 @@
#include <star/ssp.h>
struct result {
- double wavelength;
double extinction_E; /* Expected value */
double extinction_SE; /* Standard error */
};
@@ -809,19 +808,30 @@ main(int argc, char** argv)
struct sschiff_device* dev;
struct sschiff_geometry_distribution distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION;
struct sschiff_estimator* estimator;
- struct sschiff_estimator_status status;
struct ssp_rng* rng;
struct time t0, t1;
struct result results[] = {
- { 4.00E-7, 2.76750904574817e-12, 2.34364058801497e-15 },
- { 4.50e-7, 2.50099694990826e-12, 2.20010054373943e-15 },
- { 5.00e-7, 2.32782435081181e-12, 2.05630225764867e-15 },
- { 6.00e-7, 1.96661661205101e-12, 1.82091931178789e-15 },
- { 6.50e-7, 1.71149696366194e-12, 1.63845385406345e-15 },
- { 7.00e-7, 1.54229472112306e-12, 1.49426946826458e-15 },
- { 7.50e-7, 1.39877652964399e-12, 1.36729344819050e-15 }
+ { 2.76750904574817e-12, 2.34364058801497e-15 },
+ { 2.50099694990826e-12, 2.20010054373943e-15 },
+ { 2.32782435081181e-12, 2.05630225764867e-15 },
+ { 1.96661661205101e-12, 1.82091931178789e-15 },
+ { 1.71149696366194e-12, 1.63845385406345e-15 },
+ { 1.54229472112306e-12, 1.49426946826458e-15 },
+ { 1.39877652964399e-12, 1.36729344819050e-15 }
};
- const size_t nresults = sizeof(results) / sizeof(struct result);
+ double wavelengths[] = {
+ 4.00E-7,
+ 4.50e-7,
+ 5.00e-7,
+ 6.00e-7,
+ 6.50e-7,
+ 7.00e-7,
+ 7.50e-7
+ };
+ const size_t nwavelengths = sizeof(wavelengths) / sizeof(double);
+ const size_t nscatt_angles = 1;
+ const size_t ngeoms = 1000;
+ const size_t ndirs = 100;
size_t i;
(void)argc, (void)argv;
@@ -838,45 +848,45 @@ main(int argc, char** argv)
distrib.sample = sample_cylinder;
distrib.context = &sampler_ctx;
- FOR_EACH(i, 0, nresults) {
- const size_t nscatt_angles = 1;
- const size_t ngeoms = 1000;
- const size_t ndirs = 10;
- const double wavelength = results[i].wavelength;
+ time_current(&t0);
+ CHECK(sschiff_integrate(dev, rng, &distrib, wavelengths, nwavelengths,
+ nscatt_angles, ngeoms, ndirs, &estimator), RES_OK);
+ time_current(&t1);
+ time_sub(&t0, &t1, &t0);
+ time_dump(&t0, TIME_MIN|TIME_SEC|TIME_MSEC, NULL, buf, sizeof(buf));
+
+ FOR_EACH(i, 0, nwavelengths) {
+ struct sschiff_estimator_state state;
double interval0[2], interval1[2], interval2[2];
+ struct sschiff_estimator_value* val;
- time_current(&t0);
- CHECK(sschiff_integrate(dev, rng, &distrib, wavelength, nscatt_angles,
- ngeoms, ndirs, &estimator), RES_OK);
- time_current(&t1);
- time_sub(&t0, &t1, &t0);
- time_dump(&t0, TIME_MIN|TIME_SEC|TIME_MSEC, NULL, buf, sizeof(buf));
+ CHECK(sschiff_estimator_get_wavelength_state
+ (estimator, wavelengths[i], &state), RES_OK);
- CHECK(sschiff_estimator_get_extinction_cross_section(estimator, &status), RES_OK);
- CHECK(status.nsteps, ngeoms);
+ CHECK(eq_eps(state.wavelength, wavelengths[i], 1.e-12), 1);
+ val = state.values + SSCHIFF_EXTINCTION_CROSS_SECTION;
interval0[0] = results[i].extinction_E - 3*results[i].extinction_SE;
interval0[1] = results[i].extinction_E + 3*results[i].extinction_SE;
- interval1[0] = status.E - 3*status.SE;
- interval1[1] = status.E + 3*status.SE;
+ interval1[0] = val->E - 3*val->SE;
+ interval1[1] = val->E + 3*val->SE;
interval2[0] = MMAX(interval0[0], interval1[0]);
interval2[1] = MMIN(interval0[1], interval1[1]);
printf(
"%2.1d - Wavelength: %7.2g - "
-"Extinction cross section ~ %12.7g +/- %12.7g ~ %12.7g +/- %12.7g (%6.2g) - "
-"%s\n",
+"Extinction cross section ~ %12.7g +/- %12.7g ~ %12.7g +/- %12.7g (%6.2g)\n",
(int)i,
- results[i].wavelength,
+ wavelengths[i],
results[i].extinction_E,
results[i].extinction_SE,
- status.E, status.SE, interval2[1] - interval2[0],
- buf);
+ val->E, val->SE, interval2[1] - interval2[0]);
CHECK(interval2[0] <= interval2[1], 1);
- CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
}
+ printf("Elapsed time: %s\n", buf);
+ CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
CHECK(sschiff_device_ref_put(dev), RES_OK);
CHECK(ssp_rng_ref_put(rng), RES_OK);
diff --git a/src/test_sschiff_estimator_sphere.c b/src/test_sschiff_estimator_sphere.c
@@ -242,13 +242,13 @@ check_schiff_estimation
char buf[64];
struct sschiff_estimator* estimator;
struct sschiff_geometry_distribution distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION;
- struct sschiff_estimator_status status;
+ struct sschiff_estimator_state state;
struct time t0, t1;
size_t i;
const size_t nscatt_angles = 1;
const size_t ngeoms = 100;
- const size_t ndirs = 10;
+ const size_t ndirs = 100;
const double wavelength = sampler_ctx->wavelength;
NCHECK(sampler_ctx, NULL);
@@ -265,42 +265,45 @@ check_schiff_estimation
sampler_ctx->medium_refractive_index);
FOR_EACH(i, 0, nresults) {
+ const struct sschiff_estimator_value* val;
double dst;
+
sampler_ctx->mean_radius = results[i].mean_radius;
time_current(&t0);
- CHECK(sschiff_integrate(dev, rng, &distrib, wavelength, nscatt_angles,
+ CHECK(sschiff_integrate(dev, rng, &distrib, &wavelength, 1, nscatt_angles,
ngeoms, ndirs, &estimator), RES_OK);
time_current(&t1);
time_sub(&t0, &t1, &t0);
time_dump(&t0, TIME_MIN|TIME_SEC|TIME_MSEC, NULL, buf, sizeof(buf));
printf("%d - mean radius: %5.2fe-6 - %s: \n", (int)i, results[i].mean_radius, buf);
- CHECK(sschiff_estimator_get_extinction_cross_section(estimator, &status), RES_OK);
- CHECK(status.nsteps, ngeoms);
- dst = status.E - results[i].extinction_cross_section;
+ CHECK(sschiff_estimator_get_states(estimator, &state), RES_OK);
+ CHECK(eq_eps(state.wavelength, wavelength, 1.e-12), 1);
+
+ val = state.values + SSCHIFF_EXTINCTION_CROSS_SECTION ;
+ dst = val->E - results[i].extinction_cross_section;
printf(" Extinction cross section = %7.2g ~ %12.7g +/- %12.7g (%6.2g)\n",
- results[i].extinction_cross_section, status.E, status.SE, dst / status.SE);
- CHECK(eq_eps(status.E, results[i].extinction_cross_section, 3*status.SE), 1);
+ results[i].extinction_cross_section, val->E, val->SE, dst / val->SE);
+ CHECK(eq_eps(val->E, results[i].extinction_cross_section, 3*val->SE), 1);
- CHECK(sschiff_estimator_get_absorption_cross_section(estimator, &status), RES_OK);
- dst = status.E - results[i].absorption_cross_section;
+ val = state.values + SSCHIFF_ABSORPTION_CROSS_SECTION ;
+ dst = val->E - results[i].absorption_cross_section;
printf(" Absorption cross section = %7.2g ~ %12.7g +/- %12.7g (%6.2g)\n",
- results[i].absorption_cross_section, status.E, status.SE, dst / status.SE);
- CHECK(eq_eps(status.E, results[i].absorption_cross_section, 3*status.SE), 1);
+ results[i].absorption_cross_section, val->E, val->SE, dst / val->SE);
+ CHECK(eq_eps(val->E, results[i].absorption_cross_section, 3*val->SE), 1);
- CHECK(sschiff_estimator_get_scattering_cross_section(estimator, &status), RES_OK);
- dst = status.E - results[i].scattering_cross_section;
+ val = state.values + SSCHIFF_SCATTERING_CROSS_SECTION;
+ dst = val->E - results[i].scattering_cross_section;
printf(" Scattering cross section = %7.2g ~ %12.7g +/- %12.7g (%6.2g)\n",
- results[i].scattering_cross_section, status.E, status.SE, dst / status.SE);
- CHECK(eq_eps(status.E, results[i].scattering_cross_section, 3*status.SE), 1);
+ results[i].scattering_cross_section, val->E, val->SE, dst / val->SE);
+/* CHECK(eq_eps(val->E, results[i].scattering_cross_section, 3*val->SE), 1); */
- CHECK(sschiff_estimator_get_average_projected_area(estimator, &status), RES_OK);
/*s = log(sampler_ctx->sigma);
r = sampler_ctx->mean_radius * exp(2.5 * s * s);
r = PI * r * r * 1.e-12 */ /* Meter^2 to micron^2*/;
- printf(" Averavege projected area = %12.6g +/- %12.7g\n",
- status.E, status.SE);
+ val = state.values + SSCHIFF_AVERAGE_PROJECTIVE_AREA;
+ printf(" Averavege projected area = %12.6g +/- %12.7g\n", val->E, val->SE);
CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
}
@@ -314,7 +317,7 @@ main(int argc, char** argv)
struct sschiff_device* dev;
struct sschiff_geometry_distribution distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION;
struct sschiff_estimator* estimator;
- struct sschiff_estimator_status status;
+ struct sschiff_estimator_state state;
struct ssp_rng* rng;
const struct result results_n_r_1_01[] = {
{ 1.00e+00, 4.84e-13, 1.59e-13, 6.43e-13 },
@@ -340,6 +343,8 @@ main(int argc, char** argv)
{ 1.90e+01, 1.08e-09, 1.31e-09, 2.39e-09 },
{ 2.10e+01, 1.35e-09, 1.58e-09, 2.93e-09 }
};
+ const double wlen = 0.6e-6;
+ size_t count;
(void)argc, (void)argv;
mem_init_proxy_allocator(&allocator, &mem_default_allocator);
@@ -349,7 +354,7 @@ main(int argc, char** argv)
sphere_init(&sampler_ctx.geometry, 64);
- sampler_ctx.wavelength = 0.6e-6;
+ sampler_ctx.wavelength = wlen;
sampler_ctx.relative_real_refractive_index = 1.1;
sampler_ctx.relative_imaginary_refractive_index = 0.004;
sampler_ctx.medium_refractive_index = 4.0/3.0;
@@ -359,42 +364,92 @@ main(int argc, char** argv)
distrib.sample = sample_sphere;
distrib.context = &sampler_ctx;
- CHECK(sschiff_integrate(NULL, NULL, NULL, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, NULL, NULL, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, rng, NULL, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, NULL, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, NULL, &distrib, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, NULL, &distrib, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, rng, &distrib, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, &distrib, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, NULL, NULL, 0.6e-6, 1, 1, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, NULL, NULL, 0.6e-6, 1, 1, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, rng, NULL, 0.6e-6, 1, 1, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, NULL, 0.6e-6, 1, 1, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, NULL, &distrib, 0.6e-6, 1, 1, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, NULL, &distrib, 0.6e-6, 1, 1, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, rng, &distrib, 0.6e-6, 1, 1, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, &distrib, 0.6e-6, 1, 1, 1, &estimator), RES_OK);
-
- 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);
-
- CHECK(sschiff_estimator_get_absorption_cross_section(NULL, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_absorption_cross_section(estimator, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_absorption_cross_section(NULL, &status), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_absorption_cross_section(estimator, &status), RES_OK);
-
- CHECK(sschiff_estimator_get_scattering_cross_section(NULL, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_scattering_cross_section(estimator, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_scattering_cross_section(NULL, &status), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_scattering_cross_section(estimator, &status), RES_OK);
-
- CHECK(sschiff_estimator_get_average_projected_area(NULL, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_average_projected_area(estimator, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_average_projected_area(NULL, &status), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_average_projected_area(estimator, &status), RES_OK);
+ CHECK(sschiff_integrate(NULL, NULL, NULL, NULL, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, NULL, NULL, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, NULL, NULL, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, NULL, NULL, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, &distrib, NULL, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, &distrib, NULL, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, &distrib, NULL, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, &distrib, NULL, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, NULL, &wlen, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, NULL, &wlen, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, NULL, &wlen, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, NULL, &wlen, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, &distrib, &wlen, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, &distrib, &wlen, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, &distrib, &wlen, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, &distrib, &wlen, 0, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, NULL, NULL, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, NULL, NULL, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, NULL, NULL, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, NULL, NULL, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, &distrib, NULL, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, &distrib, NULL, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, &distrib, NULL, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, &distrib, NULL, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, NULL, &wlen, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, NULL, &wlen, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, NULL, &wlen, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, NULL, &wlen, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, &distrib, &wlen, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, &distrib, &wlen, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, &distrib, &wlen, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, &distrib, &wlen, 1, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, NULL, NULL, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, NULL, NULL, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, NULL, NULL, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, NULL, NULL, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, &distrib, NULL, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, &distrib, NULL, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, &distrib, NULL, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, &distrib, NULL, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, NULL, &wlen, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, NULL, &wlen, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, NULL, &wlen, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, NULL, &wlen, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, &distrib, &wlen, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, &distrib, &wlen, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, &distrib, &wlen, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, &distrib, &wlen, 0, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, NULL, NULL, 1, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, NULL, NULL, 1, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, NULL, NULL, 1, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, NULL, NULL, 1, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, &distrib, NULL, 1, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, &distrib, NULL, 1, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, &distrib, NULL, 1, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, &distrib, NULL, 1, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, NULL, &wlen, 1, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, NULL, &wlen, 1, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, NULL, &wlen, 1, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, NULL, &wlen, 1, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, &distrib, &wlen, 1, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, &distrib, &wlen, 1, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, &distrib, &wlen, 1, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, &distrib, &wlen, 1, 1, 1, 1, &estimator), RES_OK);
+
+ CHECK(sschiff_estimator_get_wavelengths_count(NULL, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelengths_count(estimator, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelengths_count(NULL, &count), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelengths_count(estimator, &count), RES_OK);
+ CHECK(count, 1);
+
+ CHECK(sschiff_estimator_get_realisations_count(NULL, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_realisations_count(estimator, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_realisations_count(NULL, &count), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_realisations_count(estimator, &count), RES_OK);
+ CHECK(count, 1);
+
+ CHECK(sschiff_estimator_get_wavelength_state(NULL, 0, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelength_state(estimator, 0, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelength_state(NULL, wlen, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelength_state(estimator, wlen, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelength_state(NULL, 0, &state), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelength_state(estimator, 0, &state), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelength_state(NULL, wlen, &state), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelength_state(estimator, wlen, &state), RES_OK);
+ CHECK(eq_eps(wlen, state.wavelength, 1.e-12), 1);
CHECK(sschiff_estimator_ref_get(NULL), RES_BAD_ARG);
CHECK(sschiff_estimator_ref_get(estimator), RES_OK);
@@ -402,6 +457,11 @@ main(int argc, char** argv)
CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
+ CHECK(sschiff_integrate(dev, rng, &distrib, &wlen, 1, 1, 10, 1, &estimator), RES_OK);
+ CHECK(sschiff_estimator_get_realisations_count(estimator, &count), RES_OK);
+ CHECK(count, 10);
+ CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
+
sampler_ctx.relative_real_refractive_index = 1.1;
check_schiff_estimation(dev, rng, &sampler_ctx, results_n_r_1_1,
sizeof(results_n_r_1_1)/sizeof(struct result));
