rnatm_properties.c (37145B)
1 /* Copyright (C) 2022, 2023, 2025 Centre National de la Recherche Scientifique 2 * Copyright (C) 2022, 2023, 2025 Institut Pierre-Simon Laplace 3 * Copyright (C) 2022, 2023, 2025 Institut de Physique du Globe de Paris 4 * Copyright (C) 2022, 2023, 2025 |Méso|Star> (contact@meso-star.com) 5 * Copyright (C) 2022, 2023, 2025 Observatoire de Paris 6 * Copyright (C) 2022, 2023, 2025 Université de Reims Champagne-Ardenne 7 * Copyright (C) 2022, 2023, 2025 Université de Versaille Saint-Quentin 8 * Copyright (C) 2022, 2023, 2025 Université Paul Sabatier 9 * 10 * This program is free software: you can redistribute it and/or modify 11 * it under the terms of the GNU General Public License as published by 12 * the Free Software Foundation, either version 3 of the License, or 13 * (at your option) any later version. 14 * 15 * This program is distributed in the hope that it will be useful, 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 * GNU General Public License for more details. 19 * 20 * You should have received a copy of the GNU General Public License 21 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 22 23 #include "rnatm.h" 24 #include "rnatm_c.h" 25 #include "rnatm_log.h" 26 27 #include <rad-net/rnsf.h> 28 #include <rad-net/rnsl.h> 29 30 #include <star/sars.h> 31 #include <star/sbuf.h> 32 #include <star/sck.h> 33 #include <star/ssf.h> 34 #include <star/suvm.h> 35 36 #include <rsys/clock_time.h> 37 #include <rsys/cstr.h> 38 39 /******************************************************************************* 40 * Helper functions 41 ******************************************************************************/ 42 static res_T 43 check_rnatm_get_radcoef_args 44 (const struct rnatm* atm, 45 const struct rnatm_get_radcoef_args* args) 46 { 47 size_t i, n; 48 49 if(!args) return RES_BAD_ARG; 50 51 /* Invalid cells */ 52 if(!args->cells) return RES_BAD_ARG; 53 54 /* Invalid band index */ 55 n = darray_band_size_get(&atm->bands) - 1; 56 if(args->iband < darray_band_cdata_get(&atm->bands)[0].index 57 || args->iband > darray_band_cdata_get(&atm->bands)[n].index) 58 return RES_BAD_ARG; 59 60 /* Invalid quadrature point index */ 61 i = args->iband - darray_band_cdata_get(&atm->bands)[0].index; 62 ASSERT(i <= n && args->iband == darray_band_cdata_get(&atm->bands)[i].index); 63 if(args->iquad >= darray_band_cdata_get(&atm->bands)[i].nquad_pts) 64 return RES_BAD_ARG; 65 66 /* Invalid radiative coefficient */ 67 if((unsigned)args->radcoef >= RNATM_RADCOEFS_COUNT__) 68 return RES_BAD_ARG; 69 70 /* Invalid K range */ 71 if(args->k_min > args->k_max) 72 return RES_BAD_ARG; 73 74 return RES_OK; 75 } 76 77 static res_T 78 check_rnatm_sample_component_args 79 (const struct rnatm* atm, 80 const struct rnatm_sample_component_args* args) 81 { 82 size_t i, n; 83 84 if(!args) return RES_BAD_ARG; 85 86 /* Invalid cells */ 87 if(!args->cells) return RES_BAD_ARG; 88 89 /* Invalid band index */ 90 n = darray_band_size_get(&atm->bands) - 1; 91 if(args->iband < darray_band_cdata_get(&atm->bands)[0].index 92 || args->iband > darray_band_cdata_get(&atm->bands)[n].index) 93 return RES_BAD_ARG; 94 95 /* Invalid quadrature point index */ 96 i = args->iband - darray_band_cdata_get(&atm->bands)[0].index; 97 ASSERT(i <= n && args->iband == darray_band_cdata_get(&atm->bands)[i].index); 98 if(args->iquad >= darray_band_cdata_get(&atm->bands)[i].nquad_pts) 99 return RES_BAD_ARG; 100 101 /* Invalid radiative coefficient */ 102 if((unsigned)args->radcoef >= RNATM_RADCOEFS_COUNT__) 103 return RES_BAD_ARG; 104 105 /* Invalid random number */ 106 if(args->r < 0 || args->r >= 1) 107 return RES_BAD_ARG; 108 109 return RES_OK; 110 } 111 112 static INLINE res_T 113 check_cell(const struct rnatm_cell_pos* cell) 114 { 115 double sum_bcoords; 116 117 if(!cell) return RES_BAD_ARG; 118 119 /* Invalid geometric primitive */ 120 if(SUVM_PRIMITIVE_NONE(&cell->prim) 121 || cell->prim.nvertices != 4) /* Expect a tetraheron */ 122 return RES_BAD_ARG; 123 124 /* Invalid barycentric coordinates */ 125 sum_bcoords = 126 cell->barycentric_coords[0] 127 + cell->barycentric_coords[1] 128 + cell->barycentric_coords[2] 129 + cell->barycentric_coords[3]; 130 if(!eq_eps(sum_bcoords, 1, 1.e-6)) 131 return RES_BAD_ARG; 132 133 return RES_OK; 134 } 135 136 static res_T 137 check_rnatm_cell_get_radcoef_args 138 (const struct rnatm* atm, 139 const struct rnatm_cell_get_radcoef_args* args) 140 { 141 size_t i, n; 142 res_T res = RES_OK; 143 ASSERT(atm && darray_band_size_get(&atm->bands)); 144 145 if(!args) return RES_BAD_ARG; 146 147 /* Invalid cell */ 148 res = check_cell(&args->cell); 149 if(res != RES_OK) return res; 150 151 /* Invalid component */ 152 if(args->cell.component != RNATM_GAS 153 && args->cell.component >= darray_aerosol_size_get(&atm->aerosols)) 154 return RES_BAD_ARG; 155 156 /* Invalid band index */ 157 n = darray_band_size_get(&atm->bands) - 1; 158 if(args->iband < darray_band_cdata_get(&atm->bands)[0].index 159 || args->iband > darray_band_cdata_get(&atm->bands)[n].index) 160 return RES_BAD_ARG; 161 162 /* Invalid quadrature point index */ 163 i = args->iband - darray_band_cdata_get(&atm->bands)[0].index; 164 ASSERT(i <= n && args->iband == darray_band_cdata_get(&atm->bands)[i].index); 165 if(args->iquad >= darray_band_cdata_get(&atm->bands)[i].nquad_pts) 166 return RES_BAD_ARG; 167 168 /* Invalid radiative coefficient */ 169 if((unsigned)args->radcoef >= RNATM_RADCOEFS_COUNT__) 170 return RES_BAD_ARG; 171 172 /* Invalid K range */ 173 if(args->k_max < 0) 174 return RES_BAD_ARG; 175 176 return RES_OK; 177 } 178 179 static res_T 180 check_rnatm_cell_create_phase_fn_args 181 (struct rnatm* atm, 182 const struct rnatm_cell_create_phase_fn_args* args) 183 { 184 res_T res = RES_OK; 185 ASSERT(atm); 186 187 if(!args) return RES_BAD_ARG; 188 189 /* Invalid cell */ 190 res = check_cell(&args->cell); 191 if(res != RES_OK) return res; 192 193 /* Invalid component */ 194 if(args->cell.component != RNATM_GAS 195 && args->cell.component >= darray_aerosol_size_get(&atm->aerosols)) 196 return RES_BAD_ARG; 197 198 /* Invalid wavelength */ 199 if(args->wavelength < 0) 200 return RES_BAD_ARG; 201 202 /* Invalid random numbers */ 203 if(args->r[0] < 0 || args->r[0] >= 1 204 || args->r[1] < 0 || args->r[1] >= 1) 205 return RES_BAD_ARG; 206 207 return RES_OK; 208 } 209 210 static INLINE void 211 reset_gas_properties(struct gas* gas) 212 { 213 if(gas->temperatures) SBUF(ref_put(gas->temperatures)); 214 if(gas->ck) SCK(ref_put(gas->ck)); 215 gas->temperatures = NULL; 216 gas->ck = NULL; 217 } 218 219 static INLINE void 220 reset_aerosol_properties(struct aerosol* aerosol) 221 { 222 if(aerosol->phase_fn_ids) SBUF(ref_put(aerosol->phase_fn_ids)); 223 if(aerosol->sars) SARS(ref_put(aerosol->sars)); 224 aerosol->phase_fn_ids = NULL; 225 aerosol->sars = NULL; 226 darray_phase_fn_clear(&aerosol->phase_fn_lst); 227 } 228 229 static INLINE res_T 230 check_gas_temperatures_desc 231 (const struct rnatm* atm, 232 const struct sbuf_desc* desc, 233 const struct rnatm_gas_args* gas_args) 234 { 235 ASSERT(atm && desc && gas_args); 236 237 if(desc->size != atm->gas.nvertices) { 238 log_err(atm, 239 "%s: no sufficient temperatures regarding the mesh %s\n", 240 gas_args->temperatures_filename, gas_args->smsh_filename); 241 return RES_BAD_ARG; 242 } 243 244 if(desc->szitem != 4 || desc->alitem != 4 || desc->pitch != 4) { 245 log_err(atm, "%s: unexpected layout of temperatures\n", 246 gas_args->temperatures_filename); 247 return RES_BAD_ARG; 248 } 249 250 return RES_OK; 251 } 252 253 static res_T 254 check_gas_temperatures(const struct rnatm* atm) 255 { 256 struct sbuf_desc sbuf_desc = SBUF_DESC_NULL; 257 size_t i; 258 res_T res = RES_OK; 259 ASSERT(atm); 260 261 /* The layout of sbuf_desc has already been checked when creating rnatm */ 262 res = sbuf_get_desc(atm->gas.temperatures, &sbuf_desc); 263 if(res != RES_OK) goto error; 264 265 FOR_EACH(i, 0, sbuf_desc.size) { 266 const float temperature = *((const float*)sbuf_desc_at(&sbuf_desc, i)); 267 if(temperature != temperature /* NaN */ || temperature < 0) { 268 log_err(atm, "%s: node %lu: invalid gas temperature `%g'\n", 269 str_cget(&atm->name), i, temperature); 270 res = RES_BAD_ARG; 271 goto error; 272 } 273 } 274 275 exit: 276 return res; 277 error: 278 goto exit; 279 } 280 281 static INLINE res_T 282 check_gas_ck_desc 283 (const struct rnatm* atm, 284 const struct rnatm_gas_args* gas_args) 285 { 286 ASSERT(atm && gas_args); 287 288 if(sck_get_nodes_count(atm->gas.ck) != atm->gas.nvertices) { 289 log_err(atm, 290 "%s: no sufficient correlated-K regarding the mesh %s\n", 291 gas_args->sck_filename, gas_args->smsh_filename); 292 return RES_BAD_ARG; 293 } 294 return RES_OK; 295 } 296 297 static INLINE res_T 298 check_aerosol_phase_fn_ids_desc 299 (const struct rnatm* atm, 300 const struct aerosol* aerosol, 301 const struct sbuf_desc* desc, 302 const struct rnatm_aerosol_args* aerosol_args) 303 { 304 ASSERT(atm && aerosol && desc && aerosol_args); 305 306 if(desc->size != aerosol->nvertices) { 307 log_err(atm, 308 "%s: no sufficient phase function ids regarding the mesh %s\n", 309 aerosol_args->phase_fn_ids_filename, aerosol_args->smsh_filename); 310 return RES_BAD_ARG; 311 } 312 313 if(desc->szitem != 4 || desc->alitem != 4 || desc->pitch != 4) { 314 log_err(atm, "%s: unexpected layout of phase function ids\n", 315 aerosol_args->phase_fn_ids_filename); 316 return RES_BAD_ARG; 317 } 318 319 return RES_OK; 320 } 321 322 static res_T 323 check_aerosol_phase_fn_ids 324 (const struct rnatm* atm, 325 const struct aerosol* aerosol) 326 { 327 struct sbuf_desc sbuf_desc = SBUF_DESC_NULL; 328 size_t i; 329 res_T res = RES_OK; 330 ASSERT(atm && aerosol); 331 332 /* The layout of sbuf_desc has already been checked when creating rnatm */ 333 res = sbuf_get_desc(aerosol->phase_fn_ids, &sbuf_desc); 334 if(res != RES_OK) goto error; 335 336 FOR_EACH(i, 0, sbuf_desc.size) { 337 const uint32_t id = *((uint32_t*)sbuf_desc_at(&sbuf_desc, i)); 338 339 if(id >= darray_phase_fn_size_get(&aerosol->phase_fn_lst)) { 340 log_err(atm, 341 "%s: node %lu: invalid phase function id `%lu'. It must be in [0, %lu[\n", 342 str_cget(&atm->name), (unsigned long)i, (unsigned long)id, 343 (unsigned long)darray_phase_fn_size_get(&aerosol->phase_fn_lst)); 344 res = RES_BAD_ARG; 345 goto error; 346 } 347 } 348 349 exit: 350 return res; 351 error: 352 goto exit; 353 } 354 355 static INLINE res_T 356 check_aerosol_sars_desc 357 (const struct rnatm* atm, 358 const struct aerosol* aerosol, 359 const struct rnatm_aerosol_args* aerosol_args) 360 { 361 ASSERT(atm && aerosol && aerosol_args); 362 363 if(sars_get_nodes_count(aerosol->sars) != aerosol->nvertices) { 364 log_err(atm, 365 "%s: no sufficient radiative properties regarding the mesh %s\n", 366 aerosol_args->sars_filename, aerosol_args->smsh_filename); 367 return RES_BAD_ARG; 368 } 369 return RES_OK; 370 } 371 372 static INLINE res_T 373 create_phase_fn_rayleigh(struct rnatm* atm, struct ssf_phase** phase_fn) 374 { 375 res_T res = RES_OK; 376 ASSERT(atm && phase_fn); 377 378 res = ssf_phase_create(atm->allocator, &ssf_phase_rayleigh, phase_fn); 379 if(res != RES_OK) goto error; 380 381 exit: 382 return res; 383 error: 384 if(*phase_fn) { SSF(phase_ref_put(*phase_fn)); *phase_fn = NULL; } 385 goto exit; 386 } 387 388 static INLINE res_T 389 create_phase_fn_hg 390 (struct rnatm* atm, 391 const struct rnsf_phase_fn_hg* hg, 392 struct ssf_phase** out_ssf_phase) 393 { 394 struct ssf_phase* ssf_phase = NULL; 395 res_T res = RES_OK; 396 ASSERT(atm && hg && out_ssf_phase); 397 398 res = ssf_phase_create(atm->allocator, &ssf_phase_hg, &ssf_phase); 399 if(res != RES_OK) goto error; 400 res = ssf_phase_hg_setup(ssf_phase, hg->g); 401 if(res != RES_OK) goto error; 402 403 exit: 404 *out_ssf_phase = ssf_phase; 405 return res; 406 error: 407 if(ssf_phase) { SSF(phase_ref_put(ssf_phase)); ssf_phase = NULL; } 408 goto exit; 409 } 410 411 static void 412 phase_fn_discrete_get_item 413 (const size_t iitem, 414 struct ssf_discrete_item* item, 415 void* context) 416 { 417 const struct rnsf_phase_fn_discrete* discrete = context; 418 ASSERT(item && context && iitem < discrete->nitems); 419 item->theta = discrete->items[iitem].theta; 420 item->value = discrete->items[iitem].value; 421 } 422 423 static INLINE res_T 424 create_phase_fn_discrete 425 (struct rnatm* atm, 426 struct rnsf_phase_fn_discrete* discrete, 427 struct ssf_phase** out_ssf_phase) 428 { 429 struct ssf_discrete_setup_args args = SSF_DISCRETE_SETUP_ARGS_NULL; 430 struct ssf_phase* ssf_phase = NULL; 431 res_T res = RES_OK; 432 ASSERT(atm && discrete && out_ssf_phase); 433 434 res = ssf_phase_create(atm->allocator, &ssf_phase_discrete, &ssf_phase); 435 if(res != RES_OK) goto error; 436 437 args.get_item = phase_fn_discrete_get_item; 438 args.context = discrete; 439 args.nitems = discrete->nitems; 440 res = ssf_phase_discrete_setup(ssf_phase, &args); 441 if(res != RES_OK) goto error; 442 443 exit: 444 *out_ssf_phase = ssf_phase; 445 return res; 446 error: 447 if(ssf_phase) { SSF(phase_ref_put(ssf_phase)); ssf_phase = NULL; } 448 goto exit; 449 } 450 451 static res_T 452 create_phase_fn 453 (struct rnatm* atm, 454 const struct rnsf_phase_fn* phase_fn, 455 struct ssf_phase** out_ssf_phase) 456 { 457 struct rnsf_phase_fn_hg hg; 458 struct rnsf_phase_fn_discrete discrete; 459 struct ssf_phase* ssf_phase = NULL; 460 res_T res = RES_OK; 461 ASSERT(atm && phase_fn && out_ssf_phase); 462 463 /* Create the scattering function */ 464 switch(rnsf_phase_fn_get_type(phase_fn)) { 465 case RNSF_PHASE_FN_HG: 466 RNSF(phase_fn_get_hg(phase_fn, &hg)); 467 res = create_phase_fn_hg(atm, &hg, &ssf_phase); 468 if(res != RES_OK) goto error; 469 break; 470 case RNSF_PHASE_FN_DISCRETE: 471 RNSF(phase_fn_get_discrete(phase_fn, &discrete)); 472 res = create_phase_fn_discrete(atm, &discrete, &ssf_phase); 473 if(res != RES_OK) goto error; 474 break; 475 default: FATAL("Unreachable code\n"); break; 476 } 477 478 exit: 479 *out_ssf_phase = ssf_phase; 480 return res; 481 error: 482 if(ssf_phase) { SSF(phase_ref_put(ssf_phase)); ssf_phase = NULL; } 483 goto exit; 484 } 485 486 static res_T 487 phase_fn_create_phase_list 488 (struct rnatm* atm, 489 const char* filename, 490 struct phase_fn* phase_fn) 491 { 492 size_t nphases = 0; 493 size_t iphase = 0; 494 res_T res = RES_OK; 495 ASSERT(atm && filename && phase_fn); 496 497 /* Reserve memory space for the list of phase functions */ 498 nphases = rnsf_get_phase_fn_count(phase_fn->rnsf); 499 res = darray_phase_resize(&phase_fn->phase_lst, nphases); 500 if(res != RES_OK) goto error; 501 502 FOR_EACH(iphase, 0, nphases) { 503 const struct rnsf_phase_fn* fn = rnsf_get_phase_fn(phase_fn->rnsf, iphase); 504 struct ssf_phase** phase = darray_phase_data_get(&phase_fn->phase_lst) + iphase; 505 506 res = create_phase_fn(atm, fn, phase); 507 if(res != RES_OK) goto error; 508 } 509 510 exit: 511 return res; 512 513 error: 514 log_err(atm, "%s: error creating the lisf of phase functions -- %s\n", 515 filename, res_to_cstr(res)); 516 darray_phase_clear(&phase_fn->phase_lst); 517 goto exit; 518 } 519 520 static res_T 521 load_phase_fn 522 (struct rnatm* atm, 523 const char* filename, 524 struct rnsf** out_phase_fn) 525 { 526 struct rnsf_create_args args = RNSF_CREATE_ARGS_DEFAULT; 527 struct rnsf* phase_fn = NULL; 528 res_T res = RES_OK; 529 ASSERT(atm && filename && out_phase_fn); 530 531 args.verbose = atm->verbose; 532 args.logger = atm->logger; 533 args.allocator = atm->allocator; 534 res = rnsf_create(&args, &phase_fn); 535 if(res != RES_OK) { 536 log_err(atm, 537 "%s: could not create the Rad-Net Scattering Function data structure\n", 538 filename); 539 goto error; 540 } 541 542 res = rnsf_load(phase_fn, filename); 543 if(res != RES_OK) goto error; 544 545 exit: 546 *out_phase_fn = phase_fn; 547 return res; 548 error: 549 if(phase_fn) { RNSF(ref_put(phase_fn)); phase_fn = NULL; } 550 goto exit; 551 } 552 553 static res_T 554 load_phase_fn_list 555 (struct rnatm* atm, 556 struct aerosol* aerosol, 557 const struct rnatm_aerosol_args* args) 558 { 559 struct rnsl_create_args rnsl_args = RNSL_CREATE_ARGS_DEFAULT; 560 struct rnsl* rnsl = NULL; 561 size_t iphase_fn, nphase_fn; 562 res_T res = RES_OK; 563 564 /* Create loader of scattefing function paths */ 565 rnsl_args.logger = atm->logger; 566 rnsl_args.allocator = atm->allocator; 567 rnsl_args.verbose = atm->verbose; 568 res = rnsl_create(&rnsl_args, &rnsl); 569 if(res != RES_OK) { 570 log_err(atm, 571 "Failed to create loader for phase function list `%s' -- %s\n", 572 args->phase_fn_lst_filename, res_to_cstr(res)); 573 goto error; 574 } 575 576 /* Load the list of phase function paths */ 577 res = rnsl_load(rnsl, args->phase_fn_lst_filename); 578 if(res != RES_OK) goto error; 579 580 /* Reserve memory space for the list of phase functions */ 581 nphase_fn = rnsl_get_strings_count(rnsl); 582 res = darray_phase_fn_resize(&aerosol->phase_fn_lst, nphase_fn); 583 if(res != RES_OK) { 584 log_err(atm, "%s: could not allocate the list of %lu phase functions -- %s\n", 585 args->phase_fn_lst_filename, nphase_fn, res_to_cstr(res)); 586 goto error; 587 } 588 589 FOR_EACH(iphase_fn, 0, nphase_fn) { 590 const char* filename = rnsl_get_string(rnsl, iphase_fn); 591 struct phase_fn* phase_fn = 592 darray_phase_fn_data_get(&aerosol->phase_fn_lst)+iphase_fn; 593 594 res = load_phase_fn(atm, filename, &phase_fn->rnsf); 595 if(res != RES_OK) goto error; 596 res = phase_fn_create_phase_list(atm, filename, phase_fn); 597 if(res != RES_OK) goto error; 598 } 599 600 exit: 601 if(rnsl) RNSL(ref_put(rnsl)); 602 return res; 603 error: 604 darray_phase_fn_clear(&aerosol->phase_fn_lst); 605 goto exit; 606 } 607 608 static res_T 609 setup_gas_properties(struct rnatm* atm, const struct rnatm_gas_args* gas_args) 610 { 611 char buf[128]; 612 struct time t0, t1; 613 struct sck_load_args sck_load_args = SCK_LOAD_ARGS_NULL; 614 struct sck_band band_low = SCK_BAND_NULL; 615 struct sck_band band_upp = SCK_BAND_NULL; 616 struct sck_create_args sck_args = SCK_CREATE_ARGS_DEFAULT; 617 struct sbuf_create_args sbuf_args = SBUF_CREATE_ARGS_DEFAULT; 618 struct sbuf_desc sbuf_desc = SBUF_DESC_NULL; 619 size_t nbands; 620 621 res_T res = RES_OK; 622 ASSERT(atm && gas_args); 623 624 /* Start time recording */ 625 log_info(atm, "load gas properties\n"); 626 time_current(&t0); 627 628 /* Create the Rayleigh phase function */ 629 res = ssf_phase_create(atm->allocator, &ssf_phase_rayleigh, &atm->gas.rayleigh); 630 if(res != RES_OK) goto error; 631 632 /* Create the Star-Buffer loader */ 633 sbuf_args.logger = atm->logger; 634 sbuf_args.allocator = atm->allocator; 635 sbuf_args.verbose = atm->verbose; 636 res = sbuf_create(&sbuf_args, &atm->gas.temperatures); 637 if(res != RES_OK) goto error; 638 639 /* Load gas temperatures */ 640 res = sbuf_load(atm->gas.temperatures, gas_args->temperatures_filename); 641 if(res != RES_OK) goto error; 642 res = sbuf_get_desc(atm->gas.temperatures, &sbuf_desc); 643 if(res != RES_OK) goto error; 644 res = check_gas_temperatures_desc(atm, &sbuf_desc, gas_args); 645 if(res != RES_OK) goto error; 646 647 /* Create the Star-CK loader */ 648 sck_args.logger = atm->logger; 649 sck_args.allocator = atm->allocator; 650 sck_args.verbose = atm->verbose; 651 res = sck_create(&sck_args, &atm->gas.ck); 652 if(res != RES_OK) goto error; 653 654 /* Load correlated-K */ 655 sck_load_args.path = gas_args->sck_filename; 656 sck_load_args.memory_mapping = 1; 657 res = sck_load(atm->gas.ck, &sck_load_args); 658 if(res != RES_OK) goto error; 659 res = check_gas_ck_desc(atm, gas_args); 660 if(res != RES_OK) goto error; 661 662 /* Print elapsed time */ 663 time_sub(&t0, time_current(&t1), &t0); 664 time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf)); 665 log_info(atm, "gas properties loaded in %s\n", buf); 666 667 /* Print gas informations */ 668 nbands = sck_get_bands_count(atm->gas.ck); 669 res = sck_get_band(atm->gas.ck, 0, &band_low); 670 if(res != RES_OK) goto error; 671 res = sck_get_band(atm->gas.ck, nbands-1, &band_upp); 672 if(res != RES_OK) goto error; 673 log_info(atm, "the gas is composed of %lu band%sin [%g, %g[ nm\n", 674 nbands, nbands > 1 ? "s " : " ", 675 band_low.lower, 676 band_upp.upper); 677 678 exit: 679 return res; 680 error: 681 reset_gas_properties(&atm->gas); 682 goto exit; 683 } 684 685 static res_T 686 setup_aerosol_properties 687 (struct rnatm* atm, 688 struct aerosol* aerosol, 689 const struct rnatm_aerosol_args* aerosol_args) 690 { 691 char buf[128]; 692 size_t ibands[2]; 693 struct time t0, t1; 694 struct sars_create_args sars_args = SARS_CREATE_ARGS_DEFAULT; 695 struct sars_load_args sars_load_args = SARS_LOAD_ARGS_NULL; 696 struct sbuf_create_args sbuf_args = SBUF_CREATE_ARGS_DEFAULT; 697 struct sbuf_desc sbuf_desc = SBUF_DESC_NULL; 698 699 res_T res = RES_OK; 700 ASSERT(atm && aerosol_args); 701 702 /* Start time recording */ 703 log_info(atm, "load %s properties\n", str_cget(&aerosol->name)); 704 time_current(&t0); 705 706 /* Create the Star-Aerosol loader */ 707 sars_args.logger = atm->logger; 708 sars_args.allocator = atm->allocator; 709 sars_args.verbose = atm->verbose; 710 res = sars_create(&sars_args, &aerosol->sars); 711 if(res != RES_OK) goto error; 712 713 /* Load the aerosol radiative properties */ 714 sars_load_args.path = aerosol_args->sars_filename; 715 sars_load_args.memory_mapping = 1; 716 res = sars_load(aerosol->sars, &sars_load_args); 717 if(res != RES_OK) goto error; 718 res = check_aerosol_sars_desc(atm, aerosol, aerosol_args); 719 if(res != RES_OK) goto error; 720 721 /* Ignore the aerosol if it has no data for the input spectral range */ 722 res = sars_find_bands(aerosol->sars, atm->spectral_range, ibands); 723 if(res != RES_OK) goto error; 724 if(ibands[0] > ibands[1]) { 725 reset_aerosol_properties(aerosol); 726 goto exit; 727 } 728 729 /* Load the aerosol phase functions */ 730 res = load_phase_fn_list(atm, aerosol, aerosol_args); 731 if(res != RES_OK) goto error; 732 733 /* Create the Star-Buffer loader */ 734 sbuf_args.logger = atm->logger; 735 sbuf_args.allocator = atm->allocator; 736 sbuf_args.verbose = atm->verbose; 737 res = sbuf_create(&sbuf_args, &aerosol->phase_fn_ids); 738 if(res != RES_OK) goto error; 739 740 /* Load phase function ids */ 741 res = sbuf_load(aerosol->phase_fn_ids, aerosol_args->phase_fn_ids_filename); 742 if(res != RES_OK) goto error; 743 res = sbuf_get_desc(aerosol->phase_fn_ids, &sbuf_desc); 744 if(res != RES_OK) goto error; 745 res = check_aerosol_phase_fn_ids_desc(atm, aerosol, &sbuf_desc, aerosol_args); 746 if(res != RES_OK) goto error; 747 748 /* Print elapsed time */ 749 time_sub(&t0, time_current(&t1), &t0); 750 time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf)); 751 log_info(atm, "%s properties loaded in %s\n", 752 str_cget(&aerosol->name), buf); 753 754 exit: 755 return res; 756 error: 757 reset_aerosol_properties(aerosol); 758 goto exit; 759 } 760 761 static res_T 762 remove_useless_aerosols(struct rnatm* atm) 763 { 764 struct aerosol* aerosols = NULL; 765 size_t i, n; 766 res_T res = RES_OK; 767 ASSERT(atm); 768 769 aerosols = darray_aerosol_data_get(&atm->aerosols); 770 771 n = 0; 772 FOR_EACH(i, 0, darray_aerosol_size_get(&atm->aerosols)) { 773 /* Discard aerosols with no radiative properties */ 774 if(!aerosols[i].sars) { 775 log_warn(atm, 776 "discard %s aerosol: no radiative properties are defined for the " 777 "spectral domain [%g, %g] nm\n", 778 str_cget(&aerosols[i].name), 779 atm->spectral_range[0], 780 atm->spectral_range[1]); 781 continue; 782 } 783 /* Compact the aerosols to consider */ 784 if(i != n) { 785 aerosol_copy_and_release(aerosols+n, aerosols+i); 786 } 787 ++n; 788 } 789 790 if(!n) { 791 /* Discard all aerosols */ 792 darray_aerosol_purge(&atm->aerosols); 793 } else { 794 /* Resize the aerosols list */ 795 res = darray_aerosol_resize(&atm->aerosols, n); 796 if(res != RES_OK) goto error; 797 } 798 799 exit: 800 return res; 801 error: 802 goto exit; 803 } 804 805 static res_T 806 find_band_range 807 (const struct rnatm* atm, 808 const double range[2], /* In nanometers */ 809 size_t bands[2]) 810 { 811 struct sck_band band_low; 812 struct sck_band band_upp; 813 size_t nbands_overlaped; 814 size_t iband; 815 res_T res = RES_OK; 816 ASSERT(atm && range && bands && range[0] <= range[1]); 817 818 res = sck_find_bands(atm->gas.ck, range, bands); 819 if(res != RES_OK) goto error; 820 821 if(bands[0] > bands[1]) { 822 log_err(atm, 823 "the spectral range [%g, %g] nm does not overlap any bands\n", 824 SPLIT2(range)); 825 res = RES_BAD_ARG; 826 goto error; 827 } 828 829 /* Check that there is no hole in the spectral data */ 830 FOR_EACH(iband, bands[0], bands[1]) { 831 struct sck_band band_curr; 832 struct sck_band band_next; 833 834 SCK(get_band(atm->gas.ck, iband+0, &band_curr)); 835 SCK(get_band(atm->gas.ck, iband+1, &band_next)); 836 837 if(band_curr.upper != band_next.lower) { 838 log_err(atm, 839 "gas spectral data is missing in [%g, %g] nm. " 840 "There is a hole in [%g, %g] nm\n", 841 atm->spectral_range[0], 842 atm->spectral_range[1], 843 band_curr.upper, 844 band_next.lower); 845 res = RES_BAD_ARG; 846 goto error; 847 } 848 } 849 850 851 SCK(get_band(atm->gas.ck, bands[0], &band_low)); 852 SCK(get_band(atm->gas.ck, bands[1], &band_upp)); 853 if(band_low.lower > range[0] 854 || band_upp.upper < range[1]) { 855 log_err(atm, 856 "gas spectral data is missing. They are defined between [%g, %g] nm " 857 "while the required spectral range is [%g, %g]\n", 858 band_low.lower, 859 band_upp.upper, 860 range[0], 861 range[1]); 862 res = RES_BAD_ARG; 863 goto error; 864 } 865 866 nbands_overlaped = bands[1] - bands[0] + 1; 867 868 log_info(atm, 869 "the spectral range [%g, %g] nm overlaps %lu band%sin [%g, %g[ nm\n", 870 SPLIT2(range), 871 (unsigned long)nbands_overlaped, 872 nbands_overlaped > 1 ? "s ": " ", 873 band_low.lower, 874 band_upp.upper); 875 876 exit: 877 return res; 878 error: 879 goto exit; 880 } 881 882 static res_T 883 setup_spectral_range(struct rnatm* atm, const struct rnatm_create_args* args) 884 { 885 size_t iband, nbands; 886 size_t offset; 887 size_t i; 888 res_T res = RES_OK; 889 ASSERT(atm && args); 890 ASSERT(args->spectral_range[0] <= args->spectral_range[1]); 891 892 atm->spectral_range[0] = args->spectral_range[0]; 893 atm->spectral_range[1] = args->spectral_range[1]; 894 895 /* Find the bands overlapped by the input spectral range */ 896 res = find_band_range(atm, atm->spectral_range, atm->ibands); 897 if(res != RES_OK) goto error; 898 899 /* Allocate the list of bands to be taken into account */ 900 nbands = atm->ibands[1] - atm->ibands[0] + 1; 901 res = darray_band_resize(&atm->bands, nbands); 902 if(res != RES_OK) { 903 log_err(atm, 904 "failed to register the bands to be taken into account -- %s\n", 905 res_to_cstr(res)); 906 goto error; 907 } 908 909 /* Register the bands */ 910 i = 0; 911 offset = 0; 912 FOR_EACH(iband, atm->ibands[0], atm->ibands[1]+1) { 913 struct band* band = NULL; 914 struct sck_band sck_band; 915 916 SCK(get_band(atm->gas.ck, iband, &sck_band)); 917 918 band = darray_band_data_get(&atm->bands) + i; 919 band->index = iband; 920 band->nquad_pts = sck_band.quad_pts_count; 921 band->offset_accel_struct = offset; 922 923 offset += band->nquad_pts; 924 ++i; 925 } 926 927 exit: 928 return res; 929 error: 930 goto exit; 931 } 932 933 static res_T 934 cell_create_phase_fn_aerosol 935 (struct rnatm* atm, 936 const struct rnatm_cell_create_phase_fn_args* args, 937 struct ssf_phase** out_ssf_phase) 938 { 939 struct sbuf_desc sbuf_desc; 940 const double* bcoords = NULL; 941 struct aerosol* aerosol = NULL; 942 struct phase_fn* phase_fn = NULL; 943 struct ssf_phase* ssf_phase = NULL; 944 uint32_t phase_fn_id = 0; 945 size_t inode = 0; 946 size_t iphase_fn = 0; 947 int icell_node = 0; 948 res_T res = RES_OK; 949 ASSERT(atm && args && out_ssf_phase); 950 ASSERT(args->cell.component < darray_aerosol_size_get(&atm->aerosols)); 951 952 /* Get the aerosol the cell belongs to */ 953 aerosol = darray_aerosol_data_get(&atm->aerosols) + args->cell.component; 954 955 /* Sample the cell node to consider */ 956 bcoords = args->cell.barycentric_coords; 957 if(args->r[0] < bcoords[0]) { 958 icell_node = 0; 959 } else if(args->r[0] < bcoords[0] + bcoords[1]) { 960 icell_node = 1; 961 } else if(args->r[0] < bcoords[0] + bcoords[1] + bcoords[2]) { 962 icell_node = 2; 963 } else { 964 icell_node = 3; 965 } 966 967 /* Retrieve the phase function id of the node */ 968 SBUF(get_desc(aerosol->phase_fn_ids, &sbuf_desc)); 969 inode = args->cell.prim.indices[icell_node]; 970 ASSERT(inode < sbuf_desc.size); 971 phase_fn_id = *((const uint32_t*)sbuf_desc_at(&sbuf_desc, inode)); 972 973 /* Retrieve the spectrally varying phase function of the node */ 974 ASSERT(phase_fn_id < darray_phase_fn_size_get(&aerosol->phase_fn_lst)); 975 phase_fn = darray_phase_fn_data_get(&aerosol->phase_fn_lst)+phase_fn_id; 976 977 /* Get the phase function based on the input wavelength */ 978 res = rnsf_fetch_phase_fn 979 (phase_fn->rnsf, args->wavelength, args->r[1], &iphase_fn); 980 if(res != RES_OK) goto error; 981 ssf_phase = darray_phase_data_get(&phase_fn->phase_lst)[iphase_fn]; 982 ASSERT(ssf_phase); 983 SSF(phase_ref_get(ssf_phase)); 984 985 exit: 986 if(out_ssf_phase) *out_ssf_phase = ssf_phase; 987 return res; 988 error: 989 log_err(atm, "error creating the %s phase function -- %s\n", 990 str_cget(&aerosol->name), res_to_cstr(res)); 991 if(ssf_phase) { SSF(phase_ref_put(ssf_phase)); ssf_phase = NULL; } 992 goto exit; 993 } 994 995 static res_T 996 compute_unnormalized_cumulative_radcoef 997 (const struct rnatm* atm, 998 const enum rnatm_radcoef radcoef, 999 const struct rnatm_cell_pos* cells, 1000 const size_t iband, 1001 const size_t iquad, 1002 float cumulative[RNATM_MAX_COMPONENTS_COUNT], 1003 /* For debug */ 1004 const double k_min, 1005 const double k_max) 1006 { 1007 struct rnatm_cell_get_radcoef_args cell_args = RNATM_CELL_GET_RADCOEF_ARGS_NULL; 1008 size_t cpnt = RNATM_GAS; 1009 size_t icumul = 0; 1010 size_t naerosols = 0; 1011 float k = 0; 1012 res_T res = RES_OK; 1013 ASSERT(atm && cells && (unsigned)radcoef < RNATM_RADCOEFS_COUNT__); 1014 ASSERT(cumulative); 1015 (void)k_min; 1016 1017 naerosols = darray_aerosol_size_get(&atm->aerosols); 1018 ASSERT(naerosols+1 < RNATM_MAX_COMPONENTS_COUNT); 1019 1020 /* Setup the arguments common to all components */ 1021 cell_args.iband = iband; 1022 cell_args.iquad = iquad; 1023 cell_args.radcoef = radcoef; 1024 cell_args.k_max = k_max; /* For Debug */ 1025 1026 do { 1027 1028 cell_args.cell = cells[cpnt+1]; 1029 1030 /* Test if the component exists here */ 1031 if(!SUVM_PRIMITIVE_NONE(&cell_args.cell.prim)) { 1032 double per_cell_k; 1033 1034 /* Add the component's contribution to the radiative coefficient */ 1035 res = rnatm_cell_get_radcoef(atm, &cell_args, &per_cell_k); 1036 if(res != RES_OK) goto error; 1037 k += (float)per_cell_k; 1038 } 1039 1040 /* Update the cumulative */ 1041 cumulative[icumul] = k; 1042 ++icumul; 1043 } while(++cpnt < naerosols); 1044 1045 ASSERT(cumulative[icumul-1] == 0 || (float)k_min <= cumulative[icumul-1]); 1046 ASSERT(cumulative[icumul-1] == 0 || (float)k_max >= cumulative[icumul-1]); 1047 1048 exit: 1049 return res; 1050 error: 1051 goto exit; 1052 } 1053 1054 /******************************************************************************* 1055 * Exported functions 1056 ******************************************************************************/ 1057 res_T 1058 rnatm_get_radcoef 1059 (const struct rnatm* atm, 1060 const struct rnatm_get_radcoef_args* args, 1061 double* out_k) 1062 { 1063 float cumul[RNATM_MAX_COMPONENTS_COUNT]; 1064 size_t ncpnts; 1065 double k = 0; 1066 res_T res = RES_OK; 1067 1068 if(!atm || !out_k) { res = RES_BAD_ARG; goto error; } 1069 res = check_rnatm_get_radcoef_args(atm, args); 1070 if(res != RES_OK) goto error; 1071 1072 ncpnts = 1/*gas*/ + darray_aerosol_size_get(&atm->aerosols); 1073 ASSERT(ncpnts <= RNATM_MAX_COMPONENTS_COUNT); 1074 1075 /* Calculate the cumulative (unnormalized) of radiative coefficients. Its last 1076 * entry is the sum of the radiative coefficients of each component, which is 1077 * the atmospheric radiative coefficient to be returned. */ 1078 res = compute_unnormalized_cumulative_radcoef(atm, args->radcoef, 1079 args->cells, args->iband, args->iquad, cumul, args->k_min, args->k_max); 1080 if(res != RES_OK) goto error; 1081 1082 if(cumul[ncpnts-1] == 0) { 1083 k = 0; /* No atmospheric data */ 1084 } else { 1085 k = cumul[ncpnts-1]; 1086 ASSERT(args->k_min <= k && k <= args->k_max); 1087 } 1088 1089 exit: 1090 *out_k = k; 1091 return res; 1092 error: 1093 goto exit; 1094 } 1095 1096 res_T 1097 rnatm_sample_component 1098 (const struct rnatm* atm, 1099 const struct rnatm_sample_component_args* args, 1100 size_t* cpnt) 1101 { 1102 float cumul[RNATM_MAX_COMPONENTS_COUNT]; 1103 float norm; 1104 size_t ncpnts; 1105 size_t i; 1106 res_T res = RES_OK; 1107 1108 if(!atm || !cpnt) { res = RES_BAD_ARG; goto error; } 1109 res = check_rnatm_sample_component_args(atm, args); 1110 if(res != RES_OK) goto error; 1111 1112 ncpnts = 1/*gas*/ + darray_aerosol_size_get(&atm->aerosols); 1113 ASSERT(ncpnts <= RNATM_MAX_COMPONENTS_COUNT); 1114 1115 /* Discard the calculation of the cumulative if there is only gas */ 1116 if(ncpnts == 1) { 1117 ASSERT(!SUVM_PRIMITIVE_NONE(&args->cells[0].prim)); 1118 *cpnt = RNATM_GAS; 1119 goto exit; 1120 } 1121 1122 res = compute_unnormalized_cumulative_radcoef(atm, args->radcoef, args->cells, 1123 args->iband, args->iquad, cumul, -DBL_MAX, DBL_MAX); 1124 if(res != RES_OK) goto error; 1125 ASSERT(cumul[ncpnts-1] > 0); 1126 1127 /* Normalize the cumulative */ 1128 norm = cumul[ncpnts-1]; 1129 FOR_EACH(i, 0, ncpnts) cumul[i] /= norm; 1130 cumul[ncpnts-1] = 1.f; /* Handle precision issues */ 1131 1132 /* Use a simple linear search to sample the component since there are 1133 * usually very few aerosols to consider */ 1134 FOR_EACH(i, 0, ncpnts) if(args->r < cumul[i]) break; 1135 1136 *cpnt = i-1; 1137 ASSERT(*cpnt == RNATM_GAS || *cpnt < darray_aerosol_size_get(&atm->aerosols)); 1138 ASSERT(!SUVM_PRIMITIVE_NONE(&args->cells[i].prim)); 1139 1140 exit: 1141 return res; 1142 error: 1143 goto exit; 1144 } 1145 1146 res_T 1147 rnatm_cell_get_radcoef 1148 (const struct rnatm* atm, 1149 const struct rnatm_cell_get_radcoef_args* args, 1150 double* out_k) 1151 { 1152 float vtx_k[4]; 1153 double k = NaN; 1154 res_T res = RES_OK; 1155 1156 if(!atm || !out_k) { res = RES_BAD_ARG; goto error; } 1157 res = check_rnatm_cell_get_radcoef_args(atm, args); 1158 if(res != RES_OK) goto error; 1159 1160 /* Get the radiative coefficients of the tetrahedron */ 1161 tetra_get_radcoef(atm, &args->cell.prim, args->cell.component, args->iband, 1162 args->iquad, args->radcoef, vtx_k); 1163 1164 if(vtx_k[0] == vtx_k[1] 1165 && vtx_k[0] == vtx_k[2] 1166 && vtx_k[0] == vtx_k[3]) { 1167 /* Avoid precision issue when iterpolating the same value */ 1168 k = vtx_k[0]; 1169 } else { 1170 float min_vtx_k; 1171 float max_vtx_k; 1172 1173 /* Calculate the radiative coefficient by linearly interpolating the 1174 * coefficients defined by tetrahedron vertex */ 1175 k = vtx_k[0] * args->cell.barycentric_coords[0] 1176 + vtx_k[1] * args->cell.barycentric_coords[1] 1177 + vtx_k[2] * args->cell.barycentric_coords[2] 1178 + vtx_k[3] * args->cell.barycentric_coords[3]; 1179 1180 /* Fix interpolation accuracy issues */ 1181 min_vtx_k = MMIN(MMIN(vtx_k[0], vtx_k[1]), MMIN(vtx_k[2], vtx_k[3])); 1182 max_vtx_k = MMAX(MMAX(vtx_k[0], vtx_k[1]), MMAX(vtx_k[2], vtx_k[3])); 1183 k = CLAMP((float)k, min_vtx_k, max_vtx_k); 1184 } 1185 ASSERT(k <= args->k_max); 1186 1187 exit: 1188 *out_k = k; 1189 return res; 1190 error: 1191 goto exit; 1192 } 1193 1194 res_T 1195 rnatm_cell_create_phase_fn 1196 (struct rnatm* atm, 1197 const struct rnatm_cell_create_phase_fn_args* args, 1198 struct ssf_phase** out_ssf_phase) 1199 { 1200 struct ssf_phase* ssf_phase = NULL; 1201 res_T res = RES_OK; 1202 1203 if(!atm || !out_ssf_phase) { res = RES_BAD_ARG; goto error; } 1204 res = check_rnatm_cell_create_phase_fn_args(atm, args); 1205 if(res != RES_OK) goto error; 1206 1207 if(args->cell.component == RNATM_GAS) { 1208 /* Get a reference on the pre-allocated Rayleigh phase function */ 1209 ssf_phase = atm->gas.rayleigh; 1210 SSF(phase_ref_get(ssf_phase)); 1211 } else { 1212 res = cell_create_phase_fn_aerosol(atm, args, &ssf_phase); 1213 if(res != RES_OK) goto error; 1214 } 1215 1216 exit: 1217 if(out_ssf_phase) *out_ssf_phase = ssf_phase; 1218 return res; 1219 error: 1220 if(ssf_phase) { SSF(phase_ref_put(ssf_phase)); ssf_phase = NULL; } 1221 goto exit; 1222 } 1223 1224 res_T 1225 rnatm_cell_get_gas_temperature 1226 (const struct rnatm* atm, 1227 const struct rnatm_cell_pos* cell, 1228 double* temperature) 1229 { 1230 struct sbuf_desc sbuf_desc; 1231 float vtx_T[4]; 1232 res_T res = RES_OK; 1233 1234 if(!atm || !temperature) return RES_OK; 1235 res = check_cell(cell); 1236 if(res != RES_OK) goto error; 1237 1238 /* Invalid cell regarding gas */ 1239 SBUF(get_desc(atm->gas.temperatures, &sbuf_desc)); 1240 if(cell->prim.indices[0] >= sbuf_desc.size 1241 || cell->prim.indices[1] >= sbuf_desc.size 1242 || cell->prim.indices[2] >= sbuf_desc.size 1243 || cell->prim.indices[3] >= sbuf_desc.size) 1244 return RES_BAD_ARG; 1245 1246 /* Get the temperature per vertex */ 1247 vtx_T[0] = *((float*)sbuf_desc_at(&sbuf_desc, cell->prim.indices[0])); 1248 vtx_T[1] = *((float*)sbuf_desc_at(&sbuf_desc, cell->prim.indices[1])); 1249 vtx_T[2] = *((float*)sbuf_desc_at(&sbuf_desc, cell->prim.indices[2])); 1250 vtx_T[3] = *((float*)sbuf_desc_at(&sbuf_desc, cell->prim.indices[3])); 1251 1252 1253 /* Calculate the temperature at the input position by linearly interpolating 1254 * the temperatures defined by tetrahedron vertex */ 1255 *temperature = 1256 vtx_T[0] * cell->barycentric_coords[0] 1257 + vtx_T[1] * cell->barycentric_coords[1] 1258 + vtx_T[2] * cell->barycentric_coords[2] 1259 + vtx_T[3] * cell->barycentric_coords[3]; 1260 1261 exit: 1262 return res; 1263 error: 1264 goto exit; 1265 } 1266 1267 /******************************************************************************* 1268 * Local functions 1269 ******************************************************************************/ 1270 res_T 1271 setup_properties(struct rnatm* atm, const struct rnatm_create_args* args) 1272 { 1273 size_t i = 0; 1274 res_T res = RES_OK; 1275 ASSERT(atm && args); 1276 1277 res = setup_gas_properties(atm, &args->gas); 1278 if(res != RES_OK) goto error; 1279 1280 res = setup_spectral_range(atm, args); 1281 if(res != RES_OK) goto error; 1282 1283 FOR_EACH(i, 0, darray_aerosol_size_get(&atm->aerosols)) { 1284 struct aerosol* aerosol = darray_aerosol_data_get(&atm->aerosols)+i; 1285 res = setup_aerosol_properties(atm, aerosol, args->aerosols+i); 1286 if(res != RES_OK) goto error; 1287 } 1288 1289 /* Remove aerosols with no radiative properties for the input spectral range */ 1290 res = remove_useless_aerosols(atm); 1291 if(res != RES_OK) goto error; 1292 1293 exit: 1294 return res; 1295 error: 1296 reset_gas_properties(&atm->gas); 1297 FOR_EACH(i, 0, darray_aerosol_size_get(&atm->aerosols)) { 1298 struct aerosol* aerosol = darray_aerosol_data_get(&atm->aerosols)+i; 1299 reset_aerosol_properties(aerosol); 1300 } 1301 goto exit; 1302 } 1303 1304 res_T 1305 check_properties(const struct rnatm* atm) 1306 { 1307 size_t i; 1308 res_T res = RES_OK; 1309 ASSERT(atm); 1310 1311 res = check_gas_temperatures(atm); 1312 if(res != RES_OK) goto error; 1313 1314 FOR_EACH(i, 0, darray_aerosol_size_get(&atm->aerosols)) { 1315 const struct aerosol* aerosol = darray_aerosol_cdata_get(&atm->aerosols)+i; 1316 res = check_aerosol_phase_fn_ids(atm, aerosol); 1317 if(res != RES_OK) goto error; 1318 } 1319 1320 exit: 1321 return res; 1322 error: 1323 goto exit; 1324 }