atrri

Refractive indices varying with wavelength
git clone git://git.meso-star.fr/atrri.git
Log | Files | Refs | README | LICENSE

atrri.c (11061B)


      1 /* Copyright (C) 2022, 2023 |Méso|Star> (contact@meso-star.com)
      2  * Copyright (C) 2021 Centre National de la Recherche Scientifique
      3  *
      4  * This program is free software: you can redistribute it and/or modify
      5  * it under the terms of the GNU General Public License as published by
      6  * the Free Software Foundation, either version 3 of the License, or
      7  * (at your option) any later version.
      8  *
      9  * This program is distributed in the hope that it will be useful,
     10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     12  * GNU General Public License for more details.
     13  *
     14  * You should have received a copy of the GNU General Public License
     15  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     16 
     17 #define _POSIX_C_SOURCE 200112L /* strtok_r support */
     18 
     19 #include "atrri.h"
     20 #include "atrri_c.h"
     21 #include "atrri_log.h"
     22 
     23 #include <rsys/algorithm.h>
     24 #include <rsys/cstr.h>
     25 #include <rsys/text_reader.h>
     26 
     27 #include <errno.h>
     28 #include <string.h>
     29 
     30 /*******************************************************************************
     31  * Local functions
     32  ******************************************************************************/
     33 static INLINE int
     34 cmp_refractive_index(const void* key, const void* entry)
     35 {
     36   const double wavelength = *((const double*)key);
     37   const struct atrri_refractive_index* id = entry;
     38 
     39   if(wavelength < id->wavelength) {
     40     return -1;
     41   } else if(wavelength > id->wavelength) {
     42     return 1;
     43   } else {
     44     return 0;
     45   }
     46 }
     47 
     48 static INLINE void
     49 reset_atrri(struct atrri* atrri)
     50 {
     51   ASSERT(atrri);
     52   darray_refractive_id_clear(&atrri->refract_ids);
     53 }
     54 
     55 static res_T
     56 parse_line(struct atrri* atrri, struct txtrdr* txtrdr)
     57 {
     58   struct atrri_refractive_index refract_id = ATRRI_REFRACTIVE_INDEX_NULL;
     59   char* tk = NULL;
     60   char* tk_ctx = NULL;
     61   size_t nrefract_ids = 0;
     62   res_T res = RES_OK;
     63   ASSERT(atrri && txtrdr);
     64 
     65   /* Parse the wavelength */
     66   tk = strtok_r(txtrdr_get_line(txtrdr), " \t", &tk_ctx);
     67   ASSERT(tk); /* A line should exist since it was parsed by txtrdr */
     68   res = cstr_to_double(tk, &refract_id.wavelength);
     69   if(res == RES_OK && refract_id.wavelength < 0) res = RES_BAD_ARG;
     70   if(res != RES_OK) {
     71     log_err(atrri, "%s:%lu: invalid wavelength `%s' -- %s.\n",
     72       txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
     73       tk, res_to_cstr(res));
     74     goto error;
     75   }
     76 
     77   nrefract_ids = darray_refractive_id_size_get(&atrri->refract_ids);
     78   if(nrefract_ids != 0) {
     79     /* Check that the indices are sorted in ascending order wrt wavelength */
     80     const struct atrri_refractive_index* prev_refract_id =
     81       darray_refractive_id_cdata_get(&atrri->refract_ids) + nrefract_ids - 1;
     82     if(prev_refract_id->wavelength >= refract_id.wavelength) {
     83       log_err(atrri,
     84         "%s:%lu: refractive indices must be sorted in ascending order "
     85         "acording to their corresponding wavelength.\n",
     86         txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr));
     87       res = RES_BAD_ARG;
     88       goto error;
     89     }
     90   }
     91 
     92   /* Parse the real part */
     93   tk = strtok_r(NULL, " \t", &tk_ctx);
     94   if(!tk) {
     95     log_err(atrri,
     96       "%s:%lu: could not read the real part of the refractive index.\n",
     97       txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr));
     98     res = RES_BAD_ARG;
     99     goto error;
    100   }
    101   res = cstr_to_double(tk, &refract_id.n);
    102   if(res == RES_OK && refract_id.n < 0) res = RES_BAD_ARG;
    103   if(res != RES_OK) {
    104     log_err(atrri,
    105       "%s:%lu: invalid real part of the refractive index `%s' -- %s.\n",
    106       txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
    107       tk, res_to_cstr(res));
    108     goto error;
    109   }
    110 
    111   /* Parse the imaginary part */
    112   tk = strtok_r(NULL, " \t", &tk_ctx);
    113   if(!tk) {
    114     log_err(atrri,
    115       "%s:%lu: could not read the imaginary part of the refractive index.\n",
    116       txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr));
    117     res = RES_BAD_ARG;
    118     goto error;
    119   }
    120   res = cstr_to_double(tk, &refract_id.kappa);
    121   if(res == RES_OK && refract_id.kappa < 0) res = RES_BAD_ARG;
    122   if(res != RES_OK) {
    123     log_err(atrri,
    124       "%s:%lu: invalid imaginary part of the refractive index `%s' -- %s.\n",
    125       txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
    126       tk, res_to_cstr(res));
    127     goto error;
    128   }
    129 
    130   tk = strtok_r(NULL, " \t", &tk_ctx);
    131   if(tk) {
    132     log_warn(atrri, "%s:%lu: unexpected text `%s'.\n",
    133       txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr), tk);
    134   }
    135 
    136   /* Register the parsed refractive index */
    137   res = darray_refractive_id_push_back(&atrri->refract_ids, &refract_id);
    138   if(res != RES_OK) {
    139     log_err(atrri, "%s:%lu: could not register the refractive index -- %s.\n",
    140       txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
    141       res_to_cstr(res));
    142     goto error;
    143   }
    144 
    145 exit:
    146   return res;
    147 error:
    148   goto exit;
    149 }
    150 
    151 static res_T
    152 load_stream(struct atrri* atrri, FILE* stream, const char* stream_name)
    153 {
    154   struct txtrdr* txtrdr = NULL;
    155   res_T res = RES_OK;
    156   ASSERT(atrri && stream && stream_name);
    157 
    158   reset_atrri(atrri);
    159 
    160   res = txtrdr_stream(atrri->allocator, stream, stream_name, '#', &txtrdr);
    161   if(res != RES_OK) {
    162     log_err(atrri, "%s: could not create the text reader -- %s.\n",
    163       stream_name, res_to_cstr(res));
    164   }
    165 
    166   for(;;) {
    167     res = txtrdr_read_line(txtrdr);
    168     if(res != RES_OK) {
    169       log_err(atrri, "%s: could not read the line `%lu' -- %s.\n",
    170         txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
    171         res_to_cstr(res));
    172       goto error;
    173     }
    174 
    175     if(!txtrdr_get_cline(txtrdr)) break; /* No more parsed line */
    176 
    177     res = parse_line(atrri, txtrdr);
    178     if(res != RES_OK) goto error;
    179   }
    180 
    181   if(darray_refractive_id_size_get(&atrri->refract_ids) == 0) {
    182     log_err(atrri, "%s: empty file.\n", stream_name);
    183     res = RES_BAD_ARG;
    184     goto error;
    185   }
    186 
    187 exit:
    188   if(txtrdr) txtrdr_ref_put(txtrdr);
    189   return res;
    190 error:
    191   goto exit;
    192 }
    193 
    194 static void
    195 release_atrri(ref_T* ref)
    196 {
    197   struct atrri* atrri;
    198   ASSERT(ref);
    199   atrri = CONTAINER_OF(ref, struct atrri, ref);
    200   darray_refractive_id_release(&atrri->refract_ids);
    201   if(atrri->logger == &atrri->logger__) logger_release(&atrri->logger__);
    202   MEM_RM(atrri->allocator, atrri);
    203 }
    204 
    205 /*******************************************************************************
    206  * Exported functions
    207  ******************************************************************************/
    208 res_T
    209 atrri_create
    210   (struct logger* logger, /* NULL <=> use default logger */
    211    struct mem_allocator* mem_allocator, /* NULL <=> use default allocator */
    212    const int verbose, /* Verbosity level */
    213    struct atrri** out_atrri)
    214 {
    215   struct atrri* atrri = NULL;
    216   struct mem_allocator* allocator = NULL;
    217   res_T res = RES_OK;
    218 
    219   if(!out_atrri) {
    220     res = RES_BAD_ARG;
    221     goto error;
    222   }
    223 
    224   allocator = mem_allocator ? mem_allocator : &mem_default_allocator;
    225   atrri = MEM_CALLOC(allocator, 1, sizeof(*atrri));
    226   if(!atrri) {
    227     if(verbose) {
    228       #define ERR_STR "Could not allocate the AtrRI device.\n"
    229       if(logger) {
    230         logger_print(logger, LOG_ERROR, ERR_STR);
    231       } else {
    232         fprintf(stderr, MSG_ERROR_PREFIX ERR_STR);
    233       }
    234       #undef ERR_STR
    235     }
    236     res = RES_MEM_ERR;
    237     goto error;
    238   }
    239   ref_init(&atrri->ref);
    240   atrri->allocator = allocator;
    241   atrri->verbose = verbose;
    242   if(logger) {
    243     atrri->logger = logger;
    244   } else {
    245     setup_log_default(atrri);
    246   }
    247 
    248   darray_refractive_id_init(atrri->allocator, &atrri->refract_ids);
    249 
    250 exit:
    251   if(out_atrri) *out_atrri = atrri;
    252   return res;
    253 error:
    254   if(atrri) {
    255     ATRRI(ref_put(atrri));
    256     atrri = NULL;
    257   }
    258   goto exit;
    259 }
    260 
    261 res_T
    262 atrri_ref_get(struct atrri* atrri)
    263 {
    264   if(!atrri) return RES_BAD_ARG;
    265   ref_get(&atrri->ref);
    266   return RES_OK;
    267 }
    268 
    269 res_T
    270 atrri_ref_put(struct atrri* atrri)
    271 {
    272   if(!atrri) return RES_BAD_ARG;
    273   ref_put(&atrri->ref, release_atrri);
    274   return RES_OK;
    275 }
    276 
    277 res_T
    278 atrri_load(struct atrri* atrri, const char* path)
    279 {
    280   FILE* file = NULL;
    281   res_T res = RES_OK;
    282 
    283   if(!atrri || !path) {
    284     res = RES_BAD_ARG;
    285     goto error;
    286   }
    287 
    288   file = fopen(path, "r");
    289   if(!file) {
    290     log_err(atrri, "%s: error opening file -- %s.\n",
    291       path, strerror(errno));
    292     res = RES_IO_ERR;
    293     goto error;
    294   }
    295 
    296   res = load_stream(atrri, file, path);
    297   if(res != RES_OK) goto error;
    298 
    299 exit:
    300   if(file) fclose(file);
    301   return res;
    302 error:
    303   goto exit;
    304 }
    305 
    306 res_T
    307 atrri_load_stream(struct atrri* atrri, FILE* stream, const char* stream_name)
    308 {
    309   res_T res = RES_OK;
    310 
    311   if(!atrri || !stream) {
    312     res = RES_BAD_ARG;
    313     goto error;
    314   }
    315 
    316   res = load_stream(atrri, stream, stream_name ? stream_name : "stream");
    317   if(res != RES_OK) goto error;
    318 
    319 exit:
    320   return res;
    321 error:
    322   goto exit;
    323 }
    324 
    325 res_T
    326 atrri_get_desc(const struct atrri* atrri, struct atrri_desc* desc)
    327 {
    328   if(!atrri || !desc) return RES_BAD_ARG;
    329   desc->indices = darray_refractive_id_cdata_get(&atrri->refract_ids);
    330   desc->nindices = darray_refractive_id_size_get(&atrri->refract_ids);
    331   return RES_OK;
    332 }
    333 
    334 res_T
    335 atrri_fetch_refractive_index
    336   (const struct atrri* atrri,
    337    const double wavelength,
    338    struct atrri_refractive_index* refract_id)
    339 {
    340   struct atrri_desc desc = ATRRI_DESC_NULL;
    341   struct atrri_refractive_index* found = NULL;
    342   res_T res = RES_OK;
    343 
    344   if(!atrri || wavelength < 0 || !refract_id) {
    345     res = RES_BAD_ARG;
    346     goto error;
    347   }
    348 
    349   ATRRI(get_desc(atrri, &desc));
    350 
    351   /* Search for the registered refractive index whose corresponding wavelength
    352    * is greater than or equal to the submitted wavelength */
    353   found = search_lower_bound(&wavelength, desc.indices, desc.nindices,
    354     sizeof(*desc.indices), cmp_refractive_index);
    355 
    356    if(!found) {
    357     /* All refractive indices are defined for wavelengths less than the
    358      * submitted one. Clamp to the last registered refractive index  */
    359     ASSERT(desc.nindices);
    360     *refract_id = desc.indices[desc.nindices-1];
    361 
    362   } else if(found == desc.indices) {
    363     /* The submitted wavelength is less than the first wavelength for which a
    364      * refractive index is registered. Clamp to the first registered refractive
    365      * index */
    366     *refract_id = desc.indices[0];
    367 
    368   } else  {
    369     const size_t inext = (size_t)(found - desc.indices);
    370     const size_t iprev = inext - 1;
    371     const struct atrri_refractive_index* next = desc.indices + inext;
    372     const struct atrri_refractive_index* prev = desc.indices + iprev;
    373     ASSERT(next->wavelength >= wavelength);
    374     ASSERT(prev->wavelength < wavelength);
    375 
    376     if(next->wavelength == wavelength) {
    377       /* The submitted wavelength strictly match a registered wavelength */
    378       *refract_id = *next;
    379 
    380     } else {
    381       /* Linearly interpolate the refractive index */
    382       const double len = next->wavelength - prev->wavelength;
    383       const double u = CLAMP((wavelength - prev->wavelength) / len, 0, 1);
    384       ASSERT(len > 0);
    385       refract_id->wavelength = wavelength;
    386       refract_id->n = u*(prev->n - next->n) + next->n;
    387       refract_id->kappa = u*(prev->kappa - next->kappa) + next->kappa;
    388     }
    389   }
    390 
    391 exit:
    392   return res;
    393 error:
    394   goto exit;
    395 }