htmie.c (21261B)
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 #include "htmie.h" 19 20 #include <rsys/algorithm.h> 21 #include <rsys/dynamic_array_double.h> 22 #include <rsys/logger.h> 23 #include <rsys/mem_allocator.h> 24 #include <rsys/ref_count.h> 25 26 #include <netcdf.h> 27 #include <stdlib.h> 28 29 #define INVALID_NC_ID -1 30 31 #define NC(Func) { \ 32 const int err__ = nc_ ## Func; \ 33 if(err__ != NC_NOERR) { \ 34 log_err(htmie, "error:%i:%s\n", __LINE__, ncerr_to_str(err__)); \ 35 abort(); \ 36 } \ 37 } (void)0 38 39 /* Context used in the sum_xsections functor */ 40 struct sum_context { 41 double prev_data; 42 double prev_wavelength; 43 double sum; 44 int first_entry; 45 }; 46 static const struct sum_context SUM_CONTEXT_NULL = 47 {DBL_MAX, DBL_MAX, 0, 1}; 48 49 struct htmie { 50 struct darray_double wavelengths; /* In nano-meters */ 51 struct darray_double xsections_absorption; /* In m^2.kg^-1 */ 52 struct darray_double xsections_scattering; /* In m^2.kg^-1 */ 53 struct darray_double g; /* Asymetric parameter */ 54 55 int verbose; /* Verbosity level */ 56 struct logger* logger; 57 struct mem_allocator* allocator; 58 ref_T ref; 59 }; 60 61 /******************************************************************************* 62 * Helper functions 63 ******************************************************************************/ 64 static INLINE void 65 log_msg 66 (const struct htmie* htmie, 67 const enum log_type stream, 68 const char* msg, 69 va_list vargs) 70 { 71 ASSERT(htmie && msg); 72 if(htmie->verbose) { 73 res_T res; (void)res; 74 res = logger_vprint(htmie->logger, stream, msg, vargs); 75 ASSERT(res == RES_OK); 76 } 77 } 78 79 static INLINE void 80 log_err(const struct htmie* htmie, const char* msg, ...) 81 { 82 va_list vargs_list; 83 ASSERT(htmie && msg); 84 va_start(vargs_list, msg); 85 log_msg(htmie, LOG_ERROR, msg, vargs_list); 86 va_end(vargs_list); 87 } 88 89 static INLINE void 90 log_warn(const struct htmie* htmie, const char* msg, ...) 91 { 92 va_list vargs_list; 93 ASSERT(htmie && msg); 94 va_start(vargs_list, msg); 95 log_msg(htmie, LOG_WARNING, msg, vargs_list); 96 va_end(vargs_list); 97 } 98 99 static INLINE int 100 cmp_dbl(const void* a, const void* b) 101 { 102 const double d0 = *((const double*)a); 103 const double d1 = *((const double*)b); 104 return d0 < d1 ? -1 : (d0 > d1 ? 1 : 0); 105 } 106 107 static INLINE const char* 108 ncerr_to_str(const int ncerror) 109 { 110 const char* str = "NC_ERR_<UNKNOWN>"; 111 switch(ncerror) { 112 case NC_EBADGRPID: str = "NC_EBADGRPID"; break; 113 case NC_EBADID: str = "NC_EBADID"; break; 114 case NC_EBADNAME: str = "NC_EBADNAME"; break; 115 case NC_ECHAR: str = "NC_ECHAR"; break; 116 case NC_EDIMMETA: str = "NC_EDIMMETA"; break; 117 case NC_EHDFERR: str = "NC_EHDFERR"; break; 118 case NC_ENOMEM: str = "NC_ENOMEM"; break; 119 case NC_ENOTATT: str = "NC_ENOTATT"; break; 120 case NC_ENOTVAR: str = "NC_ENOTVAR"; break; 121 case NC_ERANGE: str = "NC_ERANGE"; break; 122 case NC_NOERR: str = "NC_NOERR"; break; 123 } 124 return str; 125 } 126 127 static INLINE const char* 128 nctype_to_str(const nc_type type) 129 { 130 const char* str = "NC_TYPE_<UNKNOWN>"; 131 switch(type) { 132 case NC_NAT: str = "NC_NAT"; break; 133 case NC_BYTE: str = "NC_BYTE"; break; 134 case NC_CHAR: str = "NC_CHAR"; break; 135 case NC_SHORT: str = "NC_SHORT"; break; 136 case NC_LONG: str = "NC_LONG"; break; 137 case NC_FLOAT: str = "NC_FLOAT"; break; 138 case NC_DOUBLE: str = "NC_DOUBLE"; break; 139 case NC_UBYTE: str = "NC_UBYTE"; break; 140 case NC_USHORT: str = "NC_USHORT"; break; 141 case NC_UINT: str = "NC_UINT"; break; 142 case NC_INT64: str = "NC_INT64"; break; 143 case NC_UINT64: str = "NC_UINT64"; break; 144 case NC_STRING: str = "NC_STRING"; break; 145 default: FATAL("Unreachable code.\n"); break; 146 } 147 return str; 148 } 149 150 static INLINE size_t 151 sizeof_nctype(const nc_type type) 152 { 153 size_t sz; 154 switch(type) { 155 case NC_BYTE: 156 case NC_CHAR: 157 case NC_UBYTE: 158 sz = 1; break; 159 case NC_SHORT: 160 case NC_USHORT: 161 sz = 2; break; 162 case NC_FLOAT: 163 case NC_INT: 164 case NC_UINT: 165 sz = 4; break; 166 case NC_DOUBLE: 167 case NC_INT64: 168 case NC_UINT64: 169 sz = 8; break; 170 default: FATAL("Unreachable cde.\n"); break; 171 } 172 return sz; 173 } 174 175 static res_T 176 load_wavelengths(struct htmie* htmie, const int nc, const double range[2]) 177 { 178 size_t len; 179 size_t start; 180 size_t i; 181 int id; 182 int ndims; 183 int dimid; 184 int err= NC_NOERR; 185 int type; 186 res_T res = RES_OK; 187 ASSERT(htmie && nc != INVALID_NC_ID && range && range[0] <= range[1]); 188 189 err = nc_inq_varid(nc, "lambda", &id); 190 if(err != NC_NOERR) { 191 log_err(htmie, "Could not inquire the 'lambda' variable -- %s\n", 192 ncerr_to_str(err)); 193 res = RES_BAD_ARG; 194 goto error; 195 } 196 197 NC(inq_varndims(nc, id, &ndims)); 198 if(ndims != 1) { 199 log_err(htmie, "The dimension of the 'lambda' variable must be 1.\n"); 200 res = RES_BAD_ARG; 201 goto error; 202 } 203 204 NC(inq_vardimid(nc, id, &dimid)); 205 NC(inq_dimlen(nc, dimid, &len)); 206 NC(inq_vartype(nc, id, &type)); 207 208 if(type != NC_DOUBLE && type != NC_FLOAT) { 209 log_err(htmie, 210 "The type of the 'lambda' variable cannot be %s. " 211 "Expecting floating point data.\n", 212 nctype_to_str(type)); 213 res = RES_BAD_ARG; 214 goto error; 215 } 216 217 res = darray_double_resize(&htmie->wavelengths, len); 218 if(res != RES_OK) { 219 log_err(htmie, "Could not allocate the list of wavelengths.\n"); 220 goto error; 221 } 222 223 start = 0; 224 NC(get_vara_double(nc, id, &start, &len, 225 darray_double_data_get(&htmie->wavelengths))); 226 227 /* Check data validity */ 228 FOR_EACH(i, 0, len) { 229 if(darray_double_cdata_get(&htmie->wavelengths)[i] < range[0] 230 || darray_double_cdata_get(&htmie->wavelengths)[i] > range[1]) { 231 log_err(htmie, "Invalid wavelength %g. It must be in [%g, %g]\n", 232 darray_double_cdata_get(&htmie->wavelengths)[i], 233 range[0], range[1]); 234 res = RES_BAD_ARG; 235 goto error; 236 } 237 } 238 239 FOR_EACH(i, 1, len) { 240 if(darray_double_cdata_get(&htmie->wavelengths)[i-1] 241 > darray_double_cdata_get(&htmie->wavelengths)[i-0]) { 242 break; 243 } 244 } 245 246 if(i < len) { /* Wavelengths are not sorted */ 247 log_warn(htmie, "Unordered wavelengths.\n"); 248 qsort 249 (darray_double_data_get(&htmie->wavelengths), 250 darray_double_size_get(&htmie->wavelengths), 251 sizeof(double), 252 cmp_dbl); 253 } 254 255 exit: 256 return res; 257 error: 258 darray_double_clear(&htmie->wavelengths); 259 goto exit; 260 } 261 262 static res_T 263 load_distrib_x_lambda_array 264 (struct htmie* htmie, 265 const int nc, 266 const char* var_name, 267 const double range[2], /* Validity range */ 268 struct darray_double* values) 269 { 270 size_t start[2]; 271 size_t end[2]; 272 size_t len; 273 size_t i; 274 int dimids[2]; 275 int id; 276 int ndims; 277 int err; 278 int type; 279 res_T res = RES_OK; 280 ASSERT(htmie && nc != INVALID_NC_ID && var_name && values && range); 281 ASSERT(range[0] <= range[1]); 282 283 err = nc_inq_varid(nc, var_name, &id); 284 if(err != NC_NOERR) { 285 log_err(htmie, "Could not inquire the '%s' variable -- %s\n", 286 var_name, ncerr_to_str(err)); 287 res = RES_BAD_ARG; 288 goto error; 289 } 290 291 NC(inq_varndims(nc, id, &ndims)); 292 if(ndims != 2) { 293 log_err(htmie, 294 "The dimension of the '%s' variable must be 2.\n", var_name); 295 res = RES_BAD_ARG; 296 goto error; 297 } 298 299 NC(inq_vardimid(nc, id, dimids)); 300 NC(inq_vartype(nc, id, &type)); 301 302 NC(inq_dimlen(nc, dimids[0], &len)); 303 if(len != 1) { 304 log_err(htmie, 305 "Only 1 distribution is currently supported while the #distributions of " 306 "the '%s' variable is %lu\n", var_name, (unsigned long)len); 307 res = RES_BAD_ARG; 308 goto error; 309 } 310 311 NC(inq_dimlen(nc, dimids[1], &len)); 312 if(len != darray_double_size_get(&htmie->wavelengths)) { 313 log_err(htmie, 314 "Inconsistent wavelengths count. The number of defined wavelengths is %lu " 315 "while the per distribution length of the '%s' variable is %lu.\n", 316 darray_double_size_get(&htmie->wavelengths), var_name, len); 317 res = RES_BAD_ARG; 318 goto error; 319 } 320 321 if(type != NC_DOUBLE && type != NC_FLOAT) { 322 log_err(htmie, 323 "The type of the '%s' variable cannot be %s. " 324 "Expecting floating point data.\n", 325 var_name, nctype_to_str(type)); 326 res = RES_BAD_ARG; 327 goto error; 328 } 329 330 res = darray_double_resize(values, len); 331 if(res != RES_OK) { 332 log_err(htmie, 333 "Could not allocate the memory space to load the data of the '%s' variable.\n", 334 var_name); 335 res = RES_BAD_ARG; 336 goto error; 337 } 338 339 /* Read the raw data sections */ 340 start[0] = 0, start[1] = 0; 341 end[0] = 1, end[1] = len; 342 NC(get_vara_double(nc, id, start, end, darray_double_data_get(values))); 343 344 FOR_EACH(i, 0, darray_double_size_get(values)) { 345 const double* d = darray_double_cdata_get(values); 346 if(d[i] < range[0] || d[i] > range[1]) { 347 log_err(htmie, 348 "Invalid value for the %lu^th data of the '%s' variable : %g. " 349 "The data must be in [%g, %g]\n", 350 (unsigned long)i, var_name, d[i], range[0], range[1]); 351 res = RES_BAD_ARG; 352 goto error; 353 } 354 } 355 356 exit: 357 return res; 358 error: 359 darray_double_clear(values); 360 goto exit; 361 } 362 363 static FINLINE double 364 interpolate_data 365 (const struct htmie* htmie, 366 const double wavelength, /* Input wavelength */ 367 const enum htmie_filter_type type, /* Interpolation type */ 368 const size_t idata, /* Id of the upper bound data wrt `wavelength' */ 369 const double* data) /* Data to interpolate */ 370 { 371 const double* wlens; 372 double a, b, u; 373 double val; 374 ASSERT(data); 375 376 wlens = htmie_get_wavelengths(htmie); 377 378 ASSERT(idata < htmie_get_wavelengths_count(htmie)); 379 ASSERT(wlens[idata-1] < wavelength && wavelength <= wlens[idata-0]); 380 381 a = data[idata-1]; 382 b = data[idata-0]; 383 u = (wavelength - wlens[idata-1]) / (wlens[idata] - wlens[idata-1]); 384 ASSERT(eq_eps(u, 1, 1.e-6) || u < 1); 385 ASSERT(eq_eps(u, 0, 1.e-6) || u > 0); 386 387 switch(type) { 388 case HTMIE_FILTER_NEAREST: 389 val = u < 0.5 ? a : b; 390 break; 391 case HTMIE_FILTER_LINEAR: 392 u = CLAMP(u, 0, 1); /* Handle numerical inaccuracy */ 393 val = u*b + (1-u)*a; 394 break; 395 default: FATAL("Unreachable code.\n"); break; 396 } 397 return val; 398 } 399 400 static FINLINE double 401 fetch_data 402 (const struct htmie* htmie, 403 const double wavelength, 404 const enum htmie_filter_type type, 405 const double* data) /* Data to interpolate */ 406 { 407 size_t nwlens; 408 const double* wlens; 409 const double* wlen_upp; 410 double val; 411 ASSERT(htmie && (unsigned)type < HTMIE_FILTER_TYPES_COUNT__ && data); 412 413 wlens = htmie_get_wavelengths(htmie); 414 nwlens = htmie_get_wavelengths_count(htmie); 415 ASSERT(nwlens); 416 417 wlen_upp = search_lower_bound 418 (&wavelength, wlens, nwlens, sizeof(double), cmp_dbl); 419 420 if(!wlen_upp) { /* Clamp to upper */ 421 val = data[nwlens-1]; 422 } else if(wlen_upp == wlens) { /* Clamp to lower */ 423 val = data[0]; 424 } else { 425 const size_t i = (size_t)(wlen_upp - wlens); 426 val = interpolate_data(htmie, wavelength, type, i, data); 427 } 428 return val; 429 } 430 431 static FINLINE void 432 min_max(const double wavelength, const double data, void* ctx) 433 { 434 double* bounds = ctx; 435 ASSERT(ctx); 436 (void)wavelength; 437 bounds[0] = MMIN(bounds[0], data); 438 bounds[1] = MMAX(bounds[1], data); 439 } 440 441 static FINLINE void 442 sum(const double wavelength, const double data, void* context) 443 { 444 struct sum_context* ctx = context; 445 ASSERT(context); 446 447 if(ctx->first_entry) { 448 ASSERT(ctx->sum == 0); 449 ctx->first_entry = 0; 450 } else { 451 ASSERT(wavelength > ctx->prev_wavelength); 452 ctx->sum += 453 0.5 * (ctx->prev_data + data)*(wavelength - ctx->prev_wavelength); 454 } 455 ctx->prev_wavelength = wavelength; 456 ctx->prev_data = data; 457 } 458 459 static INLINE void 460 for_each_wavelength_in_spectral_band 461 (const struct htmie* htmie, 462 const double band[2], /* Boundaries of the spectral band in nanometer */ 463 const enum htmie_filter_type type, 464 const double* data, /* Input data of the spectral band*/ 465 void (*op)(const double wavelength, const double data, void* ctx), 466 void* context) /* Sent as the last argument of the 'op' functor */ 467 { 468 const double* wlens; 469 const double* wlen_low; 470 const double* wlen_upp; 471 double data_low; 472 double data_upp; 473 size_t ilow; 474 size_t iupp; 475 size_t nwlens; 476 size_t i; 477 ASSERT(htmie && band[0] <= band[1]); 478 479 wlens = htmie_get_wavelengths(htmie); 480 nwlens = htmie_get_wavelengths_count(htmie); 481 ASSERT(nwlens); 482 483 wlen_low = search_lower_bound 484 (&band[0], wlens, nwlens, sizeof(double), cmp_dbl); 485 wlen_upp = search_lower_bound 486 (&band[1], wlens, nwlens, sizeof(double), cmp_dbl); 487 488 if(!wlen_low) { 489 ilow = nwlens; 490 data_low = data[nwlens - 1]; 491 } else if(wlen_low == wlens) { 492 ilow = 1; 493 data_low = data[0]; 494 } else { 495 ilow = (size_t)(wlen_low - wlens); 496 data_low = interpolate_data(htmie, band[0], type, ilow, data); 497 } 498 499 if(!wlen_upp) { 500 iupp = nwlens; 501 data_upp = data[nwlens - 1]; 502 } else if(wlen_upp == wlens) { 503 iupp = 1; 504 data_upp = data[0]; 505 } else { 506 iupp = (size_t)(wlen_upp - wlens); 507 data_upp = interpolate_data(htmie, band[1], type, iupp, data); 508 } 509 510 if(wlens[ilow] == band[0]) ++ilow; 511 512 op(band[0], data_low, context); 513 FOR_EACH(i, ilow, iupp) { 514 op(wlens[i], data[i], context); 515 } 516 if(band[0] != band[1]) 517 op(band[1], data_upp, context); 518 } 519 520 static void 521 release_htmie(ref_T* ref) 522 { 523 struct htmie* htmie; 524 ASSERT(ref); 525 htmie = CONTAINER_OF(ref, struct htmie, ref); 526 darray_double_release(&htmie->wavelengths); 527 darray_double_release(&htmie->xsections_absorption); 528 darray_double_release(&htmie->xsections_scattering); 529 darray_double_release(&htmie->g); 530 MEM_RM(htmie->allocator, htmie); 531 } 532 533 /******************************************************************************* 534 * Exported functions 535 ******************************************************************************/ 536 res_T 537 htmie_create 538 (struct logger* log, 539 struct mem_allocator* mem_allocator, 540 const int verbose, 541 struct htmie** out_htmie) 542 { 543 struct htmie* htmie = NULL; 544 struct mem_allocator* allocator = NULL; 545 struct logger* logger = NULL; 546 res_T res = RES_OK; 547 548 if(!out_htmie) { 549 res = RES_BAD_ARG; 550 goto error; 551 } 552 553 allocator = mem_allocator ? mem_allocator : &mem_default_allocator; 554 logger = log ? log : LOGGER_DEFAULT; 555 556 htmie = MEM_CALLOC(allocator, 1, sizeof(struct htmie)); 557 if(!htmie) { 558 if(verbose) { 559 logger_print(logger, LOG_ERROR, 560 "%s: could not allocate the HTMIE handler.\n", FUNC_NAME); 561 } 562 res = RES_MEM_ERR; 563 goto error; 564 } 565 ref_init(&htmie->ref); 566 htmie->allocator = allocator; 567 htmie->logger = logger; 568 htmie->verbose = verbose; 569 darray_double_init(htmie->allocator, &htmie->wavelengths); 570 darray_double_init(htmie->allocator, &htmie->xsections_absorption); 571 darray_double_init(htmie->allocator, &htmie->xsections_scattering); 572 darray_double_init(htmie->allocator, &htmie->g); 573 574 exit: 575 if(out_htmie) *out_htmie = htmie; 576 return res; 577 error: 578 if(htmie) { 579 HTMIE(ref_put(htmie)); 580 htmie = NULL; 581 } 582 goto exit; 583 } 584 585 res_T 586 htmie_ref_get(struct htmie* htmie) 587 { 588 if(!htmie) return RES_BAD_ARG; 589 ref_get(&htmie->ref); 590 return RES_OK; 591 } 592 593 res_T 594 htmie_ref_put(struct htmie* htmie) 595 { 596 if(!htmie) return RES_BAD_ARG; 597 ref_put(&htmie->ref, release_htmie); 598 return RES_OK; 599 } 600 601 res_T 602 htmie_load(struct htmie* htmie, const char* path) 603 { 604 int err_nc = NC_NOERR; 605 int nc = INVALID_NC_ID; 606 double range[2] = {DBL_MAX, -DBL_MAX}; /* Validity range of loaded data */ 607 res_T res = RES_OK; 608 609 if(!htmie || !path) { 610 res = RES_BAD_ARG; 611 goto error; 612 } 613 614 err_nc = nc_open(path, NC_NOWRITE, &nc); 615 if(err_nc != NC_NOERR) { 616 log_err(htmie, "error opening file `%s' -- %s.\n", path, ncerr_to_str(err_nc)); 617 res = RES_BAD_ARG; 618 goto error; 619 } 620 621 #define CALL(Func) { res = Func; if(res != RES_OK) goto error; } (void)0 622 623 range[0] = DBL_EPSILON; range[1] = DBL_MAX; 624 CALL(load_wavelengths(htmie, nc, range)); 625 626 range[0] = 0; range[1] = DBL_MAX; 627 CALL(load_distrib_x_lambda_array 628 (htmie, nc, "macs", range, &htmie->xsections_absorption)); 629 630 range[0] = 0; range[1] = DBL_MAX; 631 CALL(load_distrib_x_lambda_array 632 (htmie, nc, "mscs", range, &htmie->xsections_scattering)); 633 634 range[0] = -1; range[1] = 1; 635 CALL(load_distrib_x_lambda_array(htmie, nc, "g", range, &htmie->g)); 636 637 #undef CALL 638 639 exit: 640 if(nc != INVALID_NC_ID) NC(close(nc)); 641 return res; 642 error: 643 goto exit; 644 } 645 646 const double* 647 htmie_get_wavelengths(const struct htmie* htmie) 648 { 649 ASSERT(htmie); 650 return darray_double_cdata_get(&htmie->wavelengths); 651 } 652 653 const double* 654 htmie_get_xsections_absorption(const struct htmie* htmie) 655 { 656 ASSERT(htmie); 657 return darray_double_cdata_get(&htmie->xsections_absorption); 658 } 659 660 const double* 661 htmie_get_xsections_scattering(const struct htmie* htmie) 662 { 663 ASSERT(htmie); 664 return darray_double_cdata_get(&htmie->xsections_scattering); 665 } 666 667 const double* 668 htmie_get_asymmetry_parameter(const struct htmie* htmie) 669 { 670 ASSERT(htmie); 671 return darray_double_cdata_get(&htmie->g); 672 } 673 674 size_t 675 htmie_get_wavelengths_count(const struct htmie* htmie) 676 { 677 ASSERT(htmie); 678 return darray_double_size_get(&htmie->wavelengths); 679 } 680 681 double 682 htmie_fetch_xsection_absorption 683 (const struct htmie* htmie, 684 const double wavelength, 685 const enum htmie_filter_type type) 686 { 687 ASSERT(htmie); 688 return fetch_data 689 (htmie, wavelength, type, htmie_get_xsections_absorption(htmie)); 690 } 691 692 double 693 htmie_fetch_xsection_scattering 694 (const struct htmie* htmie, 695 const double wavelength, 696 const enum htmie_filter_type type) 697 { 698 ASSERT(htmie); 699 return fetch_data 700 (htmie, wavelength, type, htmie_get_xsections_scattering(htmie)); 701 } 702 703 double 704 htmie_fetch_asymmetry_parameter 705 (const struct htmie* htmie, 706 const double wavelength, 707 const enum htmie_filter_type type) 708 { 709 ASSERT(htmie); 710 return fetch_data 711 (htmie, wavelength, type, htmie_get_asymmetry_parameter(htmie)); 712 } 713 714 void 715 htmie_compute_xsection_absorption_bounds 716 (const struct htmie* htmie, 717 const double band[2], /* Boundaries of the spectral band in nanometer */ 718 const enum htmie_filter_type type, 719 double bounds[2]) /* Min and Max scattering cross sections in m^2.kg^-1 */ 720 { 721 const double* xsections; 722 ASSERT(htmie); 723 bounds[0] = DBL_MAX; 724 bounds[1] =-DBL_MAX; 725 xsections = darray_double_cdata_get(&htmie->xsections_absorption); 726 for_each_wavelength_in_spectral_band 727 (htmie, band, type, xsections, min_max, bounds); 728 } 729 730 void 731 htmie_compute_xsection_scattering_bounds 732 (const struct htmie* htmie, 733 const double band[2], /* Boundaries of the spectral band in nanometer */ 734 const enum htmie_filter_type type, 735 double bounds[2]) /* Min and Max scattering cross sections in m^2.kg^-1 */ 736 { 737 const double* xsections; 738 ASSERT(htmie); 739 bounds[0] = DBL_MAX; 740 bounds[1] =-DBL_MAX; 741 xsections = darray_double_cdata_get(&htmie->xsections_scattering); 742 for_each_wavelength_in_spectral_band 743 (htmie, band, type, xsections, min_max, bounds); 744 } 745 746 double 747 htmie_compute_xsection_absorption_average 748 (const struct htmie* htmie, 749 const double band[2], /* Boundaries of of the spectral band in nanometer */ 750 const enum htmie_filter_type type) 751 { 752 const double* xsections; 753 struct sum_context ctx = SUM_CONTEXT_NULL; 754 ASSERT(htmie && band && band[0] < band[1]); 755 xsections = darray_double_cdata_get(&htmie->xsections_absorption); 756 for_each_wavelength_in_spectral_band(htmie, band, type, xsections, sum, &ctx); 757 return ctx.sum / (band[1] - band[0]); 758 } 759 760 double 761 htmie_compute_xsection_scattering_average 762 (const struct htmie* htmie, 763 const double band[2], /* Boundaries of of the spectral band in nanometer */ 764 const enum htmie_filter_type type) 765 { 766 const double* xsections; 767 struct sum_context ctx = SUM_CONTEXT_NULL; 768 ASSERT(htmie && band && band[0] < band[1]); 769 xsections = darray_double_cdata_get(&htmie->xsections_scattering); 770 for_each_wavelength_in_spectral_band(htmie, band, type, xsections, sum, &ctx); 771 return ctx.sum / (band[1] - band[0]); 772 } 773 774 double 775 htmie_compute_asymmetry_parameter_average 776 (const struct htmie* htmie, 777 const double band[2], /* Boundaries of of the spectral band in nanometer */ 778 const enum htmie_filter_type type) 779 { 780 const double* g; 781 struct sum_context ctx = SUM_CONTEXT_NULL; 782 ASSERT(htmie && band && band[0] < band[1]); 783 g = darray_double_cdata_get(&htmie->g); 784 for_each_wavelength_in_spectral_band(htmie, band, type, g, sum, &ctx); 785 return ctx.sum / (band[1] - band[0]); 786 } 787