htsky.c (37316B)
1 /* Copyright (C) 2018, 2019, 2020, 2021 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2018, 2019 Centre National de la Recherche Scientifique 3 * Copyright (C) 2018, 2019 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 200809L /* stat.st_time support */ 19 20 #include "htsky.h" 21 #include "htsky_c.h" 22 #include "htsky_atmosphere.h" 23 #include "htsky_cloud.h" 24 #include "htsky_log.h" 25 26 #include <high_tune/htcp.h> 27 #include <high_tune/htgop.h> 28 #include <high_tune/htmie.h> 29 30 #include <star/svx.h> 31 32 #include <rsys/clock_time.h> 33 #include <rsys/double3.h> 34 35 #include <errno.h> 36 #include <fcntl.h> /* open */ 37 #include <unistd.h> 38 #include <sys/stat.h> /* S_IRUSR & S_IWUSR */ 39 40 #include <omp.h> 41 42 /* Current version the cache structure. One should increment it and perform a 43 * version management onto serialized data when the cache data data structure 44 * is updated. */ 45 static const int CACHE_VERSION = 0; 46 47 /******************************************************************************* 48 * Helper function 49 ******************************************************************************/ 50 static INLINE int 51 check_args(const struct htsky_args* args) 52 { 53 return args 54 && args->htgop_filename 55 && args->grid_max_definition[0] 56 && args->grid_max_definition[1] 57 && args->grid_max_definition[2] 58 && args->name 59 && args->nthreads 60 && args->optical_thickness >= 0 61 && (unsigned)args->spectral_type < HTSKY_SPECTRAL_TYPES_COUNT__ 62 && args->wlen_range[0] <= args->wlen_range[1]; 63 } 64 65 static INLINE const char* 66 spectral_type_string(const enum htsky_spectral_type type) 67 { 68 const char* str = NULL; 69 switch(type) { 70 case HTSKY_SPECTRAL_LW: str = "longwave"; break; 71 case HTSKY_SPECTRAL_SW: str = "shortwave"; break; 72 default: FATAL("Unreachable code.\n"); break; 73 } 74 return str; 75 } 76 77 static res_T 78 setup_bands_properties(struct htsky* sky) 79 { 80 size_t nbands; 81 size_t iband; 82 res_T res = RES_OK; 83 ASSERT(sky); 84 85 nbands = htsky_get_spectral_bands_count(sky); 86 ASSERT(nbands); 87 88 sky->bands = MEM_CALLOC(sky->allocator, nbands, sizeof(*sky->bands)); 89 if(!sky->bands) { 90 log_err(sky, "Could not allocate the list of %s band properties.\n", 91 spectral_type_string(sky->spectral_type)); 92 res = RES_MEM_ERR; 93 goto error; 94 } 95 FOR_EACH(iband, sky->bands_range[0], sky->bands_range[1]+1) { 96 struct htgop_spectral_interval band; 97 double band_wlens[2]; 98 const size_t i = iband - sky->bands_range[0]; 99 100 switch(sky->spectral_type) { 101 case HTSKY_SPECTRAL_LW: 102 HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); 103 break; 104 case HTSKY_SPECTRAL_SW: 105 HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); 106 break; 107 default: FATAL("Unreachable code.\n"); break; 108 } 109 band_wlens[0] = wavenumber_to_wavelength(band.wave_numbers[1]); 110 band_wlens[1] = wavenumber_to_wavelength(band.wave_numbers[0]); 111 ASSERT(band_wlens[0] < band_wlens[1]); 112 113 sky->bands[i].Ca_avg = htmie_compute_xsection_absorption_average 114 (sky->htmie, band_wlens, HTMIE_FILTER_LINEAR); 115 sky->bands[i].Cs_avg = htmie_compute_xsection_scattering_average 116 (sky->htmie, band_wlens, HTMIE_FILTER_LINEAR); 117 sky->bands[i].g_avg = htmie_compute_asymmetry_parameter_average 118 (sky->htmie, band_wlens, HTMIE_FILTER_LINEAR); 119 ASSERT(sky->bands[i].Ca_avg > 0); 120 ASSERT(sky->bands[i].Cs_avg > 0); 121 ASSERT(sky->bands[i].g_avg > 0); 122 } 123 124 exit: 125 return res; 126 error: 127 if(sky->bands) { 128 MEM_RM(sky->allocator, sky->bands); 129 sky->bands = NULL; 130 } 131 goto exit; 132 } 133 134 static INLINE double 135 particle_fetch_raw_property 136 (const struct htsky* sky, 137 const enum htsky_property prop, 138 const size_t iband, 139 const size_t iquad, 140 const size_t ivox[3]) 141 { 142 double rho_da = 0; /* Dry air density */ 143 double rct = 0; /* #droplets in kg of water per kg of dry air */ 144 double ql = 0; /* Droplet density In kg.m^-3 */ 145 double Ca = 0; /* Massic absorption cross section in m^2.kg^-1 */ 146 double Cs = 0; /* Massic scattering cross section in m^2.kg^-1 */ 147 double k_particle = 0; 148 size_t i = 0; 149 (void)iquad; 150 ASSERT(sky && ivox); 151 152 /* Compute he dry air density */ 153 rho_da = cloud_dry_air_density(&sky->htcp_desc, ivox); 154 155 /* Compute the droplet density */ 156 rct = htcp_desc_RCT_at(&sky->htcp_desc, ivox[0], ivox[1], ivox[2], 0); 157 ql = rho_da * rct; 158 159 i = iband - sky->bands_range[0]; 160 161 /* Use the average cross section of the current spectral band */ 162 if(prop == HTSKY_Ka || prop == HTSKY_Kext) Ca = sky->bands[i].Ca_avg; 163 if(prop == HTSKY_Ks || prop == HTSKY_Kext) Cs = sky->bands[i].Cs_avg; 164 165 k_particle = ql*(Ca + Cs); 166 return k_particle; 167 } 168 169 static INLINE double 170 gas_fetch_raw_property 171 (const struct htsky* sky, 172 const enum htsky_property prop, 173 const size_t iband, 174 const size_t iquad, 175 const int in_clouds, 176 const double pos[3], 177 const size_t ivox[3]) 178 { 179 struct htgop_layer layer; 180 size_t ilayer = 0; 181 double k_gas = 0; 182 ASSERT(sky && pos && ivox); 183 184 /* Retrieve the quadrature point into the spectral band of the layer into 185 * which `pos' lies */ 186 HTGOP(position_to_layer_id(sky->htgop, pos[2], &ilayer)); 187 HTGOP(get_layer(sky->htgop, ilayer, &layer)); 188 189 if(sky->spectral_type == HTSKY_SPECTRAL_LW) { 190 struct htgop_layer_lw_spectral_interval band; 191 HTGOP(layer_get_lw_spectral_interval(&layer, iband, &band)); 192 193 if(!in_clouds) { 194 /* Pos is outside the clouds. Directly fetch the nominal optical 195 * properties */ 196 ASSERT(iquad < band.quadrature_length); 197 switch(prop) { 198 case HTSKY_Ka: 199 case HTSKY_Kext: 200 k_gas = band.ka_nominal[iquad]; 201 break; 202 case HTSKY_Ks: k_gas = 0; break; 203 default: FATAL("Unreachable code.\n"); break; 204 } 205 } else { 206 /* Pos is inside the clouds. Compute the water vapor molar fraction at 207 * the current voxel */ 208 const double x_h2o = cloud_water_vapor_molar_fraction(&sky->htcp_desc, ivox); 209 struct htgop_layer_lw_spectral_interval_tab tab; 210 211 /* Retrieve the tabulated data for the quadrature point */ 212 HTGOP(layer_lw_spectral_interval_get_tab(&band, iquad, &tab)); 213 214 /* Fetch the optical properties wrt x_h2o */ 215 switch(prop) { 216 case HTSKY_Ka: 217 case HTSKY_Kext: 218 HTGOP(layer_lw_spectral_interval_tab_fetch_ka(&tab, x_h2o, &k_gas)); 219 break; 220 case HTSKY_Ks: k_gas = 0; break; 221 default: FATAL("Unreachable code.\n"); break; 222 } 223 } 224 } else { 225 struct htgop_layer_sw_spectral_interval band; 226 ASSERT(sky->spectral_type == HTSKY_SPECTRAL_SW); 227 228 HTGOP(layer_get_sw_spectral_interval(&layer, iband, &band)); 229 if(!in_clouds) { 230 /* Pos is outside the clouds. Directly fetch the nominal optical 231 * properties */ 232 ASSERT(iquad < band.quadrature_length); 233 switch(prop) { 234 case HTSKY_Ka: k_gas = band.ka_nominal[iquad]; break; 235 case HTSKY_Ks: k_gas = band.ks_nominal[iquad]; break; 236 case HTSKY_Kext: 237 k_gas = band.ka_nominal[iquad] + band.ks_nominal[iquad]; 238 break; 239 default: FATAL("Unreachable code.\n"); break; 240 } 241 } else { 242 /* Pos is inside the clouds. Compute the water vapor molar fraction at 243 * the current voxel */ 244 const double x_h2o = cloud_water_vapor_molar_fraction(&sky->htcp_desc, ivox); 245 struct htgop_layer_sw_spectral_interval_tab tab; 246 247 /* Retrieve the tabulated data for the quadrature point */ 248 HTGOP(layer_sw_spectral_interval_get_tab(&band, iquad, &tab)); 249 250 /* Fetch the optical properties wrt x_h2o */ 251 switch(prop) { 252 case HTSKY_Ka: 253 HTGOP(layer_sw_spectral_interval_tab_fetch_ka(&tab, x_h2o, &k_gas)); 254 break; 255 case HTSKY_Ks: 256 HTGOP(layer_sw_spectral_interval_tab_fetch_ks(&tab, x_h2o, &k_gas)); 257 break; 258 case HTSKY_Kext: 259 HTGOP(layer_sw_spectral_interval_tab_fetch_kext(&tab, x_h2o, &k_gas)); 260 break; 261 default: FATAL("Unreachable code.\n"); break; 262 } 263 } 264 } 265 return k_gas; 266 } 267 268 static res_T 269 setup_cache_stream 270 (struct htsky* sky, 271 const char* htcp_filename, 272 const char* htgop_filename, 273 const char* htmie_filename, 274 const char* cache_filename, 275 int* out_create_cache, /* Define if the cache file was created */ 276 FILE** out_fp) 277 { 278 FILE* fp = NULL; 279 struct stat htcp_statbuf; 280 struct stat htgop_statbuf; 281 struct stat htmie_statbuf; 282 int create_cache = 0; 283 int fd = -1; 284 res_T res = RES_OK; 285 ASSERT(sky && out_create_cache && out_fp); 286 ASSERT(htcp_filename && htgop_filename && htmie_filename && cache_filename); 287 288 /* Open the cache file */ 289 fd = open(cache_filename, O_CREAT|O_EXCL|O_RDWR, S_IRUSR|S_IWUSR); 290 if(fd >= 0) { 291 create_cache = 1; 292 } else if (errno == EEXIST) { /* The cache already exists */ 293 fd = open(cache_filename, O_RDWR, 0); 294 } 295 296 if(fd < 0) { 297 log_err(sky, "Unexpected error while opening the cache file `%s'.\n", 298 cache_filename); 299 res = RES_IO_ERR; 300 goto error; 301 } 302 303 fp = fdopen(fd, "w+"); 304 if(!fp) { 305 log_err(sky, "Could not open the cache file `%s'.\n", cache_filename); 306 res = RES_IO_ERR; 307 goto error; 308 } 309 310 /* Query the status of the input */ 311 #define STAT(Filename, Statbuf) { \ 312 const int err = stat(Filename, Statbuf); \ 313 if(err) { \ 314 log_err(sky, "%s: could not stat the file -- %s\n", \ 315 Filename, strerror(errno)); \ 316 res = RES_IO_ERR; \ 317 goto error; \ 318 } \ 319 } (void)0 320 STAT(htcp_filename, &htcp_statbuf); 321 STAT(htgop_filename, &htgop_statbuf); 322 STAT(htmie_filename, &htmie_statbuf); 323 #undef STAT 324 325 if(create_cache) { 326 /* Setup the cache header, i.e. data that uniquely identify the cache 327 * regarding the input files (htcp, htmie and htgop files) */ 328 #define WRITE(Var, N) { \ 329 if(fwrite((Var), sizeof(*(Var)), (N), fp) != (N)) { \ 330 log_err(sky, "%s: could not write the cache header.\n",cache_filename);\ 331 res = RES_IO_ERR; \ 332 goto error; \ 333 } \ 334 } (void)0 335 WRITE(&CACHE_VERSION, 1); 336 WRITE(&htcp_statbuf.st_ino, 1); 337 WRITE(&htcp_statbuf.st_mtim, 1); 338 WRITE(&htgop_statbuf.st_ino, 1); 339 WRITE(&htgop_statbuf.st_mtim, 1); 340 WRITE(&htmie_statbuf.st_ino, 1); 341 WRITE(&htmie_statbuf.st_mtim, 1); 342 WRITE(&sky->spectral_type, 1); 343 WRITE(sky->bands_range, 2); 344 #undef WRITE 345 CHK(fflush(fp) == 0); 346 } else { 347 struct stat htcp_statbuf2; 348 struct stat htgop_statbuf2; 349 struct stat htmie_statbuf2; 350 int cache_version; 351 enum htsky_spectral_type spectral_type; 352 size_t bands_range[2]; 353 354 /* Read the cache header */ 355 #define READ(Var, N) { \ 356 if(fread((Var), sizeof(*(Var)), (N), fp) != (N)) { \ 357 if(feof(fp)) { \ 358 res = RES_BAD_ARG; \ 359 } else if(ferror(fp)) { \ 360 res = RES_IO_ERR; \ 361 } else { \ 362 res = RES_UNKNOWN_ERR; \ 363 } \ 364 log_err(sky, "%s: could not read the cache header.\n",cache_filename); \ 365 goto error; \ 366 } \ 367 } (void)0 368 READ(&cache_version, 1); 369 if(cache_version != CACHE_VERSION) { 370 log_err(sky, 371 "%s: invalid cache in version %d. Expecting a cache in version %d.\n", 372 cache_filename, cache_version, CACHE_VERSION); 373 res = RES_BAD_ARG; 374 goto error; 375 } 376 READ(&htcp_statbuf2.st_ino, 1); 377 READ(&htcp_statbuf2.st_mtim, 1); 378 READ(&htgop_statbuf2.st_ino, 1); 379 READ(&htgop_statbuf2.st_mtim, 1); 380 READ(&htmie_statbuf2.st_ino, 1); 381 READ(&htmie_statbuf2.st_mtim, 1); 382 READ(&spectral_type, 1); 383 READ(bands_range, 2); 384 #undef READ 385 386 /* Compare the cache header with the input file status to check that the 387 * cached data matched the input data */ 388 #define CHK_STAT(Stat0, Stat1) { \ 389 if((Stat0)->st_ino != (Stat1)->st_ino \ 390 || (Stat0)->st_mtim.tv_sec != (Stat1)->st_mtim.tv_sec \ 391 || (Stat0)->st_mtim.tv_nsec != (Stat1)->st_mtim.tv_nsec) { \ 392 log_err(sky, "%s: invalid cache regarding the input files.\n", \ 393 cache_filename); \ 394 res = RES_BAD_ARG; \ 395 goto error; \ 396 } \ 397 } (void)0 398 CHK_STAT(&htcp_statbuf, &htcp_statbuf2); 399 CHK_STAT(&htgop_statbuf, &htgop_statbuf2); 400 CHK_STAT(&htmie_statbuf, &htmie_statbuf2); 401 #undef CHK_STAT 402 403 /* Compare the handled spectral bands with the bands to handled to check 404 * that the cached octress are the expected ones */ 405 if(spectral_type != sky->spectral_type 406 || bands_range[0] != sky->bands_range[0] 407 || bands_range[1] != sky->bands_range[1]) { 408 log_err(sky, "%s: invalid cache regarding the wavelengths to handle.\n", 409 cache_filename); 410 res = RES_BAD_ARG; 411 goto error; 412 } 413 } 414 415 exit: 416 *out_fp = fp; 417 *out_create_cache = create_cache; 418 return res; 419 error: 420 if(fp) { 421 CHK(fclose(fp) == 0); 422 fp = NULL; 423 } else if(fd >= 0) { 424 CHK(close(fd) == 0); 425 } 426 goto exit; 427 } 428 429 static void 430 print_spectral_info(const struct htsky* sky) 431 { 432 struct htgop_spectral_interval band_low, band_upp; 433 size_t iband_low, iband_upp; 434 size_t nbands; 435 size_t i; 436 size_t naccels = 0; 437 ASSERT(sky); 438 439 nbands = htsky_get_spectral_bands_count(sky); 440 441 iband_low = htsky_get_spectral_band_id(sky, 0); 442 iband_upp = htsky_get_spectral_band_id(sky, nbands-1); 443 444 /* Retrieve the spectral interval boundaries */ 445 switch(sky->spectral_type) { 446 case HTSKY_SPECTRAL_LW: 447 HTGOP(get_lw_spectral_interval(sky->htgop, iband_low, &band_low)); 448 HTGOP(get_lw_spectral_interval(sky->htgop, iband_upp, &band_upp)); 449 break; 450 case HTSKY_SPECTRAL_SW: 451 HTGOP(get_sw_spectral_interval(sky->htgop, iband_low, &band_low)); 452 HTGOP(get_sw_spectral_interval(sky->htgop, iband_upp, &band_upp)); 453 break; 454 default: FATAL("Unreachable code.\n"); break; 455 } 456 457 log_info(sky, "Sky data defined in [%g, %g] nanometers over %lu %s.\n", 458 wavenumber_to_wavelength(band_upp.wave_numbers[1]), 459 wavenumber_to_wavelength(band_low.wave_numbers[0]), 460 (unsigned long)nbands, 461 nbands > 1 ? "bands" : "band"); 462 463 /* Compute the overall number of sky acceleration structures to build */ 464 FOR_EACH(i, 0, nbands) { 465 struct htgop_spectral_interval band; 466 const size_t iband = htsky_get_spectral_band_id(sky, i); 467 468 switch(sky->spectral_type) { 469 case HTSKY_SPECTRAL_LW: 470 HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); 471 break; 472 case HTSKY_SPECTRAL_SW: 473 HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); 474 break; 475 default: FATAL("Unreachable code.\n"); break; 476 } 477 naccels += band.quadrature_length; 478 } 479 480 log_info(sky, "Number of clouds partitionning structures: %lu\n", 481 (unsigned long)naccels); 482 } 483 484 static void 485 release_sky(ref_T* ref) 486 { 487 struct htsky* sky; 488 ASSERT(ref); 489 sky = CONTAINER_OF(ref, struct htsky, ref); 490 cloud_clean(sky); 491 atmosphere_clean(sky); 492 if(sky->svx) SVX(device_ref_put(sky->svx)); 493 if(sky->htcp) HTCP(ref_put(sky->htcp)); 494 if(sky->htgop) HTGOP(ref_put(sky->htgop)); 495 if(sky->htmie) HTMIE(ref_put(sky->htmie)); 496 if(sky->bands) MEM_RM(sky->allocator, sky->bands); 497 if(sky->logger == &sky->logger__) logger_release(&sky->logger__); 498 darray_split_release(&sky->svx2htcp_z); 499 str_release(&sky->name); 500 ASSERT(MEM_ALLOCATED_SIZE(&sky->svx_allocator) == 0); 501 mem_shutdown_proxy_allocator(&sky->svx_allocator); 502 MEM_RM(sky->allocator, sky); 503 } 504 505 /******************************************************************************* 506 * Local functions 507 ******************************************************************************/ 508 res_T 509 htsky_create 510 (struct logger* logger, /* NULL <=> use default logger */ 511 struct mem_allocator* mem_allocator, /* NULL <=> use default allocator */ 512 const struct htsky_args* args, 513 struct htsky** out_sky) 514 { 515 struct time t0, t1; 516 struct mem_allocator* allocator = NULL; 517 struct htsky* sky = NULL; 518 double wnums[2]; 519 char buf[128]; 520 int nthreads_max; 521 int force_cache_upd = 0; 522 FILE* cache = NULL; 523 res_T res = RES_OK; 524 525 if(!check_args(args) || !out_sky) { 526 res = RES_BAD_ARG; 527 goto error; 528 } 529 530 allocator = mem_allocator ? mem_allocator : &mem_default_allocator; 531 sky = MEM_CALLOC(allocator, 1, sizeof(*sky)); 532 if(!sky) { 533 if(args->verbose) { 534 #define ERR_STR "Could not allocate the HTSky data structure.\n" 535 if(logger) { 536 logger_print(logger, LOG_ERROR, ERR_STR); 537 } else { 538 fprintf(stderr, MSG_ERROR_PREFIX ERR_STR); 539 } 540 #undef ERR_STR 541 } 542 res = RES_MEM_ERR; 543 goto error; 544 } 545 nthreads_max = MMAX(omp_get_max_threads(), omp_get_num_procs()); 546 ref_init(&sky->ref); 547 sky->allocator = allocator; 548 sky->verbose = args->verbose; 549 sky->spectral_type = args->spectral_type ; 550 sky->repeat_clouds = args->repeat_clouds; 551 sky->is_cloudy = args->htcp_filename != NULL; 552 darray_split_init(sky->allocator, &sky->svx2htcp_z); 553 str_init(sky->allocator, &sky->name); 554 sky->bands_range[0] = 1; 555 sky->bands_range[1] = 0; 556 sky->nthreads = MMIN(args->nthreads, (unsigned)nthreads_max); 557 558 if(logger) { 559 sky->logger = logger; 560 } else { 561 setup_log_default(sky); 562 } 563 564 res = str_set(&sky->name, args->name); 565 if(res != RES_OK) { 566 log_err(sky, "Cannot setup the sky name to `%s'.\n", args->name); 567 goto error; 568 } 569 570 /* Setup an allocator specific to the SVX library */ 571 res = mem_init_proxy_allocator(&sky->svx_allocator, sky->allocator); 572 if(res != RES_OK) { 573 log_err(sky, "Cannot init the allocator used to manage the Star-VX data.\n"); 574 goto error; 575 } 576 577 /* Create the Star-VX library device */ 578 res = svx_device_create 579 (sky->logger, &sky->svx_allocator, sky->verbose, &sky->svx); 580 if(res != RES_OK) { 581 log_err(sky, "Error creating the Star-VX library device.\n"); 582 goto error; 583 } 584 585 /* Load the gas optical properties */ 586 res = htgop_create(sky->logger, sky->allocator, sky->verbose, &sky->htgop); 587 if(res != RES_OK) { 588 log_err(sky, "Could not create the gas optical properties loader.\n"); 589 goto error; 590 } 591 res = htgop_load(sky->htgop, args->htgop_filename); 592 if(res != RES_OK) { 593 log_err(sky, "Error loading the gas optical properties -- `%s'.\n", 594 args->htgop_filename); 595 goto error; 596 } 597 598 /* Retrieve the spectral bands */ 599 wnums[0] = wavelength_to_wavenumber(args->wlen_range[1]); 600 wnums[1] = wavelength_to_wavenumber(args->wlen_range[0]); 601 switch(sky->spectral_type) { 602 case HTSKY_SPECTRAL_LW: 603 res = htgop_get_lw_spectral_intervals(sky->htgop, wnums, sky->bands_range); 604 break; 605 case HTSKY_SPECTRAL_SW: 606 res = htgop_get_sw_spectral_intervals(sky->htgop, wnums, sky->bands_range); 607 break; 608 default: FATAL("Unreachable code.\n"); break; 609 } 610 if(res != RES_OK) goto error; 611 612 print_spectral_info(sky); 613 614 /* Setup the atmopshere */ 615 time_current(&t0); 616 res = atmosphere_setup(sky, args->optical_thickness); 617 if(res != RES_OK) goto error; 618 time_sub(&t0, time_current(&t1), &t0); 619 time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf)); 620 log_info(sky, "Setup atmosphere in %s\n", buf); 621 622 /* Nothing more to do */ 623 if(!sky->is_cloudy) goto exit; 624 625 if(!args->htmie_filename) { 626 log_err(sky, "Missing the HTMie filename.\n"); 627 res = RES_BAD_ARG; 628 goto error; 629 } 630 631 if(!args->htcp_filename) { 632 log_err(sky, "Missing the HTCP filename.\n"); 633 res = RES_BAD_ARG; 634 goto error; 635 } 636 637 /* Load MIE data */ 638 res = htmie_create(sky->logger, sky->allocator, sky->verbose, &sky->htmie); 639 if(res != RES_OK) { 640 log_err(sky, "Could not create the Mie's data loader.\n"); 641 goto error; 642 } 643 res = htmie_load(sky->htmie, args->htmie_filename); 644 if(res != RES_OK) { 645 log_err(sky, "Error loading the Mie's data -- `%s'.\n", args->htmie_filename); 646 goto error; 647 } 648 649 /* Setup the properties of the Short/Long wave bands */ 650 res = setup_bands_properties(sky); 651 if(res != RES_OK) goto error; 652 653 /* Load clouds properties */ 654 res = htcp_create(sky->logger, sky->allocator, sky->verbose, &sky->htcp); 655 if(res != RES_OK) { 656 log_err(sky, "Could not create the loader of cloud properties.\n"); 657 goto error; 658 } 659 res = htcp_load(sky->htcp, args->htcp_filename); 660 if(res != RES_OK) { 661 log_err(sky, "Error loading the cloud properties -- `%s'.\n", 662 args->htcp_filename); 663 goto error; 664 } 665 666 if(args->cache_filename) { 667 res = setup_cache_stream(sky, args->htcp_filename, args->htgop_filename, 668 args->htmie_filename, args->cache_filename, &force_cache_upd, &cache); 669 if(res != RES_OK) goto error; 670 } 671 672 time_current(&t0); 673 res = cloud_setup(sky, args->grid_max_definition, args->optical_thickness, 674 args->cache_filename, force_cache_upd, cache); 675 if(res != RES_OK) goto error; 676 time_sub(&t0, time_current(&t1), &t0); 677 time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf)); 678 log_info(sky, "Setup clouds in %s\n", buf); 679 680 if(sky->verbose) { 681 log_svx_memory_usage(sky); 682 } 683 684 exit: 685 if(cache) fclose(cache); 686 if(out_sky) *out_sky = sky; 687 return res; 688 error: 689 if(sky) { 690 htsky_ref_put(sky); 691 sky = NULL; 692 } 693 goto exit; 694 } 695 696 res_T 697 htsky_ref_get(struct htsky* sky) 698 { 699 if(!sky) return RES_BAD_ARG; 700 ref_get(&sky->ref); 701 return RES_OK; 702 } 703 704 res_T 705 htsky_ref_put(struct htsky* sky) 706 { 707 if(!sky) return RES_BAD_ARG; 708 ref_put(&sky->ref, release_sky); 709 return RES_OK; 710 } 711 712 const char* 713 htsky_get_name(const struct htsky* sky) 714 { 715 ASSERT(sky); 716 return str_cget(&sky->name); 717 } 718 719 double 720 htsky_fetch_particle_phase_function_asymmetry_parameter 721 (const struct htsky* sky, 722 const size_t ispectral_band, 723 const size_t iquad) 724 { 725 size_t i; 726 ASSERT(sky); 727 ASSERT(ispectral_band >= sky->bands_range[0]); 728 ASSERT(ispectral_band <= sky->bands_range[1]); 729 (void)iquad; 730 if(!sky->is_cloudy) { 731 return 0; 732 } else { 733 i = ispectral_band - sky->bands_range[0]; 734 return sky->bands[i].g_avg; 735 } 736 } 737 738 double 739 htsky_fetch_per_wavelength_particle_phase_function_asymmetry_parameter 740 (const struct htsky* sky, 741 const double wavelength) /* In nanometer */ 742 { 743 ASSERT(sky); 744 if(!sky->is_cloudy) { 745 return 0; 746 } else { 747 return htmie_fetch_asymmetry_parameter 748 (sky->htmie, wavelength, HTMIE_FILTER_LINEAR); 749 } 750 } 751 752 double 753 htsky_fetch_raw_property 754 (const struct htsky* sky, 755 const enum htsky_property prop, 756 const int components_mask, /* Combination of htsky_component_flag */ 757 const size_t iband, /* Index of the spectral band */ 758 const size_t iquad, /* Index of the quadrature point in the spectral band */ 759 const double pos[3], 760 const double k_min, 761 const double k_max) 762 { 763 size_t ivox[3]; 764 size_t i; 765 const struct svx_tree_desc* cloud_desc = NULL; 766 const struct svx_tree_desc* atmosphere_desc = NULL; 767 int comp_mask = components_mask; 768 int in_clouds; /* Defines if `pos' lies in the clouds */ 769 int in_atmosphere; /* Defines if `pos' lies in the atmosphere */ 770 double pos_cs[3]; /* Position in cloud space */ 771 double k_particle = 0; 772 double k_gas = 0; 773 double k = 0; 774 ASSERT(sky && pos); 775 ASSERT(iband >= sky->bands_range[0]); 776 ASSERT(iband <= sky->bands_range[1]); 777 ASSERT(comp_mask & HTSKY_CPNT_MASK_ALL); 778 779 i = iband - sky->bands_range[0]; 780 cloud_desc = sky->is_cloudy ? &sky->clouds[i][iquad].octree_desc : NULL; 781 atmosphere_desc = &sky->atmosphere[i][iquad].bitree_desc; 782 ASSERT(atmosphere_desc->frame[0] == SVX_AXIS_Z); 783 784 /* Is the position inside the clouds? */ 785 if(!sky->is_cloudy) { 786 in_clouds = 0; 787 } else if(sky->repeat_clouds) { 788 in_clouds = 789 pos[2] >= cloud_desc->lower[2] 790 && pos[2] <= cloud_desc->upper[2]; 791 } else { 792 in_clouds = 793 pos[0] >= cloud_desc->lower[0] 794 && pos[1] >= cloud_desc->lower[1] 795 && pos[2] >= cloud_desc->lower[2] 796 && pos[0] <= cloud_desc->upper[0] 797 && pos[1] <= cloud_desc->upper[1] 798 && pos[2] <= cloud_desc->upper[2]; 799 } 800 801 /* Is the position inside the atmosphere? */ 802 ASSERT(atmosphere_desc->frame[0] == SVX_AXIS_Z); 803 in_atmosphere = 804 pos[2] >= atmosphere_desc->lower[2] 805 && pos[2] <= atmosphere_desc->upper[2]; 806 807 if(!in_clouds) { 808 /* Make invalid the voxel index */ 809 ivox[0] = SIZE_MAX; 810 ivox[1] = SIZE_MAX; 811 ivox[2] = SIZE_MAX; 812 /* Not in clouds => No particle */ 813 comp_mask &= ~HTSKY_CPNT_FLAG_PARTICLES; 814 /* Not in atmopshere => No gas */ 815 if(!in_atmosphere) comp_mask &= ~HTSKY_CPNT_FLAG_GAS; 816 } else { 817 const double* upp; 818 const size_t* def; 819 world_to_cloud(sky, pos, pos_cs); 820 821 /* Compute the index of the voxel to fetch */ 822 ivox[0] = (size_t)((pos_cs[0] - cloud_desc->lower[0])/sky->htcp_desc.vxsz_x); 823 ivox[1] = (size_t)((pos_cs[1] - cloud_desc->lower[1])/sky->htcp_desc.vxsz_y); 824 if(!sky->htcp_desc.irregular_z) { 825 /* The voxels along the Z dimension have the same size */ 826 ivox[2] = (size_t)((pos_cs[2] - cloud_desc->lower[2])/sky->htcp_desc.vxsz_z[0]); 827 } else { 828 /* Irregular voxel size along the Z dimension. Compute the index of the Z 829 * position in the svx2htcp_z Look Up Table and use the LUT to define the 830 * voxel index into the HTCP descriptor */ 831 const struct split* splits = darray_split_cdata_get(&sky->svx2htcp_z); 832 const size_t nsplits = darray_split_size_get(&sky->svx2htcp_z); 833 const size_t ilut = (size_t) 834 ((pos_cs[2] - cloud_desc->lower[2]) / sky->lut_cell_sz); 835 if(ilut >= nsplits) FATAL("LUT index out of range\n"); 836 ivox[2] = splits[ilut].index + (pos_cs[2] > splits[ilut].height); 837 } 838 839 /* Handle numerical issues that may lead to a position lying onto the cloud 840 * upper boundaries */ 841 def = sky->htcp_desc.spatial_definition; 842 upp = cloud_desc->upper; 843 if(ivox[0] == def[0] && eq_eps(pos_cs[0], upp[0], 1.e-6)) ivox[0] = def[0]-1; 844 if(ivox[1] == def[1] && eq_eps(pos_cs[1], upp[1], 1.e-6)) ivox[1] = def[1]-1; 845 if(ivox[2] == def[2] && eq_eps(pos_cs[2], upp[2], 1.e-6)) ivox[2] = def[2]-1; 846 if(ivox[0] >= def[0]) FATAL("Out of bound X voxel coordinate\n"); 847 if(ivox[1] >= def[1]) FATAL("Out of bound Y voxel coordinate\n"); 848 if(ivox[2] >= def[2]) FATAL("Out of bound Z voxel coordinate\n"); 849 } 850 851 if(comp_mask & HTSKY_CPNT_FLAG_PARTICLES) { 852 k_particle = particle_fetch_raw_property(sky, prop, iband, iquad, ivox); 853 } 854 if(comp_mask & HTSKY_CPNT_FLAG_GAS) { 855 k_gas = gas_fetch_raw_property 856 (sky, prop, iband, iquad, in_clouds, pos, ivox); 857 } 858 k = k_particle + k_gas; 859 ASSERT(k >= k_min && k <= k_max); 860 (void)k_min, (void)k_max; 861 return k; 862 } 863 864 double 865 htsky_fetch_temperature(const struct htsky* sky, const double pos[3]) 866 { 867 double temperature = 0; 868 struct htgop_level lvl0, lvl1; 869 size_t nlvls = 0; 870 int in_clouds = 0; 871 int in_atmosphere = 0; 872 ASSERT(sky && pos); 873 874 /* Is the position inside the atmosphere */ 875 HTGOP(get_levels_count(sky->htgop, &nlvls)); 876 ASSERT(nlvls > 1); 877 HTGOP(get_level(sky->htgop, 0, &lvl0)); 878 HTGOP(get_level(sky->htgop, nlvls-1, &lvl1)); 879 in_atmosphere = pos[2] >= lvl0.height && pos[2] <= lvl1.height; 880 881 /* Is the position inside the clouds? */ 882 if(!sky->is_cloudy) { 883 in_clouds = 0; 884 } else if(sky->repeat_clouds) { 885 in_clouds = 886 pos[2] >= sky->htcp_desc.lower[2] 887 && pos[2] <= sky->htcp_desc.upper[2]; 888 } else { 889 in_clouds = 890 pos[0] >= sky->htcp_desc.lower[0] 891 && pos[1] >= sky->htcp_desc.lower[1] 892 && pos[2] >= sky->htcp_desc.lower[2] 893 && pos[0] <= sky->htcp_desc.upper[0] 894 && pos[1] <= sky->htcp_desc.upper[1] 895 && pos[2] <= sky->htcp_desc.upper[2]; 896 } 897 898 if(in_clouds) { 899 /* Clouds have priority on atmosphere */ 900 in_atmosphere = 0; 901 } 902 903 if(in_atmosphere) { 904 double u; 905 size_t ilayer; 906 907 /* Find the layer into which pos is included */ 908 HTGOP(position_to_layer_id(sky->htgop, pos[2], &ilayer)); 909 ASSERT(ilayer < nlvls-1); 910 911 /* Fetch the levels enclosing the current layer */ 912 HTGOP(get_level(sky->htgop, ilayer+0, &lvl0)); 913 HTGOP(get_level(sky->htgop, ilayer+1, &lvl1)); 914 ASSERT(lvl0.height < lvl1.height); 915 ASSERT(lvl0.height <= pos[2] && pos[2] <= lvl1.height); 916 917 /* Linearly interpolate the temperature of the levels into which pos lies */ 918 u = (pos[2] - lvl0.height) / (lvl1.height - lvl0.height); 919 temperature = u * (lvl1.temperature - lvl0.temperature) + lvl0.temperature; 920 921 } else if(in_clouds) { 922 const struct htcp_desc* desc = &sky->htcp_desc; 923 size_t ivox[3]; 924 double pos_cs[3]; 925 926 /* Transform the submitted position in local cloud space */ 927 world_to_cloud(sky, pos, pos_cs); 928 929 /* Compute the index of the voxel to fetch */ 930 ivox[0] = (size_t)((pos_cs[0] - desc->lower[0])/desc->vxsz_x); 931 ivox[1] = (size_t)((pos_cs[1] - desc->lower[1])/desc->vxsz_y); 932 if(!desc->irregular_z) { 933 /* The voxels along the Z dimension have the same size */ 934 ivox[2] = (size_t)((pos_cs[2] - desc->lower[2])/desc->vxsz_z[0]); 935 } else { 936 /* Irregular voxel size along the Z dimension. Compute the index of the Z 937 * position in the svx2htcp_z Look Up Table and use the LUT to define the 938 * voxel index into the HTCP descripptor */ 939 const struct split* splits = darray_split_cdata_get(&sky->svx2htcp_z); 940 const size_t nsplits = darray_split_size_get(&sky->svx2htcp_z); 941 const size_t ilut = (size_t)((pos_cs[2] - desc->lower[2]) / sky->lut_cell_sz); 942 if(ilut >= nsplits) FATAL("LUT index out of range\n"); 943 ivox[2] = splits[ilut].index + (pos_cs[2] > splits[ilut].height); 944 } 945 946 /* Fetch the cloud temperature */ 947 temperature = htcp_desc_T_at(desc, ivox[0], ivox[1], ivox[2], 0); 948 } 949 950 return temperature; 951 } 952 953 size_t 954 htsky_get_spectral_bands_count(const struct htsky* sky) 955 { 956 ASSERT(sky && sky->bands_range[0] <= sky->bands_range[1]); 957 return sky->bands_range[1] - sky->bands_range[0] + 1; 958 } 959 960 size_t 961 htsky_get_spectral_band_id 962 (const struct htsky* sky, const size_t i) 963 { 964 ASSERT(sky); 965 ASSERT(i < htsky_get_spectral_bands_count(sky)); 966 return sky->bands_range[0] + i; 967 } 968 969 size_t 970 htsky_get_spectral_band_quadrature_length 971 (const struct htsky* sky, const size_t iband) 972 { 973 struct htgop_spectral_interval band; 974 ASSERT(sky); 975 ASSERT(iband >= sky->bands_range[0]); 976 ASSERT(iband <= sky->bands_range[1]); 977 switch(sky->spectral_type) { 978 case HTSKY_SPECTRAL_LW: 979 HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); 980 break; 981 case HTSKY_SPECTRAL_SW: 982 HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); 983 break; 984 default: FATAL("Unreachable code.\n"); break; 985 } 986 return band.quadrature_length; 987 } 988 989 res_T 990 htsky_get_spectral_band_bounds 991 (const struct htsky* sky, 992 const size_t iband, 993 double wavelengths[2]) 994 { 995 struct htgop_spectral_interval specint; 996 res_T res = RES_OK; 997 ASSERT(sky && wavelengths); 998 999 switch(sky->spectral_type) { 1000 case HTSKY_SPECTRAL_LW: 1001 res = htgop_get_lw_spectral_interval(sky->htgop, iband, &specint); 1002 break; 1003 case HTSKY_SPECTRAL_SW: 1004 res = htgop_get_sw_spectral_interval(sky->htgop, iband, &specint); 1005 break; 1006 default: FATAL("Unreachable code.\n"); break; 1007 } 1008 if(res != RES_OK) return res; 1009 wavelengths[0] = wavenumber_to_wavelength(specint.wave_numbers[1]); 1010 wavelengths[1] = wavenumber_to_wavelength(specint.wave_numbers[0]); 1011 ASSERT(wavelengths[0] < wavelengths[1]); 1012 return RES_OK; 1013 } 1014 1015 res_T 1016 htsky_get_raw_spectral_bounds(const struct htsky* sky, double wavelengths[2]) 1017 { 1018 size_t n; 1019 double band_first_bounds[2]; 1020 double band_last_bounds[2]; 1021 size_t iband_first; 1022 size_t iband_last; 1023 ASSERT(sky && wavelengths); 1024 1025 n = htsky_get_spectral_bands_count(sky); 1026 1027 iband_first = htsky_get_spectral_band_id(sky, 0); 1028 iband_last = htsky_get_spectral_band_id(sky, n-1); 1029 1030 HTSKY(get_spectral_band_bounds(sky, iband_first, band_first_bounds)); 1031 HTSKY(get_spectral_band_bounds(sky, iband_last, band_last_bounds)); 1032 wavelengths[0] = MMIN(band_first_bounds[0], band_last_bounds[0]); 1033 wavelengths[1] = MMAX(band_first_bounds[1], band_last_bounds[1]); 1034 1035 return RES_OK; 1036 } 1037 1038 enum htsky_spectral_type 1039 htsky_get_spectral_type(const struct htsky* htsky) 1040 { 1041 ASSERT(htsky); 1042 return htsky->spectral_type; 1043 } 1044 1045 size_t 1046 htsky_find_spectral_band(const struct htsky* sky, const double wavelength) 1047 { 1048 const double wnum = wavelength_to_wavenumber(wavelength); 1049 size_t iband; 1050 ASSERT(sky); 1051 switch(sky->spectral_type) { 1052 case HTSKY_SPECTRAL_LW: 1053 HTGOP(find_lw_spectral_interval_id(sky->htgop, wnum, &iband)); 1054 break; 1055 case HTSKY_SPECTRAL_SW: 1056 HTGOP(find_sw_spectral_interval_id(sky->htgop, wnum, &iband)); 1057 break; 1058 default: FATAL("Unreachable code.\n"); break; 1059 } 1060 return iband; 1061 } 1062 1063 size_t 1064 htsky_spectral_band_sample_quadrature 1065 (const struct htsky* sky, 1066 const double r, 1067 const size_t iband) 1068 { 1069 struct htgop_spectral_interval band; 1070 size_t iquad; 1071 ASSERT(sky); 1072 ASSERT(sky->bands_range[0] <= iband || iband <= sky->bands_range[1]); 1073 1074 switch(sky->spectral_type) { 1075 case HTSKY_SPECTRAL_LW: 1076 HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); 1077 break; 1078 case HTSKY_SPECTRAL_SW: 1079 HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); 1080 break; 1081 default: FATAL("Unreachable code.\n"); break; 1082 } 1083 HTGOP(spectral_interval_sample_quadrature(&band, r, &iquad)); 1084 return iquad; 1085 } 1086 1087 /******************************************************************************* 1088 * Local functions 1089 ******************************************************************************/ 1090 double* 1091 world_to_cloud 1092 (const struct htsky* sky, 1093 const double pos_ws[3], /* World space position */ 1094 double out_pos_cs[3]) 1095 { 1096 double cloud_sz[2]; 1097 double upper[2]; 1098 double pos_cs[3]; 1099 double pos_cs_n[2]; 1100 ASSERT(sky && pos_ws && out_pos_cs); 1101 ASSERT(pos_ws[2] >= sky->htcp_desc.lower[2]); 1102 ASSERT(pos_ws[2] <= sky->htcp_desc.upper[2]); 1103 1104 if(!sky->repeat_clouds) { /* Nothing to do */ 1105 return d3_set(out_pos_cs, pos_ws); 1106 } 1107 1108 if(!sky->repeat_clouds /* Nothing to do */ 1109 || ( pos_ws[0] >= sky->htcp_desc.lower[0] 1110 && pos_ws[0] <= sky->htcp_desc.upper[0] 1111 && pos_ws[1] >= sky->htcp_desc.lower[1] 1112 && pos_ws[1] <= sky->htcp_desc.upper[1])) { 1113 return d3_set(out_pos_cs, pos_ws); 1114 } 1115 1116 /* The cloud upper bound is not inclusive. Define the inclusive upper bound 1117 * of the cloud */ 1118 upper[0] = nextafter(sky->htcp_desc.upper[0], sky->htcp_desc.lower[0]); 1119 upper[1] = nextafter(sky->htcp_desc.upper[1], sky->htcp_desc.lower[1]); 1120 cloud_sz[0] = upper[0] - sky->htcp_desc.lower[0]; 1121 cloud_sz[1] = upper[1] - sky->htcp_desc.lower[1]; 1122 1123 /* Transform pos in normalize local cloud space */ 1124 pos_cs_n[0] = (pos_ws[0] - sky->htcp_desc.lower[0]) / cloud_sz[0]; 1125 pos_cs_n[1] = (pos_ws[1] - sky->htcp_desc.lower[1]) / cloud_sz[1]; 1126 pos_cs_n[0] -= (int)pos_cs_n[0]; /* Get fractional part */ 1127 pos_cs_n[1] -= (int)pos_cs_n[1]; /* Get fractional part */ 1128 if(pos_cs_n[0] < 0) pos_cs_n[0] += 1; 1129 if(pos_cs_n[1] < 0) pos_cs_n[1] += 1; 1130 1131 /* Transform pos in local cloud space */ 1132 pos_cs[0] = sky->htcp_desc.lower[0] + pos_cs_n[0] * cloud_sz[0]; 1133 pos_cs[1] = sky->htcp_desc.lower[1] + pos_cs_n[1] * cloud_sz[1]; 1134 pos_cs[2] = pos_ws[2]; 1135 1136 ASSERT(pos_cs[0] >= sky->htcp_desc.lower[0]); 1137 ASSERT(pos_cs[0] <= sky->htcp_desc.upper[0]); 1138 ASSERT(pos_cs[1] >= sky->htcp_desc.lower[1]); 1139 ASSERT(pos_cs[1] <= sky->htcp_desc.upper[1]); 1140 1141 return d3_set(out_pos_cs, pos_cs); 1142 } 1143