htrdr_planets_compute_radiance.c (21306B)
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 "planets/htrdr_planets_c.h" 25 #include "planets/htrdr_planets_source.h" 26 27 #include <rad-net/rnatm.h> 28 #include <rad-net/rngrd.h> 29 30 #include <star/s3d.h> 31 #include <star/ssf.h> 32 #include <star/ssp.h> 33 #include <star/suvm.h> 34 #include <star/svx.h> 35 36 #include <rsys/double2.h> 37 #include <rsys/double3.h> 38 39 /* Syntactic sugar */ 40 #define DISTANCE_NONE(Dst) ((Dst) >= FLT_MAX) 41 #define SURFACE_EVENT(Event) (!S3D_HIT_NONE(&(Event)->hit)) 42 43 struct event { 44 /* Set to S3D_HIT_NULL if the event occurs in volume.*/ 45 struct s3d_hit hit; 46 47 /* The surface normal is defined only if event is on the surface. It is 48 * normalized and looks towards the incoming direction */ 49 double normal[3]; 50 51 /* Cells in which the event position is located. It makes sense only for an 52 * event in volume */ 53 struct rnatm_cell_pos cells[RNATM_MAX_COMPONENTS_COUNT]; 54 55 double distance; /* Distance from ray origin to scattering position */ 56 }; 57 #define EVENT_NULL__ { \ 58 S3D_HIT_NULL__, {0,0,0}, {RNATM_CELL_POS_NULL__}, DBL_MAX \ 59 } 60 static const struct event EVENT_NULL = EVENT_NULL__; 61 62 /* Arguments of the filtering function used to sample a position */ 63 struct sample_distance_context { 64 struct ssp_rng* rng; 65 struct rnatm* atmosphere; 66 size_t iband; 67 size_t iquad; 68 double wavelength; /* In nm */ 69 enum rnatm_radcoef radcoef; 70 double Ts; /* Sample optical thickness */ 71 72 /* Output data */ 73 struct rnatm_cell_pos* cells; 74 double distance; 75 }; 76 #define SAMPLE_DISTANCE_CONTEXT_NULL__ { \ 77 NULL, NULL, 0, 0, 0, RNATM_RADCOEFS_COUNT__, 0, NULL, DBL_MAX \ 78 } 79 static const struct sample_distance_context SAMPLE_DISTANCE_CONTEXT_NULL = 80 SAMPLE_DISTANCE_CONTEXT_NULL__; 81 82 /******************************************************************************* 83 * Helper functions 84 ******************************************************************************/ 85 static INLINE res_T 86 check_planets_compute_radiance_args 87 (const struct htrdr_planets* cmd, 88 const struct planets_compute_radiance_args* args) 89 { 90 struct rnatm_band_desc band = RNATM_BAND_DESC_NULL; 91 res_T res = RES_OK; 92 93 if(!args || !args->rng) 94 return RES_BAD_ARG; 95 96 /* Invalid thread index */ 97 if(args->ithread >= htrdr_get_threads_count(cmd->htrdr)) 98 return RES_BAD_ARG; 99 100 /* Invalid input direction */ 101 if(!d3_is_normalized(args->path_dir)) 102 return RES_BAD_ARG; 103 104 /* Invalid wavelength */ 105 if(args->wlen < cmd->spectral_domain.wlen_range[0] 106 || args->wlen > cmd->spectral_domain.wlen_range[1]) 107 return RES_BAD_ARG; 108 109 res = rnatm_band_get_desc(cmd->atmosphere, args->iband, &band); 110 if(res != RES_OK) return res; 111 112 /* Inconsistent spectral dimension */ 113 if(args->wlen < band.lower 114 || args->wlen >= band.upper /* Exclusive */ 115 || args->iquad >= band.quad_pts_count) 116 return RES_BAD_ARG; 117 118 return RES_OK; 119 } 120 121 static int 122 sample_position_hit_filter 123 (const struct svx_hit* hit, 124 const double org[3], 125 const double dir[3], 126 const double range[2], 127 void* context) 128 { 129 struct rnatm_get_radcoef_args get_k_args = RNATM_GET_RADCOEF_ARGS_NULL; 130 struct sample_distance_context* ctx = context; 131 double k_min = 0; 132 double k_max = 0; 133 double dst_travelled = 0; 134 int pursue_traversal = 1; 135 ASSERT(hit && org && range && context); 136 (void)range; 137 138 dst_travelled = hit->distance[0]; 139 140 /* Get the k_min, k_max of the voxel */ 141 k_min = rnatm_get_k_svx_voxel 142 (ctx->atmosphere, &hit->voxel, ctx->radcoef, RNATM_SVX_OP_MIN); 143 k_max = rnatm_get_k_svx_voxel 144 (ctx->atmosphere, &hit->voxel, ctx->radcoef, RNATM_SVX_OP_MAX); 145 146 /* Configure the common input arguments to retrieve the radiative coefficient 147 * of a given position */ 148 get_k_args.cells = ctx->cells; 149 get_k_args.iband = ctx->iband; 150 get_k_args.iquad = ctx->iquad; 151 get_k_args.radcoef = ctx->radcoef; 152 get_k_args.k_min = k_min; 153 get_k_args.k_max = k_max; 154 155 for(;;) { 156 /* Compute the optical thickness of the voxel */ 157 const double vox_dst = hit->distance[1] - dst_travelled; 158 const double T = vox_dst * k_max; 159 160 /* A collision occurs behind the voxel */ 161 if(ctx->Ts > T) { 162 ctx->Ts -= T; 163 pursue_traversal = 1; 164 break; 165 166 /* A collision occurs in the voxel */ 167 } else { 168 const double vox_dst_collision = ctx->Ts / k_max; 169 double pos[3]; 170 double k, r; 171 172 /* Calculate the distance travelled to the collision to be queried */ 173 dst_travelled += vox_dst_collision; 174 175 /* Retrieve the radiative coefficient at the collision position */ 176 pos[0] = org[0] + dst_travelled * dir[0]; 177 pos[1] = org[1] + dst_travelled * dir[1]; 178 pos[2] = org[2] + dst_travelled * dir[2]; 179 RNATM(fetch_cell_list(ctx->atmosphere, pos, get_k_args.cells, NULL)); 180 RNATM(get_radcoef(ctx->atmosphere, &get_k_args, &k)); 181 182 r = ssp_rng_canonical(ctx->rng); 183 184 /* Null collision */ 185 if(r > k/k_max) { 186 ctx->Ts = ssp_ran_exp(ctx->rng, 1); 187 188 /* Real collision */ 189 } else { 190 ctx->distance = dst_travelled; 191 pursue_traversal = 0; 192 break; 193 } 194 } 195 } 196 197 return pursue_traversal; 198 } 199 200 static double 201 sample_distance 202 (struct htrdr_planets* cmd, 203 const struct planets_compute_radiance_args* args, 204 struct rnatm_cell_pos* cells, 205 const enum rnatm_radcoef radcoef, 206 const double pos[3], 207 const double dir[3], 208 const double range[2]) 209 { 210 struct rnatm_trace_ray_args rt = RNATM_TRACE_RAY_ARGS_NULL; 211 struct sample_distance_context ctx = SAMPLE_DISTANCE_CONTEXT_NULL; 212 struct svx_hit hit; 213 ASSERT(cmd && args && cells && pos && dir && d3_is_normalized(dir) && range); 214 ASSERT((unsigned)radcoef < RNATM_RADCOEFS_COUNT__); 215 ASSERT(range[0] < range[1]); 216 217 /* Sample an optical thickness */ 218 ctx.Ts = ssp_ran_exp(args->rng, 1); 219 220 /* Setup the remaining arguments of the RT context */ 221 ctx.rng = args->rng; 222 ctx.atmosphere = cmd->atmosphere; 223 ctx.iband = args->iband; 224 ctx.iquad = args->iquad; 225 ctx.wavelength = args->wlen; 226 ctx.radcoef = radcoef; 227 ctx.cells = cells; 228 229 /* Trace the ray into the atmosphere */ 230 d3_set(rt.ray_org, pos); 231 d3_set(rt.ray_dir, dir); 232 rt.ray_range[0] = range[0]; 233 rt.ray_range[1] = range[1]; 234 rt.filter = sample_position_hit_filter; 235 rt.context = &ctx; 236 rt.iband = args->iband; 237 rt.iquad = args->iquad; 238 RNATM(trace_ray(cmd->atmosphere, &rt, &hit)); 239 240 if(SVX_HIT_NONE(&hit)) { /* No collision found */ 241 return INF; 242 } else { /* A (real) collision occured */ 243 return ctx.distance; 244 } 245 } 246 247 static INLINE double 248 transmissivity 249 (struct htrdr_planets* cmd, 250 const struct planets_compute_radiance_args* args, 251 const enum rnatm_radcoef radcoef, 252 const double pos[3], 253 const double dir[3], 254 const double range_max) 255 { 256 struct rnatm_cell_pos cells[RNATM_MAX_COMPONENTS_COUNT]; 257 double range[2]; 258 double dst = 0; 259 ASSERT(range_max >= 0); 260 261 range[0] = 0; 262 range[1] = range_max; 263 dst = sample_distance(cmd, args, cells, radcoef, pos, dir, range); 264 265 if(DISTANCE_NONE(dst)) { 266 return 1.0; /* No collision in the medium */ 267 } else { 268 return 0.0; /* A (real) collision occurs */ 269 } 270 } 271 272 static double 273 direct_contribution 274 (struct htrdr_planets* cmd, 275 const struct planets_compute_radiance_args* args, 276 const double pos[3], 277 const double dir[3], 278 const struct s3d_hit* hit_from) 279 { 280 struct rngrd_trace_ray_args rt = RNGRD_TRACE_RAY_ARGS_DEFAULT; 281 struct s3d_hit hit; 282 double Tr; 283 double Ld; 284 double src_dst; 285 ASSERT(cmd && args && pos && dir); 286 287 /* Is the source hidden? */ 288 d3_set(rt.ray_org, pos); 289 d3_set(rt.ray_dir, dir); 290 if(hit_from) rt.hit_from = *hit_from; 291 RNGRD(trace_ray(cmd->ground, &rt, &hit)); 292 if(!S3D_HIT_NONE(&hit)) return 0; 293 294 /* Calculate the distance between the source and `pos' */ 295 src_dst = htrdr_planets_source_distance_to(cmd->source, pos); 296 ASSERT(src_dst >= 0); 297 298 Tr = transmissivity(cmd, args, RNATM_RADCOEF_Kext, pos, dir, src_dst); 299 Ld = htrdr_planets_source_get_radiance(cmd->source, args->wlen); 300 return Ld * Tr; 301 } 302 303 static void 304 find_event 305 (struct htrdr_planets* cmd, 306 const struct planets_compute_radiance_args* args, 307 const enum rnatm_radcoef radcoef, 308 const double pos[3], 309 const double dir[3], 310 const struct s3d_hit* hit_from, 311 struct event* evt) 312 { 313 struct rngrd_trace_ray_args rt = RNGRD_TRACE_RAY_ARGS_DEFAULT; 314 struct s3d_hit hit; 315 double range[2]; 316 double dst; 317 ASSERT(cmd && args && pos && dir && hit_from && evt); 318 319 *evt = EVENT_NULL; 320 321 /* Look for a surface intersection */ 322 d3_set(rt.ray_org, pos); 323 d3_set(rt.ray_dir, dir); 324 d2(rt.ray_range, 0, INF); 325 rt.hit_from = *hit_from; 326 RNGRD(trace_ray(cmd->ground, &rt, &hit)); 327 328 /* Look for an atmospheric collision */ 329 range[0] = 0; 330 range[1] = hit.distance; 331 dst = sample_distance(cmd, args, evt->cells, radcoef, pos, dir, range); 332 333 /* Event occurs in volume */ 334 if(!DISTANCE_NONE(dst)) { 335 evt->distance = dst; 336 evt->hit = S3D_HIT_NULL; 337 338 /* Event is on surface */ 339 } else if(!S3D_HIT_NONE(&hit)) { 340 /* Normalize the normal and ensure that it points to `dir' */ 341 d3_normalize(evt->normal, d3_set_f3(evt->normal, hit.normal)); 342 if(d3_dot(evt->normal, dir) > 0) d3_minus(evt->normal, evt->normal); 343 344 evt->distance = hit.distance; 345 evt->hit = hit; 346 347 /* No event */ 348 } else { 349 evt->distance = INF; 350 evt->hit = S3D_HIT_NULL; 351 } 352 } 353 354 static INLINE struct ssf_bsdf* 355 create_bsdf 356 (struct htrdr_planets* cmd, 357 const struct planets_compute_radiance_args* args, 358 const struct s3d_hit* hit) 359 { 360 struct rngrd_create_bsdf_args bsdf_args = RNGRD_CREATE_BSDF_ARGS_NULL; 361 struct ssf_bsdf* bsdf = NULL; 362 ASSERT(!S3D_HIT_NONE(hit)); 363 364 /* Retrieve the BSDF at the intersected surface position */ 365 bsdf_args.prim = hit->prim; 366 bsdf_args.barycentric_coords[0] = hit->uv[0]; 367 bsdf_args.barycentric_coords[1] = hit->uv[1]; 368 bsdf_args.barycentric_coords[2] = 1 - hit->uv[0] - hit->uv[1]; 369 bsdf_args.wavelength = args->wlen; 370 bsdf_args.r = ssp_rng_canonical(args->rng); 371 RNGRD(create_bsdf(cmd->ground, &bsdf_args, &bsdf)); 372 373 return bsdf; 374 } 375 376 static INLINE struct ssf_phase* 377 create_phase_fn 378 (struct htrdr_planets* cmd, 379 const struct planets_compute_radiance_args* args, 380 const struct rnatm_cell_pos* cells) /* Cells in which scattering occurs */ 381 { 382 struct rnatm_sample_component_args sample_args = 383 RNATM_SAMPLE_COMPONENT_ARGS_NULL; 384 struct rnatm_cell_create_phase_fn_args phase_fn_args = 385 RNATM_CELL_CREATE_PHASE_FN_ARGS_NULL; 386 struct ssf_phase* phase = NULL; 387 size_t cpnt; 388 ASSERT(cmd && args && cells); 389 390 /* Sample the atmospheric scattering component */ 391 sample_args.cells = cells; 392 sample_args.iband = args->iband; 393 sample_args.iquad = args->iquad; 394 sample_args.radcoef = RNATM_RADCOEF_Ks; 395 sample_args.r = ssp_rng_canonical(args->rng); 396 RNATM(sample_component(cmd->atmosphere, &sample_args, &cpnt)); 397 398 /* Retrieve the component cell in which the scattering position is located */ 399 phase_fn_args.cell = RNATM_GET_COMPONENT_CELL(cells, cpnt); 400 ASSERT(!SUVM_PRIMITIVE_NONE(&phase_fn_args.cell.prim)); 401 402 /* Retrieve the component phase function */ 403 phase_fn_args.wavelength = args->wlen; 404 phase_fn_args.r[0] = ssp_rng_canonical(args->rng); 405 phase_fn_args.r[1] = ssp_rng_canonical(args->rng); 406 RNATM(cell_create_phase_fn(cmd->atmosphere, &phase_fn_args, &phase)); 407 408 return phase; 409 } 410 411 /* In shortwave, return the contribution of the external source at the bounce 412 * position. In longwave, simply return 0 */ 413 static double 414 surface_bounce 415 (struct htrdr_planets* cmd, 416 const struct planets_compute_radiance_args* args, 417 const struct event* sc, 418 const double sc_pos[3], /* Scattering position */ 419 const double in_dir[3], /* Incident direction */ 420 double sc_dir[3], /* Sampled scattering direction */ 421 double *out_refl) /* Surface reflectivity */ 422 { 423 struct ssf_bsdf* bsdf = NULL; 424 const double* N = NULL; 425 double wo[3] = {0,0,0}; 426 double reflectivity = 0; /* Surface reflectivity */ 427 double L = 0; 428 int mask = 0; 429 ASSERT(cmd && args && sc && sc_pos && in_dir && sc_dir && out_refl); 430 ASSERT(d3_dot(sc->normal, in_dir) < 0 && d3_is_normalized(sc->normal)); 431 432 bsdf = create_bsdf(cmd, args, &sc->hit); 433 N = sc->normal; 434 d3_minus(wo, in_dir); /* Match StarSF convention */ 435 ASSERT(d3_dot(wo, N) > 0); 436 437 /* Sample the scattering direction */ 438 reflectivity = ssf_bsdf_sample(bsdf, args->rng, wo, N, sc_dir, &mask, NULL); 439 440 /* Fully absorbs transmissions */ 441 if(mask & SSF_TRANSMISSION) reflectivity = 0; 442 443 /* No external source in longwave */ 444 if(cmd->spectral_domain.type == HTRDR_SPECTRAL_LW) 445 goto exit; 446 447 /* Calculate direct contribution for specular reflection */ 448 if((mask & SSF_SPECULAR) 449 && (mask & SSF_REFLECTION) 450 && htrdr_planets_source_is_targeted(cmd->source, sc_pos, sc_dir)) { 451 const double Ld = direct_contribution(cmd, args, sc_pos, sc_dir, &sc->hit); 452 L = Ld * reflectivity; 453 454 /* Calculate direct contribution in general case */ 455 } else if(!(mask & SSF_SPECULAR)) { 456 double wi[3], pdf; 457 458 /* Sample a direction toward the source */ 459 pdf = htrdr_planets_source_sample_direction(cmd->source, args->rng, sc_pos, wi); 460 461 /* The source is behind the surface */ 462 if(d3_dot(wi, N) <= 0) { 463 L = 0; 464 465 /* The source is above the surface */ 466 } else { 467 const double Ld = direct_contribution(cmd, args, sc_pos, wi, &sc->hit); 468 const double rho = ssf_bsdf_eval(bsdf, wo, N, wi); 469 const double cos_N_wi = d3_dot(N, wi); 470 ASSERT(cos_N_wi > 0); 471 472 L = Ld * rho * cos_N_wi / pdf; 473 } 474 } 475 476 exit: 477 SSF(bsdf_ref_put(bsdf)); 478 *out_refl = reflectivity; 479 return L; 480 } 481 482 /* In shortwave, return the contribution at the scattering position of the 483 * external source. Returns 0 in long wave */ 484 static INLINE double 485 volume_scattering 486 (struct htrdr_planets* cmd, 487 const struct planets_compute_radiance_args* args, 488 const struct event* sc, 489 const double sc_pos[3], /* Scattering position */ 490 const double in_dir[3], /* Incident direction */ 491 double sc_dir[3]) /* Sampled scattering direction */ 492 { 493 struct ssf_phase* phase = NULL; 494 double wo[3] = {0,0,0}; 495 double wi[3] = {0,0,0}; 496 double L = 0; 497 double pdf = 0; 498 double rho = 0; 499 double Ld = 0; 500 ASSERT(cmd && args && sc && sc_pos && in_dir && sc_dir); 501 502 phase = create_phase_fn(cmd, args, sc->cells); 503 d3_minus(wo, in_dir); /* Match StarSF convention */ 504 505 ssf_phase_sample(phase, args->rng, wo, sc_dir, NULL); 506 507 /* In short wave, manage the contribution of the external source */ 508 switch(cmd->spectral_domain.type) { 509 case HTRDR_SPECTRAL_LW: 510 /* Nothing to be done */ 511 break; 512 513 case HTRDR_SPECTRAL_SW: 514 case HTRDR_SPECTRAL_SW_CIE_XYZ: 515 /* Sample a direction toward the source */ 516 pdf = htrdr_planets_source_sample_direction 517 (cmd->source, args->rng, sc_pos, wi); 518 519 /* Calculate the direct contribution at the scattering position */ 520 Ld = direct_contribution(cmd, args, sc_pos, wi, NULL); 521 rho = ssf_phase_eval(phase, wo, wi); 522 L = Ld * rho / pdf; 523 break; 524 525 default: FATAL("Unreachable code.\n"); break; 526 } 527 528 SSF(phase_ref_put(phase)); 529 return L; 530 } 531 532 static INLINE enum rnatm_radcoef 533 sample_volume_event_type 534 (const struct htrdr_planets* cmd, 535 const struct planets_compute_radiance_args* args, 536 struct event* evt) 537 { 538 struct rnatm_get_radcoef_args get_k_args = RNATM_GET_RADCOEF_ARGS_NULL; 539 double ka, kext; 540 double r; 541 ASSERT(cmd && args && evt); 542 543 get_k_args.cells = evt->cells; 544 get_k_args.iband = args->iband; 545 get_k_args.iquad = args->iquad; 546 547 /* Retrieve the absorption coefficient */ 548 get_k_args.radcoef = RNATM_RADCOEF_Ka; 549 RNATM(get_radcoef(cmd->atmosphere, &get_k_args, &ka)); 550 551 /* Retrieve the extinction coefficient */ 552 get_k_args.radcoef = RNATM_RADCOEF_Kext; 553 RNATM(get_radcoef(cmd->atmosphere, &get_k_args, &kext)); 554 555 r = ssp_rng_canonical(args->rng); 556 if(r < ka / kext) { 557 return RNATM_RADCOEF_Ka; /* Absorption */ 558 } else { 559 return RNATM_RADCOEF_Ks; /* Scattering */ 560 } 561 } 562 563 static INLINE double 564 get_temperature 565 (const struct htrdr_planets* cmd, 566 const struct event* evt) 567 { 568 double T = 0; 569 ASSERT(cmd && evt); 570 571 if(!SURFACE_EVENT(evt)) { 572 const struct rnatm_cell_pos* cell = NULL; 573 574 /* Get gas temperature */ 575 cell = &RNATM_GET_COMPONENT_CELL(evt->cells, RNATM_GAS); 576 RNATM(cell_get_gas_temperature(cmd->atmosphere, cell, &T)); 577 578 } else { 579 struct rngrd_get_temperature_args temp_args = RNGRD_GET_TEMPERATURE_ARGS_NULL; 580 581 /* Get ground temperature */ 582 temp_args.prim = evt->hit.prim; 583 temp_args.barycentric_coords[0] = evt->hit.uv[0]; 584 temp_args.barycentric_coords[1] = evt->hit.uv[1]; 585 temp_args.barycentric_coords[2] = 1 - evt->hit.uv[0] - evt->hit.uv[1]; 586 RNGRD(get_temperature(cmd->ground, &temp_args, &T)); 587 588 } 589 return T; 590 } 591 592 /******************************************************************************* 593 * Local functions 594 ******************************************************************************/ 595 double /* Radiance in W/m²/sr/m */ 596 planets_compute_radiance 597 (struct htrdr_planets* cmd, 598 const struct planets_compute_radiance_args* args) 599 { 600 struct s3d_hit hit_from = S3D_HIT_NULL; 601 struct event ev; 602 double pos[3]; 603 double dir[3]; 604 double L = 0; /* Radiance in W/m²/sr/m */ 605 size_t nsc = 0; /* Number of surface or volume scatterings (for debug) */ 606 int longwave = 0; 607 int shortwave = 0; 608 int direct = 0; 609 int diffuse = 0; 610 ASSERT(cmd && check_planets_compute_radiance_args(cmd, args) == RES_OK); 611 612 d3_set(pos, args->path_org); 613 d3_set(dir, args->path_dir); 614 615 longwave = cmd->spectral_domain.type == HTRDR_SPECTRAL_LW; 616 shortwave = !longwave; 617 618 /* In shortwave define which components are enabled */ 619 if(shortwave) { 620 direct = (args->component & PLANETS_RADIANCE_CPNT_DIRECT) != 0; 621 diffuse = (args->component & PLANETS_RADIANCE_CPNT_DIFFUSE) != 0; 622 } 623 624 /* Handle direct shortwave contribution */ 625 if(shortwave 626 && direct 627 && htrdr_planets_source_is_targeted(cmd->source, pos, dir)) { 628 L = direct_contribution(cmd, args, pos, dir, NULL); /* In W/m²/sr/m */ 629 } 630 631 /* Nothing left to do: if only the diffuse component is required 632 * in the SW, the diffuse component should not be computed */ 633 if(shortwave && !diffuse) return L; 634 635 for(;;) { 636 double ev_pos[3]; 637 double sc_dir[3]; 638 639 find_event(cmd, args, RNATM_RADCOEF_Kext, pos, dir, &hit_from, &ev); 640 641 /* No event on surface or in volume. Stop the path */ 642 if(DISTANCE_NONE(ev.distance)) break; 643 644 /* Compute the event position */ 645 ev_pos[0] = pos[0] + dir[0] * ev.distance; 646 ev_pos[1] = pos[1] + dir[1] * ev.distance; 647 ev_pos[2] = pos[2] + dir[2] * ev.distance; 648 649 /* The event occurs on the surface */ 650 if(SURFACE_EVENT(&ev)) { 651 double refl; /* Surface reflectivity */ 652 L += surface_bounce(cmd, args, &ev, ev_pos, dir, sc_dir, &refl); 653 654 /* Check the absorption of the surface with a Russian roulette against 655 * the reflectivity of the surface */ 656 if(ssp_rng_canonical(args->rng) >= refl) break; 657 658 /* Save current intersected surface to avoid self-intersection when 659 * sampling next event */ 660 hit_from = ev.hit; 661 662 /* The event occurs in the volume */ 663 } else { 664 enum rnatm_radcoef ev_type = sample_volume_event_type(cmd, args, &ev); 665 ASSERT(ev_type == RNATM_RADCOEF_Ka || ev_type == RNATM_RADCOEF_Ks); 666 667 /* Absorption. Stop the path */ 668 if(ev_type == RNATM_RADCOEF_Ka) break; 669 670 L += volume_scattering(cmd, args, &ev, ev_pos, dir, sc_dir); 671 hit_from = S3D_HIT_NULL; /* Reset surface intersection */ 672 } 673 674 d3_set(pos, ev_pos); 675 d3_set(dir, sc_dir); 676 677 ++nsc; 678 } 679 680 /* From there, a valid event means that the path has stopped in surface or 681 * volume absorption. In longwave, add the contribution of the internal 682 * source */ 683 if(longwave && !DISTANCE_NONE(ev.distance)) { 684 const double wlen_m = args->wlen * 1.e-9; /* wavelength in meters */ 685 const double temperature = get_temperature(cmd, &ev); /* In K */ 686 L += htrdr_planck_monochromatic(wlen_m, temperature); 687 } 688 689 return L; /* In W/m²/sr/m */ 690 }