commit 690b15153d71fbdb34ceb436df39563ab749506a
parent bd69e8b92cb5ee44ee416ad205b2f030e2edb632
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 29 Oct 2015 08:29:05 +0100
Merge branch 'binary'
Diffstat:
20 files changed, 1265 insertions(+), 1338 deletions(-)
diff --git a/README.md b/README.md
@@ -1,20 +1,17 @@
-# Star Schiff
+# Schiff
-This library is a [Star-Box](https://gitlab.com/meso-star/star-box/) plugin
-that estimates the radiative properties of micro organism with an
-"Approximation Method for Short Wavelength or High-Energy Scattering"
-(L. Schiff, 1956).
+The purpose of this program is to estimate the radiative properties of micro
+organism with an "Approximation Method for Short Wavelength or High-Energy
+Scattering" (L. Schiff, 1956).
## How to build
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/)
+[RCMake](https://gitlab.com/vaplv/rcmake/) package to build. It also depends
+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-MC](https://gitlab.com/meso-star/star-mc/),
-[Star-SP](https://gitlab.com/meso-star/star-sp/) libraries, and the
-[Star-Box](https://gitlab.com/meso-schiff/star-box/) software development kit.
+[Star-Schiff](https://gitlab.com/meso-star/star-schiff/) and
+[Star-SP](https://gitlab.com/meso-star/star-sp/) libraries.
First ensure that CMake is installed on your system. Then install the RCMake
package as well as all the aforementioned prerequisites. Then generate the
@@ -24,8 +21,8 @@ the RCMake package.
## Licenses
-*Star Box Schiff* is Copyright (C) |Meso|Star> 2015 (<contact@meso-star.com>).
-It is a free software released under the [OSI](http://opensource.org)-approved
-CeCILL license. You are welcome to redistribute it under certain conditions;
-refer to the COPYING files for details.
+*Schiff* is Copyright (C) |Meso|Star> 2015 (<contact@meso-star.com>). It is a
+free software released under the [OSI](http://opensource.org)-approved CeCILL
+license. You are welcome to redistribute it under certain conditions; refer to
+the COPYING files for details.
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -27,9 +27,9 @@
# knowledge of the CeCILL license and that you accept its terms.
cmake_minimum_required(VERSION 2.8)
-project(sbox-4v_s)
+project(schiff C)
-set(SBOX_SCHIFF_SOURCE_DIR ${PROJECT_SOURCE_DIR}/../src/)
+set(SCHIFF_SOURCE_DIR ${PROJECT_SOURCE_DIR}/../src/)
################################################################################
# Check dependencies
@@ -37,7 +37,6 @@ set(SBOX_SCHIFF_SOURCE_DIR ${PROJECT_SOURCE_DIR}/../src/)
find_package(RCMake 0.2 REQUIRED)
find_package(RSys 0.2.1 REQUIRED)
find_package(Star3D 0.3 REQUIRED)
-find_package(StarBox 0.3 REQUIRED)
find_package(StarSchiff 0.3 REQUIRED)
find_package(StarMC 0.3 REQUIRED)
find_package(StarSP 0.3 REQUIRED)
@@ -48,8 +47,9 @@ endif()
include_directories(
${RSys_INCLUDE_DIR}
- ${StarBox_INCLUDE_DIR}
- ${StarSchiff_INCLUDE_DIR})
+ ${Star3D_INCLUDE_DIR}
+ ${StarSchiff_INCLUDE_DIR}
+ ${StarSP_INCLUDE_DIR})
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
include(rcmake)
@@ -63,40 +63,41 @@ set(VERSION_MINOR 3)
set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
-set(SBOX_SCHIFF_FILES_SRC
- sbox_schiff.c
- sbox_schiff_mesh.c
- sbox_schiff_optical_properties.c
- sbox_schiff_sphere.c)
-set(SBOX_SCHIFF_FILES_INC
- sbox_schiff.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)
+set(SCHIFF_FILES_SRC
+ schiff.c
+ schiff_args.c
+ schiff_mesh.c
+ schiff_optical_properties.c
+ schiff_sphere.c)
+set(SCHIFF_FILES_INC
+ schiff_args.h
+ schiff_mesh.h
+ schiff_optical_properties.h
+ schiff_sphere.h)
+set(SCHIFF_FILES_DOC COPYING.fr COPYING.en README.md)
-# Prepend each file in the `SBOX_SCHIFF_FILES_<SRC|INC>' list by the absolute
+# Prepend each file in the `SCHIFF_FILES_<SRC|INC>' list by the absolute
# path of the directory in which they lie
-rcmake_prepend_path(SBOX_SCHIFF_FILES_SRC ${SBOX_SCHIFF_SOURCE_DIR})
-rcmake_prepend_path(SBOX_SCHIFF_FILES_INC ${SBOX_SCHIFF_SOURCE_DIR})
-rcmake_prepend_path(SBOX_SCHIFF_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
+rcmake_prepend_path(SCHIFF_FILES_SRC ${SCHIFF_SOURCE_DIR})
+rcmake_prepend_path(SCHIFF_FILES_INC ${SCHIFF_SOURCE_DIR})
+rcmake_prepend_path(SCHIFF_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
if(CMAKE_COMPILER_IS_GNUCC)
set(MATH_LIB m)
endif()
-add_library(sbox-schiff SHARED ${SBOX_SCHIFF_FILES_SRC} ${SBOX_SCHIFF_FILES_INC})
-target_link_libraries(sbox-schiff RSys Star3D StarBox StarSchiff StarMC StarSP ${MATH_LIB})
-set_target_properties(sbox-schiff PROPERTIES
+add_executable(schiff ${SCHIFF_FILES_SRC} ${SCHIFF_FILES_INC})
+target_link_libraries(schiff RSys Star3D StarSchiff StarSP ${MATH_LIB})
+set_target_properties(schiff PROPERTIES
VERSION ${VERSION}
SOVERSION ${VERSION_MAJOR})
################################################################################
# Define output & install directories
################################################################################
-install(TARGETS sbox-schiff
+install(TARGETS schiff
ARCHIVE DESTINATION bin
LIBRARY DESTINATION lib
RUNTIME DESTINATION bin)
-install(FILES ${SBOX_SCHIFF_FILES_DOC} DESTINATION doc/star-box-schiff)
+install(FILES ${SCHIFF_FILES_DOC} DESTINATION doc/schiff)
diff --git a/src/sbox_schiff.c b/src/sbox_schiff.c
@@ -1,175 +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 "sbox_schiff.h"
-#include "sbox_schiff_sphere.h"
-
-#include <star/sbox.h>
-
-#include <rsys/mem_allocator.h>
-#include <rsys/rsys.h>
-
-/*******************************************************************************
- * Schiff command
- ******************************************************************************/
-static res_T
-cmd_func_schiff(struct sbox* sbox, const int argc, char** argv, void* ctx)
-{
- if(!argc) {
- logger_print(&sbox->logger, LOG_ERROR, "Missing argument.\n");
- return RES_BAD_ARG;
- } else {
- struct sbox_schiff* schiff = ctx;
- ASSERT(schiff);
- return sbox_cmd_exec(sbox, &schiff->cmd_schiff.sub_cmds, argc-1, argv+1);
- }
-}
-
-static res_T
-cmd_func_help_schiff(struct sbox* sbox, const int argc, char** argv, void* ctx)
-{
- struct sbox_schiff* schiff = ctx;
- (void)argc, (void)argv;
-
- if(argc > 1) {
- return sbox_cmd_exec(sbox, &schiff->cmd_help_schiff.sub_cmds, argc-1, argv+1);
- } else {
- static const char* help[] = {
-"Integrate the radiative properties of micro organism with respect to an\n",
-"\"Approximation Method for Short Wavelength or High-Energy Scattering\"\n",
-"(L. Schiff, 1956).\n",
-"\n",
-"schiff GEOMETRY_DISTRIBUTION [OPTION]...\n",
-"\n",
-"The geometries of the micro organisms are defined by the GEOMETRY_DISTRIBUTION\n",
-"that controls the distribution of the geometry shape.\n",
-"Type \"help schiff\" followed by a command for further documentation.\n",
-"\n",
-"List of commands:\n",
-"\n",
- };
- const size_t nlines = sizeof(help)/sizeof(const char*);
- res_T res;
-
- res = sbox_print_text(sbox, help, nlines);
- if(res != RES_OK) return res;
- sbox_cmd_print(sbox, &schiff->cmd_help_schiff.sub_cmds, nlines);
- return res;
- }
-}
-
-/*******************************************************************************
- * Exported functions
- ******************************************************************************/
-EXPORT_SYM res_T sbox_plugin_create(struct sbox* sbox, void** out_plugin);
-EXPORT_SYM res_T sbox_plugin_release(struct sbox* sbox, void* addin);
-
-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) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- 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.",
- cmd_func_schiff, schiff));
- SBOX(cmd_init(&schiff->cmd_help_schiff, "schiff",
- "Estimation of the radiative properties of micro organisms.",
- cmd_func_help_schiff, schiff));
- SBOX(cmd_init(&schiff->cmd_schiff_sphere, "sphere",
- "Estimation of the radiative properties of spherical micro organisms.",
- cmd_func_schiff_sphere, schiff));
- SBOX(cmd_init(&schiff->cmd_help_schiff_sphere, "sphere",
- "Estimation of the radiative properties of spherical micro organisms.",
- cmd_func_help_schiff_sphere, schiff));
-
- res = sbox_cmd_register(&sbox->cmds, &schiff->cmd_schiff);
- if(res != RES_OK) goto error;
- res = sbox_cmd_register(&sbox_cmd_help.sub_cmds, &schiff->cmd_help_schiff);
- if(res != RES_OK) goto error;
- res = sbox_cmd_register
- (&schiff->cmd_schiff.sub_cmds, &schiff->cmd_schiff_sphere);
- if(res != RES_OK) goto error;
- res = sbox_cmd_register
- (&schiff->cmd_help_schiff.sub_cmds, &schiff->cmd_help_schiff_sphere);
- if(res != RES_OK) goto error;
-
-exit:
- if(out_plugin) *out_plugin = schiff;
- return res;
-error:
- if(schiff) {
- SBOX(plugin_release(sbox, schiff));
- schiff = NULL;
- }
- goto exit;
-}
-
-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(&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
@@ -1,44 +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 SBOX_SCHIFF_H
-#define SBOX_SCHIFF_H
-
-#include <star/sbox_cmd.h>
-
-struct sbox_schiff {
- struct sbox_cmd cmd_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_mesh.c b/src/sbox_schiff_mesh.c
@@ -1,192 +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 "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, 0, nphis-1) {
- const double phi = -PI/2 + (double)(iphi + 1) * 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, "w");
- 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) {
- /* Obj indices start from 1 */
- fprintf(fp, "f %u %u %u\n",
- darray_uint_data_get(&mesh->indices)[i*3 + 0] + 1,
- darray_uint_data_get(&mesh->indices)[i*3 + 1] + 1,
- darray_uint_data_get(&mesh->indices)[i*3 + 2] + 1);
- }
-
-exit:
- if(fp) fclose(fp);
- return res;
-error:
- goto exit;
-}
-
diff --git a/src/sbox_schiff_mesh.h b/src/sbox_schiff_mesh.h
@@ -1,56 +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 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_optical_properties.c b/src/sbox_schiff_optical_properties.c
@@ -1,180 +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. */
-
-#define _POSIX_C_SOURCE 200112L /* wordexp support */
-
-#include "sbox_schiff_optical_properties.h"
-
-#include <star/sbox.h>
-#include <star/sbox_parse.h>
-
-#include <rsys/stretchy_array.h>
-
-#ifndef COMPILER_CL
- #include <wordexp.h>
-#endif
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static res_T
-parse_optical_properties
- (struct sbox* sbox,
- const char* filename,
- const size_t iline,
- const char* str,
- struct sbox_schiff_optical_properties* out_props)
-{
- char buf[128];
- char* tk;
- int i = 0;
- double props[4];
- res_T res = RES_OK;
- ASSERT(sbox && filename && str && out_props);
-
- if(strlen(str) >= sizeof(buf) - 1/*NULL char*/) {
- logger_print(&sbox->logger, LOG_ERROR,
- "%s:%lu: not enough memory: cannot parse the properties `%s'.\n",
- filename, (unsigned long)iline, str);
- return RES_MEM_ERR;
- }
-
- strncpy(buf, str, sizeof(buf));
- for(i=0, tk=strtok(buf, " \t"); tk && i < 4; ++i, tk=strtok(NULL, " \t")) {
- res = sbox_str_to_double(tk, props + i);
- if(res != RES_OK) {
- logger_print(&sbox->logger, LOG_ERROR,
- "%s:%lu: cannot read the optical property `%s'.\n",
- filename, (unsigned long)iline, tk);
- return res;
- }
- }
- if(i < 4) {
- logger_print(&sbox->logger, LOG_ERROR,
-"%s:%lu: not enough optical properties.\n"
-"Expect 4 parameters while %d is submitted - `%s'.\n",
- filename, (unsigned long)iline, i, str);
- return RES_BAD_ARG;
- }
- out_props->W = props[0];
- out_props->Nr = props[1];
- out_props->Kr = props[2];
- out_props->Ne = props[3];
- if(!sbox_schiff_optical_properties_check(out_props)) {
- logger_print(&sbox->logger, LOG_ERROR,
- "%s:%lu: invalid optical properties `%s'.\n",
- filename, (unsigned long)iline, str);
- return RES_BAD_ARG;
- }
- return RES_OK;
-}
-
-/*******************************************************************************
- * Local function
- ******************************************************************************/
-res_T
-sbox_schiff_optical_properties_load
- (struct sbox* sbox,
- const char* filename,
- struct sbox_schiff_optical_properties** out_props)
-{
- struct sbox_schiff_optical_properties* props = NULL;
- char* buf = NULL;
- const size_t buf_chunk = 128;
- size_t iline;
- FILE* fp;
- res_T res = RES_OK;
- ASSERT(sbox && filename && out_props);
-
-#ifdef COMPILER_CL
- fp = fopen(filename, "r");
-#else
- {
- wordexp_t wexp;
- int err = wordexp(filename, &wexp, 0);
- if(err) {
- logger_print(&sbox->logger, LOG_ERROR,
- "Error parsing the filenames `%s'\n", filename);
- res = RES_BAD_ARG;
- goto error;
- }
- ASSERT(wexp.we_wordc == 1);
- fp = fopen(wexp.we_wordv[0], "r");
- wordfree(&wexp);
- }
-#endif
- if(!fp) {
- logger_print(&sbox->logger, LOG_ERROR,
- "Cannot open the file of Schiff optical properties `%s'.\n", filename);
- res = RES_IO_ERR;
- goto error;
- }
-
- buf = sa_add(buf, buf_chunk);
- for(iline = 1; fgets(buf, (int)sa_size(buf), fp); ++iline) {
- char* line;
- size_t last_char;
- struct sbox_schiff_optical_properties tmp_props;
-
- while(!strrchr(buf, '\n')) { /* Ensure that the whole line is read */
- if(!fgets(sa_add(buf, buf_chunk), (int)buf_chunk, fp)) /* EOF */
- break;
- }
-
- /* Remove leading spaces */
- line = buf;
- while((*line == ' ' || *line == '\t') && *line != '\0') ++line;
-
- /* Remove newline character(s) */
- last_char = strlen(line);
- while(last_char-- && (line[last_char]=='\n' || line[last_char]=='\r'));
- line[last_char + 1] = '\0';
-
- if(*line == '\0' /*Empty line*/ || *line == '#' /* Comment */)
- continue;
-
- /* Read optical properties */
- res = parse_optical_properties
- (sbox, filename, iline, strtok(line, "#"), &tmp_props);
- if(res == RES_OK) { /* *NO* error */
- sa_push(props, tmp_props);
- }
- }
-
-exit:
- sa_release(buf);
- *out_props = props;
- return res;
-error:
- if(props) {
- sa_release(props);
- props = NULL;
- }
- goto exit;
-}
-
diff --git a/src/sbox_schiff_optical_properties.h b/src/sbox_schiff_optical_properties.h
@@ -1,61 +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 SBOX_SCHIFF_OPTICAL_PROPERTIES_H
-#define SBOX_SCHIFF_OPTICAL_PROPERTIES_H
-
-#include <rsys/rsys.h>
-
-struct sbox;
-
-struct sbox_schiff_optical_properties {
- double W; /* Wavelength */
- double Nr; /* Real part of the relative refractive index */
- double Kr; /* Imaginary part of the relative refractive index */
- double Ne; /* Refractive index of the medium */
-};
-
-static INLINE char
-sbox_schiff_optical_properties_check
- (const struct sbox_schiff_optical_properties* props)
-{
- ASSERT(props);
- return props->W > 0.0
- && props->Nr > 0.0
- && props->Kr > 0.0
- && props->Ne > 0.0;
-}
-
-extern LOCAL_SYM res_T
-sbox_schiff_optical_properties_load
- (struct sbox* sbox,
- const char* filename,
- struct sbox_schiff_optical_properties** props); /* Stretchy array */
-
-#endif /* SBOX_SCHIFF_OPTICAL_PROPERTIES_H */
-
diff --git a/src/sbox_schiff_sphere.c b/src/sbox_schiff_sphere.c
@@ -1,540 +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. */
-
-#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"
-
-#include <star/s3d.h>
-#include <star/sbox_parse.h>
-#include <star/ssp.h>
-#include <star/smc.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
-
-/* Argument of the schiff sphere command */
-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
-};
-
-/* Sphere geometry */
-struct sphere {
- 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 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 */
-};
-
-/*******************************************************************************
- * 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;
- return prop0->W < prop1->W ? -1 : (prop0->W > prop1->W ? 1 : 0);
-}
-
-static void
-get_indices(const unsigned itri, unsigned ids[3], void* ctx)
-{
- struct sphere* sphere = ctx;
- const size_t i = itri * 3;
- 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 void
-get_position(const unsigned ivert, float vertex[3], void* ctx)
-{
- struct sphere* sphere = ctx;
- const size_t i = ivert * 3;
- 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
-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 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;
-}
-
-/* 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 = 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);
-}
-
-static res_T
-schiff_sphere_run
- (struct sbox* sbox,
- 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 sschiff_estimator* estimator = NULL;
- struct sschiff_estimator_state* states = NULL;
- struct sschiff_device* sschiff = NULL;
- struct ssp_rng_type rng_type;
- struct ssp_rng* rng = NULL;
- struct sbox_schiff_mesh mesh;
- struct time t0, t1;
- size_t iwlen, nwlens;
- res_T res = RES_OK;
- ASSERT(sbox);
-
- /* 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);
- goto error;
- }
-
- SMC(device_get_rng_type(sbox->smc, &rng_type));
- res = ssp_rng_create(&schiff->allocator, &rng_type, &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;
- }
-
- /* 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;
-
- /* 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, "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 %g:\n", states[iwlen].wavelength);
-
- val = states[iwlen].values + SSCHIFF_EXTINCTION_CROSS_SECTION;
- logger_print(&sbox->logger, LOG_OUTPUT,
- " Extinction cross section ~ %g +/- %g\n", val->E, val->SE);
-
- val = states[iwlen].values + SSCHIFF_ABSORPTION_CROSS_SECTION;
- logger_print(&sbox->logger, LOG_OUTPUT,
- " Absorption cross section ~ %g +/- %g\n", val->E, val->SE);
-
- val = states[iwlen].values + SSCHIFF_SCATTERING_CROSS_SECTION;
- logger_print(&sbox->logger, LOG_OUTPUT,
- " Scattering cross section ~ %g +/- %g\n", val->E, val->SE);
-
- val = states[iwlen].values + SSCHIFF_AVERAGE_PROJECTED_AREA;
- logger_print(&sbox->logger, LOG_OUTPUT,
- " Average projected area ~ %g +/- %g\n", val->E, val->SE);
- }
- logger_print(&sbox->logger, LOG_OUTPUT, "Elapsed time: %s\n", buf);
-
-exit:
- if(estimator) SSCHIFF(estimator_ref_put(estimator));
- if(sschiff) SSCHIFF(device_ref_put(sschiff));
- if(rng) SSP(rng_ref_put(rng));
- sa_release(states);
- sbox_schiff_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
@@ -1,50 +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 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 */
diff --git a/src/schiff.c b/src/schiff.c
@@ -0,0 +1,51 @@
+/* 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 "schiff_args.h"
+#include "schiff_sphere.h"
+
+int
+main(int argc, char** argv)
+{
+ struct schiff_args args = SCHIFF_ARGS_NULL;
+ int err = 0;
+ res_T res;
+
+ res = schiff_args_init(&args, argc, argv);
+ if(res != RES_OK) goto error;
+ if(!args.wavelengths) goto exit;
+ res = schiff_run_sphere(&args);
+ if(res != RES_OK) goto error;
+
+exit:
+ schiff_args_release(&args);
+ return err;
+error:
+ err = 1;
+ goto exit;
+}
diff --git a/src/schiff.h b/src/schiff.h
@@ -0,0 +1,39 @@
+/* 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
@@ -0,0 +1,200 @@
+/* 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 "schiff_args.h"
+#include "schiff_optical_properties.h"
+
+#include <rsys/cstr.h>
+#include <rsys/stretchy_array.h>
+
+#ifdef COMPILER_CL
+ #include <getopt.h>
+#else
+ #include <unistd.h>
+#endif
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+static void
+print_help(const char* binary)
+{
+ printf(
+"Usage: %s [OPTION]... [FILE]\n"
+"Estimate the radiative properties of micro organisms with an \"Approximation\n"
+"Method for Short Wavelength or High Energy Scattering\" (L. Schiff, 1956).\n\n",
+ binary);
+ printf(
+"FILE lists the per wavelength optical properties of the micro organisms.\n"
+"Each line must be formatted as \"W Nr Kr Ne\" where \"W\" is the wavelength\n"
+"in meter, \"Nr\" and \"Kr\" are the real and imaginary parts, respectively, of\n"
+"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(
+" -d DIRS number of sampled directions for each geometry. Default is %u.\n",
+ SCHIFF_ARGS_NULL.ndirs);
+ printf(
+" -g GOEMS number of sampled geometries. This is actually the number of\n"
+" realisations. Default is %u.\n",
+ SCHIFF_ARGS_NULL.ngeoms);
+ printf(
+" -o OUTPUT write results to OUTPUT. If not defined, write results to\n"
+" stdout.\n");
+ printf(
+" -s R:S[:N] the micro organisms are spherical meshes discretised 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.nthetas);
+ printf(
+" -w A[:B]... list of wavelengths to integrate in meter. At least one\n"
+" wavelength must be provided.\n");
+}
+
+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
+ ******************************************************************************/
+res_T
+schiff_args_init
+ (struct schiff_args* args,
+ const int argc,
+ char** argv)
+{
+ size_t len;
+ double list[3];
+ int opt;
+ res_T res = RES_OK;
+ ASSERT(argc && argv && args);
+
+ *args = SCHIFF_ARGS_NULL;
+ while((opt = getopt(argc, argv, "d:g:ho:s:w:")) != -1) {
+ res_T res = RES_OK;
+ switch(opt) {
+ case 'd': res = cstr_to_uint(optarg, &args->ndirs); break;
+ case 'g': res = cstr_to_uint(optarg, &args->ngeoms); break;
+ case 'h':
+ print_help(argv[0]);
+ schiff_args_release(args);
+ return RES_OK;
+ case 'o': args->output_filename = optarg; break;
+ case 's':
+ res = cstr_to_list_double(optarg, list, &len, 3);
+ if(res == RES_OK) {
+ if(len < 2) {
+ res = RES_BAD_ARG;
+ } else {
+ args->distribution.sphere = SCHIFF_DISTRIBUTION_SPHERE_DEFAULT;
+ args->distribution.sphere.radius = list[0];
+ args->distribution.sphere.sigma = list[1];
+ if(len > 2) {
+ args->distribution.sphere.nthetas = (unsigned)list[2];
+ if((double)args->distribution.sphere.nthetas != list[2])
+ res = RES_BAD_ARG;
+ }
+ }
+ }
+ break;
+ case 'w':
+ sa_clear(args->wavelengths);
+ res = cstr_to_list_double(optarg, NULL, &len, 0);
+ if(res == RES_OK) {
+ args->wavelengths = sa_add(args->wavelengths, len);
+ res = cstr_to_list_double(optarg, args->wavelengths, NULL, len);
+ if(res == RES_OK) {
+ size_t i;
+ FOR_EACH(i, 0, len) {
+ if(args->wavelengths[i] < 0.0) {
+ res = RES_BAD_ARG;
+ break;
+ }
+ }
+ }
+ }
+ break;
+ default: res = RES_BAD_ARG; break;
+ }
+ if(res != RES_OK) {
+ if(optarg) {
+ fprintf(stderr, "%s: invalid option arguments '%s' -- '%c'\n",
+ argv[0], optarg, opt);
+ }
+ goto error;
+ }
+ }
+
+ if(!args->wavelengths) {
+ fprintf(stderr, "Missing wavelength(s).\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(optind == argc) {
+ res = schiff_optical_properties_load_stream(&args->properties, stdin, "stdin");
+ } else {
+ res = schiff_optical_properties_load(&args->properties, argv[optind]);
+ }
+ if(res != RES_OK) goto error;
+
+ if(!args->properties) {
+ fprintf(stderr, "Missing optical properties.\n");
+ 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;
+error:
+ schiff_args_release(args);
+ *args = SCHIFF_ARGS_NULL;
+ goto exit;
+}
+
+void
+schiff_args_release(struct schiff_args* args)
+{
+ ASSERT(args);
+ sa_release(args->properties);
+ sa_release(args->wavelengths);
+ *args = SCHIFF_ARGS_NULL;
+}
+
diff --git a/src/schiff_args.h b/src/schiff_args.h
@@ -0,0 +1,78 @@
+/* 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_ARGS_H
+#define SCHIFF_ARGS_H
+
+#include <rsys/rsys.h>
+
+struct schiff_distribution_sphere {
+ double radius;
+ double sigma;
+ unsigned nthetas;
+};
+#define SCHIFF_DISTRIBUTION_SPHERE_DEFAULT__ {1.0, 1.0, 64}
+static const struct schiff_distribution_sphere
+SCHIFF_DISTRIBUTION_SPHERE_DEFAULT = SCHIFF_DISTRIBUTION_SPHERE_DEFAULT__;
+
+union schiff_distribution {
+ struct schiff_distribution_sphere sphere;
+};
+
+#define SCHIFF_DISTRIBUTION_DEFAULT__ {SCHIFF_DISTRIBUTION_SPHERE_DEFAULT__}
+
+struct schiff_args {
+ const char* output_filename;
+ struct schiff_optical_properties* properties;
+ double* wavelengths;
+ union schiff_distribution distribution;
+ unsigned ngeoms;
+ unsigned ndirs;
+};
+
+static const struct schiff_args SCHIFF_ARGS_NULL = {
+ NULL, /* Output filename */
+ NULL, /* List of optical properties */
+ NULL, /* List of wavelength to integrate */
+ SCHIFF_DISTRIBUTION_DEFAULT__,
+ 1000, /* # Sampled geometries */
+ 100, /* # Sampled directions per gemetry */
+};
+
+extern LOCAL_SYM res_T
+schiff_args_init
+ (struct schiff_args* args,
+ const int argc,
+ char** argv);
+
+extern LOCAL_SYM void
+schiff_args_release
+ (struct schiff_args* args);
+
+#endif /* SCHIFF_ARGS_H */
+
diff --git a/src/schiff_mesh.c b/src/schiff_mesh.c
@@ -0,0 +1,267 @@
+/* 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 "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
+schiff_mesh_init_sphere
+ (struct mem_allocator* allocator,
+ struct 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, 0, nphis-1) {
+ const double phi = -PI/2 + (double)(iphi + 1) * 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:
+ schiff_mesh_release(sphere);
+ goto exit;
+}
+
+res_T
+schiff_mesh_init_cylinder
+ (struct mem_allocator* allocator,
+ struct schiff_mesh* cylinder,
+ const unsigned nsteps)
+{
+ const double step = 2*PI / (double)nsteps;
+ size_t nverts;
+ size_t ntris;
+ unsigned i;
+ res_T res = RES_OK;
+ ASSERT(allocator && cylinder && nsteps);
+
+ darray_float_init(allocator, &cylinder->vertices);
+ darray_uint_init(allocator, &cylinder->indices);
+
+ nverts = nsteps*2/* #contour verts */ + 2/* #polar verts */;
+ ntris = nsteps*2/* #contour tris */ + 2*nsteps/* #caop tris */;
+
+ res = darray_float_resize(&cylinder->vertices, nverts*3/*#coords per vert*/);
+ if(res != RES_OK) goto error;
+ res = darray_uint_resize(&cylinder->indices, ntris*3/*#indices per tri*/);
+ if(res != RES_OK) goto error;
+
+ /* Generate the vertex coordinates */
+ FOR_EACH(i, 0, nsteps) {
+ const float theta = (float)(i* step);
+ const float x = (float)cos(theta);
+ const float y = (float)sin(theta);
+ f3(darray_float_data_get(&cylinder->vertices) + (i*2 + 0)*3, x, y, 0.f);
+ f3(darray_float_data_get(&cylinder->vertices) + (i*2 + 1)*3, x, y, 1.f);
+ }
+
+ /* "Polar" vertices */
+ f3(darray_float_data_get(&cylinder->vertices) + (i*2 + 0)*3, 0.f, 0.f, 0.f);
+ f3(darray_float_data_get(&cylinder->vertices) + (i*2 + 1)*3, 0.f, 0.f, 1.f);
+
+ /* Contour primitives */
+ FOR_EACH(i, 0, nsteps) {
+ const unsigned id = i * 2;
+ unsigned* iprim = darray_uint_data_get(&cylinder->indices) + id*3;
+
+ iprim[0] = (id + 0);
+ iprim[1] = (id + 1);
+ iprim[2] = (id + 2) % (nsteps*2);
+
+ iprim += 3;
+
+ iprim[0] = (id + 2) % (nsteps*2);
+ iprim[1] = (id + 1);
+ iprim[2] = (id + 3) % (nsteps*2);
+ }
+
+ /* Cap primitives */
+ FOR_EACH(i, 0, nsteps) {
+ const unsigned id = i* 2;
+ unsigned* iprim = darray_uint_data_get(&cylinder->indices) + (nsteps + i)*6;
+
+ iprim[0] = (nsteps * 2);
+ iprim[1] = (id + 0);
+ iprim[2] = (id + 2) % (nsteps*2);
+
+ iprim += 3;
+
+ iprim[0] = (nsteps * 2) + 1;
+ iprim[1] = (id + 3) % (nsteps*2);
+ iprim[2] = (id + 1);
+ }
+exit:
+ return res;
+error:
+ schiff_mesh_release(cylinder);
+ goto exit;
+}
+
+void
+schiff_mesh_release(struct schiff_mesh* mesh)
+{
+ ASSERT(mesh);
+ darray_float_release(&mesh->vertices);
+ darray_uint_release(&mesh->indices);
+}
+
+res_T
+schiff_mesh_dump(struct 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, "w");
+ 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) {
+ /* Obj indices start from 1 */
+ fprintf(fp, "f %u %u %u\n",
+ darray_uint_data_get(&mesh->indices)[i*3 + 0] + 1,
+ darray_uint_data_get(&mesh->indices)[i*3 + 1] + 1,
+ darray_uint_data_get(&mesh->indices)[i*3 + 2] + 1);
+ }
+
+exit:
+ if(fp) fclose(fp);
+ return res;
+error:
+ goto exit;
+}
+
diff --git a/src/schiff_mesh.h b/src/schiff_mesh.h
@@ -0,0 +1,62 @@
+/* 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_MESH_H
+#define SCHIFF_MESH_H
+
+#include <rsys/dynamic_array_float.h>
+#include <rsys/dynamic_array_uint.h>
+
+struct schiff_mesh {
+ struct darray_float vertices;
+ struct darray_uint indices;
+};
+
+extern LOCAL_SYM res_T
+schiff_mesh_init_sphere
+ (struct mem_allocator* allocator,
+ struct schiff_mesh* mesh,
+ const unsigned nthetas); /* # discret points along 2PI */
+
+extern LOCAL_SYM res_T
+schiff_mesh_init_cylinder
+ (struct mem_allocator* allocator,
+ struct schiff_mesh* mesh,
+ const unsigned nslices);
+
+extern LOCAL_SYM void
+schiff_mesh_release
+ (struct schiff_mesh* mesh);
+
+extern LOCAL_SYM res_T
+schiff_mesh_dump
+ (struct schiff_mesh* mesh,
+ const char* filename);
+
+#endif /* SBOX_SCHIFF_MESH_H */
+
diff --git a/src/schiff_optical_properties.c b/src/schiff_optical_properties.c
@@ -0,0 +1,166 @@
+/* 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 "schiff_optical_properties.h"
+
+#include <rsys/cstr.h>
+#include <rsys/stretchy_array.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+parse_optical_properties
+ (struct schiff_optical_properties* out_props,
+ const char* filename,
+ const size_t iline,
+ const char* str)
+{
+ char buf[128];
+ char* tk;
+ int i = 0;
+ double props[4];
+ res_T res = RES_OK;
+ ASSERT(filename && str && out_props);
+
+ if(strlen(str) >= sizeof(buf) - 1/*NULL char*/) {
+ fprintf(stderr,
+ "%s:%lu: not enough memory: cannot parse the properties `%s'.\n",
+ filename, (unsigned long)iline, str);
+ return RES_MEM_ERR;
+ }
+
+ strncpy(buf, str, sizeof(buf));
+ for(i=0, tk=strtok(buf, " \t"); tk && i < 4; ++i, tk=strtok(NULL, " \t")) {
+ res = cstr_to_double(tk, props + i);
+ if(res != RES_OK) {
+ fprintf(stderr,
+ "%s:%lu: cannot read the optical property `%s'.\n",
+ filename, (unsigned long)iline, tk);
+ return res;
+ }
+ }
+ if(i < 4) {
+ fprintf(stderr,
+ "%s:%lu: not enough optical properties.\n"
+ "Expect 4 parameters while %d is submitted - `%s'.\n",
+ filename, (unsigned long)iline, i, str);
+ return RES_BAD_ARG;
+ }
+
+ out_props->W = props[0];
+ out_props->Nr = props[1];
+ out_props->Kr = props[2];
+ out_props->Ne = props[3];
+ if(!schiff_optical_properties_check(out_props)) {
+ fprintf(stderr, "%s:%lu: invalid optical properties `%s'.\n",
+ filename, (unsigned long)iline, str);
+ return RES_BAD_ARG;
+ }
+ return RES_OK;
+}
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
+res_T
+schiff_optical_properties_load
+ (struct schiff_optical_properties** props, const char* filename)
+{
+ FILE* fp;
+ res_T res = RES_OK;
+ ASSERT(filename && props);
+
+ fp = fopen(filename, "r");
+ if(!fp) {
+ fprintf(stderr, "Cannot open the file of optical properties `%s'.\n",
+ filename);
+ res = RES_IO_ERR;
+ goto error;
+ }
+
+ res = schiff_optical_properties_load_stream(props, fp, filename);
+ if(res != RES_OK) goto error;
+
+exit:
+ if(fp) fclose(fp);
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+schiff_optical_properties_load_stream
+ (struct schiff_optical_properties** out_props,
+ FILE* stream,
+ const char* stream_name)
+{
+ struct schiff_optical_properties* props = NULL;
+ char* buf = NULL;
+ const size_t buf_chunk = 128;
+ size_t iline;
+ res_T res = RES_OK;
+ ASSERT(out_props && stream && stream_name);
+
+ buf = sa_add(buf, buf_chunk);
+
+ for(iline = 1; fgets(buf, (int)sa_size(buf), stream); ++iline) {
+ char* line;
+ size_t last_char;
+ struct schiff_optical_properties tmp_props;
+
+ while(!strrchr(buf, '\n')) { /* Ensure that the whole line is read */
+ if(!fgets(sa_add(buf, buf_chunk), (int)buf_chunk, stream)) /* EOF */
+ break;
+ }
+
+ /* Remove leading spaces */
+ line = buf;
+ while((*line == ' ' || *line == '\t') && *line != '\0') ++line;
+
+ /* Remove newline character(s) */
+ last_char = strlen(line);
+ while(last_char-- && (line[last_char]=='\n' || line[last_char]=='\r'));
+ line[last_char + 1] = '\0';
+
+ if(*line == '\0' /*Empty line*/ || *line == '#' /* Comment */)
+ continue;
+
+ /* Read optical properties */
+ res = parse_optical_properties
+ (&tmp_props, stream_name, iline, strtok(line, "#"));
+ if(res == RES_OK) { /* *NO* error */
+ sa_push(props, tmp_props);
+ }
+ }
+
+ *out_props = props;
+ sa_release(buf);
+ return RES_OK;
+}
+
diff --git a/src/schiff_optical_properties.h b/src/schiff_optical_properties.h
@@ -0,0 +1,65 @@
+/* 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_OPTICAL_PROPERTIES_H
+#define SCHIFF_OPTICAL_PROPERTIES_H
+
+#include <rsys/dynamic_array.h>
+#include <rsys/rsys.h>
+
+struct schiff_optical_properties {
+ double W; /* Wavelength */
+ double Nr; /* Real part of the relative refractive index */
+ double Kr; /* Imaginary part of the relative refractive index */
+ double Ne; /* Refractive index of the medium */
+};
+
+static INLINE char
+schiff_optical_properties_check
+ (const struct schiff_optical_properties* props)
+{
+ ASSERT(props);
+ return props->W > 0.0
+ && props->Nr > 0.0
+ && props->Kr > 0.0
+ && props->Ne > 0.0;
+}
+
+extern LOCAL_SYM res_T
+schiff_optical_properties_load
+ (struct schiff_optical_properties** props, /* Stretchy array */
+ const char* filename);
+
+extern LOCAL_SYM res_T
+schiff_optical_properties_load_stream
+ (struct schiff_optical_properties** props, /* Stretchy array */
+ FILE* stream,
+ const char* stream_name);
+
+#endif /* SCHIFF_OPTICAL_PROPERTIES_H */
+
diff --git a/src/schiff_sphere.c b/src/schiff_sphere.c
@@ -0,0 +1,261 @@
+/* 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 "schiff_args.h"
+#include "schiff_mesh.h"
+#include "schiff_optical_properties.h"
+#include "schiff_sphere.h"
+
+#include <star/s3d.h>
+#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 {
+ struct 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 schiff_mesh* mesh; /* Triangular mesh of an unit sphere */
+ struct schiff_optical_properties* properties; /* Per wavelength properties */
+ double log_mean_radius; /* Log of the mean sphere radius */
+ double log_sigma; /* Log of the sigma argument of the lognormal distribution*/
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+get_indices(const unsigned itri, unsigned ids[3], void* ctx)
+{
+ struct sphere* sphere = ctx;
+ const size_t i = itri * 3;
+ 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 void
+get_position(const unsigned ivert, float vertex[3], void* ctx)
+{
+ struct sphere* sphere = ctx;
+ const size_t i = ivert * 3;
+ 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
+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;
+}
+
+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, ctx->log_mean_radius, ctx->log_sigma);
+
+ attrib.usage = S3D_POSITION;
+ attrib.type = S3D_FLOAT3;
+ attrib.get = get_position;
+
+ mtl->get_property = 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);
+}
+
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+schiff_run_sphere(const struct schiff_args* args)
+{
+ 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;
+ 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;
+ }
+ }
+
+ /* Generate the spherical mesh */
+ res = schiff_mesh_init_sphere
+ (&mem_default_allocator, &mesh, args->distribution.sphere.nthetas);
+ if(res != RES_OK) {
+ fprintf(stderr,
+ "Couldn't create the sphere mesh discretised in %u steps along 2PI\n",
+ args->distribution.sphere.nthetas);
+ 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");
+ 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;
+ }
+
+ /* Setup the geometry distribution context */
+ ctx.mesh = &mesh;
+ ctx.properties = args->properties;
+ ctx.log_mean_radius = log(args->distribution.sphere.radius);
+ ctx.log_sigma = log(args->distribution.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;
+ }
+
+ /* 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);
+ schiff_mesh_release(&mesh);
+ return res;
+error:
+ goto exit;
+}
+
diff --git a/src/schiff_sphere.h b/src/schiff_sphere.h
@@ -0,0 +1,38 @@
+/* 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_SPHERE_H
+#define SCHIFF_SPHERE_H
+
+#include <rsys/rsys.h>
+
+extern LOCAL_SYM res_T
+schiff_run_sphere
+ (const struct schiff_args* args);
+
+#endif /* SCHIFF_SPHERE_H */
