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