htrdr_combustion_compute_radiance_sw.c (35585B)
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 #define _POSIX_C_SOURCE 200112L /* nextafterf */ 25 26 #include "combustion/htrdr_combustion_c.h" 27 #include "combustion/htrdr_combustion_laser.h" 28 #include "combustion/htrdr_combustion_geometry_ray_filter.h" 29 30 #include "core/htrdr.h" 31 #include "core/htrdr_log.h" 32 #include "core/htrdr_geometry.h" 33 #include "core/htrdr_interface.h" 34 #include "core/htrdr_materials.h" 35 36 #include <astoria/atrstm.h> 37 38 #include <star/ssf.h> 39 #include <star/ssp.h> 40 41 #include <rsys/double2.h> 42 #include <rsys/double3.h> 43 #include <rsys/double4.h> 44 #include <rsys/dynamic_array.h> 45 46 #include <math.h> /* nextafterf */ 47 48 enum scattering_type { 49 SCATTERING_IN_VOLUME, 50 SCATTERING_AT_SURFACE, 51 SCATTERING_NONE 52 }; 53 54 /* Define a position along the ray into the semi-transparent medium */ 55 struct position { 56 double distance; /* Distance up to the position */ 57 struct suvm_primitive prim; /* Volumetric primitive of the position */ 58 double bcoords[4]; /* Local coordinate of the position into `prim' */ 59 }; 60 #define POSITION_NULL__ { \ 61 DBL_MAX, /* Distance */ \ 62 SUVM_PRIMITIVE_NULL__, /* Primitive */ \ 63 {0, 0, 0, 0} /* Barycentric coordinates */ \ 64 } 65 static const struct position POSITION_NULL = POSITION_NULL__; 66 67 /* Syntactic sugar to check if the position is valid */ 68 #define POSITION_NONE(Pos) ((Pos)->distance >= FLT_MAX) 69 70 /* Common position but preferentially sampled within a limited range. Its 71 * associated ksi variable defines the correction of the weight due to the 72 * normalization of the sampling pdf, and the recursivity associated with the 73 * null-collision technique. */ 74 struct position_limited { 75 struct position position; 76 double ksi; 77 }; 78 #define POSITION_LIMITED_NULL__ {POSITION_NULL__, 1} 79 static const struct position_limited POSITION_LIMITED_NULL = 80 POSITION_LIMITED_NULL__; 81 82 struct sample_position_context { 83 /* Input data */ 84 struct ssp_rng* rng; /* Random Number Generator */ 85 struct atrstm* medium; /* Semi-transparent medium */ 86 double wavelength; /* Wavelength to handel in nanometer */ 87 enum atrstm_radcoef radcoef; /* Radiative coef used to sample a position */ 88 double Ts; /* Sampled optical thickness */ 89 90 /* Output data */ 91 struct position position; 92 }; 93 #define SAMPLE_POSITION_CONTEXT_NULL__ { \ 94 NULL, /* RNG */ \ 95 NULL, /* Medium */ \ 96 0, /* Wavelength */ \ 97 ATRSTM_RADCOEFS_COUNT__, /* Radiative coefficient */ \ 98 0, /* Optical thickness */ \ 99 POSITION_NULL__ \ 100 } 101 static const struct sample_position_context SAMPLE_POSITION_CONTEXT_NULL = 102 SAMPLE_POSITION_CONTEXT_NULL__; 103 104 struct sample_scattering_limited_context { 105 /* Input data */ 106 struct ssp_rng* rng; /* Random number generator to use */ 107 struct atrstm* medium; /* Semi transparent medium */ 108 double wavelength; /* Wavelength to handle in nanometer */ 109 110 /* Local parameters update during ray traversal */ 111 double ks_2hat; /* Smallest scattering upper-field over the ray range */ 112 double Tmax; /* Scattering optical thickness computed using ks_2hat */ 113 double Ume; /* Normalization of the pdf for the sampled optical thickness */ 114 double sampled_vox_collision_dst; /* Scattering path length */ 115 116 /* Output data */ 117 struct position_limited position_limited; 118 }; 119 #define SAMPLE_SCATTERING_LIMITED_CONTEXT_NULL__ { \ 120 NULL, /* RNG */ \ 121 NULL, /* Medium */ \ 122 -1, /* Wavelength */ \ 123 -1, /* ks_2hat */ \ 124 -1, /* Tau max */ \ 125 -1, /* Ume */ \ 126 -1, /* Sampled collision dst */ \ 127 POSITION_LIMITED_NULL__ \ 128 } 129 static const struct sample_scattering_limited_context 130 SAMPLE_SCATTERING_LIMITED_CONTEXT_NULL = 131 SAMPLE_SCATTERING_LIMITED_CONTEXT_NULL__; 132 133 /******************************************************************************* 134 * Sample a position along a ray into the inhomogeneous medium for a given 135 * radiative coefficient 136 ******************************************************************************/ 137 static int 138 sample_position_hit_filter 139 (const struct svx_hit* hit, 140 const double org[3], 141 const double dir[3], 142 const double range[2], 143 void* context) 144 { 145 atrstm_radcoefs_svx_T radcoefs_svx; 146 struct atrstm_fetch_radcoefs_args fetch_raw_args; 147 struct atrstm_fetch_radcoefs_svx_voxel_args fetch_svx_args; 148 struct sample_position_context* ctx = context; 149 double k_min = 0; 150 double k_max = 0; 151 double traversal_dst = 0; 152 int pursue_traversal = 1; 153 ASSERT(hit && org && dir && range && ctx); 154 (void)range; 155 156 /* Fetch the K<min|max> of the current traversed voxel */ 157 fetch_svx_args.voxel = hit->voxel; 158 fetch_svx_args.radcoefs_mask = BIT(ctx->radcoef); 159 fetch_svx_args.components_mask = ATRSTM_CPNTS_MASK_ALL; 160 fetch_svx_args.operations_mask = ATRSTM_SVX_OP_FLAG_MIN | ATRSTM_SVX_OP_FLAG_MAX; 161 ATRSTM(fetch_radcoefs_svx_voxel(ctx->medium, &fetch_svx_args, radcoefs_svx)); 162 163 k_min = radcoefs_svx[ctx->radcoef][ATRSTM_SVX_OP_MIN]; 164 k_max = radcoefs_svx[ctx->radcoef][ATRSTM_SVX_OP_MAX]; 165 166 /* Setup the constants of the 'fetch' function for the current voxel */ 167 fetch_raw_args.wavelength = ctx->wavelength; 168 fetch_raw_args.radcoefs_mask = BIT(ctx->radcoef); 169 fetch_raw_args.components_mask = ATRSTM_CPNTS_MASK_ALL; 170 fetch_raw_args.k_min[ctx->radcoef] = k_min; 171 fetch_raw_args.k_max[ctx->radcoef] = k_max; 172 173 /* Initialised the already traversed distance to the distance from which the 174 * current ray enters into the current voxel */ 175 traversal_dst = hit->distance[0]; 176 177 for(;;) { 178 /* Compute the optical thickness of the current leaf */ 179 const double vox_dst = hit->distance[1] - traversal_dst; 180 const double T = vox_dst * k_max; 181 182 /* A collision occurs behind the voxel */ 183 if(ctx->Ts > T) { 184 ctx->Ts -= T; 185 pursue_traversal = 1; 186 break; 187 188 /* A collision occurs _in_ the voxel */ 189 } else { 190 const double vox_collision_dst = ctx->Ts / k_max; 191 atrstm_radcoefs_T radcoefs; 192 double x[3]; 193 double k; 194 double r; 195 196 /* Compute the traversed distance up to the challenged collision */ 197 traversal_dst += vox_collision_dst; 198 199 /* Compute the world space position where a collision may occur */ 200 x[0] = org[0] + traversal_dst * dir[0]; 201 x[1] = org[1] + traversal_dst * dir[1]; 202 x[2] = org[2] + traversal_dst * dir[2]; 203 204 /* Fetch the radiative coefficient at the collision position */ 205 ATRSTM(at(ctx->medium, x, &fetch_raw_args.prim, fetch_raw_args.bcoords)); 206 if(SUVM_PRIMITIVE_NONE(&fetch_raw_args.prim)) { 207 k = 0; 208 } else { 209 ATRSTM(fetch_radcoefs(ctx->medium, &fetch_raw_args, radcoefs)); 210 k = radcoefs[ctx->radcoef]; 211 } 212 213 r = ssp_rng_canonical(ctx->rng); 214 215 if(r > k/k_max) { /* Null collision */ 216 ctx->Ts = ssp_ran_exp(ctx->rng, 1); 217 } else { /* Real collision */ 218 /* Setup output data */ 219 ctx->position.distance = traversal_dst; 220 ctx->position.prim = fetch_raw_args.prim; 221 d4_set(ctx->position.bcoords, fetch_raw_args.bcoords); 222 223 /* Stop ray traversal */ 224 pursue_traversal = 0; 225 break; 226 } 227 } 228 } 229 return pursue_traversal; 230 } 231 232 static void 233 sample_position 234 (struct htrdr_combustion* cmd, 235 struct ssp_rng* rng, 236 const enum atrstm_radcoef radcoef, 237 const double pos[3], 238 const double dir[3], 239 const double range[2], 240 struct position* position) 241 { 242 struct atrstm_trace_ray_args args = ATRSTM_TRACE_RAY_ARGS_DEFAULT; 243 struct sample_position_context sample_pos_ctx = SAMPLE_POSITION_CONTEXT_NULL; 244 struct svx_hit svx_hit = SVX_HIT_NULL; 245 double wlen = 0; 246 ASSERT(cmd && rng && pos && dir && range); 247 248 wlen = htrdr_combustion_laser_get_wavelength(cmd->laser); 249 250 /* Sample an optical thickness */ 251 sample_pos_ctx.Ts = ssp_ran_exp(rng, 1); 252 253 /* Setup the arguments of the function invoked per voxel (i.e. the filter 254 * function) */ 255 sample_pos_ctx.rng = rng; 256 sample_pos_ctx.medium = cmd->medium; 257 sample_pos_ctx.wavelength = wlen; 258 sample_pos_ctx.radcoef = radcoef; 259 260 /* Setup ray tracing arguments */ 261 d3_set(args.ray_org, pos); 262 d3_set(args.ray_dir, dir); 263 d2_set(args.ray_range, range); 264 args.filter = sample_position_hit_filter; 265 args.context = &sample_pos_ctx; 266 267 /* Trace the ray into the heterogeneous medium */ 268 ATRSTM(trace_ray(cmd->medium, &args, &svx_hit)); 269 270 if(SVX_HIT_NONE(&svx_hit)) { 271 /* No collision with the medium was found */ 272 *position = POSITION_NULL; 273 } else { 274 /* A collision occurs into the medium */ 275 *position = sample_pos_ctx.position; 276 } 277 } 278 279 /******************************************************************************* 280 * Preferentially sample a scattering position into an inhomogeneous medium 281 * according to a limited range 282 ******************************************************************************/ 283 /* Find the constant Ks max (named ks_2hat) along the traced ray */ 284 static int 285 sample_scattering_limited_find_ks_2hat 286 (const struct svx_hit* hit, 287 const double org[3], 288 const double dir[3], 289 const double range[2], 290 void* context) 291 { 292 struct sample_scattering_limited_context* ctx = context; 293 struct atrstm_fetch_radcoefs_svx_voxel_args fetch_svx_args; 294 atrstm_radcoefs_svx_T radcoefs_svx; 295 ASSERT(hit && org && dir && range && context); 296 (void)org, (void)dir; 297 298 /* In all situations, initialise ks_2hat with the ks_max of the root node */ 299 if(hit->voxel.depth == 0) { 300 fetch_svx_args.voxel = hit->voxel; 301 fetch_svx_args.radcoefs_mask = ATRSTM_RADCOEF_FLAG_Ks; 302 fetch_svx_args.components_mask = ATRSTM_CPNTS_MASK_ALL; 303 fetch_svx_args.operations_mask = ATRSTM_SVX_OP_FLAG_MAX; 304 ATRSTM(fetch_radcoefs_svx_voxel(ctx->medium, &fetch_svx_args, radcoefs_svx)); 305 ctx->ks_2hat = radcoefs_svx[ATRSTM_RADCOEF_Ks][ATRSTM_SVX_OP_MAX]; 306 307 /* Update ks_2hat with the ks_max of the current descending node until the ray 308 * range was no more included into this node */ 309 } else if(hit->distance[0] <= range[0] && range[1] <= hit->distance[1]) { 310 fetch_svx_args.voxel = hit->voxel; 311 fetch_svx_args.radcoefs_mask = ATRSTM_RADCOEF_FLAG_Ks; 312 fetch_svx_args.components_mask = ATRSTM_CPNTS_MASK_ALL; 313 fetch_svx_args.operations_mask = ATRSTM_SVX_OP_FLAG_MAX; 314 ATRSTM(fetch_radcoefs_svx_voxel(ctx->medium, &fetch_svx_args, radcoefs_svx)); 315 ctx->ks_2hat = radcoefs_svx[ATRSTM_RADCOEF_Ks][ATRSTM_SVX_OP_MAX]; 316 } 317 318 /* Do not stop here: keep diving up to the leaves */ 319 return 0; 320 } 321 322 static int 323 sample_scattering_limited_hit_filter 324 (const struct svx_hit* hit, 325 const double org[3], 326 const double dir[3], 327 const double range[2], 328 void* context) 329 { 330 atrstm_radcoefs_svx_T radcoefs_svx; 331 struct atrstm_fetch_radcoefs_args fetch_raw_args; 332 struct atrstm_fetch_radcoefs_svx_voxel_args fetch_svx_args; 333 struct sample_scattering_limited_context* ctx = context; 334 double ks_min = 0; 335 double ks_max = 0; 336 double traversal_dst = 0; 337 int pursue_traversal = 1; 338 ASSERT(hit && org && dir && range && ctx); 339 (void)range; 340 341 /* Fetch the Ks<min|max> of the current traversed voxel */ 342 fetch_svx_args.voxel = hit->voxel; 343 fetch_svx_args.radcoefs_mask = ATRSTM_RADCOEF_FLAG_Ks; 344 fetch_svx_args.components_mask = ATRSTM_CPNTS_MASK_ALL; 345 fetch_svx_args.operations_mask = ATRSTM_SVX_OP_FLAG_MIN | ATRSTM_SVX_OP_FLAG_MAX; 346 ATRSTM(fetch_radcoefs_svx_voxel(ctx->medium, &fetch_svx_args, radcoefs_svx)); 347 348 ks_min = radcoefs_svx[ATRSTM_RADCOEF_Ks][ATRSTM_SVX_OP_MIN]; 349 ks_max = radcoefs_svx[ATRSTM_RADCOEF_Ks][ATRSTM_SVX_OP_MAX]; 350 ASSERT(ks_max <= ctx->ks_2hat); 351 352 /* Setup the constants of the 'fetch' function for the current voxel */ 353 fetch_raw_args.wavelength = ctx->wavelength; 354 fetch_raw_args.radcoefs_mask = ATRSTM_RADCOEF_FLAG_Ks; 355 fetch_raw_args.components_mask = ATRSTM_CPNTS_MASK_ALL; 356 fetch_raw_args.k_min[ATRSTM_RADCOEF_Ks] = ks_min; 357 fetch_raw_args.k_max[ATRSTM_RADCOEF_Ks] = ks_max; 358 359 /* Initialised the already traversed distance to the distance from which the 360 * current ray enters into the current voxel */ 361 traversal_dst = hit->distance[0]; 362 363 /* Tmax was not already computed */ 364 if(ctx->Tmax < 0) { 365 ctx->Tmax = (range[1] - range[0]) * ctx->ks_2hat; 366 ctx->Ume = 1 - exp(-ctx->Tmax); 367 } 368 369 /* No scattering into the whole traversed laser sheet */ 370 if(ctx->Tmax == 0) { 371 pursue_traversal = 1; 372 return pursue_traversal; 373 } 374 ASSERT(ctx->Tmax > 0); 375 376 for(;;) { 377 atrstm_radcoefs_T radcoefs; 378 double vox_dst; 379 double tau; 380 double ks; 381 double r; 382 383 /* Compute the remaining distance to traverse in the current voxel */ 384 vox_dst = hit->distance[1] - traversal_dst; 385 386 /* A collision distance was not already sampled */ 387 if(ctx->sampled_vox_collision_dst < 0) { 388 tau = ssp_ran_exp_truncated(ctx->rng, 1, ctx->Tmax); 389 ctx->sampled_vox_collision_dst = tau / ctx->ks_2hat; 390 391 /* Update the ksi output data */ 392 ctx->position_limited.ksi *= ctx->Ume; 393 } 394 395 /* Compute the traversed distance up to the challenged collision */ 396 traversal_dst = traversal_dst + ctx->sampled_vox_collision_dst; 397 398 /* The collision to challenge lies behind the current voxel */ 399 if(traversal_dst > hit->distance[1]) { 400 ctx->sampled_vox_collision_dst -= vox_dst; 401 pursue_traversal = 1; 402 break; 403 } 404 ASSERT(traversal_dst >= hit->distance[0]); 405 406 r = ssp_rng_canonical(ctx->rng); 407 if(r >= ks_max / ctx->ks_2hat) { /* Null collision */ 408 ctx->sampled_vox_collision_dst = -1; 409 } else { /* Collision with ks_max */ 410 double x[3]; 411 412 /* Compute the world space position where a collision may occur */ 413 x[0] = org[0] + traversal_dst * dir[0]; 414 x[1] = org[1] + traversal_dst * dir[1]; 415 x[2] = org[2] + traversal_dst * dir[2]; 416 417 /* Fetch the radiative coefficient at the collision position */ 418 ATRSTM(at(ctx->medium, x, &fetch_raw_args.prim, fetch_raw_args.bcoords)); 419 if(SUVM_PRIMITIVE_NONE(&fetch_raw_args.prim)) { 420 ks = 0; 421 } else { 422 ATRSTM(fetch_radcoefs(ctx->medium, &fetch_raw_args, radcoefs)); 423 ks = radcoefs[ATRSTM_RADCOEF_Ks]; 424 } 425 426 r = ssp_rng_canonical(ctx->rng); 427 428 if(r >= ks/ks_max) { /* Null collision */ 429 ctx->sampled_vox_collision_dst = -1; 430 } else { /* Real collision */ 431 /* Setup output data */ 432 ctx->position_limited.position.distance = traversal_dst; 433 ctx->position_limited.position.prim = fetch_raw_args.prim; 434 d4_set(ctx->position_limited.position.bcoords, fetch_raw_args.bcoords); 435 436 /* Stop ray traversal */ 437 pursue_traversal = 0; 438 break; 439 } 440 } 441 } 442 return pursue_traversal; 443 } 444 445 static void 446 sample_scattering_limited 447 (struct htrdr_combustion* cmd, 448 struct ssp_rng* rng, 449 const double pos[3], 450 const double dir[3], 451 const double range[2], 452 struct position_limited* position) 453 { 454 struct atrstm_trace_ray_args args = ATRSTM_TRACE_RAY_ARGS_DEFAULT; 455 struct sample_scattering_limited_context sample_scattering_limited_ctx = 456 SAMPLE_SCATTERING_LIMITED_CONTEXT_NULL; 457 struct svx_hit svx_hit = SVX_HIT_NULL; 458 double wlen = 0; 459 ASSERT(cmd && rng && pos && dir && position); 460 461 wlen = htrdr_combustion_laser_get_wavelength(cmd->laser); 462 463 /* Setup the trace ray context */ 464 sample_scattering_limited_ctx.rng = rng; 465 sample_scattering_limited_ctx.medium = cmd->medium; 466 sample_scattering_limited_ctx.wavelength = wlen; 467 sample_scattering_limited_ctx.ks_2hat = 0; 468 sample_scattering_limited_ctx.Tmax = -1; 469 sample_scattering_limited_ctx.Ume = -1; 470 sample_scattering_limited_ctx.sampled_vox_collision_dst = -1; 471 sample_scattering_limited_ctx.position_limited = POSITION_LIMITED_NULL; 472 473 /* Setup the input arguments for the ray tracing into the medium */ 474 d3_set(args.ray_org, pos); 475 d3_set(args.ray_dir, dir); 476 d2_set(args.ray_range, range); 477 args.challenge = sample_scattering_limited_find_ks_2hat; 478 args.filter = sample_scattering_limited_hit_filter; 479 args.context = &sample_scattering_limited_ctx; 480 481 /* Trace the ray into the medium to compute the transmissivity */ 482 ATRSTM(trace_ray(cmd->medium, &args, &svx_hit)); 483 484 if(SVX_HIT_NONE(&svx_hit)) { 485 /* No scattering event was found */ 486 *position = POSITION_LIMITED_NULL; 487 } else { 488 /* A scattering event occurs into the medium */ 489 *position = sample_scattering_limited_ctx.position_limited; 490 } 491 } 492 493 /******************************************************************************* 494 * Helper functions 495 ******************************************************************************/ 496 static double 497 transmissivity 498 (struct htrdr_combustion* cmd, 499 struct ssp_rng* rng, 500 const enum atrstm_radcoef radcoef, 501 const double pos[3], 502 const double dir[3], 503 const double range[2]) 504 { 505 struct position position = POSITION_NULL; 506 507 /* Degenerated range => no attenuation along dir */ 508 if(range[1] <= range[0]) return 1; 509 510 sample_position(cmd, rng, radcoef, pos, dir, range, &position); 511 512 if(POSITION_NONE(&position)) { 513 return 1; /* No collision with the medium */ 514 } else { 515 return 0; /* A collision occurs */ 516 } 517 } 518 519 /* This function computes the contribution of the in-laser scattering for the 520 * current scattering position of the reverse optical path 'pos' and current 521 * direction 'dir'. One contribution has to be computed for each scattering 522 * position 'xsc' in order to compute the value of the Monte-Carlo weight for 523 * the optical path (weight of the statistical realization). 524 * __ 525 * /\ dir 526 * xout / 527 * +---------------------x--------------- - - - 528 * lse | ---> wi / 529 * (laser surface | ---> <----x xsc 530 * emission) | ---> / 531 * +-----------------x------------------- - - - 532 * / xin 533 * / 534 * o pos */ 535 static double 536 laser_once_scattered 537 (struct htrdr_combustion* cmd, 538 const size_t ithread, 539 struct ssp_rng* rng, 540 const double pos[3], 541 const double dir[3], 542 const double range_in[2]) 543 { 544 /* Phase function */ 545 struct ssf_phase* phase = NULL; 546 547 /* Surface ray tracing */ 548 struct htrdr_geometry_trace_ray_args rt_args = 549 HTRDR_GEOMETRY_TRACE_RAY_ARGS_NULL; 550 struct geometry_ray_filter_context rt_ctx = GEOMETRY_RAY_FILTER_CONTEXT_NULL; 551 struct s3d_hit hit = S3D_HIT_NULL; 552 553 /* Scattering position into the laser sheet */ 554 struct position_limited sc_sample = POSITION_LIMITED_NULL; 555 double xsc[3] = {0, 0, 0}; /* World space position */ 556 557 /* The transmissivities to compute */ 558 double Tr_ext_pos_xin = 0; 559 double Tr_ext_xsc_lse = 0; 560 double Tr_abs_xin_xsc = 0; 561 562 /* Laser data */ 563 double laser_hit_dst[2] = {0, 0}; 564 double laser_flux_density = 0; /* In W/m^2 */ 565 double wlen = 0; /* In nm */ 566 567 /* Miscellaneous variables */ 568 double wi[3] = {0, 0, 0}; 569 double wo[3] = {0, 0, 0}; 570 double range[2] = {0, DBL_MAX}; 571 double R = 0; 572 573 /* Radiance to compute in W/m^2/sr */ 574 double L = 0; 575 576 ASSERT(cmd && rng && pos && dir && range_in); 577 ASSERT(range_in[0] <= range_in[1]); 578 579 /* Fetch the laser wavelength */ 580 wlen = htrdr_combustion_laser_get_wavelength(cmd->laser); 581 582 /* Find the intersections along dir with the laser sheet */ 583 range[0] = range_in[0]; 584 range[1] = range_in[1]; 585 htrdr_combustion_laser_trace_ray(cmd->laser, pos, dir, range, laser_hit_dst); 586 587 /* No intersection with the laser sheet => no laser contribution */ 588 if(HTRDR_COMBUSTION_LASER_HIT_NONE(laser_hit_dst)) return 0; 589 ASSERT(laser_hit_dst[0] <= laser_hit_dst[1]); 590 ASSERT(laser_hit_dst[1] <= range_in[1]); 591 592 /* Compute the transmissivity from 'pos' to 'xin' */ 593 range[0] = 0; 594 range[1] = laser_hit_dst[0]; 595 Tr_ext_pos_xin = transmissivity(cmd, rng, ATRSTM_RADCOEF_Kext, pos, dir, range); 596 if(Tr_ext_pos_xin == 0) return 0; /* no laser contribution */ 597 598 /* Sample the scattering position into the laser sheet */ 599 range[0] = laser_hit_dst[0]; 600 range[1] = laser_hit_dst[1]; 601 sample_scattering_limited(cmd, rng, pos, dir, range, &sc_sample); 602 if(POSITION_NONE(&sc_sample.position)) return 0; /* No laser contribution */ 603 ASSERT(laser_hit_dst[0] <= sc_sample.position.distance); 604 ASSERT(laser_hit_dst[1] >= sc_sample.position.distance); 605 606 /* Compute the transmissivity from 'xin' to 'xsc' */ 607 range[0] = laser_hit_dst[0]; /* <=> xin */ 608 range[1] = sc_sample.position.distance; /* <=> xsc */ 609 Tr_abs_xin_xsc = transmissivity(cmd, rng, ATRSTM_RADCOEF_Ka, pos, dir, range); 610 if(Tr_abs_xin_xsc == 0) return 0; /* No laser contribution */ 611 612 /* Compute the scattering position */ 613 xsc[0] = pos[0] + dir[0] * sc_sample.position.distance; 614 xsc[1] = pos[1] + dir[1] * sc_sample.position.distance; 615 xsc[2] = pos[2] + dir[2] * sc_sample.position.distance; 616 617 /* Retrieve the direction toward the laser surface emission */ 618 htrdr_combustion_laser_get_direction(cmd->laser, wi); 619 d3_minus(wi, wi); 620 621 /* Find the intersection with the combustion chamber */ 622 if(cmd->geom) { /* Is there a combustion chamber? */ 623 /* Setup the ray to trace */ 624 d3_set(rt_args.ray_org, xsc); 625 d3_set(rt_args.ray_dir, wi); 626 rt_args.ray_range[0] = 0; 627 rt_args.ray_range[1] = INF; 628 rt_args.hit_from = S3D_HIT_NULL; 629 630 /* Configure the "X-ray" surface ray trace filtering. This filtering 631 * function helps to avoid intersections with the side of the surfaces 632 * facing inside the combustion chamber to allow the rays to exit out. */ 633 rt_args.filter = geometry_ray_filter_discard_medium_interface; 634 rt_args.filter_context = &rt_ctx; 635 rt_ctx.geom = cmd->geom; 636 rt_ctx.medium_name = "chamber"; 637 638 HTRDR(geometry_trace_ray(cmd->geom, &rt_args, &hit)); 639 640 /* If a surface was intersected then the laser surface emissions is 641 * occluded along `wi' and thus there is no laser contribution */ 642 if(!S3D_HIT_NONE(&hit)) return 0; 643 } 644 645 /* Compute the transmissivity from xsc to the laser surface emission */ 646 range[0] = 0; 647 range[1] = htrdr_combustion_laser_compute_surface_plane_distance 648 (cmd->laser, xsc); 649 ASSERT(range[1] >= 0); 650 Tr_ext_xsc_lse = transmissivity(cmd, rng, ATRSTM_RADCOEF_Kext, xsc, wi, range); 651 if(Tr_ext_xsc_lse == 0) return 0; /* No laser contribution */ 652 653 /* Retrieve phase function */ 654 phase = combustion_fetch_phase_function 655 (cmd, wlen, &sc_sample.position.prim, sc_sample.position.bcoords, ithread); 656 657 /* Evaluate the phase function at the scattering position */ 658 d3_minus(wo, dir); /* Ensure SSF convention */ 659 R = ssf_phase_eval(phase, wo, wi); /* In sr^-1 */ 660 661 laser_flux_density = htrdr_combustion_laser_get_flux_density(cmd->laser); 662 663 L = Tr_ext_pos_xin 664 * Tr_abs_xin_xsc 665 * sc_sample.ksi * R 666 * Tr_ext_xsc_lse 667 * laser_flux_density; 668 return L; 669 } 670 671 static INLINE void 672 sample_scattering_direction 673 (struct htrdr_combustion* cmd, 674 const size_t ithread, 675 struct ssp_rng* rng, 676 const struct position* scattering, 677 const double wlen, 678 const double incoming_dir[3], 679 double scattering_dir[3]) 680 { 681 struct ssf_phase* phase = NULL; 682 double wo[3]; 683 684 ASSERT(cmd && rng && scattering && incoming_dir && scattering_dir); 685 ASSERT(!POSITION_NONE(scattering)); 686 687 /* Fetch the phase function */ 688 phase = combustion_fetch_phase_function 689 (cmd, wlen, &scattering->prim, scattering->bcoords, ithread); 690 691 /* Sample a new optical path direction from the phase function */ 692 d3_minus(wo, incoming_dir); /* Ensure SSF convention */ 693 ssf_phase_sample(phase, rng, wo, scattering_dir, NULL); 694 } 695 696 /* Return the bounce reflectivity */ 697 static double 698 sample_bounce_direction 699 (struct htrdr_combustion* cmd, 700 const size_t ithread, 701 struct ssp_rng* rng, 702 const struct s3d_hit* hit, 703 const double wlen, 704 const double incoming_dir[3], 705 double bounce_dir[3]) 706 { 707 struct htrdr_interface interf = HTRDR_INTERFACE_NULL; 708 const struct htrdr_mtl* mtl = NULL; 709 struct ssf_bsdf* bsdf = NULL; 710 double wo[3]; 711 double N[3]; 712 double bounce_reflectivity; 713 int bsdf_type = 0; 714 715 ASSERT(cmd && rng && hit && incoming_dir && bounce_dir); 716 ASSERT(!S3D_HIT_NONE(hit)); 717 718 /* Recover the bsdf of the intersected interface */ 719 htrdr_geometry_get_interface(cmd->geom, hit, &interf); 720 mtl = htrdr_interface_fetch_hit_mtl(&interf, incoming_dir, hit); 721 HTRDR(mtl_create_bsdf(cmd->htrdr, mtl, ithread, wlen, rng, &bsdf)); 722 723 /* Surface normal */ 724 d3_set_f3(N, hit->normal); 725 d3_normalize(N, N); 726 if(d3_dot(N, incoming_dir) > 0) d3_minus(N, N); /* Ensure SSF convention */ 727 728 /* Sample a new optical path direction from the brdf function */ 729 d3_minus(wo, incoming_dir); /* Ensure SSF convention */ 730 bounce_reflectivity = ssf_bsdf_sample 731 (bsdf, rng, wo, N, bounce_dir, &bsdf_type, NULL); 732 733 if(!(bsdf_type & SSF_REFLECTION)) { /* Handle only reflections */ 734 bounce_reflectivity = 0; 735 } 736 737 SSF(bsdf_ref_put(bsdf)); 738 739 return bounce_reflectivity; 740 } 741 742 static res_T 743 move_to_scattering_position 744 (struct htrdr_combustion* cmd, 745 const double pos[3], 746 const double dir[3], 747 const double sc_distance, 748 const enum scattering_type sc_type, 749 double out_pos[3]) 750 { 751 res_T res = RES_OK; 752 ASSERT(cmd && pos && dir && sc_distance >= 0 && out_pos); 753 ASSERT(sc_type == SCATTERING_IN_VOLUME || sc_type == SCATTERING_AT_SURFACE); 754 755 /* The scattering position is at a surface or in an open air volume */ 756 if(sc_type == SCATTERING_AT_SURFACE || cmd->geom == NULL) { 757 out_pos[0] = pos[0] + dir[0] * sc_distance; 758 out_pos[1] = pos[1] + dir[1] * sc_distance; 759 out_pos[2] = pos[2] + dir[2] * sc_distance; 760 761 /* The scattering position is in a volume surrounded by the geometry of the 762 * combustion chamber. Be careful when moving along 'dir'; due to numerical 763 * uncertainty, the diffusion position could be outside the combustion 764 * chamber */ 765 } else { 766 const int MAX_ATTEMPTS = 10; /* Max #attempts before reporting an error */ 767 int iattempt = 0; /* Index of the current attempt */ 768 769 double tmp[3]; /* Temporary vector */ 770 double pos_to_challenge[3]; /* Scattering position to challenge */ 771 double dst; /* Distance up to the scattering position */ 772 773 /* Search distance of a near position onto the geometry of the combustion 774 * chamber */ 775 double search_radius; 776 777 search_radius = htrdr_geometry_get_epsilon(cmd->geom) * 10.0; 778 dst = sc_distance; 779 780 while(iattempt < MAX_ATTEMPTS) { 781 struct htrdr_interface interf = HTRDR_INTERFACE_NULL; 782 struct s3d_hit hit = S3D_HIT_NULL; 783 double hit_pos[3]; 784 double N[3]; 785 double sign; 786 const struct htrdr_mtl* mtl = NULL; 787 788 /* Move to the scattering position */ 789 pos_to_challenge[0] = pos[0] + dir[0] * dst; 790 pos_to_challenge[1] = pos[1] + dir[1] * dst; 791 pos_to_challenge[2] = pos[2] + dir[2] * dst; 792 793 /* Find the geometry near the scattering position */ 794 HTRDR(geometry_find_closest_point 795 (cmd->geom, pos_to_challenge, search_radius, &hit)); 796 797 /* No geometry => the scattering position to challenge is necessarily 798 * inside the combustion chamber. */ 799 if(S3D_HIT_NONE(&hit)) break; 800 801 /* Retrieve the property of the geometry near the scattering position */ 802 htrdr_geometry_get_hit_position(cmd->geom, &hit, hit_pos); 803 htrdr_geometry_get_interface(cmd->geom, &hit, &interf); 804 805 /* Fetch the material looking toward the scattering position */ 806 d3_normalize(N, d3_set_f3(N, hit.normal)); 807 sign = d3_dot(N, d3_sub(tmp, pos_to_challenge, hit_pos)); 808 mtl = sign < 0 ? &interf.mtl_back : &interf.mtl_front; 809 810 /* The scattering position is inside the combustion chamber. Everything 811 * is fine. */ 812 if(!strcmp(mtl->name, "chamber")) break; 813 814 /* The scattering position is outside the combustion chamber! Due to 815 * numerical uncertainty, while moving along 'dir' we accidentally 816 * crossed the combustion chamber geometry. To deal with this problem, we 817 * try to move slightly less than expected. By "slightly less" we mean 818 * the next representable single precision floating point value, directly 819 * less than the previous challenge distance. */ 820 dst = (double)nextafterf((float)dst, 0); 821 822 /* Next trial */ 823 iattempt += 1; 824 } 825 826 /* Max attempts is reached. Report an error */ 827 if(iattempt == MAX_ATTEMPTS) { 828 htrdr_log_warn(cmd->htrdr, 829 "The scattering position {%g, %g, %g} " 830 "is outside the combustion chamber.\n", 831 pos_to_challenge[0], 832 pos_to_challenge[1], 833 pos_to_challenge[2]); 834 res = RES_BAD_OP; 835 goto error; 836 } 837 838 /* A valid scattering position was found */ 839 out_pos[0] = pos_to_challenge[0]; 840 out_pos[1] = pos_to_challenge[1]; 841 out_pos[2] = pos_to_challenge[2]; 842 } 843 844 exit: 845 return res; 846 error: 847 goto exit; 848 } 849 850 /******************************************************************************* 851 * Local functions 852 ******************************************************************************/ 853 extern LOCAL_SYM res_T 854 combustion_compute_radiance_sw 855 (struct htrdr_combustion* cmd, 856 const size_t ithread, 857 struct ssp_rng* rng, 858 const double pos_in[3], 859 const double dir_in[3], 860 double* out_weight) 861 { 862 /* Surface ray tracing */ 863 struct htrdr_geometry_trace_ray_args rt_args = 864 HTRDR_GEOMETRY_TRACE_RAY_ARGS_NULL; 865 struct geometry_ray_filter_context rt_ctx = GEOMETRY_RAY_FILTER_CONTEXT_NULL; 866 struct s3d_hit hit_curr = S3D_HIT_NULL; /* Current surface hit */ 867 struct s3d_hit hit_prev = S3D_HIT_NULL; /* Previous surface hit */ 868 869 /* Transmissivity between the probe position (i.e. 'pos_in') and the current 870 * scattering position over the reverse scattering path */ 871 double Tr_abs = 1; 872 873 /* Monte carlo weight of the simulated optical path */ 874 double weight = 0; 875 876 /* Miscellaneous variables */ 877 double pos[3] = {0, 0, 0}; 878 double dir[3] = {0, 0, 0}; 879 double range[2] = {0, DBL_MAX}; 880 double wlen = 0; 881 enum scattering_type sc_type = SCATTERING_NONE; 882 res_T res = RES_OK; 883 ASSERT(cmd && rng && pos_in && dir_in && out_weight); 884 885 d3_set(pos, pos_in); 886 d3_set(dir, dir_in); 887 888 wlen = htrdr_combustion_laser_get_wavelength(cmd->laser); 889 890 Tr_abs = 1; 891 weight = 0; 892 893 /* Configure the "X-ray" surface ray trace filtering. This filtering function 894 * helps to avoid intersections with the side of the surfaces facing outside 895 * the combustion chamber to allow primary rays to enter. */ 896 rt_args.filter = geometry_ray_filter_discard_medium_interface; 897 rt_args.filter_context = &rt_ctx; 898 rt_ctx.geom = cmd->geom; 899 rt_ctx.medium_name = "air"; 900 901 for(;;) { 902 struct position scattering = POSITION_NULL; 903 double laser_sheet_emissivity = 0; /* In W/m^2/sr */ 904 double Tr_abs_pos_xsc = 0; 905 double wi[3]; 906 double sc_distance; 907 908 if(cmd->geom) { /* Is there a combustion chamber? */ 909 /* Find the intersection with the combustion chamber geometry */ 910 d3_set(rt_args.ray_org, pos); 911 d3_set(rt_args.ray_dir, dir); 912 d2(rt_args.ray_range, 0, DBL_MAX); 913 rt_args.hit_from = hit_prev; /* Avoid self intersection */ 914 HTRDR(geometry_trace_ray(cmd->geom, &rt_args, &hit_curr)); 915 916 /* Deactivate the "X-ray" filtering function. It is only used during the 917 * first step of the random walk to allow primary rays (i.e. camera rays) 918 * to enter the combustion chamber. */ 919 rt_args.filter = NULL; 920 } 921 922 /* Handle the laser contribution */ 923 range[0] = 0; 924 range[1] = hit_curr.distance; 925 926 laser_sheet_emissivity = laser_once_scattered(cmd, ithread, rng, pos, dir, range); 927 weight += Tr_abs * laser_sheet_emissivity; 928 929 /* Sample a scattering position */ 930 range[0] = 0; 931 range[1] = hit_curr.distance; 932 sample_position(cmd, rng, ATRSTM_RADCOEF_Ks, pos, dir, range, &scattering); 933 934 if(!POSITION_NONE(&scattering)) { /* Scattering event in volume */ 935 sc_type = SCATTERING_IN_VOLUME; 936 sc_distance = scattering.distance; 937 } else if(!S3D_HIT_NONE(&hit_curr)) { /* Scattering event at surface */ 938 sc_type = SCATTERING_AT_SURFACE; 939 sc_distance = hit_curr.distance; 940 } else { /* No scattering event. Stop the random walk */ 941 sc_type = SCATTERING_NONE; 942 break; 943 } 944 945 /* Compute the absorption transmissivity */ 946 range[0] = 0; 947 range[1] = sc_distance; 948 Tr_abs_pos_xsc = transmissivity(cmd, rng, ATRSTM_RADCOEF_Ka, pos, dir, range); 949 if(Tr_abs_pos_xsc == 0) break; 950 951 /* Update the overall absorption transmissivity of the optical path */ 952 Tr_abs *= Tr_abs_pos_xsc; 953 954 /* Update the position of the optical path */ 955 res = move_to_scattering_position(cmd, pos, dir, sc_distance, sc_type, pos); 956 if(res != RES_OK) goto error; 957 958 if(sc_type == SCATTERING_IN_VOLUME) { 959 /* Sample a new optical path direction from the medium phase function */ 960 sample_scattering_direction(cmd, ithread, rng, &scattering, wlen, dir, wi); 961 962 } else { 963 double bounce_reflectivity; 964 double r; 965 ASSERT(sc_type == SCATTERING_AT_SURFACE); 966 967 /* Sample a new optical path direction from the surface BSDF */ 968 bounce_reflectivity = sample_bounce_direction 969 (cmd, ithread, rng, &hit_curr, wlen, dir, wi); 970 971 /* Russian roulette wrt surface scattering */ 972 r = ssp_rng_canonical(rng); 973 if(r > bounce_reflectivity) break; 974 975 /* Register the current hit to avoid a self intersection for the next 976 * optical path direction */ 977 hit_prev = hit_curr; 978 } 979 980 /* Update the optical path direction */ 981 dir[0] = wi[0]; 982 dir[1] = wi[1]; 983 dir[2] = wi[2]; 984 } 985 986 exit: 987 *out_weight = weight; 988 return res; 989 error: 990 weight = NaN; 991 goto exit; 992 }