sphin_source_surface.c (14253B)
1 /* Copyright (C) 2024-2026 Centre National de la Recherche Scientifique 2 * Copyright (C) 2024-2026 Clermont Auvergne INP 3 * Copyright (C) 2024-2026 INSA Lyon 4 * Copyright (C) 2024-2026 Institut Mines Télécom Albi-Carmaux 5 * Copyright (C) 2024-2026 Institut National Polytechnique de Toulouse 6 * Copyright (C) 2024-2026 |Méso|Star> (contact@meso-star.com) 7 * Copyright (C) 2024-2026 PhotonLyX (info@photonlyx.com) 8 * Copyright (C) 2024-2026 Université de Lorraine 9 * Copyright (C) 2024-2026 Université Paul Sabatier 10 * Copyright (C) 2024-2026 Université Toulouse - Jean Jaurès 11 * 12 * This program is free software: you can redistribute it and/or modify 13 * it under the terms of the GNU General Public License as published by 14 * the Free Software Foundation, either version 3 of the License, or 15 * (at your option) any later version. 16 * 17 * This program is distributed in the hope that it will be useful, 18 * but WITHOUT ANY WARRANTY; without even the implied warranty of 19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 20 * GNU General Public License for more details. 21 * 22 * You should have received a copy of the GNU General Public License 23 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 24 #define _POSIX_C_SOURCE 200112L /* for strtok_r support */ 25 26 #include "sphin.h" 27 #include "sphin_c.h" 28 #include "sphin_config.h" 29 #include "sphin_spectral_property.h" 30 #include "sphin_source_surface.h" 31 #include "sphin_surface.h" 32 33 #include <rsys/cstr.h> 34 #include <rsys/double3.h> 35 #include <rsys/mem_allocator.h> 36 #include <rsys/ref_count.h> 37 #include <rsys/rsys.h> 38 #include <rsys/str.h> 39 #include <rsys/text_reader.h> 40 41 struct sphin_source_surface { 42 struct str name; 43 struct sphin_source_surface_direction_distribution direction_distribution; 44 struct sphin_source_surface_flux_density flux_density; 45 46 struct sphin* sphin; 47 ref_T ref; 48 }; 49 50 /******************************************************************************* 51 * Helper functions 52 ******************************************************************************/ 53 static res_T 54 source_surface_create 55 (struct sphin* sphin, 56 const char* name, 57 struct sphin_source_surface** out_source) 58 { 59 struct sphin_source_surface* source = NULL; 60 res_T res = RES_OK; 61 62 ASSERT(NULL != sphin); 63 ASSERT(NULL != out_source); 64 ASSERT(NULL != name); 65 ASSERT('\0' != name[0]); /* Name can't be empty */ 66 67 source = MEM_CALLOC(sphin->allocator, 1, sizeof(struct sphin_source_surface)); 68 if (NULL == source) { res = RES_MEM_ERR; goto error; } 69 ref_init(&source->ref); 70 SPHIN(ref_get(sphin)); 71 source->sphin = sphin; 72 source->direction_distribution = 73 SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_NULL; 74 source->flux_density = SPHIN_SOURCE_SURFACE_FLUX_DENSITY_NULL; 75 76 str_init(sphin->allocator, &source->name); 77 res = str_set(&source->name, name); 78 if (RES_OK != res) { goto error; } 79 80 exit: 81 *out_source = source; 82 return res; 83 error: 84 if (NULL != source) { 85 SPHIN(source_surface_ref_put(source)); 86 source = NULL; 87 } 88 goto exit; 89 } 90 91 static void 92 release_source_surface 93 (ref_T* address) 94 { 95 struct sphin_source_surface* source = NULL; 96 struct sphin* sphin = NULL; 97 98 ASSERT(NULL != address); 99 100 source = CONTAINER_OF(address, struct sphin_source_surface, ref); 101 str_release(&source->name); 102 if (NULL != source->flux_density.emission_spectrum) { 103 SPHIN(spectral_property_ref_put(source->flux_density.emission_spectrum)); 104 } 105 sphin = source->sphin; 106 MEM_RM(sphin->allocator, source); 107 SPHIN(ref_put(sphin)); 108 } 109 110 static res_T 111 parse_lambertian_direction_distribution 112 (struct sphin_source_surface_direction_distribution* dir_dist, 113 struct txtrdr* txtrdr, 114 char* token) 115 { 116 char* token_ptr = NULL; 117 118 ASSERT(NULL != dir_dist); 119 ASSERT(NULL != txtrdr); 120 121 (void)txtrdr; 122 123 token = strtok_r(token, " \t", &token_ptr); 124 if (NULL != token) { return RES_BAD_ARG; } 125 126 dir_dist->type = SPHIN_SOURCE_DIRECTION_ISOTROPIC; 127 return RES_OK; 128 } 129 130 static res_T 131 parse_collim_direction_distribution 132 (struct sphin_source_surface_direction_distribution* dir_dist, 133 struct txtrdr* txtrdr, 134 char* token) 135 { 136 char* str_direction = NULL; 137 char* token_ptr = NULL; 138 double direction[3] = {0}; 139 res_T res = RES_OK; 140 141 (void)txtrdr; 142 143 ASSERT(NULL != dir_dist); 144 ASSERT(NULL != txtrdr); 145 146 if (NULL == token) { res = RES_BAD_ARG; goto error; } 147 148 dir_dist->type = SPHIN_SOURCE_DIRECTION_COLLIM; 149 str_direction = strtok_r(token, " \t", &token_ptr); 150 if (0 == strcmp(str_direction, "NORMAL")){ 151 direction[0] = direction[1] = direction[2] = 0; 152 } 153 else{ 154 res = RES_BAD_ARG; goto error; 155 } 156 157 d3_set(dir_dist->collim.direction, direction); 158 159 exit: 160 return res; 161 error: 162 goto exit; 163 } 164 165 static res_T 166 parse_cos_pow_n_direction_distribution 167 (struct sphin_source_surface_direction_distribution* dir_dist, 168 struct txtrdr* txtrdr, 169 char* token) 170 { 171 char* str_collimation_degree = NULL; 172 char* token_ptr = NULL; 173 double collimation_degree; 174 res_T res = RES_OK; 175 176 (void)txtrdr; 177 178 ASSERT(NULL != dir_dist); 179 ASSERT(NULL != txtrdr); 180 181 if (NULL == token) { res = RES_BAD_ARG; goto error; } 182 183 dir_dist->type = SPHIN_SOURCE_DIRECTION_COS_POW_N; 184 185 str_collimation_degree = strtok_r(token, " \t", &token_ptr); 186 res = cstr_to_double(str_collimation_degree, &collimation_degree); 187 if (RES_OK != res) { res = RES_BAD_ARG; goto error; } 188 if (collimation_degree < 0) { res = RES_BAD_ARG; goto error; } 189 dir_dist->cos_pow_n.collimation_degree = collimation_degree; 190 191 exit: 192 return res; 193 error: 194 goto exit; 195 } 196 197 static res_T 198 parse_direction_distribution 199 (struct sphin_source_surface* source, 200 struct txtrdr* txtrdr, 201 char* value) 202 { 203 char* direction_distribution_type = NULL; 204 char* token_ptr = NULL; 205 res_T res = RES_OK; 206 207 ASSERT(NULL != source); 208 ASSERT(NULL != txtrdr); 209 210 if (NULL == value) { res = RES_BAD_ARG; goto error; }; 211 212 /* Parse direction distribution type */ 213 direction_distribution_type = strtok_r(value, " \t", &token_ptr); 214 if (NULL == direction_distribution_type){ res = RES_BAD_ARG; goto error; } 215 216 /* Lambertian source */ 217 if (0 == strcmp(direction_distribution_type, "LAMBERT")){ 218 res = parse_lambertian_direction_distribution( 219 &source->direction_distribution, 220 txtrdr, 221 token_ptr); 222 } 223 224 /* Collimated source */ 225 else if (0 == strcmp(direction_distribution_type, "COLLIM")){ 226 res = parse_collim_direction_distribution( 227 &source->direction_distribution, 228 txtrdr, 229 token_ptr); 230 } 231 232 /* Cos^n model of direction distribution */ 233 else if (0 == strcmp(direction_distribution_type, "COS_POW_N")){ 234 res = parse_cos_pow_n_direction_distribution( 235 &source->direction_distribution, 236 txtrdr, 237 token_ptr); 238 } 239 else { res = RES_BAD_ARG; goto error; } 240 if (RES_OK != res) {goto error;} 241 242 res = txtrdr_read_line(txtrdr); 243 if (RES_OK != res) { goto error; } 244 245 exit: 246 return res; 247 error: 248 goto exit; 249 } 250 251 static res_T 252 parse_flux_density_unit 253 (struct sphin_source_surface* source, 254 struct txtrdr* txtrdr, 255 char* value) 256 { 257 res_T res = RES_OK; 258 259 ASSERT(NULL != source); 260 ASSERT(NULL != txtrdr); 261 ASSERT(NULL != value); 262 263 (void)txtrdr; 264 265 if (0 == strcmp(value, "mol/m^2/s") 266 || 0 == strcmp(value, "mol.m^-2.s^-1")) { 267 source->flux_density.unit= SPHIN_PHOTON_UNIT_MOL; 268 source->flux_density.flux_density *= 1e6; /* From mol to umol */ 269 } 270 else if (0 == strcmp(value, "umol/m^2/s") 271 || 0 == strcmp(value, "umol.m^-2.s^-1")) { 272 source->flux_density.unit= SPHIN_PHOTON_UNIT_MOL; 273 source->flux_density.flux_density *= 1; /* No conversion */ 274 } 275 276 /* Energy flux density energy */ 277 else if (0 == strcmp(value, "mW/m^2") 278 || 0 == strcmp(value, "mW.m^-2")) { 279 source->flux_density.unit= SPHIN_PHOTON_UNIT_JOULE; 280 source->flux_density.flux_density *= 1e-3; /* From mWatt to Watt*/ 281 } 282 else if (0 == strcmp(value, "W/m^2") 283 || 0 == strcmp(value, "W.m^-2") 284 || 0 == strcmp(value, "J.m^-2.s^-1") 285 || 0 == strcmp(value, "J/m^2/s")) { 286 source->flux_density.unit= SPHIN_PHOTON_UNIT_JOULE; 287 source->flux_density.flux_density *= 1; /* No conversion */ 288 } 289 else { res = RES_BAD_ARG; goto error; } 290 291 exit: 292 return res; 293 error: 294 goto exit; 295 } 296 297 static res_T 298 normalize_spectrum 299 (struct sphin_spectral_property* spectrum) 300 { 301 double sum = 0; 302 double* wavelengths = NULL; 303 double* values = NULL; 304 size_t data_count = 0; 305 size_t i = 0; 306 res_T res = RES_OK; 307 308 wavelengths = darray_double_data_get(&spectrum->wavelengths); 309 values = darray_double_data_get(&spectrum->values); 310 data_count = darray_double_size_get(&spectrum->wavelengths); 311 312 /* Compute the integral of the spectrum over the wavelengths domain using 313 * the trapeze method */ 314 315 FOR_EACH(i, 0, data_count-1) { 316 double dx = wavelengths[i+1] - wavelengths[i]; /* Interval size */ 317 sum += 0.5 * dx * (values[i] + values[i+1]); 318 } 319 320 /* Normalize spectrum */ 321 FOR_EACH(i, 0, data_count-1) { 322 values[i] = values[i] / sum; 323 } 324 325 return res; 326 } 327 328 static res_T 329 parse_flux_density 330 (struct sphin_source_surface* source, 331 struct txtrdr* txtrdr, 332 char* value) 333 { 334 char* tokens[5] = {NULL}; /* str_flux_density, flux_density_unit, [filename], 335 [spectral_unit], [property_unit] */ 336 char* filename = NULL; 337 char* flux_density_unit = NULL; 338 char* property_unit = NULL; 339 char* spectral_unit = NULL; 340 char* str_flux_density = NULL; 341 char* token = NULL; 342 char* token_ptr = NULL; 343 double flux_density = 0; 344 size_t token_count = 0; 345 res_T res = RES_OK; 346 347 ASSERT(NULL != source); 348 ASSERT(NULL != txtrdr); 349 350 if (NULL == value) { res = RES_BAD_ARG; goto error; } 351 352 token = strtok_r(value, " \t", &token_ptr); 353 while (token != NULL && token_count < 5) { 354 tokens[token_count++] = trim_keyword(token); 355 token = strtok_r(NULL, " \t", &token_ptr); 356 } 357 358 /* The number of tokens must be 5 */ 359 if (token_count != 5) { res = RES_BAD_ARG; goto error; } 360 361 /* Parse flux density value */ 362 str_flux_density = tokens[0]; 363 res = cstr_to_double(str_flux_density, &flux_density); 364 if (RES_OK != res) { res = RES_BAD_ARG; goto error; } 365 if (flux_density < 0) { res = RES_BAD_ARG; goto error; } 366 source->flux_density.flux_density = flux_density; 367 368 /* Parse flux density unit */ 369 flux_density_unit = tokens[1]; 370 if (NULL == flux_density_unit){ res = RES_BAD_ARG; goto error; } 371 372 res = parse_flux_density_unit(source, txtrdr, flux_density_unit); 373 if (RES_OK != res) { goto error; } 374 375 filename = tokens[2]; 376 if (NULL == filename){ res = RES_BAD_ARG; goto error; } 377 378 spectral_unit = tokens[3]; 379 if (NULL == spectral_unit){ res = RES_BAD_ARG; goto error; } 380 381 property_unit = tokens[4]; 382 if (NULL == property_unit){ res = RES_BAD_ARG; goto error; } 383 384 res = parse_spectral_property 385 (source->sphin, filename, &source->flux_density.emission_spectrum); 386 if (RES_OK != res) { goto error; } 387 388 /* Wavelengths will be in nm after conversion */ 389 res = convert_wavelengths 390 (&source->flux_density.emission_spectrum->wavelengths, spectral_unit); 391 if (RES_OK != res) { goto error; } 392 393 /* Normalize spectrum */ 394 res = normalize_spectrum(source->flux_density.emission_spectrum); 395 if (RES_OK != res) { goto error; } 396 397 res = txtrdr_read_line(txtrdr); 398 if (RES_OK != res) { goto error; } 399 400 exit: 401 return res; 402 error: 403 goto exit; 404 } 405 /******************************************************************************* 406 * Local functions 407 ******************************************************************************/ 408 res_T 409 parse_source_surface 410 (struct sphin* sphin, 411 struct txtrdr* txtrdr, 412 const char* name, 413 struct sphin_source_surface** out_source) 414 { 415 struct sphin_source_surface* source = NULL; 416 struct str line; 417 char* keyword = NULL; 418 char* token = NULL; 419 char* token_ptr = NULL; 420 char* value = NULL; 421 res_T res = RES_OK; 422 423 ASSERT(NULL != sphin); 424 ASSERT(NULL != out_source); 425 ASSERT(NULL != txtrdr); 426 ASSERT(NULL != name); 427 ASSERT('\0' != name[0]); /* Name can't be empty */ 428 429 str_init(sphin->allocator, &line); 430 431 res = source_surface_create(sphin, name, &source); 432 if (RES_OK != res) { goto error; } 433 434 res = txtrdr_read_line(txtrdr); 435 if (RES_OK != res) { goto error; } 436 437 while (NULL != txtrdr_get_line(txtrdr)) { 438 res = str_set(&line, txtrdr_get_cline(txtrdr)); 439 if (RES_OK != res) { goto error; } 440 441 /* parse keyword */ 442 token = strtok_r(str_get(&line), ":", &token_ptr); 443 if (NULL == token){ res = RES_BAD_ARG; goto error; } 444 keyword = trim_keyword(token); 445 if (NULL == keyword){ res = RES_BAD_ARG; goto error; } 446 447 /* parse value */ 448 token = strtok_r(NULL, "", &token_ptr); 449 if (NULL == token){ res = RES_BAD_ARG; goto error; } 450 value = token; 451 if (0 == strcmp(keyword, "flux_density")){ 452 res = parse_flux_density(source, txtrdr, value); 453 } 454 else if (0 == strcmp(keyword, "direction")){ 455 res = parse_direction_distribution(source, txtrdr, value); 456 } 457 else { 458 break; 459 } 460 if (RES_OK != res) { goto error; } 461 } 462 463 exit: 464 str_release(&line); 465 *out_source = source; 466 return res; 467 error: 468 if (NULL != source) { 469 SPHIN(source_surface_ref_put(source)); 470 source = NULL; 471 } 472 goto exit; 473 } 474 475 /******************************************************************************* 476 * Exported functions 477 ******************************************************************************/ 478 res_T 479 sphin_source_surface_ref_get 480 (struct sphin_source_surface* source) 481 { 482 if (NULL == source) { 483 return RES_BAD_ARG; 484 } 485 ref_get(&source->ref); 486 return RES_OK; 487 } 488 489 res_T 490 sphin_source_surface_ref_put 491 (struct sphin_source_surface* source) 492 { 493 if (NULL == source) { 494 return RES_BAD_ARG; 495 } 496 ref_put(&source->ref, release_source_surface); 497 return RES_OK; 498 } 499 500 res_T 501 sphin_source_surface_get_direction_distribution 502 (const struct sphin_source_surface* source, 503 struct sphin_source_surface_direction_distribution* distrib) 504 { 505 if (NULL == source || NULL == distrib) { 506 return RES_BAD_ARG; 507 } 508 509 *distrib = source->direction_distribution; 510 return RES_OK; 511 } 512 513 res_T 514 sphin_source_surface_get_flux_density 515 (const struct sphin_source_surface* source, 516 struct sphin_source_surface_flux_density* density) 517 { 518 if (NULL == source || NULL == density) { 519 return RES_BAD_ARG; 520 } 521 522 *density = source->flux_density; 523 return RES_OK; 524 }