star-line

Structure for accelerating line importance sampling
git clone git://git.meso-star.fr/star-line.git
Log | Files | Refs | README | LICENSE

sln_tree.c (21897B)


      1 /* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com)
      2  * Copyright (C) 2026 Université de Lorraine
      3  * Copyright (C) 2022 Centre National de la Recherche Scientifique
      4  * Copyright (C) 2022 Université Paul Sabatier
      5  *
      6  * This program is free software: you can redistribute it and/or modify
      7  * it under the terms of the GNU General Public License as published by
      8  * the Free Software Foundation, either version 3 of the License, or
      9  * (at your option) any later version.
     10  *
     11  * This program is distributed in the hope that it will be useful,
     12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     14  * GNU General Public License for more details.
     15  *
     16  * You should have received a copy of the GNU General Public License
     17  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     18 
     19 #include "sln.h"
     20 #include "sln_device_c.h"
     21 #include "sln_line.h"
     22 #include "sln_tree_c.h"
     23 
     24 #include <star/shtr.h>
     25 
     26 #include <rsys/algorithm.h>
     27 #include <rsys/cstr.h>
     28 #include <rsys/math.h>
     29 
     30 struct stream {
     31   const char* name;
     32   FILE* fp;
     33   int intern_fp; /* Define if the stream was internally opened */
     34 };
     35 static const struct stream STREAM_NULL = {NULL, NULL, 0};
     36 
     37 /*******************************************************************************
     38  * Helper functions
     39  ******************************************************************************/
     40 /* Check that the sum of the molecular concentrations is equal to 1 */
     41 static res_T
     42 check_molecule_concentration
     43   (struct sln_device* sln,
     44    const char* caller,
     45    const struct sln_tree_create_args* args)
     46 {
     47   int i = 0;
     48   double sum = 0;
     49   ASSERT(sln && caller && args);
     50 
     51   FOR_EACH(i, 0, SHTR_MAX_MOLECULE_COUNT) {
     52     if(i == SHTR_MOLECULE_ID_NULL) continue;
     53     sum += args->molecules[i].concentration;
     54   }
     55 
     56   /* The sum of molecular concentrations must be less than or equal to 1. It may
     57    * be less than 1 if the remaining part of the mixture is (implicitly) defined
     58    * as a radiatively inactive gas */
     59   if(sum > 1 && sum-1 > 1e-6) {
     60     ERROR(sln,
     61       "%s: the sum of molecule concentrations is greater than 1: %g\n",
     62       caller, sum);
     63     return RES_BAD_ARG;
     64   }
     65 
     66   return RES_OK;
     67 }
     68 
     69 /* Verify that the isotope abundance are valids */
     70 static res_T
     71 check_molecule_isotope_abundances
     72   (struct sln_device* sln,
     73    const char* caller,
     74    const struct sln_molecule* molecule)
     75 {
     76   int i = 0;
     77   double sum = 0;
     78   ASSERT(sln && caller && molecule);
     79 
     80   /* The isotopic abundances are the default ones. Nothing to do */
     81   if(!molecule->non_default_isotope_abundances) return RES_OK;
     82 
     83   /* The isotopic abundances are not the default ones.
     84    * Verify that they are valid ... */
     85   FOR_EACH(i, 0, SHTR_MAX_ISOTOPE_COUNT) {
     86     if(molecule->isotopes[i].abundance < 0) {
     87       const int isotope_id = i + 1; /* isotope id in [1, 10] */
     88       ERROR(sln, "%s: invalid abundance of isotopie %d of %s: %g.\n",
     89         caller, isotope_id, shtr_molecule_cstr(i),
     90         molecule->isotopes[i].abundance);
     91       return RES_BAD_ARG;
     92     }
     93 
     94     sum += molecule->isotopes[i].abundance;
     95   }
     96 
     97   /* ... and that their sum equals 1 */
     98   if(!eq_eps(sum, 1, 1e-6)) {
     99     ERROR(sln, "%s: the %s isotope abundances does not sum to 1: %g.\n",
    100       caller, shtr_molecule_cstr(i), sum);
    101     return RES_BAD_ARG;
    102   }
    103 
    104   return RES_OK;
    105 }
    106 
    107 static res_T
    108 check_molecules
    109   (struct sln_device* sln,
    110    const char* caller,
    111    const struct sln_tree_create_args* args)
    112 {
    113   char molecule_ok[SHTR_MAX_MOLECULE_COUNT] = {0};
    114 
    115   double concentrations_sum = 0;
    116   size_t iline = 0;
    117   size_t nlines = 0;
    118   res_T res = RES_OK;
    119   ASSERT(args->lines);
    120 
    121   res = check_molecule_concentration(sln, caller, args);
    122   if(res != RES_OK) return res;
    123 
    124   /* Iterate over the lines to define which molecules has to be checked, i.e.,
    125    * the ones used in the mixture */
    126   SHTR(line_list_get_size(args->lines, &nlines));
    127   FOR_EACH(iline, 0, nlines) {
    128     struct shtr_line line = SHTR_LINE_NULL;
    129     const struct sln_molecule* molecule = NULL;
    130 
    131     SHTR(line_list_at(args->lines, iline, &line));
    132 
    133     /* This molecule was already checked */
    134     if(molecule_ok[line.molecule_id]) continue;
    135 
    136     molecule = args->molecules + line.molecule_id;
    137 
    138     if(molecule->concentration == 0) {
    139       /* A molecular concentration of zero is allowed, but may be a user error,
    140        * as 0 is the default concentration in the tree creation arguments.
    141        * Therefore, warn the user about this value so that they can determine
    142        * whether or not it is an error on their part. */
    143       WARN(sln, "%s: the concentration of %s is zero.\n",
    144         caller, shtr_molecule_cstr(line.molecule_id));
    145 
    146     } else if(molecule->concentration < 0) {
    147       /* Concentration cannot be negative... */
    148       ERROR(sln, "%s: invalid %s concentration: %g.\n",
    149         FUNC_NAME, shtr_molecule_cstr(line.molecule_id),
    150         molecule->concentration);
    151       return RES_BAD_ARG;
    152     }
    153 
    154     concentrations_sum += molecule->concentration;
    155 
    156     if(molecule->cutoff <= 0) {
    157       /* ... cutoff either */
    158       ERROR(sln, "%s: invalid %s cutoff: %g.\n",
    159         caller, shtr_molecule_cstr(line.molecule_id), molecule->cutoff);
    160       return RES_BAD_ARG;
    161     }
    162 
    163     res = check_molecule_isotope_abundances(sln, caller, molecule);
    164     if(res != RES_OK) return res;
    165 
    166     molecule_ok[line.molecule_id] = 1;
    167   }
    168 
    169   /* The sum of molecular concentrations must be less than or equal to 1. It may
    170    * be less than 1 if the remaining part of the mixture is (implicitly) defined
    171    * as a radiatively inactive gas */
    172   if(concentrations_sum > 1 && (concentrations_sum - 1) > 1e-6) {
    173     ERROR(sln,
    174       "%s: the sum of molecule concentrations is greater than 1: %g\n",
    175       caller, concentrations_sum);
    176     return RES_BAD_ARG;
    177   }
    178 
    179   return RES_OK;
    180 }
    181 
    182 static res_T
    183 check_sln_tree_create_args
    184   (struct sln_device* sln,
    185    const char* caller,
    186    const struct sln_tree_create_args* args)
    187 {
    188   res_T res = RES_OK;
    189   ASSERT(sln && caller);
    190 
    191   if(!args) return RES_BAD_ARG;
    192 
    193   if(!args->metadata) {
    194     ERROR(sln, "%s: the isotope metadata are missing.\n", caller);
    195     return RES_BAD_ARG;
    196   }
    197 
    198   if(!args->lines) {
    199     ERROR(sln, "%s: the list of lines is missing.\n", caller);
    200     return RES_BAD_ARG;
    201   }
    202 
    203   if(args->nvertices_hint == 0) {
    204     ERROR(sln,
    205       "%s: invalid hint on the number of vertices around the line center %lu.\n",
    206       caller, (unsigned long)args->nvertices_hint);
    207     return RES_BAD_ARG;
    208   }
    209 
    210   if(args->mesh_decimation_err < 0) {
    211     ERROR(sln, "%s: invalid decimation error %g.\n",
    212       caller, args->mesh_decimation_err);
    213     return RES_BAD_ARG;
    214   }
    215 
    216   if((unsigned)args->mesh_type >= SLN_MESH_TYPES_COUNT__) {
    217     ERROR(sln, "%s: invalid mesh type %d.\n", caller, args->mesh_type);
    218     return RES_BAD_ARG;
    219   }
    220 
    221   if((unsigned)args->line_profile >= SLN_LINE_PROFILES_COUNT__) {
    222     ERROR(sln, "%s: invalid line profile %d.\n", caller, args->line_profile);
    223     return RES_BAD_ARG;
    224   }
    225 
    226   if(args->arity < 2 || args->arity > SLN_TREE_ARITY_MAX) {
    227     ERROR(sln, "%s: invalid arity %u. It must be in [2, %d]\n",
    228       caller, args->arity, SLN_TREE_ARITY_MAX);
    229     return RES_BAD_ARG;
    230   }
    231 
    232   res = check_molecules(sln, caller, args);
    233   if(res != RES_OK) return res;
    234 
    235   return RES_OK;
    236 }
    237 
    238 static res_T
    239 check_sln_tree_read_args
    240   (struct sln_device* sln,
    241    const char* caller,
    242    const struct sln_tree_read_args* args)
    243 {
    244   if(!args) return RES_BAD_ARG;
    245 
    246   if(!args->metadata) {
    247     ERROR(sln, "%s: the isotope metadata are missing.\n", caller);
    248     return RES_BAD_ARG;
    249   }
    250 
    251   if(!args->lines) {
    252     ERROR(sln, "%s: the list of lines is missing.\n", caller);
    253     return RES_BAD_ARG;
    254   }
    255 
    256   if(!args->file && !args->filename) {
    257     ERROR(sln,
    258       "%s: the source file is missing. No file name or stream is provided.\n",
    259       caller);
    260     return RES_BAD_ARG;
    261   }
    262 
    263   return RES_OK;
    264 }
    265 
    266 static res_T
    267 check_sln_tree_write_args
    268   (struct sln_device* sln,
    269    const char* caller,
    270    const struct sln_tree_write_args* args)
    271 {
    272   if(!args) return RES_BAD_ARG;
    273 
    274   if(!args->file && !args->filename) {
    275     ERROR(sln,
    276       "%s: the destination file is missing. "
    277       "No file name or stream is provided.\n",
    278       caller);
    279     return RES_BAD_ARG;
    280   }
    281 
    282   return RES_OK;
    283 }
    284 static INLINE void
    285 stream_release(struct stream* stream)
    286 {
    287   ASSERT(stream);
    288   if(stream->intern_fp && stream->fp) CHK(fclose(stream->fp) == 0);
    289 }
    290 
    291 static res_T
    292 stream_init
    293   (struct sln_device* sln,
    294    const char* caller,
    295    const char* name, /* NULL <=> default stream name */
    296    FILE* fp, /* NULL <=> open file "name" */
    297    const char* mode, /* mode in fopen */
    298    struct stream* stream)
    299 {
    300   res_T res = RES_OK;
    301 
    302   ASSERT(sln && caller && stream);
    303   ASSERT(fp || (name && mode));
    304 
    305   *stream = STREAM_NULL;
    306 
    307   if(fp) {
    308     stream->intern_fp = 0;
    309     stream->name = name ? name : "stream";
    310     stream->fp = fp;
    311 
    312   } else {
    313     stream->intern_fp = 1;
    314     stream->name = name;
    315     if(!(stream->fp = fopen(name, mode))) {
    316       ERROR(sln, "%s:%s: error opening file -- %s\n",
    317         caller, name, strerror(errno));
    318       res = RES_IO_ERR;
    319       goto error;
    320     }
    321   }
    322 
    323 exit:
    324   return res;
    325 error:
    326   if(stream->intern_fp && stream->fp) CHK(fclose(stream->fp) == 0);
    327   goto exit;
    328 }
    329 
    330 static res_T
    331 create_tree
    332   (struct sln_device* sln,
    333    const char* caller,
    334    struct sln_tree** out_tree)
    335 {
    336   struct sln_tree* tree = NULL;
    337   res_T res = RES_OK;
    338   ASSERT(sln && caller && out_tree);
    339 
    340   tree = MEM_CALLOC(sln->allocator, 1, sizeof(struct sln_tree));
    341   if(!tree) {
    342     ERROR(sln, "%s: could not allocate the tree data structure.\n",
    343       caller);
    344     res = RES_MEM_ERR;
    345     goto error;
    346   }
    347   ref_init(&tree->ref);
    348   SLN(device_ref_get(sln));
    349   tree->sln = sln;
    350   darray_node_init(sln->allocator, &tree->nodes);
    351   darray_vertex_init(sln->allocator, &tree->vertices);
    352 
    353 exit:
    354   *out_tree = tree;
    355   return res;
    356 error:
    357   if(tree) { SLN(tree_ref_put(tree)); tree = NULL; }
    358   goto exit;
    359 }
    360 
    361 static INLINE int
    362 cmp_nu_vtx(const void* key, const void* item)
    363 {
    364   const float nu = *((const float*)key);
    365   const struct sln_vertex* vtx = item;
    366   if(nu < vtx->wavenumber) return -1;
    367   if(nu > vtx->wavenumber) return +1;
    368   return 0;
    369 }
    370 
    371 static void
    372 release_tree(ref_T* ref)
    373 {
    374   struct sln_tree* tree = CONTAINER_OF(ref, struct sln_tree, ref);
    375   struct sln_device* sln = NULL;
    376   ASSERT(ref);
    377   sln = tree->sln;
    378   darray_node_release(&tree->nodes);
    379   darray_vertex_release(&tree->vertices);
    380   if(tree->args.lines) SHTR(line_list_ref_put(tree->args.lines));
    381   if(tree->args.metadata) SHTR(isotope_metadata_ref_put(tree->args.metadata));
    382   MEM_RM(sln->allocator, tree);
    383   SLN(device_ref_put(sln));
    384 }
    385 
    386 /*******************************************************************************
    387  * Exported symbols
    388  ******************************************************************************/
    389 res_T
    390 sln_tree_create
    391   (struct sln_device* device,
    392    const struct sln_tree_create_args* args,
    393    struct sln_tree** out_tree)
    394 {
    395   struct sln_tree* tree = NULL;
    396   res_T res = RES_OK;
    397 
    398   if(!device || !out_tree) { res = RES_BAD_ARG; goto error; }
    399   res = check_sln_tree_create_args(device, FUNC_NAME, args);
    400   if(res != RES_OK) goto error;
    401 
    402   res = create_tree(device, FUNC_NAME, &tree);
    403   if(res != RES_OK) goto error;
    404   SHTR(line_list_ref_get(args->lines));
    405   SHTR(isotope_metadata_ref_get(args->metadata));
    406   tree->args = *args;
    407 
    408   res = tree_build(tree);
    409   if(res != RES_OK) goto error;
    410 
    411 exit:
    412   if(out_tree) *out_tree = tree;
    413   return res;
    414 error:
    415   if(tree) { SLN(tree_ref_put(tree)); tree = NULL; }
    416   goto exit;
    417 }
    418 
    419 res_T
    420 sln_tree_read
    421   (struct sln_device* sln,
    422    const struct sln_tree_read_args* args,
    423    struct sln_tree** out_tree)
    424 {
    425   hash256_T hash_mdata1;
    426   hash256_T hash_mdata2;
    427   hash256_T hash_lines1;
    428   hash256_T hash_lines2;
    429 
    430   struct stream stream = STREAM_NULL;
    431   struct sln_tree* tree = NULL;
    432   size_t n = 0;
    433   int version = 0;
    434   res_T res = RES_OK;
    435 
    436   if(!sln || !out_tree) { res = RES_BAD_ARG; goto error; }
    437   res = check_sln_tree_read_args(sln, FUNC_NAME, args);
    438   if(res != RES_OK) goto error;
    439 
    440   res = create_tree(sln, FUNC_NAME, &tree);
    441   if(res != RES_OK) goto error;
    442 
    443   res = stream_init(sln, FUNC_NAME, args->filename, args->file, "r", &stream);
    444   if(res != RES_OK) goto error;
    445 
    446   #define READ(Var, Nb) {                                                      \
    447     if(fread((Var), sizeof(*(Var)), (Nb), stream.fp) != (Nb)) {                \
    448       if(feof(stream.fp)) {                                                    \
    449         res = RES_BAD_ARG;                                                     \
    450       } else if(ferror(stream.fp)) {                                           \
    451         res = RES_IO_ERR;                                                      \
    452       } else {                                                                 \
    453         res = RES_UNKNOWN_ERR;                                                 \
    454       }                                                                        \
    455       ERROR(sln, "%s: error loading the tree structure -- %s.\n",              \
    456         stream.name, res_to_cstr(res));                                        \
    457       goto error;                                                              \
    458     }                                                                          \
    459   } (void)0
    460   READ(&version, 1);
    461   if(version != SLN_TREE_VERSION) {
    462     ERROR(sln,
    463       "%s: unexpected tree version %d. Expecting a tree in version %d.\n",
    464       stream.name, version, SLN_TREE_VERSION);
    465     res = RES_BAD_ARG;
    466     goto error;
    467   }
    468 
    469   res = shtr_isotope_metadata_hash(args->metadata, hash_mdata1);
    470   if(res != RES_OK) goto error;
    471 
    472   READ(hash_mdata2, sizeof(hash256_T));
    473   if(!hash256_eq(hash_mdata1, hash_mdata2)) {
    474     ERROR(sln,
    475       "%s: the input isotopic metadata are not those used "
    476       "during tree construction.\n", stream.name);
    477     res = RES_BAD_ARG;
    478     goto error;
    479   }
    480 
    481   SHTR(isotope_metadata_ref_get(args->metadata));
    482   tree->args.metadata = args->metadata;
    483 
    484   res = shtr_line_list_hash(args->lines, hash_lines1);
    485   if(res != RES_OK) goto error;
    486 
    487   READ(hash_lines2, sizeof(hash256_T));
    488   if(!hash256_eq(hash_lines1, hash_lines2)) {
    489     ERROR(sln,
    490       "%s: the input list of lines is not the one used to build the tree.\n",
    491       stream.name);
    492     res = RES_BAD_ARG;
    493     goto error;
    494   }
    495 
    496   SHTR(line_list_ref_get(args->lines));
    497   tree->args.lines = args->lines;
    498 
    499   READ(&n, 1);
    500   if((res = darray_node_resize(&tree->nodes, n)) != RES_OK) goto error;
    501   READ(darray_node_data_get(&tree->nodes), n);
    502 
    503   READ(&n, 1);
    504   if((res = darray_vertex_resize(&tree->vertices, n)) != RES_OK) goto error;
    505   READ(darray_vertex_data_get(&tree->vertices), n);
    506 
    507   READ(&tree->args.line_profile, 1);
    508   READ(&tree->args.molecules, 1);
    509   READ(&tree->args.pressure, 1);
    510   READ(&tree->args.temperature, 1);
    511   READ(&tree->args.nvertices_hint, 1);
    512   READ(&tree->args.mesh_decimation_err, 1);
    513   READ(&tree->args.mesh_type, 1);
    514   READ(&tree->args.arity, 1);
    515   #undef READ
    516 
    517 exit:
    518   stream_release(&stream);
    519   if(out_tree) *out_tree = tree;
    520   return res;
    521 error:
    522   if(tree) { SLN(tree_ref_put(tree)); tree = NULL; }
    523   goto exit;
    524 }
    525 
    526 res_T
    527 sln_tree_ref_get(struct sln_tree* tree)
    528 {
    529   if(!tree) return RES_BAD_ARG;
    530   ref_get(&tree->ref);
    531   return RES_OK;
    532 }
    533 
    534 res_T
    535 sln_tree_ref_put(struct sln_tree* tree)
    536 {
    537   if(!tree) return RES_BAD_ARG;
    538   ref_put(&tree->ref, release_tree);
    539   return RES_OK;
    540 }
    541 
    542 res_T
    543 sln_tree_get_desc(const struct sln_tree* tree, struct sln_tree_desc* desc)
    544 {
    545   size_t nlines_adjusted = 0;
    546   unsigned depth = 0;
    547 
    548   if(!tree || !desc) return RES_BAD_ARG;
    549 
    550   desc->max_nlines_per_leaf = 1;
    551   desc->mesh_decimation_err = tree->args.mesh_decimation_err;
    552   desc->mesh_type = tree->args.mesh_type;
    553   desc->line_profile = tree->args.line_profile;
    554   desc->nnodes = darray_node_size_get(&tree->nodes);
    555   desc->nvertices = darray_vertex_size_get(&tree->vertices);
    556   desc->pressure = tree->args.pressure; /* [atm] */
    557   desc->temperature = tree->args.temperature; /* [K] */
    558   desc->arity = tree->args.arity;
    559 
    560   SHTR(line_list_get_size(tree->args.lines, &desc->nlines));
    561   nlines_adjusted = round_up_pow2(desc->nlines);
    562 
    563   for(depth=0; depth<64 && !(BIT_U64(depth) & nlines_adjusted); ++depth);
    564   desc->depth = depth;
    565 
    566 #ifndef NDEBUG
    567   {
    568     unsigned max_depth = 0;
    569     const struct sln_node* node = sln_tree_get_root(tree);
    570     while(!sln_node_is_leaf(node)) {
    571       node = sln_node_get_child(node, 0);
    572       ++max_depth;
    573     }
    574     CHK(max_depth == depth);
    575   }
    576 #endif
    577 
    578   return RES_OK;
    579 }
    580 
    581 const struct sln_node*
    582 sln_tree_get_root(const struct sln_tree* tree)
    583 {
    584   ASSERT(tree);
    585   if(darray_node_size_get(&tree->nodes)) {
    586     return darray_node_cdata_get(&tree->nodes);
    587   } else {
    588     return NULL;
    589   }
    590 }
    591 
    592 int
    593 sln_node_is_leaf(const struct sln_node* node)
    594 {
    595   ASSERT(node);
    596   return node->offset == 0;
    597 }
    598 
    599 const struct sln_node*
    600 sln_node_get_child(const struct sln_node* node, const unsigned ichild)
    601 {
    602   ASSERT(node && ichild <= 1);
    603   ASSERT(!sln_node_is_leaf(node));
    604   return node + node->offset + ichild;
    605 }
    606 
    607 size_t
    608 sln_node_get_lines_count(const struct sln_node* node)
    609 {
    610   ASSERT(node);
    611   return node->range[1] - node->range[0] + 1/*Both boundaries are inclusives*/;
    612 }
    613 
    614 res_T
    615 sln_node_get_line
    616   (const struct sln_tree* tree,
    617    const struct sln_node* node,
    618    const size_t iline,
    619    struct shtr_line* line)
    620 {
    621   if(!tree || !node || iline > sln_node_get_lines_count(node))
    622     return RES_BAD_ARG;
    623 
    624   return shtr_line_list_at
    625     (tree->args.lines, node->range[0] + iline, line);
    626 }
    627 
    628 res_T
    629 sln_node_get_mesh
    630   (const struct sln_tree* tree,
    631    const struct sln_node* node,
    632    struct sln_mesh* mesh)
    633 {
    634   if(!tree || !node || !mesh) return RES_BAD_ARG;
    635   mesh->vertices = darray_vertex_cdata_get(&tree->vertices) + node->ivertex;
    636   mesh->nvertices = node->nvertices;
    637   return RES_OK;
    638 }
    639 
    640 double
    641 sln_node_eval
    642   (const struct sln_tree* tree,
    643    const struct sln_node* node,
    644    const double nu)
    645 {
    646   double ka = 0;
    647   size_t iline;
    648   ASSERT(tree && node);
    649 
    650   FOR_EACH(iline, node->range[0], node->range[1]+1) {
    651     struct line line = LINE_NULL;
    652     res_T res = RES_OK;
    653 
    654     res = line_setup(tree, iline, &line);
    655     if(res != RES_OK) {
    656       WARN(tree->sln, "%s: could not setup the line %lu-- %s\n",
    657         FUNC_NAME, iline, res_to_cstr(res));
    658     }
    659 
    660     ka += line_value(tree, &line, nu);
    661   }
    662   return ka;
    663 }
    664 
    665 res_T
    666 sln_node_get_desc(const struct sln_node* node, struct sln_node_desc* desc)
    667 {
    668   if(!node || !desc) return RES_BAD_ARG;
    669   desc->nlines  = node->range[1] - node->range[0];
    670   desc->nlines += 1/*boundaries are inclusives*/;
    671   desc->nvertices = node->nvertices;
    672   return RES_OK;
    673 }
    674 
    675 double
    676 sln_mesh_eval(const struct sln_mesh* mesh, const double wavenumber)
    677 {
    678   const struct sln_vertex* vtx0 = NULL;
    679   const struct sln_vertex* vtx1 = NULL;
    680   const float nu = (float)wavenumber;
    681   size_t n; /* #vertices */
    682   double u; /* Linear interpolation parameter */
    683   ASSERT(mesh && mesh->nvertices);
    684 
    685   n = mesh->nvertices;
    686 
    687   /* Handle special cases */
    688   if(n == 1) return mesh->vertices[0].ka;
    689   if(nu < mesh->vertices[0].wavenumber
    690   || nu > mesh->vertices[n-1].wavenumber) {
    691     return 0;
    692   }
    693   if(nu == mesh->vertices[0].wavenumber) return mesh->vertices[0].ka;
    694   if(nu == mesh->vertices[n-1].wavenumber) return mesh->vertices[n-1].ka;
    695 
    696   /* Dichotomic search of the mesh vertex whose wavenumber is greater than or
    697    * equal to the submitted wavenumber 'nu' */
    698   vtx1 = search_lower_bound(&nu, mesh->vertices, n, sizeof(*vtx1), cmp_nu_vtx);
    699   vtx0 = vtx1 - 1;
    700   ASSERT(vtx1); /* A vertex is necessary found ...*/
    701   ASSERT(vtx1 > mesh->vertices); /* ... and it cannot be the first one */
    702   ASSERT(vtx0->wavenumber < nu && nu <= vtx1->wavenumber);
    703 
    704   /* Compute the linear interpolation parameter */
    705   u = (wavenumber - vtx0->wavenumber) / (vtx1->wavenumber - vtx0->wavenumber);
    706   u = CLAMP(u, 0, 1); /* Handle numerical imprecisions */
    707 
    708   if(u == 0) return vtx0->ka;
    709   if(u == 1) return vtx1->ka;
    710   return u*(vtx1->ka - vtx0->ka) + vtx0->ka;
    711 }
    712 
    713 res_T
    714 sln_tree_write
    715   (const struct sln_tree* tree,
    716    const struct sln_tree_write_args* args)
    717 {
    718   struct stream stream = STREAM_NULL;
    719   size_t nnodes, nverts;
    720   hash256_T hash_mdata;
    721   hash256_T hash_lines;
    722   res_T res = RES_OK;
    723 
    724   if(!tree) { res = RES_BAD_ARG; goto error; }
    725   res = check_sln_tree_write_args(tree->sln, FUNC_NAME, args);
    726   if(res != RES_OK) goto error;
    727 
    728   res = shtr_isotope_metadata_hash(tree->args.metadata, hash_mdata);
    729   if(res != RES_OK) goto error;
    730   res = shtr_line_list_hash(tree->args.lines, hash_lines);
    731   if(res != RES_OK) goto error;
    732 
    733   res = stream_init
    734     (tree->sln, FUNC_NAME, args->filename, args->file, "w", &stream);
    735   if(res != RES_OK) goto error;
    736 
    737   #define WRITE(Var, Nb) {                                                     \
    738     if(fwrite((Var), sizeof(*(Var)), (Nb), stream.fp) != (Nb)) {               \
    739       ERROR(tree->sln, "%s:%s: error writing the tree -- %s\n",                \
    740         FUNC_NAME, stream.name, strerror(errno));                              \
    741       res = RES_IO_ERR;                                                        \
    742       goto error;                                                              \
    743     }                                                                          \
    744   } (void)0
    745   WRITE(&SLN_TREE_VERSION, 1);
    746   WRITE(hash_mdata, sizeof(hash256_T));
    747   WRITE(hash_lines, sizeof(hash256_T));
    748 
    749   nnodes = darray_node_size_get(&tree->nodes);
    750   WRITE(&nnodes, 1);
    751   WRITE(darray_node_cdata_get(&tree->nodes), nnodes);
    752 
    753   nverts = darray_vertex_size_get(&tree->vertices);
    754   WRITE(&nverts, 1);
    755   WRITE(darray_vertex_cdata_get(&tree->vertices), nverts);
    756 
    757   WRITE(&tree->args.line_profile, 1);
    758   WRITE(&tree->args.molecules, 1);
    759   WRITE(&tree->args.pressure, 1);
    760   WRITE(&tree->args.temperature, 1);
    761   WRITE(&tree->args.nvertices_hint, 1);
    762   WRITE(&tree->args.mesh_decimation_err, 1);
    763   WRITE(&tree->args.mesh_type, 1);
    764   WRITE(&tree->args.arity, 1);
    765   #undef WRITE
    766 
    767 exit:
    768   stream_release(&stream);
    769   return res;
    770 error:
    771   goto exit;
    772 }