sbb_ran_planck.c (10909B)
1 /* Copyright (C) 2018-2023 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #define _POSIX_C_SOURCE 200112L /* nextafter */ 17 18 #include "sbb.h" 19 #include "sbb_c.h" 20 21 #include <rsys/algorithm.h> 22 #include <rsys/cstr.h> 23 #include <rsys/dynamic_array_double.h> 24 #include <rsys/logger.h> 25 #include <rsys/mem_allocator.h> 26 #include <rsys/ref_count.h> 27 28 struct sbb_ran_planck { 29 struct darray_double pdf; 30 struct darray_double cdf; 31 double range[2]; /* Boundaries of the spectral integration interval [m] */ 32 double band_len; /* Length of a band [m] */ 33 double ref_temperature; /* Reference temperature [K] */ 34 size_t nbands; /* #bands */ 35 36 int verbose; /* Verbosity level */ 37 struct mem_allocator* allocator; 38 struct logger* logger; 39 ref_T ref; 40 }; 41 42 /******************************************************************************* 43 * Helper functions 44 ******************************************************************************/ 45 static res_T 46 check_ranst_planck_create_args(struct sbb_ran_planck_create_args* args) 47 { 48 struct logger* logger = NULL; 49 50 if(!args) return RES_BAD_ARG; 51 52 logger = args->logger ? args->logger : LOGGER_DEFAULT; 53 54 /* Invalid spectral range */ 55 if(args->range[0] < 0 56 || args->range[1] < 0 57 || args->range[0] > args->range[1]) { 58 if(args->verbose) { 59 logger_print(logger, LOG_ERROR, MSG_ERROR_PREFIX 60 "Invalid spectral range for the Planck distribution: [%g,%g] m\n", 61 args->range[0], args->range[1]); 62 } 63 return RES_BAD_ARG; 64 } 65 66 /* Invalid reference temperature */ 67 if(args->ref_temperature < 0) { 68 if(args->verbose) { 69 logger_print(logger, LOG_ERROR, MSG_ERROR_PREFIX 70 "Invalid reference temperature for the Planck distribution: %g K\n", 71 args->ref_temperature); 72 } 73 return RES_BAD_ARG; 74 } 75 76 return RES_OK; 77 } 78 79 static res_T 80 setup_cdf(struct sbb_ran_planck* planck) 81 { 82 double* pdf = NULL; 83 double* cdf = NULL; 84 double sum = 0; 85 size_t i; 86 res_T res = RES_OK; 87 ASSERT(planck && planck->nbands != 0); 88 89 res = darray_double_resize(&planck->cdf, planck->nbands); 90 if(res != RES_OK) { 91 logger_print(planck->logger, LOG_ERROR, MSG_ERROR_PREFIX 92 "Error when allocating the CDF of the planck distribution -- %s\n", 93 res_to_cstr(res)); 94 goto error; 95 } 96 97 res = darray_double_resize(&planck->pdf, planck->nbands); 98 if(res != RES_OK) { 99 logger_print(planck->logger, LOG_ERROR, MSG_ERROR_PREFIX 100 "Error when allocating the PDF of the planck distribution -- %s\n", 101 res_to_cstr(res)); 102 goto error; 103 } 104 105 cdf = darray_double_data_get(&planck->cdf); 106 pdf = darray_double_data_get(&planck->pdf); 107 108 /* Compute the *unnormalized* probability to sample a small band */ 109 FOR_EACH(i, 0, planck->nbands) { 110 double lo = planck->range[0] + (double)i * planck->band_len; 111 double hi = MMIN(lo + planck->band_len, planck->range[1]); 112 ASSERT(lo <= hi); 113 ASSERT(lo > planck->range[0] || eq_eps(lo, planck->range[0], 1.e-6)); 114 ASSERT(lo < planck->range[1] || eq_eps(lo, planck->range[1], 1.e-6)); 115 116 /* Compute the probability of the current band */ 117 pdf[i] = sbb_blackbody_fraction(lo, hi, planck->ref_temperature); 118 119 /* Update the norm */ 120 sum += pdf[i]; 121 } 122 123 /* Compute the cumulative of the previously computed probabilities */ 124 FOR_EACH(i, 0, planck->nbands) { 125 /* Normalize the probability */ 126 pdf[i] /= sum; 127 128 /* Setup the cumulative */ 129 if(i == 0) { 130 cdf[i] = pdf[i]; 131 } else { 132 cdf[i] = pdf[i] + cdf[i-1]; 133 ASSERT(cdf[i] >= cdf[i-1]); 134 } 135 } 136 ASSERT(eq_eps(cdf[planck->nbands-1], 1, 1.e-6)); 137 cdf[planck->nbands - 1] = 1.0; /* Handle numerical issue */ 138 139 exit: 140 return res; 141 error: 142 darray_double_clear(&planck->cdf); 143 darray_double_clear(&planck->pdf); 144 goto exit; 145 } 146 147 static double 148 sample_continue 149 (const struct sbb_ran_planck* planck, 150 const double r, /* In [0, 1[ */ 151 const double range[2], /* [m]*/ 152 double* pdf) /* May be NULL */ 153 { 154 /* Numerical parameters of the dichotomy search */ 155 const size_t MAX_ITER = 100; 156 const double EPSILON_LAMBDA_M = 1e-15; 157 const double EPSILON_BF = 1e-6; 158 159 /* Local variables */ 160 double bf = 0; 161 double bf_prev = 0; 162 double bf_min_max = 0; 163 double lambda = 0; 164 double lambda_prev = 0; 165 double lambda_min = 0; 166 double lambda_max = 0; 167 size_t i; 168 169 /* Check precondition */ 170 ASSERT(planck && range); 171 ASSERT(range[0] < range[1] && range[0] >= 0 && range[1] >= 0); 172 ASSERT(0 <= r && r < 1); 173 174 /* Setup the dichotomy search */ 175 lambda_min = range[0]; 176 lambda_max = range[1]; 177 bf_min_max = sbb_blackbody_fraction 178 (range[0], range[1], planck->ref_temperature); 179 180 /* Numerically search the lambda corresponding to the submitted canonical 181 * number */ 182 FOR_EACH(i, 0, MAX_ITER) { 183 double r_test; 184 lambda = (lambda_min + lambda_max) * 0.5; 185 bf = sbb_blackbody_fraction 186 (range[0], lambda, planck->ref_temperature); 187 188 r_test = bf / bf_min_max; 189 if(r_test < r) { 190 lambda_min = lambda; 191 } else { 192 lambda_max = lambda; 193 } 194 195 if(fabs(lambda_prev - lambda) < EPSILON_LAMBDA_M 196 || fabs(bf_prev - bf) < EPSILON_BF) 197 break; 198 199 lambda_prev = lambda; 200 bf_prev = bf; 201 } 202 if(i >= MAX_ITER && planck->verbose) { 203 logger_print(planck->logger, LOG_WARNING, MSG_WARNING_PREFIX 204 "Could not sample a wavelength wrt the Planck distribution " 205 "(spectral range = [%g, %g] m; reference temperature = %g K)\n", 206 range[0], range[1], planck->ref_temperature); 207 } 208 209 if(pdf) { 210 const double Tref = planck->ref_temperature; /* K */ 211 212 /* W/m²/sr/m */ 213 const double B_lambda = sbb_planck(lambda, lambda, Tref); 214 const double B_mean = sbb_planck(range[0], range[1], Tref); 215 216 *pdf = B_lambda / (B_mean * (range[1]-range[0])); 217 } 218 219 return lambda; 220 } 221 222 static FINLINE int 223 cmp_dbl(const void* a, const void* b) 224 { 225 const double d0 = *((const double*)a); 226 const double d1 = *((const double*)b); 227 return d0 < d1 ? -1 : (d0 > d1 ? 1 : 0); 228 } 229 230 static double 231 sample_discrete 232 (const struct sbb_ran_planck* planck, 233 const double r0, /* In [0, 1[ */ 234 const double r1, /* In [0, 1[ */ 235 double* pdf) /* May be NULL */ 236 { 237 const double* cdf = NULL; 238 const double* find = NULL; 239 double r0_next = nextafter(r0, DBL_MAX); 240 double band_range[2]; 241 double lambda = 0; 242 double pdf_continue = 0; 243 double pdf_band = 0; 244 size_t cdf_length = 0; 245 size_t i; 246 ASSERT(planck && planck->nbands != 0); 247 ASSERT(0 <= r0 && r0 < 1); 248 ASSERT(0 <= r1 && r1 < 1); 249 (void)pdf_band; 250 251 cdf = darray_double_cdata_get(&planck->cdf); 252 cdf_length = darray_double_size_get(&planck->cdf); 253 ASSERT(cdf_length > 0); 254 255 /* Use r_next rather than r0 in order to find the first entry that is not less 256 * than *or equal* to r0 */ 257 find = search_lower_bound(&r0_next, cdf, cdf_length, sizeof(double), cmp_dbl); 258 ASSERT(find); 259 260 i = (size_t)(find - cdf); 261 ASSERT(i < cdf_length && cdf[i] > r0 && (!i || cdf[i-1] <= r0)); 262 263 band_range[0] = planck->range[0] + (double)i*planck->band_len; 264 band_range[1] = band_range[0] + planck->band_len; 265 266 /* Fetch the pdf of the sampled band */ 267 pdf_band = darray_double_cdata_get(&planck->pdf)[i]; 268 269 /* Uniformly sample a wavelength in the sampled band */ 270 lambda = band_range[0] + (band_range[1] - band_range[0]) * r1; 271 pdf_continue = 1.0 / (band_range[1] - band_range[0]); 272 273 if(pdf) { 274 *pdf = pdf_band * pdf_continue; 275 } 276 277 return lambda; 278 } 279 280 static void 281 release_planck(ref_T* ref) 282 { 283 struct sbb_ran_planck* planck = CONTAINER_OF(ref, struct sbb_ran_planck, ref); 284 ASSERT(ref); 285 darray_double_release(&planck->cdf); 286 darray_double_release(&planck->pdf); 287 MEM_RM(planck->allocator, planck); 288 } 289 290 /******************************************************************************* 291 * Exported functions 292 ******************************************************************************/ 293 res_T 294 sbb_ran_planck_create 295 (struct sbb_ran_planck_create_args* args, 296 struct sbb_ran_planck** out_planck) 297 { 298 struct sbb_ran_planck* planck = NULL; 299 struct mem_allocator* allocator = NULL; 300 struct logger* logger = NULL; 301 res_T res = RES_OK; 302 303 if(!out_planck) return RES_BAD_ARG; 304 res = check_ranst_planck_create_args(args); 305 if(res != RES_OK) goto error; 306 307 allocator = args->allocator ? args->allocator : &mem_default_allocator; 308 logger = args->logger ? args->logger : LOGGER_DEFAULT; 309 310 planck = MEM_CALLOC(allocator, 1, sizeof(*planck)); 311 if(!planck) { 312 res = RES_MEM_ERR; 313 if(args->verbose) { 314 logger_print(logger, LOG_ERROR, 315 "Planck distribution allocation error -- %s\n", 316 res_to_cstr(res)); 317 } 318 goto error; 319 } 320 planck->allocator = allocator; 321 planck->logger = logger; 322 ref_init(&planck->ref); 323 darray_double_init(planck->allocator, &planck->cdf); 324 darray_double_init(planck->allocator, &planck->pdf); 325 planck->range[0] = args->range[0]; 326 planck->range[1] = args->range[1]; 327 planck->ref_temperature = args->ref_temperature; 328 planck->nbands = args->nbands; 329 planck->verbose = args->verbose; 330 331 if(planck->nbands == 0) { 332 planck->band_len = 0; 333 } else { 334 planck->band_len = 335 (planck->range[1] - planck->range[0]) 336 / (double)planck->nbands; 337 338 res = setup_cdf(planck); 339 if(res != RES_OK) goto error; 340 } 341 342 exit: 343 if(out_planck) *out_planck = planck; 344 return res; 345 error: 346 if(planck) { SBB(ran_planck_ref_put(planck)); planck = NULL; } 347 goto exit; 348 } 349 350 res_T 351 sbb_ran_planck_ref_get(struct sbb_ran_planck* planck) 352 { 353 if(!planck) return RES_BAD_ARG; 354 ref_get(&planck->ref); 355 return RES_OK; 356 } 357 358 res_T 359 sbb_ran_planck_ref_put(struct sbb_ran_planck* planck) 360 { 361 if(!planck) return RES_BAD_ARG; 362 ref_put(&planck->ref, release_planck); 363 return RES_OK; 364 } 365 366 double /* Wavelength [m] */ 367 sbb_ran_planck_sample 368 (struct sbb_ran_planck* planck, 369 const double r0, /* Random number in [0, 1[ */ 370 const double r1, /* Random number in [0, 1[ */ 371 double* pdf) /* [m^-1]. May be NULL */ 372 { 373 ASSERT(planck && r0 >= 0 && r1 >= 0 && r0 < 1 && r1 < 1); 374 375 if(eq_eps(planck->range[0], planck->range[1], 1.e-12 /* <=> 1 pm */)) { 376 if(pdf) *pdf = 1; 377 return planck->range[0]; 378 379 } else if(planck->nbands == 0) { 380 return sample_continue(planck, r0, planck->range, pdf); 381 382 } else { /* Speed up wavelength sampling by presampling a band */ 383 return sample_discrete(planck, r0, r1, pdf); 384 } 385 }