commit 5c7df6fda666cc3a086cc0f23bbaa24a7f296e6f
parent bdb9e3b90ad7b1504bf5004122dbeeb54abeb9f0
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 4 Feb 2026 12:46:12 +0100
Add the mixture API
It provides functions for loading and accessing additional data about
the mixture, such as the concentration of molecules, the threshold of
their absorption lines, and, possibly, the abundance of their isotopes
if their default values are not suitable.
Only a few notes on the loaded file format have been written so far, in
other words, its man page still needs to be written. Furthermore, no
tests have been carried out yet.
Attentive readers may think that this is the reintroduction of a
previous API, but this is not the case. Previously, a mixture defined
all thermodynamic properties, molecular concentrations, isotope
abundances, etc., while selecting a subset of molecules and lines to be
taken into account. Now, a mixture only defines a restricted list of
parameters. Thus, the same mixture can now be used on different lines,
for different thermodynamic properties, etc. In short, the different
concepts are now better separated.
Diffstat:
6 files changed, 551 insertions(+), 6 deletions(-)
diff --git a/Makefile b/Makefile
@@ -35,6 +35,7 @@ SRC =\
src/sln_device.c\
src/sln_faddeeva.c\
src/sln_line.c\
+ src/sln_mixture.c\
src/sln_polyline.c\
src/sln_tree_build.c\
src/sln_tree.c
diff --git a/doc/sln-mixture.txt b/doc/sln-mixture.txt
@@ -0,0 +1,12 @@
+<mixture> ::= <molecule> <concentration> <cutoff>
+ [ <isotope-list> ]
+<isotope-list> ::= <isotope-id> <abundance>
+ ...
+
+EXAMPLE
+
+ # Molecule concentration [mol/mol] cutoff [cm^-1
+ H2O 0.3 25
+ # Isotope Abundance
+ 161 9.97E-01
+ ...
diff --git a/src/sln.h b/src/sln.h
@@ -73,14 +73,18 @@ struct sln_device_create_args {
static const struct sln_device_create_args SLN_DEVICE_CREATE_ARGS_DEFAULT =
SLN_DEVICE_CREATE_ARGS_DEFAULT__;
+struct sln_isotope {
+ double abundance; /* in [0, 1] */
+ int id; /* Identifier of the isotope */
+};
struct sln_molecule {
- double isotope_abundances[SLN_MAX_ISOTOPES_COUNT];
+ struct sln_isotope isotopes[SLN_MAX_ISOTOPES_COUNT];
double concentration;
double cutoff; /* [cm^-1] */
int non_default_isotope_abundances;
};
-#define SLN_MOLECULE_NULL__ {{0},0,0,0}
+#define SLN_MOLECULE_NULL__ {{{0}},0,0,0}
static const struct sln_molecule SLN_MOLECULE_NULL = SLN_MOLECULE_NULL__;
struct sln_tree_create_args {
@@ -143,11 +147,25 @@ struct sln_mesh {
#define SLN_MESH_NULL__ {NULL,0}
static const struct sln_mesh SLN_MESH_NULL = SLN_MESH_NULL__;
+struct sln_mixture_load_args {
+ const char* filename; /* Name of the file to load or of the provided stream */
+ FILE* file; /* Stream from where data are loaded. NULL <=> load from file */
+
+ /* Metadata from which the mix is defined */
+ struct shtr_isotope_metadata* molparam;
+};
+#define SLN_MIXTURE_LOAD_ARGS_NULL__ {NULL,NULL,NULL}
+static const struct sln_mixture_load_args SLN_MIXTURE_LOAD_ARGS_NULL =
+ SLN_MIXTURE_LOAD_ARGS_NULL__;
+
/* Forward declarations of opaque data structures */
struct sln_device;
+struct sln_mixture;
struct sln_node;
struct sln_tree;
+BEGIN_DECLS
+
/*******************************************************************************
* Device API
******************************************************************************/
@@ -164,6 +182,39 @@ SLN_API res_T
sln_device_ref_put
(struct sln_device* sln);
+
+/*******************************************************************************
+ * Mixture API
+ ******************************************************************************/
+SLN_API res_T
+sln_mixture_load
+ (struct sln_device* dev,
+ const struct sln_mixture_load_args* args,
+ struct sln_mixture** mixture);
+
+SLN_API res_T
+sln_mixture_ref_get
+ (struct sln_mixture* mixture);
+
+SLN_API res_T
+sln_mixture_ref_put
+ (struct sln_mixture* mixture);
+
+SLN_API unsigned
+sln_mixture_get_molecule_count
+ (const struct sln_mixture* mixture);
+
+SLN_API enum shtr_molecule_id
+sln_mixture_get_molecule_id
+ (const struct sln_mixture* mixture,
+ const unsigned index);
+
+SLN_API res_T
+sln_mixture_get_molecule
+ (const struct sln_mixture* mixture,
+ const unsigned index,
+ struct sln_molecule* molecule);
+
/*******************************************************************************
* Tree API
******************************************************************************/
@@ -312,4 +363,6 @@ sln_compute_voigt_profile
return k*sqrt_ln2_over_pi/gamma_d;
}
+END_DECLS
+
#endif /* SLN_H */
diff --git a/src/sln_line.c b/src/sln_line.c
@@ -131,7 +131,7 @@ line_profile_factor
if(!mol_params->non_default_isotope_abundances) { /* Use default abundance */
intensity_ref = shtr_line->intensity;
} else {
- iso_abundance = mol_params->isotope_abundances[shtr_line->isotope_id_local];
+ iso_abundance = mol_params->isotopes[shtr_line->isotope_id_local].abundance;
intensity_ref = shtr_line->intensity/isotope->abundance*iso_abundance;
}
intensity = line_intensity(intensity_ref, shtr_line->lower_state_energy, Q,
diff --git a/src/sln_mixture.c b/src/sln_mixture.c
@@ -0,0 +1,479 @@
+/* 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 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 200809L /* strtok_r support */
+
+#include "sln.h"
+#include "sln_device_c.h"
+
+#include <rsys/cstr.h>
+#include <rsys/text_reader.h>
+
+#include <ctype.h> /* isalpha support */
+#include <string.h>
+
+struct molecule {
+ struct sln_molecule param;
+ int nisotopes;
+ enum shtr_molecule_id id;
+};
+static const struct molecule MOLECULE_NULL = {
+ SLN_MOLECULE_NULL, 0, SHTR_MOLECULE_ID_NULL
+};
+
+#define MOLECULE_IS_VALID(Molecule) ((Molecule)->id >= SHTR_MOLECULE_ID_NULL)
+
+struct sln_mixture {
+ struct molecule molecules[SHTR_MAX_MOLECULE_COUNT];
+ int nmolecules;
+ struct sln_device* sln;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE res_T
+check_sln_mixture_load_args(const struct sln_mixture_load_args* args)
+{
+ if(!args || args->filename) return RES_BAD_ARG;
+ if(!args->molparam) return RES_BAD_ARG;
+ return RES_OK;
+}
+
+static INLINE res_T
+flush_molecule
+ (struct sln_mixture* mixture,
+ struct molecule* molecule,
+ const char* mixture_name)
+{
+ res_T res = RES_OK;
+ ASSERT(mixture && molecule);
+
+ mixture->molecules[mixture->nmolecules++] = *molecule;
+ *molecule = MOLECULE_NULL;
+
+ /* Check isotope abundances */
+ if(molecule->param.non_default_isotope_abundances) {
+ double sum = 0;
+ int i = 0;
+
+ FOR_EACH(i, 0, molecule->nisotopes) {
+ sum += molecule->param.isotopes[i].abundance;
+ }
+
+ if(!eq_eps(sum, 1, 1e-6)) {
+ ERROR(mixture->sln,
+ "%s: the sum of the isotopic abundances of %s is %g, "
+ "whereas it should be equal to 1\n",
+ mixture_name, shtr_molecule_cstr(molecule->id), sum);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static INLINE res_T
+parse_molecule_name
+ (struct sln_mixture* mixture,
+ struct molecule* molecule,
+ const char* name)
+{
+ size_t i = 0;
+ ASSERT(mixture && molecule && name);
+ (void)mixture;
+
+ /* Search for the molecule ID. Use a simple linear search, as there are only a
+ * few molecules. */
+ FOR_EACH(i, 1, SHTR_MAX_MOLECULE_COUNT) {
+ if(!strcmp(name, shtr_molecule_cstr(i))) break;
+ }
+ if(i >= SHTR_MAX_MOLECULE_COUNT) return RES_BAD_ARG;
+ molecule->id = i;
+
+ return RES_OK;
+}
+
+static INLINE int
+molecule_already_defined
+ (struct sln_mixture* mixture,
+ const enum shtr_molecule_id id)
+{
+ int i = 0;
+ ASSERT(mixture);
+
+ /* Check that this molecule has not already been parsed. Use a simple linear
+ * search for the molecule, as there are only a few of them. */
+ FOR_EACH(i, 0, mixture->nmolecules) {
+ if(mixture->molecules[i].id == id) return 1;
+ }
+ return 0;
+}
+
+static res_T
+parse_molecule
+ (struct sln_mixture* mixture,
+ const struct sln_mixture_load_args* args,
+ struct molecule* molecule,
+ struct txtrdr* txtrdr)
+{
+ char* line = NULL;
+ char* tk = NULL;
+ char* tk_ctx = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(mixture && args && molecule && txtrdr);
+ (void)args;
+
+ line = txtrdr_get_line(txtrdr);
+ ASSERT(line);
+
+ #define LOG(Type, Str, ...) \
+ Type(mixture->sln, "%s:%lu "Str, \
+ txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr), \
+ __VA_ARGS__)
+
+ tk = strtok_r(line, " \t", &tk_ctx);
+ res = parse_molecule_name(mixture, molecule, tk);
+ if(res != RES_OK) {
+ LOG(ERROR, "invalid molecule name `%s'\n", tk);
+ goto error;
+ }
+ if(molecule_already_defined(mixture, molecule->id)) {
+ LOG(ERROR, "duplicate molecule`%s'\n", tk);
+ goto error;
+ }
+
+ tk = strtok_r(NULL, " \t", &tk_ctx);
+ res = cstr_to_double(tk, &molecule->param.concentration);
+ if(res == RES_OK && molecule->param.concentration < 0) res = RES_BAD_ARG;
+ if(res != RES_OK) {
+ LOG(ERROR, "invalid concentration `%s'\n", tk);
+ goto error;
+ }
+
+ tk = strtok_r(NULL, " \t", &tk_ctx);
+ res = cstr_to_double(tk, &molecule->param.cutoff);
+ if(res != RES_OK && molecule->param.cutoff <= 0) res = RES_BAD_ARG;
+ if(res != RES_OK) {
+ LOG(ERROR, "invalid cutoff `%s'\n", tk);
+ goto error;
+ }
+
+ tk = strtok_r(NULL, "", &tk_ctx);
+ if(tk) LOG(WARN, "unexpected text `%s'\n", tk);
+
+ #undef LOG
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/* Initialize the isotope identifiers of molecules based on those defined in the
+ * isotope metadata */
+static res_T
+setup_molecule_isotope_ids
+ (struct shtr_isotope_metadata* molparam,
+ struct molecule* molecule)
+{
+ struct shtr_molecule shtr_molecule = SHTR_MOLECULE_NULL;
+ size_t i = 0;
+ res_T res = RES_OK;
+
+ ASSERT(molparam && molecule);
+
+ /* The molecule must be defined in the isotope metadata */
+ res = shtr_isotope_metadata_find_molecule
+ (molparam, molecule->id, &shtr_molecule);
+ if(res != RES_OK || SHTR_MOLECULE_IS_NULL(&shtr_molecule)) {
+ goto error;
+ }
+
+ FOR_EACH(i, 0, shtr_molecule.nisotopes) {
+ molecule->param.isotopes[i].id = shtr_molecule.isotopes[i].id;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_isotope
+ (struct sln_mixture* mixture,
+ const struct sln_mixture_load_args* args,
+ struct molecule* molecule,
+ struct txtrdr* txtrdr)
+{
+ char* line = NULL;
+ char* tk = NULL;
+ char* tk_ctx = NULL;
+ int iiso = 0;
+ int iso_id = 0;
+ res_T res = RES_OK;
+ ASSERT(mixture && args && molecule && txtrdr);
+
+ line = txtrdr_get_line(txtrdr);
+ ASSERT(line);
+
+ if(!molecule->param.non_default_isotope_abundances) {
+ /* No isotopes should have been parsed */
+ ASSERT(molecule->nisotopes == 0);
+
+ res = setup_molecule_isotope_ids(args->molparam, molecule);
+ if(res != RES_OK) goto error;
+
+ molecule->param.non_default_isotope_abundances = 1;
+ }
+
+ iiso = molecule->nisotopes; /* Local index of the isotope */
+
+ #define LOG(Type, Str, ...) \
+ Type(mixture->sln, "%s:%lu "Str, \
+ txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr), \
+ __VA_ARGS__)
+
+ tk = strtok_r(line, " \t", &tk_ctx);
+ res = cstr_to_int(tk, &iso_id);
+ if(res != RES_OK) {
+ LOG(ERROR, "invalid isotope index `%s'\n", tk);
+ goto error;
+ }
+ if(iso_id != molecule->param.isotopes[iiso].id) {
+ LOG(ERROR, "expecting the %s isotope %d\n",
+ shtr_molecule_cstr(molecule->id), iso_id);
+ goto error;
+ }
+
+ tk = strtok_r(NULL, " \t", &tk_ctx);
+ res = cstr_to_double(tk, &molecule->param.isotopes[iiso].abundance);
+ if(res != RES_OK || molecule->param.isotopes[iiso].abundance < 0) {
+ LOG(ERROR, "invalid isotope abundance `%s'\n", tk);
+ goto error;
+ }
+
+ tk = strtok_r(NULL, "", &tk_ctx);
+ if(tk) LOG(WARN, "unexpected text `%s'\n", tk);
+
+ #undef LOG
+
+ molecule->param.non_default_isotope_abundances = 1;
+ molecule->nisotopes += 1;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_line
+ (struct sln_mixture* mixture,
+ const struct sln_mixture_load_args* args,
+ struct molecule* molecule, /* Currently parsed molecule */
+ struct txtrdr* txtrdr)
+{
+ const char* line = NULL;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(mixture && args && molecule && txtrdr);
+
+ line = txtrdr_get_cline(txtrdr);
+ ASSERT(line);
+ i = strspn(line, " \t");
+ ASSERT(i < strlen(line));
+
+ /* New molecule */
+ if(isalpha(line[i])) {
+
+ /* Register the previous molecule if any */
+ if(MOLECULE_IS_VALID(molecule)) {
+ res = flush_molecule(mixture, molecule, txtrdr_get_name(txtrdr));
+ if(res != RES_OK) goto error;
+ }
+
+ /* Parse the molecule data, i.e., name, concentration and cutoff */
+ res = parse_molecule(mixture, args, molecule, txtrdr);
+ if(res != RES_OK) goto error;
+
+ /* If there is no molecule being analyzed, no isotopic data can be parsed */
+ } else if(!MOLECULE_IS_VALID(molecule)) {
+ ERROR(mixture->sln, "%s:%lu: missing a molecule.\n",
+ txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr));
+ res = RES_BAD_ARG;
+ goto error;
+
+ /* Parse the isotopes for the currently parsed molecule */
+ } else {
+ res = parse_isotope(mixture, args, molecule, txtrdr);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+load_stream
+ (struct sln_mixture* mixture,
+ FILE* fp,
+ const struct sln_mixture_load_args* args)
+{
+ struct molecule molecule = MOLECULE_NULL;
+ struct txtrdr* txtrdr = NULL;
+ res_T res = RES_OK;
+ ASSERT(mixture && fp && args);
+
+ res = txtrdr_stream(mixture->sln->allocator, fp, args->filename, '#', &txtrdr);
+ if(res != RES_OK) goto error;
+
+ #define READ_LINE if((res = txtrdr_read_line(txtrdr)) != RES_OK) goto error
+
+ for(;;) {
+ READ_LINE;
+ if(!txtrdr_get_cline(txtrdr)) break; /* No more parsed line */
+
+ res = parse_line(mixture, args, &molecule, txtrdr);
+ if(res != RES_OK) goto error;
+ }
+ #undef READ_LINE
+
+ if(MOLECULE_IS_VALID(&molecule)) {
+ res = flush_molecule(mixture, &molecule, txtrdr_get_name(txtrdr));
+ if(res != RES_OK) goto error;
+ }
+exit:
+ if(txtrdr) txtrdr_ref_put(txtrdr);
+ return res;
+error:
+ if(!txtrdr) {
+ ERROR(mixture->sln, "%s: loading error -- %s\n",
+ args->filename, res_to_cstr(res));
+ } else {
+ ERROR(mixture->sln, "%s:%lu: loading error -- %s\n",
+ args->filename, (unsigned long)txtrdr_get_line_num(txtrdr),
+ res_to_cstr(res));
+ }
+ goto exit;
+}
+
+static res_T
+load_mixture
+ (struct sln_mixture* mixture,
+ const struct sln_mixture_load_args* args)
+{
+ FILE* fp = NULL;
+ res_T res = RES_OK;
+ ASSERT(mixture && args);
+
+ if(args->file) { /* Load from stream */
+ fp = args->file;
+
+ } else { /* Load from file */
+ fp = fopen(args->filename, "r");
+ if(!fp) {
+ ERROR(mixture->sln, "error opening file `%s' -- %s\n",
+ args->filename, strerror(errno));
+ res = RES_IO_ERR;
+ goto error;
+ }
+ }
+
+ res = load_stream(mixture, fp, args);
+ if(res != RES_OK) goto error;
+
+exit:
+ if(!fp && fp != args->file) fclose(args->file);
+ return res;
+error:
+ goto exit;
+}
+
+static void
+release_mixture(ref_T* ref)
+{
+ struct sln_mixture* mixture = CONTAINER_OF(ref, struct sln_mixture, ref);
+ struct sln_device* sln = NULL;
+ ASSERT(ref);
+ sln = mixture->sln;
+ MEM_RM(sln->allocator, mixture);
+ SLN(device_ref_put(sln));
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+sln_mixture_load
+ (struct sln_device* sln,
+ const struct sln_mixture_load_args* args,
+ struct sln_mixture** out_mixture)
+{
+ struct sln_mixture* mixture = NULL;
+ res_T res = RES_OK;
+
+ if(!sln || !out_mixture) { res = RES_BAD_ARG; goto error; }
+ res = check_sln_mixture_load_args(args);
+ if(res != RES_OK) goto error;
+
+ mixture = MEM_CALLOC(sln->allocator, 1, sizeof(*mixture));
+ if(!mixture) {
+ ERROR(sln, "Could not allocate the mixture data structure.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&mixture->ref);
+ SLN(device_ref_get(sln));
+ mixture->sln = sln;
+
+ res = load_mixture(mixture, args);
+ if(res != RES_OK) goto error;
+
+exit:
+ if(out_mixture) *out_mixture = mixture;
+ return res;
+error:
+ if(mixture) { SLN(mixture_ref_put(mixture)); mixture = NULL; }
+ goto exit;
+}
+
+res_T
+sln_mixture_ref_get(struct sln_mixture* mixture)
+{
+ if(!mixture) return RES_BAD_ARG;
+ ref_get(&mixture->ref);
+ return RES_OK;
+}
+
+res_T
+sln_mixture_ref_put(struct sln_mixture* mixture)
+{
+ if(!mixture) return RES_BAD_ARG;
+ ref_put(&mixture->ref, release_mixture);
+ return RES_OK;
+}
diff --git a/src/sln_tree.c b/src/sln_tree.c
@@ -73,15 +73,15 @@ check_molecule_isotope_abundances
/* The isotopic abundances are not the default ones.
* Verify that they are valid ... */
FOR_EACH(i, 0, SLN_MAX_ISOTOPES_COUNT) {
- if(molecule->isotope_abundances[i] < 0) {
+ if(molecule->isotopes[i].abundance < 0) {
const int isotope_id = i + 1; /* isotope id in [1, 10] */
ERROR(sln, "%s: invalid abundance of isotopie %d of %s: %g.\n",
caller, isotope_id, shtr_molecule_cstr(i),
- molecule->isotope_abundances[i]);
+ molecule->isotopes[i].abundance);
return RES_BAD_ARG;
}
- sum += molecule->isotope_abundances[i];
+ sum += molecule->isotopes[i].abundance;
}
/* ... and that their sum equals 1 */
