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 (21382B)


      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   res = check_molecules(sln, caller, args);
    227   if(res != RES_OK) return res;
    228 
    229   return RES_OK;
    230 }
    231 
    232 static res_T
    233 check_sln_tree_read_args
    234   (struct sln_device* sln,
    235    const char* caller,
    236    const struct sln_tree_read_args* args)
    237 {
    238   if(!args) return RES_BAD_ARG;
    239 
    240   if(!args->metadata) {
    241     ERROR(sln, "%s: the isotope metadata are missing.\n", caller);
    242     return RES_BAD_ARG;
    243   }
    244 
    245   if(!args->lines) {
    246     ERROR(sln, "%s: the list of lines is missing.\n", caller);
    247     return RES_BAD_ARG;
    248   }
    249 
    250   if(!args->file && !args->filename) {
    251     ERROR(sln,
    252       "%s: the source file is missing. No file name or stream is provided.\n",
    253       caller);
    254     return RES_BAD_ARG;
    255   }
    256 
    257   return RES_OK;
    258 }
    259 
    260 static res_T
    261 check_sln_tree_write_args
    262   (struct sln_device* sln,
    263    const char* caller,
    264    const struct sln_tree_write_args* args)
    265 {
    266   if(!args) return RES_BAD_ARG;
    267 
    268   if(!args->file && !args->filename) {
    269     ERROR(sln,
    270       "%s: the destination file is missing. "
    271       "No file name or stream is provided.\n",
    272       caller);
    273     return RES_BAD_ARG;
    274   }
    275 
    276   return RES_OK;
    277 }
    278 static INLINE void
    279 stream_release(struct stream* stream)
    280 {
    281   ASSERT(stream);
    282   if(stream->intern_fp && stream->fp) CHK(fclose(stream->fp) == 0);
    283 }
    284 
    285 static res_T
    286 stream_init
    287   (struct sln_device* sln,
    288    const char* caller,
    289    const char* name, /* NULL <=> default stream name */
    290    FILE* fp, /* NULL <=> open file "name" */
    291    const char* mode, /* mode in fopen */
    292    struct stream* stream)
    293 {
    294   res_T res = RES_OK;
    295 
    296   ASSERT(sln && caller && stream);
    297   ASSERT(fp || (name && mode));
    298 
    299   *stream = STREAM_NULL;
    300 
    301   if(fp) {
    302     stream->intern_fp = 0;
    303     stream->name = name ? name : "stream";
    304     stream->fp = fp;
    305 
    306   } else {
    307     stream->intern_fp = 1;
    308     stream->name = name;
    309     if(!(stream->fp = fopen(name, mode))) {
    310       ERROR(sln, "%s:%s: error opening file -- %s\n",
    311         caller, name, strerror(errno));
    312       res = RES_IO_ERR;
    313       goto error;
    314     }
    315   }
    316 
    317 exit:
    318   return res;
    319 error:
    320   if(stream->intern_fp && stream->fp) CHK(fclose(stream->fp) == 0);
    321   goto exit;
    322 }
    323 
    324 static res_T
    325 create_tree
    326   (struct sln_device* sln,
    327    const char* caller,
    328    struct sln_tree** out_tree)
    329 {
    330   struct sln_tree* tree = NULL;
    331   res_T res = RES_OK;
    332   ASSERT(sln && caller && out_tree);
    333 
    334   tree = MEM_CALLOC(sln->allocator, 1, sizeof(struct sln_tree));
    335   if(!tree) {
    336     ERROR(sln, "%s: could not allocate the tree data structure.\n",
    337       caller);
    338     res = RES_MEM_ERR;
    339     goto error;
    340   }
    341   ref_init(&tree->ref);
    342   SLN(device_ref_get(sln));
    343   tree->sln = sln;
    344   darray_node_init(sln->allocator, &tree->nodes);
    345   darray_vertex_init(sln->allocator, &tree->vertices);
    346 
    347 exit:
    348   *out_tree = tree;
    349   return res;
    350 error:
    351   if(tree) { SLN(tree_ref_put(tree)); tree = NULL; }
    352   goto exit;
    353 }
    354 
    355 static INLINE int
    356 cmp_nu_vtx(const void* key, const void* item)
    357 {
    358   const float nu = *((const float*)key);
    359   const struct sln_vertex* vtx = item;
    360   if(nu < vtx->wavenumber) return -1;
    361   if(nu > vtx->wavenumber) return +1;
    362   return 0;
    363 }
    364 
    365 static void
    366 release_tree(ref_T* ref)
    367 {
    368   struct sln_tree* tree = CONTAINER_OF(ref, struct sln_tree, ref);
    369   struct sln_device* sln = NULL;
    370   ASSERT(ref);
    371   sln = tree->sln;
    372   darray_node_release(&tree->nodes);
    373   darray_vertex_release(&tree->vertices);
    374   if(tree->args.lines) SHTR(line_list_ref_put(tree->args.lines));
    375   if(tree->args.metadata) SHTR(isotope_metadata_ref_put(tree->args.metadata));
    376   MEM_RM(sln->allocator, tree);
    377   SLN(device_ref_put(sln));
    378 }
    379 
    380 /*******************************************************************************
    381  * Exported symbols
    382  ******************************************************************************/
    383 res_T
    384 sln_tree_create
    385   (struct sln_device* device,
    386    const struct sln_tree_create_args* args,
    387    struct sln_tree** out_tree)
    388 {
    389   struct sln_tree* tree = NULL;
    390   res_T res = RES_OK;
    391 
    392   if(!device || !out_tree) { res = RES_BAD_ARG; goto error; }
    393   res = check_sln_tree_create_args(device, FUNC_NAME, args);
    394   if(res != RES_OK) goto error;
    395 
    396   res = create_tree(device, FUNC_NAME, &tree);
    397   if(res != RES_OK) goto error;
    398   SHTR(line_list_ref_get(args->lines));
    399   SHTR(isotope_metadata_ref_get(args->metadata));
    400   tree->args = *args;
    401 
    402   res = tree_build(tree);
    403   if(res != RES_OK) goto error;
    404 
    405 exit:
    406   if(out_tree) *out_tree = tree;
    407   return res;
    408 error:
    409   if(tree) { SLN(tree_ref_put(tree)); tree = NULL; }
    410   goto exit;
    411 }
    412 
    413 res_T
    414 sln_tree_read
    415   (struct sln_device* sln,
    416    const struct sln_tree_read_args* args,
    417    struct sln_tree** out_tree)
    418 {
    419   hash256_T hash_mdata1;
    420   hash256_T hash_mdata2;
    421   hash256_T hash_lines1;
    422   hash256_T hash_lines2;
    423 
    424   struct stream stream = STREAM_NULL;
    425   struct sln_tree* tree = NULL;
    426   size_t n = 0;
    427   int version = 0;
    428   res_T res = RES_OK;
    429 
    430   if(!sln || !out_tree) { res = RES_BAD_ARG; goto error; }
    431   res = check_sln_tree_read_args(sln, FUNC_NAME, args);
    432   if(res != RES_OK) goto error;
    433 
    434   res = create_tree(sln, FUNC_NAME, &tree);
    435   if(res != RES_OK) goto error;
    436 
    437   res = stream_init(sln, FUNC_NAME, args->filename, args->file, "r", &stream);
    438   if(res != RES_OK) goto error;
    439 
    440   #define READ(Var, Nb) {                                                      \
    441     if(fread((Var), sizeof(*(Var)), (Nb), stream.fp) != (Nb)) {                \
    442       if(feof(stream.fp)) {                                                    \
    443         res = RES_BAD_ARG;                                                     \
    444       } else if(ferror(stream.fp)) {                                           \
    445         res = RES_IO_ERR;                                                      \
    446       } else {                                                                 \
    447         res = RES_UNKNOWN_ERR;                                                 \
    448       }                                                                        \
    449       ERROR(sln, "%s: error loading the tree structure -- %s.\n",              \
    450         stream.name, res_to_cstr(res));                                        \
    451       goto error;                                                              \
    452     }                                                                          \
    453   } (void)0
    454   READ(&version, 1);
    455   if(version != SLN_TREE_VERSION) {
    456     ERROR(sln,
    457       "%s: unexpected tree version %d. Expecting a tree in version %d.\n",
    458       stream.name, version, SLN_TREE_VERSION);
    459     res = RES_BAD_ARG;
    460     goto error;
    461   }
    462 
    463   res = shtr_isotope_metadata_hash(args->metadata, hash_mdata1);
    464   if(res != RES_OK) goto error;
    465 
    466   READ(hash_mdata2, sizeof(hash256_T));
    467   if(!hash256_eq(hash_mdata1, hash_mdata2)) {
    468     ERROR(sln,
    469       "%s: the input isotopic metadata are not those used "
    470       "during tree construction.\n", stream.name);
    471     res = RES_BAD_ARG;
    472     goto error;
    473   }
    474 
    475   SHTR(isotope_metadata_ref_get(args->metadata));
    476   tree->args.metadata = args->metadata;
    477 
    478   res = shtr_line_list_hash(args->lines, hash_lines1);
    479   if(res != RES_OK) goto error;
    480 
    481   READ(hash_lines2, sizeof(hash256_T));
    482   if(!hash256_eq(hash_lines1, hash_lines2)) {
    483     ERROR(sln,
    484       "%s: the input list of lines is not the one used to build the tree.\n",
    485       stream.name);
    486     res = RES_BAD_ARG;
    487     goto error;
    488   }
    489 
    490   SHTR(line_list_ref_get(args->lines));
    491   tree->args.lines = args->lines;
    492 
    493   READ(&n, 1);
    494   if((res = darray_node_resize(&tree->nodes, n)) != RES_OK) goto error;
    495   READ(darray_node_data_get(&tree->nodes), n);
    496 
    497   READ(&n, 1);
    498   if((res = darray_vertex_resize(&tree->vertices, n)) != RES_OK) goto error;
    499   READ(darray_vertex_data_get(&tree->vertices), n);
    500 
    501   READ(&tree->args.line_profile, 1);
    502   READ(&tree->args.molecules, 1);
    503   READ(&tree->args.pressure, 1);
    504   READ(&tree->args.temperature, 1);
    505   READ(&tree->args.nvertices_hint, 1);
    506   READ(&tree->args.mesh_decimation_err, 1);
    507   READ(&tree->args.mesh_type, 1);
    508   #undef READ
    509 
    510 exit:
    511   stream_release(&stream);
    512   if(out_tree) *out_tree = tree;
    513   return res;
    514 error:
    515   if(tree) { SLN(tree_ref_put(tree)); tree = NULL; }
    516   goto exit;
    517 }
    518 
    519 res_T
    520 sln_tree_ref_get(struct sln_tree* tree)
    521 {
    522   if(!tree) return RES_BAD_ARG;
    523   ref_get(&tree->ref);
    524   return RES_OK;
    525 }
    526 
    527 res_T
    528 sln_tree_ref_put(struct sln_tree* tree)
    529 {
    530   if(!tree) return RES_BAD_ARG;
    531   ref_put(&tree->ref, release_tree);
    532   return RES_OK;
    533 }
    534 
    535 res_T
    536 sln_tree_get_desc(const struct sln_tree* tree, struct sln_tree_desc* desc)
    537 {
    538   size_t nlines_adjusted = 0;
    539   unsigned depth = 0;
    540 
    541   if(!tree || !desc) return RES_BAD_ARG;
    542 
    543   desc->max_nlines_per_leaf = 1;
    544   desc->mesh_decimation_err = tree->args.mesh_decimation_err;
    545   desc->mesh_type = tree->args.mesh_type;
    546   desc->line_profile = tree->args.line_profile;
    547   desc->nnodes = darray_node_size_get(&tree->nodes);
    548   desc->nvertices = darray_vertex_size_get(&tree->vertices);
    549 
    550   SHTR(line_list_get_size(tree->args.lines, &desc->nlines));
    551   nlines_adjusted = round_up_pow2(desc->nlines);
    552 
    553   for(depth=0; depth<64 && !(BIT_U64(depth) & nlines_adjusted); ++depth);
    554   desc->depth = depth;
    555 
    556 #ifndef NDEBUG
    557   {
    558     unsigned max_depth = 0;
    559     const struct sln_node* node = sln_tree_get_root(tree);
    560     while(!sln_node_is_leaf(node)) {
    561       node = sln_node_get_child(node, 0);
    562       ++max_depth;
    563     }
    564     CHK(max_depth == depth);
    565   }
    566 #endif
    567 
    568   return RES_OK;
    569 }
    570 
    571 const struct sln_node*
    572 sln_tree_get_root(const struct sln_tree* tree)
    573 {
    574   ASSERT(tree);
    575   if(darray_node_size_get(&tree->nodes)) {
    576     return darray_node_cdata_get(&tree->nodes);
    577   } else {
    578     return NULL;
    579   }
    580 }
    581 
    582 int
    583 sln_node_is_leaf(const struct sln_node* node)
    584 {
    585   ASSERT(node);
    586   return node->offset == 0;
    587 }
    588 
    589 const struct sln_node*
    590 sln_node_get_child(const struct sln_node* node, const unsigned ichild)
    591 {
    592   ASSERT(node && ichild <= 1);
    593   ASSERT(!sln_node_is_leaf(node));
    594   return node + node->offset + ichild;
    595 }
    596 
    597 size_t
    598 sln_node_get_lines_count(const struct sln_node* node)
    599 {
    600   ASSERT(node);
    601   return node->range[1] - node->range[0] + 1/*Both boundaries are inclusives*/;
    602 }
    603 
    604 res_T
    605 sln_node_get_line
    606   (const struct sln_tree* tree,
    607    const struct sln_node* node,
    608    const size_t iline,
    609    struct shtr_line* line)
    610 {
    611   if(!tree || !node || iline > sln_node_get_lines_count(node))
    612     return RES_BAD_ARG;
    613 
    614   return shtr_line_list_at
    615     (tree->args.lines, node->range[0] + iline, line);
    616 }
    617 
    618 res_T
    619 sln_node_get_mesh
    620   (const struct sln_tree* tree,
    621    const struct sln_node* node,
    622    struct sln_mesh* mesh)
    623 {
    624   if(!tree || !node || !mesh) return RES_BAD_ARG;
    625   mesh->vertices = darray_vertex_cdata_get(&tree->vertices) + node->ivertex;
    626   mesh->nvertices = node->nvertices;
    627   return RES_OK;
    628 }
    629 
    630 double
    631 sln_node_eval
    632   (const struct sln_tree* tree,
    633    const struct sln_node* node,
    634    const double nu)
    635 {
    636   double ka = 0;
    637   size_t iline;
    638   ASSERT(tree && node);
    639 
    640   FOR_EACH(iline, node->range[0], node->range[1]+1) {
    641     struct line line = LINE_NULL;
    642     res_T res = RES_OK;
    643 
    644     res = line_setup(tree, iline, &line);
    645     if(res != RES_OK) {
    646       WARN(tree->sln, "%s: could not setup the line %lu-- %s\n",
    647         FUNC_NAME, iline, res_to_cstr(res));
    648     }
    649 
    650     ka += line_value(tree, &line, nu);
    651   }
    652   return ka;
    653 }
    654 
    655 res_T
    656 sln_node_get_desc(const struct sln_node* node, struct sln_node_desc* desc)
    657 {
    658   if(!node || !desc) return RES_BAD_ARG;
    659   desc->nlines  = node->range[1] - node->range[0];
    660   desc->nlines += 1/*boundaries are inclusives*/;
    661   desc->nvertices = node->nvertices;
    662   return RES_OK;
    663 }
    664 
    665 double
    666 sln_mesh_eval(const struct sln_mesh* mesh, const double wavenumber)
    667 {
    668   const struct sln_vertex* vtx0 = NULL;
    669   const struct sln_vertex* vtx1 = NULL;
    670   const float nu = (float)wavenumber;
    671   size_t n; /* #vertices */
    672   float u; /* Linear interpolation paraeter */
    673   ASSERT(mesh && mesh->nvertices);
    674 
    675   n = mesh->nvertices;
    676 
    677   /* Handle special cases */
    678   if(n == 1) return mesh->vertices[0].ka;
    679   if(nu <= mesh->vertices[0].wavenumber) return mesh->vertices[0].ka;
    680   if(nu >= mesh->vertices[n-1].wavenumber) return mesh->vertices[n-1].ka;
    681 
    682   /* Dichotomic search of the mesh vertex whose wavenumber is greater than or
    683    * equal to the submitted wavenumber 'nu' */
    684   vtx1 = search_lower_bound(&nu, mesh->vertices, n, sizeof(*vtx1), cmp_nu_vtx);
    685   vtx0 = vtx1 - 1;
    686   ASSERT(vtx1); /* A vertex is necessary found ...*/
    687   ASSERT(vtx1 > mesh->vertices); /* ... and it cannot be the first one */
    688   ASSERT(vtx0->wavenumber < nu && nu <= vtx1->wavenumber);
    689 
    690   /* Compute the linear interpolation parameter */
    691   u = (nu - vtx0->wavenumber) / (vtx1->wavenumber - vtx0->wavenumber);
    692   u = CLAMP(u, 0, 1); /* Handle numerical imprecisions */
    693 
    694   if(u == 0) return vtx0->ka;
    695   if(u == 1) return vtx1->ka;
    696   return u*(vtx1->ka - vtx0->ka) + vtx0->ka;
    697 }
    698 
    699 res_T
    700 sln_tree_write
    701   (const struct sln_tree* tree,
    702    const struct sln_tree_write_args* args)
    703 {
    704   struct stream stream = STREAM_NULL;
    705   size_t nnodes, nverts;
    706   hash256_T hash_mdata;
    707   hash256_T hash_lines;
    708   res_T res = RES_OK;
    709 
    710   if(!tree) { res = RES_BAD_ARG; goto error; }
    711   res = check_sln_tree_write_args(tree->sln, FUNC_NAME, args);
    712   if(res != RES_OK) goto error;
    713 
    714   res = shtr_isotope_metadata_hash(tree->args.metadata, hash_mdata);
    715   if(res != RES_OK) goto error;
    716   res = shtr_line_list_hash(tree->args.lines, hash_lines);
    717   if(res != RES_OK) goto error;
    718 
    719   res = stream_init
    720     (tree->sln, FUNC_NAME, args->filename, args->file, "w", &stream);
    721   if(res != RES_OK) goto error;
    722 
    723   #define WRITE(Var, Nb) {                                                     \
    724     if(fwrite((Var), sizeof(*(Var)), (Nb), stream.fp) != (Nb)) {               \
    725       ERROR(tree->sln, "%s:%s: error writing the tree -- %s\n",                \
    726         FUNC_NAME, stream.name, strerror(errno));                              \
    727       res = RES_IO_ERR;                                                        \
    728       goto error;                                                              \
    729     }                                                                          \
    730   } (void)0
    731   WRITE(&SLN_TREE_VERSION, 1);
    732   WRITE(hash_mdata, sizeof(hash256_T));
    733   WRITE(hash_lines, sizeof(hash256_T));
    734 
    735   nnodes = darray_node_size_get(&tree->nodes);
    736   WRITE(&nnodes, 1);
    737   WRITE(darray_node_cdata_get(&tree->nodes), nnodes);
    738 
    739   nverts = darray_vertex_size_get(&tree->vertices);
    740   WRITE(&nverts, 1);
    741   WRITE(darray_vertex_cdata_get(&tree->vertices), nverts);
    742 
    743   WRITE(&tree->args.line_profile, 1);
    744   WRITE(&tree->args.molecules, 1);
    745   WRITE(&tree->args.pressure, 1);
    746   WRITE(&tree->args.temperature, 1);
    747   WRITE(&tree->args.nvertices_hint, 1);
    748   WRITE(&tree->args.mesh_decimation_err, 1);
    749   WRITE(&tree->args.mesh_type, 1);
    750   #undef WRITE
    751 
    752 exit:
    753   stream_release(&stream);
    754   return res;
    755 error:
    756   goto exit;
    757 }