commit d95149b1c8a679f6726611c39c987bc42def30e8
parent ce30fcf5c1a2aad59475e447d9ba0c245221f807
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 21 Oct 2015 14:56:59 +0200
First draft of the "schiff sphere" computation
Implement a sub-set of the "schiff sphere" command. Output file is not
handled yet.
Diffstat:
6 files changed, 685 insertions(+), 290 deletions(-)
diff --git a/README.md b/README.md
@@ -9,8 +9,10 @@ that estimates the radiative properties of micro organism with an
The library relies on the [CMake](http://www.cmake.org) and the
[RCMake](https://gitlab.com/vaplv/rcmake/) package to build. It also depends
-on the [RSys](https://gitlab.com/vaplv/rsys/) and the
-[Star-Schiff](https://gitlab.com/meso-star/star-schiff/) libraries and the
+on the [RSys](https://gitlab.com/vaplv/rsys/)
+[Star-3D](https://gitlab.com/meso-star/star-3d/),
+[Star-Schiff](https://gitlab.com/meso-star/star-schiff/),
+[Star-SP](https://gitlab.com/meso-star/star-sp/) libraries, and the
[Star-Box](https://gitlab.com/meso-schiff/star-box/) software development kit.
First ensure that CMake is installed on your system. Then install the RCMake
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -36,8 +36,10 @@ set(SBOX_SCHIFF_SOURCE_DIR ${PROJECT_SOURCE_DIR}/../src/)
################################################################################
find_package(RCMake 0.2 REQUIRED)
find_package(RSys 0.2.1 REQUIRED)
-find_package(StarSchiff 0.3 REQUIRED)
+find_package(Star3D 0.3 REQUIRED)
find_package(StarBox 0.3 REQUIRED)
+find_package(StarSchiff 0.3 REQUIRED)
+find_package(StarSP 0.3 REQUIRED)
if(MSVC)
find_package(MuslGetopt REQUIRED)
include_directories(${MuslGetopt_INCLUDE_DIR})
@@ -62,9 +64,11 @@ 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)
set(SBOX_SCHIFF_FILES_INC
sbox_schiff.h
+ sbox_schiff_sphere.h
sbox_schiff_optical_properties.h)
set(SBOX_SCHIFF_FILES_DOC COPYING.fr COPYING.en README.md)
@@ -79,7 +83,7 @@ if(CMAKE_COMPILER_IS_GNUCC)
endif()
add_library(sbox-schiff SHARED ${SBOX_SCHIFF_FILES_SRC} ${SBOX_SCHIFF_FILES_INC})
-target_link_libraries(sbox-schiff RSys StarBox StarSchiff ${MATH_LIB})
+target_link_libraries(sbox-schiff RSys Star3D StarBox StarSchiff StarSP ${MATH_LIB})
set_target_properties(sbox-schiff PROPERTIES
VERSION ${VERSION}
SOVERSION ${VERSION_MAJOR})
diff --git a/src/sbox_schiff.c b/src/sbox_schiff.c
@@ -26,194 +26,14 @@
* The fact that you are presently reading this means that you have had
* knowledge of the CeCILL license and that you accept its terms. */
-#define _POSIX_C_SOURCE 2 /* getopt support */
#include "sbox_schiff.h"
-#include "sbox_schiff_optical_properties.h"
+#include "sbox_schiff_sphere.h"
#include <star/sbox.h>
-#include <star/sbox_parse.h>
#include <rsys/mem_allocator.h>
#include <rsys/rsys.h>
-#include <rsys/stretchy_array.h>
-
-#ifdef COMPER_CL
- #include <getopt.h>
-#else
- #include <unistd.h>
-#endif
-
-#define NDIRS 100
-#define NGEOMS 1000
-#define RADIUS 1.0
-#define SIGMA 1.0
-#define NTHETAS 64
-
-struct schiff_sphere_args {
- const char* material;
- const char* output;
- struct sbox_schiff_optical_properties* props;
- double* wavelengths;
- double radius;
- double sigma;
- unsigned ndirs;
- unsigned ngeoms;
- unsigned nthetas;
-};
-
-static const struct schiff_sphere_args SCHIFF_SPHERE_ARGS_NULL = {
- NULL, NULL, NULL, NULL, RADIUS, SIGMA, NDIRS, NGEOMS, NTHETAS
-};
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static char
-check_wavelength(const double wavelength)
-{
- ASSERT(wavelength);
- return wavelength > 0.0;
-}
-
-static int
-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;
-}
-
-static res_T
-schiff_sphere_parse_args
- (struct sbox* sbox,
- const int argc,
- char** argv,
- struct schiff_sphere_args* args)
-{
- double props[4];
- double* wavelength;
- int opt;
- res_T res = RES_OK;
- ASSERT(sbox && argc && argv && args);
-
- *args = SCHIFF_SPHERE_ARGS_NULL;
- optind = 0;
- opterr = 0;
-
- while((opt = getopt(argc, argv, ":d:g:i:n:p:r:s:w:")) != -1) {
- res_T res = RES_OK;
- switch(opt) {
- case 'd': res = sbox_str_to_uint(optarg, &args->ndirs); break;
- case 'g': res = sbox_str_to_uint(optarg, &args->ngeoms); break;
- case 'i': args->material = optarg; break;
- case 'n': res = sbox_str_to_uint(optarg, &args->nthetas); break;
- case '0': args->output = optarg; break;
- case 'p':
- res = sbox_parse_list_double(optarg, props, 4);
- if(res == RES_OK) {
- struct sbox_schiff_optical_properties* opt_props = sa_add(args->props, 1);
- opt_props->W = props[0];
- opt_props->Nr = props[1];
- opt_props->Kr = props[2];
- opt_props->Ne = props[3];
- if(!sbox_schiff_optical_properties_check(opt_props))
- res = RES_BAD_ARG;
- }
- break;
- case 'r':
- res = sbox_str_to_double(optarg, &args->radius);
- if(res == RES_OK && args->radius < 0.0) res = RES_BAD_ARG;
- break;
- case 's':
- res = sbox_str_to_double(optarg, &args->sigma);
- if(res == RES_OK && args->sigma < 0.0) res = RES_BAD_ARG;
- break;
- case 'w':
- wavelength = sa_add(args->wavelengths, 1);
- res = sbox_str_to_double(optarg, wavelength);
- if(res == RES_OK && !check_wavelength(*wavelength)) res = RES_BAD_ARG;
- break;
- case ':':
- logger_print(&sbox->logger, LOG_ERROR,
- "%s: option requires an argument -- '%c'\n",
- argv[0], argv[optind-1][1]);
- res = RES_BAD_ARG;
- break;
- case '?':
- logger_print(&sbox->logger, LOG_ERROR,
- "%s: invalid option -- '%c'\n",
- argv[0], argv[optind-1][1]);
- res = RES_BAD_ARG;
- break;
- default: res = RES_BAD_ARG; break;
- }
- if(res != RES_OK) {
- if(optarg) {
- logger_print(&sbox->logger, LOG_ERROR,
- "%s: invalid option arguments '%s' -- '%c'\n",
- argv[0], optarg, opt);
- }
- goto error;
- }
- }
-exit:
- return res;
-error:
- sa_release(args->props);
- sa_release(args->wavelengths);
- *args = SCHIFF_SPHERE_ARGS_NULL;
- goto exit;
-}
-
-/* Assume that props0 and props1 are sorted. If two entries have the same
- * wavelength the entry of props1 is skipped. */
-static void
-merge_properties
- (struct sbox_schiff_optical_properties** props0,
- const struct sbox_schiff_optical_properties* props1)
-{
- const size_t sizeof_prop = sizeof(struct sbox_schiff_optical_properties);
- const struct sbox_schiff_optical_properties* src0;
- const struct sbox_schiff_optical_properties* src1;
- struct sbox_schiff_optical_properties* dst = NULL;
- size_t nsrc0, nsrc1;
- size_t isrc0 = 0, isrc1 =0;
- ASSERT(props0 && props1);
-
- src0 = *props0;
- src1 = props1;
- nsrc0 = sa_size(src0);
- nsrc1 = sa_size(src1);
-
- while(isrc0 < nsrc0 && isrc1 < nsrc1) {
- if(src0[isrc0].W < src1[isrc1].W) {
- sa_push(dst, src0[isrc0]);
- ++isrc0;
- } else if(src1[isrc1].W < src0[isrc0].W) {
- sa_push(dst, src1[isrc1]);
- ++isrc1;
- } else { /* src0[isrc0].W == src1[isrc1].W */
- sa_push(dst, src0[isrc0]);
- ++isrc0;
- ++isrc1;
- }
-
- if(isrc0 >= nsrc0 && isrc1 < nsrc1) {
- const size_t nremains = nsrc1 - isrc1;
- memcpy(sa_add(dst, nremains), src1 + isrc1, sizeof_prop * nremains);
- isrc1 = nsrc1;
- } else if(isrc1 >= nsrc1 && isrc0 < nsrc0) {
- const size_t nremains = nsrc0 - isrc0;
- memcpy(sa_add(dst, nremains), src0 + isrc0, sizeof_prop * nremains);
- isrc0 = nsrc0;
- }
- }
- sa_release(*props0);
- *props0 = dst;
-}
/*******************************************************************************
* Schiff command
@@ -232,71 +52,6 @@ cmd_func_schiff(struct sbox* sbox, const int argc, char** argv, void* ctx)
}
static res_T
-cmd_func_schiff_sphere(struct sbox* sbox, const int argc, char** argv, void* ctx)
-{
- struct schiff_sphere_args args = SCHIFF_SPHERE_ARGS_NULL;
- struct sbox_schiff_optical_properties* mtl_props = NULL;
- const size_t sizeof_prop = sizeof(struct sbox_schiff_optical_properties);
- size_t i;
- res_T res = RES_OK;
- (void)ctx;
-
- res = schiff_sphere_parse_args(sbox, argc, argv, &args);
- if(res != RES_OK) goto error;
-
- if(!args.wavelengths) {
- logger_print(&sbox->logger, LOG_ERROR,
- "Missing the wavelength argument(s).\n");
- res = RES_BAD_ARG;
- goto error;
- }
-
- /* Sort the submitted properties in ascending order with respect to the
- * wavelength */
- if(args.props) {
- qsort(args.props, sa_size(args.props), sizeof_prop, cmp_properties);
- }
- /* Load the material file and sort the its properties in ascending order with
- * respect to the wavelength */
- if(args.material) {
- res = sbox_schiff_optical_properties_load(sbox, args.material, &mtl_props);
- if(res != RES_OK) goto error;
- qsort(mtl_props, sa_size(mtl_props), sizeof_prop, cmp_properties);
- }
-
- /* Define the optical properties to use */
- if(mtl_props && !args.props) {
- args.props = mtl_props;
- mtl_props = NULL;
- } else if(args.props && mtl_props) {
- merge_properties(&args.props, mtl_props);
- sa_release(mtl_props);
- mtl_props = NULL;
- }
-
- if(!args.props) {
- logger_print(&sbox->logger, LOG_ERROR,
- "Radiative properties are not defined.\n");
- res = RES_BAD_ARG;
- goto error;
- }
-
- FOR_EACH(i, 0, sa_size(args.props)) {
- printf("%g %g %g %g\n",
- args.props[i].W, args.props[i].Nr, args.props[i].Kr, args.props[i].Ne);
- }
-
-exit:
- sa_release(args.props);
- sa_release(args.wavelengths);
- sa_release(mtl_props);
- return res;
-error:
- goto exit;
-}
-
-
-static res_T
cmd_func_help_schiff(struct sbox* sbox, const int argc, char** argv, void* ctx)
{
struct sbox_schiff* schiff = ctx;
@@ -329,43 +84,6 @@ cmd_func_help_schiff(struct sbox* sbox, const int argc, char** argv, void* ctx)
}
}
-static res_T
-cmd_func_help_schiff_sphere
- (struct sbox* sbox, const int argc, char** argv, void* ctx)
-{
- static const char* help[] = {
-"The geometry of the micro organisms are distributed with respect to a lognormal\n",
-"distribution of spherical meshes:\n",
-"P(x) dx = 1/(sigma*x*sqrt(2*PI)) * exp(-(ln(x)-radius)^2 / (2*sigma^2)) dx\n",
-"\n",
-"schiff sphere [OPTIONS]...\n",
-"\n",
-"-d DIRECTIONS Number of directions sampled for each spherical geometry.\n",
-" Default is "STR(NDIRS)".\n",
-"-g GEOMETRIES Number of spherical geometries sampled with the lognormal\n",
-" distribution. This is actually the number of realisations of\n",
-" the integration. Default is "STR(NGEOMS)".\n",
-"-i PROPERTIES File that lists the radiative properties of the geometries for\n",
-" a set of wavelengths. Each line must be formatted as \"W Nr Kr Ne\"\n",
-" where \"W\" is the wavelength in meter, \"Nr\" and \"Kr\" are\n",
-" respectively the real and imaginary parts of the relative\n",
-" refractive index, and \"Ne\" the refractive index of the medium.\n",
-"-n STEPS Number of points used to discretised the spherical mesh along 2PI.\n",
-" Default is "STR(NTHETAS)".\n",
-"-o FILENAME Write output to FILENAME.\n",
-"-p W:Nr:Kr:Ne Optical properties of the geometries for a given wavelength. \"W\"\n",
-" is the wavelength in meter, \"Nr\" and \"Kr\" are respectively the\n",
-" real and imaginary parts of the relative refractive index and\n",
-" \"Ne\" is the refractive index of the medium. Overwrite the\n",
-" properties loaded with the -i option.\n",
-"-r RADIUS \"radius\" argument of the lognormal distribution Default is "STR(RADIUS)".\n",
-"-s SIGMA \"sigma\" argument of the lognormal distribution. Default is "STR(SIGMA)".\n",
-"-w WAVELENGTH Wavelength to integrate.\n"
- };
- (void)argc, (void)argv, (void)ctx;
- return sbox_print_text(sbox, help, sizeof(help)/sizeof(const char*));
-}
-
/*******************************************************************************
* Exported functions
******************************************************************************/
@@ -376,6 +94,7 @@ res_T
sbox_plugin_create(struct sbox* sbox, void** out_plugin)
{
struct sbox_schiff* schiff = NULL;
+ struct mem_allocator allocator;
res_T res = RES_OK;
if(!sbox || !out_plugin) {
@@ -383,13 +102,16 @@ sbox_plugin_create(struct sbox* sbox, void** out_plugin)
goto error;
}
- schiff = MEM_CALLOC(sbox->allocator, 1, sizeof(struct sbox_schiff));
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ schiff = MEM_CALLOC(&allocator, 1, sizeof(struct sbox_schiff));
if(!schiff) {
logger_print(&sbox->logger, LOG_ERROR,
"Not enough memory: couldn't create the Schiff Star-Box plugin.\n");
res = RES_MEM_ERR;
goto error;
}
+ schiff->allocator = allocator;
SBOX(cmd_init(&schiff->cmd_schiff, "schiff",
"Estimation of the radiative properties of micro organisms.",
@@ -430,12 +152,24 @@ res_T
sbox_plugin_release(struct sbox* sbox, void* addin)
{
struct sbox_schiff* schiff = addin;
+ struct mem_allocator allocator = schiff->allocator;
+ res_T res = RES_OK;
+
if(!sbox || !schiff) return RES_BAD_ARG;
+
SBOX(cmd_unregister(&schiff->cmd_schiff));
SBOX(cmd_unregister(&schiff->cmd_schiff_sphere));
SBOX(cmd_unregister(&schiff->cmd_help_schiff));
SBOX(cmd_unregister(&schiff->cmd_help_schiff_sphere));
- MEM_RM(sbox->allocator, schiff);
- return RES_OK;
+ MEM_RM(&allocator, schiff);
+
+ if(MEM_ALLOCATED_SIZE(&allocator)) {
+ char dump[512];
+ MEM_DUMP(&allocator, dump, sizeof(dump)/sizeof(char));
+ logger_print(&sbox->logger, LOG_ERROR, "Memory leaks:\n%s\n", dump);
+ res = RES_MEM_ERR;
+ }
+ mem_shutdown_proxy_allocator(&allocator);
+ return res;
}
diff --git a/src/sbox_schiff.h b/src/sbox_schiff.h
@@ -36,6 +36,8 @@ struct sbox_schiff {
struct sbox_cmd cmd_schiff_sphere;
struct sbox_cmd cmd_help_schiff;
struct sbox_cmd cmd_help_schiff_sphere;
+
+ struct mem_allocator allocator;
};
#endif /* SBOX_SCHIFF_H */
diff --git a/src/sbox_schiff_sphere.c b/src/sbox_schiff_sphere.c
@@ -0,0 +1,603 @@
+/* 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. */
+
+#define _POSIX_C_SOURCE 2 /* getopt support */
+
+#include "sbox_schiff.h"
+#include "sbox_schiff_optical_properties.h"
+#include "sbox_schiff_sphere.h"
+
+#include <star/s3d.h>
+#include <star/sbox_parse.h>
+#include <star/ssp.h>
+#include <star/sschiff.h>
+
+#include <rsys/clock_time.h>
+#include <rsys/float3.h>
+#include <rsys/stretchy_array.h>
+
+#ifdef COMPILER_CL
+ #include <getopt.h>
+#else
+ #include <unistd.h>
+#endif
+
+
+/* Default values */
+#define NDIRS 100
+#define NGEOMS 1000
+#define RADIUS 1.0
+#define SIGMA 1.0
+#define NTHETAS 64
+
+struct schiff_sphere_args {
+ const char* material;
+ const char* output;
+ struct sbox_schiff_optical_properties* props;
+ double* wavelengths;
+ double radius;
+ double sigma;
+ unsigned ndirs;
+ unsigned ngeoms;
+ unsigned nthetas;
+};
+
+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 };
+
+struct sphere {
+ struct mesh* mesh;
+ float radius;
+};
+
+struct schiff_sphere_geometry_distribution_context {
+ struct mesh* mesh;
+ struct sbox_schiff_optical_properties* properties;
+ double mean_radius;
+ double sigma;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static char
+check_wavelength(const double wavelength)
+{
+ ASSERT(wavelength);
+ return wavelength > 0.0;
+}
+
+static int
+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;
+}
+
+static void
+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];
+}
+
+static INLINE 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;
+}
+
+static void
+get_material_property
+ (void* mtl,
+ const double wavelength,
+ struct sschiff_material_properties* props)
+{
+ struct sbox_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 */
+ 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;
+ 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;
+ }
+
+ NCHECK(props, NULL);
+ props->medium_refractive_index = Ne;
+ props->relative_imaginary_refractive_index = Kr;
+ 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
+merge_properties
+ (struct sbox_schiff_optical_properties** props0,
+ const struct sbox_schiff_optical_properties* props1)
+{
+ const size_t sizeof_prop = sizeof(struct sbox_schiff_optical_properties);
+ const struct sbox_schiff_optical_properties* src0;
+ const struct sbox_schiff_optical_properties* src1;
+ struct sbox_schiff_optical_properties* dst = NULL;
+ size_t nsrc0, nsrc1;
+ size_t isrc0 = 0, isrc1 =0;
+ ASSERT(props0 && props1);
+
+ src0 = *props0;
+ src1 = props1;
+ nsrc0 = sa_size(src0);
+ nsrc1 = sa_size(src1);
+
+ while(isrc0 < nsrc0 && isrc1 < nsrc1) {
+ if(src0[isrc0].W < src1[isrc1].W) {
+ sa_push(dst, src0[isrc0]);
+ ++isrc0;
+ } else if(src1[isrc1].W < src0[isrc0].W) {
+ sa_push(dst, src1[isrc1]);
+ ++isrc1;
+ } else { /* src0[isrc0].W == src1[isrc1].W */
+ sa_push(dst, src0[isrc0]);
+ ++isrc0;
+ ++isrc1;
+ }
+
+ if(isrc0 >= nsrc0 && isrc1 < nsrc1) {
+ const size_t nremains = nsrc1 - isrc1;
+ memcpy(sa_add(dst, nremains), src1 + isrc1, sizeof_prop * nremains);
+ isrc1 = nsrc1;
+ } else if(isrc1 >= nsrc1 && isrc0 < nsrc0) {
+ const size_t nremains = nsrc0 - isrc0;
+ memcpy(sa_add(dst, nremains), src0 + isrc0, sizeof_prop * nremains);
+ isrc0 = nsrc0;
+ }
+ }
+ sa_release(*props0);
+ *props0 = dst;
+}
+
+
+static res_T
+schiff_sphere_parse_args
+ (struct sbox* sbox,
+ const int argc,
+ char** argv,
+ struct schiff_sphere_args* args)
+{
+ double props[4];
+ double* wavelength;
+ int opt;
+ res_T res = RES_OK;
+ ASSERT(sbox && argc && argv && args);
+
+ *args = SCHIFF_SPHERE_ARGS_NULL;
+ optind = 0;
+ opterr = 0;
+
+ while((opt = getopt(argc, argv, ":d:g:i:n:p:r:s:w:")) != -1) {
+ res_T res = RES_OK;
+ switch(opt) {
+ case 'd': res = sbox_str_to_uint(optarg, &args->ndirs); break;
+ case 'g': res = sbox_str_to_uint(optarg, &args->ngeoms); break;
+ case 'i': args->material = optarg; break;
+ case 'n': res = sbox_str_to_uint(optarg, &args->nthetas); break;
+ case '0': args->output = optarg; break;
+ case 'p':
+ res = sbox_parse_list_double(optarg, props, 4);
+ if(res == RES_OK) {
+ struct sbox_schiff_optical_properties* opt_props = sa_add(args->props, 1);
+ opt_props->W = props[0];
+ opt_props->Nr = props[1];
+ opt_props->Kr = props[2];
+ opt_props->Ne = props[3];
+ if(!sbox_schiff_optical_properties_check(opt_props))
+ res = RES_BAD_ARG;
+ }
+ break;
+ case 'r':
+ res = sbox_str_to_double(optarg, &args->radius);
+ if(res == RES_OK && args->radius < 0.0) res = RES_BAD_ARG;
+ break;
+ case 's':
+ res = sbox_str_to_double(optarg, &args->sigma);
+ if(res == RES_OK && args->sigma < 0.0) res = RES_BAD_ARG;
+ break;
+ case 'w':
+ wavelength = sa_add(args->wavelengths, 1);
+ res = sbox_str_to_double(optarg, wavelength);
+ if(res == RES_OK && !check_wavelength(*wavelength)) res = RES_BAD_ARG;
+ break;
+ case ':':
+ logger_print(&sbox->logger, LOG_ERROR,
+ "%s: option requires an argument -- '%c'\n",
+ argv[0], argv[optind-1][1]);
+ res = RES_BAD_ARG;
+ break;
+ case '?':
+ logger_print(&sbox->logger, LOG_ERROR,
+ "%s: invalid option -- '%c'\n",
+ argv[0], argv[optind-1][1]);
+ res = RES_BAD_ARG;
+ break;
+ default: res = RES_BAD_ARG; break;
+ }
+ if(res != RES_OK) {
+ if(optarg) {
+ logger_print(&sbox->logger, LOG_ERROR,
+ "%s: invalid option arguments '%s' -- '%c'\n",
+ argv[0], optarg, opt);
+ }
+ goto error;
+ }
+ }
+exit:
+ return res;
+error:
+ sa_release(args->props);
+ sa_release(args->wavelengths);
+ *args = SCHIFF_SPHERE_ARGS_NULL;
+ goto exit;
+}
+
+static res_T
+schiff_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 sphere;
+ struct s3d_vertex_data attrib;
+ size_t nverts, nprims;
+ ASSERT(rng && mtl && shape && context);
+
+ sphere.mesh = ctx->mesh;
+ sphere.radius = (float)ssp_ran_lognormal
+ (rng, log(ctx->mean_radius), log(ctx->sigma));
+
+ attrib.usage = S3D_POSITION;
+ attrib.type = S3D_FLOAT3;
+ attrib.get = get_position;
+
+ 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*/;
+
+ return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims,
+ get_indices, (unsigned)nverts, &attrib, 1, &sphere);
+}
+
+static res_T
+schiff_sphere_run
+ (struct sbox* sbox,
+ struct sbox_schiff* schiff,
+ struct schiff_sphere_args* args)
+{
+ struct sschiff_geometry_distribution distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION;
+ struct schiff_sphere_geometry_distribution_context ctx;
+ struct mesh mesh = MESH_NULL;
+ struct sschiff_device* sschiff = NULL;
+ struct ssp_rng* rng = NULL;
+ size_t iwave, nwaves;
+ res_T res = RES_OK;
+ ASSERT(sbox);
+
+ mesh_init_sphere(&mesh, args->nthetas);
+
+ /* FIXME use the rng type defined in sbox */
+ res = ssp_rng_create(&schiff->allocator, &ssp_rng_threefry, &rng);
+ if(res != RES_OK) {
+ logger_print(&sbox->logger, LOG_ERROR,
+ "Couldn't create the Star Schiff Random Number Generator.\n");
+ goto error;
+ }
+
+ res = sschiff_device_create
+ (&sbox->logger, &schiff->allocator, 0, sbox->s3d, &sschiff);
+ if(res != RES_OK) {
+ logger_print(&sbox->logger, LOG_ERROR,
+ "Couldn't create the Star Schiff device.\n");
+ goto error;
+ }
+
+ ctx.mesh = &mesh;
+ ctx.properties = args->props;
+ ctx.mean_radius = args->radius;
+ ctx.sigma = args->sigma;
+
+ 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,
+ 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_sub(&t0, &t1, &t0);
+ time_dump(&t0, TIME_MIN|TIME_SEC|TIME_MSEC, NULL, buf, sizeof(buf));
+
+ SSCHIFF(estimator_get_extinction_cross_section(estimator, &status));
+ 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);
+
+ SSCHIFF(estimator_ref_put(estimator));
+ }
+
+exit:
+ if(sschiff) SSCHIFF(device_ref_put(sschiff));
+ if(rng) SSP(rng_ref_put(rng));
+ mesh_release(&mesh);
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+cmd_func_schiff_sphere(struct sbox* sbox, const int argc, char** argv, void* ctx)
+{
+ struct sbox_schiff* schiff = ctx;
+ struct schiff_sphere_args args = SCHIFF_SPHERE_ARGS_NULL;
+ struct sbox_schiff_optical_properties* mtl_props = NULL;
+ const size_t sizeof_prop = sizeof(struct sbox_schiff_optical_properties);
+ res_T res = RES_OK;
+ (void)ctx;
+
+ res = schiff_sphere_parse_args(sbox, argc, argv, &args);
+ if(res != RES_OK) goto error;
+
+ if(!args.wavelengths) {
+ logger_print(&sbox->logger, LOG_ERROR,
+ "Missing the wavelength argument(s).\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Sort the submitted properties in ascending order with respect to the
+ * wavelength */
+ if(args.props) {
+ qsort(args.props, sa_size(args.props), sizeof_prop, cmp_properties);
+ }
+ /* Load the material file and sort the its properties in ascending order with
+ * respect to the wavelength */
+ if(args.material) {
+ res = sbox_schiff_optical_properties_load(sbox, args.material, &mtl_props);
+ if(res != RES_OK) goto error;
+ qsort(mtl_props, sa_size(mtl_props), sizeof_prop, cmp_properties);
+ }
+
+ /* Define the optical properties to use */
+ if(mtl_props && !args.props) {
+ args.props = mtl_props;
+ mtl_props = NULL;
+ } else if(args.props && mtl_props) {
+ merge_properties(&args.props, mtl_props);
+ sa_release(mtl_props);
+ mtl_props = NULL;
+ }
+
+ if(!args.props) {
+ logger_print(&sbox->logger, LOG_ERROR,
+ "Radiative properties are not defined.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = schiff_sphere_run(sbox, schiff, &args);
+ if(res != RES_OK) goto error;
+
+exit:
+ sa_release(args.props);
+ sa_release(args.wavelengths);
+ sa_release(mtl_props);
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+cmd_func_help_schiff_sphere
+ (struct sbox* sbox, const int argc, char** argv, void* ctx)
+{
+ static const char* help[] = {
+"The geometry of the micro organisms are distributed with respect to a lognormal\n",
+"distribution of spherical meshes:\n",
+"P(x) dx = 1/(sigma*x*sqrt(2*PI)) * exp(-(ln(x)-ln(radius))^2 / (2*ln(sigma)^2)) dx\n",
+"\n",
+"schiff sphere [OPTIONS]...\n",
+"\n",
+"-d DIRECTIONS Number of directions sampled for each spherical geometry.\n",
+" Default is "STR(NDIRS)".\n",
+"-g GEOMETRIES Number of spherical geometries sampled with the lognormal\n",
+" distribution. This is actually the number of realisations of\n",
+" the integration. Default is "STR(NGEOMS)".\n",
+"-i PROPERTIES File that lists the radiative properties of the geometries for\n",
+" a set of wavelengths. Each line must be formatted as \"W Nr Kr Ne\"\n",
+" where \"W\" is the wavelength in meter, \"Nr\" and \"Kr\" are\n",
+" respectively the real and imaginary parts of the relative\n",
+" refractive index, and \"Ne\" the refractive index of the medium.\n",
+"-n STEPS Number of points used to discretised the spherical mesh along 2PI.\n",
+" Default is "STR(NTHETAS)".\n",
+"-o FILENAME Write output to FILENAME.\n",
+"-p W:Nr:Kr:Ne Optical properties of the geometries for a given wavelength. \"W\"\n",
+" is the wavelength in meter, \"Nr\" and \"Kr\" are respectively the\n",
+" real and imaginary parts of the relative refractive index and\n",
+" \"Ne\" is the refractive index of the medium. Overwrite the\n",
+" properties loaded with the -i option.\n",
+"-r RADIUS \"radius\" argument of the lognormal distribution Default is "STR(RADIUS)".\n",
+"-s SIGMA \"sigma\" argument of the lognormal distribution. Default is "STR(SIGMA)".\n",
+"-w WAVELENGTH Wavelength to integrate.\n"
+ };
+ (void)argc, (void)argv, (void)ctx;
+ return sbox_print_text(sbox, help, sizeof(help)/sizeof(const char*));
+}
+
diff --git a/src/sbox_schiff_sphere.h b/src/sbox_schiff_sphere.h
@@ -0,0 +1,50 @@
+/* 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_SPHERE_H
+#define SBOX_SCHOFF_SPHERE_H
+
+#include <rsys/rsys.h>
+
+struct sbox;
+
+extern LOCAL_SYM res_T
+cmd_func_schiff_sphere
+ (struct sbox* sbox,
+ const int argc,
+ char** argv,
+ void* ctx);
+
+extern LOCAL_SYM res_T
+cmd_func_help_schiff_sphere
+ (struct sbox* sbox,
+ const int argc,
+ char** argv,
+ void* ctx);
+
+#endif /* SBOX_SCHIFF_SPHERE_H */
