commit eeee76acf3d587adfdcadd36c08444def80b684a
parent a57cf548c99f5ea0af24d4df231b45f3d14d2df9
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 13 Oct 2015 12:15:42 +0200
Refactor the Schiff estimator test
Improve the performances of the geometry sampler
Diffstat:
4 files changed, 389 insertions(+), 396 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -93,7 +93,7 @@ if(NOT NO_TEST)
add_test(${_name} ${_name})
endfunction(new_test)
new_test(test_sschiff_device)
- new_test(test_sschiff_estimator)
+ new_test(test_sschiff_estimator_sphere)
endif()
################################################################################
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -147,7 +147,7 @@ compute_path_length
if(S3D_HIT_NONE(&hit)) return RES_BAD_ARG;
- len += hit.distance;
+ len += hit.distance - range[0];
range[0] = hit.distance;
prim_from = hit.prim;
@@ -161,7 +161,7 @@ compute_path_length
} while(!S3D_HIT_NONE(&hit));
- *length = (len - first_hit->distance) * 1.e-6/* Micron to meter */;
+ *length = len * 1.e-6/* Meter to micron */;
return RES_OK;
}
diff --git a/src/test_sschiff_estimator.c b/src/test_sschiff_estimator.c
@@ -1,393 +0,0 @@
-/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
- *
- * This software is governed by the CeCILL license under French law and
- * abiding by the rules of distribution of free software. You can use,
- * modify and/or redistribute the software under the terms of the CeCILL
- * license as circulated by CEA, CNRS and INRIA at the following URL
- * "http://www.cecill.info".
- *
- * As a counterpart to the access to the source code and rights to copy,
- * modify and redistribute granted by the license, users are provided only
- * with a limited warranty and the software's author, the holder of the
- * economic rights, and the successive licensors have only limited
- * liability.
- *
- * In this respect, the user's attention is drawn to the risks associated
- * with loading, using, modifying and/or developing or reproducing the
- * software by the user in light of its specific status of free software,
- * that may mean that it is complicated to manipulate, and that also
- * therefore means that it is reserved for developers and experienced
- * professionals having in-depth computer knowledge. Users are therefore
- * encouraged to load and test the software's suitability as regards their
- * requirements in conditions enabling the security of their systems and/or
- * data to be ensured and, more generally, to use and operate it in the
- * same conditions as regards security.
- *
- * The fact that you are presently reading this means that you have had
- * knowledge of the CeCILL license and that you accept its terms. */
-
-#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.1f
-#define NREALISATIONS 1000
-
-struct optical_properties { double m, n_r, kappa_r, n_medium; };
-
-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
-get_material_property
- (void* mtl,
- const double wavelength,
- struct sschiff_material_properties* properties)
-{
- 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 / (sphere->nphi - 1);
- const size_t j = iquad % (sphere->nphi - 1);
- size_t A, B;
-
- 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;
- ids[2] = (unsigned)B;
- } else { ASSERT(iquad_tri == 1);
- ids[0] = (unsigned)B;
- ids[1] = (unsigned)A;
- ids[2] = (unsigned)A + 1;
- }
-}
-
-static INLINE void
-sphere_get_vertex_position(const unsigned ivert, float vertex[3], void* ctx)
-{
- 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];
-
- if(iphi == sphere->nphi - 1) { /* Polar vertices: Force polar coordinates */
- angle[0] = PI;
- angle[1] = PI/2;
- } else {
- angle[0] = -PI + (double)itheta * STEP;
- angle[1] = -PI/2 + (double)iphi * STEP;
- }
-
- 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)(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_sphere
- (struct ssp_rng* rng,
- struct sschiff_material* mtl,
- struct s3d_shape* shape,
- void* sampler_context)
-{
- struct s3d_vertex_data attrib;
- 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*/;
-
- 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);
- NCHECK(shape, NULL);
-
- attrib.usage = S3D_POSITION;
- attrib.type = S3D_FLOAT3;
- attrib.get = sphere_get_vertex_position;
-
- mtl->get_property = get_material_property;
- mtl->material = ctx;
-
- return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims,
- 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);
-
- 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_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);
- CHECK(sschiff_integrate(NULL, rng, NULL, 0, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, NULL, 0, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, NULL, &desc, 0, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, NULL, &desc, 0, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, rng, &desc, 0, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, &desc, 0, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, NULL, NULL, 0, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, NULL, NULL, 0, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, rng, NULL, 0, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, NULL, 0, 1, &estimator), RES_BAD_ARG);
- 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);
-
- 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);
- CHECK(sschiff_estimator_ref_put(NULL), RES_BAD_ARG);
- 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);
-
- check_memory_allocator(&allocator);
- mem_shutdown_proxy_allocator(&allocator);
- CHECK(mem_allocated_size(), 0);
- return 0;
-}
-
diff --git a/src/test_sschiff_estimator_sphere.c b/src/test_sschiff_estimator_sphere.c
@@ -0,0 +1,386 @@
+/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
+ *
+ * This software is governed by the CeCILL license under French law and
+ * abiding by the rules of distribution of free software. You can use,
+ * modify and/or redistribute the software under the terms of the CeCILL
+ * license as circulated by CEA, CNRS and INRIA at the following URL
+ * "http://www.cecill.info".
+ *
+ * As a counterpart to the access to the source code and rights to copy,
+ * modify and redistribute granted by the license, users are provided only
+ * with a limited warranty and the software's author, the holder of the
+ * economic rights, and the successive licensors have only limited
+ * liability.
+ *
+ * In this respect, the user's attention is drawn to the risks associated
+ * with loading, using, modifying and/or developing or reproducing the
+ * software by the user in light of its specific status of free software,
+ * that may mean that it is complicated to manipulate, and that also
+ * therefore means that it is reserved for developers and experienced
+ * professionals having in-depth computer knowledge. Users are therefore
+ * encouraged to load and test the software's suitability as regards their
+ * requirements in conditions enabling the security of their systems and/or
+ * data to be ensured and, more generally, to use and operate it in the
+ * same conditions as regards security.
+ *
+ * The fact that you are presently reading this means that you have had
+ * knowledge of the CeCILL license and that you accept its terms. */
+
+#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 NSTEPS 1000
+
+struct result {
+ double mean_radius;
+ double extinction_cross_section;
+};
+
+struct sphere {
+ float* vertices; /* List of float[3] */
+ unsigned* indices; /* List of unsigned[3] */
+};
+
+struct sampler_context {
+ struct sphere unit_sphere;
+ double wavelength;
+ double medium_refractive_index;
+ double relative_real_refractive_index;
+ double relative_imaginary_refractive_index;
+ double mean_radius;
+};
+
+struct sphere_context {
+ struct sphere* unit_sphere;
+ float radius;
+};
+
+static void
+get_material_property
+ (void* mtl,
+ const double wavelength,
+ struct sschiff_material_properties* properties)
+{
+ struct sampler_context* ctx = mtl;
+
+ NCHECK(mtl, NULL);
+ NCHECK(properties, NULL);
+ CHECK(eq_eps(ctx->wavelength, wavelength, 1.e-8), 1);
+
+ properties->medium_refractive_index = ctx->medium_refractive_index;
+ properties->relative_real_refractive_index =
+ ctx->relative_real_refractive_index;
+ properties->relative_imaginary_refractive_index =
+ ctx->relative_imaginary_refractive_index;
+}
+
+static void
+get_indices(const unsigned itri, unsigned ids[3], void* ctx)
+{
+ struct sphere_context* sphere = ctx;
+ const size_t i = itri * 3;
+
+ CHECK(sa_size(sphere->unit_sphere->indices) % 3, 0);
+ CHECK(itri < sa_size(sphere->unit_sphere->indices) / 3, 1);
+ ids[0] = sphere->unit_sphere->indices[i + 0];
+ ids[1] = sphere->unit_sphere->indices[i + 1];
+ ids[2] = sphere->unit_sphere->indices[i + 2];
+}
+
+static INLINE void
+get_position(const unsigned ivert, float vertex[3], void* ctx)
+{
+ struct sphere_context* sphere = ctx;
+ const size_t i = ivert * 3;
+
+ CHECK(sa_size(sphere->unit_sphere->vertices) % 3, 0);
+ CHECK(ivert < sa_size(sphere->unit_sphere->vertices) / 3, 1);
+ vertex[0] = sphere->unit_sphere->vertices[i + 0] * sphere->radius;
+ vertex[1] = sphere->unit_sphere->vertices[i + 1] * sphere->radius;
+ vertex[2] = sphere->unit_sphere->vertices[i + 2] * sphere->radius;
+}
+
+static res_T
+sample_sphere
+ (struct ssp_rng* rng,
+ struct sschiff_material* mtl,
+ struct s3d_shape* shape,
+ void* sampler_context)
+{
+ struct s3d_vertex_data attrib;
+ struct sampler_context* sampler_ctx = sampler_context;
+ struct sphere_context sphere;
+ size_t nverts, nprims;
+
+ NCHECK(rng, NULL);
+ NCHECK(mtl, NULL);
+ NCHECK(shape, NULL);
+
+ sphere.unit_sphere = &sampler_ctx->unit_sphere;
+ sphere.radius = (float)ssp_ran_lognormal
+ (rng, log(sampler_ctx->mean_radius), log(1.18));
+
+ attrib.usage = S3D_POSITION;
+ attrib.type = S3D_FLOAT3;
+ attrib.get = get_position;
+
+ mtl->get_property = get_material_property;
+ mtl->material = sampler_ctx;
+
+ nverts = sa_size(sphere.unit_sphere->vertices) / 3/*#coords*/;
+ nprims = sa_size(sphere.unit_sphere->indices) / 3/*#indices per prim*/;
+
+ return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims,
+ get_indices, (unsigned)nverts, &attrib, 1, &sphere);
+}
+
+static void
+unit_sphere_init(struct sphere* sphere, const unsigned ntheta)
+{
+ const unsigned nphi = (unsigned)(((double)ntheta + 0.5) / 2.0);
+ const double step_theta = 2*PI / (double)ntheta;
+ const double step_phi = PI / (double)nphi;
+ double* cos_theta = NULL;
+ double* sin_theta = NULL;
+ double* cos_phi = NULL;
+ double* sin_phi = NULL;
+ unsigned itheta, iphi;
+
+ NCHECK(sphere, NULL);
+ NCHECK(ntheta, 0);
+
+ sphere->vertices = NULL;
+ sphere->indices = NULL;
+
+ /* Precompute the cosine/sinus of the theta/phi angles */
+ FOR_EACH(itheta, 0, ntheta) {
+ const double theta = -PI + (double)itheta * step_theta;
+ sa_push(cos_theta, cos(theta));
+ sa_push(sin_theta, sin(theta));
+ }
+ FOR_EACH(iphi, 1, nphi) {
+ const double phi = -PI/2 + (double)iphi * step_phi;
+ sa_push(cos_phi, cos(phi));
+ sa_push(sin_phi, sin(phi));
+ }
+ /* Compute the contour vertices */
+ FOR_EACH(itheta, 0, ntheta) {
+ FOR_EACH(iphi, 0, nphi-1) {
+ sa_push(sphere->vertices, (float)(cos_theta[itheta]*cos_phi[iphi]));
+ sa_push(sphere->vertices, (float)(sin_theta[itheta]*cos_phi[iphi]));
+ sa_push(sphere->vertices, (float)sin_phi[iphi]);
+ }
+ }
+ /* Compute the polar vertices */
+ f3(sa_add(sphere->vertices, 3), 0.f, 0.f,-1.f);
+ f3(sa_add(sphere->vertices, 3), 0.f, 0.f, 1.f);
+
+ /* Define the indices of the contour primitives */
+ FOR_EACH(itheta, 0, ntheta) {
+ const unsigned itheta0 = itheta * (nphi - 1);
+ const unsigned itheta1 = ((itheta + 1) % ntheta) * (nphi - 1);
+ FOR_EACH(iphi, 0, nphi-2) {
+ const unsigned iphi0 = iphi + 0;
+ const unsigned iphi1 = iphi + 1;
+
+ sa_push(sphere->indices, itheta0 + iphi0);
+ sa_push(sphere->indices, itheta0 + iphi1);
+ sa_push(sphere->indices, itheta1 + iphi0);
+
+ sa_push(sphere->indices, itheta1 + iphi0);
+ sa_push(sphere->indices, itheta0 + iphi1);
+ sa_push(sphere->indices, itheta1 + iphi1);
+ }
+ }
+
+ /* Define the indices of the polar primitives */
+ FOR_EACH(itheta, 0, ntheta) {
+ const unsigned itheta0 = itheta * (nphi - 1);
+ const unsigned itheta1 = ((itheta + 1) % ntheta) * (nphi - 1);
+
+ sa_push(sphere->indices, ntheta * (nphi - 1));
+ sa_push(sphere->indices, itheta0);
+ sa_push(sphere->indices, itheta1);
+
+ sa_push(sphere->indices, ntheta * (nphi - 1) + 1);
+ sa_push(sphere->indices, itheta1 + (nphi - 2));
+ sa_push(sphere->indices, itheta0 + (nphi - 2));
+ }
+
+ /* Release the intermediary data structure */
+ sa_release(cos_theta);
+ sa_release(sin_theta);
+ sa_release(cos_phi);
+ sa_release(sin_phi);
+}
+
+static void
+unit_sphere_release(struct sphere* sphere)
+{
+ sa_release(sphere->vertices);
+ sa_release(sphere->indices);
+ sphere->vertices = NULL;
+ sphere->indices = NULL;
+}
+
+static void
+check_schiff_estimation
+ (struct sschiff_device* dev,
+ struct ssp_rng* rng,
+ struct sampler_context* sampler_ctx,
+ const struct result* results,
+ const size_t nresults)
+{
+ char buf[64];
+ struct sschiff_estimator* estimator;
+ struct sschiff_integrator_desc desc = SSCHIFF_NULL_INTEGRATOR_DESC;
+ struct sschiff_estimator_status status;
+ struct time t0, t1;
+ size_t i;
+
+ NCHECK(sampler_ctx, NULL);
+ NCHECK(results, NULL);
+ NCHECK(results, 0);
+
+ desc.sample_geometry = sample_sphere;
+ desc.scattering_angles_count = 1;
+ desc.sampler_context = sampler_ctx;
+
+ printf("Schiff sphere estimation - m: %g; n_r: %g; kappa_r: %g; n_e: %g\n",
+ sampler_ctx->wavelength,
+ sampler_ctx->relative_real_refractive_index,
+ sampler_ctx->relative_imaginary_refractive_index,
+ sampler_ctx->medium_refractive_index);
+
+ FOR_EACH(i, 0, nresults) {
+ sampler_ctx->mean_radius = results[i].mean_radius;
+
+ time_current(&t0);
+ CHECK(sschiff_integrate(dev, rng, &desc, 6.e-7, NSTEPS, &estimator), RES_OK);
+ time_current(&t1);
+
+ CHECK(sschiff_estimator_get_extinction_cross_section(estimator, &status), RES_OK);
+
+ 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 - "
+"Extinction cross section = %7.2g ~ %12.7g +/- %12.7g - "
+"%s\n",
+ (int)i,
+ results[i].mean_radius,
+ results[i].extinction_cross_section,
+ status.E, status.SE,
+ buf);
+
+ CHECK(eq_eps(status.E, results[i].extinction_cross_section, 2*status.SE), 1);
+ CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
+ }
+}
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sampler_context sampler_ctx;
+ struct sschiff_device* dev;
+ struct sschiff_integrator_desc desc = SSCHIFF_NULL_INTEGRATOR_DESC;
+ struct sschiff_estimator* estimator;
+ struct sschiff_estimator_status status;
+ struct ssp_rng* rng;
+ const struct result results_n_r_1_01[] = {
+ { 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 struct result results_n_r_1_1[] = {
+ { 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);
+
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng), RES_OK);
+ CHECK(sschiff_device_create(NULL, &allocator, 0, NULL, &dev), RES_OK);
+
+ unit_sphere_init(&sampler_ctx.unit_sphere, 64);
+
+ sampler_ctx.wavelength = 0.6e-6;
+ 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;
+ sampler_ctx.mean_radius = 1.0;
+
+ desc.sample_geometry = sample_sphere;
+ desc.scattering_angles_count = 1;
+ desc.sampler_context = &sampler_ctx;
+
+ CHECK(sschiff_integrate(NULL, NULL, NULL, 6e-7, 0, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, NULL, 6e-7, 0, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, NULL, 6e-7, 0, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, NULL, 6e-7, 0, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, &desc, 6e-7, 0, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, &desc, 6e-7, 0, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, &desc, 6e-7, 0, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, &desc, 6e-7, 0, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, NULL, 6e-7, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, NULL, 6e-7, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, NULL, 6e-7, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, NULL, 6e-7, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, &desc, 6e-7, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, &desc, 6e-7, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, &desc, 6e-7, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, &desc, 6e-7, 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_ref_get(NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_ref_get(estimator), RES_OK);
+ CHECK(sschiff_estimator_ref_put(NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
+ 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));
+
+ sampler_ctx.relative_real_refractive_index = 1.01;
+ check_schiff_estimation(dev, rng, &sampler_ctx, results_n_r_1_01,
+ sizeof(results_n_r_1_01)/sizeof(struct result));
+
+ CHECK(sschiff_device_ref_put(dev), RES_OK);
+ CHECK(ssp_rng_ref_put(rng), RES_OK);
+
+ unit_sphere_release(&sampler_ctx.unit_sphere);
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHECK(mem_allocated_size(), 0);
+ return 0;
+}
+
