commit e88aeea9206c33e54cd184a8876046af19fd8ced
parent facc7d874dcb11552c33cc175ac7469b8ff2190b
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 23 Sep 2026 12:32:56 +0200
Merge branch 'release_0.1'
Diffstat:
32 files changed, 1513 insertions(+), 125 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -10,6 +10,7 @@ mixture.txt
sln-build
sln-get
sln-slab
+sln-stat
tags
tags
test_*
diff --git a/Makefile b/Makefile
@@ -3,6 +3,8 @@
# Copyright (C) 2022 Centre National de la Recherche Scientifique
# Copyright (C) 2022 Université Paul Sabatier
#
+# This file is part of Star-Line.
+#
# 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
@@ -84,7 +86,7 @@ libsln.o: $(OBJ)
################################################################################
# Utils
################################################################################
-UTIL_SRC = src/sln_build.c src/sln_get.c src/sln_slab.c
+UTIL_SRC = src/sln_build.c src/sln_get.c src/sln_slab.c src/sln_stat.c
UTIL_OBJ = $(UTIL_SRC:.c=.o)
UTIL_DEP = $(UTIL_SRC:.c=.d)
@@ -103,14 +105,14 @@ LDFLAGS_SLAB = $(LDFLAGS_EXE) $(LIBS_SLAB)
utils: library $(UTIL_DEP) .config
@$(MAKE) -fMakefile \
$$(for i in $(UTIL_DEP); do printf -- '-f%s\n' "$${i}"; done) \
- sln-build sln-get sln-slab
+ sln-build sln-get sln-slab sln-stat
$(UTIL_DEP) $(UTIL_OBJ): config.mk sln-local.pc
-
src/sln_build.d src/sln_build.o: src/sln_build.c
src/sln_get.d src/sln_get.o: src/sln_get.c
src/sln_slab.d src/sln_slab.o: src/sln_slab.c
+src/sln_stat.d src/sln_stat.o: src/sln_stat.c
sln-build: config.mk sln-local.pc src/sln_build.o $(LIBNAME)
$(CC) $(CFLAGS_UTIL) -o $@ src/sln_build.o $(LDFLAGS_UTIL)
@@ -121,20 +123,23 @@ sln-get: config.mk sln-local.pc src/sln_get.o $(LIBNAME)
sln-slab: config.mk sln-local.pc src/sln_slab.o $(LIBNAME)
$(CC) $(CFLAGS_SLAB) -o $@ src/sln_slab.o $(LDFLAGS_SLAB)
+sln-stat: config.mk sln-local.pc src/sln_stat.o $(LIBNAME)
+ $(CC) $(CFLAGS_SLAB) -o $@ src/sln_stat.o $(LDFLAGS_SLAB)
+
src/sln_build.d src/sln_get.d:
@$(CC) $(CFLAGS_UTIL) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
-src/sln_slab.d:
+src/sln_slab.d src/sln_stat.d:
@$(CC) $(CFLAGS_SLAB) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
src/sln_build.o src/sln_get.o:
$(CC) $(CFLAGS_UTIL) -c $(@:.o=.c) -o $@
-src/sln_slab.o:
+src/sln_slab.o src/sln_stat.o:
$(CC) $(CFLAGS_SLAB) -c $(@:.o=.c) -o $@
clean_utils:
- rm -f $(UTIL_OBJ) $(UTIL_DEP) sln-build sln-get sln-slab
+ rm -f $(UTIL_OBJ) $(UTIL_DEP) sln-build sln-get sln-slab sln-stat
################################################################################
# Installation
@@ -172,11 +177,13 @@ install: library pkg utils
install 755 "$(DESTDIR)$(BINPREFIX)" sln-build; \
install 755 "$(DESTDIR)$(BINPREFIX)" sln-get; \
install 755 "$(DESTDIR)$(BINPREFIX)" sln-slab; \
+ install 755 "$(DESTDIR)$(BINPREFIX)" sln-stat; \
install 644 "$(DESTDIR)$(LIBPREFIX)/pkgconfig" sln.pc; \
install 644 "$(DESTDIR)$(INCPREFIX)/star" src/sln.h; \
install 644 "$(DESTDIR)$(MANPREFIX)/man1" doc/sln-build.1; \
install 644 "$(DESTDIR)$(MANPREFIX)/man1" doc/sln-get.1; \
install 644 "$(DESTDIR)$(MANPREFIX)/man1" doc/sln-slab.1; \
+ install 644 "$(DESTDIR)$(MANPREFIX)/man1" doc/sln-stat.1; \
install 644 "$(DESTDIR)$(MANPREFIX)/man5" doc/sln-mixture.5; \
install 644 "$(DESTDIR)$(PREFIX)/share/doc/star-line" COPYING README.md
@@ -185,12 +192,14 @@ uninstall:
rm -f "$(DESTDIR)$(BINPREFIX)/sln-build"
rm -f "$(DESTDIR)$(BINPREFIX)/sln-get"
rm -f "$(DESTDIR)$(BINPREFIX)/sln-slab"
+ rm -f "$(DESTDIR)$(BINPREFIX)/sln-stat"
rm -f "$(DESTDIR)$(LIBPREFIX)/pkgconfig/sln.pc"
rm -f "$(DESTDIR)$(BINPREFIX)/sln"
rm -f "$(DESTDIR)$(INCPREFIX)/star/sln.h"
rm -f "$(DESTDIR)$(MANPREFIX)/man1/sln-build.1"
rm -f "$(DESTDIR)$(MANPREFIX)/man1/sln-get.1"
rm -f "$(DESTDIR)$(MANPREFIX)/man1/sln-slab.1"
+ rm -f "$(DESTDIR)$(MANPREFIX)/man1/sln-stat.1"
rm -f "$(DESTDIR)$(MANPREFIX)/man5/sln-mixture.5"
rm -f "$(DESTDIR)$(PREFIX)/share/doc/star-line/COPYING"
rm -f "$(DESTDIR)$(PREFIX)/share/doc/star-line/README.md"
@@ -203,6 +212,7 @@ lint:
mandoc -Tlint -Wwarning doc/sln-build.1
mandoc -Tlint -Wwarning doc/sln-get.1
mandoc -Tlint -Wwarning doc/sln-slab.1
+ mandoc -Tlint -Wwarning doc/sln-stat.1
mandoc -Tlint doc/sln-mixture.5
################################################################################
@@ -212,6 +222,7 @@ TEST_SRC =\
src/test_sln_device.c\
src/test_sln_mesh.c\
src/test_sln_mixture.c\
+ src/test_sln_thermo_props.c\
src/test_sln_tree.c\
src/test_sln_tree_sample.c
TEST_OBJ = $(TEST_SRC:.c=.o)
@@ -246,6 +257,7 @@ $(TEST_OBJ): config.mk sln-local.pc
test_sln_device \
test_sln_mesh \
test_sln_mixture \
+test_sln_thermo_props\
test_sln_tree_sample \
: config.mk sln-local.pc $(LIBNAME)
$(CC) $(CFLAGS_TEST) -o $@ src/$@.o $(LDFLAGS_TEST)
diff --git a/README.md b/README.md
@@ -21,6 +21,19 @@ Edit config.mk as needed, then run:
## Release notes
+### Version 0.1
+
+- Allow the caller to modify the thermodynamic properties at which it
+ evaluates the node of a tree or the lines it partitions, i.e., to use
+ others than those used to build the tree.
+ The same acceleration structure can therefore be used to sample lines
+ with different thermodynamic properties *as long as the sampling
+ criteria are met*;
+ such as the polylines of nodes that must be an upper limit of the
+ lines they encompass.
+- Provide the `sln-stat` utility, which uses line importance sampling to
+ estimate the absorption coefficient of a spectrum, and its square.
+
### Version 0.0
- Initial version of the library that structures a set of lines into an
diff --git a/config.mk b/config.mk
@@ -1,4 +1,4 @@
-VERSION = 0.0
+VERSION = 0.1
PREFIX = /usr/local
BINPREFIX = $(PREFIX)/bin
diff --git a/doc/sln-build.1 b/doc/sln-build.1
@@ -3,6 +3,8 @@
.\" Copyright (C) 2022 Centre National de la Recherche Scientifique
.\" Copyright (C) 2022 Université Paul Sabatier
.\"
+.\" This file is part of Star-Line.
+.\"
.\" 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
@@ -15,7 +17,7 @@
.\"
.\" You should have received a copy of the GNU General Public License
.\" along with this program. If not, see <http://www.gnu.org/licenses/>.
-.Dd April 30, 2026
+.Dd August 17, 2026
.Dt SLN-BUILD 1
.Os
.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
@@ -204,7 +206,7 @@ global radiative cooling
.Rs
.%A L.S. Rothman et al.
.%T HITEMP, the high-temperature molecular spectroscopic database
-.%J Journal of Quantitative Spectroscopu & Radiative Transfer
+.%J Journal of Quantitative Spectroscopy & Radiative Transfer
.%V 111
.%P pp. 2139\(en2150
.%D 2010
diff --git a/doc/sln-get.1 b/doc/sln-get.1
@@ -3,6 +3,8 @@
.\" Copyright (C) 2022 Centre National de la Recherche Scientifique
.\" Copyright (C) 2022 Université Paul Sabatier
.\"
+.\" This file is part of Star-Line.
+.\"
.\" 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
@@ -15,7 +17,7 @@
.\"
.\" You should have received a copy of the GNU General Public License
.\" along with this program. If not, see <http://www.gnu.org/licenses/>.
-.Dd April 10, 2026
+.Dd August 17, 2026
.Dt SLN-GET 1
.Os
.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
@@ -229,7 +231,7 @@ sln-get -i lines.par -p molparams.txt -c0:3 -w 50 tree.sln
.Rs
.%A L.S. Rothman et al.
.%T HITEMP, the high-temperature molecular spectroscopic database
-.%J Journal of Quantitative Spectroscopu & Radiative Transfer
+.%J Journal of Quantitative Spectroscopy & Radiative Transfer
.%V 111
.%P pp. 2139\(en2150
.%D 2010
diff --git a/doc/sln-mixture.5 b/doc/sln-mixture.5
@@ -3,6 +3,8 @@
.\" Copyright (C) 2022 Centre National de la Recherche Scientifique
.\" Copyright (C) 2022 Université Paul Sabatier
.\"
+.\" This file is part of Star-Line.
+.\"
.\" 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
diff --git a/doc/sln-slab.1 b/doc/sln-slab.1
@@ -3,6 +3,8 @@
.\" Copyright (C) 2022 Centre National de la Recherche Scientifique
.\" Copyright (C) 2022 Université Paul Sabatier
.\"
+.\" This file is part of Star-Line.
+.\"
.\" 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
diff --git a/doc/sln-stat.1 b/doc/sln-stat.1
@@ -0,0 +1,170 @@
+.\" Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com)
+.\" Copyright (C) 2026 Université de Lorraine
+.\" Copyright (C) 2022 Centre National de la Recherche Scientifique
+.\" Copyright (C) 2022 Université Paul Sabatier
+.\"
+.\" This file is part of Star-Line.
+.\"
+.\" 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/>.
+.Dd May 5, 2026
+.Dt SLN-STAT 1
+.Os
+.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+.Sh NAME
+.Nm sln-stat
+.Nd computations of basic statistics over k spectrum
+.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+.Sh SYNOPSIS
+.Nm
+.Op Fl dhsv
+.Op Fl n Ar nrealisations
+.Op Fl t Ar threads
+.Fl S Ar nu_min , Ns Ar nu_max
+.Fl a Ar accel_struct
+.Fl m Ar molparams
+.Fl l Ar lines
+.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+.Sh DESCRIPTION
+.Nm
+calculates the mean and mean of squares of k absorption values using a Monte
+Carlo algorithm that samples the spectral lines that make up the gas mixture.
+These computations are accelerated by sampling the lines based on the magnitude
+of their contribution to the mixture’s spectrum, so that few Monte Carlo runs
+are required to estimate the spectrum statistics with a high degree of
+confidence.
+The core of the proposal rests on this sampling strategy, made possible
+by constructing an acceleration structure from the set of lines in the
+mixture.
+A structure built using the
+.Xr sln-build 1
+utility and provided as input to the program.
+.Pp
+The output of
+.Nm
+displays the estimated mean and mean of squares, their standard
+deviation, and the number of Monte Carlo realisations rejected due to
+issues encountered during the computation, such as numerical
+uncertainty.
+Each estimate is displayed on a line formatted as follows:
+.Bd -literal -offset Ds
+"%-16s: %e +/- %e; %lu\en", name, estimate, std_err, rejects_count
+.Ed
+.Pp
+The options are as follows:
+.Bl -tag -width Ds
+.\""""""""""""""""""""""""""""""""""
+.It Fl a Ar accel_struct
+An acceleration structure corresponding to the input
+.Ar lines ,
+used to accelerate their sampling based on their importance.
+This structure is generated by the
+.Xr sln-build 1
+tool.
+.\""""""""""""""""""""""""""""""""""
+.It Fl d
+Disables verification of the correspondence between the lines provided
+by the
+.Fl l
+option and those used to construct
+the acceleration structure defined by the
+.Fl a
+option.
+.Pp
+Warning!
+It is always recommended to verify that the data is correct, even though
+this verification can take a significant amount of time when there are a
+large number of lines.
+Anyway, a user who is
+.Em certain
+of the data’s consistency may nevertheless use this option
+.Pq at their own risk
+to disable this verification and thus speed up the execution.
+.\""""""""""""""""""""""""""""""""""
+.It Fl h
+Display short help and exit.
+.It Fl l Ar lines
+List of lines from which the tree was built.
+This list is in binary format as generated by the
+.Xr shtr 1
+binary, or in plain text HITRAN format, depending on whether the
+.Fl s
+option is set or not, respectively.
+.\""""""""""""""""""""""""""""""""""
+.It Fl m Ar molparams
+Isotopologue metadata in HITRAN format.
+.\""""""""""""""""""""""""""""""""""
+.It Fl n Ar nrealisations
+Number of Monte Carlo realisations.
+By default the number of realisations is 10000.
+.\""""""""""""""""""""""""""""""""""
+.It Fl S Ar nu_min , Ns Ar nu_max
+The spectral range, in cm^-1, over which the computations are performed.
+The default spectral range is from 0 to infinity.
+.\""""""""""""""""""""""""""""""""""
+.It Fl s
+Specifies that input lines are formatted according to the binary format
+as written by the
+.Xr shtr 1
+utility, and not according to the HITRAN format.
+This format is more compact, allowing for faster loading of line data.
+.\""""""""""""""""""""""""""""""""""
+.It Fl t Ar threads
+Advice on the number of threads to use.
+By default,
+.Nm
+uses as many threads as processor cores.
+.\""""""""""""""""""""""""""""""""""
+.It Fl v
+Make
+.Nm
+verbose.
+Multiple
+.Fl v
+options increase the verbosity.
+The maximum is 3.
+.El
+.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+.Sh EXIT STATUS
+.Ex -std
+.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+.Sh EXAMPLES
+Estimate the mean k and mean square k between 100 and 2500 cm^-1
+for a gaz mixture made of H2O, CO2 and CO molecules.
+The thermodynamic properties of the mixture, such as its pressure,
+temperature and molecular concentrations, correspond to those used to
+construct the acceleration structures with sln-build, provided as input
+arguments
+.Pq option Fl a .
+The isotopic metadata
+.Pq option Fl m
+and the list of lines
+.Pq option Fl l
+partitioned by the acceleration structure, complete the list of input
+data.
+The latter is encoded in the format generated by the
+.Xr shtr 1
+tool
+.Pq option Fl s .
+The isotopes are in HITRAN format.
+Finally, make the program as verbose as possible
+.Pq options Fl vvv .
+.Bd -literal -offset Ds
+sln-stat -S 100,2500 -a tree_H2O_CO2_CO_1atm_600K.sln \e
+ -m molparam.txt -sl H2O_CO2_CO_100-2500cm-1.shtr -vvv
+.Ed
+.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+.Sh SEE ALSO
+.Xr shtr 1 ,
+.Xr sln-build 1 ,
+.Xr sln-slab 1
diff --git a/src/sln.h b/src/sln.h
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
@@ -210,6 +212,15 @@ struct sln_tree_desc {
}
static const struct sln_tree_desc SLN_TREE_DESC_NULL = SLN_TREE_DESC_NULL__;
+struct sln_thermo_props {
+ double concentrations[SHTR_MAX_MOLECULE_COUNT];
+ double pressure; /* [atm] */
+ double temperature; /* [K] */
+};
+#define SLN_THERMO_PROPS_NULL__ {{0},0,0}
+static const struct sln_thermo_props SLN_THERMO_PROPS_NULL =
+ SLN_THERMO_PROPS_NULL__;
+
struct sln_node_desc {
/* Range of lines belonging to the node. The endpoints are included */
size_t ilines[2];
@@ -351,6 +362,9 @@ SLN_API res_T
sln_tree_get_line
(const struct sln_tree* tree,
const size_t iline,
+ /* Thermodynamic properties to which the line is recovered.
+ * Can be NULL, so these properties are those used to build the tree */
+ const struct sln_thermo_props* props,
struct sln_line* line);
SLN_API res_T
@@ -381,6 +395,9 @@ SLN_API double
sln_node_eval
(const struct sln_tree* tree,
const struct sln_node* node,
+ /* Thermodynamic properties to which the node lines are evaluated.
+ * Can be NULL, so these properties are those used to build the tree */
+ const struct sln_thermo_props* props,
const double wavenumber); /* In cm^-1 */
SLN_API res_T
diff --git a/src/sln_build.c b/src/sln_build.c
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
diff --git a/src/sln_device.c b/src/sln_device.c
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
diff --git a/src/sln_device_c.h b/src/sln_device_c.h
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
diff --git a/src/sln_faddeeva.c b/src/sln_faddeeva.c
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
diff --git a/src/sln_get.c b/src/sln_get.c
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
@@ -427,7 +429,7 @@ print_node_value(const struct cmd* cmd)
if(res != RES_OK) goto error;
val_mesh = sln_mesh_eval(&mesh, cmd->args.wavenumber);
- val_node = sln_node_eval(cmd->tree, node, cmd->args.wavenumber);
+ val_node = sln_node_eval(cmd->tree, node, NULL, cmd->args.wavenumber);
printf("ka(%e) = %e ~ %e\n", cmd->args.wavenumber, val_node, val_mesh);
diff --git a/src/sln_line.c b/src/sln_line.c
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
@@ -70,6 +72,9 @@ static res_T
line_profile_factor
(const struct sln_tree* tree,
const struct shtr_line* shtr_line,
+ const double concentration,
+ const double pressure,
+ const double temperature,
double* out_profile_factor)
{
/* Star-HITRAN data */
@@ -103,19 +108,19 @@ line_profile_factor
ASSERT(molecule.nisotopes > (size_t)shtr_line->isotope_id_local);
isotope = molecule.isotopes + shtr_line->isotope_id_local;
- nu_c = line_center(shtr_line, tree->args.pressure);
+ nu_c = line_center(shtr_line, pressure);
/* Compute the intensity */
- Ps = tree->args.pressure * mol_params->concentration;
+ Ps = pressure * concentration;
density = (AVOGADRO_NUMBER * Ps);
- density = density / (PERFECT_GAZ_CONSTANT * tree->args.temperature);
+ density = density / (PERFECT_GAZ_CONSTANT * temperature);
density = density * 1e-6; /* Convert in molec.cm^-3 */
- /* Compute the partition function. TODO precompute it for molid/isoid */
+ /* Compute the partition function */
Q_Tref = isotope->Q296K;
molid = shtr_line->molecule_id;
isoid = shtr_line->isotope_id_local+1/*Local indices start at 1 in BD_TIPS*/;
- T = tree->args.temperature;
+ T = temperature;
BD_TIPS_2017(&molid, &T, &isoid, &gj, &Q_T);
if(Q_T <= 0) {
ERROR(tree->sln,
@@ -135,7 +140,7 @@ line_profile_factor
intensity_ref = shtr_line->intensity/isotope->abundance*iso_abundance;
}
intensity = line_intensity(intensity_ref, shtr_line->lower_state_energy, Q,
- tree->args.temperature, T_REF, nu_c);
+ temperature, T_REF, nu_c);
profile_factor = 1.e2 * density * intensity; /* In m^-1.cm^-1 */
@@ -427,12 +432,15 @@ res_T
line_setup
(const struct sln_tree* tree,
const size_t iline,
+ const struct sln_thermo_props* props,
struct sln_line* line)
{
struct shtr_molecule molecule = SHTR_MOLECULE_NULL;
struct shtr_line shtr_line = SHTR_LINE_NULL;
- double molar_mass = 0; /* In kg.mol^-1 */
- const struct sln_molecule* mol_params = NULL;
+ double concentration = 0;
+ double molar_mass = 0; /*[kg.mol^-1]*/
+ double pressure = 0; /*[atm]*/
+ double temperature = 0; /*[K]*/
res_T res = RES_OK;
ASSERT(tree && line);
@@ -443,21 +451,30 @@ line_setup
ASSERT(!SHTR_MOLECULE_IS_NULL(&molecule));
ASSERT(molecule.nisotopes > (size_t)shtr_line.isotope_id_local);
- mol_params = tree->args.molecules + shtr_line.molecule_id;
+ if(!props) { /* Use thermo properties used to build the tree */
+ concentration = tree->args.molecules[shtr_line.molecule_id].concentration;
+ pressure = tree->args.pressure; /*[atm]*/
+ temperature = tree->args.temperature; /*[K]*/
+ } else {
+ concentration = props->concentrations[shtr_line.molecule_id];
+ pressure = props->pressure; /*[atm]*/
+ temperature = props->temperature; /*[K]*/
+ }
/* Convert the molar mass of the line from g.mol^-1 to kg.mol^-1 */
molar_mass = molecule.isotopes[shtr_line.isotope_id_local].molar_mass*1e-3;
/* Setup the line */
- res = line_profile_factor(tree, &shtr_line, &line->profile_factor);
+ res = line_profile_factor(tree, &shtr_line, concentration, pressure,
+ temperature, &line->profile_factor);
if(res != RES_OK) goto error;
- line->wavenumber = line_center(&shtr_line, tree->args.pressure);
+ line->wavenumber = line_center(&shtr_line, pressure);
line->gamma_d = sln_compute_line_half_width_doppler
- (line->wavenumber, molar_mass, tree->args.temperature);
+ (line->wavenumber, molar_mass, temperature);
line->gamma_l = sln_compute_line_half_width_lorentz
- (shtr_line.gamma_air, shtr_line.gamma_self, tree->args.temperature,
- tree->args.pressure, shtr_line.n_air, mol_params->concentration);
+ (shtr_line.gamma_air, shtr_line.gamma_self, temperature,
+ pressure, shtr_line.n_air, concentration);
line->molecule_id = shtr_line.molecule_id;
exit:
@@ -493,7 +510,7 @@ line_mesh
/* Setup the line wrt molecule concentration, isotope abundance, temperature
* and pressure */
- res = line_setup(tree, iline, &line);
+ res = line_setup(tree, iline, NULL/*default thermo props*/, &line);
if(res != RES_OK) goto error;
/* Adjust the hint on the number of vertices. This is not actually the real
diff --git a/src/sln_line.h b/src/sln_line.h
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
@@ -43,6 +45,9 @@ extern LOCAL_SYM res_T
line_setup
(const struct sln_tree* tree,
const size_t iline,
+ /* Thermodynamic properties to which the line is recovered.
+ * Can be NULL, so these properties are those used to build the tree */
+ const struct sln_thermo_props* props,
struct sln_line* line);
extern LOCAL_SYM res_T
diff --git a/src/sln_mixture.c b/src/sln_mixture.c
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
diff --git a/src/sln_polyline.c b/src/sln_polyline.c
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
diff --git a/src/sln_polyline.h b/src/sln_polyline.h
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
diff --git a/src/sln_slab.c b/src/sln_slab.c
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
@@ -411,7 +413,7 @@ realisation
/* Evaluate the value of the line and compute the probability of being
* absorbed by it */
- leaf_ka = sln_node_eval(cmd->tree, leaf, nu);
+ leaf_ka = sln_node_eval(cmd->tree, leaf, NULL, nu);
proba_abs = leaf_ka / (leaf_proba*ka_max);
if((res = check_proba(cmd, proba_abs)) != RES_OK) goto error;
diff --git a/src/sln_stat.c b/src/sln_stat.c
@@ -0,0 +1,488 @@
+/* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com)
+ * Copyright (C) 2026 Université de Lorraine
+ * Copyright (C) 2022 Centre National de la Recherche Scientifique
+ * Copyright (C) 2022 Université Paul Sabatier
+ *
+ * This file is part of Star-Line.
+ *
+ * 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 /* getopt */
+
+#include "sln.h"
+
+#include <star/shtr.h>
+#include <star/sbb.h>
+#include <star/ssp.h>
+
+#include <rsys/cstr.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/str.h>
+
+#include <omp.h>
+
+#include <unistd.h> /* getopt */
+
+enum estimate {
+ MEAN,
+ SQMEAN,
+ ESTIMATE_COUNT__
+};
+
+#define WAVENUMBER_TO_WAVELENGTH(Nu/* [cm^-1] */) (1.e-2/(Nu))/*[m]*/
+
+struct args {
+ const char* tree; /* Acceleration structure */
+ const char* molparams;
+ const char* lines;
+
+ double spectral_range[2]; /* [cm^-1]^2 */
+
+ unsigned long nrealisations; /* Number of Monte Carlo realisations */
+
+ /* Miscellaneous */
+ unsigned nthreads_hint; /* Hint on the number of threads to use */
+ int disable_line_hash_check;
+ int lines_in_shtr_format;
+ int verbose;
+ int quit;
+};
+#define ARGS_DEFAULT__ {NULL,NULL,NULL,{0,DBL_MAX},10000,UINT_MAX,0,0,0,0}
+static const struct args ARGS_DEFAULT = ARGS_DEFAULT__;
+
+struct cmd {
+ struct args args;
+
+ struct sln_tree* tree;
+ unsigned nthreads;
+};
+#define CMD_NULL__ {0}
+static const struct cmd CMD_NULL = CMD_NULL__;
+
+struct accum {
+ double sum;
+ double sum2;
+ size_t count;
+};
+#define ACCUM_NULL__ {0}
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+usage(FILE* stream)
+{
+ fprintf(stream,
+"usage: sln-stat [-dhsv] [-n nrealisations] [-t threads]\n"
+" -S nu_min,nu_max -a accel_struct -m molparams -l lines\n");
+}
+
+static res_T
+parse_spectral_range(const char* str, double spectral_range[2])
+{
+ size_t len = 0;
+ res_T res = RES_OK;
+ ASSERT(str && spectral_range);
+
+ res = cstr_to_list_double(str, ',', spectral_range, &len, 2);
+ if(res == RES_OK && len < 2) res = RES_BAD_ARG;
+
+ return res;
+}
+
+static res_T
+args_init(struct args* args, int argc, char** argv)
+{
+ int opt = 0;
+ res_T res = RES_OK;
+
+ ASSERT(args);
+
+ *args = ARGS_DEFAULT;
+
+ while((opt = getopt(argc, argv, "a:dhl:m:n:S:st:v")) != -1) {
+ switch(opt) {
+ case 'a': args->tree = optarg; break;
+ case 'd': args->disable_line_hash_check = 1; break;
+ case 'h':
+ usage(stdout);
+ args->quit = 1;
+ goto exit;
+ case 'l': args->lines = optarg; break;
+ case 'm': args->molparams = optarg; break;
+ case 'n': res = cstr_to_ulong(optarg, &args->nrealisations); break;
+ case 'S': res = parse_spectral_range(optarg, args->spectral_range); break;
+ case 's': args->lines_in_shtr_format = 1; break;
+ case 't':
+ res = cstr_to_uint(optarg, &args->nthreads_hint);
+ if(res == RES_OK && args->nthreads_hint == 0) res = RES_BAD_ARG;
+ break;
+ case 'v': args->verbose += (args->verbose < 3); break;
+ default: res = RES_BAD_ARG; break;
+ }
+ if(res != RES_OK) {
+ if(optarg) {
+ fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n",
+ argv[0], optarg, opt);
+ }
+ goto error;
+ }
+ }
+
+ #define MANDATORY(Cond, Name, Opt) { \
+ if(!(Cond)) { \
+ fprintf(stderr, "%s: %s missing -- option '-%c'\n", argv[0], (Name), (Opt)); \
+ res = RES_BAD_ARG; \
+ goto error; \
+ } \
+ } (void)0
+ MANDATORY(args->molparams, "molparams", 'm');
+ MANDATORY(args->lines, "line list", 'l');
+ MANDATORY(args->tree, "acceleration structure", 'a');
+ #undef MANDATORY
+
+exit:
+ return res;
+error:
+ usage(stderr);
+ goto exit;
+}
+
+static res_T
+load_lines
+ (struct shtr* shtr,
+ const struct args* args,
+ struct shtr_line_list** out_lines)
+{
+ struct shtr_line_list* lines = NULL;
+ res_T res = RES_OK;
+ ASSERT(shtr && args && out_lines);
+
+ if(args->lines_in_shtr_format) {
+ struct shtr_line_list_read_args read_args = SHTR_LINE_LIST_READ_ARGS_NULL;
+
+ /* Loads lines from data serialized by the Star-HITRAN library */
+ read_args.filename = args->lines;
+ res = shtr_line_list_read(shtr, &read_args, &lines);
+ if(res != RES_OK) goto error;
+
+ } else {
+ struct shtr_line_list_load_args load_args = SHTR_LINE_LIST_LOAD_ARGS_NULL;
+
+ /* Loads lines from a file in HITRAN format */
+ load_args.filename = args->lines;
+ res = shtr_line_list_load(shtr, &load_args, &lines);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ *out_lines = lines;
+ return res;
+error:
+ if(lines) { SHTR(line_list_ref_put(lines)); lines = NULL; }
+ goto exit;
+}
+
+static void
+delete_per_thread_rngs(const struct cmd* cmd, struct ssp_rng* rngs[])
+{
+ unsigned i = 0;
+ ASSERT(cmd && rngs);
+
+ FOR_EACH(i, 0, cmd->nthreads) {
+ if(rngs[i]) SSP(rng_ref_put(rngs[i]));
+ }
+ mem_rm(rngs);
+}
+
+static res_T
+create_per_thread_rngs(const struct cmd* cmd, struct ssp_rng** out_rngs[])
+{
+ struct ssp_rng_proxy* proxy = NULL;
+ struct ssp_rng** rngs = NULL;
+ size_t i = 0;
+ res_T res = RES_OK;
+ ASSERT(cmd);
+
+ rngs = mem_calloc(cmd->nthreads, sizeof(*rngs));
+ if(!rngs) { res = RES_MEM_ERR; goto error; }
+
+ res = ssp_rng_proxy_create(NULL, SSP_RNG_THREEFRY, cmd->nthreads, &proxy);
+ if(res != RES_OK) goto error;
+
+ FOR_EACH(i, 0, cmd->nthreads) {
+ res = ssp_rng_proxy_create_rng(proxy, i, &rngs[i]);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ *out_rngs = rngs;
+ if(proxy) SSP(rng_proxy_ref_put(proxy));
+ return res;
+error:
+ if(cmd->args.verbose >= 1) {
+ fprintf(stderr,
+ "Error creating the list of per thread RNG -- %s\n",
+ res_to_cstr(res));
+ }
+ if(rngs) delete_per_thread_rngs(cmd, rngs);
+ rngs = NULL;
+ goto exit;
+}
+
+static void
+cmd_release(struct cmd* cmd)
+{
+ ASSERT(cmd);
+ if(cmd->tree) SLN(tree_ref_put(cmd->tree));
+}
+
+static res_T
+cmd_init(struct cmd* cmd, const struct args* args)
+{
+ /* Star Line */
+ struct sln_device_create_args sln_args = SLN_DEVICE_CREATE_ARGS_DEFAULT;
+ struct sln_tree_read_args tree_args = SLN_TREE_READ_ARGS_NULL;
+ struct sln_device* sln = NULL;
+
+ /* Star HITRAN */
+ struct shtr_create_args shtr_args = SHTR_CREATE_ARGS_DEFAULT;
+ struct shtr* shtr = NULL;
+ struct shtr_isotope_metadata* molparams = NULL;
+ struct shtr_line_list* lines = NULL;
+
+ /* Miscellaneous */
+ unsigned nthreads_max = 0;
+ res_T res = RES_OK;
+
+ ASSERT(cmd && args);
+
+ *cmd = CMD_NULL;
+
+ shtr_args.verbose = args->verbose;
+ res = shtr_create(&shtr_args, &shtr);
+ if(res != RES_OK) goto error;
+
+ res = shtr_isotope_metadata_load(shtr, args->molparams, &molparams);
+ if(res != RES_OK) goto error;
+
+ res = load_lines(shtr, args, &lines);
+ if(res != RES_OK) goto error;
+
+ sln_args.verbose = args->verbose;
+ res = sln_device_create(&sln_args, &sln);
+ if(res != RES_OK) goto error;
+
+ tree_args.metadata = molparams;
+ tree_args.lines = lines;
+ tree_args.filename = args->tree;
+ tree_args.disable_line_hash_check = args->disable_line_hash_check;
+ res = sln_tree_read(sln, &tree_args, &cmd->tree);
+ if(res != RES_OK) goto error;
+
+ nthreads_max = (unsigned)MMAX(omp_get_max_threads(), omp_get_num_procs());
+ cmd->args = *args;
+ cmd->nthreads = MMIN(cmd->args.nthreads_hint, nthreads_max);
+
+exit:
+ if(sln) SLN(device_ref_put(sln));
+ if(shtr) SHTR(ref_put(shtr));
+ if(molparams) SHTR(isotope_metadata_ref_put(molparams));
+ if(lines) SHTR(line_list_ref_put(lines));
+ return res;
+error:
+ cmd_release(cmd);
+ *cmd = CMD_NULL;
+ goto exit;
+}
+
+static INLINE const char*
+estimate_cstr(const enum estimate estimate)
+{
+ const char* cstr = NULL;
+ switch(estimate) {
+ case MEAN: cstr="mean"; break;
+ case SQMEAN: cstr="mean_of_squares"; break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ return cstr;
+}
+
+static res_T
+realisation
+ (const struct cmd* cmd,
+ struct ssp_rng* rng,
+ double out_weights[ESTIMATE_COUNT__])
+{
+ /* Acceleration structure */
+ struct sln_tree_desc tree_desc = SLN_TREE_DESC_NULL;
+ const struct sln_node* root = NULL;
+
+ /* Variables to sample k */
+ const struct sln_node* leaf1 = NULL;
+ const struct sln_node* leaf2 = NULL;
+ double leaf_proba1 = 0; /* Probability of sampling a line */
+ double leaf_proba2 = 0; /* Probability of sampling a line */
+ double leaf_ka1 = 0; /* Value of a line */
+ double leaf_ka2 = 0; /* Value of a line */
+
+ /* Miscellaneous */
+ double w[ESTIMATE_COUNT__] = {0, 0}; /* Monte Carlo weight */
+ double nu = 0; /* Sampled wavenumber [cm^-1] */
+ int i = 0;
+ res_T res = RES_OK;
+
+ ASSERT(cmd && rng && out_weights); /* Check pre-conditions */
+
+ /* Uniformly sample the spectral dimension */
+ nu = ssp_rng_uniform_double
+ (rng, cmd->args.spectral_range[0], cmd->args.spectral_range[1]);
+
+ SLN(tree_get_desc(cmd->tree, &tree_desc));
+
+ /* Store the root node of the tree */
+ root = sln_tree_get_root(cmd->tree);
+
+ /* Importance sampling of a line and evaluation of the line contribution */
+ leaf1 = sln_node_sample_leaf(cmd->tree, root, nu, rng, &leaf_proba1);
+ if(!leaf1) { res = RES_BAD_ARG; goto error; }
+ leaf_ka1 = sln_node_eval(cmd->tree, leaf1, NULL, nu);
+
+ /* Importance sampling of a line and evaluation of the line contribution */
+ leaf2 = sln_node_sample_leaf(cmd->tree, root, nu, rng, &leaf_proba2);
+ if(!leaf2) { res = RES_BAD_ARG; goto error; }
+ leaf_ka2 = sln_node_eval(cmd->tree, leaf2, NULL, nu);
+
+ w[MEAN] = leaf_ka1 / leaf_proba1;
+ w[SQMEAN] = (leaf_ka1/leaf_proba1)*(leaf_ka2/leaf_proba2);
+
+exit:
+ FOR_EACH(i, 0, ESTIMATE_COUNT__) out_weights[i] = w[i];
+ return res;
+error:
+ FOR_EACH(i, 0, ESTIMATE_COUNT__) w[i] = NaN;
+ goto exit;
+}
+
+static res_T
+cmd_run(const struct cmd* cmd)
+{
+ /* Random Number Generator */
+ struct ssp_rng** rngs = NULL;
+
+ /* Monte Carlo */
+ struct accum accum[ESTIMATE_COUNT__] = {0};
+ int64_t i = 0; /* Index of the realisation */
+ size_t nrejects = 0; /* Number of rejected realisations */
+
+ /* Progress */
+ size_t nrealisations = 0;
+ size_t realisation_done = 0;
+ int progress = 0;
+ int progress_pcent = 10;
+
+ res_T res = RES_OK;
+ ASSERT(cmd);
+
+ res = create_per_thread_rngs(cmd, &rngs);
+ if(res != RES_OK) goto error;
+
+ #define PROGRESS_MSG "Solving: %3d%%\n"
+ if(cmd->args.verbose >= 3) fprintf(stderr, PROGRESS_MSG, progress);
+
+ nrealisations = cmd->args.nrealisations;
+
+ omp_set_num_threads((int)cmd->nthreads);
+
+ #pragma omp parallel for schedule(static)
+ for(i = 0; i < (int64_t)nrealisations; ++i) {
+ double w[ESTIMATE_COUNT__] = {0}; /* Monte Carlo weights */
+ const int ithread = omp_get_thread_num();
+ int pcent = 0;
+ res_T res_realisation = RES_OK;
+
+ res_realisation = realisation(cmd, rngs[ithread], w);
+
+ #pragma omp critical
+ {
+ /* Update the Monte Carlo accumulator */
+ if(res_realisation == RES_OK) {
+ int iestim = 0;
+ FOR_EACH(iestim, 0, ESTIMATE_COUNT__) {
+ accum[iestim].sum += w[iestim];
+ accum[iestim].sum2 += w[iestim]*w[iestim];
+ accum[iestim].count += 1;
+ }
+ }
+
+ if(cmd->args.verbose >= 3) {
+ /* Update progress */
+ realisation_done += 1;
+ pcent = (int)((double)realisation_done*100.0/(double)nrealisations+0.5);
+ if(pcent/progress_pcent > progress/progress_pcent) {
+ progress = pcent;
+ fprintf(stderr, PROGRESS_MSG, progress);
+ }
+ }
+ }
+ }
+
+ #undef PROGRESS_MSG
+
+ nrejects = nrealisations - accum[0].count;
+
+ FOR_EACH(i, 0, ESTIMATE_COUNT__) {
+ const double E = accum[i].sum / (double)accum[i].count;
+ const double V = accum[i].sum2 / (double)accum[i].count - E*E;
+ const double SE = sqrt(V/(double)accum[i].count);
+
+ /* Assume that the number of realisations is the same for all estimates */
+ ASSERT(accum[i].count == accum[0].count);
+
+ printf("%-16s: %e +/- %e; %lu\n",
+ estimate_cstr(i), E, SE, (unsigned long)nrejects);
+ }
+
+exit:
+ delete_per_thread_rngs(cmd, rngs);
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Main function
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct args args = ARGS_DEFAULT;
+ struct cmd cmd = CMD_NULL;
+ int err = 0;
+ res_T res = RES_OK;
+
+ if((res = args_init(&args, argc, argv)) != RES_OK) goto error;
+ if(args.quit) goto exit;
+
+ if((res = cmd_init(&cmd, &args)) != RES_OK) goto error;
+ if((res = cmd_run(&cmd)) != RES_OK) goto error;
+
+exit:
+ cmd_release(&cmd);
+ CHK(mem_allocated_size() == 0);
+ return err;
+error:
+ err = 1;
+ goto exit;
+}
diff --git a/src/sln_tree.c b/src/sln_tree.c
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
@@ -40,20 +42,52 @@ static const struct stream STREAM_NULL = {NULL, NULL, 0};
/*******************************************************************************
* Helper functions
******************************************************************************/
-/* Check that the sum of the molecular concentrations is equal to 1 */
-static res_T
+static INLINE res_T
check_molecule_concentration
- (struct sln_device* sln,
+ (const struct sln_device* sln,
const char* caller,
- const struct sln_tree_create_args* args)
+ const enum shtr_molecule_id molecule_id,
+ const double concentration)
+{
+ ASSERT(sln && caller);
+
+ if(concentration == 0) {
+ /* A molecular concentration of zero is allowed, but may be a user error,
+ * as 0 is the default concentration in the tree creation arguments.
+ * Therefore, warn the user about this value so that they can determine
+ * whether or not it is an error on their part. */
+ WARN(sln, "%s: the concentration of %s is zero.\n",
+ caller, shtr_molecule_cstr(molecule_id));
+
+ } else if(concentration < 0) {
+ /* Concentration cannot be negative... */
+ ERROR(sln, "%s: invalid %s concentration: %g.\n",
+ FUNC_NAME, shtr_molecule_cstr(molecule_id),
+ concentration);
+ return RES_BAD_ARG;
+ }
+
+ return RES_OK;
+}
+
+static res_T
+check_concentrations_list
+ (const struct sln_device* sln,
+ const char* caller,
+ const double concentrations[SHTR_MAX_MOLECULE_COUNT])
{
- int i = 0;
double sum = 0;
- ASSERT(sln && caller && args);
+ int i = 0;
+ res_T res = RES_OK;
+ ASSERT(sln && caller && concentrations);
FOR_EACH(i, 0, SHTR_MAX_MOLECULE_COUNT) {
if(i == SHTR_MOLECULE_ID_NULL) continue;
- sum += args->molecules[i].concentration;
+
+ res = check_molecule_concentration(sln, caller, i, concentrations[i]);
+ if(res != RES_OK) goto error;
+
+ sum += concentrations[i];
}
/* The sum of molecular concentrations must be less than or equal to 1. It may
@@ -63,10 +97,32 @@ check_molecule_concentration
ERROR(sln,
"%s: the sum of molecule concentrations is greater than 1: %g\n",
caller, sum);
- return RES_BAD_ARG;
+ res = RES_BAD_ARG;
+ goto error;
}
- return RES_OK;
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/* Check the consistency of the molecular concentrations */
+static INLINE res_T
+check_mixture_concentrations
+ (struct sln_device* sln,
+ const char* caller,
+ const struct sln_tree_create_args* args)
+{
+ double concentrations[SHTR_MAX_MOLECULE_COUNT] = {0};
+ int i = 0;
+ ASSERT(sln && caller && args);
+
+ FOR_EACH(i, 0, SHTR_MAX_MOLECULE_COUNT) {
+ concentrations[i] = args->molecules[i].concentration;
+ }
+
+ return check_concentrations_list(sln, caller, concentrations);
}
/* Verify that the isotope abundance are valids */
@@ -115,14 +171,13 @@ check_molecules
{
char molecule_ok[SHTR_MAX_MOLECULE_COUNT] = {0};
- double concentrations_sum = 0;
size_t iline = 0;
size_t nlines = 0;
res_T res = RES_OK;
ASSERT(args->lines);
- res = check_molecule_concentration(sln, caller, args);
- if(res != RES_OK) return res;
+ res = check_mixture_concentrations(sln, caller, args);
+ if(res != RES_OK) goto error;
/* Iterate over the lines to define which molecules has to be checked, i.e.,
* the ones used in the mixture */
@@ -138,24 +193,6 @@ check_molecules
molecule = args->molecules + line.molecule_id;
- if(molecule->concentration == 0) {
- /* A molecular concentration of zero is allowed, but may be a user error,
- * as 0 is the default concentration in the tree creation arguments.
- * Therefore, warn the user about this value so that they can determine
- * whether or not it is an error on their part. */
- WARN(sln, "%s: the concentration of %s is zero.\n",
- caller, shtr_molecule_cstr(line.molecule_id));
-
- } else if(molecule->concentration < 0) {
- /* Concentration cannot be negative... */
- ERROR(sln, "%s: invalid %s concentration: %g.\n",
- FUNC_NAME, shtr_molecule_cstr(line.molecule_id),
- molecule->concentration);
- return RES_BAD_ARG;
- }
-
- concentrations_sum += molecule->concentration;
-
if(molecule->cutoff <= 0) {
/* ... cutoff either */
ERROR(sln, "%s: invalid %s cutoff: %g.\n",
@@ -164,21 +201,40 @@ check_molecules
}
res = check_molecule_isotope_abundances(sln, caller, molecule);
- if(res != RES_OK) return res;
+ if(res != RES_OK) goto error;
molecule_ok[line.molecule_id] = 1;
}
- /* The sum of molecular concentrations must be less than or equal to 1. It may
- * be less than 1 if the remaining part of the mixture is (implicitly) defined
- * as a radiatively inactive gas */
- if(concentrations_sum > 1 && (concentrations_sum - 1) > 1e-6) {
- ERROR(sln,
- "%s: the sum of molecule concentrations is greater than 1: %g\n",
- caller, concentrations_sum);
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static INLINE res_T
+check_pressure
+ (const struct sln_device* sln,
+ const char* caller,
+ const double pressure /*[atm]*/)
+{
+ if(pressure < 0) {
+ ERROR(sln, "%s: invalid negative pressure %g atm\n", caller, pressure);
return RES_BAD_ARG;
}
+ return RES_OK;
+}
+static INLINE res_T
+check_temperature
+ (const struct sln_device* sln,
+ const char* caller,
+ const double temperature /*[K]*/)
+{
+ if(temperature < 0) {
+ ERROR(sln, "%s: invalid negative temperature %g K\n", caller, temperature);
+ return RES_BAD_ARG;
+ }
return RES_OK;
}
@@ -203,6 +259,14 @@ check_sln_tree_create_args
return RES_BAD_ARG;
}
+ if((res = check_pressure(sln, caller, args->pressure)) != RES_OK) {
+ return res;
+ }
+
+ if((res = check_temperature(sln, caller, args->temperature)) != RES_OK) {
+ return res;
+ }
+
if(args->nvertices_hint == 0) {
ERROR(sln,
"%s: invalid hint on the number of vertices around the line center %lu.\n",
@@ -296,6 +360,64 @@ check_sln_tree_write_args
return RES_OK;
}
+
+static res_T
+check_line_thermo_props
+ (const struct sln_tree* tree,
+ const char* caller,
+ const size_t iline,
+ const struct sln_thermo_props* props)
+{
+ struct shtr_line line = SHTR_LINE_NULL;
+ res_T res = RES_OK;
+ ASSERT(tree && caller);
+
+ if(!props) goto exit; /* Default thermo props */
+
+ SHTR(line_list_at(tree->args.lines, iline, &line));
+
+ res = check_molecule_concentration(tree->sln, caller, line.molecule_id,
+ props->concentrations[line.molecule_id]);
+ if(res != RES_OK) goto error;
+
+ res = check_pressure(tree->sln, caller, props->pressure);
+ if(res != RES_OK) goto error;
+
+ res = check_temperature(tree->sln, caller, props->temperature);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static INLINE res_T
+check_sln_thermo_props
+ (const struct sln_device* sln,
+ const char* caller,
+ const struct sln_thermo_props* props)
+{
+ res_T res = RES_OK;
+ ASSERT(sln && caller);
+
+ if(!props) goto exit; /* Default thermo props */
+
+ res = check_concentrations_list(sln, caller, props->concentrations);
+ if(res != RES_OK) goto error;
+
+ res = check_pressure(sln, caller, props->pressure);
+ if(res != RES_OK) goto error;
+
+ res = check_temperature(sln, caller, props->temperature);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
static INLINE void
stream_release(struct stream* stream)
{
@@ -703,6 +825,7 @@ res_T
sln_tree_get_line
(const struct sln_tree* tree,
const size_t iline,
+ const struct sln_thermo_props* props,
struct sln_line* line)
{
size_t nlines = 0;
@@ -713,7 +836,10 @@ sln_tree_get_line
SHTR(line_list_get_size(tree->args.lines, &nlines));
if(iline >= nlines) { res = RES_BAD_ARG; goto error; }
- res = line_setup(tree, iline, line);
+ res = check_line_thermo_props(tree, FUNC_NAME, iline, props);
+ if(res != RES_OK) goto error;
+
+ res = line_setup(tree, iline, props, line);
if(res != RES_OK) {
ERROR(tree->sln, "%s: could not setup the line %lu-- %s\n",
FUNC_NAME, iline, res_to_cstr(res));
@@ -775,17 +901,19 @@ double
sln_node_eval
(const struct sln_tree* tree,
const struct sln_node* node,
+ const struct sln_thermo_props* props,
const double nu)
{
double ka = 0;
size_t iline;
ASSERT(tree && node);
+ ASSERT(check_sln_thermo_props(tree->sln, FUNC_NAME, props) == RES_OK);
FOR_EACH(iline, node->range[0], node->range[1]+1) {
struct sln_line line = SLN_LINE_NULL;
res_T res = RES_OK;
- res = line_setup(tree, iline, &line);
+ res = line_setup(tree, iline, props, &line);
if(res != RES_OK) {
WARN(tree->sln, "%s: could not setup the line %lu-- %s\n",
FUNC_NAME, iline, res_to_cstr(res));
diff --git a/src/sln_tree_build.c b/src/sln_tree_build.c
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
@@ -270,8 +272,9 @@ build_leaf_polyline_from_Nlines
const size_t iline = leaf->range[0] + i;
/* Mesh the line in the temporary vertex buffer */
- res = line_mesh(tree, iline, tree->args.nvertices_hint, &scratch->vertices,
- vertices_range);
+ res = line_mesh
+ (/* in */ tree, iline, tree->args.nvertices_hint,
+ /* out */ &scratch->vertices, vertices_range);
if(res != RES_OK) goto error;
/* Decimate the line mesh */
diff --git a/src/sln_tree_c.h b/src/sln_tree_c.h
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
diff --git a/src/test_sln_device.c b/src/test_sln_device.c
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
diff --git a/src/test_sln_lines.h b/src/test_sln_lines.h
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
diff --git a/src/test_sln_mesh.c b/src/test_sln_mesh.c
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
diff --git a/src/test_sln_mixture.c b/src/test_sln_mixture.c
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
diff --git a/src/test_sln_thermo_props.c b/src/test_sln_thermo_props.c
@@ -0,0 +1,379 @@
+/* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com)
+ * Copyright (C) 2026 Université de Lorraine
+ * Copyright (C) 2022 Centre National de la Recherche Scientifique
+ * Copyright (C) 2022 Université Paul Sabatier
+ *
+ * This file is part of Star-Line.
+ *
+ * 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 <test_sln_lines.h>
+
+#include "sln.h"
+
+#include <rsys/mem_allocator.h>
+#include <rsys/rsys.h>
+
+struct thermo_props {
+ double xH2O;
+ double xCO2;
+ double xO3;
+ double pressure; /*[atm]*/
+ double temperature; /*[K]*/
+};
+
+static const struct thermo_props thermo_props1 = {0.15, 0.10, 0.05, 10, 600};
+static const struct thermo_props thermo_props2 = {0.10, 0.15, 0.25, 20, 450};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static struct shtr_isotope_metadata*
+setup_isotopes(struct shtr* shtr)
+{
+ struct shtr_isotope_metadata* metadata = NULL;
+ FILE* fp = NULL;
+
+ CHK(fp = tmpfile());
+ fprintf(fp, "Molecule # Iso Abundance Q(296K) gj Molar Mass(g)\n");
+ write_shtr_molecule(fp, &g_H2O);
+ write_shtr_molecule(fp, &g_CO2);
+ write_shtr_molecule(fp, &g_O3);
+ rewind(fp);
+
+ CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &metadata) == RES_OK);
+
+ CHK(fclose(fp) == 0);
+
+ return metadata;
+}
+
+static struct shtr_line_list*
+setup_lines(struct shtr* shtr)
+{
+ struct shtr_line_list_load_args args = SHTR_LINE_LIST_LOAD_ARGS_NULL;
+ struct shtr_line_list* lines = NULL;
+ FILE* fp = NULL;
+
+ CHK(fp = tmpfile());
+ write_shtr_lines(fp, g_lines, g_nlines);
+ rewind(fp);
+
+ args.filename = "stream";
+ args.file = fp;
+ CHK(shtr_line_list_load(shtr, &args, &lines) == RES_OK);
+
+ CHK(fclose(fp) == 0);
+
+ return lines;
+}
+
+static struct sln_tree*
+create_tree
+ (struct sln_device* sln,
+ struct shtr_isotope_metadata* mdata,
+ struct shtr_line_list* lines,
+ const struct thermo_props* props)
+{
+ struct sln_tree_create_args args = SLN_TREE_CREATE_ARGS_DEFAULT;
+ struct sln_tree* tree = NULL;
+
+ args.metadata = mdata;
+ args.lines = lines;
+
+ args.molecules[SHTR_H2O].concentration = props->xH2O;
+ args.molecules[SHTR_H2O].cutoff = 25; /* [cm^-1] */
+ args.molecules[SHTR_CO2].concentration = props->xCO2;
+ args.molecules[SHTR_CO2].cutoff = 50; /* [cm^-1] */
+ args.molecules[SHTR_O3].concentration = props->xO3;
+ args.molecules[SHTR_O3].cutoff = 25; /* [cm^-1] */
+
+ args.pressure = props->pressure; /*[atm]*/
+ args.temperature = props->temperature; /*[K]*/
+
+ CHK(sln_tree_create(sln, &args, &tree) == RES_OK);
+ return tree;
+}
+
+static INLINE double /* in [0,1[ */
+rand_canonic(void)
+{
+ return (double)rand() / (double)((long)RAND_MAX+1);
+}
+
+static double /* [cm^-1] */
+line_sample_nu(const struct sln_tree* tree, const size_t iline)
+{
+ struct sln_line line = SLN_LINE_NULL;
+ double nu_range[2] = {0,0}; /* [cm^-1] */
+ double nu = 0; /* [cm^-1] */
+
+ CHK(sln_tree_get_line(tree, iline, NULL, &line) == RES_OK);
+ nu_range[0] = line.wavenumber - 20;
+ nu_range[1] = line.wavenumber + 20;
+ nu = nu_range[0] + rand_canonic() * (nu_range[1] - nu_range[0]);
+ return nu;
+}
+
+static double /* [cm^-1] */
+node_sample_nu(const struct sln_tree* tree, const struct sln_node* node)
+{
+ struct sln_line line = SLN_LINE_NULL;
+ struct sln_node_desc desc = SLN_NODE_DESC_NULL;
+ double nu_range[2] = {0,0}; /* [cm^-1] */
+ double nu = 0; /* [cm^-1] */
+
+ CHK(sln_node_get_desc(tree, node, &desc) == RES_OK);
+
+ CHK(sln_tree_get_line(tree, desc.ilines[0], NULL, &line) == RES_OK);
+ nu_range[0] = line.wavenumber - 10;
+ CHK(sln_tree_get_line(tree, desc.ilines[1], NULL, &line) == RES_OK);
+ nu_range[1] = line.wavenumber + 10;
+
+ nu = nu_range[0] + rand_canonic() * (nu_range[1] - nu_range[0]);
+ return nu;
+}
+
+static size_t
+node_sample_line(const struct sln_tree* tree, const struct sln_node* node)
+{
+ struct sln_node_desc desc = SLN_NODE_DESC_NULL;
+ const double r = rand_canonic();
+ size_t iline = 0;
+
+ CHK(sln_node_get_desc(tree, node, &desc) == RES_OK);
+ iline = desc.ilines[0];
+ iline += (size_t)(r * (double)(desc.ilines[1] - desc.ilines[0] + 1));
+ return iline;
+}
+
+/* Check that, even belonging to 2 trees built from different thermodynamic
+ * properties, a line has the _exact_ same value when queried with the same
+ * thermodynamic properties. */
+static void
+cmp_lines_values
+ (const struct sln_tree* tree1,
+ const struct sln_tree* tree2,
+ const size_t iline,
+ const double nu/*[cm^-1]*/)
+{
+ struct sln_line line1 = SLN_LINE_NULL;
+ struct sln_line line2 = SLN_LINE_NULL;
+ struct sln_thermo_props props = SLN_THERMO_PROPS_NULL;
+ double ka1 = 0;
+ double ka2 = 0;
+
+ CHK(sln_tree_get_line(tree1, iline, NULL, &line1) == RES_OK);
+ CHK(sln_tree_get_line(tree2, iline, NULL, &line2) == RES_OK);
+ ka1 = sln_line_eval(tree1, &line1, nu);
+ ka2 = sln_line_eval(tree2, &line2, nu);
+ CHK(ka1 != ka2);
+
+ props.concentrations[SHTR_H2O] = thermo_props1.xH2O;
+ props.concentrations[SHTR_CO2] = thermo_props1.xCO2;
+ props.concentrations[SHTR_O3] = thermo_props1.xO3;
+ props.pressure = thermo_props1.pressure;
+ props.temperature = thermo_props1.temperature;
+ CHK(sln_tree_get_line(tree1, iline, NULL, &line1) == RES_OK);
+ CHK(sln_tree_get_line(tree2, iline, &props, &line2) == RES_OK);
+ ka1 = sln_line_eval(tree1, &line1, nu);
+ ka2 = sln_line_eval(tree2, &line2, nu);
+ CHK(ka1 == ka2);
+
+ props.concentrations[SHTR_H2O] = thermo_props2.xH2O;
+ props.concentrations[SHTR_CO2] = thermo_props2.xCO2;
+ props.concentrations[SHTR_O3] = thermo_props2.xO3;
+ props.pressure = thermo_props2.pressure;
+ props.temperature = thermo_props2.temperature;
+ CHK(sln_tree_get_line(tree1, iline, &props, &line1) == RES_OK);
+ CHK(sln_tree_get_line(tree2, iline, NULL, &line2) == RES_OK);
+ ka1 = sln_line_eval(tree1, &line1, nu);
+ ka2 = sln_line_eval(tree2, &line2, nu);
+ CHK(ka1 == ka2);
+
+ props.concentrations[SHTR_H2O] = thermo_props1.xH2O;
+ props.concentrations[SHTR_CO2] = thermo_props2.xCO2;
+ props.concentrations[SHTR_O3] = thermo_props1.xO3;
+ props.pressure = thermo_props2.pressure;
+ props.temperature = thermo_props1.temperature;
+ CHK(sln_tree_get_line(tree1, iline, &props, &line1) == RES_OK);
+ CHK(sln_tree_get_line(tree2, iline, &props, &line2) == RES_OK);
+ ka1 = sln_line_eval(tree1, &line1, nu);
+ ka2 = sln_line_eval(tree2, &line2, nu);
+ CHK(ka1 == ka2);
+}
+
+static void
+cmp_lines
+ (const struct sln_tree* tree1,
+ const struct sln_tree* tree2,
+ const size_t iline)
+{
+ const size_t N = 50;
+ size_t i = 0;
+
+ FOR_EACH(i, 0, N) {
+ const double nu = line_sample_nu(tree1, iline);
+ cmp_lines_values(tree1, tree2, iline, nu);
+ }
+}
+
+static void
+cmp_nodes_lines
+ (const struct sln_tree* tree1,
+ const struct sln_node* node1,
+ const struct sln_tree* tree2,
+ const struct sln_node* node2)
+{
+ const size_t N = 50;
+ size_t i = 0;
+ (void)node2;
+
+ FOR_EACH(i, 0, N) {
+ const size_t iline = node_sample_line(tree1, node1);
+ cmp_lines(tree1, tree2, iline);
+ }
+}
+
+/* Check that, even belonging to 2 trees built from different thermodynamic
+ * properties, a node has the _exact_ same value when evaluated with the same
+ * thermodynamic properties. */
+static void
+cmp_nodes_values
+ (const struct sln_tree* tree1,
+ const struct sln_node* node1,
+ const struct sln_tree* tree2,
+ const struct sln_node* node2,
+ const double nu/*[cm^-1]*/)
+{
+ struct sln_thermo_props props = SLN_THERMO_PROPS_NULL;
+ double ka1 = 0;
+ double ka2 = 0;
+
+ ka1 = sln_node_eval(tree1, node1, NULL, nu);
+ ka2 = sln_node_eval(tree2, node2, NULL, nu);
+ CHK(ka1 != ka2);
+
+ props.concentrations[SHTR_H2O] = thermo_props1.xH2O;
+ props.concentrations[SHTR_CO2] = thermo_props1.xCO2;
+ props.concentrations[SHTR_O3] = thermo_props1.xO3;
+ props.pressure = thermo_props1.pressure;
+ props.temperature = thermo_props1.temperature;
+ ka1 = sln_node_eval(tree1, node1, NULL, nu);
+ ka2 = sln_node_eval(tree2, node2, &props, nu);
+ CHK(ka1 == ka2);
+
+ props.concentrations[SHTR_H2O] = thermo_props2.xH2O;
+ props.concentrations[SHTR_CO2] = thermo_props2.xCO2;
+ props.concentrations[SHTR_O3] = thermo_props2.xO3;
+ props.pressure = thermo_props2.pressure;
+ props.temperature = thermo_props2.temperature;
+ ka1 = sln_node_eval(tree1, node1, &props, nu);
+ ka2 = sln_node_eval(tree2, node2, NULL, nu);
+ CHK(ka1 == ka2);
+
+ props.concentrations[SHTR_H2O] = thermo_props1.xH2O;
+ props.concentrations[SHTR_CO2] = thermo_props2.xCO2;
+ props.concentrations[SHTR_O3] = thermo_props1.xO3;
+ props.pressure = thermo_props2.pressure;
+ props.temperature = thermo_props1.temperature;
+ ka1 = sln_node_eval(tree1, node1, &props, nu);
+ ka2 = sln_node_eval(tree2, node2, &props, nu);
+ CHK(ka1 == ka2);
+}
+
+static void
+cmp_nodes
+ (const struct sln_tree* tree1,
+ const struct sln_node* node1,
+ const struct sln_tree* tree2,
+ const struct sln_node* node2)
+{
+ const size_t N = 50;
+ size_t i = 0;
+
+ FOR_EACH(i, 0, N) {
+ const double nu = node_sample_nu(tree1, node1);
+ cmp_nodes_values(tree1, node1, tree2, node2, nu);
+ }
+
+ cmp_nodes_lines(tree1, node1, tree2, node2);
+}
+
+static void
+cmp_trees(const struct sln_tree* tree1, const struct sln_tree* tree2)
+{
+ const struct sln_node* node1 = NULL;
+ const struct sln_node* node2 = NULL;
+ struct sln_tree_desc desc = SLN_TREE_DESC_NULL;
+
+ node1 = sln_tree_get_root(tree1);
+ node2 = sln_tree_get_root(tree2);
+
+ CHK(sln_tree_get_desc(tree1, &desc) == RES_OK);
+ CHK(desc.arity == 2); /* Assume that the arity of the tree is 2 */
+
+ for(;;) {
+ unsigned ichild = 0;
+
+ cmp_nodes(tree1, node1, tree2, node2);
+ if(sln_node_is_leaf(node1)) break;
+
+ /* Randomly choose one node child */
+ ichild = rand_canonic() < 0.5 ? 0 : 1;
+ node1 = sln_node_get_child(tree1, node1, ichild);
+ node2 = sln_node_get_child(tree2, node2, ichild);
+ }
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(void)
+{
+ struct shtr_create_args shtr_args = SHTR_CREATE_ARGS_DEFAULT;
+ struct shtr* shtr = NULL;
+ struct shtr_isotope_metadata* mdata = NULL;
+ struct shtr_line_list* lines = NULL;
+
+ struct sln_device_create_args sln_args = SLN_DEVICE_CREATE_ARGS_DEFAULT;
+ struct sln_device* sln = NULL;
+ struct sln_tree* tree1 = NULL;
+ struct sln_tree* tree2 = NULL;
+
+ shtr_args.verbose = 1;
+ CHK(shtr_create(&shtr_args, &shtr) == RES_OK);
+ sln_args.verbose = 1;
+ CHK(sln_device_create(&sln_args, &sln) == RES_OK);
+
+ mdata = setup_isotopes(shtr);
+ lines = setup_lines(shtr);
+
+ tree1 = create_tree(sln, mdata, lines, &thermo_props1);
+ tree2 = create_tree(sln, mdata, lines, &thermo_props2);
+
+ cmp_trees(tree1, tree2);
+
+ CHK(shtr_ref_put(shtr) == RES_OK);
+ CHK(shtr_line_list_ref_put(lines) == RES_OK);
+ CHK(shtr_isotope_metadata_ref_put(mdata) == RES_OK);
+
+ CHK(sln_device_ref_put(sln) == RES_OK);
+ CHK(sln_tree_ref_put(tree1) == RES_OK);
+ CHK(sln_tree_ref_put(tree2) == RES_OK);
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sln_tree.c b/src/test_sln_tree.c
@@ -3,6 +3,8 @@
* Copyright (C) 2022 Centre National de la Recherche Scientifique
* Copyright (C) 2022 Université Paul Sabatier
*
+ * This file is part of Star-Line.
+ *
* 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
@@ -28,50 +30,6 @@
/*******************************************************************************
* Helper function
******************************************************************************/
-/* Return the index of the line in the line_list or SIZE_MAX if the line does
- * not exist */
-static INLINE size_t
-find_line
- (struct shtr_line_list* line_list,
- const struct shtr_line* line)
-{
- struct shtr_line ln = SHTR_LINE_NULL;
- size_t lo, hi, mid;
- size_t iline;
- size_t nlines;
-
- CHK(shtr_line_list_get_size(line_list, &nlines) == RES_OK);
-
- /* Dichotomic search */
- lo = 0; hi = nlines -1;
- while(lo < hi) {
- mid = (lo+hi)/2;
-
- CHK(shtr_line_list_at(line_list, mid, &ln) == RES_OK);
- if(line->wavenumber > ln.wavenumber) {
- lo = mid + 1;
- } else {
- hi = mid;
- }
- }
- iline = lo;
-
- CHK(shtr_line_list_at(line_list, iline, &ln) == RES_OK);
- if(ln.wavenumber != line->wavenumber) return SIZE_MAX;
-
-
- /* Find a line with the same wavenumber as the one searched for and whose
- * other member variables are also equal to those of the line searched for */
- while(ln.wavenumber == line->wavenumber
- && !shtr_line_eq(&ln, line)
- && iline < nlines) {
- iline += 1;
- CHK(shtr_line_list_at(line_list, iline, &ln) == RES_OK);
- }
-
- return shtr_line_eq(&ln, line) ? iline : SIZE_MAX;
-}
-
/* This test assumes that all the lines contained into the list are
* partitionned in the tree */
static void
@@ -179,8 +137,8 @@ check_node_equality
FOR_EACH(iline, desc1.ilines[0], desc1.ilines[1]+1) {
struct sln_line line1 = SLN_LINE_NULL;
struct sln_line line2 = SLN_LINE_NULL;
- CHK(sln_tree_get_line(tree1, iline, &line1) == RES_OK);
- CHK(sln_tree_get_line(tree2, iline, &line2) == RES_OK);
+ CHK(sln_tree_get_line(tree1, iline, NULL, &line1) == RES_OK);
+ CHK(sln_tree_get_line(tree2, iline, NULL, &line2) == RES_OK);
CHK(line1.wavenumber == line2.wavenumber);
CHK(line1.profile_factor == line2.profile_factor);
@@ -267,15 +225,15 @@ check_node_value(const struct sln_tree* tree, const struct sln_node* node)
CHK(sln_node_get_desc(tree, node, &desc) == RES_OK);
- CHK(sln_tree_get_line(tree, desc.ilines[0], &line) == RES_OK);
+ CHK(sln_tree_get_line(tree, desc.ilines[0], NULL, &line) == RES_OK);
nu = line.wavenumber + 10;
- CHK(sln_tree_get_line(tree, desc.ilines[1], &line) == RES_OK);
+ CHK(sln_tree_get_line(tree, desc.ilines[1], NULL, &line) == RES_OK);
nu += line.wavenumber - 10;
nu *= 0.5;
- ka_node = sln_node_eval(tree, node, nu);
+ ka_node = sln_node_eval(tree, node, NULL, nu);
FOR_EACH(iline, desc.ilines[0], desc.ilines[1]+1/*inclusive*/) {
- CHK(sln_tree_get_line(tree, iline, &line) == RES_OK);
+ CHK(sln_tree_get_line(tree, iline, NULL, &line) == RES_OK);
ka_ref += sln_line_eval(tree, &line, nu);
}
@@ -329,11 +287,11 @@ test_tree
CHK(desc.temperature == tree_args.temperature);
CHK(desc.arity == tree_args.arity);
- CHK(sln_tree_get_line(NULL, 0, &line) == RES_BAD_ARG);
- CHK(sln_tree_get_line(tree, nlines, &line) == RES_BAD_ARG);
- CHK(sln_tree_get_line(tree, 0, NULL) == RES_BAD_ARG);
- CHK(sln_tree_get_line(tree, 0, &line) == RES_OK);
- CHK(sln_tree_get_line(tree, nlines-1, &line) == RES_OK);
+ CHK(sln_tree_get_line(NULL, 0, NULL, &line) == RES_BAD_ARG);
+ CHK(sln_tree_get_line(tree, nlines, NULL, &line) == RES_BAD_ARG);
+ CHK(sln_tree_get_line(tree, 0, NULL, NULL) == RES_BAD_ARG);
+ CHK(sln_tree_get_line(tree, 0, NULL, &line) == RES_OK);
+ CHK(sln_tree_get_line(tree, nlines-1, NULL, &line) == RES_OK);
CHK(node = sln_tree_get_root(tree));
CHK(node != NULL);
diff --git a/src/test_sln_tree_sample.c b/src/test_sln_tree_sample.c
@@ -0,0 +1,161 @@
+/* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com)
+ * Copyright (C) 2026 Université de Lorraine
+ * Copyright (C) 2022 Centre National de la Recherche Scientifique
+ * Copyright (C) 2022 Université Paul Sabatier
+ *
+ * This file is part of Star-Line.
+ *
+ * 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 "test_sln_lines.h"
+
+#include "sln.h"
+
+#include <rsys/mem_allocator.h>
+#include <rsys/rsys.h>
+
+#include <star/ssp.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static struct shtr_isotope_metadata*
+setup_isotopes
+ (struct shtr* shtr,
+ struct sln_molecule molecules[SHTR_MAX_MOLECULE_COUNT])
+{
+ struct shtr_isotope_metadata* metadata = NULL;
+ FILE* fp = NULL;
+
+ CHK(fp = tmpfile());
+ fprintf(fp, "Molecule # Iso Abundance Q(296K) gj Molar Mass(g)\n");
+ write_shtr_molecule(fp, &g_H2O);
+ write_shtr_molecule(fp, &g_CO2);
+ write_shtr_molecule(fp, &g_O3);
+ rewind(fp);
+
+ CHK(molecules);
+ molecules[SHTR_H2O].concentration = 0.15;
+ molecules[SHTR_H2O].cutoff = 25; /* [cm^-1] */
+ molecules[SHTR_CO2].concentration = 0.10;
+ molecules[SHTR_CO2].cutoff = 50; /* [cm^-1] */
+ molecules[SHTR_O3].concentration = 0.05;
+ molecules[SHTR_O3].cutoff = 25; /* [cm^-1] */
+
+ CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &metadata) == RES_OK);
+
+ CHK(fclose(fp) == 0);
+
+ return metadata;
+}
+
+static struct shtr_line_list*
+setup_lines(struct shtr* shtr)
+{
+ struct shtr_line_list_load_args args = SHTR_LINE_LIST_LOAD_ARGS_NULL;
+ struct shtr_line_list* lines = NULL;
+ FILE* fp = NULL;
+
+ CHK(fp = tmpfile());
+ write_shtr_lines(fp, g_lines, g_nlines);
+ rewind(fp);
+
+ args.filename = "stream";
+ args.file = fp;
+ CHK(shtr_line_list_load(shtr, &args, &lines) == RES_OK);
+
+ CHK(fclose(fp) == 0);
+
+ return lines;
+}
+
+static void
+test_sample(struct sln_tree* tree)
+{
+ struct ssp_rng* rng = NULL;
+ const struct sln_node* node = NULL;
+ const struct sln_node* leaf = NULL;
+ double nu = 0; /* [cm^-2] */
+ double proba = 0;
+
+ CHK(ssp_rng_create(NULL, SSP_RNG_MT19937_64, &rng) == RES_OK);
+
+ CHK(node = sln_tree_get_root(tree));
+
+ /* Set an arbitrary wave number within the range of the lines */
+ nu = (g_lines[g_nlines-1].wavenumber + g_lines[0].wavenumber) / 3.0;
+
+ CHK(sln_node_sample_leaf(NULL, node, nu, rng, NULL) == NULL);
+ CHK(sln_node_sample_leaf(tree, NULL, nu, rng, NULL) == NULL);
+ CHK(sln_node_sample_leaf(tree, node, nu, NULL, NULL) == NULL);
+ CHK(sln_node_sample_leaf(tree, node, nu, rng, NULL) != NULL);
+
+ CHK(leaf = sln_node_sample_leaf(tree, node, nu, rng, &proba));
+ CHK(proba > 0 && proba < 1);
+
+ CHK(sln_node_sample_leaf(tree, leaf, nu, rng, &proba));
+ CHK(proba == 1);
+
+ /* Attempt to sample a line outside the spectral range. There are no lines
+ * with a non-zero value at the wavelength in question. In this case, the
+ * library assumes that no line can be sampled.
+ *
+ * To ensure that the nu value corresponds to a wave number whose value at the
+ * node is 0, take the last line and add 51 cm^-1 to it, which is 1 cm^-1 more
+ * than the maximum cutoff defined for the molecules in the mixture */
+ nu = g_lines[g_nlines-1].wavenumber + 51 /* [cm^-1] */;
+ CHK(sln_node_sample_leaf(tree, node, nu, rng, &proba) == NULL);
+ CHK(sln_node_sample_leaf(tree, node, INF, rng, &proba) == NULL);
+
+ CHK(ssp_rng_ref_put(rng) == RES_OK);
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(void)
+{
+ struct sln_device_create_args dev_args = SLN_DEVICE_CREATE_ARGS_DEFAULT;
+ struct sln_tree_create_args tree_args = SLN_TREE_CREATE_ARGS_DEFAULT;
+ struct sln_device* sln = NULL;
+ struct sln_tree* tree = NULL;
+
+ struct shtr_create_args shtr_args = SHTR_CREATE_ARGS_DEFAULT;
+ struct shtr* shtr = NULL;
+
+ shtr_args.verbose = 1;
+ CHK(shtr_create(&shtr_args, &shtr) == RES_OK);
+
+ dev_args.verbose = 1;
+ CHK(sln_device_create(&dev_args, &sln) == RES_OK);
+
+ tree_args.metadata = setup_isotopes(shtr, tree_args.molecules);
+ tree_args.lines = setup_lines(shtr);
+ tree_args.pressure = 10; /* [atm] */
+ tree_args.temperature = 600; /* [K] */
+ CHK(sln_tree_create(sln, &tree_args, &tree) == RES_OK);
+
+ test_sample(tree);
+
+ CHK(shtr_ref_put(shtr) == RES_OK);
+ CHK(shtr_line_list_ref_put(tree_args.lines) == RES_OK);
+ CHK(shtr_isotope_metadata_ref_put(tree_args.metadata) == RES_OK);
+
+ CHK(sln_tree_ref_put(tree) == RES_OK);
+ CHK(sln_device_ref_put(sln) == RES_OK);
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
