sln_tree.c (24470B)
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 #include "sln.h" 20 #include "sln_device_c.h" 21 #include "sln_line.h" 22 #include "sln_tree_c.h" 23 24 #include <star/shtr.h> 25 26 #include <rsys/algorithm.h> 27 #include <rsys/cstr.h> 28 #include <rsys/math.h> 29 30 #include <omp.h> 31 32 struct stream { 33 const char* name; 34 FILE* fp; 35 int intern_fp; /* Define if the stream was internally opened */ 36 }; 37 static const struct stream STREAM_NULL = {NULL, NULL, 0}; 38 39 /******************************************************************************* 40 * Helper functions 41 ******************************************************************************/ 42 /* Check that the sum of the molecular concentrations is equal to 1 */ 43 static res_T 44 check_molecule_concentration 45 (struct sln_device* sln, 46 const char* caller, 47 const struct sln_tree_create_args* args) 48 { 49 int i = 0; 50 double sum = 0; 51 ASSERT(sln && caller && args); 52 53 FOR_EACH(i, 0, SHTR_MAX_MOLECULE_COUNT) { 54 if(i == SHTR_MOLECULE_ID_NULL) continue; 55 sum += args->molecules[i].concentration; 56 } 57 58 /* The sum of molecular concentrations must be less than or equal to 1. It may 59 * be less than 1 if the remaining part of the mixture is (implicitly) defined 60 * as a radiatively inactive gas */ 61 if(sum > 1 && sum-1 > 1e-6) { 62 ERROR(sln, 63 "%s: the sum of molecule concentrations is greater than 1: %g\n", 64 caller, sum); 65 return RES_BAD_ARG; 66 } 67 68 return RES_OK; 69 } 70 71 /* Verify that the isotope abundance are valids */ 72 static res_T 73 check_molecule_isotope_abundances 74 (struct sln_device* sln, 75 const char* caller, 76 const struct sln_molecule* molecule) 77 { 78 int i = 0; 79 double sum = 0; 80 ASSERT(sln && caller && molecule); 81 82 /* The isotopic abundances are the default ones. Nothing to do */ 83 if(!molecule->non_default_isotope_abundances) return RES_OK; 84 85 /* The isotopic abundances are not the default ones. 86 * Verify that they are valid ... */ 87 FOR_EACH(i, 0, SHTR_MAX_ISOTOPE_COUNT) { 88 if(molecule->isotopes[i].abundance < 0) { 89 const int isotope_id = i + 1; /* isotope id in [1, 10] */ 90 ERROR(sln, "%s: invalid abundance of isotopie %d of %s: %g.\n", 91 caller, isotope_id, shtr_molecule_cstr(i), 92 molecule->isotopes[i].abundance); 93 return RES_BAD_ARG; 94 } 95 96 sum += molecule->isotopes[i].abundance; 97 } 98 99 /* ... and that their sum equals 1 */ 100 if(!eq_eps(sum, 1, 1e-6)) { 101 ERROR(sln, "%s: the %s isotope abundances does not sum to 1: %g.\n", 102 caller, shtr_molecule_cstr(i), sum); 103 return RES_BAD_ARG; 104 } 105 106 return RES_OK; 107 } 108 109 static res_T 110 check_molecules 111 (struct sln_device* sln, 112 const char* caller, 113 const struct sln_tree_create_args* args) 114 { 115 char molecule_ok[SHTR_MAX_MOLECULE_COUNT] = {0}; 116 117 double concentrations_sum = 0; 118 size_t iline = 0; 119 size_t nlines = 0; 120 res_T res = RES_OK; 121 ASSERT(args->lines); 122 123 res = check_molecule_concentration(sln, caller, args); 124 if(res != RES_OK) return res; 125 126 /* Iterate over the lines to define which molecules has to be checked, i.e., 127 * the ones used in the mixture */ 128 SHTR(line_list_get_size(args->lines, &nlines)); 129 FOR_EACH(iline, 0, nlines) { 130 struct shtr_line line = SHTR_LINE_NULL; 131 const struct sln_molecule* molecule = NULL; 132 133 SHTR(line_list_at(args->lines, iline, &line)); 134 135 /* This molecule was already checked */ 136 if(molecule_ok[line.molecule_id]) continue; 137 138 molecule = args->molecules + line.molecule_id; 139 140 if(molecule->concentration == 0) { 141 /* A molecular concentration of zero is allowed, but may be a user error, 142 * as 0 is the default concentration in the tree creation arguments. 143 * Therefore, warn the user about this value so that they can determine 144 * whether or not it is an error on their part. */ 145 WARN(sln, "%s: the concentration of %s is zero.\n", 146 caller, shtr_molecule_cstr(line.molecule_id)); 147 148 } else if(molecule->concentration < 0) { 149 /* Concentration cannot be negative... */ 150 ERROR(sln, "%s: invalid %s concentration: %g.\n", 151 FUNC_NAME, shtr_molecule_cstr(line.molecule_id), 152 molecule->concentration); 153 return RES_BAD_ARG; 154 } 155 156 concentrations_sum += molecule->concentration; 157 158 if(molecule->cutoff <= 0) { 159 /* ... cutoff either */ 160 ERROR(sln, "%s: invalid %s cutoff: %g.\n", 161 caller, shtr_molecule_cstr(line.molecule_id), molecule->cutoff); 162 return RES_BAD_ARG; 163 } 164 165 res = check_molecule_isotope_abundances(sln, caller, molecule); 166 if(res != RES_OK) return res; 167 168 molecule_ok[line.molecule_id] = 1; 169 } 170 171 /* The sum of molecular concentrations must be less than or equal to 1. It may 172 * be less than 1 if the remaining part of the mixture is (implicitly) defined 173 * as a radiatively inactive gas */ 174 if(concentrations_sum > 1 && (concentrations_sum - 1) > 1e-6) { 175 ERROR(sln, 176 "%s: the sum of molecule concentrations is greater than 1: %g\n", 177 caller, concentrations_sum); 178 return RES_BAD_ARG; 179 } 180 181 return RES_OK; 182 } 183 184 static res_T 185 check_sln_tree_create_args 186 (struct sln_device* sln, 187 const char* caller, 188 const struct sln_tree_create_args* args) 189 { 190 res_T res = RES_OK; 191 ASSERT(sln && caller); 192 193 if(!args) return RES_BAD_ARG; 194 195 if(!args->metadata) { 196 ERROR(sln, "%s: the isotope metadata are missing.\n", caller); 197 return RES_BAD_ARG; 198 } 199 200 if(!args->lines) { 201 ERROR(sln, "%s: the list of lines is missing.\n", caller); 202 return RES_BAD_ARG; 203 } 204 205 if(args->nvertices_hint == 0) { 206 ERROR(sln, 207 "%s: invalid hint on the number of vertices around the line center %lu.\n", 208 caller, (unsigned long)args->nvertices_hint); 209 return RES_BAD_ARG; 210 } 211 212 if(args->mesh_decimation_err < 0) { 213 ERROR(sln, "%s: invalid decimation error %g.\n", 214 caller, args->mesh_decimation_err); 215 return RES_BAD_ARG; 216 } 217 218 if((unsigned)args->mesh_type >= SLN_MESH_TYPES_COUNT__) { 219 ERROR(sln, "%s: invalid mesh type %d.\n", caller, args->mesh_type); 220 return RES_BAD_ARG; 221 } 222 223 if((unsigned)args->line_profile >= SLN_LINE_PROFILES_COUNT__) { 224 ERROR(sln, "%s: invalid line profile %d.\n", caller, args->line_profile); 225 return RES_BAD_ARG; 226 } 227 228 if(args->arity < 2 || args->arity > SLN_TREE_ARITY_MAX) { 229 ERROR(sln, "%s: invalid arity %u. It must be in [2, %d]\n", 230 caller, args->arity, SLN_TREE_ARITY_MAX); 231 return RES_BAD_ARG; 232 } 233 234 if(args->leaf_nlines < 1 || args->leaf_nlines > SLN_LEAF_NLINES_MAX) { 235 ERROR(sln, "%s: invalid number of lines per leaf %u. It must be in [1, %d]\n", 236 caller, args->leaf_nlines, SLN_LEAF_NLINES_MAX); 237 return RES_BAD_ARG; 238 } 239 240 if(args->nthreads_hint == 0) { 241 ERROR(sln, "%s: invalid number of threads %u\n", 242 caller, args->nthreads_hint); 243 return RES_BAD_ARG; 244 } 245 246 res = check_molecules(sln, caller, args); 247 if(res != RES_OK) return res; 248 249 return RES_OK; 250 } 251 252 static res_T 253 check_sln_tree_read_args 254 (struct sln_device* sln, 255 const char* caller, 256 const struct sln_tree_read_args* args) 257 { 258 if(!args) return RES_BAD_ARG; 259 260 if(!args->metadata) { 261 ERROR(sln, "%s: the isotope metadata are missing.\n", caller); 262 return RES_BAD_ARG; 263 } 264 265 if(!args->lines) { 266 ERROR(sln, "%s: the list of lines is missing.\n", caller); 267 return RES_BAD_ARG; 268 } 269 270 if(!args->file && !args->filename) { 271 ERROR(sln, 272 "%s: the source file is missing. No file name or stream is provided.\n", 273 caller); 274 return RES_BAD_ARG; 275 } 276 277 return RES_OK; 278 } 279 280 static res_T 281 check_sln_tree_write_args 282 (struct sln_device* sln, 283 const char* caller, 284 const struct sln_tree_write_args* args) 285 { 286 if(!args) return RES_BAD_ARG; 287 288 if(!args->file && !args->filename) { 289 ERROR(sln, 290 "%s: the destination file is missing. " 291 "No file name or stream is provided.\n", 292 caller); 293 return RES_BAD_ARG; 294 } 295 296 return RES_OK; 297 } 298 static INLINE void 299 stream_release(struct stream* stream) 300 { 301 ASSERT(stream); 302 if(stream->intern_fp && stream->fp) CHK(fclose(stream->fp) == 0); 303 } 304 305 static res_T 306 stream_init 307 (struct sln_device* sln, 308 const char* caller, 309 const char* name, /* NULL <=> default stream name */ 310 FILE* fp, /* NULL <=> open file "name" */ 311 const char* mode, /* mode in fopen */ 312 struct stream* stream) 313 { 314 res_T res = RES_OK; 315 316 ASSERT(sln && caller && stream); 317 ASSERT(fp || (name && mode)); 318 319 *stream = STREAM_NULL; 320 321 if(fp) { 322 stream->intern_fp = 0; 323 stream->name = name ? name : "stream"; 324 stream->fp = fp; 325 326 } else { 327 stream->intern_fp = 1; 328 stream->name = name; 329 if(!(stream->fp = fopen(name, mode))) { 330 ERROR(sln, "%s:%s: error opening file -- %s\n", 331 caller, name, strerror(errno)); 332 res = RES_IO_ERR; 333 goto error; 334 } 335 } 336 337 exit: 338 return res; 339 error: 340 if(stream->intern_fp && stream->fp) CHK(fclose(stream->fp) == 0); 341 goto exit; 342 } 343 344 static res_T 345 create_tree 346 (struct sln_device* sln, 347 const char* caller, 348 struct sln_tree** out_tree) 349 { 350 struct sln_tree* tree = NULL; 351 res_T res = RES_OK; 352 ASSERT(sln && caller && out_tree); 353 354 tree = MEM_CALLOC(sln->allocator, 1, sizeof(struct sln_tree)); 355 if(!tree) { 356 ERROR(sln, "%s: could not allocate the tree data structure.\n", 357 caller); 358 res = RES_MEM_ERR; 359 goto error; 360 } 361 ref_init(&tree->ref); 362 SLN(device_ref_get(sln)); 363 tree->sln = sln; 364 darray_node_init(sln->allocator, &tree->nodes); 365 darray_vertex_init(sln->allocator, &tree->vertices); 366 367 exit: 368 *out_tree = tree; 369 return res; 370 error: 371 if(tree) { SLN(tree_ref_put(tree)); tree = NULL; } 372 goto exit; 373 } 374 375 static INLINE int 376 cmp_nu_vtx(const void* key, const void* item) 377 { 378 const float nu = *((const float*)key); 379 const struct sln_vertex* vtx = item; 380 if(nu < vtx->wavenumber) return -1; 381 if(nu > vtx->wavenumber) return +1; 382 return 0; 383 } 384 385 static void 386 release_tree(ref_T* ref) 387 { 388 struct sln_tree* tree = CONTAINER_OF(ref, struct sln_tree, ref); 389 struct sln_device* sln = NULL; 390 ASSERT(ref); 391 sln = tree->sln; 392 darray_node_release(&tree->nodes); 393 darray_vertex_release(&tree->vertices); 394 if(tree->args.lines) SHTR(line_list_ref_put(tree->args.lines)); 395 if(tree->args.metadata) SHTR(isotope_metadata_ref_put(tree->args.metadata)); 396 MEM_RM(sln->allocator, tree); 397 SLN(device_ref_put(sln)); 398 } 399 400 /******************************************************************************* 401 * Local function 402 ******************************************************************************/ 403 unsigned 404 node_child_count(const struct sln_node* node, const unsigned tree_arity) 405 { 406 size_t nlines = 0; /* #lines in the node */ 407 size_t nlines_per_child = 0; /* Max #lines in a child */ 408 size_t nchildren = 0; 409 410 /* Pre-conditions */ 411 ASSERT(node && tree_arity >= 2); 412 413 /* Retrieve the node data and compute the #lines it partitions */ 414 nlines = node->range[1] - node->range[0] + 1; 415 ASSERT(nlines); 416 417 /* Based on the arity of the tree, calculate how the lines of the node are 418 * distributed among its children. For low lines count, i.e. when the minimum 419 * number of lines par child is less than the tree arity, the policy below 420 * prioritizes an equal distribution of lines among the children over 421 * maintaining the tree's arity. Thus, if a smaller number of children 422 * results in a more equitable distribution, this option is preferred over 423 * ensuring a number of children equal to the tree's arity. In other words, 424 * the tree's balance is prioritized. */ 425 nlines_per_child = (nlines + tree_arity-1/*ceil*/)/tree_arity; 426 427 /* From the previous line repartition, compute the number of children */ 428 nchildren = (nlines + nlines_per_child-1/*ceil*/)/nlines_per_child; 429 ASSERT(nchildren >= 2); 430 431 ASSERT(nchildren <= UINT_MAX); 432 return (unsigned)nchildren; 433 } 434 435 /******************************************************************************* 436 * Exported symbols 437 ******************************************************************************/ 438 res_T 439 sln_tree_create 440 (struct sln_device* device, 441 const struct sln_tree_create_args* args, 442 struct sln_tree** out_tree) 443 { 444 struct sln_tree* tree = NULL; 445 unsigned nthreads_max = 0; 446 res_T res = RES_OK; 447 448 if(!device || !out_tree) { res = RES_BAD_ARG; goto error; } 449 res = check_sln_tree_create_args(device, FUNC_NAME, args); 450 if(res != RES_OK) goto error; 451 452 res = create_tree(device, FUNC_NAME, &tree); 453 if(res != RES_OK) goto error; 454 SHTR(line_list_ref_get(args->lines)); 455 SHTR(isotope_metadata_ref_get(args->metadata)); 456 tree->args = *args; 457 458 /* Set the #threads to match the maximum number of available threads */ 459 nthreads_max = (unsigned)MMAX(omp_get_max_threads(), omp_get_num_procs()); 460 tree->args.nthreads_hint = MMIN(tree->args.nthreads_hint, nthreads_max); 461 462 res = tree_build(tree); 463 if(res != RES_OK) goto error; 464 465 exit: 466 if(out_tree) *out_tree = tree; 467 return res; 468 error: 469 if(tree) { SLN(tree_ref_put(tree)); tree = NULL; } 470 goto exit; 471 } 472 473 res_T 474 sln_tree_read 475 (struct sln_device* sln, 476 const struct sln_tree_read_args* args, 477 struct sln_tree** out_tree) 478 { 479 hash256_T hash_mdata1; 480 hash256_T hash_mdata2; 481 hash256_T hash_lines1; 482 hash256_T hash_lines2; 483 484 struct stream stream = STREAM_NULL; 485 struct sln_tree* tree = NULL; 486 size_t n = 0; 487 int version = 0; 488 res_T res = RES_OK; 489 490 if(!sln || !out_tree) { res = RES_BAD_ARG; goto error; } 491 res = check_sln_tree_read_args(sln, FUNC_NAME, args); 492 if(res != RES_OK) goto error; 493 494 res = create_tree(sln, FUNC_NAME, &tree); 495 if(res != RES_OK) goto error; 496 497 res = stream_init(sln, FUNC_NAME, args->filename, args->file, "r", &stream); 498 if(res != RES_OK) goto error; 499 500 #define READ(Var, Nb) { \ 501 if(fread((Var), sizeof(*(Var)), (Nb), stream.fp) != (Nb)) { \ 502 if(feof(stream.fp)) { \ 503 res = RES_BAD_ARG; \ 504 } else if(ferror(stream.fp)) { \ 505 res = RES_IO_ERR; \ 506 } else { \ 507 res = RES_UNKNOWN_ERR; \ 508 } \ 509 ERROR(sln, "%s: error loading the tree structure -- %s.\n", \ 510 stream.name, res_to_cstr(res)); \ 511 goto error; \ 512 } \ 513 } (void)0 514 READ(&version, 1); 515 if(version != SLN_TREE_VERSION) { 516 ERROR(sln, 517 "%s: unexpected tree version %d. Expecting a tree in version %d.\n", 518 stream.name, version, SLN_TREE_VERSION); 519 res = RES_BAD_ARG; 520 goto error; 521 } 522 523 res = shtr_isotope_metadata_hash(args->metadata, hash_mdata1); 524 if(res != RES_OK) goto error; 525 526 READ(hash_mdata2, sizeof(hash256_T)); 527 if(!hash256_eq(hash_mdata1, hash_mdata2)) { 528 ERROR(sln, 529 "%s: the input isotopic metadata are not those used " 530 "during tree construction.\n", stream.name); 531 res = RES_BAD_ARG; 532 goto error; 533 } 534 535 SHTR(isotope_metadata_ref_get(args->metadata)); 536 tree->args.metadata = args->metadata; 537 538 READ(hash_lines1, sizeof(hash256_T)); 539 if(!args->disable_line_hash_check) { 540 res = shtr_line_list_hash(args->lines, hash_lines2); 541 if(res != RES_OK) goto error; 542 543 if(!hash256_eq(hash_lines1, hash_lines2)) { 544 ERROR(sln, 545 "%s: the input list of lines is not the one used to build the tree.\n", 546 stream.name); 547 res = RES_BAD_ARG; 548 goto error; 549 } 550 } 551 552 SHTR(line_list_ref_get(args->lines)); 553 tree->args.lines = args->lines; 554 555 READ(&n, 1); 556 if((res = darray_node_resize(&tree->nodes, n)) != RES_OK) goto error; 557 READ(darray_node_data_get(&tree->nodes), n); 558 559 READ(&n, 1); 560 if((res = darray_vertex_resize(&tree->vertices, n)) != RES_OK) goto error; 561 READ(darray_vertex_data_get(&tree->vertices), n); 562 563 READ(&tree->args.line_profile, 1); 564 READ(&tree->args.molecules, 1); 565 READ(&tree->args.pressure, 1); 566 READ(&tree->args.temperature, 1); 567 READ(&tree->args.nvertices_hint, 1); 568 READ(&tree->args.mesh_decimation_err, 1); 569 READ(&tree->args.mesh_type, 1); 570 READ(&tree->args.arity, 1); 571 READ(&tree->args.leaf_nlines, 1); 572 #undef READ 573 574 exit: 575 stream_release(&stream); 576 if(out_tree) *out_tree = tree; 577 return res; 578 error: 579 if(tree) { SLN(tree_ref_put(tree)); tree = NULL; } 580 goto exit; 581 } 582 583 res_T 584 sln_tree_ref_get(struct sln_tree* tree) 585 { 586 if(!tree) return RES_BAD_ARG; 587 ref_get(&tree->ref); 588 return RES_OK; 589 } 590 591 res_T 592 sln_tree_ref_put(struct sln_tree* tree) 593 { 594 if(!tree) return RES_BAD_ARG; 595 ref_put(&tree->ref, release_tree); 596 return RES_OK; 597 } 598 599 res_T 600 sln_tree_get_desc(const struct sln_tree* tree, struct sln_tree_desc* desc) 601 { 602 const struct sln_node* node = NULL; 603 unsigned depth = 0; 604 605 if(!tree || !desc) return RES_BAD_ARG; 606 607 desc->mesh_decimation_err = tree->args.mesh_decimation_err; 608 desc->mesh_type = tree->args.mesh_type; 609 desc->line_profile = tree->args.line_profile; 610 desc->nnodes = darray_node_size_get(&tree->nodes); 611 desc->nvertices = darray_vertex_size_get(&tree->vertices); 612 desc->pressure = tree->args.pressure; /* [atm] */ 613 desc->temperature = tree->args.temperature; /* [K] */ 614 desc->arity = tree->args.arity; 615 desc->leaf_nlines = tree->args.leaf_nlines; 616 617 SHTR(line_list_get_size(tree->args.lines, &desc->nlines)); 618 619 node = sln_tree_get_root(tree); 620 while(!sln_node_is_leaf(node)) { 621 node = sln_node_get_child(tree, node, 0); 622 ++depth; 623 } 624 desc->depth = depth; 625 626 return RES_OK; 627 } 628 629 const struct sln_node* 630 sln_tree_get_root(const struct sln_tree* tree) 631 { 632 ASSERT(tree); 633 if(darray_node_size_get(&tree->nodes)) { 634 return darray_node_cdata_get(&tree->nodes); 635 } else { 636 return NULL; 637 } 638 } 639 640 res_T 641 sln_tree_get_line 642 (const struct sln_tree* tree, 643 const size_t iline, 644 struct sln_line* line) 645 { 646 size_t nlines = 0; 647 res_T res = RES_OK; 648 649 if(!tree || !line) { res = RES_BAD_ARG; goto error; } 650 651 SHTR(line_list_get_size(tree->args.lines, &nlines)); 652 if(iline >= nlines) { res = RES_BAD_ARG; goto error; } 653 654 res = line_setup(tree, iline, line); 655 if(res != RES_OK) { 656 ERROR(tree->sln, "%s: could not setup the line %lu-- %s\n", 657 FUNC_NAME, iline, res_to_cstr(res)); 658 goto error; 659 } 660 661 exit: 662 return res; 663 error: 664 goto exit; 665 } 666 667 int 668 sln_node_is_leaf(const struct sln_node* node) 669 { 670 ASSERT(node); 671 return node->offset == 0; 672 } 673 674 unsigned 675 sln_node_get_child_count 676 (const struct sln_tree* tree, 677 const struct sln_node* node) 678 { 679 ASSERT(tree && node); 680 681 if(sln_node_is_leaf(node)) { 682 return 0; /* No child */ 683 } else { 684 return node_child_count(node, tree->args.arity); 685 } 686 } 687 688 const struct sln_node* 689 sln_node_get_child 690 (const struct sln_tree* tree, 691 const struct sln_node* node, 692 const unsigned ichild) 693 { 694 ASSERT(node && ichild < sln_node_get_child_count(tree, node)); 695 ASSERT(!sln_node_is_leaf(node)); 696 (void)tree; 697 return node + node->offset + ichild; 698 } 699 700 res_T 701 sln_node_get_mesh 702 (const struct sln_tree* tree, 703 const struct sln_node* node, 704 struct sln_mesh* mesh) 705 { 706 if(!tree || !node || !mesh) return RES_BAD_ARG; 707 mesh->vertices = darray_vertex_cdata_get(&tree->vertices) + node->ivertex; 708 mesh->nvertices = node->nvertices; 709 return RES_OK; 710 } 711 712 double 713 sln_node_eval 714 (const struct sln_tree* tree, 715 const struct sln_node* node, 716 const double nu) 717 { 718 double ka = 0; 719 size_t iline; 720 ASSERT(tree && node); 721 722 FOR_EACH(iline, node->range[0], node->range[1]+1) { 723 struct sln_line line = SLN_LINE_NULL; 724 res_T res = RES_OK; 725 726 res = line_setup(tree, iline, &line); 727 if(res != RES_OK) { 728 WARN(tree->sln, "%s: could not setup the line %lu-- %s\n", 729 FUNC_NAME, iline, res_to_cstr(res)); 730 continue; 731 } 732 733 ka += sln_line_eval(tree, &line, nu); 734 } 735 return ka; 736 } 737 738 res_T 739 sln_node_get_desc 740 (const struct sln_tree* tree, 741 const struct sln_node* node, 742 struct sln_node_desc* desc) 743 { 744 if(!tree || !node || !desc) return RES_BAD_ARG; 745 desc->ilines[0] = node->range[0]; 746 desc->ilines[1] = node->range[1]; 747 desc->nvertices = node->nvertices; 748 desc->nchildren = sln_node_get_child_count(tree, node); 749 return RES_OK; 750 } 751 752 double 753 sln_mesh_eval(const struct sln_mesh* mesh, const double wavenumber) 754 { 755 const struct sln_vertex* vtx0 = NULL; 756 const struct sln_vertex* vtx1 = NULL; 757 const float nu = (float)wavenumber; 758 size_t n; /* #vertices */ 759 double u; /* Linear interpolation parameter */ 760 ASSERT(mesh && mesh->nvertices); 761 762 n = mesh->nvertices; 763 764 /* Handle special cases */ 765 if(n == 1) return mesh->vertices[0].ka; 766 if(nu < mesh->vertices[0].wavenumber 767 || nu > mesh->vertices[n-1].wavenumber) { 768 return 0; 769 } 770 if(nu == mesh->vertices[0].wavenumber) return mesh->vertices[0].ka; 771 if(nu == mesh->vertices[n-1].wavenumber) return mesh->vertices[n-1].ka; 772 773 /* Dichotomic search of the mesh vertex whose wavenumber is greater than or 774 * equal to the submitted wavenumber 'nu' */ 775 vtx1 = search_lower_bound(&nu, mesh->vertices, n, sizeof(*vtx1), cmp_nu_vtx); 776 vtx0 = vtx1 - 1; 777 ASSERT(vtx1); /* A vertex is necessary found ...*/ 778 ASSERT(vtx1 > mesh->vertices); /* ... and it cannot be the first one */ 779 ASSERT(vtx0->wavenumber < nu && nu <= vtx1->wavenumber); 780 781 /* Compute the linear interpolation parameter */ 782 u = (wavenumber - vtx0->wavenumber) / (vtx1->wavenumber - vtx0->wavenumber); 783 u = CLAMP(u, 0, 1); /* Handle numerical imprecisions */ 784 785 if(u == 0) return vtx0->ka; 786 if(u == 1) return vtx1->ka; 787 return u*(vtx1->ka - vtx0->ka) + vtx0->ka; 788 } 789 790 res_T 791 sln_tree_write 792 (const struct sln_tree* tree, 793 const struct sln_tree_write_args* args) 794 { 795 struct stream stream = STREAM_NULL; 796 size_t nnodes, nverts; 797 hash256_T hash_mdata; 798 hash256_T hash_lines; 799 res_T res = RES_OK; 800 801 if(!tree) { res = RES_BAD_ARG; goto error; } 802 res = check_sln_tree_write_args(tree->sln, FUNC_NAME, args); 803 if(res != RES_OK) goto error; 804 805 res = shtr_isotope_metadata_hash(tree->args.metadata, hash_mdata); 806 if(res != RES_OK) goto error; 807 res = shtr_line_list_hash(tree->args.lines, hash_lines); 808 if(res != RES_OK) goto error; 809 810 res = stream_init 811 (tree->sln, FUNC_NAME, args->filename, args->file, "w", &stream); 812 if(res != RES_OK) goto error; 813 814 #define WRITE(Var, Nb) { \ 815 if(fwrite((Var), sizeof(*(Var)), (Nb), stream.fp) != (Nb)) { \ 816 ERROR(tree->sln, "%s:%s: error writing the tree -- %s\n", \ 817 FUNC_NAME, stream.name, strerror(errno)); \ 818 res = RES_IO_ERR; \ 819 goto error; \ 820 } \ 821 } (void)0 822 WRITE(&SLN_TREE_VERSION, 1); 823 WRITE(hash_mdata, sizeof(hash256_T)); 824 WRITE(hash_lines, sizeof(hash256_T)); 825 826 nnodes = darray_node_size_get(&tree->nodes); 827 WRITE(&nnodes, 1); 828 WRITE(darray_node_cdata_get(&tree->nodes), nnodes); 829 830 nverts = darray_vertex_size_get(&tree->vertices); 831 WRITE(&nverts, 1); 832 WRITE(darray_vertex_cdata_get(&tree->vertices), nverts); 833 834 WRITE(&tree->args.line_profile, 1); 835 WRITE(&tree->args.molecules, 1); 836 WRITE(&tree->args.pressure, 1); 837 WRITE(&tree->args.temperature, 1); 838 WRITE(&tree->args.nvertices_hint, 1); 839 WRITE(&tree->args.mesh_decimation_err, 1); 840 WRITE(&tree->args.mesh_type, 1); 841 WRITE(&tree->args.arity, 1); 842 WRITE(&tree->args.leaf_nlines, 1); 843 #undef WRITE 844 845 exit: 846 stream_release(&stream); 847 return res; 848 error: 849 goto exit; 850 }