sln_mixture.c (14505B)
1 /* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2026 Université de Lorraine 3 * Copyright (C) 2022 Centre National de la Recherche Scientifique 4 * Copyright (C) 2022 Université Paul Sabatier 5 * 6 * This program is free software: you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation, either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 18 19 #define _POSIX_C_SOURCE 200809L /* strtok_r support */ 20 21 #include "sln.h" 22 #include "sln_device_c.h" 23 24 #include <rsys/cstr.h> 25 #include <rsys/text_reader.h> 26 27 #include <ctype.h> /* isalpha support */ 28 #include <string.h> 29 30 struct molecule { 31 struct sln_molecule param; 32 int nisotopes; 33 enum shtr_molecule_id id; 34 }; 35 static const struct molecule MOLECULE_NULL = { 36 SLN_MOLECULE_NULL, 0, SHTR_MOLECULE_ID_NULL 37 }; 38 39 #define MOLECULE_IS_VALID(Molecule) ((Molecule)->id > SHTR_MOLECULE_ID_NULL) 40 41 struct sln_mixture { 42 struct molecule molecules[SHTR_MAX_MOLECULE_COUNT]; 43 int nmolecules; 44 struct sln_device* sln; 45 ref_T ref; 46 }; 47 48 /******************************************************************************* 49 * Helper functions 50 ******************************************************************************/ 51 static INLINE res_T 52 check_sln_mixture_load_args(const struct sln_mixture_load_args* args) 53 { 54 if(!args || !args->filename) return RES_BAD_ARG; 55 if(!args->molparam) return RES_BAD_ARG; 56 return RES_OK; 57 } 58 59 static res_T 60 check_molecules_cocentration 61 (const struct sln_mixture* mixture, 62 const char* mixture_name) 63 { 64 double sum = 0; 65 int i = 0; 66 res_T res = RES_OK; 67 ASSERT(mixture); 68 69 FOR_EACH(i, 0, mixture->nmolecules) { 70 sum += mixture->molecules[i].param.concentration; 71 } 72 73 /* The sum of molecular concentrations must be less than or equal to 1. It may 74 * be less than 1 if the remaining part of the mixture is (implicitly) defined 75 * as a radiatively inactive gas */ 76 if(sum > 1 && sum-1 > 1e-6) { 77 ERROR(mixture->sln, 78 "%s: the sum of molecule concentrations is greater than 1: %g\n", 79 mixture_name, sum); 80 res = RES_BAD_ARG; 81 goto error; 82 } 83 84 exit: 85 return res; 86 error: 87 goto exit; 88 } 89 90 static INLINE res_T 91 flush_molecule 92 (struct sln_mixture* mixture, 93 struct molecule* molecule, 94 const char* mixture_name) 95 { 96 res_T res = RES_OK; 97 ASSERT(mixture && molecule); 98 99 mixture->molecules[mixture->nmolecules++] = *molecule; 100 101 /* Check isotope abundances */ 102 if(molecule->param.non_default_isotope_abundances) { 103 double sum = 0; 104 int i = 0; 105 106 FOR_EACH(i, 0, molecule->nisotopes) { 107 sum += molecule->param.isotopes[i].abundance; 108 } 109 110 if(!eq_eps(sum, 1, 1e-6)) { 111 ERROR(mixture->sln, 112 "%s: the sum of the isotopic abundances of %s is %g, " 113 "whereas it should be equal to 1\n", 114 mixture_name, shtr_molecule_cstr(molecule->id), sum); 115 res = RES_BAD_ARG; 116 goto error; 117 } 118 } 119 120 *molecule = MOLECULE_NULL; 121 122 exit: 123 return res; 124 error: 125 goto exit; 126 } 127 128 static INLINE res_T 129 parse_molecule_name 130 (struct sln_mixture* mixture, 131 struct molecule* molecule, 132 const char* name) 133 { 134 size_t i = 0; 135 ASSERT(mixture && molecule && name); 136 (void)mixture; 137 138 /* Search for the molecule ID. Use a simple linear search, as there are only a 139 * few molecules. */ 140 FOR_EACH(i, 1, SHTR_MAX_MOLECULE_COUNT) { 141 if(!strcmp(name, shtr_molecule_cstr(i))) break; 142 } 143 if(i >= SHTR_MAX_MOLECULE_COUNT) return RES_BAD_ARG; 144 molecule->id = i; 145 146 return RES_OK; 147 } 148 149 static INLINE int 150 molecule_has_metadata 151 (const enum shtr_molecule_id id, 152 struct shtr_isotope_metadata* molparam) 153 { 154 struct shtr_molecule shtr_molecule = SHTR_MOLECULE_NULL; 155 ASSERT(molparam); 156 157 SHTR(isotope_metadata_find_molecule(molparam, id, &shtr_molecule)); 158 return !SHTR_MOLECULE_IS_NULL(&shtr_molecule); 159 } 160 161 static INLINE int 162 molecule_already_defined 163 (struct sln_mixture* mixture, 164 const enum shtr_molecule_id id) 165 { 166 int i = 0; 167 ASSERT(mixture); 168 169 /* Check that this molecule has not already been parsed. Use a simple linear 170 * search for the molecule, as there are only a few of them. */ 171 FOR_EACH(i, 0, mixture->nmolecules) { 172 if(mixture->molecules[i].id == id) return 1; 173 } 174 return 0; 175 } 176 177 static res_T 178 parse_molecule 179 (struct sln_mixture* mixture, 180 const struct sln_mixture_load_args* args, 181 struct molecule* molecule, 182 struct txtrdr* txtrdr) 183 { 184 char* line = NULL; 185 char* tk = NULL; 186 char* tk_ctx = NULL; 187 res_T res = RES_OK; 188 189 ASSERT(mixture && args && molecule && txtrdr); 190 (void)args; 191 192 line = txtrdr_get_line(txtrdr); 193 ASSERT(line); 194 195 #define LOG(Type, Str, ...) \ 196 Type(mixture->sln, "%s:%lu: "Str, \ 197 txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr), \ 198 __VA_ARGS__) 199 200 tk = strtok_r(line, " \t", &tk_ctx); 201 res = parse_molecule_name(mixture, molecule, tk); 202 if(res != RES_OK) { 203 LOG(ERROR, "invalid molecule name `%s'\n", tk); 204 goto error; 205 } 206 if(molecule_already_defined(mixture, molecule->id)) { 207 LOG(ERROR, "duplicate molecule `%s'\n", tk); 208 res = RES_BAD_ARG; 209 goto error; 210 } 211 if(!molecule_has_metadata(molecule->id, args->molparam)) { 212 LOG(ERROR, "`%s' does not have isotope metadata\n", tk); 213 res = RES_BAD_ARG; 214 goto error; 215 } 216 217 tk = strtok_r(NULL, " \t", &tk_ctx); 218 res = cstr_to_double(tk, &molecule->param.concentration); 219 if(res == RES_OK) { 220 if(molecule->param.concentration < 0 221 || molecule->param.concentration > 1) { 222 res = RES_BAD_ARG; 223 } 224 } 225 if(res != RES_OK) { 226 LOG(ERROR, "invalid concentration `%s'\n", tk ? tk : "(null)"); 227 goto error; 228 } 229 230 tk = strtok_r(NULL, " \t", &tk_ctx); 231 res = cstr_to_double(tk, &molecule->param.cutoff); 232 if(res == RES_OK && molecule->param.cutoff <= 0) res = RES_BAD_ARG; 233 if(res != RES_OK) { 234 LOG(ERROR, "invalid cutoff `%s'\n", tk ? tk : "(null)"); 235 goto error; 236 } 237 238 tk = strtok_r(NULL, "", &tk_ctx); 239 if(tk) LOG(WARN, "unexpected text `%s'\n", tk); 240 241 #undef LOG 242 243 exit: 244 return res; 245 error: 246 goto exit; 247 } 248 249 /* Initialize the isotope identifiers of molecules based on those defined in the 250 * isotope metadata */ 251 static res_T 252 setup_molecule_isotope_ids 253 (struct shtr_isotope_metadata* molparam, 254 struct molecule* molecule) 255 { 256 struct shtr_molecule shtr_molecule = SHTR_MOLECULE_NULL; 257 size_t i = 0; 258 res_T res = RES_OK; 259 260 ASSERT(molparam && molecule); 261 262 /* The molecule must be defined in the isotope metadata */ 263 res = shtr_isotope_metadata_find_molecule 264 (molparam, molecule->id, &shtr_molecule); 265 if(res == RES_OK && SHTR_MOLECULE_IS_NULL(&shtr_molecule)) res = RES_BAD_ARG; 266 if(res != RES_OK) goto error; 267 268 FOR_EACH(i, 0, shtr_molecule.nisotopes) { 269 molecule->param.isotopes[i].id = shtr_molecule.isotopes[i].id; 270 } 271 272 exit: 273 return res; 274 error: 275 goto exit; 276 } 277 278 static res_T 279 parse_isotope 280 (struct sln_mixture* mixture, 281 const struct sln_mixture_load_args* args, 282 struct molecule* molecule, 283 struct txtrdr* txtrdr) 284 { 285 char* line = NULL; 286 char* tk = NULL; 287 char* tk_ctx = NULL; 288 int iiso = 0; 289 int iso_id = 0; 290 res_T res = RES_OK; 291 ASSERT(mixture && args && molecule && txtrdr); 292 293 line = txtrdr_get_line(txtrdr); 294 ASSERT(line); 295 296 if(!molecule->param.non_default_isotope_abundances) { 297 /* No isotopes should have been parsed */ 298 ASSERT(molecule->nisotopes == 0); 299 300 res = setup_molecule_isotope_ids(args->molparam, molecule); 301 if(res != RES_OK) goto error; 302 303 molecule->param.non_default_isotope_abundances = 1; 304 } 305 306 iiso = molecule->nisotopes; /* Local index of the isotope */ 307 308 #define LOG(Type, Str, ...) \ 309 Type(mixture->sln, "%s:%lu: "Str, \ 310 txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr), \ 311 __VA_ARGS__) 312 313 tk = strtok_r(line, " \t", &tk_ctx); 314 res = cstr_to_int(tk, &iso_id); 315 if(res != RES_OK) { 316 LOG(ERROR, "invalid isotope index `%s'\n", tk ? tk : "(null)"); 317 goto error; 318 } 319 if(iso_id != molecule->param.isotopes[iiso].id) { 320 LOG(ERROR, 321 "expecting isotope %d of %s, whereas the actual isotope is %d\n", 322 molecule->param.isotopes[iiso].id, shtr_molecule_cstr(molecule->id), 323 iso_id); 324 res = RES_BAD_ARG; 325 goto error; 326 } 327 328 tk = strtok_r(NULL, " \t", &tk_ctx); 329 res = cstr_to_double(tk, &molecule->param.isotopes[iiso].abundance); 330 if(res != RES_OK || molecule->param.isotopes[iiso].abundance < 0) { 331 LOG(ERROR, "invalid isotope abundance `%s'\n", tk ? tk : "(null)"); 332 res = RES_BAD_ARG; 333 goto error; 334 } 335 336 tk = strtok_r(NULL, "", &tk_ctx); 337 if(tk) LOG(WARN, "unexpected text `%s'\n", tk); 338 339 #undef LOG 340 341 molecule->param.non_default_isotope_abundances = 1; 342 molecule->nisotopes += 1; 343 344 exit: 345 return res; 346 error: 347 goto exit; 348 } 349 350 static res_T 351 parse_line 352 (struct sln_mixture* mixture, 353 const struct sln_mixture_load_args* args, 354 struct molecule* molecule, /* Currently parsed molecule */ 355 struct txtrdr* txtrdr) 356 { 357 const char* line = NULL; 358 size_t i; 359 res_T res = RES_OK; 360 ASSERT(mixture && args && molecule && txtrdr); 361 362 line = txtrdr_get_cline(txtrdr); 363 ASSERT(line); 364 i = strspn(line, " \t"); 365 ASSERT(i < strlen(line)); 366 367 /* New molecule */ 368 if(isalpha(line[i])) { 369 370 /* Register the previous molecule if any */ 371 if(MOLECULE_IS_VALID(molecule)) { 372 res = flush_molecule(mixture, molecule, txtrdr_get_name(txtrdr)); 373 if(res != RES_OK) goto error; 374 } 375 376 /* Parse the molecule data, i.e., name, concentration and cutoff */ 377 res = parse_molecule(mixture, args, molecule, txtrdr); 378 if(res != RES_OK) goto error; 379 380 /* If there is no molecule being analyzed, no isotopic data can be parsed */ 381 } else if(!MOLECULE_IS_VALID(molecule)) { 382 ERROR(mixture->sln, "%s:%lu: missing a molecule\n", 383 txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr)); 384 res = RES_BAD_ARG; 385 goto error; 386 387 /* Parse the isotopes for the currently parsed molecule */ 388 } else { 389 res = parse_isotope(mixture, args, molecule, txtrdr); 390 if(res != RES_OK) goto error; 391 } 392 393 exit: 394 return res; 395 error: 396 goto exit; 397 } 398 399 static res_T 400 load_stream 401 (struct sln_mixture* mixture, 402 FILE* fp, 403 const struct sln_mixture_load_args* args) 404 { 405 struct molecule molecule = MOLECULE_NULL; 406 struct txtrdr* txtrdr = NULL; 407 res_T res = RES_OK; 408 ASSERT(mixture && fp && args); 409 410 res = txtrdr_stream(mixture->sln->allocator, fp, args->filename, '#', &txtrdr); 411 if(res != RES_OK) goto error; 412 413 #define READ_LINE if((res = txtrdr_read_line(txtrdr)) != RES_OK) goto error 414 415 for(;;) { 416 READ_LINE; 417 if(!txtrdr_get_cline(txtrdr)) break; /* No more parsed line */ 418 419 res = parse_line(mixture, args, &molecule, txtrdr); 420 if(res != RES_OK) goto error; 421 } 422 #undef READ_LINE 423 424 if(MOLECULE_IS_VALID(&molecule)) { 425 res = flush_molecule(mixture, &molecule, txtrdr_get_name(txtrdr)); 426 if(res != RES_OK) goto error; 427 } 428 exit: 429 if(txtrdr) txtrdr_ref_put(txtrdr); 430 return res; 431 error: 432 if(!txtrdr) { 433 ERROR(mixture->sln, "%s: loading error -- %s\n", 434 args->filename, res_to_cstr(res)); 435 } else { 436 ERROR(mixture->sln, "%s:%lu: loading error -- %s\n", 437 args->filename, (unsigned long)txtrdr_get_line_num(txtrdr), 438 res_to_cstr(res)); 439 } 440 goto exit; 441 } 442 443 static res_T 444 load_mixture 445 (struct sln_mixture* mixture, 446 const struct sln_mixture_load_args* args) 447 { 448 FILE* fp = NULL; 449 res_T res = RES_OK; 450 ASSERT(mixture && args); 451 452 if(args->file) { /* Load from stream */ 453 fp = args->file; 454 455 } else { /* Load from file */ 456 fp = fopen(args->filename, "r"); 457 if(!fp) { 458 ERROR(mixture->sln, "error opening file `%s' -- %s\n", 459 args->filename, strerror(errno)); 460 res = RES_IO_ERR; 461 goto error; 462 } 463 } 464 465 res = load_stream(mixture, fp, args); 466 if(res != RES_OK) goto error; 467 468 res = check_molecules_cocentration(mixture, args->filename); 469 if(res != RES_OK) goto error; 470 471 exit: 472 if(fp && fp != args->file) fclose(fp); 473 return res; 474 error: 475 goto exit; 476 } 477 478 static void 479 release_mixture(ref_T* ref) 480 { 481 struct sln_mixture* mixture = CONTAINER_OF(ref, struct sln_mixture, ref); 482 struct sln_device* sln = NULL; 483 ASSERT(ref); 484 sln = mixture->sln; 485 MEM_RM(sln->allocator, mixture); 486 SLN(device_ref_put(sln)); 487 } 488 489 /******************************************************************************* 490 * Exported functions 491 ******************************************************************************/ 492 res_T 493 sln_mixture_load 494 (struct sln_device* sln, 495 const struct sln_mixture_load_args* args, 496 struct sln_mixture** out_mixture) 497 { 498 struct sln_mixture* mixture = NULL; 499 res_T res = RES_OK; 500 501 if(!sln || !out_mixture) { res = RES_BAD_ARG; goto error; } 502 res = check_sln_mixture_load_args(args); 503 if(res != RES_OK) goto error; 504 505 mixture = MEM_CALLOC(sln->allocator, 1, sizeof(*mixture)); 506 if(!mixture) { 507 ERROR(sln, "Could not allocate the mixture data structure.\n"); 508 res = RES_MEM_ERR; 509 goto error; 510 } 511 ref_init(&mixture->ref); 512 SLN(device_ref_get(sln)); 513 mixture->sln = sln; 514 515 res = load_mixture(mixture, args); 516 if(res != RES_OK) goto error; 517 518 exit: 519 if(out_mixture) *out_mixture = mixture; 520 return res; 521 error: 522 if(mixture) { SLN(mixture_ref_put(mixture)); mixture = NULL; } 523 goto exit; 524 } 525 526 res_T 527 sln_mixture_ref_get(struct sln_mixture* mixture) 528 { 529 if(!mixture) return RES_BAD_ARG; 530 ref_get(&mixture->ref); 531 return RES_OK; 532 } 533 534 res_T 535 sln_mixture_ref_put(struct sln_mixture* mixture) 536 { 537 if(!mixture) return RES_BAD_ARG; 538 ref_put(&mixture->ref, release_mixture); 539 return RES_OK; 540 } 541 542 int 543 sln_mixture_get_molecule_count(const struct sln_mixture* mixture) 544 { 545 ASSERT(mixture); 546 return mixture->nmolecules; 547 } 548 549 enum shtr_molecule_id 550 sln_mixture_get_molecule_id(const struct sln_mixture* mixture, const int index) 551 { 552 ASSERT(mixture && 0 <= index && index < mixture->nmolecules); 553 return mixture->molecules[index].id; 554 } 555 556 res_T 557 sln_mixture_get_molecule 558 (const struct sln_mixture* mixture, 559 const int index, 560 struct sln_molecule* molecule) 561 { 562 if(!mixture || index < 0 || index >= mixture->nmolecules || !molecule) { 563 return RES_BAD_ARG; 564 } 565 566 *molecule = mixture->molecules[index].param; 567 return RES_OK; 568 }