htrdr_planets_source.c (13677B)
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 "core/htrdr.h" 28 #include "core/htrdr_log.h" 29 30 #include <star/sbuf.h> 31 #include <star/ssp.h> 32 33 #include <rsys/algorithm.h> 34 #include <rsys/cstr.h> 35 #include <rsys/double3.h> 36 #include <rsys/ref_count.h> 37 38 typedef struct ALIGN(16) { 39 double wavelength; /* in nm */ 40 double radiance; /* in W/m²/sr/m */ 41 } source_radiance_T; 42 43 struct htrdr_planets_source { 44 double position[3]; /* In m */ 45 46 double radius; /* In m */ 47 48 /* In Kelvin. Defined if the radiances by wavelength is no set */ 49 double temperature; 50 51 struct sbuf* per_wlen_radiances; /* List of radiances by wavelength */ 52 53 ref_T ref; 54 struct htrdr* htrdr; 55 }; 56 57 /******************************************************************************* 58 * Helper functions 59 ******************************************************************************/ 60 static INLINE res_T 61 check_per_wlen_radiance_sbuf_desc 62 (const struct htrdr_planets_source* src, 63 const struct sbuf_desc* desc) 64 { 65 const source_radiance_T* spectrum = NULL; 66 size_t i; 67 ASSERT(src && desc); 68 69 /* Invalid size */ 70 if(desc->size == 0) { 71 htrdr_log_err(src->htrdr, "invalid empty source spectrum\n"); 72 return RES_BAD_ARG; 73 } 74 75 /* Invalid memory layout */ 76 if(desc->szitem != 16 || desc->alitem != 16 || desc->pitch != 16) { 77 htrdr_log_err(src->htrdr, "unexpected layout of source spectrum\n"); 78 return RES_BAD_ARG; 79 } 80 81 /* Data must be sorted */ 82 spectrum = desc->buffer; 83 FOR_EACH(i, 1, desc->size) { 84 if(spectrum[i-1].wavelength >= spectrum[i].wavelength) { 85 htrdr_log_err(src->htrdr, 86 "the source spectrum is not sorted in ascending order " 87 "with respect to wavelengths\n"); 88 return RES_BAD_ARG; 89 } 90 } 91 92 return RES_OK; 93 } 94 95 static res_T 96 setup_per_wavelength_radiances 97 (struct htrdr_planets_source* src, 98 const struct htrdr_planets_source_args* args) 99 { 100 struct sbuf_create_args sbuf_args; 101 struct sbuf_desc desc; 102 res_T res = RES_OK; 103 ASSERT(src && args && args->rnrl_filename && args->temperature < 0); 104 105 sbuf_args.logger = htrdr_get_logger(src->htrdr); 106 sbuf_args.allocator = htrdr_get_allocator(src->htrdr); 107 sbuf_args.verbose = htrdr_get_verbosity_level(src->htrdr); 108 res = sbuf_create(&sbuf_args, &src->per_wlen_radiances); 109 if(res != RES_OK) goto error; 110 111 res = sbuf_load(src->per_wlen_radiances, args->rnrl_filename); 112 if(res != RES_OK) goto error; 113 res = sbuf_get_desc(src->per_wlen_radiances, &desc); 114 if(res != RES_OK) goto error; 115 res = check_per_wlen_radiance_sbuf_desc(src, &desc); 116 if(res != RES_OK) goto error; 117 118 exit: 119 return res; 120 error: 121 htrdr_log_err(src->htrdr, "error loading %s -- %s\n", 122 args->rnrl_filename, res_to_cstr(res)); 123 goto exit; 124 } 125 126 static INLINE int 127 cmp_wlen(const void* a, const void* b) 128 { 129 const double wlen = *((double*)a); 130 const source_radiance_T* src_rad1 = b; 131 ASSERT(a && b); 132 133 if(wlen < src_rad1->wavelength) { 134 return -1; 135 } else if(wlen > src_rad1->wavelength) { 136 return +1; 137 } else { 138 return 0; 139 } 140 } 141 142 static double 143 get_radiance 144 (const struct htrdr_planets_source* src, 145 const double wlen) 146 { 147 struct sbuf_desc desc; 148 const source_radiance_T* spectrum; 149 const source_radiance_T* find; 150 size_t id; 151 ASSERT(src && src->per_wlen_radiances); 152 153 SBUF(get_desc(src->per_wlen_radiances, &desc)); 154 spectrum = desc.buffer; 155 156 if(wlen < spectrum[0].wavelength) { 157 htrdr_log_warn(src->htrdr, 158 "The wavelength %g nm is before the spectrum of the source\n", wlen); 159 return spectrum[0].radiance; 160 } 161 if(wlen > spectrum[desc.size-1].wavelength) { 162 htrdr_log_warn(src->htrdr, 163 "The wavelength %g nm is above the spectrum of the source\n", wlen); 164 return spectrum[desc.size-1].radiance; 165 } 166 167 /* Look for the first item whose wavelength is not less than 'wlen' */ 168 find = search_lower_bound(&wlen, spectrum, desc.size, desc.pitch, cmp_wlen); 169 ASSERT(find); 170 id = (size_t)(find - spectrum); 171 ASSERT(id < desc.size); 172 173 if(id == 0) { 174 ASSERT(wlen == spectrum[0].wavelength); 175 return spectrum[0].radiance; 176 } else { 177 const double w0 = spectrum[id-1].wavelength; 178 const double w1 = spectrum[id-0].wavelength; 179 const double L0 = spectrum[id-1].radiance; 180 const double L1 = spectrum[id-0].radiance; 181 const double u = (wlen - w0) / (w1 - w0); 182 const double L = L0 + u*(L1 - L0); /* Linear interpolation */ 183 return L; 184 } 185 } 186 187 static void 188 release_source(ref_T* ref) 189 { 190 struct htrdr_planets_source* source; 191 struct htrdr* htrdr; 192 ASSERT(ref); 193 194 source = CONTAINER_OF(ref, struct htrdr_planets_source, ref); 195 htrdr = source->htrdr; 196 if(source->per_wlen_radiances) SBUF(ref_put(source->per_wlen_radiances)); 197 MEM_RM(htrdr_get_allocator(htrdr), source); 198 htrdr_ref_put(htrdr); 199 } 200 201 /******************************************************************************* 202 * Local functions 203 ******************************************************************************/ 204 res_T 205 htrdr_planets_source_create 206 (struct htrdr* htrdr, 207 const struct htrdr_planets_source_args* args, 208 struct htrdr_planets_source** out_source) 209 { 210 struct htrdr_planets_source* src = NULL; 211 double dst; /* In m */ 212 double lat; /* In radians */ 213 double lon; /* In radians */ 214 res_T res = RES_OK; 215 ASSERT(htrdr && out_source); 216 ASSERT(htrdr_planets_source_args_check(args) == RES_OK); 217 218 src = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*src)); 219 if(!src) { 220 htrdr_log_err(htrdr, "error allocating source\n"); 221 res = RES_MEM_ERR; 222 goto error; 223 } 224 ref_init(&src->ref); 225 htrdr_ref_get(htrdr); 226 src->htrdr = htrdr; 227 src->radius = args->radius * 1e3/*From km to m*/; 228 229 if(!args->rnrl_filename) { 230 src->temperature = args->temperature; 231 } else { 232 res = setup_per_wavelength_radiances(src, args); 233 if(res != RES_OK) goto error; 234 src->temperature = -1; /* Not used */ 235 } 236 237 /* Convert latitude and longitude to radians and distance in m */ 238 lat = MDEG2RAD(args->latitude); 239 lon = MDEG2RAD(args->longitude); 240 dst = args->distance * 1e3/*From km to m*/; 241 242 /* Compute the position of the source */ 243 src->position[0] = dst * cos(lat) * cos(lon); 244 src->position[1] = dst * cos(lat) * sin(lon); 245 src->position[2] = dst * sin(lat); 246 247 exit: 248 *out_source = src; 249 return res; 250 error: 251 if(src) { htrdr_planets_source_ref_put(src); src = NULL; } 252 goto exit; 253 } 254 255 void 256 htrdr_planets_source_ref_get(struct htrdr_planets_source* source) 257 { 258 ASSERT(source); 259 ref_get(&source->ref); 260 } 261 262 void htrdr_planets_source_ref_put(struct htrdr_planets_source* source) 263 { 264 ASSERT(source); 265 ref_put(&source->ref, release_source); 266 } 267 268 double 269 htrdr_planets_source_sample_direction 270 (const struct htrdr_planets_source* source, 271 struct ssp_rng* rng, 272 const double pos[3], 273 double dir[3]) 274 { 275 double main_dir[3]; 276 double half_angle; /* In radians */ 277 double cos_half_angle; 278 double dst; /* In m */ 279 double pdf; 280 ASSERT(source && rng && pos && dir); 281 282 /* compute the direction of `pos' toward the center of the source */ 283 d3_sub(main_dir, source->position, pos); 284 285 /* Normalize the direction and keep the distance from `pos' to the center of 286 * the source */ 287 dst = d3_normalize(main_dir, main_dir); 288 CHK(dst > source->radius); 289 290 /* Sample the source according to its solid angle, 291 * i.e. 2*PI*(1 - cos(half_angle)) */ 292 half_angle = asin(source->radius/dst); 293 cos_half_angle = cos(half_angle); 294 ssp_ran_sphere_cap_uniform(rng, main_dir, cos_half_angle, dir, &pdf); 295 296 return pdf; 297 } 298 299 double /* In W/m²/sr/m */ 300 htrdr_planets_source_get_radiance 301 (const struct htrdr_planets_source* source, 302 const double wlen) 303 { 304 if(source->per_wlen_radiances) { 305 return get_radiance(source, wlen); 306 } else { 307 return htrdr_planck_monochromatic 308 (wlen*1e-9/*From nm to m*/, source->temperature); 309 } 310 } 311 312 double 313 htrdr_planets_source_distance_to 314 (const struct htrdr_planets_source* source, 315 const double pos[3]) 316 { 317 double vec[3]; 318 double dst; 319 ASSERT(source && pos); 320 321 d3_sub(vec, source->position, pos); 322 dst = d3_len(vec); 323 return dst - source->radius; 324 } 325 326 int 327 htrdr_planets_source_is_targeted 328 (const struct htrdr_planets_source* source, 329 const double pos[3], 330 const double dir[3]) 331 { 332 double main_dir[3]; 333 double half_angle; /* In radians */ 334 double dst; /* In m */ 335 ASSERT(source && dir && d3_is_normalized(dir)); 336 337 /* compute the direction of `pos' toward the center of the source */ 338 d3_sub(main_dir, source->position, pos); 339 340 /* Normalize the direction and keep the distance from `pos' to the center of 341 * the source */ 342 dst = d3_normalize(main_dir, main_dir); 343 CHK(dst > source->radius); 344 345 /* Compute the the half angle of the source as seen from pos */ 346 half_angle = asin(source->radius/dst); 347 348 return d3_dot(dir, main_dir) >= cos(half_angle); 349 } 350 351 res_T 352 htrdr_planets_source_get_spectral_range 353 (const struct htrdr_planets_source* source, 354 double range[2]) 355 { 356 res_T res = RES_OK; 357 ASSERT(source && range); 358 359 if(!source->per_wlen_radiances) { 360 range[0] = 0; 361 range[1] = INF; 362 } else { 363 struct sbuf_desc desc = SBUF_DESC_NULL; 364 const source_radiance_T* spectrum = NULL; 365 366 res = sbuf_get_desc(source->per_wlen_radiances, &desc); 367 if(res != RES_OK) goto error; 368 369 spectrum = desc.buffer; 370 range[0] = spectrum[0].wavelength; 371 range[1] = spectrum[desc.size-1].wavelength; 372 } 373 374 exit: 375 return res; 376 error: 377 goto exit; 378 } 379 380 int 381 htrdr_planets_source_does_radiance_vary_spectrally 382 (const struct htrdr_planets_source* source) 383 { 384 ASSERT(source); 385 return source->per_wlen_radiances != NULL; 386 } 387 388 res_T 389 htrdr_planets_source_get_spectrum 390 (const struct htrdr_planets_source* source, 391 const double range[2], /* In nm. Limits are inclusive */ 392 struct htrdr_planets_source_spectrum* source_spectrum) 393 { 394 double full_range[2]; 395 res_T res = RES_OK; 396 ASSERT(source && range && source_spectrum && range[0] <= range[1]); 397 398 if(!htrdr_planets_source_does_radiance_vary_spectrally(source)) { 399 res = RES_BAD_ARG; 400 goto error; 401 } 402 403 res = htrdr_planets_source_get_spectral_range(source, full_range); 404 if(res != RES_OK) goto error; 405 406 if(range[0] < full_range[0] || full_range[1] < range[1]) { 407 res = RES_BAD_ARG; 408 goto error; 409 } 410 411 source_spectrum->source = source; 412 source_spectrum->range[0] = range[0]; 413 source_spectrum->range[1] = range[1]; 414 415 if(range[0] == range[1]) { 416 /* Degenerated spectral range */ 417 source_spectrum->size = 1; 418 source_spectrum->buffer = NULL; 419 420 } else { 421 const source_radiance_T* spectrum; 422 const source_radiance_T* low; 423 const source_radiance_T* upp; 424 struct sbuf_desc desc; 425 426 res = sbuf_get_desc(source->per_wlen_radiances, &desc); 427 if(res != RES_OK) goto error; 428 429 spectrum = desc.buffer; 430 low = search_lower_bound(&range[0], spectrum, desc.size, desc.pitch, cmp_wlen); 431 upp = search_lower_bound(&range[1], spectrum, desc.size, desc.pitch, cmp_wlen); 432 ASSERT(low && upp); 433 434 if(low == upp) { 435 /* The range is fully included in a band */ 436 ASSERT(low->wavelength > range[0] && upp->wavelength >= range[1]); 437 source_spectrum->size = 2; 438 source_spectrum->buffer = NULL; 439 440 } else { 441 source_spectrum->size = 442 2/* Boundaries */ + (size_t)(upp - low)/*discrete items*/; 443 444 if(low->wavelength == range[0]) { 445 /* The lower limit coincide with a discrete element. 446 * Remove the discrete element */ 447 source_spectrum->size -= 1; 448 source_spectrum->buffer = low + 1; 449 } else { 450 source_spectrum->buffer = low; 451 } 452 453 } 454 } 455 456 exit: 457 return res; 458 error: 459 goto exit; 460 } 461 462 void 463 htrdr_planets_source_spectrum_at 464 (void* source_spectrum, 465 const size_t i, /* between [0, spectrum->size[ */ 466 double* wavelength, /* In nm */ 467 double* radiance) /* In W/m²/sr/m */ 468 { 469 struct htrdr_planets_source_spectrum* spectrum = source_spectrum; 470 ASSERT(spectrum && i < spectrum->size && wavelength && radiance); 471 472 /* Lower limit */ 473 if(i == 0) { 474 *wavelength = spectrum->range[0]; 475 *radiance = htrdr_planets_source_get_radiance 476 (spectrum->source, spectrum->range[0]); 477 478 /* Upper limit */ 479 } else if(i == spectrum->size-1) { 480 *wavelength = spectrum->range[1]; 481 *radiance = htrdr_planets_source_get_radiance 482 (spectrum->source, spectrum->range[1]); 483 484 /* Discrete element */ 485 } else { 486 const source_radiance_T* item = 487 (const source_radiance_T*)spectrum->buffer + (i-1); 488 *wavelength = item->wavelength; 489 *radiance = item->radiance; 490 } 491 }