htgop

htgop.c (30706B)

---

 /* Copyright (C) 2018-2021, 2023 |Méso|Star> (contact@meso-star.com)
  *
  * 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 /* nextafter support */
 
 #include "htgop.h"
 #include "htgop_c.h"
 #include "htgop_reader.h"
 
 #include <rsys/algorithm.h>
 #include <rsys/cstr.h>
 #include <rsys/logger.h>
 #include <rsys/mem_allocator.h>
 
 #include <math.h>
 
 /*******************************************************************************
  * Helper functions
  ******************************************************************************/
 static void
 log_msg
   (const struct htgop* htgop,
    const enum log_type stream,
    const char* msg,
    va_list vargs)
 {
   ASSERT(htgop && msg);
   if(htgop->verbose) {
     res_T res; (void)res;
     res = logger_vprint(htgop->logger, stream, msg, vargs);
     ASSERT(res == RES_OK);
   }
 }
 
 static INLINE int
 cmp_lvl(const void* key, const void* data)
 {
   const double height = *((const double*)key);
   const struct htgop_level* lvl = (const struct htgop_level*)data;
   return height < lvl->height ? -1 : (height > lvl->height ? 1 : 0);
 }
 
 static res_T
 parse_spectral_intervals
   (struct htgop* htgop,
    struct reader* rdr,
    const unsigned long nspecints,
    struct spectral_intervals* specints)
 {
   unsigned long ispecint;
   double* wave_numbers = NULL;
   struct darray_double* pdfs = NULL;
   struct darray_double* cdfs = NULL;
   res_T res = RES_OK;
   (void)htgop;
   ASSERT(htgop && rdr && nspecints && specints);
 
   #define CALL(Func) { if((res = Func) != RES_OK) goto error; } (void)0
 
   /* Allocate the memory space for the wave numbers and the quadratures */
   CALL(darray_double_resize(&specints->wave_numbers, nspecints + 1));
   CALL(darray_dbllst_resize(&specints->quadrature_pdfs, nspecints));
   CALL(darray_dbllst_resize(&specints->quadrature_cdfs, nspecints));
   wave_numbers = darray_double_data_get(&specints->wave_numbers);
   pdfs = darray_dbllst_data_get(&specints->quadrature_pdfs);
   cdfs = darray_dbllst_data_get(&specints->quadrature_cdfs);
 
   FOR_EACH(ispecint, 0, nspecints) {
     double* pdf = NULL;
     double* cdf = NULL;
     double wave_number_low;
     double wave_number_upp;
     double sum = 0;
     unsigned long quad_len;
     unsigned long iquad;
 
     /* Parse the interval bounds */
     CALL(cstr_to_double(read_line(rdr), &wave_number_low));
     CALL(cstr_to_double(read_line(rdr), &wave_number_upp));
 
     /* Check and register the interval bounds */
     if(wave_number_low >= wave_number_upp) {
       res = RES_BAD_ARG;
       goto error;
     }
     if(ispecint == 0) {
       wave_numbers[ispecint] = wave_number_low;
     } else if(wave_number_low != wave_numbers[ispecint]) {
       res = RES_BAD_ARG;
       goto error;
     }
     wave_numbers[ispecint + 1] = wave_number_upp;
 
     /* Parse and allocate the quadrature length */
     CALL(cstr_to_ulong(read_line(rdr), &quad_len));
 
     /* Allocate the points of the quadrature */
     CALL(darray_double_resize(&pdfs[ispecint], quad_len));
     CALL(darray_double_resize(&cdfs[ispecint], quad_len));
     pdf = darray_double_data_get(&pdfs[ispecint]);
     cdf = darray_double_data_get(&cdfs[ispecint]);
 
     /* Read the weight of the quadrature points */
     sum = 0;
     FOR_EACH(iquad, 0, quad_len) {
       CALL(cstr_to_double(read_line(rdr), &pdf[iquad]));
       cdf[iquad] = iquad == 0 ? pdf[iquad] : pdf[iquad]+cdf[iquad-1];
       sum += pdf[iquad];
     }
     ASSERT(eq_eps(sum, 1.0, 1.e6));
     ASSERT(eq_eps(cdf[quad_len-1], 1.0, 1.e6));
     cdf[quad_len-1] = 1.0; /* Handle numerical imprecision */
   }
 
   #undef CALL
 
 exit:
   return res;
 error:
   darray_double_clear(&specints->wave_numbers);
   darray_dbllst_clear(&specints->quadrature_pdfs);
   darray_dbllst_clear(&specints->quadrature_cdfs);
   goto exit;
 }
 
 /* Generate the parse_layers_spectral_intervals_ka_lw function */
 #define DOMAIN lw
 #define DATA ka
 #include "htgop_parse_layers_spectral_intervals_data.h"
 /* Generate the parse_layers_spectral_intervals_ka_sw function */
 #define DOMAIN sw
 #define DATA ka
 #include "htgop_parse_layers_spectral_intervals_data.h"
 /* Generate the parse_layers_spectral_intervals_ks_sw function */
 #define DOMAIN sw
 #define DATA ks
 #include "htgop_parse_layers_spectral_intervals_data.h"
 
 static res_T
 parse_layers_spectral_intervals(struct htgop* htgop, struct reader* rdr)
 {
   struct layer* layers = NULL;
   size_t lw_nspecints, sw_nspecints, nlays;
   size_t ilay;
   res_T res = RES_OK;
   ASSERT(htgop && rdr);
 
   layers = darray_layer_data_get(&htgop->layers);
   nlays = darray_layer_size_get(&htgop->layers);
   ASSERT(nlays > 0);
   lw_nspecints = darray_double_size_get(&htgop->lw_specints.wave_numbers) - 1;
   sw_nspecints = darray_double_size_get(&htgop->sw_specints.wave_numbers) - 1;
   ASSERT(lw_nspecints > 0 && sw_nspecints);
 
   #define CALL(Func) { if((res = Func) != RES_OK) goto error; } (void)0
   /* Allocate the per layer spectral intervals */
   FOR_EACH(ilay, 0, nlays) {
     CALL(darray_lay_lw_specint_resize(&layers[ilay].lw_specints, lw_nspecints));
     CALL(darray_lay_sw_specint_resize(&layers[ilay].sw_specints, sw_nspecints));
   }
   /* Parse the spectral data of the layers */
   CALL(parse_layers_spectral_intervals_ka_lw(htgop, rdr));
   CALL(parse_layers_spectral_intervals_ka_sw(htgop, rdr));
   CALL(parse_layers_spectral_intervals_ks_sw(htgop, rdr));
   #undef CALL
 
 exit:
   return res;
 error:
   FOR_EACH(ilay, 0, nlays) {
     darray_lay_lw_specint_clear(&layers[ilay].lw_specints);
     darray_lay_sw_specint_clear(&layers[ilay].sw_specints);
   }
   goto exit;
 }
 
 static res_T
 load_stream(struct htgop* htgop, FILE* stream, const char* stream_name)
 {
   struct reader rdr;
   struct htgop_level* levels = NULL;
   struct layer* layers = NULL;
   unsigned long nlvls, nlays, tab_len, nspecints_lw, nspecints_sw;
   unsigned long ilvl, ilay, itab;
   unsigned long iline_saved;
   res_T res = RES_OK;
   ASSERT(htgop && stream && stream_name);
   reader_init(&rdr, stream, stream_name);
 
   #define CALL(Func) {                                                         \
     if((res = Func) != RES_OK) {                                               \
       log_err(htgop, "%s:%lu: Parsing error.\n", rdr.name, rdr.iline);         \
       goto error;                                                              \
     }                                                                          \
   } (void)0
 
   /* Parse the number of levels/layers */
   CALL(cstr_to_ulong(read_line(&rdr), &nlvls));
   iline_saved = rdr.iline;
   CALL(cstr_to_ulong(read_line(&rdr), &nlays));
   if(!nlays) {
     log_err(htgop, "%s:%lu: The number of layers cannot be null.\n",
       rdr.name, rdr.iline);
     res = RES_BAD_ARG;
     goto error;
   }
   if(nlvls != nlays + 1) {
     log_err(htgop,
       "%s:%lu: Invalid number of levels `%lu' (#layers = `%lu').\n",
       rdr.name, iline_saved, nlvls, nlays);
     res = RES_BAD_ARG;
     goto error;
   }
 
   /* Parse the ground temperature */
   CALL(cstr_to_double(read_line(&rdr), &htgop->ground.temperature));
 
   /* Allocate the per level data */
   CALL(darray_level_resize(&htgop->levels, nlvls));
   CALL(darray_layer_resize(&htgop->layers, nlays));
   levels = darray_level_data_get(&htgop->levels);
   layers = darray_layer_data_get(&htgop->layers);
 
   /* Per level data */
   FOR_EACH(ilvl, 0, nlvls) { /* Pressure */
     CALL(cstr_to_double(read_line(&rdr), &levels[ilvl].pressure));
   }
   FOR_EACH(ilvl, 0, nlvls) { /* Temperature */
     CALL(cstr_to_double(read_line(&rdr), &levels[ilvl].temperature));
   }
   FOR_EACH(ilvl, 0, nlvls) { /* Height */
     CALL(cstr_to_double(read_line(&rdr), &levels[ilvl].height));
     if(ilvl && levels[ilvl].height <= levels[ilvl-1].height) {
       log_err(htgop,
         "%s:%lu: The levels must be sorted in strict ascending order "
         "wrt their height.\n", rdr.name, rdr.iline);
       res = RES_BAD_ARG;
       goto error;
     }
   }
 
   /* Per layer x H2O nominal */
   FOR_EACH(ilay, 0, nlays) {
     CALL(cstr_to_double(read_line(&rdr), &layers[ilay].x_h2o_nominal));
   }
 
   /* Parse the length of the tabulation */
   CALL(cstr_to_ulong(read_line(&rdr), &tab_len));
   if(!tab_len) {
     log_err(htgop,
       "%s:%lu: The tabulation length cannot be null.\n", rdr.name, rdr.iline);
     res = RES_BAD_ARG;
     goto error;
   }
 
   /* Allocate the tabulated xH2O of each layer */
   FOR_EACH(ilay, 0, nlays) {
     CALL(darray_double_resize(&layers[ilay].x_h2o_tab, tab_len));
   }
 
   /* Parse the tabulated xH2O of each layer */
   FOR_EACH(itab, 0, tab_len) {
     FOR_EACH(ilay, 0, nlays) {
       double* x_h2o_tab = darray_double_data_get(&layers[ilay].x_h2o_tab);
       CALL(cstr_to_double(read_line(&rdr), &x_h2o_tab[itab]));
       if(x_h2o_tab[itab] < 0) {
         log_err(htgop,
           "%s:%lu: the tabulated water vapor molar fractions must be in [0, 1].\n",
           rdr.name, rdr.iline);
         res = RES_BAD_ARG;
         goto error;
       }
       if(itab && x_h2o_tab[itab] < x_h2o_tab[itab-1]) {
         log_err(htgop,
           "%s:%lu: the tabulated water vapor molar fractions must be sorted in "
           "strict ascending order.\n", rdr.name, rdr.iline);
         res = RES_BAD_ARG;
         goto error;
       }
       if(itab == tab_len-1 && x_h2o_tab[itab] != 1) {
         log_err(htgop,
         "%s:%lu: the tabulated water vapor molar fractions is not normalized.\n",
         rdr.name, rdr.iline);
         res = RES_BAD_ARG;
         goto error;
       }
     }
   }
 
   /* Parse the long/short wave emissivity of the hround */
   CALL(cstr_to_double(read_line(&rdr), &htgop->ground.lw_emissivity));
   CALL(cstr_to_double(read_line(&rdr), &htgop->ground.sw_emissivity));
 
   /* Parse the number of long/short wave spectral intervals */
   CALL(cstr_to_ulong(read_line(&rdr), &nspecints_lw));
   CALL(cstr_to_ulong(read_line(&rdr), &nspecints_sw));
 
   /* Parse the data of the long/short wave spectral intervals */
   CALL(parse_spectral_intervals(htgop, &rdr, nspecints_lw, &htgop->lw_specints));
   CALL(parse_spectral_intervals(htgop, &rdr, nspecints_sw, &htgop->sw_specints));
 
   /* Parse the spectral data of the layers */
   CALL(parse_layers_spectral_intervals(htgop, &rdr));
 
   #undef CALL
 
 exit:
   reader_release(&rdr);
   return res;
 error:
   goto exit;
 }
 
 static void
 release_htgop(ref_T* ref)
 {
   struct htgop* htgop;
   ASSERT(ref);
   htgop = CONTAINER_OF(ref, struct htgop, ref);
   spectral_intervals_release(&htgop->lw_specints);
   spectral_intervals_release(&htgop->sw_specints);
   darray_double_release(&htgop->sw_X_cdf);
   darray_double_release(&htgop->sw_Y_cdf);
   darray_double_release(&htgop->sw_Z_cdf);
   darray_level_release(&htgop->levels);
   darray_layer_release(&htgop->layers);
   MEM_RM(htgop->allocator, htgop);
 }
 
 /*******************************************************************************
  * Exported functions
  ******************************************************************************/
 res_T
 htgop_create
   (struct logger* log,
    struct mem_allocator* mem_allocator,
    const int verbose,
    struct htgop** out_htgop)
 {
   struct mem_allocator* allocator = NULL;
   struct logger* logger = NULL;
   struct htgop* htgop = NULL;
   res_T res = RES_OK;
 
   if(!out_htgop) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   allocator = mem_allocator ? mem_allocator : &mem_default_allocator;
   logger = log ? log : LOGGER_DEFAULT;
 
   htgop = MEM_CALLOC(allocator, 1, sizeof(struct htgop));
   if(!htgop) {
     if(verbose) {
       logger_print(logger, LOG_ERROR,
         "Could not allocate the HTGOP handler.\n");
     }
     res = RES_MEM_ERR;
     goto error;
   }
 
   ref_init(&htgop->ref);
   htgop->allocator = allocator;
   htgop->logger = logger;
   htgop->verbose = verbose;
   spectral_intervals_init(allocator, &htgop->lw_specints);
   spectral_intervals_init(allocator, &htgop->sw_specints);
   darray_double_init(allocator, &htgop->sw_X_cdf);
   darray_double_init(allocator, &htgop->sw_Y_cdf);
   darray_double_init(allocator, &htgop->sw_Z_cdf);
   darray_level_init(allocator, &htgop->levels);
   darray_layer_init(allocator, &htgop->layers);
 
 exit:
   if(out_htgop) *out_htgop = htgop;
   return res;
 error:
   if(htgop) {
     HTGOP(ref_put(htgop));
     htgop = NULL;
   }
   goto exit;
 }
 
 res_T
 htgop_ref_get(struct htgop* htgop)
 {
   if(!htgop) return RES_BAD_ARG;
   ref_get(&htgop->ref);
   return RES_OK;
 }
 
 res_T
 htgop_ref_put(struct htgop* htgop)
 {
   if(!htgop) return RES_BAD_ARG;
   ref_put(&htgop->ref, release_htgop);
   return RES_OK;
 }
 
 res_T
 htgop_load(struct htgop* htgop, const char* filename)
 {
   FILE* file = NULL;
   res_T res = RES_OK;
 
   if(!htgop || !filename) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   file = fopen(filename, "r");
   if(!file) {
     log_err(htgop, "%s: error opening file `%s'.\n", FUNC_NAME, filename);
     res = RES_IO_ERR;
     goto error;
   }
 
   res = load_stream(htgop, file, filename);
   if(res != RES_OK) goto error;
 
 exit:
   if(file) fclose(file);
   return res;
 error:
   goto exit;
 }
 
 res_T
 htgop_load_stream(struct htgop* htgop, FILE* stream)
 {
   if(!htgop || !stream) return RES_BAD_ARG;
   return load_stream(htgop, stream, "<stream>");
 }
 
 res_T
 htgop_get_ground(const struct htgop* htgop, struct htgop_ground* ground)
 {
   if(!htgop || !ground) return RES_BAD_ARG;
   *ground = htgop->ground;
   return RES_OK;
 }
 
 res_T
 htgop_get_layers_count(const struct htgop* htgop, size_t* nlayers)
 {
   if(!htgop || !nlayers) return RES_BAD_ARG;
   *nlayers = darray_layer_size_get(&htgop->layers);
   return RES_OK;
 }
 
 res_T
 htgop_get_levels_count(const struct htgop* htgop, size_t* nlevels)
 {
   if(!htgop || !nlevels) return RES_BAD_ARG;
   *nlevels = darray_level_size_get(&htgop->levels);
   return RES_OK;
 }
 
 res_T
 htgop_get_level
   (const struct htgop* htgop, const size_t ilvl, struct htgop_level* lvl)
 {
   size_t n;
   if(!htgop || !lvl) return RES_BAD_ARG;
   HTGOP(get_levels_count(htgop, &n));
   if(ilvl >= n) {
     log_err(htgop, "%s: invalid level `%lu'.\n", FUNC_NAME, ilvl);
     return RES_BAD_ARG;
   }
   *lvl = darray_level_cdata_get(&htgop->levels)[ilvl];
   return RES_OK;
 }
 
 res_T
 htgop_get_layer
   (const struct htgop* htgop,
    const size_t ilayer,
    struct htgop_layer* layer)
 {
   const struct layer* l;
   size_t n;
   if(!htgop || !layer) return RES_BAD_ARG;
   HTGOP(get_layers_count(htgop, &n));
   if(ilayer >= n) {
     log_err(htgop, "%s: invalid layer `%lu'.\n",
       FUNC_NAME, (unsigned long)ilayer);
     return RES_BAD_ARG;
   }
   l = darray_layer_cdata_get(&htgop->layers) + ilayer;
   layer->x_h2o_nominal = l->x_h2o_nominal;
   layer->x_h2o_tab = darray_double_cdata_get(&l->x_h2o_tab);
   layer->tab_length = darray_double_size_get(&l->x_h2o_tab);
   layer->lw_spectral_intervals_count =
     darray_lay_lw_specint_size_get(&l->lw_specints);
   layer->sw_spectral_intervals_count =
     darray_lay_sw_specint_size_get(&l->sw_specints);
   layer->data__ = l;
   layer->htgop = htgop;
   return RES_OK;
 }
 
 res_T
 htgop_layer_get_lw_spectral_interval
   (const struct htgop_layer* layer,
    const size_t ispecint,
    struct htgop_layer_lw_spectral_interval* specint)
 {
   const struct layer* l;
   const struct layer_lw_spectral_interval* lspecint;
   if(!layer || !specint) return RES_BAD_ARG;
   if(ispecint >= layer->lw_spectral_intervals_count) {
     log_err(layer->htgop, "%s: invalid spectral interval `%lu'.\n",
       FUNC_NAME, (unsigned long)ispecint);
     return RES_BAD_ARG;
   }
   l = layer->data__;
   lspecint = darray_lay_lw_specint_cdata_get(&l->lw_specints) + ispecint;
   specint->ka_nominal = darray_double_cdata_get(&lspecint->ka_nominal);
   specint->quadrature_length = darray_double_size_get(&lspecint->ka_nominal);
   specint->htgop = layer->htgop;
   specint->data__ = lspecint;
   specint->layer__ = l;
   return RES_OK;
 }
 
 res_T
 htgop_layer_get_sw_spectral_interval
   (const struct htgop_layer* layer,
    const size_t ispecint,
    struct htgop_layer_sw_spectral_interval* specint)
 {
   const struct layer* l;
   const struct layer_sw_spectral_interval* lspecint;
   if(!layer || !specint) return RES_BAD_ARG;
   if(ispecint >= layer->sw_spectral_intervals_count) {
     log_err(layer->htgop, "%s: invalid spectral interval `%lu'.\n",
       FUNC_NAME, (unsigned long)ispecint);
     return RES_BAD_ARG;
   }
   l = layer->data__;
   lspecint = darray_lay_sw_specint_cdata_get(&l->sw_specints) + ispecint;
   specint->ka_nominal = darray_double_cdata_get(&lspecint->ka_nominal);
   specint->ks_nominal = darray_double_cdata_get(&lspecint->ks_nominal);
   specint->quadrature_length = darray_double_size_get(&lspecint->ka_nominal);
   specint->htgop = layer->htgop;
   specint->data__ = lspecint;
   specint->layer__ = l;
   return RES_OK;
 }
 
 res_T
 htgop_layer_lw_spectral_interval_get_tab
   (const struct htgop_layer_lw_spectral_interval* specint,
    const size_t iquad,
    struct htgop_layer_lw_spectral_interval_tab* tab)
 {
   const struct layer_lw_spectral_interval* lspecint;
   const struct layer* l;
   const struct darray_double* t;
   if(!specint || !tab) return RES_BAD_ARG;
   if(iquad >= specint->quadrature_length) {
     log_err(specint->htgop, "%s: invalid quadrature point `%lu'.\n",
       FUNC_NAME, (unsigned long)iquad);
     return RES_BAD_ARG;
   }
   lspecint = specint->data__;
   l = specint->layer__;
   t = darray_dbllst_cdata_get(&lspecint->ka_tab) + iquad;
   tab->x_h2o_tab = darray_double_cdata_get(&l->x_h2o_tab);
   tab->ka_tab = darray_double_cdata_get(t);
   tab->tab_length = darray_double_size_get(t);
   tab->htgop = specint->htgop;
   return RES_OK;
 }
 
 res_T
 htgop_layer_sw_spectral_interval_get_tab
   (const struct htgop_layer_sw_spectral_interval* specint,
    const size_t iquad,
    struct htgop_layer_sw_spectral_interval_tab* tab)
 {
   const struct layer_sw_spectral_interval* lspecint;
   const struct layer* l;
   const struct darray_double* ka_tab;
   const struct darray_double* ks_tab;
   if(!specint || !tab) return RES_BAD_ARG;
   if(iquad >= specint->quadrature_length) {
     log_err(specint->htgop, "%s: invalid quadrature point `%lu'.\n",
       FUNC_NAME, (unsigned long)iquad);
     return RES_BAD_ARG;
   }
   lspecint = specint->data__;
   l = specint->layer__;
 
   ka_tab = darray_dbllst_cdata_get(&lspecint->ka_tab) + iquad;
   ks_tab = darray_dbllst_cdata_get(&lspecint->ks_tab) + iquad;
   tab->x_h2o_tab = darray_double_cdata_get(&l->x_h2o_tab);
   tab->ka_tab = darray_double_cdata_get(ka_tab);
   tab->ks_tab = darray_double_cdata_get(ks_tab);
   tab->tab_length = darray_double_size_get(ka_tab);
   tab->htgop = specint->htgop;
   return RES_OK;
 }
 
 res_T
 htgop_get_lw_spectral_intervals_count(const struct htgop* htgop, size_t* nspecints)
 {
   if(!htgop || !nspecints) return RES_BAD_ARG;
   *nspecints = darray_dbllst_size_get(&htgop->lw_specints.quadrature_pdfs);
   return RES_OK;
 }
 
 res_T
 htgop_get_sw_spectral_intervals_count(const struct htgop* htgop, size_t* nspecints)
 {
   if(!htgop || !nspecints) return RES_BAD_ARG;
   *nspecints = darray_dbllst_size_get(&htgop->sw_specints.quadrature_pdfs);
   return RES_OK;
 }
 
 res_T
 htgop_get_lw_spectral_intervals_wave_numbers
   (const struct htgop* htgop, const double* wave_numbers[])
 {
   if(!htgop || !wave_numbers) return RES_BAD_ARG;
   *wave_numbers = darray_double_cdata_get(&htgop->lw_specints.wave_numbers);
   return RES_OK;
 }
 
 res_T
 htgop_get_sw_spectral_intervals_wave_numbers
   (const struct htgop* htgop, const double* wave_numbers[])
 {
   if(!htgop || !wave_numbers) return RES_BAD_ARG;
   *wave_numbers = darray_double_cdata_get(&htgop->sw_specints.wave_numbers);
   return RES_OK;
 }
 
 res_T
 htgop_get_lw_spectral_interval
   (const struct htgop* htgop,
    const size_t ispecint,
    struct htgop_spectral_interval* interval)
 {
   const double* wave_numbers;
   const struct darray_double* quad_cdfs;
   const struct darray_double* quad_pdfs;
   size_t n;
   if(!htgop || !interval) return RES_BAD_ARG;
   HTGOP(get_lw_spectral_intervals_count(htgop, &n));
   if(ispecint >= n) {
     log_err(htgop, "%s: invalid long wave spectral interval `%lu'.\n",
       FUNC_NAME, (unsigned long)ispecint);
     return RES_BAD_ARG;
   }
   wave_numbers = darray_double_cdata_get(&htgop->lw_specints.wave_numbers);
   quad_pdfs = darray_dbllst_cdata_get(&htgop->lw_specints.quadrature_pdfs);
   quad_cdfs = darray_dbllst_cdata_get(&htgop->lw_specints.quadrature_cdfs);
   interval->wave_numbers[0] = wave_numbers[ispecint+0];
   interval->wave_numbers[1] = wave_numbers[ispecint+1];
   interval->quadrature_pdf = darray_double_cdata_get(&quad_pdfs[ispecint]);
   interval->quadrature_cdf = darray_double_cdata_get(&quad_cdfs[ispecint]);
   interval->quadrature_length = darray_double_size_get(&quad_pdfs[ispecint]);
   interval->htgop = htgop;
   return RES_OK;
 }
 
 res_T
 htgop_get_sw_spectral_interval
   (const struct htgop* htgop,
    const size_t ispecint,
    struct htgop_spectral_interval* interval)
 {
   const double* wave_numbers;
   const struct darray_double* quad_pdfs;
   const struct darray_double* quad_cdfs;
   size_t n;
   if(!htgop || !interval) return RES_BAD_ARG;
   HTGOP(get_sw_spectral_intervals_count(htgop, &n));
   if(ispecint >= n) {
     log_err(htgop, "%s: invalid short wave spectral interval `%lu'.\n",
       FUNC_NAME, (unsigned long)ispecint);
     return RES_BAD_ARG;
   }
   wave_numbers = darray_double_cdata_get(&htgop->sw_specints.wave_numbers);
   quad_pdfs = darray_dbllst_cdata_get(&htgop->sw_specints.quadrature_pdfs);
   quad_cdfs = darray_dbllst_cdata_get(&htgop->sw_specints.quadrature_cdfs);
   interval->wave_numbers[0] = wave_numbers[ispecint+0];
   interval->wave_numbers[1] = wave_numbers[ispecint+1];
   interval->quadrature_pdf = darray_double_cdata_get(&quad_pdfs[ispecint]);
   interval->quadrature_cdf = darray_double_cdata_get(&quad_cdfs[ispecint]);
   interval->quadrature_length = darray_double_size_get(&quad_pdfs[ispecint]);
   interval->htgop = htgop;
   return RES_OK;
 }
 
 res_T
 htgop_find_lw_spectral_interval_id
   (const struct htgop* htgop,
    const double wnum,
    size_t* ispecint)
 {
   const double* find = NULL;
   const double* wnums = NULL;
   size_t nwnums = 0;
   if(!htgop || !ispecint) return RES_BAD_ARG;
 
   wnums = darray_double_cdata_get(&htgop->lw_specints.wave_numbers);
   nwnums = darray_double_size_get(&htgop->lw_specints.wave_numbers);
 
   find = search_lower_bound(&wnum, wnums, nwnums, sizeof(*wnums), cmp_dbl);
   if(!find || find == wnums) {
     log_warn(htgop, "%s: the wavenumber `%g' is not included in any band.\n",
       FUNC_NAME, wnum);
     *ispecint = SIZE_MAX;
   } else {
     *ispecint = (size_t)(find - wnums - 1);
 #ifndef NDEBUG
     {
       size_t n;
       HTGOP(get_lw_spectral_intervals_count(htgop, &n));
       ASSERT(*ispecint < n);
     }
 #endif
   }
   return RES_OK;
 }
 
 res_T
 htgop_find_sw_spectral_interval_id
   (const struct htgop* htgop,
    const double wnum,
    size_t* ispecint)
 {
   const double* find = NULL;
   const double* wnums = NULL;
   size_t nwnums = 0;
   if(!htgop || !ispecint) return RES_BAD_ARG;
 
   wnums = darray_double_cdata_get(&htgop->sw_specints.wave_numbers);
   nwnums = darray_double_size_get(&htgop->sw_specints.wave_numbers);
 
   find = search_lower_bound(&wnum, wnums, nwnums, sizeof(*wnums), cmp_dbl);
   if(!find || find == wnums) {
     log_warn(htgop, "%s: the wavenumber `%g' is not included in any band.\n",
       FUNC_NAME, wnum);
     *ispecint = SIZE_MAX;
   } else {
     *ispecint = (size_t)(find - wnums - 1);
 #ifndef NDEBUG
     {
       size_t n;
       HTGOP(get_sw_spectral_intervals_count(htgop, &n));
       ASSERT(*ispecint < n);
     }
 #endif
 
   }
   return RES_OK;
 }
 
 res_T
 htgop_spectral_interval_sample_quadrature
   (const struct htgop_spectral_interval* specint,
    const double r, /* Canonical number in [0, 1[ */
    size_t* iquad_point) /* Id of the sample quadrature point */
 {
   double r_next = nextafter(r, DBL_MAX);
   double* find;
   size_t i;
 
   if(!specint || !iquad_point) return RES_BAD_ARG;
 
   if(r < 0 || r >= 1) {
     log_err(specint->htgop,
       "%s: invalid canonical random number `%g'.\n", FUNC_NAME, r);
     return RES_BAD_ARG;
   }
 
   /* Use r_next rather than r in order to find the first entry that is not less
    * than *or equal* to r */
   find = search_lower_bound(&r_next, specint->quadrature_cdf,
     specint->quadrature_length, sizeof(double), cmp_dbl);
   ASSERT(find);
   i = (size_t)(find - specint->quadrature_cdf);
   ASSERT(i < specint->quadrature_length);
   ASSERT(specint->quadrature_cdf[i] > r);
   ASSERT(!i || specint->quadrature_cdf[i-1] <= r);
   *iquad_point = i;
   return RES_OK;
 }
 
 res_T
 htgop_position_to_layer_id
   (const struct htgop* htgop, const double pos, size_t* ilayer)
 {
   const struct htgop_level* lvls;
   size_t nlvls;
 
   if(!htgop || !ilayer) return RES_BAD_ARG;
 
   lvls = darray_level_cdata_get(&htgop->levels);
   nlvls = darray_level_size_get(&htgop->levels);
   ASSERT(nlvls);
 
   if(pos == lvls[0].height) {
     *ilayer = 0;
   } else if(pos == lvls[nlvls-1].height) {
     *ilayer = nlvls - 2;
   } else {
     const struct htgop_level* find;
     size_t i;
 
     find = search_lower_bound(&pos, lvls, nlvls, sizeof(*lvls), cmp_lvl);
     if(!find || find == lvls) {
       log_err(htgop, "%s: the position `%g' is outside the atmospheric slices.\n",
         FUNC_NAME, pos);
       return RES_BAD_ARG;
     }
 
     i = (size_t)(find - lvls);
     ASSERT(i < nlvls && i > 0);
     ASSERT(lvls[i].height > pos && (!i || lvls[i-1].height <= pos));
     *ilayer = i - 1;
   }
   return RES_OK;
 }
 
 res_T
 htgop_get_sw_spectral_intervals
   (const struct htgop* htgop,
    const double wnum_range[2],
    size_t specint_range[2])
 {
   const double* wnums = NULL;
   size_t nwnums = 0;
   res_T res = RES_OK;
 
   if(!htgop || !wnum_range || wnum_range[0]>wnum_range[1] || !specint_range) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   wnums = darray_double_cdata_get(&htgop->sw_specints.wave_numbers);
   nwnums = darray_double_size_get(&htgop->sw_specints.wave_numbers);
 
   res = get_spectral_intervals
     (htgop, FUNC_NAME, wnum_range, wnums, nwnums, specint_range);
   if(res != RES_OK) goto error;
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 res_T
 htgop_get_lw_spectral_intervals
   (const struct htgop* htgop,
    const double wnum_range[2],
    size_t specint_range[2])
 {
   const double* wnums = NULL;
   size_t nwnums = 0;
   res_T res = RES_OK;
 
   if(!htgop || !wnum_range || wnum_range[0]>wnum_range[1] || !specint_range) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   wnums = darray_double_cdata_get(&htgop->lw_specints.wave_numbers);
   nwnums = darray_double_size_get(&htgop->lw_specints.wave_numbers);
 
   res = get_spectral_intervals
     (htgop, FUNC_NAME, wnum_range, wnums, nwnums, specint_range);
   if(res != RES_OK) goto error;
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 /* Generate the htgop_layer_fetch_lw_ka function */
 #define DOMAIN lw
 #define DATA ka
 #include "htgop_fetch_radiative_properties.h"
 
 /* Generate the htgop_layer_fetch_sw_ka function */
 #define DOMAIN sw
 #define DATA ka
 #include "htgop_fetch_radiative_properties.h"
 
 /* Generate the htgop_layer_fetch_sw_ks function */
 #define DOMAIN sw
 #define DATA ks
 #include "htgop_fetch_radiative_properties.h"
 
 /* Generate the htgop_layer_get_sw_kext function */
 #define GET_DATA(Tab, Id) ((Tab)->ka_tab[Id] + (Tab)->ks_tab[Id])
 #define DOMAIN sw
 #define DATA kext
 #include "htgop_fetch_radiative_properties.h"
 
 /* Generate the functions that get the boundaries of ka in LW */
 #define DOMAIN lw
 #define DATA ka
 #include "htgop_get_radiative_properties_bounds.h"
 
 /* Generate the functions that get the boundaries of ka in SW */
 #define DOMAIN sw
 #define DATA ka
 #include "htgop_get_radiative_properties_bounds.h"
 
 /* Generate the functions that get the boundaries of ks in SW */
 #define DOMAIN sw
 #define DATA ks
 #include "htgop_get_radiative_properties_bounds.h"
 
 /* Generate the functions that get the boundaries of kext in SW */
 #define GET_DATA(Tab, Id) ((Tab)->ka_tab[Id] + (Tab)->ks_tab[Id])
 #define DOMAIN sw
 #define DATA kext
 #include "htgop_get_radiative_properties_bounds.h"
 
 /*******************************************************************************
  * Local functions
  ******************************************************************************/
 void
 log_err(const struct htgop* htgop, const char* msg, ...)
 {
   va_list vargs_list;
   ASSERT(htgop && msg);
 
   va_start(vargs_list, msg);
   log_msg(htgop, LOG_ERROR, msg, vargs_list);
   va_end(vargs_list);
 }
 
 void
 log_warn(const struct htgop* htgop, const char* msg, ...)
 {
   va_list vargs_list;
   ASSERT(htgop && msg);
 
   va_start(vargs_list, msg);
   log_msg(htgop, LOG_WARNING, msg, vargs_list);
   va_end(vargs_list);
 }
 
 res_T
 get_spectral_intervals
   (const struct htgop* htgop,
    const char* func_name,
    const double wnum_range[2],
    const double* wnums,
    const size_t nwnums,
    size_t specint_range[2])
 {
   double wnum_min = 0;
   double wnum_max = 0;
   double* low = NULL;
   double* upp = NULL;
   res_T res = RES_OK;
   ASSERT(htgop && wnum_range && wnums && nwnums && specint_range);
   ASSERT(wnum_range[0] <= wnum_range[1]);
 
   wnum_min = nextafter(wnum_range[0], DBL_MAX);
   wnum_max = wnum_range[1];
   low = search_lower_bound(&wnum_min, wnums, nwnums, sizeof(double), cmp_dbl);
   upp = search_lower_bound(&wnum_max, wnums, nwnums, sizeof(double), cmp_dbl);
 
   if(!low || upp == wnums) {
     log_err(htgop,
       "%s: the loaded data do not overlap the wave numbers in [%g, %g].\n",
       func_name, wnum_range[0], wnum_range[1]);
     res = RES_BAD_OP;
     goto error;
   }
   ASSERT(low[0] >= wnum_min && (low == wnums || low[-1] < wnum_min));
   ASSERT(!upp || (upp[0] >= wnum_max && upp[-1] < wnum_max));
 
   specint_range[0] = low == wnums ? 0 : (size_t)(low - wnums) - 1;
   specint_range[1] = !upp ? nwnums - 1 : (size_t)(upp - wnums);
   /* Transform the upper wnum bound in spectral interval id */
   specint_range[1] = specint_range[1] - 1;
 
 
   ASSERT(specint_range[0] <= specint_range[1]);
   ASSERT(wnums[specint_range[0]+1] >  wnum_range[0]);
   ASSERT(wnums[specint_range[1]+0] <  wnum_range[1]);
 
   ASSERT(wnums[specint_range[0]+0] <= wnum_range[0]
     || specint_range[0]==0); /* The data do not include `wnum_range' */
   ASSERT(wnums[specint_range[1]+1] >= wnum_range[1]
     || specint_range[1] + 1 == nwnums-1);/* The data do not include `wnum_range' */
 
 exit:
   return res;
 error:
   if(specint_range) {
     specint_range[0] = SIZE_MAX;
     specint_range[1] = 0;
   }
   goto exit;
 }