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