htgop

Optical properties of a gas mixture
git clone git://git.meso-star.fr/htgop.git
Log | Files | Refs | README | LICENSE

htgop_fetch_radiative_properties.h (3539B)

      1 /* Copyright (C) 2018-2021, 2023 |Méso|Star> (contact@meso-star.com)
      2  *
      3  * This program is free software: you can redistribute it and/or modify
      4  * it under the terms of the GNU General Public License as published by
      5  * the Free Software Foundation, either version 3 of the License, or
      6  * (at your option) any later version.
      7  *
      8  * This program is distributed in the hope that it will be useful,
      9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     11  * GNU General Public License for more details.
     12  *
     13  * You should have received a copy of the GNU General Public License
     14  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     15 
     16 #include "htgop.h"
     17 #include "htgop_c.h"
     18 
     19 #include <rsys/algorithm.h>
     20 
     21 #if !defined(DATA) || !defined(DOMAIN)
     22   #error "Missing the <DATA|DOMAIN> macro."
     23 #endif
     24 
     25 /*
     26  * Generate the functions that fetch and interpolate a radiative property in a
     27  * given spectral domain
     28  */
     29 
     30 #ifndef GET_DATA
     31   #define GET_DATA(Tab, Id) ((Tab)->CONCAT(DATA, _tab)[Id])
     32 #endif
     33 
     34 /*******************************************************************************
     35  * Exported functions
     36  ******************************************************************************/
     37 res_T
     38 CONCAT(CONCAT(CONCAT(htgop_layer_,DOMAIN),_spectral_interval_tab_fetch_),DATA)
     39   (const struct CONCAT(CONCAT(htgop_layer_,DOMAIN),_spectral_interval_tab)* tab,
     40    const double x_h2o, /* H2O mol / mixture mol */
     41    double* out_k)
     42 {
     43   double* find;
     44 
     45   if(!tab || !out_k) return RES_BAD_ARG;
     46 
     47   find = search_lower_bound(&x_h2o, tab->x_h2o_tab, tab->tab_length,
     48     sizeof(double), cmp_dbl);
     49   if(!find) return RES_BAD_ARG;
     50 
     51   *out_k = CONCAT(CONCAT(CONCAT(
     52     layer_,DOMAIN),_spectral_interval_tab_interpolate_),DATA)(tab, x_h2o, find);
     53   return RES_OK;
     54 }
     55 
     56 /*******************************************************************************
     57  * Local functions
     58  ******************************************************************************/
     59 double
     60 CONCAT(CONCAT(CONCAT(layer_,DOMAIN),_spectral_interval_tab_interpolate_),DATA)
     61   (const struct CONCAT(CONCAT(htgop_layer_,DOMAIN),_spectral_interval_tab)* tab,
     62    const double x_h2o,
     63    const double* x_h2o_upp) /* Pointer toward the upper bound used to interpolate */
     64 {
     65   double x1, x2, k1, k2;
     66   double k;
     67 
     68   ASSERT(tab && x_h2o_upp);
     69 
     70   ASSERT(x_h2o_upp >= tab->x_h2o_tab);
     71   ASSERT(x_h2o_upp < tab->x_h2o_tab + tab->tab_length);
     72   ASSERT(*x_h2o_upp >= x_h2o);
     73   ASSERT(x_h2o_upp == tab->x_h2o_tab || x_h2o_upp[-1] < x_h2o);
     74 
     75   if(x_h2o_upp == tab->x_h2o_tab) {
     76     /* The submitted x_h2o is smaller that the first tabulated value of the
     77      * water vapor molar fraction. Use a simple linear interpolation to define
     78      * K by assuming that k1 and x1 is null */
     79     ASSERT(tab->x_h2o_tab[0] >= x_h2o);
     80     x2 = tab->x_h2o_tab[0];
     81     k2 = GET_DATA(tab, 0);
     82     k = k2 / x2*x_h2o;
     83   } else {
     84     const size_t i = (size_t)(x_h2o_upp - tab->x_h2o_tab);
     85     ASSERT(i < tab->tab_length);
     86     ASSERT(tab->x_h2o_tab[i-0] >= x_h2o);
     87     ASSERT(tab->x_h2o_tab[i-1] <  x_h2o);
     88     x1 = tab->x_h2o_tab[i-1];
     89     x2 = tab->x_h2o_tab[i-0];
     90     k1 = GET_DATA(tab, i-1);
     91     k2 = GET_DATA(tab, i-0);
     92 
     93     if(!k1 || !k2) { /* Linear interpolation */
     94       k = k1 + (k2-k1)/(x2-x1) * (x_h2o-x1);
     95     } else {
     96       const double alpha = (log(k2) - log(k1))/(log(x2) - log(x1));
     97       const double beta = log(k1) - alpha * log(x1);
     98       k = exp(alpha * log(x_h2o) + beta);
     99     }
    100   }
    101   return k;
    102 }
    103 
    104 #undef GET_DATA
    105 #undef DATA
    106 #undef DOMAIN