schiff_args.c (51638B)
1 /* Copyright (C) 2015, 2016, 2026 Centre National de la Recherche Scientifique 2 * Copyright (C) 2026 Clermont Auvergne INP 3 * Copyright (C) 2026 Institut Mines Télécom Albi-Carmaux 4 * Copyright (C) 2017, 2019-2021, 2026 |Méso|Star> (contact@meso-star.com) 5 * Copyright (C) 2026 Université de Lorraine 6 * Copyright (C) 2026 Université de Toulouse 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 2 /* getopt support */ 22 23 #include "schiff_args.h" 24 #include "schiff_version.h" 25 #include "schiff_optical_properties.h" 26 #include "schiff_geometry.h" 27 28 #include <rsys/dynamic_array_char.h> 29 #include <rsys/cstr.h> 30 #include <rsys/str.h> 31 #include <rsys/stretchy_array.h> 32 33 #include <stdarg.h> 34 #include <yaml.h> 35 36 #define MIN_NANGLES 2 37 #define MIN_NANGLES_INV 2 38 39 #include <unistd.h> 40 41 /******************************************************************************* 42 * Helper function 43 ******************************************************************************/ 44 static void 45 usage(FILE* stream) 46 { 47 fprintf(stream, 48 "schiff [-Dhqv] [-a nangles_phase_func] [-A nangles_phase_func_invcumul]\n" 49 " [-d ninner_samples] [-g nrealisations] [-G nparticles]\n" 50 " [-i distribution] [-l particles_length] [-n nthreads] [-o output]\n" 51 " [-w wavelelength[:wavelelength ...]] [optical_props]\n"); 52 } 53 54 static INLINE void 55 print_version(void) 56 { 57 printf("Schiff %d.%d.%d\n", 58 SCHIFF_VERSION_MAJOR, 59 SCHIFF_VERSION_MINOR, 60 SCHIFF_VERSION_PATCH); 61 } 62 63 static int 64 cmp_double(const void* op0, const void* op1) 65 { 66 const double* a = op0; 67 const double* b = op1; 68 return *a < *b ? -1 : (*a > *b ? 1 : 0); 69 } 70 71 static INLINE void 72 param_release(struct schiff_param* param) 73 { 74 ASSERT(param); 75 76 if(param->distribution == SCHIFF_PARAM_HISTOGRAM 77 && param->data.histogram.entries) { 78 sa_release(param->data.histogram.entries); 79 } 80 param->distribution = SCHIFF_PARAM_NONE; 81 } 82 83 static void 84 geometry_release(struct schiff_geometry* geom) 85 { 86 int i; 87 ASSERT(geom); 88 switch(geom->type) { 89 case SCHIFF_ELLIPSOID: 90 case SCHIFF_ELLIPSOID_AS_SPHERE: 91 param_release(&geom->data.ellipsoid.a); 92 param_release(&geom->data.ellipsoid.c); 93 param_release(&geom->data.ellipsoid.radius_sphere); 94 break; 95 case SCHIFF_CYLINDER_AS_SPHERE: 96 case SCHIFF_CYLINDER: 97 param_release(&geom->data.cylinder.radius); 98 param_release(&geom->data.cylinder.height); 99 param_release(&geom->data.cylinder.radius_sphere); 100 break; 101 case SCHIFF_HELICAL_PIPE: 102 case SCHIFF_HELICAL_PIPE_AS_SPHERE: 103 param_release(&geom->data.helical_pipe.pitch); 104 param_release(&geom->data.helical_pipe.height); 105 param_release(&geom->data.helical_pipe.radius_helicoid); 106 param_release(&geom->data.helical_pipe.radius_circle); 107 param_release(&geom->data.helical_pipe.radius_sphere); 108 break; 109 case SCHIFF_SPHERE: 110 param_release(&geom->data.sphere.radius); 111 break; 112 case SCHIFF_SUPERSHAPE: 113 case SCHIFF_SUPERSHAPE_AS_SPHERE: 114 FOR_EACH(i, 0, 6) { 115 param_release(&geom->data.supershape.formulas[0][i]); 116 param_release(&geom->data.supershape.formulas[1][i]); 117 } 118 param_release(&geom->data.supershape.radius_sphere); 119 break; 120 case SCHIFF_NONE: /* Do nothing */ break; 121 default: FATAL("Unreachable code\n"); break; 122 } 123 geom->type = SCHIFF_NONE; 124 } 125 126 static res_T 127 parse_wavelengths(const char* str, struct schiff_args* args) 128 { 129 size_t len; 130 size_t i; 131 res_T res = RES_OK; 132 ASSERT(args && str); 133 134 /* How many wavelengths are submitted */ 135 res = cstr_to_list_double(str, ':', NULL, &len, 0); 136 if(res != RES_OK) goto error; 137 138 /* Reserve the wavelengths memory space */ 139 sa_clear(args->wavelengths); 140 args->wavelengths = sa_add(args->wavelengths, len); 141 142 /* Read the wavelengths */ 143 res = cstr_to_list_double(optarg, ':', args->wavelengths, NULL, len); 144 if(res != RES_OK) goto error; 145 146 /* Check the validity of read wavelengths */ 147 FOR_EACH(i, 0, len) { 148 if(args->wavelengths[i] < 0.0) { 149 res = RES_BAD_ARG; 150 goto error; 151 } 152 } 153 exit: 154 return res; 155 error: 156 goto exit; 157 } 158 159 static INLINE void 160 log_err 161 (const char* filename, 162 const yaml_node_t* node, 163 const char* fmt, 164 ...) 165 { 166 va_list vargs_list; 167 ASSERT(node && fmt); 168 169 fprintf(stderr, "%s:%lu:%lu: ", 170 filename, 171 (unsigned long)node->start_mark.line+1, 172 (unsigned long)node->start_mark.column+1); 173 va_start(vargs_list, fmt); 174 vfprintf(stderr, fmt, vargs_list); 175 va_end(vargs_list); 176 } 177 178 static res_T 179 parse_yaml_double 180 (const char* filename, 181 const yaml_node_t* node, 182 const double min_val, /* Minimum valid value */ 183 const double max_val, /* Maximum valid value */ 184 double* value) 185 { 186 res_T res = RES_OK; 187 ASSERT(filename && node && value && min_val < max_val); 188 189 if(node->type != YAML_SCALAR_NODE) { 190 log_err(filename, node, "expecting a floating point number.\n"); 191 return RES_BAD_ARG; 192 } 193 res = cstr_to_double((char*)node->data.scalar.value, value); 194 if(res != RES_OK) { 195 log_err(filename, node, "invalid floating point number `%s'.\n", 196 node->data.scalar.value); 197 return RES_BAD_ARG; 198 } 199 if(*value < min_val || *value > max_val) { 200 log_err(filename, node, 201 "the floating point number %g is not in [%g, %g].\n", 202 *value, min_val, max_val); 203 return RES_BAD_ARG; 204 } 205 return RES_OK; 206 } 207 208 static res_T 209 parse_yaml_uint 210 (const char* filename, 211 const yaml_node_t* node, 212 const unsigned min_val, /* Minimum valid value */ 213 const unsigned max_val, /* Maximum valid value */ 214 unsigned* value) 215 { 216 res_T res = RES_OK; 217 ASSERT(filename && node && value); 218 ASSERT(min_val < max_val); 219 220 if(node->type != YAML_SCALAR_NODE) { 221 log_err(filename, node, "expecting an unsigned integer.\n"); 222 return RES_BAD_ARG; 223 } 224 res = cstr_to_uint((char*)node->data.scalar.value, value); 225 if(res != RES_OK) { 226 log_err(filename, node, "invalid unsigned integer `%s'.\n", 227 node->data.scalar.value); 228 return RES_BAD_ARG; 229 } 230 if(*value < min_val || *value > max_val) { 231 log_err(filename, node, 232 "the unsigned integer %u is not in [%u, %u].\n", 233 *value, min_val, max_val); 234 return RES_BAD_ARG; 235 } 236 237 return RES_OK; 238 } 239 240 static res_T 241 parse_yaml_param_mu_sigma 242 (const char* filename, 243 yaml_document_t* doc, 244 const yaml_node_t* distrib, 245 const char* distrib_name, 246 const double min_val, /* Minimum valid value fo the parameter */ 247 const double max_val, /* Maximum valid value of the parameter */ 248 double* mu, 249 double* sigma) 250 { 251 enum { 252 MU = BIT(0), 253 SIGMA = BIT(1) 254 }; 255 size_t nattrs; 256 size_t i; 257 int mask = 0; /* Register the parsed attributes */ 258 res_T res = RES_OK; 259 ASSERT(filename && doc && distrib && distrib_name); 260 ASSERT(min_val < max_val && mu && sigma); 261 262 if(distrib->type != YAML_MAPPING_NODE) { 263 log_err(filename, distrib, 264 "expecting a mapping of the %s attributes.\n", distrib_name); 265 return RES_BAD_ARG; 266 } 267 268 /* Parse the log normal attributes */ 269 nattrs = (size_t) 270 (distrib->data.mapping.pairs.top - distrib->data.mapping.pairs.start); 271 FOR_EACH(i, 0, nattrs) { 272 yaml_node_t* key, *val; 273 274 key=yaml_document_get_node(doc, distrib->data.mapping.pairs.start[i].key); 275 val=yaml_document_get_node(doc, distrib->data.mapping.pairs.start[i].value); 276 ASSERT(key->type == YAML_SCALAR_NODE); 277 278 #define SETUP_MASK(Flag, Name) { \ 279 if(mask & Flag) { \ 280 log_err(filename, key, "the "Name" %s attribute is already defined.\n",\ 281 distrib_name); \ 282 return RES_BAD_ARG; \ 283 } \ 284 mask |= Flag; \ 285 } (void)0 286 287 /* mu attribute */ 288 if(!strcmp((char*)key->data.scalar.value, "mu")) { 289 SETUP_MASK(MU, "`mu'"); 290 res = parse_yaml_double(filename, val, min_val, max_val, mu); 291 292 /* sigma attribute */ 293 } else if(!strcmp((char*)key->data.scalar.value, "sigma")) { 294 SETUP_MASK(SIGMA, "`sigma'"); 295 res = parse_yaml_double(filename, val, DBL_MIN, DBL_MAX, sigma); 296 297 /* unknown attribute */ 298 } else { 299 log_err(filename, key, "unknown %s attribute `%s'.\n", 300 distrib_name, key->data.scalar.value); 301 return RES_BAD_ARG; 302 } 303 if(res != RES_OK) return res; 304 305 #undef SETUP_MASK 306 } 307 308 /* Ensure that the attributes are all parsed */ 309 if(!(mask & MU)) { 310 log_err(filename, distrib, "missing the `mu' %s attribute.\n", 311 distrib_name); 312 return RES_BAD_ARG; 313 } else if(!(mask & SIGMA)) { 314 log_err(filename, distrib, "missing the `sigma' %s attribute.\n", 315 distrib_name); 316 return RES_BAD_ARG; 317 } 318 319 return RES_OK; 320 } 321 322 static res_T 323 parse_yaml_param_histogram 324 (const char* filename, 325 yaml_document_t* doc, 326 const yaml_node_t* histo, 327 const double min_val, /* Minimum valid value fo the parameter */ 328 const double max_val, /* Maximum valid value of the parameter */ 329 struct schiff_param* param) 330 { 331 enum { 332 LOWER = BIT(0), 333 UPPER = BIT(1), 334 PROBAS = BIT(2) 335 }; 336 size_t ientry, nentries, nattrs; 337 size_t i; 338 int mask = 0; /* Register the parsed histogram attributes */ 339 double accum_proba; 340 res_T res = RES_OK; 341 ASSERT(filename && doc && histo && param && min_val < max_val); 342 343 param->distribution = SCHIFF_PARAM_HISTOGRAM; 344 param->data.histogram.entries = NULL; 345 346 if(histo->type != YAML_MAPPING_NODE) { 347 log_err(filename, histo, "expecting a mapping of the histogram data.\n"); 348 res = RES_BAD_ARG; 349 goto error; 350 } 351 352 /* Parse the histogram data */ 353 nattrs = (size_t) 354 (histo->data.mapping.pairs.top - histo->data.mapping.pairs.start); 355 FOR_EACH(i, 0, nattrs) { 356 yaml_node_t *key, *val; 357 358 key = yaml_document_get_node(doc, histo->data.mapping.pairs.start[i].key); 359 val = yaml_document_get_node(doc, histo->data.mapping.pairs.start[i].value); 360 ASSERT(key->type == YAML_SCALAR_NODE); 361 362 #define SETUP_MASK(Flag, Name) { \ 363 if(mask & Flag) { \ 364 log_err(filename, key, "the "Name" is already defined.\n"); \ 365 return RES_BAD_ARG; \ 366 } \ 367 mask |= Flag; \ 368 } (void)0 369 370 /* Histogram lower bound */ 371 if(!strcmp((char*)key->data.scalar.value, "lower")) { 372 SETUP_MASK(LOWER, "histogram lower bound"); 373 res = parse_yaml_double 374 (filename, val, min_val, max_val, ¶m->data.histogram.lower); 375 if(res != RES_OK) goto error; 376 377 /* Histogram upper bound */ 378 } else if(!strcmp((char*)key->data.scalar.value, "upper")) { 379 SETUP_MASK(UPPER, "histogram upper bound"); 380 res = parse_yaml_double 381 (filename, val, min_val, max_val, ¶m->data.histogram.upper); 382 if(res != RES_OK) goto error; 383 384 /* Histogram entries */ 385 } else if(!strcmp((char*)key->data.scalar.value, "probabilities")) { 386 SETUP_MASK(PROBAS, "histogram data"); 387 388 if(val->type != YAML_SEQUENCE_NODE) { 389 log_err(filename, val, 390 "expecting a sequence of floating point numbers.\n"); 391 res = RES_BAD_ARG; 392 goto error; 393 } 394 395 nentries = (size_t) 396 (val->data.sequence.items.top - val->data.sequence.items.start); 397 if(!sa_add(param->data.histogram.entries, nentries)) { 398 log_err(filename, val, 399 "couldn't allocate an histogram with %lu entries.\n", 400 (unsigned long)nentries); 401 res = RES_MEM_ERR; 402 goto error; 403 } 404 405 /* Parse histogram entries */ 406 accum_proba = 0; 407 FOR_EACH(ientry, 0, nentries) { 408 double proba; 409 yaml_node_t* node; 410 node = yaml_document_get_node(doc, val->data.sequence.items.start[ientry]); 411 412 res = parse_yaml_double(filename, node, DBL_MIN, DBL_MAX, &proba); 413 if(res != RES_OK) goto error; 414 415 accum_proba += proba; 416 param->data.histogram.entries[ientry] = accum_proba; 417 } 418 419 /* Normalize the histogram entries */ 420 FOR_EACH(ientry, 0, nentries) 421 param->data.histogram.entries[ientry] /= accum_proba; 422 423 } else { 424 log_err(filename, key, "unknown histogram data `%s'.\n", 425 key->data.scalar.value); 426 res = RES_BAD_ARG; 427 goto error; 428 } 429 430 #undef SETUP_MASK 431 } 432 433 /* Ensure that the histogram data are all parsed. */ 434 if(!(mask & LOWER)) { 435 log_err(filename, histo, "missing the histogram lower parameter.\n"); 436 res = RES_BAD_ARG; 437 goto error; 438 } else if(!(mask & UPPER)) { 439 log_err(filename, histo, "missing the histogram upper parameter.\n"); 440 res = RES_BAD_ARG; 441 goto error; 442 } else if(!(mask & PROBAS)) { 443 log_err(filename, histo, "missing the histogram probabilities.\n"); 444 res = RES_BAD_ARG; 445 goto error; 446 } 447 448 /* Check the histogram interval */ 449 if(param->data.histogram.upper <= param->data.histogram.lower) { 450 log_err(filename, histo, "invalid histogram interval [%g, %g].\n", 451 param->data.histogram.lower, param->data.histogram.upper); 452 res = RES_BAD_ARG; 453 goto error; 454 } 455 456 exit: 457 return res; 458 error: 459 if(param->data.histogram.entries) { 460 sa_release(param->data.histogram.entries); 461 param->data.histogram.entries = NULL; 462 } 463 goto exit; 464 } 465 466 static res_T 467 parse_yaml_param_distribution 468 (const char* filename, 469 yaml_document_t* doc, 470 const yaml_node_t* node, 471 const double min_val, /* Minimum valid value fo the parameter */ 472 const double max_val, /* Maximum valid value of the parameter */ 473 struct schiff_param* param) 474 { 475 res_T res = RES_OK; 476 ASSERT(filename && doc && node && param && min_val < max_val); 477 478 if(node->type == YAML_SCALAR_NODE) { /* Floating point constant */ 479 param->distribution = SCHIFF_PARAM_CONSTANT; 480 res = parse_yaml_double 481 (filename, node, min_val, max_val, ¶m->data.constant); 482 if(res != RES_OK) return res; 483 484 } else if(node->type == YAML_MAPPING_NODE) { 485 yaml_node_t* key, *val; 486 487 if(node->data.mapping.pairs.top - node->data.mapping.pairs.start != 1) { 488 log_err(filename, node, 489 "expecting a mapping from a parameter distribution to its attributes.\n"); 490 return RES_BAD_ARG; 491 } 492 493 key = yaml_document_get_node(doc, node->data.mapping.pairs.start[0].key); 494 val = yaml_document_get_node(doc, node->data.mapping.pairs.start[0].value); 495 ASSERT(key->type == YAML_SCALAR_NODE); 496 497 /* Lognormal distribution */ 498 if(!strcmp((char*)key->data.scalar.value, "lognormal")) { 499 res = parse_yaml_param_mu_sigma(filename, doc, val, "lognormal", min_val, 500 max_val, ¶m->data.lognormal.mu, ¶m->data.lognormal.sigma); 501 if(res != RES_OK) return res; 502 param->distribution = SCHIFF_PARAM_LOGNORMAL; 503 504 /* Gaussian distribution */ 505 } else if(!strcmp((char*)key->data.scalar.value, "gaussian")) { 506 res = parse_yaml_param_mu_sigma(filename, doc, val, "gaussian", min_val, 507 max_val, ¶m->data.gaussian.mu, ¶m->data.gaussian.sigma); 508 if(res != RES_OK) return res; 509 param->distribution = SCHIFF_PARAM_GAUSSIAN; 510 param->data.gaussian.range[0] = min_val; 511 param->data.gaussian.range[1] = max_val; 512 513 /* Histogram distribution */ 514 } else if(!strcmp((char*)key->data.scalar.value, "histogram")) { 515 res = parse_yaml_param_histogram 516 (filename, doc, val, min_val, max_val, param); 517 if(res != RES_OK) return res; 518 519 /* Unknown distribution */ 520 } else { 521 log_err(filename, key, "unknown parameter distribution `%s'.\n", 522 key->data.scalar.value); 523 return RES_BAD_ARG; 524 } 525 526 } else { 527 log_err(filename, node, "unexpected YAML node.\n"); 528 return RES_BAD_ARG; 529 } 530 531 return RES_OK; 532 } 533 534 static res_T 535 parse_yaml_superformula 536 (const char* filename, 537 yaml_document_t* doc, 538 const yaml_node_t* node, 539 struct schiff_param formula[6]) 540 { 541 int mask = 0; /* Register the parsed histogram attributes */ 542 size_t nattrs; 543 size_t i; 544 res_T res = RES_OK; 545 ASSERT(filename && doc && node && formula); 546 547 if(node->type != YAML_MAPPING_NODE) { 548 log_err(filename, node, 549 "expecting a mapping of superformula parameters.\n"); 550 return RES_BAD_ARG; 551 } 552 553 nattrs = (size_t) 554 (node->data.mapping.pairs.top - node->data.mapping.pairs.start); 555 556 FOR_EACH(i, 0, nattrs) { 557 yaml_node_t* key, *val; 558 559 key = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].key); 560 val = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].value); 561 ASSERT(key->type == YAML_SCALAR_NODE); 562 563 #define PARSE_SUPER_SHAPE_PARAM(Param) \ 564 if(!strcmp((char*)key->data.scalar.value, STR(Param))) { \ 565 if(mask & BIT(Param)) { \ 566 log_err(filename, key, \ 567 "the "STR(Param)" superformula parameter is already defined.\n"); \ 568 return RES_BAD_ARG; \ 569 } \ 570 mask |= BIT(Param); \ 571 res = parse_yaml_param_distribution \ 572 (filename, doc, val, DBL_MIN, DBL_MAX, formula + Param); \ 573 if(res != RES_OK) return res; \ 574 continue; \ 575 } (void)0 576 PARSE_SUPER_SHAPE_PARAM(A); 577 PARSE_SUPER_SHAPE_PARAM(B); 578 PARSE_SUPER_SHAPE_PARAM(M); 579 PARSE_SUPER_SHAPE_PARAM(N0); 580 PARSE_SUPER_SHAPE_PARAM(N1); 581 PARSE_SUPER_SHAPE_PARAM(N2); 582 #undef PARSE_SUPER_SHAPE_PARAM 583 584 log_err(filename, key, "unknown superformula parameter `%s'.\n", 585 key->data.scalar.value); 586 return RES_BAD_ARG; 587 } 588 #define CHECK_SUPER_SHAPE_PARAM(Param) \ 589 if(!(mask & BIT(Param))) { \ 590 log_err(filename, node, \ 591 "missing the "STR(Param)" superformula parameter.\n"); \ 592 return RES_BAD_ARG; \ 593 } (void)0 594 CHECK_SUPER_SHAPE_PARAM(A); 595 CHECK_SUPER_SHAPE_PARAM(B); 596 CHECK_SUPER_SHAPE_PARAM(M); 597 CHECK_SUPER_SHAPE_PARAM(N0); 598 CHECK_SUPER_SHAPE_PARAM(N1); 599 CHECK_SUPER_SHAPE_PARAM(N2); 600 #undef CHECK_SUPER_SHAPE_PARAM 601 return RES_OK; 602 } 603 604 static res_T 605 parse_yaml_ellipsoid 606 (const char* filename, 607 yaml_document_t* doc, 608 const yaml_node_t* node, 609 struct schiff_geometry* geom, 610 double* geom_proba) 611 { 612 enum { 613 PROBA = BIT(0), 614 LENGTH_a = BIT(1), 615 LENGTH_c = BIT(2), 616 RADIUS_SPHERE = BIT(3), 617 SLICES = BIT(4) 618 }; 619 int mask = 0; /* Register the parsed histogram attributes */ 620 size_t nattrs; 621 size_t i; 622 res_T res = RES_OK; 623 ASSERT(filename && doc && node && geom && geom_proba); 624 625 /* Setup the type at the beginning in order to define what arguments should 626 * be released if a parsing error occurs. Note that one can define the main 627 * type or its "equivalent sphere" variation. */ 628 geom->type = SCHIFF_ELLIPSOID; 629 630 if(node->type != YAML_MAPPING_NODE) { 631 log_err(filename, node, "expecting a mapping of ellipsoid attributes.\n"); 632 return RES_BAD_ARG; 633 } 634 635 nattrs = (size_t) 636 (node->data.mapping.pairs.top - node->data.mapping.pairs.start); 637 638 FOR_EACH(i, 0, nattrs) { 639 yaml_node_t* key; 640 yaml_node_t* val; 641 642 key = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].key); 643 val = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].value); 644 ASSERT(key->type == YAML_SCALAR_NODE); 645 646 #define SETUP_MASK(Flag, Name) { \ 647 if(mask & Flag) { \ 648 log_err(filename, key, "the "Name" is already defined.\n"); \ 649 return RES_BAD_ARG; \ 650 } \ 651 mask |= Flag; \ 652 } (void)0 653 654 /* Probability of the distribution */ 655 if(!strcmp((char*)key->data.scalar.value, "proba")) { 656 SETUP_MASK(PROBA, "sphere proba"); 657 res = parse_yaml_double(filename, val, DBL_MIN, DBL_MAX, geom_proba); 658 659 /* # slices used to discretized the triangular ellipsoid */ 660 } else if(!strcmp((char*)key->data.scalar.value, "slices")) { 661 SETUP_MASK(SLICES, "ellipsoid number of slices"); 662 res = parse_yaml_uint 663 (filename, val, 4, 32768, &geom->data.ellipsoid.nslices); 664 665 /* equivalent sphere radius */ 666 } else if(!strcmp((char*)key->data.scalar.value, "radius_sphere")) { 667 SETUP_MASK(RADIUS_SPHERE, "equivalent sphere radius of the ellipsoid"); 668 res = parse_yaml_param_distribution(filename, doc, val, DBL_MIN, DBL_MAX, 669 &geom->data.ellipsoid.radius_sphere); 670 671 /* Length of the ellipsoid "a" semi-principal axis */ 672 } else if(!strcmp((char*)key->data.scalar.value, "a")) { 673 SETUP_MASK(LENGTH_a, "ellipsoid \"a\" parameter"); 674 res = parse_yaml_param_distribution 675 (filename, doc, val, DBL_MIN, DBL_MAX, &geom->data.ellipsoid.a); 676 677 /* Length of the ellipsoid "c" semi-principal axis */ 678 } else if(!strcmp((char*)key->data.scalar.value, "c")) { 679 SETUP_MASK(LENGTH_c, "ellipsoid \"c\" parameter"); 680 res = parse_yaml_param_distribution 681 (filename, doc, val, DBL_MIN, DBL_MAX, &geom->data.ellipsoid.c); 682 683 /* Error */ 684 } else { 685 log_err(filename, key, "unknown sphere attribute `%s'.\n", 686 key->data.scalar.value); 687 return RES_BAD_ARG; 688 } 689 if(res != RES_OK) return res; 690 691 #undef SETUP_MASK 692 } 693 694 /* Ensure the completeness of the ellipsoid distribution */ 695 if(!(mask & LENGTH_a)) { 696 log_err(filename, node, 697 "missing the length of the semi principal axis \"a\".\n"); 698 return RES_BAD_ARG; 699 } else if(!(mask & LENGTH_c)) { 700 log_err(filename, node, 701 "missing the length of the semi principal axis \"c\".\n"); 702 return RES_BAD_ARG; 703 } 704 /* Setup the default values if required */ 705 if(!(mask & PROBA)) { /* Default proba */ 706 *geom_proba = 1.0; 707 } 708 if(!(mask & SLICES)) { /* Default number of slices */ 709 geom->data.ellipsoid.nslices = SCHIFF_ELLIPSOID_DEFAULT.nslices; 710 } 711 /* Define the geometry type */ 712 if(!(mask & RADIUS_SPHERE)) { 713 geom->type = SCHIFF_ELLIPSOID; 714 } else { 715 if(geom->data.ellipsoid.a.distribution != SCHIFF_PARAM_CONSTANT 716 || geom->data.ellipsoid.c.distribution != SCHIFF_PARAM_CONSTANT) { 717 log_err(filename, node, 718 "the radius_sphere parameter cannot be defined with non constant " 719 "ellipsoid attributes.\n"); 720 return RES_BAD_ARG; 721 } 722 geom->type = SCHIFF_ELLIPSOID_AS_SPHERE; 723 } 724 return RES_OK; 725 } 726 727 static res_T 728 parse_yaml_cylinder 729 (const char* filename, 730 yaml_document_t* doc, 731 const yaml_node_t* node, 732 struct schiff_geometry* geom, 733 double* geom_proba) 734 { 735 enum { 736 PROBA = BIT(0), 737 RADIUS = BIT(1), 738 HEIGHT = BIT(2), 739 RADIUS_SPHERE = BIT(3), 740 SLICES = BIT(4) 741 }; 742 int mask = 0; /* Register the parsed attributes */ 743 size_t nattrs; 744 size_t i; 745 res_T res = RES_OK; 746 ASSERT(filename && doc && node && geom && geom_proba); 747 748 /* Setup the type at the beginning in order to define what arguments should 749 * be released if a parsing error occurs. Note that one can define the main 750 * type or its "equivalent sphere" variation. */ 751 geom->type = SCHIFF_CYLINDER; 752 753 if(node->type != YAML_MAPPING_NODE) { 754 log_err(filename, node, "expecting a mapping of cylinder attributes.\n"); 755 return RES_BAD_ARG; 756 } 757 758 nattrs = (size_t) 759 (node->data.mapping.pairs.top - node->data.mapping.pairs.start); 760 761 FOR_EACH(i, 0, nattrs) { 762 yaml_node_t* key; 763 yaml_node_t* val; 764 765 key = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].key); 766 val = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].value); 767 ASSERT(key->type == YAML_SCALAR_NODE); 768 769 #define SETUP_MASK(Flag, Name) { \ 770 if(mask & Flag) { \ 771 log_err(filename, key, "the "Name" is already defined.\n"); \ 772 return RES_BAD_ARG; \ 773 } \ 774 mask |= Flag; \ 775 } (void)0 776 777 /* Distribution probability */ 778 if(!strcmp((char*)key->data.scalar.value, "proba")) { 779 SETUP_MASK(PROBA, "cylinder proba"); 780 res = parse_yaml_double(filename, val, DBL_MIN, DBL_MAX, geom_proba); 781 782 /* Cylinder radius */ 783 } else if(!strcmp((char*)key->data.scalar.value, "radius")) { 784 SETUP_MASK(RADIUS, "cylinder radius"); 785 res = parse_yaml_param_distribution 786 (filename, doc, val, DBL_MIN, DBL_MAX, &geom->data.cylinder.radius); 787 788 /* Cylinder height */ 789 } else if(!strcmp((char*)key->data.scalar.value, "height")) { 790 SETUP_MASK(HEIGHT, "cylinder height"); 791 res = parse_yaml_param_distribution 792 (filename, doc, val, DBL_MIN, DBL_MAX, &geom->data.cylinder.height); 793 794 /* Equivalent sphere radius */ 795 } else if(!strcmp((char*)key->data.scalar.value, "radius_sphere")) { 796 SETUP_MASK(RADIUS_SPHERE, "equivalent sphere radius of the cylinder"); 797 res = parse_yaml_param_distribution(filename, doc, val, DBL_MIN, DBL_MAX, 798 &geom->data.cylinder.radius_sphere); 799 800 /* # slices used to discretized the cylinder */ 801 } else if(!strcmp((char*)key->data.scalar.value, "slices")) { 802 SETUP_MASK(SLICES, "cylinder number of slices"); 803 res = parse_yaml_uint 804 (filename, val, 4, 32768, &geom->data.cylinder.nslices); 805 806 /* Error */ 807 } else { 808 log_err(filename, key, "unknown cylinder attribute `%s'.\n", 809 key->data.scalar.value); 810 res = RES_BAD_ARG; 811 } 812 if(res != RES_OK) return res; 813 814 #undef SETUP_MASK 815 } 816 817 /* Ensure the completness of the cylinder distribution */ 818 if(!(mask & RADIUS)) { 819 log_err(filename, node, "missing the radius attribute.\n"); 820 return RES_BAD_ARG; 821 } else if(!(mask & HEIGHT)) { 822 log_err(filename, node, "missing the height attribute.\n"); 823 return RES_BAD_ARG; 824 } 825 /* Setup the default values if required */ 826 if(!(mask & PROBA)) { /* Default proba */ 827 *geom_proba = 1.0; 828 } 829 if(!(mask & SLICES)) { /* Default number of slices */ 830 geom->data.cylinder.nslices = SCHIFF_CYLINDER_DEFAULT.nslices; 831 } 832 /* Define the geometry type */ 833 if(!(mask & RADIUS_SPHERE)) { 834 geom->type = SCHIFF_CYLINDER; 835 } else { 836 if(geom->data.cylinder.radius.distribution != SCHIFF_PARAM_CONSTANT 837 || geom->data.cylinder.height.distribution != SCHIFF_PARAM_CONSTANT) { 838 log_err(filename, node, 839 "the radius_sphere parameter cannot be defined with non constant " 840 "cylinder attributes.\n"); 841 return RES_BAD_ARG; 842 } 843 geom->type = SCHIFF_CYLINDER_AS_SPHERE; 844 } 845 return RES_OK; 846 } 847 848 static res_T 849 parse_yaml_helical_pipe 850 (const char* filename, 851 yaml_document_t* doc, 852 const yaml_node_t* node, 853 struct schiff_geometry* geom, 854 double* geom_proba) 855 { 856 enum { 857 PROBA = BIT(0), 858 PITCH = BIT(1), 859 HEIGHT = BIT(2), 860 RADIUS_HELICOID = BIT(3), 861 RADIUS_CIRCLE = BIT(4), 862 SLICES_HELICOID = BIT(5), 863 SLICES_CIRCLE = BIT(6), 864 RADIUS_SPHERE = BIT(7) 865 }; 866 int mask = 0; /* Register the parsed attributes */ 867 size_t nattrs; 868 size_t i; 869 res_T res = RES_OK; 870 ASSERT(filename && doc && node && geom && geom_proba); 871 872 /* Setup the type at the beginning in order to define what arguments should 873 * be released if a parsing error occurs. Note that one can define the main 874 * type or its "equivalent sphere" variation. */ 875 geom->type = SCHIFF_HELICAL_PIPE; 876 877 if(node->type != YAML_MAPPING_NODE) { 878 log_err(filename, node, "expecting a mapping of helical pipe attributes.\n"); 879 return RES_BAD_ARG; 880 } 881 882 nattrs = (size_t) 883 (node->data.mapping.pairs.top - node->data.mapping.pairs.start); 884 885 FOR_EACH(i, 0, nattrs) { 886 yaml_node_t* key; 887 yaml_node_t* val; 888 889 key = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].key); 890 val = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].value); 891 ASSERT(key->type == YAML_SCALAR_NODE); 892 893 #define SETUP_MASK(Flag, Name) { \ 894 if(mask & Flag) { \ 895 log_err(filename, key, "the "Name" is already defined.\n"); \ 896 return RES_BAD_ARG; \ 897 } \ 898 mask |= Flag; \ 899 } (void)0 900 901 /* Probability of the distribution */ 902 if(!strcmp((char*)key->data.scalar.value, "proba")) { 903 SETUP_MASK(PROBA, "helical pipe proba"); 904 res = parse_yaml_double(filename, val, DBL_MIN, DBL_MAX, geom_proba); 905 906 /* Helicoid pitch */ 907 } else if(!strcmp((char*)key->data.scalar.value, "pitch")) { 908 SETUP_MASK(PITCH, "helical pipe pitch"); 909 res = parse_yaml_param_distribution 910 (filename, doc, val, DBL_MIN, DBL_MAX, &geom->data.helical_pipe.pitch); 911 912 /* Helicoid height */ 913 } else if(!strcmp((char*)key->data.scalar.value, "height")) { 914 SETUP_MASK(HEIGHT, "helical pipe height"); 915 res = parse_yaml_param_distribution 916 (filename, doc, val, DBL_MIN, DBL_MAX, &geom->data.helical_pipe.height); 917 918 /* Radius of the helicoid */ 919 } else if(!strcmp((char*)key->data.scalar.value, "radius_helicoid")) { 920 SETUP_MASK(RADIUS_HELICOID, "helicoid radius"); 921 res = parse_yaml_param_distribution(filename, doc, val, DBL_MIN, DBL_MAX, 922 &geom->data.helical_pipe.radius_helicoid); 923 924 /* Radius of the meridian circle */ 925 } else if(!strcmp((char*)key->data.scalar.value, "radius_circle")) { 926 SETUP_MASK(RADIUS_CIRCLE, "circle radius of the helical pipe"); 927 res = parse_yaml_param_distribution(filename, doc, val, DBL_MIN, DBL_MAX, 928 &geom->data.helical_pipe.radius_circle); 929 930 /* # slices of the helicoid */ 931 } else if(!strcmp((char*)key->data.scalar.value, "slices_helicoid")) { 932 SETUP_MASK(SLICES_HELICOID, "helicoid number of slices"); 933 res = parse_yaml_uint 934 (filename, val, 4, 32768, &geom->data.helical_pipe.nslices_helicoid); 935 936 /* # slices of the circle */ 937 } else if(!strcmp((char*)key->data.scalar.value, "slices_circle")) { 938 SETUP_MASK(SLICES_CIRCLE, "helicoid meridian circle number of slices"); 939 res = parse_yaml_uint 940 (filename, val, 4, 32768, &geom->data.helical_pipe.nslices_circle); 941 942 /* Equivalent sphere radius */ 943 } else if(!strcmp((char*)key->data.scalar.value, "radius_sphere")) { 944 SETUP_MASK(RADIUS_SPHERE, "equivalent sphere radius of the helical pipe"); 945 res = parse_yaml_param_distribution(filename, doc, val, DBL_MIN, DBL_MAX, 946 &geom->data.helical_pipe.radius_sphere); 947 948 /* Error */ 949 } else { 950 log_err(filename, key, "unknown helical pipe attribute `%s'.\n", 951 key->data.scalar.value); 952 return RES_BAD_ARG; 953 } 954 if(res != RES_OK) return res; 955 #undef SETUP_MASK 956 } 957 958 /* Ensure the completeness of the helical pipe distribution */ 959 if(!(mask & PITCH)) { 960 log_err(filename, node, "missing the pitch of the helical pipe.\n"); 961 return RES_BAD_ARG; 962 } else if(!(mask & HEIGHT)) { 963 log_err(filename, node, "missing the height of the helical pipe.\n"); 964 return RES_BAD_ARG; 965 } else if(!(mask & RADIUS_HELICOID)) { 966 log_err(filename, node, "missing the radius of the helicoid.\n"); 967 return RES_BAD_ARG; 968 } else if(!(mask & RADIUS_CIRCLE)) { 969 log_err(filename, node, "missing the radius of the meridian circle.\n"); 970 return RES_BAD_ARG; 971 } 972 /* Setup the default values if required */ 973 if(!(mask & PROBA)) { 974 *geom_proba = 1.0; 975 } 976 if(!(mask & SLICES_HELICOID)) { 977 geom->data.helical_pipe.nslices_helicoid = 978 SCHIFF_HELICAL_PIPE_DEFAULT.nslices_helicoid; 979 } 980 if(!(mask & SLICES_CIRCLE)) { 981 geom->data.helical_pipe.nslices_circle = 982 SCHIFF_HELICAL_PIPE_DEFAULT.nslices_circle; 983 } 984 /* Define the geometry type */ 985 if(!(mask & RADIUS_SPHERE)) { 986 geom->type = SCHIFF_HELICAL_PIPE; 987 } else { 988 const struct schiff_helical_pipe* hpipe = &geom->data.helical_pipe; 989 if(hpipe->pitch.distribution != SCHIFF_PARAM_CONSTANT 990 || hpipe->height.distribution != SCHIFF_PARAM_CONSTANT 991 || hpipe->radius_helicoid.distribution != SCHIFF_PARAM_CONSTANT 992 || hpipe->radius_circle.distribution != SCHIFF_PARAM_CONSTANT) { 993 log_err(filename, node, 994 "the radius_sphere parameter cannot be defined with non constant " 995 "helical pipe attributes.\n"); 996 return RES_BAD_ARG; 997 } 998 geom->type = SCHIFF_HELICAL_PIPE_AS_SPHERE; 999 } 1000 return RES_OK; 1001 } 1002 1003 static res_T 1004 parse_yaml_sphere 1005 (const char* filename, 1006 yaml_document_t* doc, 1007 const yaml_node_t* node, 1008 struct schiff_geometry* geom, 1009 double* geom_proba) 1010 { 1011 enum { 1012 PROBA = BIT(0), 1013 RADIUS = BIT(1), 1014 SLICES = BIT(2) 1015 }; 1016 int mask = 0; /* Register the parsed attributes */ 1017 size_t nattrs; 1018 size_t i; 1019 res_T res = RES_OK; 1020 ASSERT(filename && doc && node && geom && geom_proba); 1021 1022 /* Setup the type at the beginning in order to define what arguments should 1023 * be released if a parsing error occurs. */ 1024 geom->type = SCHIFF_SPHERE; 1025 1026 if(node->type != YAML_MAPPING_NODE) { 1027 log_err(filename, node, "expecting a mapping of sphere attributes.\n"); 1028 return RES_BAD_ARG; 1029 } 1030 1031 nattrs = (size_t) 1032 (node->data.mapping.pairs.top - node->data.mapping.pairs.start); 1033 1034 FOR_EACH(i, 0, nattrs) { 1035 yaml_node_t* key; 1036 yaml_node_t* val; 1037 1038 key = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].key); 1039 val = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].value); 1040 ASSERT(key->type == YAML_SCALAR_NODE); 1041 1042 #define SETUP_MASK(Flag, Name) { \ 1043 if(mask & Flag) { \ 1044 log_err(filename, key, "the "Name" is already defined.\n"); \ 1045 return RES_BAD_ARG; \ 1046 } \ 1047 mask |= Flag; \ 1048 } (void)0 1049 1050 /* Probality to sample this geometry */ 1051 if(!strcmp((char*)key->data.scalar.value, "proba")) { 1052 SETUP_MASK(PROBA, "sphere proba"); 1053 res = parse_yaml_double(filename, val, DBL_MIN, DBL_MAX, geom_proba); 1054 1055 /* Sphere radius */ 1056 } else if(!strcmp((char*)key->data.scalar.value, "radius")) { 1057 SETUP_MASK(RADIUS, "sphere radius"); 1058 res = parse_yaml_param_distribution 1059 (filename, doc, val, DBL_MIN, DBL_MAX, &geom->data.sphere.radius); 1060 1061 /* # slices used to discretized the triangular sphere */ 1062 } else if(!strcmp((char*)key->data.scalar.value, "slices")) { 1063 SETUP_MASK(SLICES, "sphere number of slices"); 1064 res = parse_yaml_uint 1065 (filename, val, 4, 32768, &geom->data.sphere.nslices); 1066 1067 /* Error */ 1068 } else { 1069 log_err(filename, key, "unkown sphere attribute `%s'.\n", 1070 key->data.scalar.value); 1071 return RES_BAD_ARG; 1072 } 1073 if(res != RES_OK) return res; 1074 1075 #undef SETUP_MASK 1076 } 1077 1078 /* Ensure the completness of the spherical distribution */ 1079 if(!(mask & RADIUS)) { 1080 log_err(filename, node, "missing the radius attribute.\n"); 1081 return RES_BAD_ARG; 1082 } 1083 if(!(mask & PROBA)) { /* Default proba */ 1084 *geom_proba = 1.0; 1085 } 1086 if(!(mask & SLICES)) { /* Default number of slices */ 1087 geom->data.sphere.nslices = SCHIFF_SPHERE_DEFAULT.nslices; 1088 } 1089 return RES_OK; 1090 } 1091 1092 static res_T 1093 parse_yaml_supershape 1094 (const char* filename, 1095 yaml_document_t* doc, 1096 const yaml_node_t* node, 1097 struct schiff_geometry* geom, 1098 double* geom_proba) 1099 { 1100 enum { 1101 FORMULA0 = BIT(0), 1102 FORMULA1 = BIT(1), 1103 RADIUS_SPHERE = BIT(2), 1104 PROBA = BIT(3), 1105 SLICES = BIT(4) 1106 }; 1107 int mask = 0; /* Register the parsed attributes */ 1108 size_t nattrs; 1109 size_t i; 1110 res_T res = RES_OK; 1111 ASSERT(filename && doc && node && geom && geom_proba); 1112 1113 /* Setup the type at the beginning in order to define what arguments should 1114 * be released if a parsing error occurs. Note that one can define the main 1115 * type or its "equivalent sphere" variation. */ 1116 geom->type = SCHIFF_SUPERSHAPE; 1117 1118 if(node->type != YAML_MAPPING_NODE) { 1119 log_err(filename, node, "expecting a mapping of supershape attributes.\n"); 1120 return RES_BAD_ARG; 1121 } 1122 1123 nattrs = (size_t) 1124 (node->data.mapping.pairs.top - node->data.mapping.pairs.start); 1125 1126 FOR_EACH(i, 0, nattrs) { 1127 yaml_node_t* key; 1128 yaml_node_t* val; 1129 1130 key = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].key); 1131 val = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].value); 1132 ASSERT(key->type == YAML_SCALAR_NODE); 1133 1134 #define SETUP_MASK(Flag, Name) { \ 1135 if(mask & Flag) { \ 1136 log_err(filename, key, "the "Name" is already defined.\n"); \ 1137 return RES_BAD_ARG; \ 1138 } \ 1139 mask |= Flag; \ 1140 } (void)0 1141 1142 /* Geometry probability */ 1143 if(!strcmp((char*)key->data.scalar.value, "proba")) { 1144 SETUP_MASK(PROBA, "supershape proba"); 1145 res = parse_yaml_double(filename, val, DBL_MIN, DBL_MAX, geom_proba); 1146 1147 /* Super shape formula0 */ 1148 } else if(!strcmp((char*)key->data.scalar.value, "formula0")) { 1149 SETUP_MASK(FORMULA0, "supershape formula0"); 1150 res = parse_yaml_superformula 1151 (filename, doc, val, geom->data.supershape.formulas[0]); 1152 1153 /* Super shape formula1 */ 1154 } else if(!strcmp((char*)key->data.scalar.value, "formula1")) { 1155 SETUP_MASK(FORMULA1, "supershape formula1"); 1156 res = parse_yaml_superformula 1157 (filename, doc, val, geom->data.supershape.formulas[1]); 1158 1159 /* Equivalent sphere radius */ 1160 } else if(!strcmp((char*)key->data.scalar.value, "radius_sphere")) { 1161 SETUP_MASK(RADIUS_SPHERE, "equivalent sphere radius of the supershape"); 1162 res = parse_yaml_param_distribution(filename, doc, val, DBL_MIN, DBL_MAX, 1163 &geom->data.supershape.radius_sphere); 1164 1165 /* # slices used to discretized the super shape */ 1166 } else if(!strcmp((char*)key->data.scalar.value, "slices")) { 1167 SETUP_MASK(SLICES, "supershape number of slices"); 1168 res = parse_yaml_uint 1169 (filename, val, 4, 32768, &geom->data.supershape.nslices); 1170 1171 /* Error */ 1172 } else { 1173 log_err(filename, key, "unknown supershape attribute `%s'.\n", 1174 key->data.scalar.value); 1175 res = RES_BAD_ARG; 1176 } 1177 if(res != RES_OK) return res; 1178 1179 #undef SETUP_MASK 1180 } 1181 1182 /* Ensure the completness of the cylinder distribution */ 1183 if(!(mask & FORMULA0)) { 1184 log_err(filename, node, "missing the formula0 attribute.\n"); 1185 return RES_BAD_ARG; 1186 } else if(!(mask & FORMULA1)) { 1187 log_err(filename, node, "missing the formula1 attribute.\n"); 1188 return RES_BAD_ARG; 1189 } 1190 /* Setup the default values if required */ 1191 if(!(mask & PROBA)) { /* Default proba */ 1192 *geom_proba = 1.0; 1193 } 1194 if(!(mask & SLICES)) { /* Default number of slices */ 1195 geom->data.supershape.nslices = SCHIFF_SUPERSHAPE_DEFAULT.nslices; 1196 } 1197 /* Define the geometry type */ 1198 if(!(mask & RADIUS_SPHERE)) { 1199 geom->type = SCHIFF_SUPERSHAPE; 1200 } else { 1201 const struct schiff_supershape* sshape = &geom->data.supershape; 1202 FOR_EACH(i, 0, 6) { 1203 if(sshape->formulas[0][i].distribution != SCHIFF_PARAM_CONSTANT 1204 || sshape->formulas[1][i].distribution != SCHIFF_PARAM_CONSTANT) { 1205 log_err(filename, node, 1206 "the radius_sphere parameter cannot be defined with non constant " 1207 "supershape attributes.\n"); 1208 return RES_BAD_ARG; 1209 } 1210 } 1211 geom->type = SCHIFF_SUPERSHAPE_AS_SPHERE; 1212 } 1213 return RES_OK; 1214 } 1215 1216 static res_T 1217 parse_yaml_geom_distrib 1218 (const char* filename, 1219 yaml_document_t* doc, 1220 yaml_node_t* node, 1221 struct schiff_geometry* geom, 1222 double* proba) 1223 { 1224 res_T res; 1225 yaml_node_t *key, *val; 1226 ASSERT(filename && doc && node && geom && proba); 1227 1228 if(node->type != YAML_MAPPING_NODE 1229 || node->data.mapping.pairs.top - node->data.mapping.pairs.start > 1) { 1230 log_err(filename, node, 1231 "expecting a mapping of the geometry distribution to its parameters\n"); 1232 return RES_BAD_ARG; 1233 } 1234 1235 key = yaml_document_get_node(doc, node->data.mapping.pairs.start[0].key); 1236 val = yaml_document_get_node(doc, node->data.mapping.pairs.start[0].value); 1237 ASSERT(key->type == YAML_SCALAR_NODE); 1238 1239 if(!strcmp((char*)key->data.scalar.value, "ellipsoid")) { 1240 res = parse_yaml_ellipsoid(filename, doc, val, geom, proba); 1241 } else if(!strcmp((char*)key->data.scalar.value, "cylinder")) { 1242 res = parse_yaml_cylinder(filename, doc, val, geom, proba); 1243 } else if(!strcmp((char*)key->data.scalar.value, "sphere")) { 1244 res = parse_yaml_sphere(filename, doc, val, geom, proba); 1245 } else if(!strcmp((char*)key->data.scalar.value, "helical_pipe")) { 1246 res = parse_yaml_helical_pipe(filename, doc, val, geom, proba); 1247 } else if(!strcmp((char*)key->data.scalar.value, "supershape")) { 1248 res = parse_yaml_supershape(filename, doc, val, geom, proba); 1249 } else { 1250 log_err(filename, key, "unknown distribution `%s'.\n", 1251 key->data.scalar.value); 1252 return RES_BAD_ARG; 1253 } 1254 if(res != RES_OK) return res; 1255 return RES_OK; 1256 } 1257 1258 static res_T 1259 parse_yaml 1260 (const char* filename, 1261 struct schiff_geometry** out_geoms, 1262 struct ssp_ranst_discrete** out_ran) 1263 { 1264 yaml_parser_t parser; 1265 yaml_document_t doc; 1266 yaml_node_t* root; 1267 size_t ndistribs = 0; 1268 size_t idistrib; 1269 struct schiff_geometry* geoms = NULL; 1270 double* probas = NULL; 1271 struct ssp_ranst_discrete* ran = NULL; 1272 FILE* file = NULL; 1273 int doc_is_init = 0; 1274 res_T res = RES_OK; 1275 ASSERT(filename && out_geoms && out_ran); 1276 1277 if(!yaml_parser_initialize(&parser)) { 1278 fprintf(stderr, "Couldn't intialise the YAML parser.\n"); 1279 res = RES_UNKNOWN_ERR; 1280 goto exit; 1281 } 1282 1283 file = fopen(filename, "rb"); 1284 if(!file) { 1285 fprintf(stderr, "Couldn't open the YAML file `%s'.\n", filename); 1286 res = RES_IO_ERR; 1287 goto error; 1288 } 1289 1290 yaml_parser_set_input_file(&parser, file); 1291 1292 if(!yaml_parser_load(&parser, &doc)) { 1293 fprintf(stderr, "%s:%lu:%lu: %s.\n", 1294 filename, 1295 (unsigned long)parser.problem_mark.line+1, 1296 (unsigned long)parser.problem_mark.column+1, 1297 parser.problem); 1298 res = RES_IO_ERR; 1299 goto error; 1300 } 1301 doc_is_init = 1; 1302 1303 root = yaml_document_get_root_node(&doc); 1304 if(!root) { 1305 fprintf(stderr, "Unexpected empty file `%s'.\n", filename); 1306 res = RES_BAD_ARG; 1307 goto error; 1308 } 1309 if(root->type == YAML_MAPPING_NODE) { 1310 ndistribs = (size_t) 1311 (root->data.mapping.pairs.top - root->data.mapping.pairs.start); 1312 } else if(root->type == YAML_SEQUENCE_NODE) { 1313 /* Define the number of submitted distributions */ 1314 ndistribs = (size_t) 1315 (root->data.sequence.items.top - root->data.sequence.items.start); 1316 } 1317 1318 if((root->type == YAML_MAPPING_NODE && ndistribs > 1) 1319 || (root->type != YAML_MAPPING_NODE && root->type != YAML_SEQUENCE_NODE)) { 1320 fprintf(stderr, 1321 "Expecting either one or a list of geometry distributions.\n"); 1322 res = RES_BAD_ARG; 1323 goto error; 1324 } 1325 1326 /* Allocate the list geometry distributions */ 1327 if(!sa_add(geoms, ndistribs)) { 1328 log_err(filename, root, 1329 "couldn't allocate up to %lu geometry distributions.\n", 1330 (unsigned long)ndistribs); 1331 res = RES_MEM_ERR; 1332 goto error; 1333 } 1334 FOR_EACH(idistrib, 0, ndistribs) 1335 geoms[idistrib] = SCHIFF_GEOMETRY_NULL; 1336 1337 /* Allocate the per geometry distribution proba */ 1338 if(!sa_add(probas, ndistribs)) { 1339 log_err(filename, root, 1340 "couldn't allocate the list of geometry distribution probabilities.\n"); 1341 res = RES_MEM_ERR; 1342 goto error; 1343 } 1344 1345 /* Create the geometry distribution random variate */ 1346 res = ssp_ranst_discrete_create(&mem_default_allocator, &ran); 1347 if(res != RES_OK) { 1348 log_err(filename, root, 1349 "couldn't allocate the random variate of geometry distributions.\n"); 1350 goto error; 1351 } 1352 1353 /* Only one distribution */ 1354 if(root->type == YAML_MAPPING_NODE) { 1355 res = parse_yaml_geom_distrib(filename, &doc, root, &geoms[0], &probas[0]); 1356 if(res != RES_OK) goto error; 1357 1358 /* List of geometry distributions */ 1359 } else { 1360 double accum_proba = 0; 1361 ASSERT(root->type == YAML_SEQUENCE_NODE); 1362 1363 FOR_EACH(idistrib, 0, ndistribs) { 1364 yaml_node_t* distrib; 1365 1366 distrib = yaml_document_get_node 1367 (&doc, root->data.sequence.items.start[idistrib]); 1368 res = parse_yaml_geom_distrib 1369 (filename, &doc, distrib, &geoms[idistrib], &probas[idistrib]); 1370 if(res != RES_OK) goto error; 1371 1372 accum_proba += probas[idistrib]; 1373 } 1374 /* Normalized the geometry distribution probabilities */ 1375 FOR_EACH(idistrib, 0, ndistribs-1) probas[idistrib] /= accum_proba; 1376 probas[ndistribs-1] = 1.f; /* Handle precision issues */ 1377 } 1378 1379 /* Setup the geometry distributions random variate */ 1380 res = ssp_ranst_discrete_setup(ran, probas, ndistribs); 1381 if(res != RES_OK) { 1382 log_err(filename, root, 1383 "couldn't setup the discrete geometry distributions.\n"); 1384 goto error; 1385 } 1386 1387 exit: 1388 yaml_parser_delete(&parser); 1389 if(doc_is_init) yaml_document_delete(&doc); 1390 if(file) fclose(file); 1391 if(probas) sa_release(probas); 1392 *out_geoms = geoms; 1393 *out_ran = ran; 1394 return res; 1395 error: 1396 if(ran) { 1397 SSP(ranst_discrete_ref_put(ran)); 1398 ran = NULL; 1399 } 1400 if(geoms) { 1401 FOR_EACH(idistrib, 0, ndistribs) { 1402 geometry_release(&geoms[idistrib]); 1403 } 1404 sa_release(geoms); 1405 geoms = NULL; 1406 } 1407 goto exit; 1408 } 1409 1410 /******************************************************************************* 1411 * Local function 1412 ******************************************************************************/ 1413 res_T 1414 schiff_args_init 1415 (struct schiff_args* args, 1416 const int argc, 1417 char** argv) 1418 { 1419 int quiet = 0; 1420 int opt; 1421 res_T res = RES_OK; 1422 ASSERT(argc && argv && args); 1423 1424 *args = SCHIFF_ARGS_NULL; 1425 1426 while((opt = getopt(argc, argv, "a:A:d:Dg:G:hi:l:n:o:qvw:")) != -1) { 1427 switch(opt) { 1428 case 'a': 1429 res = cstr_to_uint(optarg, &args->nangles); 1430 if(res == RES_OK && args->nangles < MIN_NANGLES) { 1431 fprintf(stderr, 1432 "%s: expecting at least "STR(MIN_NANGLES)" scattering angles.\n", 1433 argv[0]); 1434 res = RES_BAD_ARG; 1435 } 1436 break; 1437 case 'A': 1438 res = cstr_to_uint(optarg, &args->nangles_inv); 1439 if(res == RES_OK && args->nangles_inv < MIN_NANGLES_INV) { 1440 fprintf(stderr, 1441 "%s: expecting at least "STR(MIN_NANGLES_INV) 1442 " inverse cumulative phase function values.\n", 1443 argv[0]); 1444 res = RES_BAD_ARG; 1445 } 1446 break; 1447 case 'd': 1448 res = cstr_to_uint(optarg, &args->ninsamps); 1449 if(res == RES_OK && !args->ninsamps) { 1450 fprintf(stderr, "%s: the number of inner samples cannot be null.\n", 1451 argv[0]); 1452 res = RES_BAD_ARG; 1453 } 1454 break; 1455 case 'D': args->discard_large_angles = 1; break; 1456 case 'g': 1457 res = cstr_to_uint(optarg, &args->nrealisations); 1458 if(res == RES_OK && !args->nrealisations) { 1459 fprintf(stderr, "%s: the number of realisations cannot be null.\n", 1460 argv[0]); 1461 res = RES_BAD_ARG; 1462 } 1463 break; 1464 case 'G': res = cstr_to_uint(optarg, &args->ngeoms_dump); break; 1465 case 'h': 1466 usage(stderr); 1467 schiff_args_release(args); 1468 return RES_OK; 1469 case 'i': 1470 res = parse_yaml(optarg, &args->geoms, &args->ran_geoms); 1471 break; 1472 case 'l': 1473 res = cstr_to_double(optarg, &args->characteristic_length); 1474 if(res == RES_OK && args->characteristic_length <= 0.0) 1475 res = RES_BAD_ARG; 1476 break; 1477 case 'n': 1478 res = cstr_to_uint(optarg, &args->nthreads); 1479 if(res == RES_OK && args->nthreads == 0) 1480 res = RES_BAD_ARG; 1481 break; 1482 case 'o': args->output_filename = optarg; break; 1483 case 'q': quiet = 1; break; 1484 case 'w': res = parse_wavelengths(optarg, args); break; 1485 case 'v': 1486 print_version(); 1487 schiff_args_release(args); 1488 return RES_OK; 1489 default: res = RES_BAD_ARG; break; 1490 } 1491 if(res != RES_OK) { 1492 if(optarg) { 1493 fprintf(stderr, "%s: invalid option arguments '%s' -- '%c'\n", 1494 argv[0], optarg, opt); 1495 } 1496 goto error; 1497 } 1498 } 1499 1500 #define MANDATORY(Cond, Name, Opt) { \ 1501 if(!(Cond)) { \ 1502 fprintf(stderr, "%s: %s missing -- option '-%c'\n", argv[0], (Name), (Opt)); \ 1503 res = RES_BAD_ARG; \ 1504 goto error; \ 1505 } \ 1506 } (void)0 1507 MANDATORY(args->geoms != NULL, "geometry distribution", 'i'); 1508 if(args->ngeoms_dump) goto exit; 1509 MANDATORY(args->wavelengths != NULL, "wavelengths", 'w'); 1510 MANDATORY(args->characteristic_length >= 0, "characteristic length", 'l'); 1511 #undef MANDATORY 1512 1513 /* Sort the submitted wavelengths in ascending order */ 1514 qsort(args->wavelengths, sa_size(args->wavelengths), 1515 sizeof(args->wavelengths[0]), cmp_double); 1516 1517 if(optind == argc) { 1518 if(!quiet) { 1519 fprintf(stderr, 1520 "Enter the optical properties. Type ^D (i.e. CTRL+d) to stop:\n"); 1521 } 1522 res = schiff_optical_properties_load_stream(&args->properties, stdin, "stdin"); 1523 } else { 1524 res = schiff_optical_properties_load(&args->properties, argv[optind]); 1525 } 1526 if(res != RES_OK) goto error; 1527 1528 if(!args->properties) { 1529 fprintf(stderr, "%s: missing optical properties.\n", argv[0]); 1530 res = RES_BAD_ARG; 1531 goto error; 1532 } 1533 1534 exit: 1535 return res; 1536 error: 1537 usage(stderr); 1538 schiff_args_release(args); 1539 *args = SCHIFF_ARGS_NULL; 1540 goto exit; 1541 } 1542 1543 void 1544 schiff_args_release(struct schiff_args* args) 1545 { 1546 size_t i, count; 1547 ASSERT(args); 1548 sa_release(args->properties); 1549 sa_release(args->wavelengths); 1550 1551 count = sa_size(args->geoms); 1552 FOR_EACH(i, 0, count) geometry_release(&args->geoms[i]); 1553 sa_release(args->geoms); 1554 if(args->ran_geoms) SSP(ranst_discrete_ref_put(args->ran_geoms)); 1555 args->geoms = NULL; 1556 *args = SCHIFF_ARGS_NULL; 1557 } 1558