htrdr_atmosphere_draw_map.c (20895B)
1 /* Copyright (C) 2018-2019, 2022-2025 Centre National de la Recherche Scientifique 2 * Copyright (C) 2020-2022 Institut Mines Télécom Albi-Carmaux 3 * Copyright (C) 2022-2025 Institut Pierre-Simon Laplace 4 * Copyright (C) 2022-2025 Institut de Physique du Globe de Paris 5 * Copyright (C) 2018-2025 |Méso|Star> (contact@meso-star.com) 6 * Copyright (C) 2022-2025 Observatoire de Paris 7 * Copyright (C) 2022-2025 Université de Reims Champagne-Ardenne 8 * Copyright (C) 2022-2025 Université de Versaille Saint-Quentin 9 * Copyright (C) 2018-2019, 2022-2025 Université Paul Sabatier 10 * 11 * This program is free software: you can redistribute it and/or modify 12 * it under the terms of the GNU General Public License as published by 13 * the Free Software Foundation, either version 3 of the License, or 14 * (at your option) any later version. 15 * 16 * This program is distributed in the hope that it will be useful, 17 * but WITHOUT ANY WARRANTY; without even the implied warranty of 18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 * GNU General Public License for more details. 20 * 21 * You should have received a copy of the GNU General Public License 22 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 23 24 #include "atmosphere/htrdr_atmosphere_c.h" 25 #include "atmosphere/htrdr_atmosphere_ground.h" 26 #include "atmosphere/htrdr_atmosphere_sun.h" 27 28 #include "core/htrdr.h" 29 #include "core/htrdr_buffer.h" 30 #include "core/htrdr_draw_map.h" 31 #include "core/htrdr_interface.h" 32 #include "core/htrdr_log.h" 33 #include "core/htrdr_ran_wlen_cie_xyz.h" 34 #include "core/htrdr_ran_wlen_planck.h" 35 #include "core/htrdr_rectangle.h" 36 37 #include <high_tune/htsky.h> 38 39 #include <star/s3d.h> 40 #include <star/scam.h> 41 #include <star/ssp.h> 42 43 #include <rsys/clock_time.h> 44 #include <rsys/str.h> 45 46 #include <string.h> 47 48 /******************************************************************************* 49 * Helper functions 50 ******************************************************************************/ 51 static res_T 52 sample_rectangle_ray 53 (struct htrdr_atmosphere* cmd, 54 struct htrdr_rectangle* rect, 55 const size_t ipix[2], 56 const double pix_sz[2], 57 struct ssp_rng* rng, 58 double ray_org[3], 59 double ray_dir[3]) 60 { 61 struct s3d_hit hit = S3D_HIT_NULL; 62 double pix_samp[2]; 63 const double up_dir[3] = {0,0,1}; 64 const double range[2] = {0, DBL_MAX}; 65 double normal[3]; 66 ASSERT(cmd && rect && ipix && pix_sz && rng && ray_org && ray_dir); 67 68 /* Sample a position into the pixel, in the normalized image plane */ 69 pix_samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0]; 70 pix_samp[1] = ((double)ipix[1] + ssp_rng_canonical(rng)) * pix_sz[1]; 71 72 /* Retrieve the world space position of pix_samp */ 73 htrdr_rectangle_sample_pos(rect, pix_samp, ray_org); 74 75 /* Check that `ray_org' is not included into a geometry */ 76 HTRDR(atmosphere_ground_trace_ray 77 (cmd->ground, ray_org, up_dir, range, NULL, &hit)); 78 79 /* Up direction is occluded. Check if the sample must be rejected, i.e. does it 80 * lies inside a geometry? */ 81 if(!S3D_HIT_NONE(&hit)) { 82 struct htrdr_interface interf = HTRDR_INTERFACE_NULL; 83 const struct htrdr_mtl* mtl = NULL; 84 float N[3] = {0,0,0}; /* Normalized normal of the hit */ 85 float wi[3]; 86 float cos_wi_N; 87 88 /* Compute the cosine between the up direction and the hit normal */ 89 f3_set_d3(wi, up_dir); 90 f3_normalize(N, hit.normal); 91 cos_wi_N = f3_dot(wi, N); 92 93 /* Fetch the hit interface and retrieve the material into which the ray was 94 * traced */ 95 htrdr_atmosphere_ground_get_interface(cmd->ground, &hit, &interf); 96 mtl = cos_wi_N < 0 ? &interf.mtl_front : &interf.mtl_back; 97 98 /* Reject the sample if the incident direction do not travel into the sky */ 99 if(strcmp(mtl->name, cmd->sky_mtl_name) != 0) return RES_BAD_OP; 100 } 101 102 /* Sample a ray direction */ 103 htrdr_rectangle_get_normal(rect, normal); 104 ssp_ran_hemisphere_cos(rng, normal, ray_dir, NULL); 105 106 return RES_OK; 107 } 108 109 static void 110 draw_pixel_image 111 (struct htrdr* htrdr, 112 const struct htrdr_draw_pixel_args* args, 113 void* data) 114 { 115 struct htrdr_accum XYZ[3]; /* X, Y, and Z */ 116 struct htrdr_accum time; 117 struct htrdr_atmosphere* cmd; 118 struct atmosphere_pixel_image* pixel = data; 119 size_t ichannel; 120 ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data); 121 (void)htrdr; 122 123 cmd = args->context; 124 ASSERT(cmd); 125 ASSERT(cmd->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ); 126 ASSERT(cmd->output_type == HTRDR_ATMOSPHERE_ARGS_OUTPUT_IMAGE); 127 128 /* Reset accumulators */ 129 XYZ[0] = HTRDR_ACCUM_NULL; 130 XYZ[1] = HTRDR_ACCUM_NULL; 131 XYZ[2] = HTRDR_ACCUM_NULL; 132 time = HTRDR_ACCUM_NULL; 133 134 FOR_EACH(ichannel, 0, 3) { 135 size_t isamp; 136 137 FOR_EACH(isamp, 0, args->spp) { 138 struct time t0, t1; 139 struct scam_sample sample = SCAM_SAMPLE_NULL; 140 struct scam_ray ray = SCAM_RAY_NULL; 141 double weight; 142 double r0, r1, r2; 143 double wlen; /* Sampled wavelength into the spectral band */ 144 double pdf; 145 size_t iband; /* Sampled spectral band */ 146 size_t iquad; /* Sampled quadrature point into the spectral band */ 147 double usec; 148 149 /* Begin the registration of the time spent to in the realisation */ 150 time_current(&t0); 151 152 /* Sample a position into the pixel, in the normalized image plane */ 153 sample.film[0] = (double)args->pixel_coord[0]+ssp_rng_canonical(args->rng); 154 sample.film[1] = (double)args->pixel_coord[1]+ssp_rng_canonical(args->rng); 155 sample.film[0] *= args->pixel_normalized_size[0]; 156 sample.film[1] *= args->pixel_normalized_size[1]; 157 sample.lens[0] = ssp_rng_canonical(args->rng); 158 sample.lens[1] = ssp_rng_canonical(args->rng); 159 160 /* Generate a camera ray */ 161 scam_generate_ray(cmd->camera, &sample, &ray); 162 163 r0 = ssp_rng_canonical(args->rng); 164 r1 = ssp_rng_canonical(args->rng); 165 r2 = ssp_rng_canonical(args->rng); 166 167 /* Sample a spectral band and a quadrature point */ 168 switch(ichannel) { 169 case 0: wlen = htrdr_ran_wlen_cie_xyz_sample_X(cmd->cie, r0, r1, &pdf); break; 170 case 1: wlen = htrdr_ran_wlen_cie_xyz_sample_Y(cmd->cie, r0, r1, &pdf); break; 171 case 2: wlen = htrdr_ran_wlen_cie_xyz_sample_Z(cmd->cie, r0, r1, &pdf); break; 172 default: FATAL("Unreachable code.\n"); break; 173 } 174 pdf *= 1.e9; /* Transform the pdf from nm^-1 to m^-1 */ 175 176 iband = htsky_find_spectral_band(cmd->sky, wlen); 177 iquad = htsky_spectral_band_sample_quadrature(cmd->sky, r2, iband); 178 179 /* Compute the radiance in W/m^2/sr/m */ 180 weight = atmosphere_compute_radiance_sw(cmd, args->ithread, args->rng, 181 ATMOSPHERE_RADIANCE_ALL, ray.org, ray.dir, wlen, iband, iquad); 182 ASSERT(weight >= 0); 183 184 weight /= pdf; /* In W/m^2/sr */ 185 186 /* End the registration of the per realisation time */ 187 time_sub(&t0, time_current(&t1), &t0); 188 usec = (double)time_val(&t0, TIME_NSEC) * 0.001; 189 190 /* Update the pixel accumulator of the current channel */ 191 XYZ[ichannel].sum_weights += weight; 192 XYZ[ichannel].sum_weights_sqr += weight*weight; 193 XYZ[ichannel].nweights += 1; 194 195 /* Update the pixel accumulator of per realisation time */ 196 time.sum_weights += usec; 197 time.sum_weights_sqr += usec*usec; 198 time.nweights += 1; 199 } 200 } 201 202 /* Flush pixel data */ 203 htrdr_accum_get_estimation(XYZ+0, &pixel->X); 204 htrdr_accum_get_estimation(XYZ+1, &pixel->Y); 205 htrdr_accum_get_estimation(XYZ+2, &pixel->Z); 206 pixel->time = time; 207 } 208 209 static void 210 draw_pixel_flux 211 (struct htrdr* htrdr, 212 const struct htrdr_draw_pixel_args* args, 213 void* data) 214 { 215 struct htrdr_accum flux; 216 struct htrdr_accum time; 217 struct htrdr_atmosphere* cmd; 218 struct atmosphere_pixel_flux* pixel = data; 219 size_t isamp; 220 ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data); 221 (void)htrdr; 222 223 cmd = args->context; 224 ASSERT(cmd); 225 ASSERT(cmd->output_type == HTRDR_ATMOSPHERE_ARGS_OUTPUT_FLUX_MAP); 226 ASSERT(cmd->spectral_type == HTRDR_SPECTRAL_LW 227 || cmd->spectral_type == HTRDR_SPECTRAL_SW); 228 229 /* Reset the pixel accumulators */ 230 flux = HTRDR_ACCUM_NULL; 231 time = HTRDR_ACCUM_NULL; 232 233 FOR_EACH(isamp, 0, args->spp) { 234 struct time t0, t1; 235 double ray_org[3]; 236 double ray_dir[3]; 237 double weight; 238 double r0, r1, r2; 239 double wlen; 240 size_t iband; 241 size_t iquad; 242 double usec; 243 double band_pdf; 244 res_T res = RES_OK; 245 246 /* Begin the registration of the time spent in the realisation */ 247 time_current(&t0); 248 249 res = sample_rectangle_ray(cmd, cmd->flux_map, args->pixel_coord, 250 args->pixel_normalized_size, args->rng, ray_org, ray_dir); 251 if(res != RES_OK) continue; /* Reject the current sample */ 252 253 r0 = ssp_rng_canonical(args->rng); 254 r1 = ssp_rng_canonical(args->rng); 255 r2 = ssp_rng_canonical(args->rng); 256 257 /* Sample a wavelength */ 258 wlen = htrdr_ran_wlen_planck_sample(cmd->planck, r0, r1, &band_pdf); 259 band_pdf *= 1.e9; /* Transform the pdf from nm^-1 to m^-1 */ 260 261 /* Select the associated band and sample a quadrature point */ 262 iband = htsky_find_spectral_band(cmd->sky, wlen); 263 iquad = htsky_spectral_band_sample_quadrature(cmd->sky, r2, iband); 264 265 if(cmd->spectral_type == HTRDR_SPECTRAL_LW) { 266 weight = atmosphere_compute_radiance_lw(cmd, args->ithread, args->rng, 267 ray_org, ray_dir, wlen, iband, iquad); 268 weight *= PI / band_pdf; /* Transform weight from W/m^2/sr/m to W/m^2 */ 269 } else { 270 double sun_dir[3]; 271 double N[3]; 272 double L_direct; 273 double L_diffuse; 274 double cos_N_sun_dir; 275 double sun_solid_angle; 276 ASSERT(cmd->spectral_type == HTRDR_SPECTRAL_SW); 277 278 /* Compute direct contribution if necessary */ 279 htrdr_atmosphere_sun_sample_direction(cmd->sun, args->rng, sun_dir); 280 htrdr_rectangle_get_normal(cmd->flux_map, N); 281 cos_N_sun_dir = d3_dot(N, sun_dir); 282 283 if(cos_N_sun_dir <= 0) { 284 L_direct = 0; 285 } else { 286 L_direct = atmosphere_compute_radiance_sw(cmd, args->ithread, 287 args->rng, ATMOSPHERE_RADIANCE_DIRECT, ray_org, sun_dir, wlen, iband, 288 iquad); 289 } 290 291 /* Compute diffuse contribution */ 292 L_diffuse = atmosphere_compute_radiance_sw(cmd, args->ithread, args->rng, 293 ATMOSPHERE_RADIANCE_DIFFUSE, ray_org, ray_dir, wlen, iband, iquad); 294 295 sun_solid_angle = htrdr_atmosphere_sun_get_solid_angle(cmd->sun); 296 297 /* Compute the weight in W/m^2/m */ 298 weight = cos_N_sun_dir * sun_solid_angle * L_direct + PI * L_diffuse; 299 300 /* Importance sampling: correct weight with pdf */ 301 weight /= band_pdf; /* In W/m^2 */ 302 } 303 304 /* End the registration of the per realisation time */ 305 time_sub(&t0, time_current(&t1), &t0); 306 usec = (double)time_val(&t0, TIME_NSEC) * 0.001; 307 308 /* Update the pixel accumulator of the flux */ 309 flux.sum_weights += weight; 310 flux.sum_weights_sqr += weight*weight; 311 flux.nweights += 1; 312 313 /* Update the pixel accumulator of per realisation time */ 314 time.sum_weights += usec; 315 time.sum_weights_sqr += usec*usec; 316 time.nweights += 1; 317 } 318 319 /* Write the accumulators */ 320 pixel->flux = flux; 321 pixel->time = time; 322 } 323 324 static void 325 draw_pixel_xwave 326 (struct htrdr* htrdr, 327 const struct htrdr_draw_pixel_args* args, 328 void* data) 329 { 330 struct htrdr_accum radiance; 331 struct htrdr_accum time; 332 struct htrdr_atmosphere* cmd; 333 struct atmosphere_pixel_xwave* pixel = data; 334 size_t isamp; 335 double temp_min, temp_max; 336 ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data); 337 (void)htrdr; 338 339 cmd = args->context; 340 ASSERT(cmd->output_type == HTRDR_ATMOSPHERE_ARGS_OUTPUT_IMAGE); 341 ASSERT(cmd->spectral_type == HTRDR_SPECTRAL_LW 342 || cmd->spectral_type == HTRDR_SPECTRAL_SW); 343 344 /* Reset the pixel accumulators */ 345 radiance = HTRDR_ACCUM_NULL; 346 time = HTRDR_ACCUM_NULL; 347 348 FOR_EACH(isamp, 0, args->spp) { 349 struct time t0, t1; 350 struct scam_sample sample = SCAM_SAMPLE_NULL; 351 struct scam_ray ray = SCAM_RAY_NULL; 352 double weight; 353 double r0, r1, r2; 354 double wlen; 355 size_t iband; 356 size_t iquad; 357 double usec; 358 double band_pdf; 359 360 /* Begin the registration of the time spent in the realisation */ 361 time_current(&t0); 362 363 /* Sample a position into the pixel, in the normalized image plane */ 364 sample.film[0] = (double)args->pixel_coord[0]+ssp_rng_canonical(args->rng); 365 sample.film[1] = (double)args->pixel_coord[1]+ssp_rng_canonical(args->rng); 366 sample.film[0] *= args->pixel_normalized_size[0]; 367 sample.film[1] *= args->pixel_normalized_size[1]; 368 sample.lens[0] = ssp_rng_canonical(args->rng); 369 sample.lens[1] = ssp_rng_canonical(args->rng); 370 371 /* Generate a camera ray */ 372 scam_generate_ray(cmd->camera, &sample, &ray); 373 374 r0 = ssp_rng_canonical(args->rng); 375 r1 = ssp_rng_canonical(args->rng); 376 r2 = ssp_rng_canonical(args->rng); 377 378 /* Sample a wavelength */ 379 wlen = htrdr_ran_wlen_planck_sample(cmd->planck, r0, r1, &band_pdf); 380 band_pdf *= 1.e9; /* Transform the pdf from nm^-1 to m^-1 */ 381 382 /* Select the associated band and sample a quadrature point */ 383 iband = htsky_find_spectral_band(cmd->sky, wlen); 384 iquad = htsky_spectral_band_sample_quadrature(cmd->sky, r2, iband); 385 386 /* Compute the spectral radiance in W/m^2/sr/m */ 387 switch(cmd->spectral_type) { 388 case HTRDR_SPECTRAL_LW: 389 weight = atmosphere_compute_radiance_lw(cmd, args->ithread, args->rng, 390 ray.org, ray.dir, wlen, iband, iquad); 391 break; 392 case HTRDR_SPECTRAL_SW: 393 weight = atmosphere_compute_radiance_sw(cmd, args->ithread, args->rng, 394 ATMOSPHERE_RADIANCE_ALL, ray.org, ray.dir, wlen, iband, iquad); 395 break; 396 default: FATAL("Unreachable code.\n"); break; 397 } 398 ASSERT(weight >= 0); 399 /* Importance sampling: correct weight with pdf */ 400 weight /= band_pdf; /* In W/m^2/sr */ 401 402 /* End the registration of the per realisation time */ 403 time_sub(&t0, time_current(&t1), &t0); 404 usec = (double)time_val(&t0, TIME_NSEC) * 0.001; 405 406 /* Update the pixel accumulator */ 407 radiance.sum_weights += weight; 408 radiance.sum_weights_sqr += weight*weight; 409 radiance.nweights += 1; 410 411 /* Update the pixel accumulator of per realisation time */ 412 time.sum_weights += usec; 413 time.sum_weights_sqr += usec*usec; 414 time.nweights += 1; 415 } 416 417 /* Compute the estimation of the pixel radiance */ 418 htrdr_accum_get_estimation(&radiance, &pixel->radiance); 419 420 /* Save the per realisation integration time */ 421 pixel->time = time; 422 423 /* Compute the brightness_temperature of the pixel and estimate its standard 424 * error if the sources were in the medium (<=> longwave) */ 425 if(cmd->spectral_type == HTRDR_SPECTRAL_LW) { 426 const double wlen_min = cmd->wlen_range_m[0]; 427 const double wlen_max = cmd->wlen_range_m[1]; 428 pixel->radiance_temperature.E = htrdr_radiance_temperature 429 (cmd->htrdr, wlen_min, wlen_max, pixel->radiance.E); 430 temp_min = htrdr_radiance_temperature 431 (cmd->htrdr, wlen_min, wlen_max, pixel->radiance.E - pixel->radiance.SE); 432 temp_max = htrdr_radiance_temperature 433 (cmd->htrdr, wlen_min, wlen_max, pixel->radiance.E + pixel->radiance.SE); 434 pixel->radiance_temperature.SE = temp_max - temp_min; 435 } 436 } 437 438 static INLINE void 439 setup_draw_map_args_rectangle 440 (struct htrdr_atmosphere* cmd, 441 struct htrdr_draw_map_args* args) 442 { 443 ASSERT(cmd && args); 444 ASSERT(cmd->output_type == HTRDR_ATMOSPHERE_ARGS_OUTPUT_FLUX_MAP); 445 *args = HTRDR_DRAW_MAP_ARGS_NULL; 446 args->draw_pixel = draw_pixel_flux; 447 args->buffer_layout = cmd->buf_layout; 448 args->spp = cmd->spp; 449 args->context = cmd; 450 } 451 452 static INLINE void 453 setup_draw_map_args_camera 454 (struct htrdr_atmosphere* cmd, 455 struct htrdr_draw_map_args* args) 456 { 457 ASSERT(cmd && args); 458 ASSERT(cmd->output_type == HTRDR_ATMOSPHERE_ARGS_OUTPUT_IMAGE); 459 460 *args = HTRDR_DRAW_MAP_ARGS_NULL; 461 args->buffer_layout = cmd->buf_layout; 462 args->spp = cmd->spp; 463 args->context = cmd; 464 465 switch(cmd->spectral_type) { 466 case HTRDR_SPECTRAL_LW: 467 case HTRDR_SPECTRAL_SW: 468 args->draw_pixel = draw_pixel_xwave; 469 break; 470 case HTRDR_SPECTRAL_SW_CIE_XYZ: 471 args->draw_pixel = draw_pixel_image; 472 break; 473 default: FATAL("Unreachable code.\n"); break; 474 } 475 } 476 477 static INLINE void 478 dump_accum 479 (const struct htrdr_accum* acc, /* Accum to dump */ 480 struct htrdr_accum* out_acc, /* May be NULL */ 481 FILE* stream) 482 { 483 ASSERT(acc && stream); 484 485 if(acc->nweights == 0) { 486 fprintf(stream, "0 0 "); 487 } else { 488 struct htrdr_estimate estimate = HTRDR_ESTIMATE_NULL; 489 490 htrdr_accum_get_estimation(acc, &estimate); 491 fprintf(stream, "%g %g ", estimate.E, estimate.SE); 492 493 if(out_acc) { 494 out_acc->sum_weights += acc->sum_weights; 495 out_acc->sum_weights_sqr += acc->sum_weights_sqr; 496 out_acc->nweights += acc->nweights; 497 } 498 } 499 } 500 501 static INLINE void 502 dump_pixel_flux 503 (const struct atmosphere_pixel_flux* pix, 504 struct htrdr_accum* time_acc, /* May be NULL */ 505 struct htrdr_accum* flux_acc, /* May be NULL */ 506 FILE* stream) 507 { 508 ASSERT(pix && stream); 509 dump_accum(&pix->flux, flux_acc, stream); 510 fprintf(stream, "0 0 0 0 "); 511 dump_accum(&pix->time, time_acc, stream); 512 fprintf(stream, "\n"); 513 } 514 515 static INLINE void 516 dump_pixel_image 517 (const struct atmosphere_pixel_image* pix, 518 struct htrdr_accum* time_acc, /* May be NULL */ 519 FILE* stream) 520 { 521 ASSERT(pix && stream); 522 fprintf(stream, "%g %g ", pix->X.E, pix->X.SE); 523 fprintf(stream, "%g %g ", pix->Y.E, pix->Y.SE); 524 fprintf(stream, "%g %g ", pix->Z.E, pix->Z.SE); 525 dump_accum(&pix->time, time_acc, stream); 526 fprintf(stream, "\n"); 527 } 528 529 static INLINE void 530 dump_pixel_xwave 531 (const struct atmosphere_pixel_xwave* pix, 532 struct htrdr_accum* time_acc, /* May be NULL */ 533 FILE* stream) 534 { 535 ASSERT(pix && stream); 536 fprintf(stream, "%g %g %f %f 0 0 ", 537 pix->radiance_temperature.E, 538 pix->radiance_temperature.SE, 539 pix->radiance.E, 540 pix->radiance.SE); 541 dump_accum(&pix->time, time_acc, stream); 542 fprintf(stream, "\n"); 543 } 544 545 static res_T 546 dump_buffer 547 (struct htrdr_atmosphere* cmd, 548 struct htrdr_buffer* buf, 549 struct htrdr_accum* time_acc, /* May be NULL */ 550 struct htrdr_accum* flux_acc, /* May be NULL */ 551 FILE* stream) 552 { 553 struct htrdr_pixel_format pixfmt; 554 struct htrdr_buffer_layout layout; 555 size_t x, y; 556 ASSERT(cmd && buf && stream); 557 558 atmosphere_get_pixel_format(cmd, &pixfmt); 559 htrdr_buffer_get_layout(buf, &layout); 560 CHK(pixfmt.size == layout.elmt_size); 561 562 fprintf(stream, "%lu %lu\n", layout.width, layout.height); 563 564 if(time_acc) *time_acc = HTRDR_ACCUM_NULL; 565 if(flux_acc) *flux_acc = HTRDR_ACCUM_NULL; 566 567 FOR_EACH(y, 0, layout.height) { 568 FOR_EACH(x, 0, layout.width) { 569 void* pix_raw = htrdr_buffer_at(buf, x, y); 570 ASSERT(IS_ALIGNED(pix_raw, pixfmt.alignment)); 571 572 if(cmd->output_type == HTRDR_ATMOSPHERE_ARGS_OUTPUT_FLUX_MAP) { 573 const struct atmosphere_pixel_flux* pix = pix_raw; 574 dump_pixel_flux(pix, time_acc, flux_acc, stream); 575 } else if(cmd->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ) { 576 const struct atmosphere_pixel_image* pix = pix_raw; 577 dump_pixel_image(pix, time_acc, stream); 578 } else { 579 const struct atmosphere_pixel_xwave* pix = pix_raw; 580 dump_pixel_xwave(pix, time_acc, stream); 581 } 582 } 583 fprintf(stream, "\n"); 584 } 585 586 return RES_OK; 587 } 588 589 /******************************************************************************* 590 * Local functions 591 ******************************************************************************/ 592 res_T 593 atmosphere_draw_map(struct htrdr_atmosphere* cmd) 594 { 595 struct htrdr_draw_map_args args = HTRDR_DRAW_MAP_ARGS_NULL; 596 struct htrdr_accum path_time_acc = HTRDR_ACCUM_NULL; 597 struct htrdr_accum flux_acc = HTRDR_ACCUM_NULL; 598 struct htrdr_estimate path_time; 599 struct htrdr_estimate flux; 600 res_T res = RES_OK; 601 ASSERT(cmd); 602 603 args.spp = cmd->spp; 604 605 switch(cmd->output_type) { 606 case HTRDR_ATMOSPHERE_ARGS_OUTPUT_FLUX_MAP: 607 setup_draw_map_args_rectangle(cmd, &args); 608 break; 609 case HTRDR_ATMOSPHERE_ARGS_OUTPUT_IMAGE: 610 setup_draw_map_args_camera(cmd, &args); 611 break; 612 default: FATAL("Unreachable code.\n"); break; 613 } 614 615 res = htrdr_draw_map(cmd->htrdr, &args, cmd->buf); 616 if(res != RES_OK) goto error; 617 618 /* No more to do on non master processes */ 619 if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit; 620 621 /* Write buffer to output */ 622 res = dump_buffer(cmd, cmd->buf, &path_time_acc, &flux_acc, cmd->output); 623 if(res != RES_OK) goto error; 624 625 htrdr_accum_get_estimation(&path_time_acc, &path_time); 626 htrdr_log(cmd->htrdr, 627 "Time per radiative path (in micro seconds): %g +/- %g\n", 628 path_time.E, path_time.SE); 629 630 if(cmd->output_type == HTRDR_ATMOSPHERE_ARGS_OUTPUT_FLUX_MAP) { 631 htrdr_accum_get_estimation(&flux_acc, &flux); 632 htrdr_log(cmd->htrdr, 633 "Radiative flux density (in W/(external m^2)): %g +/- %g\n", 634 flux.E, flux.SE); 635 } 636 637 exit: 638 return res; 639 error: 640 goto exit; 641 }