commit 1279c65ed70567198a98237457e91319247a1951
parent d95149b1c8a679f6726611c39c987bc42def30e8
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 23 Oct 2015 16:18:09 +0200
Refactor the Schiff mesh generation
Replace the ad-hock mesh data structure by the sbox_schiff_mesh
structure. Store mesh data in dynamic array rather than stretchy arrays
in order to handle memory allocation errors.
Diffstat:
4 files changed, 333 insertions(+), 151 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -64,12 +64,14 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SBOX_SCHIFF_FILES_SRC
sbox_schiff.c
- sbox_schiff_sphere.c
- sbox_schiff_optical_properties.c)
+ sbox_schiff_mesh.c
+ sbox_schiff_optical_properties.c
+ sbox_schiff_sphere.c)
set(SBOX_SCHIFF_FILES_INC
sbox_schiff.h
- sbox_schiff_sphere.h
- sbox_schiff_optical_properties.h)
+ sbox_schiff_mesh.h
+ sbox_schiff_optical_properties.h
+ sbox_schiff_sphere.h)
set(SBOX_SCHIFF_FILES_DOC COPYING.fr COPYING.en README.md)
# Prepend each file in the `SBOX_SCHIFF_FILES_<SRC|INC>' list by the absolute
diff --git a/src/sbox_schiff_mesh.c b/src/sbox_schiff_mesh.c
@@ -0,0 +1,191 @@
+/* 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 "sbox_schiff_mesh.h"
+#include <rsys/float3.h>
+
+struct sin_cos { double sinus, cosine; };
+
+#define DARRAY_NAME sincos
+#define DARRAY_DATA struct sin_cos
+#include <rsys/dynamic_array.h>
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+sbox_schiff_mesh_init_sphere
+ (struct mem_allocator* allocator,
+ struct sbox_schiff_mesh* sphere,
+ const unsigned nthetas)
+{
+ /* Theta in [0, 2PI[ and phi in [0, PI[ */
+ const unsigned nphis = (unsigned)(((double)nthetas + 0.5) / 2.0);
+ const double step_theta = 2*PI / (double)nthetas;
+ const double step_phi = PI / (double)nphis;
+ struct darray_sincos sincos_thetas;
+ struct darray_sincos sincos_phis;
+ size_t nverts;
+ size_t ntris;
+ size_t i;
+ unsigned itheta, iphi;
+ res_T res = RES_OK;
+ ASSERT(allocator && sphere && nthetas);
+
+ darray_float_init(allocator, &sphere->vertices);
+ darray_uint_init(allocator, &sphere->indices);
+ darray_sincos_init(allocator, &sincos_thetas);
+ darray_sincos_init(allocator, &sincos_phis);
+
+ nverts = nthetas*(nphis-1)/* #contour verts */ + 2/* polar verts */;
+ ntris = 2*nthetas*(nphis-2)/* #contour tris */ + 2*nthetas/* #polar tris */;
+
+ res = darray_float_resize(&sphere->vertices, nverts * 3/*#coords per vert*/);
+ if(res != RES_OK) goto error;
+ res = darray_uint_resize(&sphere->indices, ntris * 3/*#indices per tri*/);
+ if(res != RES_OK) goto error;
+ res = darray_sincos_resize(&sincos_thetas, nthetas);
+ if(res != RES_OK) goto error;
+ res = darray_sincos_resize(&sincos_phis, nphis);
+ if(res != RES_OK) goto error;
+
+ /* Precompute the cosine/sinus of the theta/phi angles */
+ FOR_EACH(itheta, 0, nthetas) {
+ const double theta = -PI + (double)itheta * step_theta;
+ darray_sincos_data_get(&sincos_thetas)[itheta].sinus = sin(theta);
+ darray_sincos_data_get(&sincos_thetas)[itheta].cosine = cos(theta);
+ }
+ FOR_EACH(iphi, 1, nphis) {
+ const double phi = -PI/2 + (double)iphi * step_phi;
+ darray_sincos_data_get(&sincos_phis)[iphi].sinus = sin(phi);
+ darray_sincos_data_get(&sincos_phis)[iphi].cosine = cos(phi);
+ }
+
+ /* Build the contour vertices */
+ i = 0;
+ FOR_EACH(itheta, 0, nthetas) {
+ const struct sin_cos* theta = darray_sincos_data_get(&sincos_thetas) + itheta;
+ FOR_EACH(iphi, 0, nphis-1) {
+ const struct sin_cos* phi = darray_sincos_data_get(&sincos_phis) + iphi;
+ darray_float_data_get(&sphere->vertices)[i++] = (float)(theta->cosine * phi->cosine);
+ darray_float_data_get(&sphere->vertices)[i++] = (float)(theta->sinus * phi->cosine);
+ darray_float_data_get(&sphere->vertices)[i++] = (float)phi->sinus;
+ }
+ }
+ /* Setup polar vertices */
+ f3(darray_float_data_get(&sphere->vertices) + i + 0, 0.f, 0.f,-1.f);
+ f3(darray_float_data_get(&sphere->vertices) + i + 3, 0.f, 0.f, 1.f);
+
+ /* Define the indices of the contour primitives */
+ i = 0;
+ FOR_EACH(itheta, 0, nthetas) {
+ const unsigned itheta0 = itheta * (nphis - 1);
+ const unsigned itheta1 = ((itheta + 1) % nthetas) * (nphis - 1);
+ FOR_EACH(iphi, 0, nphis-2) {
+ const unsigned iphi0 = iphi + 0;
+ const unsigned iphi1 = iphi + 1;
+
+ darray_uint_data_get(&sphere->indices)[i++] = itheta0 + iphi0;
+ darray_uint_data_get(&sphere->indices)[i++] = itheta0 + iphi1;
+ darray_uint_data_get(&sphere->indices)[i++] = itheta1 + iphi0;
+
+ darray_uint_data_get(&sphere->indices)[i++] = itheta1 + iphi0;
+ darray_uint_data_get(&sphere->indices)[i++] = itheta0 + iphi1;
+ darray_uint_data_get(&sphere->indices)[i++] = itheta1 + iphi1;
+ }
+ }
+
+ /* Define the indices of the polar primitives */
+ FOR_EACH(itheta, 0, nthetas) {
+ const unsigned itheta0 = itheta * (nphis - 1);
+ const unsigned itheta1 = ((itheta + 1) % nthetas) * (nphis - 1);
+
+ darray_uint_data_get(&sphere->indices)[i++] = nthetas * (nphis - 1);
+ darray_uint_data_get(&sphere->indices)[i++] = itheta0;
+ darray_uint_data_get(&sphere->indices)[i++] = itheta1;
+
+ darray_uint_data_get(&sphere->indices)[i++] = nthetas * (nphis - 1) + 1;
+ darray_uint_data_get(&sphere->indices)[i++] = itheta1 + (nphis - 2);
+ darray_uint_data_get(&sphere->indices)[i++] = itheta0 + (nphis - 2);
+ }
+
+exit:
+ darray_sincos_release(&sincos_thetas);
+ darray_sincos_release(&sincos_phis);
+ return res;
+error:
+ sbox_schiff_mesh_release(sphere);
+ goto exit;
+}
+
+void
+sbox_schiff_mesh_release(struct sbox_schiff_mesh* mesh)
+{
+ ASSERT(mesh);
+ darray_float_release(&mesh->vertices);
+ darray_uint_release(&mesh->indices);
+}
+
+res_T
+sbox_schiff_mesh_dump(struct sbox_schiff_mesh* mesh, const char* filename)
+{
+ FILE* fp = NULL;
+ size_t nverts, ntris;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(mesh && filename);
+ ASSERT(darray_float_size_get(&mesh->vertices) % 3 == 0); /* 3D */
+ ASSERT(darray_uint_size_get(&mesh->indices) % 3 == 0); /* Triangles */
+
+ fp = fopen(filename, "r");
+ if(!fp) {
+ res = RES_IO_ERR;
+ goto error;
+ }
+
+ nverts = darray_float_size_get(&mesh->vertices) / 3;
+ FOR_EACH(i, 0, nverts) {
+ fprintf(fp, "v %f %f %f\n",
+ SPLIT3(darray_float_data_get(&mesh->vertices) + i*3));
+ }
+
+ ntris = darray_uint_size_get(&mesh->indices) / 3;
+ FOR_EACH(i, 0, ntris) {
+ fprintf(fp, "f %u %u %u\n",
+ darray_uint_data_get(&mesh->indices)[i*3 + 0],
+ darray_uint_data_get(&mesh->indices)[i*3 + 1],
+ darray_uint_data_get(&mesh->indices)[i*3 + 2]);
+ }
+
+exit:
+ if(fp) fclose(fp);
+ return res;
+error:
+ goto exit;
+}
+
diff --git a/src/sbox_schiff_mesh.h b/src/sbox_schiff_mesh.h
@@ -0,0 +1,56 @@
+/* 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. */
+
+#ifndef SBOX_SCHIFF_MESH_H
+#define SBOX_SCHIFF_MESH_H
+
+#include <rsys/dynamic_array_float.h>
+#include <rsys/dynamic_array_uint.h>
+
+struct sbox_schiff_mesh {
+ struct darray_float vertices;
+ struct darray_uint indices;
+};
+
+extern LOCAL_SYM res_T
+sbox_schiff_mesh_init_sphere
+ (struct mem_allocator* alocator,
+ struct sbox_schiff_mesh* mesh,
+ const unsigned nthetas); /* # discret points along 2PI */
+
+extern LOCAL_SYM void
+sbox_schiff_mesh_release
+ (struct sbox_schiff_mesh* mesh);
+
+extern LOCAL_SYM res_T
+sbox_schiff_mesh_dump
+ (struct sbox_schiff_mesh* mesh,
+ const char* filename);
+
+#endif /* SBOX_SCHIFF_MESH_H */
+
diff --git a/src/sbox_schiff_sphere.c b/src/sbox_schiff_sphere.c
@@ -29,6 +29,7 @@
#define _POSIX_C_SOURCE 2 /* getopt support */
#include "sbox_schiff.h"
+#include "sbox_schiff_mesh.h"
#include "sbox_schiff_optical_properties.h"
#include "sbox_schiff_sphere.h"
@@ -47,7 +48,6 @@
#include <unistd.h>
#endif
-
/* Default values */
#define NDIRS 100
#define NGEOMS 1000
@@ -55,6 +55,7 @@
#define SIGMA 1.0
#define NTHETAS 64
+/* Argument of the schiff sphere command */
struct schiff_sphere_args {
const char* material;
const char* output;
@@ -71,23 +72,18 @@ static const struct schiff_sphere_args SCHIFF_SPHERE_ARGS_NULL = {
NULL, NULL, NULL, NULL, RADIUS, SIGMA, NDIRS, NGEOMS, NTHETAS
};
-struct mesh {
- unsigned* indices;
- float* vertices;
-};
-
-static const struct mesh MESH_NULL = { NULL, NULL };
-
+/* Sphere geometry */
struct sphere {
- struct mesh* mesh;
- float radius;
+ struct sbox_schiff_mesh* mesh; /* Triangular mesh of an unit sphere */
+ float radius; /* Sphere radius */
};
+/* Geometry distribution of a sphere */
struct schiff_sphere_geometry_distribution_context {
- struct mesh* mesh;
- struct sbox_schiff_optical_properties* properties;
- double mean_radius;
- double sigma;
+ struct sbox_schiff_mesh* mesh; /* Triangular mesh of an unit sphere */
+ struct sbox_schiff_optical_properties* properties; /* Per wavelength properties */
+ double mean_radius; /* Mean radius of the log normal sphere distribution */
+ double sigma; /* Sigma argument of the log normal sphere distribution */
};
/*******************************************************************************
@@ -105,9 +101,7 @@ cmp_properties(const void* op0, const void* op1)
{
const struct sbox_schiff_optical_properties* prop0 = op0;
const struct sbox_schiff_optical_properties* prop1 = op1;
- if (prop0->W < prop1->W) return -1;
- else if(prop0->W > prop1->W) return 1;
- else return 0;
+ return prop0->W < prop1->W ? -1 : (prop0->W > prop1->W ? 1 : 0);
}
static void
@@ -115,23 +109,22 @@ get_indices(const unsigned itri, unsigned ids[3], void* ctx)
{
struct sphere* sphere = ctx;
const size_t i = itri * 3;
- ASSERT(sa_size(sphere->mesh->indices) % 3 == 0);
- ASSERT(itri < sa_size(sphere->mesh->indices) / 3);
- ids[0] = sphere->mesh->indices[i + 0];
- ids[1] = sphere->mesh->indices[i + 1];
- ids[2] = sphere->mesh->indices[i + 2];
+ ASSERT(darray_uint_size_get(&sphere->mesh->indices) % 3 == 0);
+ ASSERT(itri < darray_uint_size_get(&sphere->mesh->indices) / 3);
+ ids[0] = darray_uint_data_get(&sphere->mesh->indices)[i + 0];
+ ids[1] = darray_uint_data_get(&sphere->mesh->indices)[i + 1];
+ ids[2] = darray_uint_data_get(&sphere->mesh->indices)[i + 2];
}
-static INLINE void
+static void
get_position(const unsigned ivert, float vertex[3], void* ctx)
{
struct sphere* sphere = ctx;
const size_t i = ivert * 3;
- ASSERT(sa_size(sphere->mesh->vertices) % 3 == 0);
- ASSERT(ivert < sa_size(sphere->mesh->vertices) / 3);
- vertex[0] = sphere->mesh->vertices[i + 0] * sphere->radius;
- vertex[1] = sphere->mesh->vertices[i + 1] * sphere->radius;
- vertex[2] = sphere->mesh->vertices[i + 2] * sphere->radius;
+ ASSERT(darray_float_size_get(&sphere->mesh->vertices) % 3 == 0);
+ ASSERT(ivert < darray_float_size_get(&sphere->mesh->vertices) / 3);
+ f3_set(vertex, darray_float_data_get(&sphere->mesh->vertices) + i);
+ f3_mulf(vertex, vertex, sphere->radius);
}
static void
@@ -146,7 +139,7 @@ get_material_property
size_t i;
/* Assume that properties are sorted in ascending order with respect to the
- * wavelength. TODO use a binary search */
+ * wavelength. TODO use a binary search of the upper bound */
FOR_EACH(i, 0, nproperties) {
if(properties[i].W >= wavelength)
break;
@@ -161,9 +154,10 @@ get_material_property
const size_t i_prev = i - 1;
const double d = properties[i].W - properties[i_prev].W;
const double u = (wavelength - properties[i_prev].W) / d;
- Ne = u * properties[i].Ne + (1-u) * properties[i_prev].Ne;
- Kr = u * properties[i].Kr + (1-u) * properties[i_prev].Kr;
- Nr = u * properties[i].Nr + (1-u) * properties[i_prev].Nr;
+ const double v = 1.0 - u;
+ Ne = u * properties[i].Ne + v * properties[i_prev].Ne;
+ Kr = u * properties[i].Kr + v * properties[i_prev].Kr;
+ Nr = u * properties[i].Nr + v * properties[i_prev].Nr;
}
NCHECK(props, NULL);
@@ -172,93 +166,6 @@ get_material_property
props->relative_real_refractive_index = Nr;
}
-static void
-mesh_init_sphere(struct mesh* 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;
- ASSERT(sphere && ntheta);
-
- 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
-mesh_release(struct mesh* sphere)
-{
- sa_release(sphere->vertices);
- sa_release(sphere->indices);
- sphere->vertices = NULL;
- sphere->indices = NULL;
-}
-
/* Assume that props0 and props1 are sorted. If two entries have the same
* wavelength the entry of props1 is skipped. */
static void
@@ -413,8 +320,8 @@ schiff_sphere_sample_geometry
mtl->get_property = get_material_property;
mtl->material = ctx->properties;
- nverts = sa_size(ctx->mesh->vertices) / 3/*#coords*/;
- nprims = sa_size(ctx->mesh->indices) / 3/*#indices per prim*/;
+ nverts = darray_float_size_get(&ctx->mesh->vertices) / 3/*#coords*/;
+ nprims = darray_uint_size_get(&ctx->mesh->indices) / 3/*#indices per prim*/;
return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims,
get_indices, (unsigned)nverts, &attrib, 1, &sphere);
@@ -426,16 +333,26 @@ schiff_sphere_run
struct sbox_schiff* schiff,
struct schiff_sphere_args* args)
{
+ char buf[64];
struct sschiff_geometry_distribution distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION;
struct schiff_sphere_geometry_distribution_context ctx;
- struct mesh mesh = MESH_NULL;
+ struct sschiff_estimator* estimator = NULL;
+ struct sschiff_estimator_state* states = NULL;
struct sschiff_device* sschiff = NULL;
struct ssp_rng* rng = NULL;
- size_t iwave, nwaves;
+ struct sbox_schiff_mesh mesh;
+ struct time t0, t1;
+ size_t iwlen, nwlens;
res_T res = RES_OK;
ASSERT(sbox);
- mesh_init_sphere(&mesh, args->nthetas);
+ /* Generate the spherical mesh */
+ res = sbox_schiff_mesh_init_sphere(&schiff->allocator, &mesh, args->nthetas);
+ if(res != RES_OK) {
+ logger_print(&sbox->logger, LOG_ERROR,
+ "Couldn't create the Schiff sphere mesh discretised in %u steps along 2PI\n",
+ args->nthetas);
+ }
/* FIXME use the rng type defined in sbox */
res = ssp_rng_create(&schiff->allocator, &ssp_rng_threefry, &rng);
@@ -453,47 +370,63 @@ schiff_sphere_run
goto error;
}
+ /* Setup the geometry distribution context */
ctx.mesh = &mesh;
ctx.properties = args->props;
ctx.mean_radius = args->radius;
ctx.sigma = args->sigma;
+ /* Setup the geometry distribution */
distrib.sample = schiff_sphere_sample_geometry;
distrib.context = &ctx;
- nwaves = sa_size(args->wavelengths);
- FOR_EACH(iwave, 0, nwaves) {
- struct sschiff_estimator* estimator;
- struct sschiff_estimator_status status;
- struct time t0, t1;
- char buf[64];
-
- time_current(&t0);
- res = sschiff_integrate(sschiff, rng, &distrib, args->wavelengths[iwave], 1,
+ /* Invoke the Schiff integration */
+ nwlens = sa_size(args->wavelengths);
+ time_current(&t0);
+ res = sschiff_integrate(sschiff, rng, &distrib, args->wavelengths, nwlens, 1,
args->ngeoms, args->ndirs, &estimator);
- time_current(&t1);
- if(res != RES_OK) {
- logger_print(&sbox->logger, LOG_ERROR,
- "Error in Schiff integration of the wavelength `%g'.\n",
- args->wavelengths[iwave]);
- continue;
- }
+ time_current(&t1);
+ if(res != RES_OK) {
+ logger_print(&sbox->logger, LOG_ERROR, "Schiff integration error.\n");
+ goto error;
+ }
+ time_sub(&t0, &t1, &t0);
+ time_dump(&t0, TIME_MIN|TIME_SEC|TIME_MSEC, NULL, buf, sizeof(buf));
+
+ /* Retrieve estimation results */
+ states = sa_add(states, nwlens);
+ SSCHIFF(estimator_get_states(estimator, states));
+
+ /* Print the estimation results */
+ FOR_EACH(iwlen, 0, nwlens) {
+ const struct sschiff_estimator_value* val;
+
+ logger_print(&sbox->logger, LOG_OUTPUT,
+ "Wavelength %7.2g:\n", states[iwlen].wavelength);
- time_sub(&t0, &t1, &t0);
- time_dump(&t0, TIME_MIN|TIME_SEC|TIME_MSEC, NULL, buf, sizeof(buf));
+ val = states[iwlen].values + SSCHIFF_EXTINCTION_CROSS_SECTION;
+ logger_print(&sbox->logger, LOG_OUTPUT,
+ " Extinction cross section ~ %12.7g +/- %12.7g\n", val->E, val->SE);
- SSCHIFF(estimator_get_extinction_cross_section(estimator, &status));
+ val = states[iwlen].values + SSCHIFF_ABSORPTION_CROSS_SECTION;
logger_print(&sbox->logger, LOG_OUTPUT,
-"Wavelength: %7.2g - Extinction cross section ~ %12.7g +/- %12.7g - %s\n",
- args->wavelengths[iwave], status.E, status.SE, buf);
+ " Absorption cross section ~ %12.7g +/- %12.7g\n", val->E, val->SE);
- SSCHIFF(estimator_ref_put(estimator));
+ val = states[iwlen].values + SSCHIFF_SCATTERING_CROSS_SECTION;
+ logger_print(&sbox->logger, LOG_OUTPUT,
+ " Scattering cross section ~ %12.7g +/- %12.7g\n", val->E, val->SE);
+
+ val = states[iwlen].values + SSCHIFF_AVERAGE_PROJECTED_AREA;
+ logger_print(&sbox->logger, LOG_OUTPUT,
+ " Average projected area ~ %12.7g +/- %12.7g\n", val->E, val->SE);
}
exit:
+ if(estimator) SSCHIFF(estimator_ref_put(estimator));
if(sschiff) SSCHIFF(device_ref_put(sschiff));
if(rng) SSP(rng_ref_put(rng));
- mesh_release(&mesh);
+ sa_release(states);
+ sbox_schiff_mesh_release(&mesh);
return res;
error:
goto exit;
