htrdr_atmosphere_compute_radiance_lw.c (9169B)
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 27 #include "core/htrdr.h" 28 #include "core/htrdr_interface.h" 29 30 #include <high_tune/htsky.h> 31 32 #include <star/s3d.h> 33 #include <star/ssf.h> 34 #include <star/ssp.h> 35 #include <star/svx.h> 36 37 #include <rsys/double2.h> 38 #include <rsys/double3.h> 39 40 enum event { 41 EVENT_ABSORPTION, 42 EVENT_SCATTERING, 43 EVENT_NONE 44 }; 45 46 struct filter_context { 47 struct ssp_rng* rng; 48 const struct htsky* htsky; 49 size_t iband; /* Index of the spectral band */ 50 size_t iquad; /* Index of the quadrature point into the band */ 51 52 double Ts; /* Sampled optical thickness */ 53 double traversal_dst; /* Distance traversed along the ray */ 54 55 enum event event_type; 56 }; 57 static const struct filter_context FILTER_CONTEXT_NULL = { 58 NULL, NULL, 0, 0, 0.0, 0.0, EVENT_NONE 59 }; 60 61 /******************************************************************************* 62 * Helper functions 63 ******************************************************************************/ 64 static int 65 hit_filter 66 (const struct svx_hit* hit, 67 const double org[3], 68 const double dir[3], 69 const double range[2], 70 void* context) 71 { 72 struct filter_context* ctx = context; 73 double kext_max; 74 int pursue_traversal = 1; 75 ASSERT(hit && ctx && !SVX_HIT_NONE(hit) && org && dir && range); 76 (void)range; 77 78 kext_max = htsky_fetch_svx_voxel_property(ctx->htsky, HTSKY_Kext, 79 HTSKY_SVX_MAX, HTSKY_CPNT_MASK_ALL, ctx->iband, ctx->iquad, &hit->voxel); 80 81 ctx->traversal_dst = hit->distance[0]; 82 83 for(;;) { 84 double vox_dst = hit->distance[1] - ctx->traversal_dst; 85 const double T = vox_dst * kext_max; 86 87 /* A collision occurs behind `vox_dst' */ 88 if(ctx->Ts > T) { 89 ctx->Ts -= T; 90 ctx->traversal_dst = hit->distance[1]; 91 pursue_traversal = 1; 92 break; 93 94 /* A real/null collision occurs before `vox_dst' */ 95 } else { 96 const double collision_dst = ctx->Ts / kext_max; 97 double pos[3]; 98 double ks, ka; 99 double r; 100 double proba_abs; 101 double proba_sca; 102 103 /* Compute the traversed distance up to the challenged collision */ 104 ctx->traversal_dst += collision_dst; 105 ASSERT(ctx->traversal_dst >= hit->distance[0]); 106 ASSERT(ctx->traversal_dst <= hit->distance[1]); 107 108 /* Compute the world space position where a collision may occur */ 109 pos[0] = org[0] + ctx->traversal_dst * dir[0]; 110 pos[1] = org[1] + ctx->traversal_dst * dir[1]; 111 pos[2] = org[2] + ctx->traversal_dst * dir[2]; 112 113 ka = htsky_fetch_raw_property(ctx->htsky, HTSKY_Ka, HTSKY_CPNT_MASK_ALL, 114 ctx->iband, ctx->iquad, pos, -DBL_MAX, DBL_MAX); 115 ks = htsky_fetch_raw_property(ctx->htsky, HTSKY_Ks, HTSKY_CPNT_MASK_ALL, 116 ctx->iband, ctx->iquad, pos, -DBL_MAX, DBL_MAX); 117 118 r = ssp_rng_canonical(ctx->rng); 119 proba_abs = ka / kext_max; 120 proba_sca = ks / kext_max; 121 if(r < proba_abs) { /* Absorption event */ 122 pursue_traversal = 0; 123 ctx->event_type = EVENT_ABSORPTION; 124 break; 125 } else if(r < proba_abs + proba_sca) { /* Scattering event */ 126 pursue_traversal = 0; 127 ctx->event_type = EVENT_SCATTERING; 128 break; 129 } else { /* Null collision */ 130 ctx->Ts = ssp_ran_exp(ctx->rng, 1); /* Sample a new optical thickness */ 131 } 132 } 133 } 134 return pursue_traversal; 135 } 136 137 /******************************************************************************* 138 * Local functions 139 ******************************************************************************/ 140 double 141 atmosphere_compute_radiance_lw 142 (struct htrdr_atmosphere* cmd, 143 const size_t ithread, 144 struct ssp_rng* rng, 145 const double pos_in[3], 146 const double dir_in[3], 147 const double wlen, /* In nanometer */ 148 const size_t iband, 149 const size_t iquad) 150 { 151 struct s3d_hit s3d_hit = S3D_HIT_NULL; 152 struct s3d_hit s3d_hit_prev = S3D_HIT_NULL; 153 struct svx_hit svx_hit = SVX_HIT_NULL; 154 struct ssf_phase* phase_hg = NULL; 155 156 double wo[3]; 157 double pos[3]; 158 double dir[3]; 159 double range[2]; 160 double pos_next[3]; 161 double dir_next[3]; 162 double temperature; 163 double wlen_m = wlen * 1.e-9; 164 double g; 165 double w = 0; /* Weight */ 166 167 ASSERT(cmd && rng && pos_in && dir_in); 168 ASSERT(ithread < htrdr_get_threads_count(cmd->htrdr)); 169 170 /* Setup the phase function for this spectral band & quadrature point */ 171 CHK(RES_OK == ssf_phase_create 172 (htrdr_get_thread_allocator(cmd->htrdr, ithread), 173 &ssf_phase_hg, 174 &phase_hg)); 175 g = htsky_fetch_per_wavelength_particle_phase_function_asymmetry_parameter 176 (cmd->sky, wlen); 177 SSF(phase_hg_setup(phase_hg, g)); 178 179 /* Initialise the random walk */ 180 d3_set(pos, pos_in); 181 d3_set(dir, dir_in); 182 d2(range, 0, INF); 183 184 for(;;) { 185 struct filter_context ctx = FILTER_CONTEXT_NULL; 186 187 /* Find the first intersection with the surface geometry */ 188 d2(range, 0, DBL_MAX); 189 HTRDR(atmosphere_ground_trace_ray 190 (cmd->ground, pos, dir, range, &s3d_hit_prev, &s3d_hit)); 191 192 /* Sample an optical thickness */ 193 ctx.Ts = ssp_ran_exp(rng, 1); 194 195 /* Setup the remaining fields of the hit filter context */ 196 ctx.rng = rng; 197 ctx.htsky = cmd->sky; 198 ctx.iband = iband; 199 ctx.iquad = iquad; 200 201 /* Fit the ray range to the surface distance along the ray */ 202 d2(range, 0, s3d_hit.distance); 203 204 /* Trace a ray into the participating media */ 205 HTSKY(trace_ray(cmd->sky, pos, dir, range, NULL, 206 hit_filter, &ctx, iband, iquad, &svx_hit)); 207 208 /* No scattering and no surface reflection. 209 * Congratulation !! You are in space. */ 210 if(S3D_HIT_NONE(&s3d_hit) && SVX_HIT_NONE(&svx_hit)) { 211 w = 0; 212 break; 213 } 214 215 /* Compute the next position */ 216 pos_next[0] = pos[0] + dir[0]*ctx.traversal_dst; 217 pos_next[1] = pos[1] + dir[1]*ctx.traversal_dst; 218 pos_next[2] = pos[2] + dir[2]*ctx.traversal_dst; 219 220 /* Absorption event. Stop the realisation */ 221 if(ctx.event_type == EVENT_ABSORPTION) { 222 ASSERT(!SVX_HIT_NONE(&svx_hit)); 223 temperature = htsky_fetch_temperature(cmd->sky, pos_next); 224 /* weight is planck integrated over the spectral sub-interval */ 225 w = htrdr_planck_monochromatic(wlen_m, temperature); 226 break; 227 } 228 229 /* Negate the incoming dir to match the convention of the SSF library */ 230 d3_minus(wo, dir); 231 232 /* Scattering in the volume */ 233 if(ctx.event_type == EVENT_SCATTERING) { 234 ASSERT(!SVX_HIT_NONE(&svx_hit)); 235 ssf_phase_sample(phase_hg, rng, wo, dir_next, NULL); 236 s3d_hit_prev = S3D_HIT_NULL; 237 238 /* Scattering at a surface */ 239 } else { 240 struct htrdr_interface interf = HTRDR_INTERFACE_NULL; 241 const struct htrdr_mtl* mtl = NULL; 242 struct ssf_bsdf* bsdf = NULL; 243 double bounce_reflectivity = 0; 244 double N[3]; 245 int type; 246 ASSERT(ctx.event_type == EVENT_NONE); 247 ASSERT(!S3D_HIT_NONE(&s3d_hit)); 248 249 /* Fetch the hit interface materal and build its BSDF */ 250 htrdr_atmosphere_ground_get_interface(cmd->ground, &s3d_hit, &interf); 251 mtl = htrdr_interface_fetch_hit_mtl(&interf, dir, &s3d_hit); 252 HTRDR(mtl_create_bsdf(cmd->htrdr, mtl, ithread, wlen, rng, &bsdf)); 253 254 d3_normalize(N, d3_set_f3(N, s3d_hit.normal)); 255 if(d3_dot(N, wo) < 0) d3_minus(N, N); 256 257 bounce_reflectivity = ssf_bsdf_sample 258 (bsdf, rng, wo, N, dir_next, &type, NULL); 259 if(!(type & SSF_REFLECTION)) { /* Handle only reflections */ 260 bounce_reflectivity = 0; 261 } 262 263 /* Release the BSDF */ 264 SSF(bsdf_ref_put(bsdf)); 265 266 if(ssp_rng_canonical(rng) >= bounce_reflectivity) { /* Absorbed at boundary */ 267 temperature = mtl->temperature; /* Fetch mtl temperature */ 268 /* Weight is planck integrated over the spectral sub-interval */ 269 w = temperature>0 ? htrdr_planck_monochromatic(wlen_m, temperature) : 0; 270 break; 271 } 272 s3d_hit_prev = s3d_hit; 273 } 274 d3_set(pos, pos_next); 275 d3_set(dir, dir_next); 276 } 277 SSF(phase_ref_put(phase_hg)); 278 return w; 279 } 280