commit ed1549d80066aabc1b21e3dc5818e66331019651
parent 70bb72288ba6d6d31e8229807c09e6c3592b1b27
Author: Eduardo Fontana Lazzari <edufonlaz@gmail.com>
Date: Mon, 25 Aug 2025 17:07:29 +0200
Introduces spectral properties framework
Up to this point, all physical properties in sphin have been supported
as 'doubles'. However, in photoreactive systems, spectral (i.e.,
per-wavelength) properties are often necessary to fully describe the
system and the underlying physics. This is the case for the radiative
properties of particles and the emission properties of sources, which
may vary significantly across the electromagnetic spectrum. These
properties can also be expressed in different ways depending on the
situation. The absorption properties of a particle, for instance, may be
represented as an absorption coefficient or as absorption cross
sections. They may also be computed using Mie theory.
From an implementation perspective, this translates into a structure
that can easily become complex. This commit introduces the numerical
structure that allows applications using sphin to handle spectral
properties, covering a broad range of the possibilities currently used
in photoreactive systems engineering. A generic structure,
sphin_spectral_property, is added. It contains two lists of equal
length: one with the wavelengths and the other with the corresponding
property values. Physical properties can be constructed upon this
structure. Similar to the representation of geometries, it uses rsys
internally to simplify operations and is exposed with a descriptor
containing simpler types to streamline interactions through the sphin
API.
This commit also begins implementing some of the derived structures that
use sphin_spectral_property, such as absorption cross sections and
emission spectra.
For now, backward compatibility with scalar properties (such as ka) is
preserved.
Diffstat:
16 files changed, 1845 insertions(+), 60 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -7,12 +7,15 @@ file.txt
*.so
*.swp
*.stl
+abs_cross_sec**.txt
+source_spec.txt
tags
star-phor-input.5
sphin.1
sphin
test_sphin
test_sphin_load_geometry
+test_sphin_load_prop_rad
test_sphin_load_source
test_sphin_load_surface
test_sphin_load_volume
diff --git a/Makefile b/Makefile
@@ -41,10 +41,13 @@ SRC =\
src/sphin_brdf.c\
src/sphin_config.c\
src/sphin_geometry.c\
+ src/sphin_prop_rad.c\
+ src/sphin_scatterer.c\
src/sphin_sensor.c\
src/sphin_sensor_surface.c\
src/sphin_sensor_volume.c\
src/sphin_source_surface.c\
+ src/sphin_spectral_property.c\
src/sphin_surface.c\
src/sphin_volume.c
OBJ = $(SRC:.c=.o)
@@ -168,6 +171,7 @@ clean: clean_test clean_tool
TEST_SRC =\
src/test_sphin.c\
src/test_sphin_load_geometry.c\
+ src/test_sphin_load_prop_rad.c\
src/test_sphin_load_source.c\
src/test_sphin_load_surface.c\
src/test_sphin_load_volume.c
@@ -218,6 +222,7 @@ $(TEST_OBJ): config.mk sphin-local.pc
test_sphin \
test_sphin_load_geometry \
+test_sphin_load_prop_rad \
test_sphin_load_source \
test_sphin_load_surface \
test_sphin_load_volume \
@@ -226,5 +231,8 @@ test_sphin_load_volume \
clean_test:
rm -f $(TEST_DEP) $(TEST_OBJ) $(TEST_TGT)
- rm -f file.txt test_0.stl test_1.stl test_2.stl test_3.stl
+ rm -f file.txt test_0.stl test_1.stl test_2.stl test_3.stl\
+ source_spec.txt abs_cross_sec.txt abs_cross_sec_0.txt\
+ abs_cross_sec_1.txt abs_cross_sec_2.txt abs_cross_sec_3.txt\
+ abs_cross_sec_4.txt
for i in $(TEST_SRC); do rm -f "$$(basename "$${i}" ".c")"; done
diff --git a/doc/star-phor-input.scd b/doc/star-phor-input.scd
@@ -85,6 +85,7 @@ The file format describing a photoreactive system is as follows:
<volume-props> ::= <ka>
| <geometry>
| <sensor>
+ | <prop-rad>
<surface> ::= 'surface:' <name>
[<surface-props> ...]
@@ -113,12 +114,45 @@ The file format describing a photoreactive system is as follows:
| 'm^-1'
| '1/m'
+<prop-rad> ::= 'prop_rad:' <name> <prop-rad-type>
+ <concentration>
+ <cross-sections>
+
+<prop-rad-type> ::= 'SCATTERER'
+
+<concentration> ::= 'concentration:' real <concentration-unit> % conc > 0
+<concentration-unit> ::= 'kg/L'
+ | 'kg.L^-1'
+ | 'kg.m^-3'
+ | 'kg/m^3'
+ | 'mol/L'
+ | 'mol.L^-1'
+ | 'mol.m^-3'
+ | 'mol/m^3'
+ | 'part/L'
+ | 'part.L^-1'
+ | 'part.m^-3'
+ | 'part/m^3'
+<cross-section> ::= 'cross_sections:'
+ <abs-cross-sections>
+<abs-cross-sections> ::= 'abs_cross_sec:' <prop-file> <spec-unit> <cross-sec-unit>
+<prop-file> ::= path % no spaces allowed
+<spec-unit> ::= %TODO
+<cross-sec-unit> ::= % should be compatible with concentration-unit
+ | 'm^2/kg'
+ | 'm^2.kg^-1'
+ | 'm^2/mol'
+ | 'm^2.mol^-1'
+ | 'm^2/part'
+ | 'm^2.part^-1'
+
+
<brdf> ::= 'brdf:' <brdf-type> <reflectivity>
<brdf-type> ::= 'LAMBERT' | 'SPECULAR'
<reflectivity> ::= real % In [0, 1]
<surface-source> ::= 'source:'
- 'flux_density:' <flux-val> <flux-unit>
+ 'flux_density:' <flux-val> <flux-unit> <prop-file> <spec-unit> %TODO <pdf-unit>
'direction:' <direction-distrib>
<flux-val> ::= real % Flux density value > 0
diff --git a/src/sphin.h b/src/sphin.h
@@ -69,6 +69,12 @@ enum sphin_source_surface_direction_distribution_type {
SPHIN_SOURCE_DIRECTION_NONE__
};
+enum sphin_prop_rad_type {
+ SPHIN_PROP_RAD_BOLTZMANN,
+ SPHIN_PROP_RAD_SCATTERER,
+ SPHIN_PROP_RAD_NONE__
+};
+
struct sphin_create_args {
struct logger* logger; /* May be NULL <=> default logger */
struct mem_allocator* allocator; /* NULL <=> use default allocator */
@@ -106,9 +112,10 @@ static const struct sphin_brdf_specular SPHIN_BRDF_SPECULAR_NULL =
struct sphin_source_surface_flux_density {
double flux_density; /* either [umol/m^2/s] or [W/m^2] */
+ struct sphin_spectral_property* emission_spectrum;
enum sphin_photon_unit unit;
};
-#define SPHIN_SOURCE_SURFACE_FLUX_DENSITY_NULL__ {0, SPHIN_PHOTON_UNIT_NONE__}
+#define SPHIN_SOURCE_SURFACE_FLUX_DENSITY_NULL__ {0, NULL, SPHIN_PHOTON_UNIT_NONE__}
static const struct sphin_source_surface_flux_density SPHIN_SOURCE_SURFACE_FLUX_DENSITY_NULL =
SPHIN_SOURCE_SURFACE_FLUX_DENSITY_NULL__;
@@ -181,6 +188,17 @@ static const struct sphin_geometry_descriptor
SPHIN_GEOMETRY_DESCRIPTOR_NULL =
SPHIN_GEOMETRY_DESCRIPTOR_NULL__;
+struct sphin_spectral_property_descriptor {
+ const char* filename;
+ const double* wavelengths; /* In nm */
+ const double* values; /* Corresponding value to each of the wavelengths */
+ size_t data_count; /* Number of elements in the arrays */
+};
+#define SPHIN_SPECTRAL_PROPERTY_DESCRIPTOR_NULL__ {NULL, NULL, NULL, 0}
+static const struct sphin_spectral_property_descriptor
+ SPHIN_SPECTRAL_PROPERTY_DESCRIPTOR_NULL =
+ SPHIN_SPECTRAL_PROPERTY_DESCRIPTOR_NULL__;
+
/* Forward declarations
*
* These structures are declared but remain opaque to users of the
@@ -193,9 +211,12 @@ struct sphin; /* Library handler */
struct sphin_config; /* Physical configuration */
struct sphin_brdf;
struct sphin_geometry;
+struct sphin_prop_rad;
+struct sphin_scatterer;
struct sphin_sensor_surface;
struct sphin_sensor_volume;
struct sphin_source_surface;
+struct sphin_spectral_property;
struct sphin_surface;
struct sphin_volume;
@@ -295,6 +316,54 @@ sphin_volume_get_geometry
size_t igeometry,
struct sphin_geometry** geometry);
+SPHIN_API res_T
+sphin_volume_get_prop_rad_count
+ (struct sphin_volume* volume,
+ size_t* prop_rad_count);
+
+SPHIN_API res_T
+sphin_volume_get_prop_rad
+ (struct sphin_volume* volume,
+ size_t iprop_rad,
+ struct sphin_prop_rad** prop_rad);
+
+/*******************************************************************************
+ * API of the radiative properties
+ ******************************************************************************/
+SPHIN_API res_T
+sphin_prop_rad_ref_get
+ (struct sphin_prop_rad* prop_rad);
+
+SPHIN_API res_T
+sphin_prop_rad_ref_put
+ (struct sphin_prop_rad* prop_rad);
+
+SPHIN_API res_T
+sphin_prop_rad_get_scatterer
+ (struct sphin_prop_rad* prop_rad,
+ struct sphin_scatterer** scatterer);
+
+/*******************************************************************************
+ * API of the scatterer
+ ******************************************************************************/
+SPHIN_API res_T
+sphin_scatterer_ref_get
+ (struct sphin_scatterer* scatterer);
+
+SPHIN_API res_T
+sphin_scatterer_ref_put
+ (struct sphin_scatterer* scatterer);
+
+SPHIN_API res_T
+sphin_scatterer_get_concentration
+ (struct sphin_scatterer* scatterer,
+ double* concentration);
+
+SPHIN_API res_T
+sphin_scatterer_get_abs_cross_sec
+ (struct sphin_scatterer* scatterer,
+ struct sphin_spectral_property** abs_cross_sec);
+
/*******************************************************************************
* API of the surface
******************************************************************************/
@@ -427,6 +496,27 @@ sphin_geometry_get_desc
(const struct sphin_geometry* geometry,
struct sphin_geometry_descriptor* desc);
-END_DECLS
+/*******************************************************************************
+ * API of the spectral_property
+ ******************************************************************************/
+SPHIN_API res_T
+sphin_spectral_property_ref_get
+ (struct sphin_spectral_property* property);
+
+SPHIN_API res_T
+sphin_spectral_property_ref_put
+ (struct sphin_spectral_property* property);
+SPHIN_API res_T
+sphin_spectral_property_get_desc
+ (const struct sphin_spectral_property* property,
+ struct sphin_spectral_property_descriptor* desc);
+
+SPHIN_API res_T
+sphin_spectral_property_interpolate_at_wavelength
+ (const struct sphin_spectral_property* property,
+ double wavelength,
+ double* value);
+
+END_DECLS
#endif /* SPHIN_H */
diff --git a/src/sphin_prop_rad.c b/src/sphin_prop_rad.c
@@ -0,0 +1,210 @@
+/* Copyright (C) 2024-2025 Centre National de la Recherche Scientifique
+ * Copyright (C) 2024-2025 Clermont Auvergne INP
+ * Copyright (C) 2024-2025 INSA Lyon
+ * Copyright (C) 2024-2025 Institut Mines Télécom Albi-Carmaux
+ * Copyright (C) 2024-2025 Institut National Polytechnique de Toulouse
+ * Copyright (C) 2024-2025 |Méso|Star> (contact@meso-star.com)
+ * Copyright (C) 2024-2025 PhotonLyX (info@photonlyx.com)
+ * Copyright (C) 2024-2025 Université de Lorraine
+ * Copyright (C) 2024-2025 Université Paul Sabatier
+ * Copyright (C) 2024-2025 Université Toulouse - Jean Jaurès
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+#define _POSIX_C_SOURCE 200112L /* for strtok_r support */
+
+#include "sphin.h"
+#include "sphin_c.h"
+#include "sphin_config.h"
+#include "sphin_prop_rad.h"
+#include "sphin_scatterer.h"
+#include "sphin_volume.h"
+
+#include <rsys/cstr.h> /* str_to_double */
+#include <rsys/ref_count.h>
+#include <rsys/str.h>
+#include <rsys/text_reader.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+prop_rad_create
+ (struct sphin* sphin,
+ const char* name,
+ struct sphin_prop_rad** out_prop_rad)
+{
+ struct sphin_prop_rad* prop_rad = NULL;
+ res_T res = RES_OK;
+
+ ASSERT('\0' != name[0]); /* Name can't be empty */
+ ASSERT(NULL != name);
+ ASSERT(NULL != out_prop_rad);
+ ASSERT(NULL != sphin);
+
+ prop_rad = MEM_CALLOC(sphin->allocator, 1, sizeof(struct sphin_prop_rad));
+ if (NULL == prop_rad) { res = RES_MEM_ERR; goto error; }
+ ref_init(&prop_rad->ref);
+ SPHIN(ref_get(sphin));
+ prop_rad->sphin = sphin;
+ prop_rad->scatterer = NULL;
+
+ str_init(sphin->allocator, &prop_rad->name);
+ res = str_set(&prop_rad->name, name);
+ if (RES_OK != res) { goto error; }
+
+exit:
+ *out_prop_rad = prop_rad;
+ return res;
+error:
+ if (NULL != prop_rad) {
+ SPHIN(prop_rad_ref_put(prop_rad));
+ prop_rad = NULL;
+ }
+ goto exit;
+}
+
+static void
+release_prop_rad
+ (ref_T* address)
+{
+ struct sphin_prop_rad* prop_rad = NULL;
+ struct sphin* sphin = NULL;
+
+ ASSERT(NULL != address);
+
+ prop_rad = CONTAINER_OF(address, struct sphin_prop_rad, ref);
+ str_release(&prop_rad->name);
+ if (NULL != prop_rad->scatterer) {
+ SPHIN(scatterer_ref_put(prop_rad->scatterer));
+ }
+ sphin = prop_rad->sphin;
+ MEM_RM(sphin->allocator, prop_rad);
+ SPHIN(ref_put(sphin));
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+parse_prop_rad
+ (struct sphin_volume* volume,
+ struct txtrdr* txtrdr,
+ char* value)
+{
+ char* tokens[2] = {NULL};
+ char* name = NULL;
+ char* prop_rad_type_str = NULL;
+ char* token = NULL;
+ char* token_ptr = NULL;
+ size_t token_count = 0;
+ enum sphin_prop_rad_type prop_rad_type = SPHIN_PROP_RAD_NONE__;
+ struct sphin_prop_rad* prop_rad = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != volume);
+ ASSERT(NULL != txtrdr);
+
+ if (NULL == value) { res = RES_BAD_ARG; goto error; }
+
+ token = strtok_r(value, " \t", &token_ptr);
+ while (token != NULL && token_count < 5) {
+ tokens[token_count++] = trim_keyword(token);
+ token = strtok_r(NULL, " \t", &token_ptr);
+ }
+
+ name = tokens[0];
+ if (NULL == name) { res = RES_BAD_ARG; goto error; }
+ prop_rad_type_str = tokens[1];
+ if (NULL == prop_rad_type_str) { res = RES_BAD_ARG; goto error; }
+
+ if (NULL == prop_rad_type_str) {
+ prop_rad_type = SPHIN_PROP_RAD_BOLTZMANN;
+ } else if ( 0 == strcmp(prop_rad_type_str, "SCATTERER")) {
+ prop_rad_type = SPHIN_PROP_RAD_SCATTERER;
+ } else {
+ res = RES_BAD_ARG; goto error;
+ }
+
+ res = prop_rad_create(volume->sphin, name, &prop_rad);
+ if (RES_OK != res) { goto error; }
+
+ switch(prop_rad_type) {
+ case SPHIN_PROP_RAD_BOLTZMANN:
+ break;
+ /* TODO Parse prop rad coefficients */
+
+ case SPHIN_PROP_RAD_SCATTERER:
+ res = parse_scatterer(prop_rad, txtrdr);
+ if (RES_OK != res) { goto error; }
+ break;
+
+ case SPHIN_PROP_RAD_NONE__:
+ res = RES_BAD_ARG;
+ goto error;
+
+ default:
+ FATAL("Unreachable code\n"); break;
+ }
+
+ res = darray_sphin_prop_rad_ptr_push_back(&volume->prop_rads, &prop_rad);
+ if (RES_OK != res) { goto error; }
+
+exit:
+ return res;
+error:
+ if (NULL != prop_rad) {
+ SPHIN(prop_rad_ref_put(prop_rad));
+ prop_rad = NULL;
+ }
+ goto exit;
+}
+
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+sphin_prop_rad_ref_get
+ (struct sphin_prop_rad* prop_rad)
+{
+ if (NULL == prop_rad) {
+ return RES_BAD_ARG;
+ }
+ ref_get(&prop_rad->ref);
+ return RES_OK;
+}
+
+res_T
+sphin_prop_rad_ref_put
+ (struct sphin_prop_rad* prop_rad)
+{
+ if (NULL == prop_rad) {
+ return RES_BAD_ARG;
+ }
+ ref_put(&prop_rad->ref, release_prop_rad);
+ return RES_OK;
+}
+
+res_T
+sphin_prop_rad_get_scatterer
+ (struct sphin_prop_rad* prop_rad,
+ struct sphin_scatterer** scatterer)
+{
+ if (NULL == prop_rad
+ || NULL == scatterer) {
+ return RES_BAD_ARG;
+ }
+ *scatterer = prop_rad->scatterer;
+ return RES_OK;
+}
diff --git a/src/sphin_prop_rad.h b/src/sphin_prop_rad.h
@@ -0,0 +1,62 @@
+/* Copyright (C) 2024-2025 Centre National de la Recherche Scientifique
+ * Copyright (C) 2024-2025 Clermont Auvergne INP
+ * Copyright (C) 2024-2025 INSA Lyon
+ * Copyright (C) 2024-2025 Institut Mines Télécom Albi-Carmaux
+ * Copyright (C) 2024-2025 Institut National Polytechnique de Toulouse
+ * Copyright (C) 2024-2025 |Méso|Star> (contact@meso-star.com)
+ * Copyright (C) 2024-2025 PhotonLyX (info@photonlyx.com)
+ * Copyright (C) 2024-2025 Université de Lorraine
+ * Copyright (C) 2024-2025 Université Paul Sabatier
+ * Copyright (C) 2024-2025 Université Toulouse - Jean Jaurès
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SPHIN_PROP_RAD_H
+#define SPHIN_PROP_RAD_H
+
+#include "sphin.h"
+
+#include <rsys/str.h>
+#include <rsys/dynamic_array.h>
+#include <rsys/ref_count.h>
+#include <rsys/rsys.h>
+
+struct txtrdr;
+
+struct sphin_prop_rad{
+ struct str name;
+ /* TODO struct coeffs* coeffs */
+ struct sphin_scatterer* scatterer;
+
+ struct sphin* sphin;
+ ref_T ref;
+};
+
+extern LOCAL_SYM res_T
+parse_prop_rad
+ (struct sphin_volume* volume,
+ struct txtrdr* txtrdr,
+ char* value);
+
+/* Generate the dynamic array of pointers for the structure sphin_prop_rad
+ *
+ * Refer to sphin_geometry.h to a more detailed explanation on this */
+#define DARRAY_NAME sphin_prop_rad_ptr/* Prefix for api functions and structures:
+ darray_sphin_prop_rad_ptr */
+#define DARRAY_DATA struct sphin_prop_rad*
+
+/* Generate the code */
+#include <rsys/dynamic_array.h>
+
+#endif /* SPHIN_PROP_RAD_H */
diff --git a/src/sphin_scatterer.c b/src/sphin_scatterer.c
@@ -0,0 +1,407 @@
+/* Copyright (C) 2024-2025 Centre National de la Recherche Scientifique
+ * Copyright (C) 2024-2025 Clermont Auvergne INP
+ * Copyright (C) 2024-2025 INSA Lyon
+ * Copyright (C) 2024-2025 Institut Mines Télécom Albi-Carmaux
+ * Copyright (C) 2024-2025 Institut National Polytechnique de Toulouse
+ * Copyright (C) 2024-2025 |Méso|Star> (contact@meso-star.com)
+ * Copyright (C) 2024-2025 PhotonLyX (info@photonlyx.com)
+ * Copyright (C) 2024-2025 Université de Lorraine
+ * Copyright (C) 2024-2025 Université Paul Sabatier
+ * Copyright (C) 2024-2025 Université Toulouse - Jean Jaurès
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+#define _POSIX_C_SOURCE 200112L /* for strtok_r support */
+
+#include "sphin.h"
+#include "sphin_c.h"
+#include "sphin_config.h"
+#include "sphin_prop_rad.h"
+#include "sphin_scatterer.h"
+#include "sphin_spectral_property.h"
+
+#include <rsys/cstr.h> /* str_to_double */
+#include <rsys/ref_count.h>
+#include <rsys/str.h>
+#include <rsys/text_reader.h>
+
+struct sphin_scatterer{
+ double concentration;
+ struct sphin_spectral_property* abs_cross_sec;
+
+ /* TODO
+ * struct sca_cross_sec sca_cross_sec;
+ * struct mie* mie;
+ * struct phase_fn* phase_fn;
+ * */
+
+ struct sphin* sphin;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+scatterer_create
+ (struct sphin* sphin,
+ struct sphin_scatterer** out_scatterer)
+{
+ struct sphin_scatterer* scatterer = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != out_scatterer);
+ ASSERT(NULL != sphin);
+
+ scatterer = MEM_CALLOC(sphin->allocator, 1, sizeof(struct sphin_scatterer));
+ if (NULL == scatterer) { res = RES_MEM_ERR; goto error; }
+ ref_init(&scatterer->ref);
+ SPHIN(ref_get(sphin));
+ scatterer->sphin = sphin;
+ scatterer->abs_cross_sec = NULL;
+
+exit:
+ *out_scatterer = scatterer;
+ return res;
+error:
+ if (NULL != scatterer) {
+ SPHIN(scatterer_ref_put(scatterer));
+ scatterer = NULL;
+ }
+ goto exit;
+}
+
+static void
+release_scatterer
+ (ref_T* address)
+{
+ struct sphin_scatterer* scatterer = NULL;
+ struct sphin* sphin = NULL;
+
+ ASSERT(NULL != address);
+
+ scatterer = CONTAINER_OF(address, struct sphin_scatterer, ref);
+ sphin = scatterer->sphin;
+ if (NULL != scatterer->abs_cross_sec) {
+ SPHIN(spectral_property_ref_put(scatterer->abs_cross_sec));
+ }
+ MEM_RM(sphin->allocator, scatterer);
+ SPHIN(ref_put(sphin));
+}
+
+static res_T
+parse_concentration
+ (struct sphin_scatterer* scatterer,
+ struct txtrdr* txtrdr,
+ char* value)
+{
+ char* concentration_val = NULL;
+ char* concentration_unit = NULL;
+ char* token = NULL;
+ char* token_ptr = NULL;
+ res_T res = RES_OK;
+ double concentration;
+
+ ASSERT(NULL != scatterer);
+ ASSERT(NULL != txtrdr);
+
+ if (NULL == value){ res = RES_BAD_ARG; goto error; }
+
+ /* Parse concentration value */
+ concentration_val = strtok_r(value, " \t", &token_ptr);
+ res = cstr_to_double(concentration_val, &concentration);
+ if (RES_OK != res) { goto error; }
+ if (concentration < 0) { res = RES_BAD_ARG; goto error; }
+
+ /* Parse concentration unit and convert its value to unit/m^-3 accordingly,
+ * where unit is kg, mol or part */
+ concentration_unit = strtok_r(NULL, " \t", &token_ptr);
+
+ if (0 == strcmp(concentration_unit, "kg/m^3")
+ || 0 == strcmp(concentration_unit, "kg.m^-3")
+ || 0 == strcmp(concentration_unit, "mol/m^3")
+ || 0 == strcmp(concentration_unit, "mol.m^-3")
+ || 0 == strcmp(concentration_unit, "part/m^3")
+ || 0 == strcmp(concentration_unit, "part.m^-3")
+ ) {
+ scatterer->concentration = concentration * 1;
+ } else
+ if (0 == strcmp(concentration_unit, "kg/L")
+ || 0 == strcmp(concentration_unit, "kg.L^-1")
+ || 0 == strcmp(concentration_unit, "mol/L")
+ || 0 == strcmp(concentration_unit, "mol.L^-1")
+ || 0 == strcmp(concentration_unit, "part/L")
+ || 0 == strcmp(concentration_unit, "part.L^-1")
+ ) {
+ scatterer->concentration = concentration * 1000;
+ }
+ else { res = RES_BAD_ARG; goto error; }
+
+ /* Assert that there isnt anything but the concentration value and its unit
+ * in the line */
+ token = strtok_r(NULL, " \t", &token_ptr);
+ if (NULL != token) { res = RES_BAD_ARG; goto error; }
+
+ res = txtrdr_read_line(txtrdr);
+ if (RES_OK != res) { goto error; }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_abs_cross_sec
+ (struct sphin_scatterer* scatterer,
+ struct txtrdr* txtrdr,
+ char* value)
+{
+ char* tokens[3] = {NULL}; /* filename, spectral_unit, property_unit */
+ char* filename = NULL;
+ char* spectral_unit = NULL;
+ char* property_unit = NULL;
+ char* token = NULL;
+ char* token_ptr = NULL;
+ size_t token_count = 0;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != scatterer);
+ ASSERT(NULL != txtrdr);
+
+ if (NULL == value) { res = RES_BAD_ARG; goto error; }
+
+ token = strtok_r(value, " \t", &token_ptr);
+ while (token != NULL && token_count < 5) {
+ tokens[token_count++] = trim_keyword(token);
+ token = strtok_r(NULL, " \t", &token_ptr);
+ }
+
+ /* The number of tokens must be 3 (filename and each one of the units) */
+ if (token_count != 3) { res = RES_BAD_ARG; goto error; }
+
+ filename = tokens[0];
+ if (NULL == filename) { res = RES_BAD_ARG; goto error; }
+ spectral_unit = tokens[1];
+ if (NULL == spectral_unit) { res = RES_BAD_ARG; goto error; }
+ property_unit = tokens[2];
+ if (NULL == property_unit) { res = RES_BAD_ARG; goto error; }
+
+ res = parse_spectral_property
+ (scatterer->sphin, filename, &scatterer->abs_cross_sec);
+ if (RES_OK != res) { goto error; }
+
+ res = convert_wavelengths
+ (&scatterer->abs_cross_sec->wavelengths, spectral_unit);
+ if (RES_OK != res) { goto error; }
+
+ /* TODO Convert property according to unit */
+
+ res = txtrdr_read_line(txtrdr);
+ if (RES_OK != res) { goto error; }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_tabulated_cross_sections
+ (struct sphin_scatterer* scatterer,
+ struct txtrdr* txtrdr)
+{
+ char* keyword = NULL;
+ char* token = NULL;
+ char* token_ptr = NULL;
+ struct str line;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != scatterer);
+ ASSERT(NULL != txtrdr);
+
+ str_init(scatterer->sphin->allocator, &line);
+
+ res = txtrdr_read_line(txtrdr);
+ if (RES_OK != res) { goto error; }
+
+ while (NULL != txtrdr_get_line(txtrdr)) {
+ res = str_set(&line, txtrdr_get_cline(txtrdr));
+ if (RES_OK != res) { goto error; }
+ /* parse keyword */
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ if (NULL == token){ res = RES_BAD_ARG; goto error; }
+ keyword = trim_keyword(token);
+ if (NULL == keyword){ res = RES_BAD_ARG; goto error; }
+ if (0 == strcmp(keyword, "abs_cross_sec")){
+ res = parse_abs_cross_sec(scatterer, txtrdr, token_ptr);
+ }
+ /* TODO
+ * if (0 == strcmp(keyword, "sca_cross_sec")){
+ * res = parse_sca_cross_sec(scatterer, txtrdr, token_ptr);
+ * }
+ * if (0 == strcmp(keyword, "phase_fn")){
+ * res = parse_phase_fn(scatterer, txtrdr, token_ptr);
+ * }
+ */
+ else {
+ break;
+ }
+ if (RES_OK != res) { goto error; }
+ }
+
+exit:
+ str_release(&line);
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_cross_sections
+ (struct sphin_scatterer* scatterer,
+ struct txtrdr* txtrdr,
+ char* value)
+{
+ char* type = NULL;
+ char* token_ptr = NULL;
+ res_T res = RES_OK;
+
+ type = strtok_r(value, "\t ", &token_ptr);
+
+ if (NULL == type) {
+ /* Parse the cross sections as imposed in the input file */
+ res = parse_tabulated_cross_sections(scatterer, txtrdr);
+ if (RES_OK != res) { goto error; }
+ } else if (0 == strcmp(type, "MIE")) {
+ /* TODO Scatterer type is mie */
+ } else {
+ res = RES_BAD_ARG; goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+parse_scatterer
+ (struct sphin_prop_rad* prop_rad,
+ struct txtrdr* txtrdr)
+{
+ struct sphin_scatterer* scatterer = NULL;
+ struct str line;
+ char* keyword = NULL;
+ char* token = NULL;
+ char* token_ptr = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != prop_rad);
+ ASSERT(NULL != txtrdr);
+
+ str_init(prop_rad->sphin->allocator, &line);
+
+ res = scatterer_create(prop_rad->sphin, &scatterer);
+ if (RES_OK != res) { goto error; }
+
+ res = txtrdr_read_line(txtrdr);
+ if (RES_OK != res) { goto error; }
+
+ while (NULL != txtrdr_get_line(txtrdr)) {
+ res = str_set(&line, txtrdr_get_cline(txtrdr));
+ if (RES_OK != res) { goto error; }
+ /* parse keyword */
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ if (NULL == token){ res = RES_BAD_ARG; goto error; }
+ keyword = trim_keyword(token);
+ if (NULL == keyword){ res = RES_BAD_ARG; goto error; }
+
+ /* parse value */
+ token = token_ptr;
+ if (0 == strcmp(keyword, "concentration")){
+ res = parse_concentration(scatterer, txtrdr, token);
+ }
+ else if (0 == strcmp(keyword, "cross_sections")){
+ /* TODO Parse cross sections */
+ res = parse_cross_sections(scatterer, txtrdr, token);
+ }
+ else {
+ break;
+ }
+ if (RES_OK != res) { goto error; }
+ }
+
+exit:
+ prop_rad->scatterer = scatterer;
+ str_release(&line);
+ return res;
+error:
+ if (NULL != scatterer) {
+ SPHIN(scatterer_ref_put(scatterer));
+ scatterer = NULL;
+ }
+ goto exit;
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+sphin_scatterer_ref_get
+ (struct sphin_scatterer* scatterer)
+{
+ if (NULL == scatterer) {
+ return RES_BAD_ARG;
+ }
+ ref_get(&scatterer->ref);
+ return RES_OK;
+}
+
+res_T
+sphin_scatterer_ref_put
+ (struct sphin_scatterer* scatterer)
+{
+ if (NULL == scatterer) {
+ return RES_BAD_ARG;
+ }
+ ref_put(&scatterer->ref, release_scatterer);
+ return RES_OK;
+}
+
+res_T
+sphin_scatterer_get_concentration
+ (struct sphin_scatterer* scatterer,
+ double* concentration)
+{
+ if (NULL == scatterer || NULL == concentration) {
+ return RES_BAD_ARG;
+ }
+ *concentration = scatterer->concentration;
+ return RES_OK;
+}
+
+res_T
+sphin_scatterer_get_abs_cross_sec
+ (struct sphin_scatterer* scatterer,
+ struct sphin_spectral_property** abs_cross_sec)
+{
+ if (NULL == scatterer || NULL == abs_cross_sec) {
+ return RES_BAD_ARG;
+ }
+ *abs_cross_sec = scatterer->abs_cross_sec;
+ return RES_OK;
+}
diff --git a/src/sphin_scatterer.h b/src/sphin_scatterer.h
@@ -0,0 +1,41 @@
+/* Copyright (C) 2024-2025 Centre National de la Recherche Scientifique
+ * Copyright (C) 2024-2025 Clermont Auvergne INP
+ * Copyright (C) 2024-2025 INSA Lyon
+ * Copyright (C) 2024-2025 Institut Mines Télécom Albi-Carmaux
+ * Copyright (C) 2024-2025 Institut National Polytechnique de Toulouse
+ * Copyright (C) 2024-2025 |Méso|Star> (contact@meso-star.com)
+ * Copyright (C) 2024-2025 PhotonLyX (info@photonlyx.com)
+ * Copyright (C) 2024-2025 Université de Lorraine
+ * Copyright (C) 2024-2025 Université Paul Sabatier
+ * Copyright (C) 2024-2025 Université Toulouse - Jean Jaurès
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SPHIN_SCATTERER_H
+#define SPHIN_SCATTERER_H
+
+#include "sphin.h"
+
+#include <rsys/dynamic_array.h>
+#include <rsys/rsys.h>
+
+struct sphin_prop_rad;
+struct txtrdr;
+
+extern LOCAL_SYM res_T
+parse_scatterer
+ (struct sphin_prop_rad* prop_rad,
+ struct txtrdr* txtrdr);
+
+#endif /* SPHIN_SCATTERER_H */
diff --git a/src/sphin_source_surface.c b/src/sphin_source_surface.c
@@ -26,6 +26,7 @@
#include "sphin.h"
#include "sphin_c.h"
#include "sphin_config.h"
+#include "sphin_spectral_property.h"
#include "sphin_source_surface.h"
#include "sphin_surface.h"
@@ -95,6 +96,9 @@ release_source_surface
source = CONTAINER_OF(address, struct sphin_source_surface, ref);
str_release(&source->name);
+ if (NULL != source->flux_density.emission_spectrum) {
+ SPHIN(spectral_property_ref_put(source->flux_density.emission_spectrum));
+ }
sphin = source->sphin;
MEM_RM(sphin->allocator, source);
SPHIN(ref_put(sphin));
@@ -244,61 +248,119 @@ error:
}
static res_T
-parse_flux_density
+parse_flux_density_unit
(struct sphin_source_surface* source,
struct txtrdr* txtrdr,
char* value)
{
- char* flux_density_unit = NULL;
- char* str_flux_density = NULL;
- char* token_ptr = NULL;
- double flux_density;
res_T res = RES_OK;
ASSERT(NULL != source);
ASSERT(NULL != txtrdr);
+ ASSERT(NULL != value);
- if (NULL == value) { res = RES_BAD_ARG; goto error; }
-
- /* Parse flux density value */
- str_flux_density = strtok_r(value, " \t", &token_ptr);
- res = cstr_to_double(str_flux_density, &flux_density);
- if (RES_OK != res) { res = RES_BAD_ARG; goto error; }
- if (flux_density < 0) { res = RES_BAD_ARG; goto error; }
-
- source->flux_density.flux_density = flux_density;
-
- /* Parse unit */
- flux_density_unit = strtok_r(NULL, " \t", &token_ptr);
- if (NULL == flux_density_unit){ res = RES_BAD_ARG; goto error; }
+ (void)txtrdr;
- /* Parse kinetic flux density unit */
- if (0 == strcmp(flux_density_unit, "mol/m^2/s")
- || 0 == strcmp(flux_density_unit, "mol.m^-2.s^-1")) {
+ if (0 == strcmp(value, "mol/m^2/s")
+ || 0 == strcmp(value, "mol.m^-2.s^-1")) {
source->flux_density.unit= SPHIN_PHOTON_UNIT_MOL;
source->flux_density.flux_density *= 1e6; /* From mol to umol */
}
- else if (0 == strcmp(flux_density_unit, "umol/m^2/s")
- || 0 == strcmp(flux_density_unit, "umol.m^-2.s^-1")) {
+ else if (0 == strcmp(value, "umol/m^2/s")
+ || 0 == strcmp(value, "umol.m^-2.s^-1")) {
source->flux_density.unit= SPHIN_PHOTON_UNIT_MOL;
source->flux_density.flux_density *= 1; /* No conversion */
}
/* Energy flux density energy */
- else if (0 == strcmp(flux_density_unit, "mW/m^2")
- || 0 == strcmp(flux_density_unit, "mW.m^-2")) {
+ else if (0 == strcmp(value, "mW/m^2")
+ || 0 == strcmp(value, "mW.m^-2")) {
source->flux_density.unit= SPHIN_PHOTON_UNIT_JOULE;
source->flux_density.flux_density *= 1e-3; /* From mWatt to Watt*/
}
- else if (0 == strcmp(flux_density_unit, "W/m^2")
- || 0 == strcmp(flux_density_unit, "W.m^-2")
- || 0 == strcmp(flux_density_unit, "J.m^-2.s^-1")
- || 0 == strcmp(flux_density_unit, "J/m^2/s")) {
+ else if (0 == strcmp(value, "W/m^2")
+ || 0 == strcmp(value, "W.m^-2")
+ || 0 == strcmp(value, "J.m^-2.s^-1")
+ || 0 == strcmp(value, "J/m^2/s")) {
source->flux_density.unit= SPHIN_PHOTON_UNIT_JOULE;
source->flux_density.flux_density *= 1; /* No conversion */
}
else { res = RES_BAD_ARG; goto error; }
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_flux_density
+ (struct sphin_source_surface* source,
+ struct txtrdr* txtrdr,
+ char* value)
+{
+ char* tokens[5] = {NULL}; /* str_flux_density, flux_density_unit, [filename],
+ [spectral_unit], [property_unit] */
+ char* filename = NULL;
+ char* flux_density_unit = NULL;
+ char* property_unit = NULL;
+ char* spectral_unit = NULL;
+ char* str_flux_density = NULL;
+ char* token = NULL;
+ char* token_ptr = NULL;
+ double flux_density = 0;
+ size_t token_count = 0;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != source);
+ ASSERT(NULL != txtrdr);
+
+ if (NULL == value) { res = RES_BAD_ARG; goto error; }
+
+ token = strtok_r(value, " \t", &token_ptr);
+ while (token != NULL && token_count < 5) {
+ tokens[token_count++] = trim_keyword(token);
+ token = strtok_r(NULL, " \t", &token_ptr);
+ }
+
+ /* The number of tokens must be either 2 (gray source) or 5 (read from file) */
+ if (token_count != 2 && token_count != 5) { res = RES_BAD_ARG; goto error; }
+
+ /* Parse flux density value */
+ str_flux_density = tokens[0];
+ res = cstr_to_double(str_flux_density, &flux_density);
+ if (RES_OK != res) { res = RES_BAD_ARG; goto error; }
+ if (flux_density < 0) { res = RES_BAD_ARG; goto error; }
+ source->flux_density.flux_density = flux_density;
+
+ /* Parse flux density unit */
+ flux_density_unit = tokens[1];
+ if (NULL == flux_density_unit){ res = RES_BAD_ARG; goto error; }
+
+ res = parse_flux_density_unit(source, txtrdr, flux_density_unit);
+ if (RES_OK != res){ goto error; }
+
+ if (token_count == 5) {
+ filename = tokens[2];
+ if (NULL == filename){ res = RES_BAD_ARG; goto error; }
+
+ spectral_unit = tokens[3];
+ if (NULL == spectral_unit){ res = RES_BAD_ARG; goto error; }
+
+ property_unit = tokens[4];
+ if (NULL == property_unit){ res = RES_BAD_ARG; goto error; }
+
+ res = parse_spectral_property
+ (source->sphin, filename, &source->flux_density.emission_spectrum);
+ if (RES_OK != res) { goto error; }
+
+ /* wavelengths will be in nm after conversion */
+ res = convert_wavelengths(&source->flux_density.emission_spectrum->wavelengths, spectral_unit);
+ if (RES_OK != res) { goto error; }
+
+ /* TODO Convert and property according to units */
+ }
+
res = txtrdr_read_line(txtrdr);
if (RES_OK != res) { goto error; }
diff --git a/src/sphin_spectral_property.c b/src/sphin_spectral_property.c
@@ -0,0 +1,353 @@
+/* Copyright (C) 2024-2025 Centre National de la Recherche Scientifique
+ * Copyright (C) 2024-2025 Clermont Auvergne INP
+ * Copyright (C) 2024-2025 INSA Lyon
+ * Copyright (C) 2024-2025 Institut Mines Télécom Albi-Carmaux
+ * Copyright (C) 2024-2025 Institut National Polytechnique de Toulouse
+ * Copyright (C) 2024-2025 |Méso|Star> (contact@meso-star.com)
+ * Copyright (C) 2024-2025 PhotonLyX (info@photonlyx.com)
+ * Copyright (C) 2024-2025 Université de Lorraine
+ * Copyright (C) 2024-2025 Université Paul Sabatier
+ * Copyright (C) 2024-2025 Université Toulouse - Jean Jaurès
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+#define _POSIX_C_SOURCE 200112L /* for strtok_r support */
+
+#include "sphin.h"
+#include "sphin_c.h"
+#include "sphin_spectral_property.h"
+
+#include <rsys/algorithm.h>
+#include <rsys/cstr.h>
+#include <rsys/double3.h>
+#include <rsys/text_reader.h>
+#include <star/sstl.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+spectral_property_create
+ (struct sphin* sphin,
+ const char* property_filename,
+ struct sphin_spectral_property** out_property)
+{
+ struct sphin_spectral_property* property = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != sphin);
+ ASSERT(NULL != out_property);
+ ASSERT(NULL != property_filename);
+ ASSERT('\0' != property_filename[0]); /* Name can't be empty */
+
+ property = MEM_CALLOC(sphin->allocator, 1,
+ sizeof(struct sphin_spectral_property));
+ if (NULL == property) { res = RES_MEM_ERR; goto error; }
+
+ /* Init geometry ref counter and init member variables */
+ ref_init(&property->ref);
+ SPHIN(ref_get(sphin));
+ property->sphin = sphin;
+ darray_double_init(sphin->allocator, &property->wavelengths);
+ darray_double_init(sphin->allocator, &property->values);
+
+ str_init(sphin->allocator, &property->property_filename);
+ res = str_set(&property->property_filename, property_filename);
+ if (RES_OK != res) { goto error; }
+
+exit:
+ *out_property = property;
+ return res;
+error:
+ if (NULL != property) {
+ SPHIN(spectral_property_ref_put(property));
+ property = NULL;
+ }
+ goto exit;
+}
+
+static void
+release_spectral_property
+ (ref_T* address)
+{
+ struct sphin_spectral_property* property = NULL;
+ struct sphin* sphin = NULL;
+
+ ASSERT(NULL != address);
+
+ property = CONTAINER_OF(address, struct sphin_spectral_property, ref);
+ str_release(&property->property_filename);
+ darray_double_release(&property->wavelengths);
+ darray_double_release(&property->values);
+ sphin = property->sphin;
+ MEM_RM(sphin->allocator, property);
+ SPHIN(ref_put(sphin));
+}
+
+/* As needed in search_lower_bound, see rsys/algorithm.h */
+static int
+compare_wavelengths
+ (const void* key,
+ const void* element)
+{
+ double target = *(const double*) key;
+ double array_element = *(const double*) element;
+
+ return (target > array_element) - (target < array_element);
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+parse_spectral_property
+ (struct sphin* sphin,
+ char* filename,
+ struct sphin_spectral_property** out_property)
+{
+ char* token = NULL;
+ char* token_ptr = NULL;
+ double wavelength = 0;
+ double wavelength_prev = 0; /* verify that the wavelengths are sorted */
+ double value = 0;
+
+ FILE* stream = NULL;
+ struct sphin_spectral_property* property = NULL;
+ struct txtrdr* txtrdr = NULL;
+ struct str line;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != sphin);
+ ASSERT(NULL != out_property);
+ ASSERT(NULL != filename);
+ ASSERT('\0' != filename[0]); /* filename can't be empty */
+
+ stream = fopen(filename, "r");
+ if (NULL == stream) {
+ ERROR(sphin, "Not possible to open file %s -- %s\n",
+ filename, strerror(errno));
+ res = RES_IO_ERR;
+ goto error;
+ }
+
+ str_init(sphin->allocator, &line);
+
+ res = spectral_property_create(sphin, filename, &property);
+ if (RES_OK != res) { goto error; }
+
+ res = txtrdr_stream(sphin->allocator, stream, filename, '#', &txtrdr);
+ if (RES_OK != res) { goto error; }
+
+ res = txtrdr_read_line(txtrdr);
+ if (RES_OK != res) { goto error; }
+
+ while (NULL != txtrdr_get_line(txtrdr)) {
+ res = str_set(&line, txtrdr_get_cline(txtrdr));
+ if (RES_OK != res) { goto error; }
+
+ /* Parse wavelength */
+ token = strtok_r(str_get(&line), "\t ", &token_ptr);
+ if (NULL == token) { res = RES_BAD_ARG; goto error; }
+ if (NULL == token_ptr) { res = RES_BAD_ARG; goto error; }
+
+ res = cstr_to_double(token, &wavelength);
+ if (RES_OK != res) {
+ ERROR
+ (sphin,
+ "%s: %lu: could not read numerical value of wavelenght\n",
+ filename, txtrdr_get_line_num(txtrdr));
+ goto error; }
+
+ if (wavelength_prev >= wavelength) {
+ res = RES_BAD_ARG;
+ ERROR
+ (sphin,
+ "%s: %lu: Wavelengths in spectral properties must be sorted\n",
+ filename, txtrdr_get_line_num(txtrdr));
+ goto error;
+ }
+
+ res = darray_double_push_back(&property->wavelengths, &wavelength);
+ if (RES_OK != res) { goto error; }
+
+ res = cstr_to_double(token_ptr, &value);
+ if (RES_OK != res) {
+ ERROR
+ (sphin,
+ "%s: %lu: could not read numerical value of property\n",
+ filename, txtrdr_get_line_num(txtrdr));
+ goto error; }
+
+ res = darray_double_push_back(&property->values, &value);
+ if (RES_OK != res) { goto error; }
+
+ wavelength_prev = wavelength;
+
+ res = txtrdr_read_line(txtrdr);
+ if (RES_OK != res) { goto error; }
+ }
+
+exit:
+ fclose(stream);
+ str_release(&line);
+ if (NULL != txtrdr) {
+ txtrdr_ref_put(txtrdr);
+ }
+ *out_property = property;
+ return res;
+error:
+ if( NULL != property) {
+ SPHIN(spectral_property_ref_put(property));
+ property = NULL;
+ }
+ goto exit;
+}
+
+res_T
+convert_wavelengths
+ (struct darray_double* wavelengths,
+ char* spectral_unit)
+{
+ double scaling_factor = 1;
+ int invert = 0;
+ double* wls = NULL;
+ size_t i, wl_count = 0;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != wavelengths);
+ ASSERT(NULL != spectral_unit);
+
+ if (0 == strcmp(spectral_unit, "nm")) {scaling_factor = 1e+00;}
+ else if (0 == strcmp(spectral_unit, "cm")) {scaling_factor = 1e07;}
+ else if (0 == strcmp(spectral_unit, "m")) {scaling_factor = 1e09;}
+ else if (0 == strcmp(spectral_unit, "cm^-1")) {
+ scaling_factor = 1e07;
+ invert = 1;}
+ else if (0 == strcmp(spectral_unit, "1/cm")) {
+ scaling_factor = 1e07;
+ invert = 1;}
+ else { res = RES_BAD_ARG; goto error; }
+
+ wl_count = darray_double_size_get(wavelengths);
+ wls = darray_double_data_get(wavelengths);
+
+ FOR_EACH(i, 0, wl_count) {
+ if (invert) {
+ if (0.0 == wls[i]) { res = RES_BAD_ARG; goto error; }
+ wls[i] = scaling_factor / wls[i];
+ }
+ else {wls[i] = wls[i] * scaling_factor;}
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+sphin_spectral_property_ref_get
+ (struct sphin_spectral_property* property)
+{
+ if (NULL == property) {
+ return RES_BAD_ARG;
+ }
+ ref_get(&property->ref);
+ return RES_OK;
+}
+
+res_T
+sphin_spectral_property_ref_put
+ (struct sphin_spectral_property* property)
+{
+ if (NULL == property) {
+ return RES_BAD_ARG;
+ }
+ ref_put(&property->ref, release_spectral_property);
+ return RES_OK;
+}
+
+res_T
+sphin_spectral_property_get_desc
+ (const struct sphin_spectral_property* property,
+ struct sphin_spectral_property_descriptor* desc)
+{
+ if (NULL == property || NULL == desc) {
+ return RES_BAD_ARG;
+ }
+
+ desc->wavelengths = darray_double_data_get
+ ((struct darray_double*)&property->wavelengths);
+ desc->values = darray_double_data_get
+ ((struct darray_double*)&property->values);
+ desc->filename = (char*)str_get((struct str*)&property->property_filename);
+
+ desc->data_count = darray_double_size_get(&property->wavelengths);
+
+ return RES_OK;
+}
+
+res_T
+sphin_spectral_property_interpolate_at_wavelength
+ (const struct sphin_spectral_property* property,
+ double wavelength,
+ double* value)
+{
+ double wl_min, wl_max; /* Min and max wl values in the whole table */
+ double wl_lower, wl_upper; /* Wl values used in interpolation */
+ double v_lower, v_upper; /* Property values used in interpolation */
+ const double* wavelengths = NULL;
+ const double* values = NULL;
+ const double* upper = NULL;
+
+ size_t data_count, upper_index, lower_index;
+
+ if (NULL == property || NULL == value) {
+ return RES_BAD_ARG;
+ }
+
+ data_count = darray_double_size_get(&property->wavelengths);
+ wavelengths = darray_double_cdata_get(&property->wavelengths);
+ values = darray_double_cdata_get(&property->values);
+ wl_min = wavelengths[0];
+ wl_max = wavelengths[data_count - 1];
+
+ if (wavelength < wl_min || wavelength > wl_max) {
+ return RES_BAD_ARG;
+ }
+
+ upper = search_lower_bound
+ (&wavelength, wavelengths, data_count, sizeof(double), compare_wavelengths);
+
+ if (upper == wavelengths) {
+ /* Exact match is the first element of the wavelengths */
+ *value = values[0];
+
+ return RES_OK;
+ }
+ upper_index = (size_t)(upper - wavelengths);
+ lower_index = upper_index - 1;
+
+ wl_lower = wavelengths[lower_index];
+ wl_upper = wavelengths[upper_index];
+ v_lower = values[lower_index];
+ v_upper = values[upper_index];
+
+ *value = v_lower + (wavelength - wl_lower) *
+ (v_upper - v_lower) / (wl_upper - wl_lower);
+
+ return RES_OK;
+}
diff --git a/src/sphin_spectral_property.h b/src/sphin_spectral_property.h
@@ -0,0 +1,61 @@
+/* Copyright (C) 2024-2025 Centre National de la Recherche Scientifique
+ * Copyright (C) 2024-2025 Clermont Auvergne INP
+ * Copyright (C) 2024-2025 INSA Lyon
+ * Copyright (C) 2024-2025 Institut Mines Télécom Albi-Carmaux
+ * Copyright (C) 2024-2025 Institut National Polytechnique de Toulouse
+ * Copyright (C) 2024-2025 |Méso|Star> (contact@meso-star.com)
+ * Copyright (C) 2024-2025 PhotonLyX (info@photonlyx.com)
+ * Copyright (C) 2024-2025 Université de Lorraine
+ * Copyright (C) 2024-2025 Université Paul Sabatier
+ * Copyright (C) 2024-2025 Université Toulouse - Jean Jaurès
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SPHIN_SPECTRAL_PROPERTY_H
+#define SPHIN_SPECTRAL_PROPERTY_H
+
+#include "sphin.h"
+
+#include <rsys/dynamic_array.h>
+#include <rsys/dynamic_array_double.h>
+#include <rsys/dynamic_array_size_t.h>
+#include <rsys/ref_count.h>
+#include <rsys/str.h>
+
+struct mem_allocator;
+
+struct sphin_spectral_property {
+ /* Property data */
+ struct darray_double wavelengths;
+ struct darray_double values;
+
+ struct str property_filename;
+
+ struct sphin* sphin;
+ ref_T ref;
+};
+
+extern LOCAL_SYM res_T
+parse_spectral_property
+ (struct sphin* sphin,
+ char* filename,
+ struct sphin_spectral_property** out_property);
+
+extern LOCAL_SYM res_T
+convert_wavelengths
+ (struct darray_double* wavelengths,
+ char* spectral_unit);
+
+
+#endif /* SPHIN_SPECTRAL_PROPERTY_H */
diff --git a/src/sphin_volume.c b/src/sphin_volume.c
@@ -27,6 +27,7 @@
#include "sphin_c.h"
#include "sphin_config.h"
#include "sphin_geometry.h"
+#include "sphin_prop_rad.h"
#include "sphin_sensor_volume.h"
#include "sphin_volume.h"
@@ -35,15 +36,6 @@
#include <rsys/str.h>
#include <rsys/text_reader.h>
-struct sphin_volume {
- struct str name;
- double ka;
- struct darray_sphin_geometry_ptr geometries; /* dynamic array of struct geometry. see
- rsys/dynamic_array.h */
- struct sphin_sensor_volume* sensor_volume;
- struct sphin* sphin;
- ref_T ref;
-};
/*******************************************************************************
* Helper functions
@@ -73,6 +65,7 @@ volume_create
str_init(sphin->allocator, &volume->name);
darray_sphin_geometry_ptr_init(sphin->allocator, &volume->geometries);
+ darray_sphin_prop_rad_ptr_init(sphin->allocator, &volume->prop_rads);
res = str_set(&volume->name, name);
if (RES_OK != res) { goto error; }
@@ -91,8 +84,9 @@ static void
release_volume
(ref_T* address)
{
- size_t i, ngeometries;
+ size_t i, ngeometries, nprop_rads;
struct sphin_geometry** geometries = NULL;
+ struct sphin_prop_rad** prop_rads = NULL;
struct sphin* sphin = NULL;
struct sphin_sensor_volume* sensor_volume = NULL;
struct sphin_volume* volume = NULL;
@@ -101,11 +95,14 @@ release_volume
volume = CONTAINER_OF(address, struct sphin_volume, ref);
str_release(&volume->name);
- /* Retrieve the number of geometries associated with the volume and decrease
- * the reference counter of each one of them */
+ /* Retrieve the number of geometries and prop rads associated with the volume
+ * and decrease the reference counter of each one of them */
ngeometries = darray_sphin_geometry_ptr_size_get(&volume->geometries);
geometries = darray_sphin_geometry_ptr_data_get(&volume->geometries);
+ nprop_rads = darray_sphin_prop_rad_ptr_size_get(&volume->prop_rads);
+ prop_rads = darray_sphin_prop_rad_ptr_data_get(&volume->prop_rads);
+
/* Put references for each one of the geometries */
FOR_EACH(i, 0, ngeometries) {
SPHIN(geometry_ref_put(geometries[i]));
@@ -114,6 +111,12 @@ release_volume
darray_sphin_geometry_ptr_release(&volume->geometries);
}
+ FOR_EACH(i, 0, nprop_rads) {
+ SPHIN(prop_rad_ref_put(prop_rads[i]));
+ }
+ if (NULL != prop_rads) {
+ darray_sphin_prop_rad_ptr_release(&volume->prop_rads);
+ }
sphin = volume->sphin;
sensor_volume = volume->sensor_volume;
MEM_RM(sphin->allocator, volume);
@@ -251,9 +254,13 @@ parse_volume
if (0 == strcmp(keyword, "geometry")){
res = parse_geometry(volume, txtrdr, token);
}
+ /* TODO Remove ka from volumes */
else if (0 == strcmp(keyword, "ka")){
res = parse_ka(volume, txtrdr, token);
}
+ else if (0 == strcmp(keyword, "prop_rad")){
+ res = parse_prop_rad(volume, txtrdr, token);
+ }
else if (0 == strcmp(keyword, "sensor")) {
/* for the sensor_volume,
* strings are not allowed in the input file after the :
@@ -368,3 +375,38 @@ sphin_volume_get_geometry
return RES_OK;
}
+
+res_T
+sphin_volume_get_prop_rad_count
+ (struct sphin_volume* volume,
+ size_t* prop_rad_count)
+{
+ if (NULL == volume || NULL == prop_rad_count) {
+ return RES_BAD_ARG;
+ }
+
+ *prop_rad_count = darray_sphin_prop_rad_ptr_size_get(&volume->prop_rads);
+
+ return RES_OK;
+}
+
+res_T
+sphin_volume_get_prop_rad
+ (struct sphin_volume* volume,
+ size_t iprop_rad,
+ struct sphin_prop_rad** prop_rad)
+{
+ size_t prop_rad_count;
+ if (NULL == volume || NULL == prop_rad) {
+ return RES_BAD_ARG;
+ }
+
+ prop_rad_count = darray_sphin_prop_rad_ptr_size_get(&volume->prop_rads);
+ if (iprop_rad > prop_rad_count) {
+ return RES_BAD_ARG;
+ }
+
+ *prop_rad = darray_sphin_prop_rad_ptr_data_get(&volume->prop_rads)[iprop_rad];
+
+ return RES_OK;
+}
diff --git a/src/sphin_volume.h b/src/sphin_volume.h
@@ -25,12 +25,28 @@
#ifndef SPHIN_VOLUME_H
#define SPHIN_VOLUME_H
+#include "sphin_geometry.h"
+#include "sphin_prop_rad.h"
+
+#include <rsys/ref_count.h>
+#include <rsys/str.h>
#include <rsys/dynamic_array.h>
#include <rsys/rsys.h>
struct sphin_config;
struct txtrdr;
+struct sphin_volume {
+ struct str name;
+ double ka;
+ struct darray_sphin_geometry_ptr geometries; /* dynamic array of struct geometry. see
+ rsys/dynamic_array.h */
+ struct darray_sphin_prop_rad_ptr prop_rads;
+ struct sphin_sensor_volume* sensor_volume;
+ struct sphin* sphin;
+ ref_T ref;
+};
+
extern LOCAL_SYM res_T
parse_volume
(struct sphin_config* config,
diff --git a/src/test_sphin_load_prop_rad.c b/src/test_sphin_load_prop_rad.c
@@ -0,0 +1,318 @@
+/* Copyright (C) 2024-2025 Centre National de la Recherche Scientifique
+ * Copyright (C) 2024-2025 Clermont Auvergne INP
+ * Copyright (C) 2024-2025 INSA Lyon
+ * Copyright (C) 2024-2025 Institut Mines Télécom Albi-Carmaux
+ * Copyright (C) 2024-2025 Institut National Polytechnique de Toulouse
+ * Copyright (C) 2024-2025 |Méso|Star> (contact@meso-star.com)
+ * Copyright (C) 2024-2025 PhotonLyX (info@photonlyx.com)
+ * Copyright (C) 2024-2025 Université de Lorraine
+ * Copyright (C) 2024-2025 Université Paul Sabatier
+ * Copyright (C) 2024-2025 Université Toulouse - Jean Jaurès
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sphin.h"
+
+#include <rsys/math.h>
+#include <rsys/mem_allocator.h>
+#include <string.h>
+
+
+static void
+write_prop_rad_test_files
+ (void)
+{
+ FILE* file;
+ static const char* test0 =
+ "#wv abs_cross_sec\n"
+ "100 1\n"
+ "200 2\n"
+ "300 3\n"
+ "400 4\n";
+ static const char* test1 = /* Invalid first line */
+ "wv abs_cross_sec\n"
+ "100 1\n"
+ "200 2\n"
+ "300 3\n"
+ "400 4\n";
+ static const char* test2 = /* Unsorted wavelengths */
+ "100 1\n"
+ "300 3\n"
+ "200 2\n"
+ "400 4\n";
+ static const char* test3 = /* Missing data */
+ "100 1\n"
+ "200 2\n"
+ "300 \n"
+ "400 4\n";
+ static const char* test4 = /* Bad formating */
+ "100 1\n"
+ "200 2\n"
+ "300 2,5\n"
+ "400 4\n";
+ file = fopen("abs_cross_sec_0.txt", "w");
+ CHK(file != NULL);
+ fwrite(test0, sizeof(char), strlen(test0), file);
+ fclose(file);
+
+ file = fopen("abs_cross_sec_1.txt", "w");
+ CHK(file != NULL);
+ fwrite(test1, sizeof(char), strlen(test1), file);
+ fclose(file);
+
+ file = fopen("abs_cross_sec_2.txt", "w");
+ CHK(file != NULL);
+ fwrite(test2, sizeof(char), strlen(test2), file);
+ fclose(file);
+
+ file = fopen("abs_cross_sec_3.txt", "w");
+ CHK(file != NULL);
+ fwrite(test3, sizeof(char), strlen(test3), file);
+ fclose(file);
+
+ file = fopen("abs_cross_sec_4.txt", "w");
+ CHK(file != NULL);
+ fwrite(test4, sizeof(char), strlen(test4), file);
+ fclose(file);
+}
+
+static void
+test_prop_rad_api
+ (struct sphin* sphin)
+{
+ struct sphin_config* config = NULL;
+ const char* path = "file.txt";
+ struct sphin_volume* volume = NULL;
+ struct sphin_prop_rad* prop_rad = NULL;
+ struct sphin_scatterer* scatterer = NULL;
+ struct sphin_spectral_property* abs_cross_sec = NULL;
+ struct sphin_spectral_property_descriptor abs_cs_desc =
+ SPHIN_SPECTRAL_PROPERTY_DESCRIPTOR_NULL;
+ size_t nprop_rad;
+ double concentration, value;
+ FILE* fp = NULL;
+ CHK(fp = fopen(path, "w+"));
+ fprintf(fp, "#Mot Clé Nom\n");
+ fprintf(fp, "volume: \"reaction volume\"\n");
+ fprintf(fp, "\tprop_rad: \"name\" SCATTERER\n");
+ fprintf(fp, "\t\tconcentration: 1.935 mol/m^3\n");
+ fprintf(fp, "\t\tcross_sections:\n");
+ fprintf(fp, "\t\t\tabs_cross_sec: abs_cross_sec_0.txt nm m^2/part\n");
+ fclose(fp);
+
+ CHK(sphin_load(sphin, path, &config) == RES_OK);
+ CHK(sphin_config_get_volume(config, 0, &volume) == RES_OK);
+ CHK(sphin_volume_get_prop_rad_count(NULL, NULL) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad_count(volume, NULL) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad_count(NULL, &nprop_rad) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad_count(volume, &nprop_rad) == RES_OK);
+ CHK(nprop_rad == 1);
+ CHK(sphin_volume_get_prop_rad(NULL, 0, NULL) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad(volume, 0, NULL) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad(NULL, 0, &prop_rad) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad(volume, (size_t)-1, &prop_rad) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad(volume, 10, &prop_rad) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad(volume, 0, &prop_rad) == RES_OK);
+ CHK(sphin_prop_rad_ref_get(prop_rad) == RES_OK);
+ CHK(sphin_prop_rad_ref_get(NULL) == RES_BAD_ARG);
+ CHK(sphin_prop_rad_ref_put(NULL) == RES_BAD_ARG);
+ CHK(sphin_prop_rad_ref_put(prop_rad) == RES_OK);
+ CHK(sphin_prop_rad_get_scatterer(NULL, NULL) == RES_BAD_ARG);
+ CHK(sphin_prop_rad_get_scatterer(prop_rad, NULL) == RES_BAD_ARG);
+ CHK(sphin_prop_rad_get_scatterer(NULL, &scatterer) == RES_BAD_ARG);
+ CHK(sphin_prop_rad_get_scatterer(prop_rad, &scatterer) == RES_OK);
+ CHK(sphin_scatterer_ref_get(NULL) == RES_BAD_ARG);
+ CHK(sphin_scatterer_ref_get(scatterer) == RES_OK);
+ CHK(sphin_scatterer_ref_put(NULL) == RES_BAD_ARG);
+ CHK(sphin_scatterer_ref_put(scatterer) == RES_OK);
+ CHK(sphin_scatterer_get_concentration(NULL, NULL) == RES_BAD_ARG);
+ CHK(sphin_scatterer_get_concentration(scatterer, NULL) == RES_BAD_ARG);
+ CHK(sphin_scatterer_get_concentration(NULL, &concentration) == RES_BAD_ARG);
+ CHK(sphin_scatterer_get_concentration(scatterer, &concentration) == RES_OK);
+ CHK(eq_eps(concentration, 1.935, 1e-15));
+ CHK(sphin_scatterer_get_abs_cross_sec(NULL, NULL) == RES_BAD_ARG);
+ CHK(sphin_scatterer_get_abs_cross_sec(scatterer, NULL) == RES_BAD_ARG);
+ CHK(sphin_scatterer_get_abs_cross_sec(NULL, &abs_cross_sec) == RES_BAD_ARG);
+ CHK(sphin_scatterer_get_abs_cross_sec(scatterer, &abs_cross_sec) == RES_OK);
+ CHK(sphin_spectral_property_ref_get(NULL) == RES_BAD_ARG);
+ CHK(sphin_spectral_property_ref_get(abs_cross_sec) == RES_OK);
+ CHK(sphin_spectral_property_ref_put(NULL) == RES_BAD_ARG);
+ CHK(sphin_spectral_property_ref_put(abs_cross_sec) == RES_OK);
+ CHK(sphin_spectral_property_get_desc(NULL, NULL) == RES_BAD_ARG);
+ CHK(sphin_spectral_property_get_desc(abs_cross_sec, NULL) == RES_BAD_ARG);
+ CHK(sphin_spectral_property_get_desc(NULL, &abs_cs_desc) == RES_BAD_ARG);
+ CHK(sphin_spectral_property_get_desc(abs_cross_sec, &abs_cs_desc) == RES_OK);
+ CHK(sphin_spectral_property_interpolate_at_wavelength
+ (abs_cross_sec, 100, &value) == RES_OK);
+ CHK(eq_eps(value, 1., 1e-15));
+ CHK(sphin_spectral_property_interpolate_at_wavelength
+ (abs_cross_sec, 175, &value) == RES_OK);
+ CHK(eq_eps(value, 1.75, 1e-15));
+ CHK(sphin_spectral_property_interpolate_at_wavelength
+ (abs_cross_sec, 400, &value) == RES_OK);
+ CHK(eq_eps(value, 4., 1e-15));
+ CHK(eq_eps(abs_cs_desc.wavelengths[0], 100, 1e-15));
+ CHK(eq_eps(abs_cs_desc.wavelengths[2], 300, 1e-15));
+ CHK(eq_eps(abs_cs_desc.values[0], 1, 1e-15));
+ CHK(eq_eps(abs_cs_desc.values[2], 3, 1e-15));
+ CHK(abs_cs_desc.data_count == 4);
+ CHK(sphin_config_ref_put(config) == RES_OK);
+}
+
+static void
+test_prop_rad_api_units
+ (struct sphin* sphin)
+{
+ struct sphin_config* config = NULL;
+ const char* path = "file.txt";
+ struct sphin_volume* volume = NULL;
+ struct sphin_prop_rad* prop_rad = NULL;
+ struct sphin_scatterer* scatterer = NULL;
+ struct sphin_spectral_property* abs_cross_sec = NULL;
+ struct sphin_spectral_property_descriptor abs_cs_desc =
+ SPHIN_SPECTRAL_PROPERTY_DESCRIPTOR_NULL;
+ double value;
+ FILE* fp = NULL;
+ CHK(fp = fopen(path, "w+"));
+ fprintf(fp, "#Mot Clé Nom\n");
+ fprintf(fp, "volume: \"reaction volume\"\n");
+ fprintf(fp, "\tprop_rad: \"name\" SCATTERER\n");
+ fprintf(fp, "\t\tconcentration: 1.935 mol/m^3\n");
+ fprintf(fp, "\t\tcross_sections:\n");
+ fprintf(fp, "\t\t\tabs_cross_sec: abs_cross_sec_0.txt cm m^2/part\n");
+ fclose(fp);
+
+ CHK(sphin_load(sphin, path, &config) == RES_OK);
+ CHK(sphin_config_get_volume(config, 0, &volume) == RES_OK);
+ CHK(sphin_volume_get_prop_rad(volume, 0, &prop_rad) == RES_OK);
+ CHK(sphin_prop_rad_get_scatterer(prop_rad, &scatterer) == RES_OK);
+ CHK(sphin_scatterer_get_abs_cross_sec(scatterer, &abs_cross_sec) == RES_OK);
+ CHK(sphin_spectral_property_get_desc(abs_cross_sec, &abs_cs_desc) == RES_OK);
+ CHK(sphin_spectral_property_interpolate_at_wavelength
+ (abs_cross_sec, 100e7, &value) == RES_OK);
+ CHK(eq_eps(value, 1., 1e-15));
+ CHK(sphin_spectral_property_interpolate_at_wavelength
+ (abs_cross_sec, 175e7, &value) == RES_OK);
+ CHK(eq_eps(value, 1.75, 1e-15));
+ CHK(sphin_spectral_property_interpolate_at_wavelength
+ (abs_cross_sec, 400e7, &value) == RES_OK);
+ CHK(eq_eps(value, 4., 1e-15));
+ CHK(eq_eps(abs_cs_desc.wavelengths[0], 100e7, 1e-15));
+ CHK(eq_eps(abs_cs_desc.wavelengths[2], 300e7, 1e-15));
+ CHK(sphin_config_ref_put(config) == RES_OK);
+}
+static void
+test_prop_rad_api_bad_comment
+ (struct sphin* sphin)
+{
+ struct sphin_config* config = NULL;
+ const char* path = "file.txt";
+ FILE* fp = NULL;
+
+ CHK(fp = fopen(path, "w+"));
+ fprintf(fp, "#Mot Clé Nom\n");
+ fprintf(fp, "volume: \"reaction volume\"\n");
+ fprintf(fp, "\tprop_rad: SCATTERER \"name\"\n");
+ fprintf(fp, "\t\tconcentration: 1.935 mol/m^3\n");
+ fprintf(fp, "\t\tcross_sections:\n");
+ fprintf(fp, "\t\t\tabs_cross_sec: abs_cross_sec_1.txt nm m^2/part\n");
+ fclose(fp);
+
+ CHK(sphin_load(sphin, path, &config) == RES_BAD_ARG);
+}
+
+static void
+test_prop_rad_api_unsorted_wavelengths
+ (struct sphin* sphin)
+{
+ struct sphin_config* config = NULL;
+ const char* path = "file.txt";
+ FILE* fp = NULL;
+
+ CHK(fp = fopen(path, "w+"));
+ fprintf(fp, "#Mot Clé Nom\n");
+ fprintf(fp, "volume: \"reaction volume\"\n");
+ fprintf(fp, "\tprop_rad: SCATTERER \"name\"\n");
+ fprintf(fp, "\t\tconcentration: 1.935 mol/m^3\n");
+ fprintf(fp, "\t\tcross_sections:\n");
+ fprintf(fp, "\t\t\tabs_cross_sec: abs_cross_sec_2.txt nm m^2/part\n");
+ fclose(fp);
+
+ CHK(sphin_load(sphin, path, &config) == RES_BAD_ARG);
+}
+
+static void
+test_prop_rad_api_missing_data
+ (struct sphin* sphin)
+{
+ struct sphin_config* config = NULL;
+ const char* path = "file.txt";
+ FILE* fp = NULL;
+
+ CHK(fp = fopen(path, "w+"));
+ fprintf(fp, "#Mot Clé Nom\n");
+ fprintf(fp, "volume: \"reaction volume\"\n");
+ fprintf(fp, "\tprop_rad: SCATTERER \"name\"\n");
+ fprintf(fp, "\t\tconcentration: 1.935 mol/m^3\n");
+ fprintf(fp, "\t\tcross_sections:\n");
+ fprintf(fp, "\t\t\tabs_cross_sec: abs_cross_sec_3.txt nm m^2/part\n");
+ fclose(fp);
+
+ CHK(sphin_load(sphin, path, &config) == RES_BAD_ARG);
+}
+
+static void
+test_prop_rad_api_bad_formatting
+ (struct sphin* sphin)
+{
+ struct sphin_config* config = NULL;
+ const char* path = "file.txt";
+ FILE* fp = NULL;
+
+ CHK(fp = fopen(path, "w+"));
+ fprintf(fp, "#Mot Clé Nom\n");
+ fprintf(fp, "volume: \"reaction volume\"\n");
+ fprintf(fp, "\tprop_rad: SCATTERER \"name\"\n");
+ fprintf(fp, "\t\tconcentration: 1.935 mol/m^3\n");
+ fprintf(fp, "\t\tcross_sections:\n");
+ fprintf(fp, "\t\t\tabs_cross_sec: abs_cross_sec_4.txt nm m^2/part\n");
+ fclose(fp);
+
+ CHK(sphin_load(sphin, path, &config) == RES_BAD_ARG);
+}
+
+int
+main(int argc, char** argv)
+{
+ struct sphin* sphin = NULL;
+ struct sphin_create_args args = SPHIN_CREATE_ARGS_DEFAULT;
+
+ (void)argc;
+ (void)argv;
+
+ args.verbose = 1;
+ sphin_create(&args, &sphin);
+
+ write_prop_rad_test_files();
+ test_prop_rad_api(sphin);
+ test_prop_rad_api_units(sphin);
+ test_prop_rad_api_bad_comment(sphin);
+ test_prop_rad_api_unsorted_wavelengths(sphin);
+ test_prop_rad_api_missing_data(sphin);
+ test_prop_rad_api_bad_formatting(sphin);
+ CHK(sphin_ref_put(sphin) == RES_OK);
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sphin_load_source.c b/src/test_sphin_load_source.c
@@ -50,32 +50,54 @@ write_stl_test_files
}
static void
+write_emission_spectrum_test_files
+ (void)
+{
+ FILE* file;
+ static const char* test0 =
+ "#wv,abs_cross_sec\n"
+ "100 1\n"
+ "200 2\n"
+ "300 3\n"
+ "400 4\n";
+ file = fopen("source_spec.txt", "w");
+ CHK(file != NULL);
+ fwrite(test0, sizeof(char), strlen(test0), file);
+ fclose(file);
+}
+
+static void
test_source_api
(struct sphin* sphin)
{
const char* path = "file.txt";
double power;
+ double value;
size_t nsurfaces, ngeometries;
struct sphin_config* config = NULL;
struct sphin_surface* surface = NULL;
struct sphin_source_surface* source = NULL;
- struct sphin_source_surface_direction_distribution dir_dist = SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_NULL;
- struct sphin_source_surface_flux_density flux_density = SPHIN_SOURCE_SURFACE_FLUX_DENSITY_NULL;
+ struct sphin_source_surface_direction_distribution dir_dist =
+ SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_NULL;
+ struct sphin_source_surface_flux_density flux_density =
+ SPHIN_SOURCE_SURFACE_FLUX_DENSITY_NULL;
+ struct sphin_spectral_property_descriptor emission_desc =
+ SPHIN_SPECTRAL_PROPERTY_DESCRIPTOR_NULL;
FILE* fp = NULL;
CHK(fp = fopen(path, "w+"));
- fprintf(fp, "\t\t surface : \"surface name\"\n");
+ fprintf(fp, "\t\t surface : \"surface name\"\n"); /* Source 0 */
fprintf(fp, "\tgeometry: FRONT test_0.stl\t\n");
fprintf(fp, "\tsource:#commentaire\n");
- fprintf(fp, "\tflux_density: 200e-6 mol/m^2/s\t\n");
+ fprintf(fp, "\tflux_density: 200e-6 mol/m^2/s \t\n");
fprintf(fp, "\tdirection: LAMBERT\n");
- fprintf(fp, "\t\t surface : \"surface name 2\"\n");
+ fprintf(fp, "\t\t surface : \"surface name 2\"\n"); /* Source 1 */
fprintf(fp, "\tgeometry: BACK test_0.stl\t\n");
fprintf(fp, "\tsource:\t #commentaire\n");
- fprintf(fp, "\tflux_density: 1e3 umol/m^2/s\t\n");
+ fprintf(fp, "\tflux_density: 1e3 umol/m^2/s source_spec.txt nm unit\t\n");
fprintf(fp, "\tdirection: COLLIM NORMAL\n");
- fprintf(fp, "\t\t surface : \"surface name 2\"\n");
+ fprintf(fp, "\t\t surface : \"surface name 3\"\n"); /* Source 2 */
fprintf(fp, "\tgeometry: BACK test_0.stl\t\n");
fprintf(fp, "\tsource:\t #commentaire\n");
fprintf(fp, "\tflux_density: 200e3 mW/m^2\t\n");
@@ -94,10 +116,14 @@ test_source_api
CHK(sphin_source_surface_get_direction_distribution(source, NULL) == RES_BAD_ARG);
CHK(sphin_source_surface_get_direction_distribution(source, &dir_dist) == RES_OK);
CHK(dir_dist.type == SPHIN_SOURCE_DIRECTION_ISOTROPIC);
- CHK(dir_dist.cos_pow_n.collimation_degree == SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COS_POW_N_NULL.collimation_degree);
- CHK(dir_dist.collim.direction[0] == SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COLLIM_NULL.direction[0]);
- CHK(dir_dist.collim.direction[1] == SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COLLIM_NULL.direction[1]);
- CHK(dir_dist.collim.direction[2] == SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COLLIM_NULL.direction[2]);
+ CHK(dir_dist.cos_pow_n.collimation_degree ==
+ SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COS_POW_N_NULL.collimation_degree);
+ CHK(dir_dist.collim.direction[0] ==
+ SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COLLIM_NULL.direction[0]);
+ CHK(dir_dist.collim.direction[1] ==
+ SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COLLIM_NULL.direction[1]);
+ CHK(dir_dist.collim.direction[2] ==
+ SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COLLIM_NULL.direction[2]);
CHK(sphin_source_surface_get_flux_density(NULL, &flux_density) == RES_BAD_ARG);
CHK(sphin_source_surface_get_flux_density(source, NULL) == RES_BAD_ARG);
CHK(sphin_source_surface_get_flux_density(source, &flux_density) == RES_OK);
@@ -112,22 +138,38 @@ test_source_api
CHK(sphin_surface_get_source(surface, &source) == RES_OK);
CHK(sphin_source_surface_get_direction_distribution(source, &dir_dist) == RES_OK);
CHK(dir_dist.type == SPHIN_SOURCE_DIRECTION_COLLIM);
- CHK(dir_dist.cos_pow_n.collimation_degree == SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COS_POW_N_NULL.collimation_degree);
+ CHK(dir_dist.cos_pow_n.collimation_degree ==
+ SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COS_POW_N_NULL.collimation_degree);
CHK(dir_dist.collim.direction[0] == 0);
CHK(dir_dist.collim.direction[1] == 0);
CHK(dir_dist.collim.direction[2] == 0);
CHK(sphin_source_surface_get_flux_density(source, &flux_density) == RES_OK);
CHK(eq_eps(flux_density.flux_density, 1000, 1e-15));
CHK(flux_density.unit == SPHIN_PHOTON_UNIT_MOL);
+ CHK(sphin_spectral_property_ref_get(flux_density.emission_spectrum) == RES_OK);
+ CHK(sphin_spectral_property_ref_put(flux_density.emission_spectrum) == RES_OK);
+ CHK(sphin_spectral_property_get_desc(flux_density.emission_spectrum, &emission_desc) == RES_OK);
+ CHK(sphin_spectral_property_interpolate_at_wavelength
+ (flux_density.emission_spectrum, 100, &value) == RES_OK);
+ CHK(eq_eps(value, 1., 1e-15));
+ CHK(sphin_spectral_property_interpolate_at_wavelength
+ (flux_density.emission_spectrum, 175, &value) == RES_OK);
+ CHK(eq_eps(value, 1.75, 1e-15));
+ CHK(eq_eps(emission_desc.wavelengths[2], 300, 1e-15));
+ CHK(eq_eps(emission_desc.values[2], 3, 1e-15));
+ CHK(emission_desc.data_count == 4);
CHK(sphin_config_get_surface(config, 2, &surface) == RES_OK);
CHK(sphin_surface_get_source(surface, &source) == RES_OK);
CHK(sphin_source_surface_get_direction_distribution(source, &dir_dist) == RES_OK);
CHK(dir_dist.type == SPHIN_SOURCE_DIRECTION_COS_POW_N);
CHK(eq_eps(dir_dist.cos_pow_n.collimation_degree, 1, 1e-15));
- CHK(dir_dist.collim.direction[0] == SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COLLIM_NULL.direction[0]);
- CHK(dir_dist.collim.direction[1] == SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COLLIM_NULL.direction[1]);
- CHK(dir_dist.collim.direction[2] == SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COLLIM_NULL.direction[2]);
+ CHK(dir_dist.collim.direction[0] ==
+ SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COLLIM_NULL.direction[0]);
+ CHK(dir_dist.collim.direction[1] ==
+ SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COLLIM_NULL.direction[1]);
+ CHK(dir_dist.collim.direction[2] ==
+ SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_COLLIM_NULL.direction[2]);
CHK(sphin_source_surface_get_flux_density(source, NULL) == RES_BAD_ARG);
CHK(sphin_source_surface_get_flux_density(NULL, &flux_density) == RES_BAD_ARG);
CHK(sphin_source_surface_get_flux_density(source, &flux_density) == RES_OK);
@@ -249,6 +291,7 @@ main(int argc, char** argv)
sphin_create(&args, &sphin);
write_stl_test_files();
+ write_emission_spectrum_test_files();
test_source_api(sphin);
test_source_api_bad_flux_density(sphin);
test_source_api_bad_flux_density_value(sphin);
diff --git a/src/test_sphin_load_volume.c b/src/test_sphin_load_volume.c
@@ -86,6 +86,23 @@ write_stl_test_files
}
static void
+write_prop_rad_test_files
+ (void)
+{
+ FILE* file;
+ static const char* test0 =
+ "#wv,abs_cross_sec\n"
+ "100 1\n"
+ "200 2\n"
+ "300 3\n"
+ "400 4\n";
+ file = fopen("abs_cross_sec_0.txt", "w");
+ CHK(file != NULL);
+ fwrite(test0, sizeof(char), strlen(test0), file);
+ fclose(file);
+}
+
+static void
test_volume_api
(struct sphin* sphin)
{
@@ -93,13 +110,19 @@ test_volume_api
const char* path = "file.txt";
struct sphin_volume* volume = NULL;
struct sphin_sensor_volume* sensor_volume = NULL;
+ struct sphin_prop_rad* prop_rad = NULL;
size_t nvolumes;
+ size_t nprop_rad;
double ka, response_function;
FILE* fp = NULL;
CHK(fp = fopen(path, "w+"));
fprintf(fp, "#Mot Clé Nom\n");
fprintf(fp, "volume: \"reaction volume\"\n");
fprintf(fp, "\tka: 2.11 m^-1\n");
+ fprintf(fp, "\tprop_rad: \"name\" SCATTERER\n");
+ fprintf(fp, "\t\tconcentration: 1.935 mol/m^3\n");
+ fprintf(fp, "\t\tcross_sections:\n");
+ fprintf(fp, "\t\t\tabs_cross_sec: abs_cross_sec_0.txt nm m^2/part\n");
fprintf(fp, "\tsensor:\n");
fprintf(fp, "\t\tresponse_function: 1\n");
fprintf(fp, "\tgeometry: FRONT test_0.stl\n");
@@ -131,6 +154,17 @@ test_volume_api
CHK(sphin_volume_get_ka(volume, NULL) == RES_BAD_ARG);
CHK(sphin_volume_get_ka(volume, &ka) == RES_OK);
CHK(eq_eps(ka, 2.11, 1e-15));
+ CHK(sphin_volume_get_prop_rad_count(NULL, NULL) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad_count(volume, NULL) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad_count(NULL, &nprop_rad) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad_count(volume, &nprop_rad) == RES_OK);
+ CHK(nprop_rad == 1);
+ CHK(sphin_volume_get_prop_rad(NULL, 0, NULL) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad(volume, 0, NULL) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad(NULL, 0, &prop_rad) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad(volume, (size_t)-1, &prop_rad) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad(volume, 10, &prop_rad) == RES_BAD_ARG);
+ CHK(sphin_volume_get_prop_rad(volume, 0, &prop_rad) == RES_OK);
CHK(sphin_volume_get_sensor(NULL, NULL) == RES_BAD_ARG);
CHK(sphin_volume_get_sensor(NULL, &sensor_volume) == RES_BAD_ARG);
CHK(sphin_volume_get_sensor(volume, NULL) == RES_BAD_ARG);
@@ -228,6 +262,7 @@ main(int argc, char** argv)
sphin_create(&args, &sphin);
write_stl_test_files();
+ write_prop_rad_test_files();
test_volume_api(sphin);
test_volume_api_empty_file(sphin);
test_volume_api_bad_geometry(sphin);
