schiff

Estimate the radiative properties of soft particless
git clone git://git.meso-star.com/schiff.git
Log | Files | Refs | README | LICENSE

commit b24b7eb2b0e983fc7b1e56a64375b9623db8ae2e
parent 025f2d37528e7a0797d01a5e7b27742034ba392b
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 20 Oct 2015 15:28:53 +0200

Implement the loading of the optical properties

Merge the loaded properties with the properties submitted by the command
line.

Diffstat:
Mcmake/CMakeLists.txt | 8++++++--
Msrc/sbox_schiff.c | 254++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------
Asrc/sbox_schiff_optical_properties.c | 161+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/sbox_schiff_optical_properties.h | 61+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4 files changed, 414 insertions(+), 70 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -60,8 +60,12 @@ set(VERSION_MINOR 3) set(VERSION_PATCH 0) set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) -set(SBOX_SCHIFF_FILES_SRC sbox_schiff.c) -set(SBOX_SCHIFF_FILES_INC sbox_schiff.h) +set(SBOX_SCHIFF_FILES_SRC + sbox_schiff.c + sbox_schiff_optical_properties.c) +set(SBOX_SCHIFF_FILES_INC + sbox_schiff.h + sbox_schiff_optical_properties.h) set(SBOX_SCHIFF_FILES_DOC COPYING.fr COPYING.en README.md) # Prepend each file in the `SBOX_SCHIFF_FILES_<SRC|INC>' list by the absolute diff --git a/src/sbox_schiff.c b/src/sbox_schiff.c @@ -29,6 +29,7 @@ #define _POSIX_C_SOURCE 2 /* getopt support */ #include "sbox_schiff.h" +#include "sbox_schiff_optical_properties.h" #include <star/sbox.h> #include <star/sbox_parse.h> @@ -49,85 +50,89 @@ #define SIGMA 1.0 #define NTHETAS 64 +struct schiff_sphere_args { + const char* material; + const char* output; + struct sbox_schiff_optical_properties* props; + double* wavelengths; + double radius; + double sigma; + unsigned ndirs; + unsigned ngeoms; + unsigned nthetas; +}; + +static const struct schiff_sphere_args SCHIFF_SPHERE_ARGS_NULL = { + NULL, NULL, NULL, NULL, RADIUS, SIGMA, NDIRS, NGEOMS, NTHETAS +}; + /******************************************************************************* * Helper functions ******************************************************************************/ static char -check_optical_properties(const double propperties[4]) -{ - ASSERT(propperties); - return propperties[0] > 0.0 - && propperties[1] > 0.0 - && propperties[2] > 0.0 - && propperties[3] > 0.0; -} - -static char check_wavelength(const double wavelength) { ASSERT(wavelength); return wavelength > 0.0; } -/******************************************************************************* - * Schiff command - ******************************************************************************/ -static res_T -cmd_func_schiff(struct sbox* sbox, const int argc, char** argv, void* ctx) +static int +cmp_properties(const void* op0, const void* op1) { - if(!argc) { - logger_print(&sbox->logger, LOG_ERROR, "Missing argument.\n"); - return RES_BAD_ARG; - } else { - struct sbox_schiff* schiff = ctx; - ASSERT(schiff); - return sbox_cmd_exec(sbox, &schiff->cmd_schiff.sub_cmds, argc-1, argv+1); - } + const struct sbox_schiff_optical_properties* prop0 = op0; + const struct sbox_schiff_optical_properties* prop1 = op1; + if (prop0->W < prop1->W) return -1; + else if(prop0->W > prop1->W) return 1; + else return 0; } static res_T -cmd_func_schiff_sphere(struct sbox* sbox, const int argc, char** argv, void* ctx) +schiff_sphere_parse_args + (struct sbox* sbox, + const int argc, + char** argv, + struct schiff_sphere_args* args) { - const char* material = NULL; - const char* output = NULL; - double* properties = NULL; - double* wavelengths = NULL; - double* prop; + double props[4]; double* wavelength; - double radius = RADIUS; - double sigma = SIGMA; - unsigned ndirs = NDIRS; - unsigned ngeoms = NGEOMS; - unsigned nthetas = NTHETAS; int opt; - size_t i; - (void)ctx; + res_T res = RES_OK; + ASSERT(sbox && argc && argv && args); + *args = SCHIFF_SPHERE_ARGS_NULL; optind = 0; opterr = 0; + while((opt = getopt(argc, argv, ":d:g:i:n:p:r:s:w:")) != -1) { res_T res = RES_OK; switch(opt) { - case 'd': res = sbox_str_to_uint(optarg, &ndirs); break; - case 'g': res = sbox_str_to_uint(optarg, &ngeoms); break; - case 'i': material = optarg; break; - case 'n': res = sbox_str_to_uint(optarg, &nthetas); break; - case '0': output = optarg; break; + case 'd': res = sbox_str_to_uint(optarg, &args->ndirs); break; + case 'g': res = sbox_str_to_uint(optarg, &args->ngeoms); break; + case 'i': args->material = optarg; break; + case 'n': res = sbox_str_to_uint(optarg, &args->nthetas); break; + case '0': args->output = optarg; break; case 'p': - prop = sa_add(properties, 4); - res = sbox_parse_list_double(optarg, prop, 4); - if(res == RES_OK && !check_optical_properties(prop)) res = RES_BAD_ARG; + res = sbox_parse_list_double(optarg, props, 4); + if(res == RES_OK) { + struct sbox_schiff_optical_properties* opt_props = sa_add(args->props, 1); + opt_props->W = props[0]; + opt_props->Nr = props[1]; + opt_props->Kr = props[2]; + opt_props->Ne = props[3]; + if(!sbox_schiff_optical_properties_check(opt_props)) + res = RES_BAD_ARG; + } break; case 'r': - res = sbox_str_to_double(optarg, &radius); - if(res == RES_OK && radius < 0.0) res = RES_BAD_ARG; + res = sbox_str_to_double(optarg, &args->radius); + if(res == RES_OK && args->radius < 0.0) res = RES_BAD_ARG; break; case 's': - res = sbox_str_to_double(optarg, &sigma); - if(res == RES_OK && sigma < 0.0) res = RES_BAD_ARG; + res = sbox_str_to_double(optarg, &args->sigma); + if(res == RES_OK && args->sigma < 0.0) res = RES_BAD_ARG; break; case 'w': - wavelength = sa_add(wavelengths, 1); + wavelength = sa_add(args->wavelengths, 1); res = sbox_str_to_double(optarg, wavelength); if(res == RES_OK && !check_wavelength(*wavelength)) res = RES_BAD_ARG; break; @@ -135,12 +140,14 @@ cmd_func_schiff_sphere(struct sbox* sbox, const int argc, char** argv, void* ctx logger_print(&sbox->logger, LOG_ERROR, "%s: option requires an argument -- '%c'\n", argv[0], argv[optind-1][1]); - return RES_BAD_ARG; + res = RES_BAD_ARG; + break; case '?': logger_print(&sbox->logger, LOG_ERROR, "%s: invalid option -- '%c'\n", argv[0], argv[optind-1][1]); - return RES_BAD_ARG; + res = RES_BAD_ARG; + break; default: res = RES_BAD_ARG; break; } if(res != RES_OK) { @@ -149,29 +156,140 @@ cmd_func_schiff_sphere(struct sbox* sbox, const int argc, char** argv, void* ctx "%s: invalid option arguments '%s' -- '%c'\n", argv[0], optarg, opt); } - return res; + goto error; } } - logger_print(&sbox->logger, LOG_OUTPUT, -"ndirs = %u\nmaterial = %s\nngeoms = %u\noutput = %s\nradius = %g\nsigma = %g\n" -"nsteps = %u\n", - ndirs, material, ngeoms, output, radius, sigma, nthetas); - - logger_print(&sbox->logger, LOG_OUTPUT, "Properties:\n"); - for(i=0; i<sa_size(properties); i += 4) { - logger_print(&sbox->logger, LOG_OUTPUT, " %g %g %g %g\n", - properties[i+0], properties[i+1], properties[i+2], properties[i+3]); +exit: + return res; +error: + sa_release(args->props); + sa_release(args->wavelengths); + *args = SCHIFF_SPHERE_ARGS_NULL; + goto exit; +} + +/* Assume that props0 and props1 are sorted. If two entries have the same + * wavelength the entry of props1 is skipped. */ +static void +merge_properties + (struct sbox_schiff_optical_properties** props0, + const struct sbox_schiff_optical_properties* props1) +{ + const size_t sizeof_prop = sizeof(struct sbox_schiff_optical_properties); + const struct sbox_schiff_optical_properties* src0; + const struct sbox_schiff_optical_properties* src1; + struct sbox_schiff_optical_properties* dst = NULL; + size_t nsrc0, nsrc1; + size_t isrc0 = 0, isrc1 =0; + ASSERT(props0 && props1); + + src0 = *props0; + src1 = props1; + nsrc0 = sa_size(src0); + nsrc1 = sa_size(src1); + + while(isrc0 < nsrc0 && isrc1 < nsrc1) { + if(src0[isrc0].W < src1[isrc1].W) { + memcpy(sa_add(dst, 4), src0 + isrc0, sizeof_prop); + isrc0 += 4; + } else if(src1[isrc1].W < src0[isrc0].W) { + memcpy(sa_add(dst, 4), src1 + isrc1, sizeof_prop); + isrc1 += 4; + } else { /* src0[isrc0].W == src1[isrc1].W */ + memcpy(sa_add(dst, 4), src0 + isrc0, sizeof_prop); + isrc0 += 4; + isrc1 += 4; + } + + if(isrc0 >= nsrc0 && isrc1 < nsrc1) { + const size_t nremains = nsrc1 - isrc1; + memcpy(sa_add(dst, nremains), src1 + isrc1, sizeof_prop * nremains); + isrc1 = nsrc1; + } else if(isrc1 >= nsrc1 && isrc0 < nsrc0) { + const size_t nremains = nsrc0 - isrc0; + memcpy(sa_add(dst, nremains), src0 + isrc0, sizeof_prop * nremains); + isrc0 = nsrc0; + } } - logger_print(&sbox->logger, LOG_OUTPUT, "Wavelengths:\n"); - FOR_EACH(i, 0, sa_size(wavelengths)) { - logger_print(&sbox->logger, LOG_OUTPUT, " %g\n", wavelengths[i]); + sa_release(*props0); + *props0 = dst; +} + +/******************************************************************************* + * Schiff command + ******************************************************************************/ +static res_T +cmd_func_schiff(struct sbox* sbox, const int argc, char** argv, void* ctx) +{ + if(!argc) { + logger_print(&sbox->logger, LOG_ERROR, "Missing argument.\n"); + return RES_BAD_ARG; + } else { + struct sbox_schiff* schiff = ctx; + ASSERT(schiff); + return sbox_cmd_exec(sbox, &schiff->cmd_schiff.sub_cmds, argc-1, argv+1); } - sa_release(properties); - sa_release(wavelengths); - return RES_OK; } static res_T +cmd_func_schiff_sphere(struct sbox* sbox, const int argc, char** argv, void* ctx) +{ + struct schiff_sphere_args args = SCHIFF_SPHERE_ARGS_NULL; + struct sbox_schiff_optical_properties* mtl_props = NULL; + const size_t sizeof_prop = sizeof(struct sbox_schiff_optical_properties); + res_T res = RES_OK; + (void)ctx; + + res = schiff_sphere_parse_args(sbox, argc, argv, &args); + if(res != RES_OK) goto error; + + if(!args.wavelengths) { + logger_print(&sbox->logger, LOG_ERROR, + "Missing the wavelength argument(s).\n"); + res = RES_BAD_ARG; + goto error; + } + + /* Sort the submitted properties in ascending order with respect to the + * wavelength */ + if(args.props) { + qsort(args.props, sa_size(args.props), sizeof_prop, cmp_properties); + } + /* Load the material file and sort the its properties in ascending order with + * respect to the wavelength */ + if(args.material) { + res = sbox_schiff_optical_properties_load(sbox, args.material, &mtl_props); + if(res != RES_OK) goto error; + qsort(mtl_props, sa_size(mtl_props), sizeof_prop, cmp_properties); + } + + /* Define the optical properties to use */ + if(mtl_props && !args.props) { + args.props = mtl_props; + mtl_props = NULL; + } else if(args.props && mtl_props) { + merge_properties(&args.props, mtl_props); + mtl_props = NULL; + } + + if(!args.props) { + logger_print(&sbox->logger, LOG_ERROR, + "Radiative properties are not defined.\n"); + res = RES_BAD_ARG; + goto error; + } + +exit: + sa_release(args.props); + sa_release(args.wavelengths); + sa_release(mtl_props); + return res; +error: + goto exit; +} + + +static res_T cmd_func_help_schiff(struct sbox* sbox, const int argc, char** argv, void* ctx) { struct sbox_schiff* schiff = ctx; @@ -232,7 +350,7 @@ cmd_func_help_schiff_sphere " is the wavelength in meter, \"Nr\" and \"Kr\" are respectively the\n", " real and imaginary parts of the relative refractive index and\n", " \"Ne\" is the refractive index of the medium. Overwrite the\n", -" properties loaded with the -i option \n", +" properties loaded with the -i option.\n", "-r RADIUS \"radius\" argument of the lognormal distribution Default is "STR(RADIUS)".\n", "-s SIGMA \"sigma\" argument of the lognormal distribution. Default is "STR(SIGMA)".\n", "-w WAVELENGTH Wavelength to integrate.\n" diff --git a/src/sbox_schiff_optical_properties.c b/src/sbox_schiff_optical_properties.c @@ -0,0 +1,161 @@ +/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com) + * + * This software is governed by the CeCILL license under French law and + * abiding by the rules of distribution of free software. You can use, + * modify and/or redistribute the software under the terms of the CeCILL + * license as circulated by CEA, CNRS and INRIA at the following URL + * "http://www.cecill.info". + * + * As a counterpart to the access to the source code and rights to copy, + * modify and redistribute granted by the license, users are provided only + * with a limited warranty and the software's author, the holder of the + * economic rights, and the successive licensors have only limited + * liability. + * + * In this respect, the user's attention is drawn to the risks associated + * with loading, using, modifying and/or developing or reproducing the + * software by the user in light of its specific status of free software, + * that may mean that it is complicated to manipulate, and that also + * therefore means that it is reserved for developers and experienced + * professionals having in-depth computer knowledge. Users are therefore + * encouraged to load and test the software's suitability as regards their + * requirements in conditions enabling the security of their systems and/or + * data to be ensured and, more generally, to use and operate it in the + * same conditions as regards security. + * + * The fact that you are presently reading this means that you have had + * knowledge of the CeCILL license and that you accept its terms. */ + +#include "sbox_schiff_optical_properties.h" + +#include <star/sbox.h> +#include <star/sbox_parse.h> + +#include <rsys/stretchy_array.h> + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static res_T +parse_optical_properties + (struct sbox* sbox, + const char* filename, + const size_t iline, + const char* str, + struct sbox_schiff_optical_properties* out_props) +{ + char buf[128]; + char* tk; + int i = 0; + double props[4]; + res_T res = RES_OK; + ASSERT(sbox && filename && str && out_props); + + if(strlen(str) >= sizeof(buf) - 1/*NULL char*/) { + logger_print(&sbox->logger, LOG_WARNING, + "%s:%lu: not enough memory: cannot parse the properties `%s'.\n", + filename, (unsigned long)iline, str); + return RES_MEM_ERR; + } + + strncpy(buf, str, sizeof(buf)); + tk = strtok(buf, " \t"); + while(tk) { + res = sbox_str_to_double(tk, props + i); + if(res != RES_OK) { + logger_print(&sbox->logger, LOG_WARNING, + "%s:%lu: cannot read the optical property `%s'.\n", + filename, (unsigned long)iline, tk); + return res; + } + tk = strtok(NULL, " \t"); + ++i; + } + if(i < 4) { + logger_print(&sbox->logger, LOG_WARNING, +"%s:%lu: not enough optical properties.\n" +"Expect 4 parameters while %d is submitted - `%s'.\n", + filename, (unsigned long)iline, i, str); + return RES_BAD_ARG; + } + out_props->W = props[0]; + out_props->Nr = props[1]; + out_props->Kr = props[2]; + out_props->Ne = props[3]; + if(!sbox_schiff_optical_properties_check(out_props)) { + logger_print(&sbox->logger, LOG_WARNING, + "%s:%lu: invalid optical properties `%s'.\n", + filename, (unsigned long)iline, str); + return RES_BAD_ARG; + } + return RES_OK; +} + +/******************************************************************************* + * Local function + ******************************************************************************/ +res_T +sbox_schiff_optical_properties_load + (struct sbox* sbox, + const char* filename, + struct sbox_schiff_optical_properties** out_props) +{ + struct sbox_schiff_optical_properties* props = NULL; + char* buf = NULL; + const size_t buf_chunk = 128; + size_t iline; + FILE* fp; + res_T res = RES_OK; + ASSERT(sbox && filename && out_props); + + fp = fopen(filename, "r"); + if(!fp) { + logger_print(&sbox->logger, LOG_ERROR, + "Cannot open the file of Schiff optical properties `%s'.\n", filename); + res = RES_IO_ERR; + goto error; + } + + buf = sa_add(buf, buf_chunk); + iline = 1; + while(fgets(buf, (int)sa_size(buf), fp)) { + char* line; + size_t last_char; + struct sbox_schiff_optical_properties tmp_props; + + while(!strrchr(buf, '\n')) { /* Ensure that the whole line is read */ + if(!fgets(sa_add(buf, buf_chunk), (int)buf_chunk, fp)) /* EOF */ + break; + } + + /* Remove leading spaces */ + line = buf; + while((*line == ' ' || *line == '\t') && *line != '\0') ++line; + + /* Remove newline character(s) */ + last_char = strlen(line); + while(last_char-- && (line[last_char]=='\n' || line[last_char]=='\r')); + line[last_char + 1] = '\0'; + + if(*line == '\0' /*Empty line*/ && *line == '#' /* Comment */) + continue; + + /* Read optical properties */ + parse_optical_properties(sbox, filename, iline, line, &tmp_props); + if(res == RES_OK) { /* *NO* error */ + sa_push(props, tmp_props); + } + } + +exit: + sa_release(buf); + *out_props = props; + return res; +error: + if(props) { + sa_release(props); + props = NULL; + } + goto exit; +} + diff --git a/src/sbox_schiff_optical_properties.h b/src/sbox_schiff_optical_properties.h @@ -0,0 +1,61 @@ +/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com) + * + * This software is governed by the CeCILL license under French law and + * abiding by the rules of distribution of free software. You can use, + * modify and/or redistribute the software under the terms of the CeCILL + * license as circulated by CEA, CNRS and INRIA at the following URL + * "http://www.cecill.info". + * + * As a counterpart to the access to the source code and rights to copy, + * modify and redistribute granted by the license, users are provided only + * with a limited warranty and the software's author, the holder of the + * economic rights, and the successive licensors have only limited + * liability. + * + * In this respect, the user's attention is drawn to the risks associated + * with loading, using, modifying and/or developing or reproducing the + * software by the user in light of its specific status of free software, + * that may mean that it is complicated to manipulate, and that also + * therefore means that it is reserved for developers and experienced + * professionals having in-depth computer knowledge. Users are therefore + * encouraged to load and test the software's suitability as regards their + * requirements in conditions enabling the security of their systems and/or + * data to be ensured and, more generally, to use and operate it in the + * same conditions as regards security. + * + * The fact that you are presently reading this means that you have had + * knowledge of the CeCILL license and that you accept its terms. */ + +#ifndef SBOX_SCHIFF_OPTICAL_PROPERTIES_H +#define SBOX_SCHIFF_OPTICAL_PROPERTIES_H + +#include <rsys/rsys.h> + +struct sbox; + +struct sbox_schiff_optical_properties { + double W; /* Wavelength */ + double Nr; /* Real part of the relative refractive index */ + double Kr; /* Imaginary part of the relative refractive index */ + double Ne; /* Refractive index of the medium */ +}; + +static INLINE char +sbox_schiff_optical_properties_check + (const struct sbox_schiff_optical_properties* props) +{ + ASSERT(props); + return props->W > 0.0 + && props->Nr > 0.0 + && props->Kr > 0.0 + && props->Ne > 0.0; +} + +extern LOCAL_SYM res_T +sbox_schiff_optical_properties_load + (struct sbox* sbox, + const char* filename, + struct sbox_schiff_optical_properties** props); /* Stretchy array */ + +#endif /* SBOX_SCHIFF_OPTICAL_PROPERTIES_H */ +