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