commit 40b3d9d737008dc8497b135421b0f48894061112
parent 134ae95a64866f8b5668b022635b4e5e60906c78
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 29 Oct 2015 15:35:41 +0100
Refactor the integration function
Make the main integration function more agnostic of the geometric
distribution.
Diffstat:
8 files changed, 317 insertions(+), 210 deletions(-)
diff --git a/src/schiff.c b/src/schiff.c
@@ -28,7 +28,111 @@
#include "schiff_args.h"
#include "schiff_sphere.h"
+#include "schiff_optical_properties.h"
+#include <rsys/stretchy_array.h>
+
+#include <star/sschiff.h>
+#include <star/ssp.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+run(const struct schiff_args* args)
+{
+ FILE* fp = stdout;
+ struct sschiff_geometry_distribution distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION;
+ struct sschiff_estimator* estimator = NULL;
+ struct sschiff_estimator_state* states = NULL;
+ struct sschiff_device* sschiff = NULL;
+ struct ssp_rng* rng = NULL;
+ size_t iwlen, nwlens;
+ res_T res = RES_OK;
+ ASSERT(args);
+
+ if(args->output_filename) {
+ fp = fopen(args->output_filename, "w");
+ if(!fp) {
+ fprintf(stderr,
+ "Couldn't open the output file `%s'.", args->output_filename);
+ res = RES_IO_ERR;
+ goto error;
+ }
+ }
+
+ res = ssp_rng_create(&mem_default_allocator, &ssp_rng_threefry, &rng);
+ if(res != RES_OK) {
+ fprintf(stderr, "Couldn't create the Random Number Generator.\n");
+ goto error;
+ }
+
+ res = sschiff_device_create(NULL, NULL, 0, NULL, &sschiff);
+ if(res != RES_OK) {
+ fprintf(stderr, "Couldn't create the Star Schiff device.\n");
+ goto error;
+ }
+
+ res = schiff_geometry_distribution_sphere_init
+ (&distrib, &args->geom.sphere, args->properties);
+ if(res != RES_OK) {
+ fprintf(stderr,
+ "Couldn't initialise the Star Schiff geometry distribution.\n");
+ goto error;
+ }
+
+ /* Invoke the Schiff integration */
+ nwlens = sa_size(args->wavelengths);
+ res = sschiff_integrate(sschiff, rng, &distrib, args->wavelengths, nwlens, 1,
+ args->ngeoms, args->ndirs, &estimator);
+ schiff_geometry_distribution_sphere_release(&distrib);
+
+ if(res != RES_OK) {
+ fprintf(stderr, "Schiff integration error.\n");
+ goto error;
+ }
+
+ /* Retrieve estimation results */
+ states = sa_add(states, nwlens);
+ SSCHIFF(estimator_get_states(estimator, states));
+
+ fprintf(fp,
+"# Line format \"W E E' A A' S S' P P'\" with \"W\" the wavelength in meter,\n"
+"# \"E\", \"A\", and \"S\" the estimation of the extinction, absorption and\n"
+"# scattering cross section, respectively, and \"P\" the estimation of the\n"
+"# average projected area. The `prime' values, i.e. \"E'\", \"A'\", \"S'\" and \"P'\"\n"
+"# are the standard error of the aforementionned estimations.\n\n");
+
+ /* Print the estimation results */
+ FOR_EACH(iwlen, 0, nwlens) {
+ const struct sschiff_estimator_value* val;
+
+ fprintf(fp, "%g ", states[iwlen].wavelength);
+
+ val = states[iwlen].values + SSCHIFF_EXTINCTION_CROSS_SECTION;
+ fprintf(fp, "%g %g ", val->E, val->SE);
+ val = states[iwlen].values + SSCHIFF_ABSORPTION_CROSS_SECTION;
+ fprintf(fp, "%g %g ", val->E, val->SE);
+ val = states[iwlen].values + SSCHIFF_SCATTERING_CROSS_SECTION;
+ fprintf(fp, "%g %g ", val->E, val->SE);
+ val = states[iwlen].values + SSCHIFF_AVERAGE_PROJECTED_AREA;
+ fprintf(fp, "%g %g\n", val->E, val->SE);
+ }
+
+exit:
+ if(fp && fp != stdout) fclose(fp);
+ if(estimator) SSCHIFF(estimator_ref_put(estimator));
+ if(sschiff) SSCHIFF(device_ref_put(sschiff));
+ if(rng) SSP(rng_ref_put(rng));
+ sa_release(states);
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Entry point
+ ******************************************************************************/
int
main(int argc, char** argv)
{
@@ -39,11 +143,14 @@ main(int argc, char** argv)
res = schiff_args_init(&args, argc, argv);
if(res != RES_OK) goto error;
if(!args.wavelengths) goto exit;
- res = schiff_run_sphere(&args);
+ res = run(&args);
if(res != RES_OK) goto error;
exit:
schiff_args_release(&args);
+ if(mem_allocated_size() != 0) {
+ fprintf(stderr, "Memory leaks: %lu Bytes\n", mem_allocated_size());
+ }
return err;
error:
err = 1;
diff --git a/src/schiff.h b/src/schiff.h
@@ -1,39 +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. */
-
-#ifndef SCHIFF_H
-#define SCHIFF_H
-
-#include <star/sbox_cmd.h>
-
-struct sbox_schiff {
- struct mem_allocator allocator;
-};
-
-#endif /* SBOX_SCHIFF_H */
-
diff --git a/src/schiff_args.c b/src/schiff_args.c
@@ -58,6 +58,14 @@ print_help(const char* binary)
"the relative refractive index, and \"Ne\" the refractive index of the medium.\n"
"With no FILE, read optical properties from standard input.\n\n");
printf(
+" -c R:S:A[:N] the micro organisms are cylindrical meshes discretized in N\n"
+" slices. By default N is %u. The radius `r' and the height `h'\n"
+" of the cylinders are distributed as follow:\n"
+" y = 1/(ln(S)*x*sqrt(2PI)) * exp(-(ln(x)-ln(R))^2 / (2*ln(S)^2))\n"
+" r = y / (3 / (2*A))^1/3\n"
+" h = 2*r/A\n",
+ SCHIFF_CYLINDER_DEFAULT.nslices);
+ printf(
" -d DIRS number of sampled directions for each geometry. Default is %u.\n",
SCHIFF_ARGS_NULL.ndirs);
printf(
@@ -69,17 +77,63 @@ print_help(const char* binary)
" stdout.\n");
printf(
" -s R:S[:N] the micro organisms are spherical meshes discretized in N\n"
-" slices along 2PI. Their radius is distributed with respect to a\n"
-" lognormal distribution:\n"
-" 1/(ln(S)*x*sqrt(2PI)) * exp(-(ln(x)-ln(R))^2 / (2*ln(S)^2))\n"
-" By default N is %u.\n",
- SCHIFF_DISTRIBUTION_SPHERE_DEFAULT.nslices);
+" slices along 2PI. By default N is %u. Their radius `r' is\n"
+" distributed with respect to a lognormal distribution:\n"
+" r = 1/(ln(S)*x*sqrt(2PI)) * exp(-(ln(x)-ln(R))^2 / (2*ln(S)^2))\n",
+ SCHIFF_SPHERE_DEFAULT.nslices);
printf(
" -w A[:B]... list of wavelengths (in micron) to integrate. At least one\n"
" wavelength must be provided.\n");
}
static res_T
+parse_cylinder_distribution(const char* str, struct schiff_args* args)
+{
+ double list[4];
+ size_t len;
+ unsigned nslices;
+ res_T res = RES_OK;
+ ASSERT(args && str);
+
+ /* The distribution was already set to a non cylinder distribution. */
+ if(args->geom_type!= SCHIFF_NONE && args->geom_type != SCHIFF_CYLINDER) {
+ res = RES_BAD_ARG;
+ goto error;
+
+ }
+ res = cstr_to_list_double(str, list, &len, 4);
+ if(res != RES_OK) goto error; /* Wrong argument formating */
+
+ /* At least the radius the sigma and the aspect_ratio must be provided */
+ if(len < 3) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(len == 3) {
+ nslices = SCHIFF_CYLINDER_DEFAULT.nslices;
+ } else {
+ nslices = (unsigned)list[3];
+ if((double)nslices != list[3]) { /* The #slices arg is not an integer */
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+
+ /* Update the distribution */
+ args->geom_type = SCHIFF_CYLINDER;
+ args->geom.cylinder.radius = list[0];
+ args->geom.cylinder.sigma = list[1];
+ args->geom.cylinder.aspect_ratio = list[2];
+ args->geom.cylinder.nslices = nslices;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
parse_sphere_distribution(const char* str, struct schiff_args* args)
{
double list[3];
@@ -89,7 +143,7 @@ parse_sphere_distribution(const char* str, struct schiff_args* args)
ASSERT(args && str);
/* The distribution was already set to a non sphere distribution. */
- if(args->distrib_type != SCHIFF_NONE && args->distrib_type != SCHIFF_SPHERE) {
+ if(args->geom_type != SCHIFF_NONE && args->geom_type != SCHIFF_SPHERE) {
res = RES_BAD_ARG;
goto error;
@@ -103,7 +157,7 @@ parse_sphere_distribution(const char* str, struct schiff_args* args)
}
if(len == 2) {
- nslices = SCHIFF_DISTRIBUTION_SPHERE_DEFAULT.nslices;
+ nslices = SCHIFF_SPHERE_DEFAULT.nslices;
} else {
nslices = (unsigned)list[2];
if((double)nslices != list[2]) { /* The #slices arg is not an integer */
@@ -113,10 +167,10 @@ parse_sphere_distribution(const char* str, struct schiff_args* args)
}
/* Update the distribution */
- args->distrib_type = SCHIFF_SPHERE;
- args->distrib.sphere.radius = list[0];
- args->distrib.sphere.sigma = list[1];
- args->distrib.sphere.nslices = nslices;
+ args->geom_type = SCHIFF_SPHERE;
+ args->geom.sphere.radius = list[0];
+ args->geom.sphere.sigma = list[1];
+ args->geom.sphere.nslices = nslices;
exit:
return res;
@@ -157,14 +211,6 @@ error:
goto exit;
}
-static int
-cmp_properties(const void* op0, const void* op1)
-{
- const struct schiff_optical_properties* prop0 = op0;
- const struct schiff_optical_properties* prop1 = op1;
- return prop0->W < prop1->W ? -1 : (prop0->W > prop1->W ? 1 : 0);
-}
-
/*******************************************************************************
* Local function
******************************************************************************/
@@ -179,9 +225,10 @@ schiff_args_init
ASSERT(argc && argv && args);
*args = SCHIFF_ARGS_NULL;
- while((opt = getopt(argc, argv, "d:g:ho:s:w:")) != -1) {
+ while((opt = getopt(argc, argv, "c:d:g:ho:s:w:")) != -1) {
res_T res = RES_OK;
switch(opt) {
+ case 'c': res = parse_cylinder_distribution(optarg, args); break;
case 'd': res = cstr_to_uint(optarg, &args->ndirs); break;
case 'g': res = cstr_to_uint(optarg, &args->ngeoms); break;
case 'h':
@@ -220,10 +267,6 @@ schiff_args_init
res = RES_BAD_ARG;
goto error;
}
- /* Sort the optical properties in ascending order with respect to the
- * wavelength */
- qsort(args->properties, sa_size(args->properties),
- sizeof(struct schiff_optical_properties), cmp_properties);
exit:
return res;
diff --git a/src/schiff_args.h b/src/schiff_args.h
@@ -31,49 +31,48 @@
#include <rsys/rsys.h>
-/* List of supported distribution */
-enum schiff_distribution_type {
+/* List of supported geometry */
+enum schiff_type {
SCHIFF_CYLINDER,
SCHIFF_SPHERE,
SCHIFF_NONE
};
-/* Sphere distribution */
-struct schiff_distribution_sphere {
+struct schiff_sphere {
double radius;
double sigma;
unsigned nslices;
};
-#define SCHIFF_DISTRIBUTION_SPHERE_DEFAULT__ {1.0, 1.0, 64}
-static const struct schiff_distribution_sphere
-SCHIFF_DISTRIBUTION_SPHERE_DEFAULT = SCHIFF_DISTRIBUTION_SPHERE_DEFAULT__;
+#define SCHIFF_SPHERE_DEFAULT__ {1.0, 1.0, 64}
+static const struct schiff_sphere SCHIFF_SPHERE_DEFAULT =
+ SCHIFF_SPHERE_DEFAULT__;
-/* Cylinder distribution */
-struct schiff_distribution_cylinder {
+struct schiff_cylinder {
double radius;
+ double sigma;
double aspect_ratio;
unsigned nslices;
};
-#define SCHIFF_DISTRIBUTION_CYLINDER_DEFAULT__ {1.0, 1.0, 64}
-static const struct schiff_distribution_cylinder
-SCHIFF_DISTRIBUTION_CYLINDER_DEFAULT = SCHIFF_DISTRIBUTION_SPHERE_DEFAULT__;
+#define SCHIFF_CYLINDER_DEFAULT__ {1.0, 1.0, 1.0, 64}
+static const struct schiff_cylinder SCHIFF_CYLINDER_DEFAULT =
+ SCHIFF_CYLINDER_DEFAULT__;
/* Generic distribution */
-union schiff_distribution {
- struct schiff_distribution_sphere sphere;
- struct schiff_distribution_cylinder cylinder;
+union schiff_geometry {
+ struct schiff_sphere sphere;
+ struct schiff_cylinder cylinder;
};
-#define SCHIFF_DISTRIBUTION_DEFAULT__ {SCHIFF_DISTRIBUTION_SPHERE_DEFAULT__}
+#define SCHIFF_GEOMETRY_DEFAULT__ {SCHIFF_SPHERE_DEFAULT__}
struct schiff_args {
const char* output_filename;
struct schiff_optical_properties* properties;
double* wavelengths;
- enum schiff_distribution_type distrib_type;
- union schiff_distribution distrib;
+ enum schiff_type geom_type;
+ union schiff_geometry geom;
unsigned ngeoms;
unsigned ndirs;
};
@@ -83,7 +82,7 @@ static const struct schiff_args SCHIFF_ARGS_NULL = {
NULL, /* List of optical properties */
NULL, /* List of wavelength to integrate */
SCHIFF_NONE,
- { SCHIFF_DISTRIBUTION_SPHERE_DEFAULT__ },
+ SCHIFF_GEOMETRY_DEFAULT__,
1000, /* # Sampled geometries */
100, /* # Sampled directions per gemetry */
};
diff --git a/src/schiff_optical_properties.c b/src/schiff_optical_properties.c
@@ -31,6 +31,8 @@
#include <rsys/cstr.h>
#include <rsys/stretchy_array.h>
+#include <star/sschiff.h>
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -85,6 +87,14 @@ parse_optical_properties
return RES_OK;
}
+static int
+cmp_properties(const void* op0, const void* op1)
+{
+ const struct schiff_optical_properties* prop0 = op0;
+ const struct schiff_optical_properties* prop1 = op1;
+ return prop0->W < prop1->W ? -1 : (prop0->W > prop1->W ? 1 : 0);
+}
+
/*******************************************************************************
* Local function
******************************************************************************/
@@ -159,8 +169,50 @@ schiff_optical_properties_load_stream
}
}
+ /* Sort the optical properties in ascending order with respect to the
+ * wavelength */
+ qsort(props, sa_size(props),
+ sizeof(struct schiff_optical_properties), cmp_properties);
+
*out_props = props;
sa_release(buf);
return RES_OK;
}
+void
+schiff_optical_properties_fetch
+ (struct schiff_optical_properties* properties, /* Stetchy array */
+ const double wavelength,
+ struct sschiff_material_properties* props)
+{
+ const size_t nproperties = sa_size(properties);
+ double Ne, Kr, Nr;
+ size_t i;
+ ASSERT(properties && props && wavelength > 0.0);
+
+ /* Assume that properties are sorted in ascending order with respect to the
+ * wavelength. TODO use a binary search of the upper bound */
+ FOR_EACH(i, 0, nproperties) {
+ if(properties[i].W >= wavelength) break;
+ }
+
+ if(eq_eps(properties[i].W, wavelength, 1.e-12)
+ || i==0 || i==nproperties) { /* Clamp to wavelength */
+ Ne = properties[i].Ne;
+ Kr = properties[i].Kr;
+ Nr = properties[i].Nr;
+ } else { /* Linear interpolation of optical properties */
+ 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;
+ 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;
+ }
+
+ props->medium_refractive_index = Ne;
+ props->relative_imaginary_refractive_index = Kr;
+ props->relative_real_refractive_index = Nr;
+}
+
diff --git a/src/schiff_optical_properties.h b/src/schiff_optical_properties.h
@@ -32,6 +32,8 @@
#include <rsys/dynamic_array.h>
#include <rsys/rsys.h>
+struct sschiff_material_properties;
+
struct schiff_optical_properties {
double W; /* Wavelength */
double Nr; /* Real part of the relative refractive index */
@@ -61,5 +63,11 @@ schiff_optical_properties_load_stream
FILE* stream,
const char* stream_name);
+extern LOCAL_SYM void
+schiff_optical_properties_fetch
+ (struct schiff_optical_properties* properties, /* Stetchy array */
+ const double wavelength,
+ struct sschiff_material_properties* props);
+
#endif /* SCHIFF_OPTICAL_PROPERTIES_H */
diff --git a/src/schiff_sphere.c b/src/schiff_sphere.c
@@ -35,9 +35,7 @@
#include <star/ssp.h>
#include <star/sschiff.h>
-#include <rsys/clock_time.h>
#include <rsys/float3.h>
-#include <rsys/stretchy_array.h>
/* Sphere geometry */
struct sphere {
@@ -46,7 +44,7 @@ struct sphere {
};
/* Geometry distribution of a sphere */
-struct schiff_sphere_geometry_distribution_context {
+struct sphere_geometry_distribution_context {
struct schiff_mesh* mesh; /* Triangular mesh of the template mesh */
struct schiff_optical_properties* properties; /* Per wavelength properties */
double log_mean_radius; /* Log of the mean radius */
@@ -57,14 +55,14 @@ struct schiff_sphere_geometry_distribution_context {
* Helper functions
******************************************************************************/
static void
-get_indices(const unsigned itri, unsigned ids[3], void* ctx)
+sphere_get_indices(const unsigned itri, unsigned ids[3], void* ctx)
{
struct sphere* sphere = ctx;
schiff_mesh_get_indices(sphere->mesh, itri, ids);
}
static void
-get_position(const unsigned ivert, float vertex[3], void* ctx)
+sphere_get_position(const unsigned ivert, float vertex[3], void* ctx)
{
struct sphere* sphere = ctx;
schiff_mesh_get_position(sphere->mesh, ivert, vertex);
@@ -72,52 +70,23 @@ get_position(const unsigned ivert, float vertex[3], void* ctx)
}
static void
-get_material_property
+sphere_get_material_property
(void* mtl,
const double wavelength,
struct sschiff_material_properties* props)
{
struct schiff_optical_properties* properties = mtl;
- const size_t nproperties = sa_size(properties);
- double Ne, Kr, Nr;
- size_t i;
-
- /* Assume that properties are sorted in ascending order with respect to the
- * wavelength. TODO use a binary search of the upper bound */
- FOR_EACH(i, 0, nproperties) {
- if(properties[i].W >= wavelength)
- break;
- }
-
- if(eq_eps(properties[i].W, wavelength, 1.e-12)
- || i==0 || i==nproperties) { /* Clamp to wavelength */
- Ne = properties[i].Ne;
- Kr = properties[i].Kr;
- Nr = properties[i].Nr;
- } else { /* Linear interpolation of optical properties */
- 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;
- 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);
- props->medium_refractive_index = Ne;
- props->relative_imaginary_refractive_index = Kr;
- props->relative_real_refractive_index = Nr;
+ schiff_optical_properties_fetch(properties, wavelength, props);
}
static res_T
-schiff_sphere_sample_geometry
+sphere_sample_geometry
(struct ssp_rng* rng,
struct sschiff_material* mtl,
struct s3d_shape* shape,
void* context)
{
- struct schiff_sphere_geometry_distribution_context* ctx = context;
+ struct sphere_geometry_distribution_context* ctx = context;
struct sphere sphere;
struct s3d_vertex_data attrib;
size_t nverts, nprims;
@@ -129,124 +98,82 @@ schiff_sphere_sample_geometry
attrib.usage = S3D_POSITION;
attrib.type = S3D_FLOAT3;
- attrib.get = get_position;
+ attrib.get = sphere_get_position;
- mtl->get_property = get_material_property;
+ mtl->get_property = sphere_get_material_property;
mtl->material = ctx->properties;
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);
+ sphere_get_indices, (unsigned)nverts, &attrib, 1, &sphere);
}
/*******************************************************************************
* Local functions
******************************************************************************/
-
res_T
-schiff_run_sphere
- (const struct schiff_args* args)
+schiff_geometry_distribution_sphere_init
+ (struct sschiff_geometry_distribution* distrib,
+ const struct schiff_sphere* sphere,
+ struct schiff_optical_properties* properties)
{
- FILE* fp = stdout;
- struct sschiff_geometry_distribution distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION;
- struct schiff_sphere_geometry_distribution_context ctx;
- struct sschiff_estimator* estimator = NULL;
- struct sschiff_estimator_state* states = NULL;
- struct sschiff_device* sschiff = NULL;
- struct ssp_rng* rng = NULL;
- struct schiff_mesh mesh;
- size_t iwlen, nwlens;
+ struct sphere_geometry_distribution_context* ctx = NULL;
res_T res = RES_OK;
- ASSERT(args);
-
- if(args->output_filename) {
- fp = fopen(args->output_filename, "w");
- if(!fp) {
- fprintf(stderr,
- "Couldn't open the output file `%s'.", args->output_filename);
- res = RES_IO_ERR;
- goto error;
- }
- }
+ ASSERT(sphere);
- /* Generate the mesh template of the distribution */
- res = schiff_mesh_init_sphere
- (&mem_default_allocator, &mesh, args->distrib.sphere.nslices);
- if(res != RES_OK) {
- fprintf(stderr, "Couldn't create the mesh template of the distribution.\n");
+ ctx = mem_calloc(1, sizeof(struct sphere_geometry_distribution_context));
+ if(!ctx) {
+ fprintf(stderr,
+"Not enough memory: couldn't allocate the sphere geometry distribution context\n");
+ res = RES_MEM_ERR;
goto error;
}
-
- /* TODO this can be parametrised */
- res = ssp_rng_create(&mem_default_allocator, &ssp_rng_threefry, &rng);
- if(res != RES_OK) {
- fprintf(stderr, "Couldn't create the Random Number Generator.\n");
+ ctx->mesh = mem_calloc(1, sizeof(struct schiff_mesh));
+ if(!ctx->mesh) {
+ fprintf(stderr,
+"Not enough memory: couldn't allocate the sphere mesh template.\n");
+ res = RES_MEM_ERR;
goto error;
}
- res = sschiff_device_create(NULL, NULL, 0, NULL, &sschiff);
+ /* Generate the mesh template of the distribution */
+ res = schiff_mesh_init_sphere
+ (&mem_default_allocator, ctx->mesh, sphere->nslices);
if(res != RES_OK) {
- fprintf(stderr, "Couldn't create the Star Schiff device.\n");
+ fprintf(stderr, "Couldn't initialise the sphere mesh template.\n");
goto error;
}
/* Setup the geometry distribution context */
- ctx.mesh = &mesh;
- ctx.properties = args->properties;
- ctx.log_mean_radius = log(args->distrib.sphere.radius);
- ctx.log_sigma = log(args->distrib.sphere.sigma);
-
- /* Setup the geometry distribution */
- distrib.sample = schiff_sphere_sample_geometry;
- distrib.context = &ctx;
-
- /* Invoke the Schiff integration */
- nwlens = sa_size(args->wavelengths);
- res = sschiff_integrate(sschiff, rng, &distrib, args->wavelengths, nwlens, 1,
- args->ngeoms, args->ndirs, &estimator);
- if(res != RES_OK) {
- fprintf(stderr, "Schiff integration error.\n");
- goto error;
- }
+ ctx->properties = properties;
+ ctx->log_mean_radius = log(sphere->radius);
+ ctx->log_sigma = log(sphere->sigma);
- /* Retrieve estimation results */
- states = sa_add(states, nwlens);
- SSCHIFF(estimator_get_states(estimator, states));
-
- fprintf(fp,
-"# Line format \"W E E' A A' S S' P P'\" with \"W\" the wavelength in meter,\n"
-"# \"E\", \"A\", and \"S\" the estimation of the extinction, absorption and\n"
-"# scattering cross section, respectively, and \"P\" the estimation of the\n"
-"# average projected area. The `prime' values, i.e. \"E'\", \"A'\", \"S'\" and \"P'\"\n"
-"# are the standard error of the aforementionned estimations.\n\n");
-
- /* Print the estimation results */
- FOR_EACH(iwlen, 0, nwlens) {
- const struct sschiff_estimator_value* val;
-
- fprintf(fp, "%g ", states[iwlen].wavelength);
-
- val = states[iwlen].values + SSCHIFF_EXTINCTION_CROSS_SECTION;
- fprintf(fp, "%g %g ", val->E, val->SE);
- val = states[iwlen].values + SSCHIFF_ABSORPTION_CROSS_SECTION;
- fprintf(fp, "%g %g ", val->E, val->SE);
- val = states[iwlen].values + SSCHIFF_SCATTERING_CROSS_SECTION;
- fprintf(fp, "%g %g ", val->E, val->SE);
- val = states[iwlen].values + SSCHIFF_AVERAGE_PROJECTED_AREA;
- fprintf(fp, "%g %g\n", val->E, val->SE);
- }
+ distrib->sample = sphere_sample_geometry;
+ distrib->context = ctx;
exit:
- if(fp && fp != stdout) fclose(fp);
- if(estimator) SSCHIFF(estimator_ref_put(estimator));
- if(sschiff) SSCHIFF(device_ref_put(sschiff));
- if(rng) SSP(rng_ref_put(rng));
- sa_release(states);
- schiff_mesh_release(&mesh);
return res;
error:
+ if(ctx) {
+ if(ctx->mesh) mem_rm(ctx->mesh);
+ mem_rm(ctx);
+ ctx = NULL;
+ }
goto exit;
}
+void
+schiff_geometry_distribution_sphere_release
+ (struct sschiff_geometry_distribution* distrib)
+{
+ struct sphere_geometry_distribution_context* ctx;
+ ASSERT(distrib);
+ ctx = distrib->context;
+ schiff_mesh_release(ctx->mesh);
+ mem_rm(ctx->mesh);
+ mem_rm(ctx);
+}
+
diff --git a/src/schiff_sphere.h b/src/schiff_sphere.h
@@ -31,8 +31,18 @@
#include <rsys/rsys.h>
+struct schiff_distribution_sphere;
+struct schiff_optical_properties;
+struct sschiff_geometry_distribution;
+
extern LOCAL_SYM res_T
-schiff_run_sphere
- (const struct schiff_args* args);
+schiff_geometry_distribution_sphere_init
+ (struct sschiff_geometry_distribution* distrib,
+ const struct schiff_sphere* sphere,
+ struct schiff_optical_properties* properties);
+
+extern LOCAL_SYM void
+schiff_geometry_distribution_sphere_release
+ (struct sschiff_geometry_distribution* distrib);
#endif /* SCHIFF_SPHERE_H */
