test_htgop_fetch_radiative_properties.h (3979B)
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 <math.h> 18 19 #if !defined(DOMAIN) || !defined(DATA) 20 #error "Missing the <DATA|DOMAIN> macro." 21 #endif 22 23 /* Helper macros */ 24 #define GET_NSPECINTS \ 25 CONCAT(CONCAT(htgop_get_, DOMAIN), _spectral_intervals_count) 26 #define GET_SPECINT \ 27 CONCAT(CONCAT(htgop_get_, DOMAIN), _spectral_interval) 28 #define LAY_SPECINT \ 29 CONCAT(CONCAT(htgop_layer_, DOMAIN),_spectral_interval) 30 #define LAY_SPECINT_TAB \ 31 CONCAT(CONCAT(htgop_layer_, DOMAIN),_spectral_interval_tab) 32 #define GET_LAY_SPECINT \ 33 CONCAT(CONCAT(htgop_layer_get_, DOMAIN), _spectral_interval) 34 #define GET_LAY_SPECINT_TAB \ 35 CONCAT(CONCAT(htgop_layer_, DOMAIN), _spectral_interval_get_tab) 36 #define FETCH_K \ 37 CONCAT(CONCAT(CONCAT(htgop_layer_, DOMAIN),_spectral_interval_tab_fetch_), DATA) 38 39 #ifndef GET_K 40 #define GET_K(Tab, Id) ((Tab)->CONCAT(DATA, _tab)[Id]) 41 #endif 42 43 static void 44 CONCAT(CONCAT(CONCAT(check_layer_fetch_, DOMAIN),_), DATA) 45 (const struct htgop* htgop, 46 const struct htgop_layer* layer) 47 { 48 struct htgop_spectral_interval specint; 49 size_t ispecint; 50 size_t nspecints; 51 double k; 52 53 CHK(htgop && layer); 54 CHK(GET_NSPECINTS(htgop, &nspecints) == RES_OK); 55 CHK(nspecints); 56 57 CHK(GET_SPECINT(htgop, 0, &specint) == RES_OK); 58 CHK(specint.quadrature_length); 59 CHK(layer->tab_length); 60 61 FOR_EACH(ispecint, 0, nspecints) { 62 struct LAY_SPECINT lay_specint; 63 struct LAY_SPECINT_TAB lay_tab; 64 size_t iquad; 65 66 CHK(GET_SPECINT(htgop, ispecint, &specint) == RES_OK); 67 CHK(specint.quadrature_length); 68 69 FOR_EACH(iquad, 0, specint.quadrature_length) { 70 size_t itab; 71 72 CHK(GET_LAY_SPECINT(layer, ispecint, &lay_specint) == RES_OK); 73 CHK(GET_LAY_SPECINT_TAB(&lay_specint, iquad, &lay_tab) == RES_OK); 74 CHK(lay_tab.tab_length == layer->tab_length); 75 76 FOR_EACH(itab, 0, layer->tab_length) { 77 double x_h2o; 78 const size_t N = 10; 79 size_t i; 80 81 /* Check boundaries */ 82 x_h2o = layer->x_h2o_tab[itab]; 83 CHK(FETCH_K(&lay_tab, x_h2o, &k) == RES_OK); 84 CHK(eq_eps(k, GET_K(&lay_tab,itab), 1.e-6)); 85 86 FOR_EACH(i, 0, N) { 87 const double r = rand_canonic(); 88 const double x1 = itab ? layer->x_h2o_tab[itab-1] : 0; 89 const double x2 = layer->x_h2o_tab[itab]; 90 const double k1 = itab ? GET_K(&lay_tab, itab-1) : 0; 91 const double k2 = GET_K(&lay_tab, itab); 92 double ref; 93 94 x_h2o = x1 + (x2 - x1) * r; 95 if(!itab || !k1 || !k2) { 96 ref = k1 + (k2-k1) / (x2 - x1) * x_h2o; 97 } else { 98 const double alpha = (log(k2) - log(k1)) / (log(x2) - log(x1)); 99 const double beta = log(k1) - alpha*log(x1); 100 ref = exp(alpha*log(x_h2o) + beta); 101 } 102 103 CHK(FETCH_K(&lay_tab, x_h2o, &k) == RES_OK); 104 CHK(eq_eps(k, ref, 1.e-6)); 105 #ifdef Kext 106 { 107 double ka, ks; 108 CHK(FETCH_Ka(&lay_tab, x_h2o, &ka) == RES_OK); 109 CHK(FETCH_Ks(&lay_tab ,x_h2o, &ks) == RES_OK); 110 CHK(ka <= k && ks <= k); 111 } 112 #endif 113 } 114 } 115 } 116 } 117 } 118 119 #undef GET_NSPECINTS 120 #undef GET_SPECINT 121 #undef LAY_SPECINT 122 #undef LAY_SPECINT_TAB 123 #undef GET_LAY_SPECINT 124 #undef GET_LAY_SPECINT_TAB 125 #undef GET_K 126 #undef DOMAIN 127 #undef DATA 128