test_htmie_load.c (13178B)
1 /* Copyright (C) 2018, 2020-2023 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2018 Centre National de la Recherche Scientifique 3 * Copyright (C) 2018 Université Paul Sabatier 4 * 5 * This program is free software: you can redistribute it and/or modify 6 * it under the terms of the GNU General Public License as published by 7 * the Free Software Foundation, either version 3 of the License, or 8 * (at your option) any later version. 9 * 10 * This program is distributed in the hope that it will be useful, 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 * GNU General Public License for more details. 14 * 15 * You should have received a copy of the GNU General Public License 16 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 17 18 #define _POSIX_C_SOURCE 200112L /* snprintf support */ 19 20 #include "htmie.h" 21 #include "test_htmie_utils.h" 22 23 #include <rsys/math.h> 24 25 #include <stdlib.h> 26 #include <string.h> 27 28 29 #define NEAREST HTMIE_FILTER_NEAREST 30 #define LINEAR HTMIE_FILTER_LINEAR 31 32 static INLINE double 33 rand_canonic(void) 34 { 35 return rand() / (double)(RAND_MAX - 1); 36 } 37 38 static void 39 test_fetch(struct htmie* htmie) 40 { 41 size_t nwlens; 42 size_t i; 43 const double* wlens; 44 const double* absor; 45 const double* scatt; 46 const double* asym; 47 double wlen; 48 double xabs; 49 double xsca; 50 double u; 51 double g; 52 53 CHK(htmie); 54 55 nwlens = htmie_get_wavelengths_count(htmie); 56 CHK(nwlens > 1); 57 58 wlens = htmie_get_wavelengths(htmie); 59 scatt = htmie_get_xsections_scattering(htmie); 60 absor = htmie_get_xsections_absorption(htmie); 61 asym = htmie_get_asymmetry_parameter(htmie); 62 63 /* Test clamp to border */ 64 CHK(absor[0] == htmie_fetch_xsection_absorption(htmie, wlens[0]-1, NEAREST)); 65 CHK(scatt[0] == htmie_fetch_xsection_scattering(htmie, wlens[0]-1, NEAREST)); 66 CHK(asym[0] == htmie_fetch_asymmetry_parameter(htmie, wlens[0]-1, NEAREST)); 67 CHK(absor[0] == htmie_fetch_xsection_absorption(htmie, wlens[0]-1, LINEAR)); 68 CHK(scatt[0] == htmie_fetch_xsection_scattering(htmie, wlens[0]-1, LINEAR)); 69 CHK(asym[0] == htmie_fetch_asymmetry_parameter(htmie, wlens[0]-1, LINEAR)); 70 71 i = nwlens - 1; 72 CHK(absor[i] == htmie_fetch_xsection_absorption(htmie, wlens[i]+1, NEAREST)); 73 CHK(scatt[i] == htmie_fetch_xsection_scattering(htmie, wlens[i]+1, NEAREST)); 74 CHK(asym[i] == htmie_fetch_asymmetry_parameter(htmie, wlens[i]+1, NEAREST)); 75 CHK(absor[i] == htmie_fetch_xsection_absorption(htmie, wlens[i]+1, LINEAR)); 76 CHK(scatt[i] == htmie_fetch_xsection_scattering(htmie, wlens[i]+1, LINEAR)); 77 CHK(asym[i] == htmie_fetch_asymmetry_parameter(htmie, wlens[i]+1, LINEAR)); 78 79 80 FOR_EACH(i, 0, nwlens-1) { 81 CHK(absor[i+0] == htmie_fetch_xsection_absorption(htmie, wlens[i+0], NEAREST)); 82 CHK(absor[i+1] == htmie_fetch_xsection_absorption(htmie, wlens[i+1], NEAREST)); 83 CHK(scatt[i+0] == htmie_fetch_xsection_scattering(htmie, wlens[i+0], NEAREST)); 84 CHK(scatt[i+1] == htmie_fetch_xsection_scattering(htmie, wlens[i+1], NEAREST)); 85 CHK(asym[i+1] == htmie_fetch_asymmetry_parameter(htmie, wlens[i+1], NEAREST)); 86 87 u = 0.25; 88 wlen = wlens[i+0] + u * (wlens[i+1] - wlens[i+0]); 89 xabs = absor[i+0] + u * (absor[i+1] - absor[i+0]); 90 xsca = scatt[i+0] + u * (scatt[i+1] - scatt[i+0]); 91 g = asym[i+0] + u *(asym[i+1] - asym[i+0]); 92 CHK(absor[i+0] == htmie_fetch_xsection_absorption(htmie, wlen, NEAREST)); 93 CHK(scatt[i+0] == htmie_fetch_xsection_scattering(htmie, wlen, NEAREST)); 94 CHK(asym[i+0] == htmie_fetch_asymmetry_parameter(htmie, wlen, NEAREST)); 95 CHK(eq_eps(xabs, htmie_fetch_xsection_absorption(htmie, wlen, LINEAR), 1.e-6)); 96 CHK(eq_eps(xsca, htmie_fetch_xsection_scattering(htmie, wlen, LINEAR), 1.e-6)); 97 CHK(eq_eps(g, htmie_fetch_asymmetry_parameter(htmie, wlen, LINEAR), 1.e-6)); 98 99 u = 0.51; 100 wlen = wlens[i+0] + u * (wlens[i+1] - wlens[i+0]); 101 xabs = absor[i+0] + u * (absor[i+1] - absor[i+0]); 102 xsca = scatt[i+0] + u * (scatt[i+1] - scatt[i+0]); 103 g = asym[i+0] + u *(asym[i+1] - asym[i+0]); 104 CHK(absor[i+1] == htmie_fetch_xsection_absorption(htmie, wlen, NEAREST)); 105 CHK(scatt[i+1] == htmie_fetch_xsection_scattering(htmie, wlen, NEAREST)); 106 CHK(asym[i+1] == htmie_fetch_asymmetry_parameter(htmie, wlen, NEAREST)); 107 CHK(eq_eps(xabs, htmie_fetch_xsection_absorption(htmie, wlen, LINEAR), 1.e-6)); 108 CHK(eq_eps(xsca, htmie_fetch_xsection_scattering(htmie, wlen, LINEAR), 1.e-6)); 109 CHK(eq_eps(g, htmie_fetch_asymmetry_parameter(htmie, wlen, LINEAR), 1.e-6)); 110 } 111 } 112 113 static void 114 test_bounds(struct htmie* htmie) 115 { 116 double bounds[2]; 117 double band[2]; 118 double min_val = DBL_MAX; 119 double max_val =-DBL_MAX; 120 const double* wlens; 121 const size_t ntests = 128; 122 size_t ilow; 123 size_t iupp; 124 size_t nwlens; 125 size_t i; 126 size_t itest; 127 CHK(htmie); 128 129 wlens = htmie_get_wavelengths(htmie); 130 nwlens = htmie_get_wavelengths_count(htmie); 131 132 CHK(nwlens > 1); 133 134 FOR_EACH(itest, 0, ntests) { 135 do { 136 ilow = (size_t)(rand_canonic() * (double)nwlens); 137 iupp = (size_t)(rand_canonic() * (double)nwlens); 138 } while(ilow == iupp); 139 140 if(ilow > iupp) SWAP(size_t, ilow, iupp); 141 band[0] = wlens[ilow]; 142 band[1] = wlens[iupp]; 143 144 #define TEST(Name, Filter) { \ 145 htmie_compute_xsection_## Name ## _bounds(htmie, band, Filter, bounds); \ 146 min_val = htmie_fetch_xsection_ ## Name(htmie, band[0], Filter); \ 147 max_val = htmie_fetch_xsection_ ## Name(htmie, band[1], Filter); \ 148 if(min_val > max_val) SWAP(double, min_val, max_val); \ 149 FOR_EACH(i, ilow+1, iupp) { \ 150 double tmp = htmie_fetch_xsection_ ## Name(htmie, wlens[i], Filter); \ 151 min_val = MMIN(min_val, tmp); \ 152 max_val = MMAX(max_val, tmp); \ 153 } \ 154 CHK(eq_eps(bounds[0], min_val, 1.e-6)); \ 155 CHK(eq_eps(bounds[1], max_val, 1.e-6)); \ 156 } (void)0 157 158 TEST(absorption, NEAREST); 159 TEST(scattering, NEAREST); 160 TEST(absorption, LINEAR); 161 TEST(scattering, LINEAR); 162 163 band[0] = wlens[ilow] + 0.25*(wlens[ilow+1] - wlens[ilow+0]); 164 band[1] = wlens[iupp] - 0.25*(wlens[iupp+0] - wlens[iupp-1]); 165 CHK(band[0] < band[1]); 166 167 TEST(absorption, NEAREST); 168 TEST(scattering, NEAREST); 169 TEST(absorption, LINEAR); 170 TEST(scattering, LINEAR); 171 172 #undef TEST 173 174 band[0] = band[1]; 175 176 htmie_compute_xsection_absorption_bounds(htmie, band, NEAREST, bounds); 177 min_val = htmie_fetch_xsection_absorption(htmie, band[0], NEAREST); 178 max_val = htmie_fetch_xsection_absorption(htmie, band[1], NEAREST); 179 CHK(eq_eps(bounds[0], min_val, 1.e-6)); 180 CHK(eq_eps(bounds[1], max_val, 1.e-6)); 181 182 htmie_compute_xsection_absorption_bounds(htmie, band, LINEAR, bounds); 183 min_val = htmie_fetch_xsection_absorption(htmie, band[0], LINEAR); 184 max_val = htmie_fetch_xsection_absorption(htmie, band[1], LINEAR); 185 CHK(eq_eps(bounds[0], min_val, 1.e-6)); 186 CHK(eq_eps(bounds[1], max_val, 1.e-6)); 187 188 htmie_compute_xsection_scattering_bounds(htmie, band, NEAREST, bounds); 189 min_val = htmie_fetch_xsection_scattering(htmie, band[0], NEAREST); 190 max_val = htmie_fetch_xsection_scattering(htmie, band[1], NEAREST); 191 CHK(eq_eps(bounds[0], min_val, 1.e-6)); 192 CHK(eq_eps(bounds[1], max_val, 1.e-6)); 193 194 htmie_compute_xsection_scattering_bounds(htmie, band, LINEAR, bounds); 195 min_val = htmie_fetch_xsection_scattering(htmie, band[0], LINEAR); 196 max_val = htmie_fetch_xsection_scattering(htmie, band[1], LINEAR); 197 CHK(eq_eps(bounds[0], min_val, 1.e-6)); 198 CHK(eq_eps(bounds[1], max_val, 1.e-6)); 199 } 200 } 201 202 static void 203 test_avg(struct htmie* htmie) 204 { 205 double band[2]; 206 const double* wlens; 207 const size_t ntests = 128; 208 size_t itest; 209 size_t nwlens; 210 size_t ilow; 211 size_t iupp; 212 size_t i; 213 CHK(htmie); 214 215 wlens = htmie_get_wavelengths(htmie); 216 nwlens = htmie_get_wavelengths_count(htmie); 217 218 CHK(nwlens > 1); 219 220 FOR_EACH(itest, 0, ntests) { 221 do { 222 ilow = (size_t)(rand_canonic() * (double)nwlens); 223 iupp = (size_t)(rand_canonic() * (double)nwlens); 224 } while(ilow == iupp); 225 226 if(ilow > iupp) SWAP(size_t, ilow, iupp); 227 band[0] = wlens[ilow]; 228 band[1] = wlens[iupp]; 229 230 #define TEST(Name, Filter) { \ 231 double avg = 0; \ 232 double sum = 0; \ 233 double data = 0; \ 234 double prev_wlen = band[0]; \ 235 double prev_data = htmie_fetch_ ## Name(htmie, band[0], Filter);\ 236 FOR_EACH(i, ilow+1, iupp) { \ 237 data = htmie_fetch_ ## Name(htmie, wlens[i], Filter); \ 238 sum += 0.5*(prev_data + data) * (wlens[i] - prev_wlen); \ 239 prev_wlen = wlens[i]; \ 240 prev_data = data; \ 241 } \ 242 data = htmie_fetch_ ## Name(htmie, band[1], Filter); \ 243 sum += 0.5*(prev_data + data) * (band[1] - prev_wlen); \ 244 sum /= (band[1] - band[0]); \ 245 avg = htmie_compute_## Name ## _average(htmie, band, Filter); \ 246 CHK(eq_eps(avg, sum, 1.e-6)); \ 247 } (void)0 248 249 TEST(xsection_absorption, NEAREST); 250 TEST(xsection_scattering, NEAREST); 251 TEST(asymmetry_parameter, NEAREST); 252 TEST(xsection_absorption, LINEAR); 253 TEST(xsection_scattering, LINEAR); 254 TEST(asymmetry_parameter, LINEAR); 255 256 band[0] = wlens[ilow] + 0.25*(wlens[ilow+1] - wlens[ilow+0]); 257 band[1] = wlens[iupp] - 0.25*(wlens[iupp+0] - wlens[iupp-1]); 258 CHK(band[0] < band[1]); 259 260 TEST(xsection_absorption, NEAREST); 261 TEST(xsection_scattering, NEAREST); 262 TEST(asymmetry_parameter, NEAREST); 263 TEST(xsection_absorption, LINEAR); 264 TEST(xsection_scattering, LINEAR); 265 TEST(asymmetry_parameter, LINEAR); 266 267 #undef TEST 268 } 269 } 270 271 int 272 main(int argc, char** argv) 273 { 274 char buf[128]; 275 struct mem_allocator allocator; 276 char* filename = NULL; 277 char* path = NULL; 278 char* base = NULL; 279 char* p = NULL; 280 FILE* fp = NULL; 281 struct htmie* htmie = NULL; 282 size_t i; 283 284 if(argc < 3) { 285 fprintf(stderr, "Usage: %s <netcdf> <ref-data-path>\n", argv[0]); 286 return -1; 287 } 288 289 filename = argv[1]; 290 path = argv[2]; 291 292 CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK); 293 CHK(htmie_create(NULL, &allocator, 1, &htmie) == RES_OK); 294 295 CHK(htmie_load(NULL, NULL) == RES_BAD_ARG); 296 CHK(htmie_load(htmie, NULL) == RES_BAD_ARG); 297 CHK(htmie_load(NULL, filename) == RES_BAD_ARG); 298 CHK(htmie_load(htmie, filename) == RES_OK); 299 300 p = strrchr(filename, '/'); 301 base = p ? p+1 : filename; 302 p = strrchr(base, '.'); 303 if(p) *p = '\0'; 304 305 /* Check the wavelengths list */ 306 CHK((size_t)snprintf(buf, sizeof(buf), "%s/%s.lambda", path, base)<sizeof(buf)); 307 CHK(fp = fopen(buf, "r")); 308 FOR_EACH(i, 0, htmie_get_wavelengths_count(htmie)) { 309 double lambda; 310 CHK(fscanf(fp, "%lg", &lambda) == 1); 311 CHK(eq_eps(lambda, htmie_get_wavelengths(htmie)[i], 1.e-6)); 312 } 313 CHK(fscanf(fp, "%*g") == EOF); 314 CHK(fclose(fp) == 0); 315 316 /* Check absorption cross sections */ 317 CHK((size_t)snprintf(buf, sizeof(buf), "%s/%s.macs", path, base)<sizeof(buf)); 318 CHK(fp = fopen(buf, "r")); 319 FOR_EACH(i, 0, htmie_get_wavelengths_count(htmie)) { 320 const double Cabs = htmie_get_xsections_absorption(htmie)[i]; 321 double macs; 322 CHK(fscanf(fp, "%lg", &macs) == 1); 323 CHK(eq_eps(macs, Cabs, 1.e-6)); 324 } 325 CHK(fscanf(fp, "%*g") == EOF); 326 CHK(fclose(fp) == 0); 327 328 /* Check scattering cross sections */ 329 CHK((size_t)snprintf(buf, sizeof(buf), "%s/%s.mscs", path, base)<sizeof(buf)); 330 CHK(fp = fopen(buf, "r")); 331 FOR_EACH(i, 0, htmie_get_wavelengths_count(htmie)) { 332 const double Csca = htmie_get_xsections_scattering(htmie)[i]; 333 double mscs; 334 CHK(fscanf(fp, "%lg", &mscs) == 1); 335 CHK(eq_eps(mscs, Csca, 1.e-6)); 336 } 337 CHK(fscanf(fp, "%*g") == EOF); 338 CHK(fclose(fp) == 0); 339 340 /* Check scattering asymmetry parameter */ 341 CHK((size_t)snprintf(buf, sizeof(buf), "%s/%s.g", path, base)<sizeof(buf)); 342 CHK(fp = fopen(buf, "r")); 343 FOR_EACH(i, 0, htmie_get_wavelengths_count(htmie)) { 344 double g; 345 CHK(fscanf(fp, "%lg", &g) == 1); 346 CHK(eq_eps(g, htmie_get_asymmetry_parameter(htmie)[i], 1.e-6)); 347 } 348 CHK(fscanf(fp, "%*g") == EOF); 349 CHK(fclose(fp) == 0); 350 351 test_fetch(htmie); 352 test_bounds(htmie); 353 test_avg(htmie); 354 355 CHK(htmie_ref_put(htmie) == RES_OK); 356 357 check_memory_allocator(&allocator); 358 mem_shutdown_proxy_allocator(&allocator); 359 CHK(mem_allocated_size() == 0); 360 return 0; 361 } 362