schiff

schiff_optical_properties.c (6121B)

---

 /* Copyright (C) 2015-2016 CNRS
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation, either version 3 of the License, or
  * (at your option) any later version.
  *
  * This program is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  * GNU General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
 
 #include "schiff_optical_properties.h"
 #include "schiff_streamline.h"
 
 #include <rsys/algorithm.h>
 #include <rsys/cstr.h>
 #include <rsys/stretchy_array.h>
 
 #include <star/sschiff.h>
 
 /*******************************************************************************
  * Helper functions
  ******************************************************************************/
 static res_T
 parse_optical_properties
   (struct schiff_optical_properties* out_props,
    const char* filename,
    const size_t iline,
    const char* str)
 {
   char buf[128];
   char* tk;
   int i = 0;
   double props[4];
   res_T res = RES_OK;
   ASSERT(filename && str && out_props);
 
   if(strlen(str) >= sizeof(buf) - 1/*NULL char*/) {
     fprintf(stderr,
       "%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));
   for(i=0, tk=strtok(buf, " \t"); tk && i < 4; ++i, tk=strtok(NULL, " \t")) {
     res = cstr_to_double(tk, props + i);
     if(res != RES_OK) {
       fprintf(stderr,
         "%s:%lu: cannot read the optical property `%s'.\n",
         filename, (unsigned long)iline, tk);
       return res;
     }
   }
   if(i < 4) {
     fprintf(stderr,
       "%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->Ne = props[3];
   out_props->Nr = props[1] / out_props->Ne; /* Absolute to relative */
   out_props->Kr = props[2] / out_props->Ne; /* Absolute to relative */
   if(!schiff_optical_properties_check(out_props)) {
     fprintf(stderr, "%s:%lu: invalid optical properties `%s'.\n",
       filename, (unsigned long)iline, str);
     return RES_BAD_ARG;
   }
   return RES_OK;
 }
 
 static FINLINE int
 cmp_properties(const void* op0, const void* op1)
 {
   const struct schiff_optical_properties* prop0 = op0;
   const struct schiff_optical_properties* prop1 = op1;
   return prop0->W < prop1->W ? -1 : (prop0->W > prop1->W ? 1 : 0);
 }
 
 static FINLINE int
 cmp_wavelength_to_properties(const void* wavelength, const void* op)
 {
   const double* W = wavelength;
   const struct schiff_optical_properties* prop = op;
   return *W < prop->W ? -1 : (*W > prop->W ? 1 : 0);
 }
 
 /*******************************************************************************
  * Local function
  ******************************************************************************/
 res_T
 schiff_optical_properties_load
   (struct schiff_optical_properties** props, const char* filename)
 {
   FILE* fp = NULL;
   res_T res = RES_OK;
   ASSERT(filename && props);
 
   fp = fopen(filename, "r");
   if(!fp) {
     fprintf(stderr, "Cannot open the file of optical properties `%s'.\n",
       filename);
     res = RES_IO_ERR;
     goto error;
   }
 
   res = schiff_optical_properties_load_stream(props, fp, filename);
   if(res != RES_OK) goto error;
 
 exit:
   if(fp) fclose(fp);
   return res;
 error:
   goto exit;
 }
 
 res_T
 schiff_optical_properties_load_stream
   (struct schiff_optical_properties** out_props,
    FILE* stream,
    const char* stream_name)
 {
   struct schiff_optical_properties* props = NULL;
   struct schiff_streamline streamline;
   size_t iline;
   char* line;
   res_T res = RES_OK;
   ASSERT(out_props && stream && stream_name);
 
   schiff_streamline_init(&streamline);
 
   for(iline=1;
     (res = schiff_streamline_read(&streamline, stream, &line)) != RES_EOF;
     ++iline) {
     struct schiff_optical_properties tmp_props;
 
     if(strcspn(line, "# \t") == 0) /* Empty line */
       continue;
 
     /* Read optical properties */
     line = strtok(line, "#");
     res = parse_optical_properties(&tmp_props, stream_name, iline, line);
     if(res == RES_OK) { /* *NO* error */
       sa_push(props, tmp_props);
     }
   }
 
   /* Sort the optical properties in ascending order with respect to the
    * wavelength */
   qsort(props, sa_size(props),
     sizeof(struct schiff_optical_properties), cmp_properties);
 
   *out_props = props;
   schiff_streamline_release(&streamline);
   return RES_OK;
 }
 
 void
 schiff_optical_properties_fetch
   (struct schiff_optical_properties* properties, /* Stetchy array */
    const double wavelength,
    struct sschiff_material_properties* props)
 {
   const size_t nproperties = sa_size(properties);
   struct schiff_optical_properties* find;
   double Ne, Kr, Nr;
   ASSERT(properties && props && wavelength > 0.0);
 
   /* Assume that properties are sorted in ascending order with respect to the
    * wavelength. */
   find = search_lower_bound(&wavelength, properties, nproperties,
     sizeof(struct schiff_optical_properties), cmp_wavelength_to_properties);
   ASSERT(!find || find->W >= wavelength);
 
   if(!find) /* Clamp to upper wavelength */
     find = properties + (nproperties - 1);
 
   if(find == properties || find == (properties + (nproperties - 1))
   || eq_eps(find->W, wavelength, 1.e-12)) {
     Ne = find[0].Ne;
     Kr = find[0].Kr;
     Nr = find[0].Nr;
   } else { /* Linear interpolation of optical properties */
     const double d = find[0].W - find[-1].W;
     const double u = (wavelength - find[-1].W) / d;
     const double v = 1.0 - u;
     Ne = u * find[0].Ne + v * find[-1].Ne;
     Kr = u * find[0].Kr + v * find[-1].Kr;
     Nr = u * find[0].Nr + v * find[-1].Nr;
   }
 
   props->medium_refractive_index = Ne;
   props->relative_imaginary_refractive_index = Kr;
   props->relative_real_refractive_index = Nr;
 }