htrdr_ran_wlen_planck.c (11228B)
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 /* nextafter */ 25 26 #include "core/htrdr.h" 27 #include "core/htrdr_c.h" 28 #include "core/htrdr_log.h" 29 #include "core/htrdr_ran_wlen_planck.h" 30 #include "core/htrdr_spectral.h" 31 32 #include <rsys/algorithm.h> 33 #include <rsys/cstr.h> 34 #include <rsys/dynamic_array_double.h> 35 #include <rsys/mem_allocator.h> 36 #include <rsys/ref_count.h> 37 38 #include <math.h> /* nextafter */ 39 40 struct htrdr_ran_wlen_planck { 41 struct darray_double pdf; 42 struct darray_double cdf; 43 double range[2]; /* Boundaries of the spectral integration interval */ 44 double band_len; /* Length in nanometers of a band */ 45 double ref_temperature; /* In Kelvin */ 46 size_t nbands; /* # bands */ 47 48 struct htrdr* htrdr; 49 ref_T ref; 50 }; 51 52 /******************************************************************************* 53 * Helper functions 54 ******************************************************************************/ 55 static res_T 56 setup_wlen_planck_ran_cdf 57 (struct htrdr_ran_wlen_planck* planck, 58 const char* func_name) 59 { 60 double* pdf = NULL; 61 double* cdf = NULL; 62 double sum = 0; 63 size_t i; 64 res_T res = RES_OK; 65 ASSERT(planck && func_name); 66 ASSERT(planck->nbands != HTRDR_WLEN_RAN_PLANCK_CONTINUE); 67 68 res = darray_double_resize(&planck->cdf, planck->nbands); 69 if(res != RES_OK) { 70 htrdr_log_err(planck->htrdr, 71 "%s: Error allocating the CDF of the spectral bands -- %s.\n", 72 func_name, res_to_cstr(res)); 73 goto error; 74 } 75 res = darray_double_resize(&planck->pdf, planck->nbands); 76 if(res != RES_OK) { 77 htrdr_log_err(planck->htrdr, 78 "%s: Error allocating the pdf of the spectral bands -- %s.\n", 79 func_name, res_to_cstr(res)); 80 goto error; 81 } 82 83 cdf = darray_double_data_get(&planck->cdf); 84 pdf = darray_double_data_get(&planck->pdf); /* Now save pdf to correct weight */ 85 86 htrdr_log(planck->htrdr, 87 "Number of bands used to speed up Planck distribution: %lu\n", 88 (unsigned long)planck->nbands); 89 90 /* Compute the *unnormalized* probability to sample a small band */ 91 FOR_EACH(i, 0, planck->nbands) { 92 double lambda_lo = planck->range[0] + (double)i * planck->band_len; 93 double lambda_hi = MMIN(lambda_lo + planck->band_len, planck->range[1]); 94 ASSERT(lambda_lo<= lambda_hi); 95 ASSERT(lambda_lo > planck->range[0] 96 || eq_eps(lambda_lo, planck->range[0], 1.e-6)); 97 ASSERT(lambda_lo < planck->range[1] 98 || eq_eps(lambda_lo, planck->range[1], 1.e-6)); 99 100 /* Convert from nanometer to meter */ 101 lambda_lo *= 1.e-9; 102 lambda_hi *= 1.e-9; 103 104 /* Compute the probability of the current band */ 105 pdf[i] = htrdr_blackbody_fraction 106 (lambda_lo, lambda_hi, planck->ref_temperature); 107 108 /* Update the norm */ 109 sum += pdf[i]; 110 } 111 112 /* Compute the cumulative of the previously computed probabilities */ 113 FOR_EACH(i, 0, planck->nbands) { 114 /* Normalize the probability */ 115 pdf[i] /= sum; 116 117 /* Setup the cumulative */ 118 if(i == 0) { 119 cdf[i] = pdf[i]; 120 } else { 121 cdf[i] = pdf[i] + cdf[i-1]; 122 ASSERT(cdf[i] >= cdf[i-1]); 123 } 124 } 125 ASSERT(eq_eps(cdf[planck->nbands-1], 1, 1.e-6)); 126 cdf[planck->nbands - 1] = 1.0; /* Handle numerical issue */ 127 128 exit: 129 return res; 130 error: 131 darray_double_clear(&planck->cdf); 132 darray_double_clear(&planck->pdf); 133 goto exit; 134 } 135 136 static double 137 wlen_ran_sample_continue 138 (const struct htrdr_ran_wlen_planck* planck, 139 const double r, 140 const double range[2], /* In nanometer */ 141 const char* func_name, 142 double* pdf) 143 { 144 /* Numerical parameters of the dichotomy search */ 145 const size_t MAX_ITER = 100; 146 const double EPSILON_LAMBDA_M = 1e-15; 147 const double EPSILON_BF = 1e-6; 148 149 /* Local variables */ 150 double bf = 0; 151 double bf_prev = 0; 152 double bf_min_max = 0; 153 double lambda_m = 0; 154 double lambda_m_prev = 0; 155 double lambda_m_min = 0; 156 double lambda_m_max = 0; 157 double range_m[2] = {0, 0}; 158 size_t i; 159 160 /* Check precondition */ 161 ASSERT(planck && func_name); 162 ASSERT(range && range[0] < range[1]); 163 ASSERT(0 <= r && r < 1); 164 165 /* Convert the wavelength range in meters */ 166 range_m[0] = range[0] * 1.0e-9; 167 range_m[1] = range[1] * 1.0e-9; 168 169 /* Setup the dichotomy search */ 170 lambda_m_min = range_m[0]; 171 lambda_m_max = range_m[1]; 172 bf_min_max = htrdr_blackbody_fraction 173 (range_m[0], range_m[1], planck->ref_temperature); 174 175 /* Numerically search the lambda corresponding to the submitted canonical 176 * number */ 177 FOR_EACH(i, 0, MAX_ITER) { 178 double r_test; 179 lambda_m = (lambda_m_min + lambda_m_max) * 0.5; 180 bf = htrdr_blackbody_fraction 181 (range_m[0], lambda_m, planck->ref_temperature); 182 183 r_test = bf / bf_min_max; 184 if(r_test < r) { 185 lambda_m_min = lambda_m; 186 } else { 187 lambda_m_max = lambda_m; 188 } 189 190 if(fabs(lambda_m_prev - lambda_m) < EPSILON_LAMBDA_M 191 || fabs(bf_prev - bf) < EPSILON_BF) 192 break; 193 194 lambda_m_prev = lambda_m; 195 bf_prev = bf; 196 } 197 if(i >= MAX_ITER) { 198 htrdr_log_warn(planck->htrdr, 199 "%s: could not sample a wavelength in the range [%g, %g] nanometers " 200 "for the reference temperature %g Kelvin.\n", 201 func_name, SPLIT2(range), planck->ref_temperature); 202 } 203 204 if(pdf) { 205 const double Tref = planck->ref_temperature; /* K */ 206 207 /* W/m²/sr/m */ 208 const double B_lambda = htrdr_planck(lambda_m, lambda_m, Tref); 209 const double B_mean = htrdr_planck(range_m[0], range_m[1], Tref); 210 211 *pdf = B_lambda / (B_mean * (range_m[1]-range_m[0])); 212 *pdf *= 1.e-9; /* Transform from m⁻¹ to nm⁻¹ */ 213 } 214 215 return lambda_m * 1.e9; /* Convert in nanometers */ 216 } 217 218 static double 219 wlen_ran_sample_discrete 220 (const struct htrdr_ran_wlen_planck* planck, 221 const double r0, 222 const double r1, 223 const char* func_name, 224 double* pdf) 225 { 226 const double* cdf = NULL; 227 const double* find = NULL; 228 double r0_next = nextafter(r0, DBL_MAX); 229 double band_range[2]; 230 double lambda = 0; 231 double pdf_continue = 0; 232 double pdf_band = 0; 233 size_t cdf_length = 0; 234 size_t i; 235 ASSERT(planck && planck->nbands != HTRDR_WLEN_RAN_PLANCK_CONTINUE); 236 ASSERT(0 <= r0 && r0 < 1); 237 ASSERT(0 <= r1 && r1 < 1); 238 (void)func_name; 239 (void)pdf_band; 240 241 cdf = darray_double_cdata_get(&planck->cdf); 242 cdf_length = darray_double_size_get(&planck->cdf); 243 ASSERT(cdf_length > 0); 244 245 /* Use r_next rather than r0 in order to find the first entry that is not less 246 * than *or equal* to r0 */ 247 find = search_lower_bound(&r0_next, cdf, cdf_length, sizeof(double), cmp_dbl); 248 ASSERT(find); 249 250 i = (size_t)(find - cdf); 251 ASSERT(i < cdf_length && cdf[i] > r0 && (!i || cdf[i-1] <= r0)); 252 253 band_range[0] = planck->range[0] + (double)i*planck->band_len; 254 band_range[1] = band_range[0] + planck->band_len; 255 256 /* Fetch the pdf of the sampled band */ 257 pdf_band = darray_double_cdata_get(&planck->pdf)[i]; 258 259 /* Uniformly sample a wavelength in the sampled band */ 260 lambda = band_range[0] + (band_range[1] - band_range[0]) * r1; 261 pdf_continue = 1.0 / (band_range[1] - band_range[0]); 262 263 if(pdf) { 264 *pdf = pdf_band * pdf_continue; 265 } 266 267 return lambda; 268 } 269 270 static void 271 release_wlen_planck_ran(ref_T* ref) 272 { 273 struct htrdr_ran_wlen_planck* planck = NULL; 274 struct htrdr* htrdr = NULL; 275 ASSERT(ref); 276 planck = CONTAINER_OF(ref, struct htrdr_ran_wlen_planck, ref); 277 darray_double_release(&planck->cdf); 278 darray_double_release(&planck->pdf); 279 htrdr = planck->htrdr; 280 MEM_RM(htrdr_get_allocator(htrdr), planck); 281 htrdr_ref_put(planck->htrdr); 282 } 283 284 /******************************************************************************* 285 * Local functions 286 ******************************************************************************/ 287 res_T 288 htrdr_ran_wlen_planck_create 289 (struct htrdr* htrdr, 290 /* range must be included in [200,1000] nm for shortwave or in [1000,100000] 291 * nanometers for longwave (thermal) */ 292 const double range[2], 293 const size_t nbands, /* # bands used to discretized CDF */ 294 const double ref_temperature, 295 struct htrdr_ran_wlen_planck** out_wlen_planck_ran) 296 { 297 struct htrdr_ran_wlen_planck* planck = NULL; 298 res_T res = RES_OK; 299 ASSERT(htrdr && range && out_wlen_planck_ran && ref_temperature > 0); 300 ASSERT(ref_temperature > 0); 301 ASSERT(range[0] <= range[1]); 302 303 planck = MEM_CALLOC(htrdr->allocator, 1, sizeof(*planck)); 304 if(!planck) { 305 res = RES_MEM_ERR; 306 htrdr_log_err(htrdr, 307 "%s: could not allocate Planck distribution -- %s.\n", 308 FUNC_NAME, res_to_cstr(res)); 309 goto error; 310 } 311 ref_init(&planck->ref); 312 darray_double_init(htrdr->allocator, &planck->cdf); 313 darray_double_init(htrdr->allocator, &planck->pdf); 314 htrdr_ref_get(htrdr); 315 planck->htrdr = htrdr; 316 317 planck->range[0] = range[0]; 318 planck->range[1] = range[1]; 319 planck->ref_temperature = ref_temperature; 320 planck->nbands = nbands; 321 322 if(nbands == HTRDR_WLEN_RAN_PLANCK_CONTINUE) { 323 planck->band_len = 0; 324 } else { 325 planck->band_len = (range[1] - range[0]) / (double)planck->nbands; 326 327 res = setup_wlen_planck_ran_cdf(planck, FUNC_NAME); 328 if(res != RES_OK) goto error; 329 } 330 331 htrdr_log(htrdr, "Spectral interval defined on [%g, %g] nanometers.\n", 332 range[0], range[1]); 333 334 exit: 335 *out_wlen_planck_ran = planck; 336 return res; 337 error: 338 if(planck) { 339 htrdr_ran_wlen_planck_ref_put(planck); 340 planck = NULL; 341 } 342 goto exit; 343 } 344 345 void 346 htrdr_ran_wlen_planck_ref_get(struct htrdr_ran_wlen_planck* planck) 347 { 348 ASSERT(planck); 349 ref_get(&planck->ref); 350 } 351 352 void 353 htrdr_ran_wlen_planck_ref_put(struct htrdr_ran_wlen_planck* planck) 354 { 355 ASSERT(planck); 356 ref_put(&planck->ref, release_wlen_planck_ran); 357 } 358 359 double 360 htrdr_ran_wlen_planck_sample 361 (const struct htrdr_ran_wlen_planck* planck, 362 const double r0, 363 const double r1, 364 double* pdf) 365 { 366 ASSERT(planck); 367 if(planck->nbands != HTRDR_WLEN_RAN_PLANCK_CONTINUE) { /* Discrete */ 368 return wlen_ran_sample_discrete(planck, r0, r1, FUNC_NAME, pdf); 369 } else if(eq_eps(planck->range[0], planck->range[1], 1.e-6)) { 370 if(pdf) *pdf = 1; 371 return planck->range[0]; 372 } else { /* Continue */ 373 return wlen_ran_sample_continue 374 (planck, r0, planck->range, FUNC_NAME, pdf); 375 } 376 } 377