star-line

sln.h (17012B)

---

 /* 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/>. */
 
 #ifndef SLN_H
 #define SLN_H
 
 #include <star/shtr.h>
 #include <rsys/rsys.h>
 
 #include <float.h>
 #include <math.h>
 
 /* Library symbol management */
 #if defined(SLN_SHARED_BUILD)  /* Build shared library */
   #define SLN_API extern EXPORT_SYM
 #elif defined(SLN_STATIC)  /* Use/build static library */
   #define SLN_API extern LOCAL_SYM
 #else
   #define SLN_API extern IMPORT_SYM
 #endif
 
 /* Helper macro that asserts if the invocation of the sln function `Func'
  * returns an error. One should use this macro on sln calls for which no
  * explicit error checking is performed */
 #ifndef NDEBUG
   #define SLN(Func) ASSERT(sln_ ## Func == RES_OK)
 #else
   #define SLN(Func) sln_ ## Func
 #endif
 
 #define SLN_TREE_DEPTH_MAX 64 /* Maximum depth of a tree */
 #define SLN_TREE_ARITY_MAX 256  /* Maximum arity of a tree */
 #define SLN_LEAF_NLINES_MAX 16384 /* Maximum number of lines per leaf */
 
 /* Forward declaration of external data structures */
 struct logger;
 struct mem_allocator;
 struct shtr;
 struct shtr_line;
 struct shtr_isotope_metadata;
 struct shtr_line_list;
 
 enum sln_mesh_type {
   SLN_MESH_FIT, /* Fit the spectrum */
   SLN_MESH_UPPER, /* Upper limit of the spectrum */
   SLN_MESH_TYPES_COUNT__
 };
 
 enum sln_line_profile {
   SLN_LINE_PROFILE_VOIGT,
   SLN_LINE_PROFILES_COUNT__
 };
 
 struct sln_device_create_args {
   struct logger* logger; /* May be NULL <=> default logger */
   struct mem_allocator* allocator; /* NULL <=> use default allocator */
   int verbose; /* Verbosity level */
 };
 #define SLN_DEVICE_CREATE_ARGS_DEFAULT__ {NULL,NULL,0}
 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 {
   struct sln_isotope isotopes[SHTR_MAX_ISOTOPE_COUNT];
   double concentration;
   double cutoff; /* [cm^-1] */
   int non_default_isotope_abundances;
 };
 #define SLN_MOLECULE_NULL__ {{{0}},0,0,0}
 static const struct sln_molecule SLN_MOLECULE_NULL = SLN_MOLECULE_NULL__;
 
 struct sln_tree_create_args {
   /* Isotope metadata and list of spectral lines */
   struct shtr_isotope_metadata* metadata;
   struct shtr_line_list* lines;
 
   enum sln_line_profile line_profile;
   /* Mixture description */
   struct sln_molecule molecules[SHTR_MAX_MOLECULE_COUNT];
 
     /* Thermo dynamic properties */
   double pressure; /* [atm] */
   double temperature; /* [K] */
 
   /* Hint on the number of vertices around the line center */
   size_t nvertices_hint;
 
   /* Relative error used to simplify the spectrum mesh. The larger it is, the
    * coarser the mesh */
   double mesh_decimation_err; /* > 0 */
   enum sln_mesh_type mesh_type; /* Type of mesh to generate */
 
   /* Maximum number of children per node */
   unsigned arity;
 
   /* Maximum number of lines per leaf */
   unsigned leaf_nlines;
 
   /* When this option is enabled, the polylines of internal nodes are
    * constructed by merging their children's polylines in pairs (and then
    * simplifying the result), and repeating the process until only a single
    * polyline remains, which becomes the internal node's polyline.
    *
    * If this option is disabled, all child polylines are merged in a single step
    * before being simplified.
    *
    * Enabling this option only makes sense for trees with an arity greater than
    * two. For a binary tree, both methods should produce exactly the same tree,
    * down to the bit */
   int collapse_polylines;
 
   /* Advice on the number of threads to use */
   unsigned nthreads_hint;
 };
 #define SLN_TREE_CREATE_ARGS_DEFAULT__ { \
   NULL, /* metadata */ \
   NULL, /* line list */ \
   SLN_LINE_PROFILE_VOIGT, /* Profile */ \
   {SLN_MOLECULE_NULL__}, /* Molecules */ \
   0, /* Pressure [atm] */ \
   0, /* Temperature [K] */ \
   16, /* #vertices hint */ \
   0.01f, /* Mesh decimation error */ \
   SLN_MESH_UPPER, /* Mesh type */ \
   2, /* Arity */ \
   1, /* Number of lines per leaf */ \
   0, /* Collapse polylines */ \
   (unsigned)(-1), /* #threads hint */ \
 }
 static const struct sln_tree_create_args SLN_TREE_CREATE_ARGS_DEFAULT =
   SLN_TREE_CREATE_ARGS_DEFAULT__;
 
 struct sln_tree_read_args {
   /* Metadata and list of spectral lines from which the tree was constructed */
   struct shtr_isotope_metadata* metadata;
   struct shtr_line_list* lines;
 
   /* Name of the file to read or of the provided stream.
    * NULL <=> uses a default name for the stream to be read, which must
    * therefore be defined. */
   const char* filename; /* Name of the file to read */
   FILE* file; /* Stream from where data are read. NULL <=> read from file */
 
   /* Verify that the digital signature of the input lines matches the one stored
    * in the tree. In other words, ensure that this list of lines is indeed the
    * one used to construct the tree. An error is returned if the signatures do
    * not match.
    *
    * Although it is always advisable to verify that the data matches what is
    * expected, calculating the signatures of the lines can be time-consuming.
    * Therefore, a user who is _certain_ that the data matches can disable this
    * verification */
   int disable_line_hash_check;
 };
 #define SLN_TREE_READ_ARGS_NULL__ {NULL,NULL,NULL,NULL,0}
 static const struct sln_tree_read_args SLN_TREE_READ_ARGS_NULL =
   SLN_TREE_READ_ARGS_NULL__;
 
 struct sln_tree_write_args {
   /* Name of the file in which the tree is serialized.
    * NULL <=> uses a default name for the stream to be written, which must
    * therefore be defined. */
   const char* filename; /* Name of the file to read */
 
   /* Stream where data is written.
    * NULL <=> write to the file defined by "filename" */
   FILE* file;
 };
 #define SLN_TREE_WRITE_ARGS_NULL__ {NULL,NULL}
 static const struct sln_tree_write_args SLN_TREE_WRITE_ARGS_NULL =
   SLN_TREE_WRITE_ARGS_NULL__;
 
 struct sln_tree_desc {
   double mesh_decimation_err;
   enum sln_mesh_type mesh_type;
   enum sln_line_profile line_profile;
 
   double pressure; /* [atm] */
   double temperature; /* [K] */
 
   unsigned depth; /* #edges from the root to the deepest leaf */
   size_t nlines;
   size_t nvertices;
   size_t nnodes;
   unsigned arity;
   unsigned leaf_nlines;
 };
 #define SLN_TREE_DESC_NULL__ { \
   0,SLN_MESH_TYPES_COUNT__,SLN_LINE_PROFILES_COUNT__,0,0,0,0,0,0,0,0 \
 }
 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];
   size_t nvertices;
   unsigned nchildren;
 };
 #define SLN_NODE_DESC_NULL__ {{0,0},0,0}
 static const struct sln_node_desc SLN_NODE_DESC_NULL = SLN_NODE_DESC_NULL__;
 
 struct sln_vertex { /* 8 Bytes */
   float wavenumber; /* in cm^-1 */
   float ka;
 };
 #define SLN_VERTEX_NULL__ {0,0}
 static const struct sln_vertex SLN_VERTEX_NULL = SLN_VERTEX_NULL__;
 
 struct sln_mesh {
   const struct sln_vertex* vertices;
   size_t nvertices;
 };
 #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__;
 
 struct sln_line {
   double wavenumber; /* Line center wrt pressure in cm^-1 */
   double profile_factor; /* m^-1.cm^-1 (1e2*density*intensity) */
   double gamma_d; /* Doppler half width */
   double gamma_l; /* Lorentz half width */
   enum shtr_molecule_id molecule_id;
 };
 #define SLN_LINE_NULL__ {0,0,0,0,SHTR_MOLECULE_ID_NULL}
 static const struct sln_line SLN_LINE_NULL = SLN_LINE_NULL__;
 
 /* External data structure */
 struct ssp_rng;
 
 /* Forward declarations of opaque data structures */
 struct sln_device;
 struct sln_mixture;
 struct sln_node;
 struct sln_tree;
 
 BEGIN_DECLS
 
 /*******************************************************************************
  * Device API
  ******************************************************************************/
 SLN_API res_T
 sln_device_create
   (const struct sln_device_create_args* args,
    struct sln_device** sln);
 
 SLN_API res_T
 sln_device_ref_get
   (struct sln_device* sln);
 
 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 int
 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 int index);
 
 SLN_API res_T
 sln_mixture_get_molecule
   (const struct sln_mixture* mixture,
    const int index,
    struct sln_molecule* molecule);
 
 /*******************************************************************************
  * Tree API
  ******************************************************************************/
 SLN_API res_T
 sln_tree_create
   (struct sln_device* dev,
    const struct sln_tree_create_args* args,
    struct sln_tree** tree);
 
 /* Read a tree serialized with the "sln_tree_write" function */
 SLN_API res_T
 sln_tree_read
   (struct sln_device* sln,
    const struct sln_tree_read_args* args,
    struct sln_tree** tree);
 
 SLN_API res_T
 sln_tree_ref_get
   (struct sln_tree* tree);
 
 SLN_API res_T
 sln_tree_ref_put
   (struct sln_tree* tree);
 
 SLN_API res_T
 sln_tree_get_desc
   (const struct sln_tree* tree,
    struct sln_tree_desc* desc);
 
 SLN_API const struct sln_node* /* NULL <=> No node */
 sln_tree_get_root
   (const struct sln_tree* tree);
 
 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
 sln_tree_write
   (const struct sln_tree* tree,
    const struct sln_tree_write_args* args);
 
 /*******************************************************************************
  * Node API
  ******************************************************************************/
 SLN_API int
 sln_node_is_leaf
   (const struct sln_node* node);
 
 SLN_API unsigned
 sln_node_get_child_count
   (const struct sln_tree* tree,
    const struct sln_node* node);
 
 /* The node must not be a leaf */
 SLN_API const struct sln_node*
 sln_node_get_child
   (const struct sln_tree* tree,
    const struct sln_node* node,
    const unsigned ichild); /* 0 or #children */
 
 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
 sln_node_get_desc
   (const struct sln_tree* tree,
    const struct sln_node* node,
    struct sln_node_desc* desc);
 
 SLN_API res_T
 sln_node_get_mesh
   (const struct sln_tree* tree,
    const struct sln_node* node,
    struct sln_mesh* mesh);
 
 /* Sample a leaf based on its importance, that is, its contribution to the
  * node's absorption spectrum at the given wave number */
 SLN_API const struct sln_node* /* NULL <=> an error occurs */
 sln_node_sample_leaf
   (const struct sln_tree* tree,
    const struct sln_node* node,
    const double nu, /* [cm^-1] */
    struct ssp_rng* rng,
    double* proba); /* May be NULL */
 
 /*******************************************************************************
  * Miscellaneous
  ******************************************************************************/
 SLN_API double
 sln_line_eval
   (const struct sln_tree* tree,
    const struct sln_line* line,
    const double wavenumber); /* In cm^-1 */
 
 SLN_API double
 sln_mesh_eval
   (const struct sln_mesh* mesh,
    const double wavenumber); /* In cm^-1 */
 
 /*******************************************************************************
  * Helper functions
  ******************************************************************************/
 /* Purpose: to calculate the Faddeeva function with relative error less than
  * 10^(-4).
  *
  * Inputs: x and y, parameters for the Voigt function :
  * - x is defined as x=(nu-nu_c)/gamma_D*sqrt(ln(2)) with nu the current
  *   wavenumber, nu_c the wavenumber at line center, gamma_D the Doppler
  *   linewidth.
  * - y is defined as y=gamma_L/gamma_D*sqrt(ln(2)) with gamma_L the Lorentz
  *   linewith and gamma_D the Doppler linewidth
  *
  * Output: k, the Voigt function; it has to be multiplied by
  * sqrt(ln(2)/pi)*1/gamma_D so that the result may be interpretable in terms of
  * line profile.
  *
  * TODO check the copyright */
 SLN_API double
 sln_faddeeva
   (const double x,
    const double y);
 
 static INLINE double
 sln_compute_line_half_width_doppler
   (const double nu, /* Line center wrt pressure in cm^-1 */ /* TODO check this */
    const double molar_mass, /* In kg.mol^-1 */
    const double temperature) /* In K */
 {
   /* kb = 1.3806e-23
    * Na = 6.02214076e23
    * c = 299792458
    * sqrt(2*log(2)*kb*Na)/c */
   const double sqrt_two_ln2_kb_Na_over_c = 1.1324431552553545042e-08;
   const double gamma_d = nu * sqrt_two_ln2_kb_Na_over_c * sqrt(temperature/molar_mass);
   ASSERT(temperature >= 0 && molar_mass > 0);
   return gamma_d;
 }
 
 static INLINE double
 sln_compute_line_half_width_lorentz
   (const double gamma_air, /* Air broadening half width [cm^-1.atm^-1] */
    const double gamma_self, /* Air broadening half width [cm^-1.atm^-1] */
    const double temperature, /* [K] */
    const double pressure, /* [atm^-1] */
    const double n_air,
    const double concentration)
 {
   const double TREF=296; /* Ref temperature [K] for HITRAN/HITEMP database */
   const double Ps = pressure * concentration;
   const double n_self = n_air; /* In HITRAN n_air == n_self */
   const double gamma_l =
     pow(TREF/temperature, n_air)  * (pressure - Ps) * gamma_air
   + pow(TREF/temperature, n_self) * Ps * gamma_self;
   ASSERT(gamma_air > 0 && gamma_self > 0);
   ASSERT(pressure > 0 && concentration >= 0 && concentration <= 1);
 
   return gamma_l;
 }
 
 static INLINE double
 sln_compute_voigt_profile
   (const double wavenumber, /* In cm^-1 */
    const double nu, /* Line center in cm^-1 */
    const double gamma_d, /* Doppler line half width in cm^-1 */
    const double gamma_l) /* Lorentz line half width in cm^-1 */
 {
   /* Constants */
   const double sqrt_ln2 = 0.83255461115769768821; /* sqrt(log(2)) */
   const double sqrt_ln2_over_pi = 0.46971863934982566180; /* sqrt(log(2)/M_PI) */
   const double sqrt_ln2_over_gamma_d = sqrt_ln2 / gamma_d;
 
   const double x = (wavenumber - nu) * sqrt_ln2_over_gamma_d;
   const double y = gamma_l * sqrt_ln2_over_gamma_d;
   const double k = sln_faddeeva(x, y);
   return k*sqrt_ln2_over_pi/gamma_d;
 }
 
 static INLINE const char*
 sln_mesh_type_cstr(const enum sln_mesh_type type)
 {
   const char* cstr = NULL;
 
   switch(type) {
     case SLN_MESH_FIT:   cstr = "fit"; break;
     case SLN_MESH_UPPER: cstr = "upper"; break;
     default: FATAL("Unreachable code\n"); break;
   }
   return cstr;
 }
 
 static INLINE const char*
 sln_line_profile_cstr(const enum sln_line_profile profile)
 {
   const char* cstr = NULL;
 
   switch(profile) {
     case SLN_LINE_PROFILE_VOIGT: cstr = "voigt"; break;
     default: FATAL("Unreachable code\n"); break;
   }
   return cstr;
 }
 
 END_DECLS
 
 #endif /* SLN_H */