rnsf

Define and load a phase function data format
git clone git://git.meso-star.fr/rnsf.git
Log | Files | Refs | README | LICENSE

rnsf_fetch_phase_fn.c (4338B)

      1 /* Copyright (C) 2022, 2023 Centre National de la Recherche Scientifique
      2  * Copyright (C) 2022, 2023 Institut Pierre-Simon Laplace
      3  * Copyright (C) 2022, 2023 Institut de Physique du Globe de Paris
      4  * Copyright (C) 2022, 2023 |Méso|Star> (contact@meso-star.com)
      5  * Copyright (C) 2022, 2023 Observatoire de Paris
      6  * Copyright (C) 2022, 2023 Université de Reims Champagne-Ardenne
      7  * Copyright (C) 2022, 2023 Université de Versaille Saint-Quentin
      8  * Copyright (C) 2022, 2023 Université Paul Sabatier
      9  *
     10  * This program is free software: you can redistribute it and/or modify
     11  * it under the terms of the GNU General Public License as published by
     12  * the Free Software Foundation, either version 3 of the License, or
     13  * (at your option) any later version.
     14  *
     15  * This program is distributed in the hope that it will be useful,
     16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     18  * GNU General Public License for more details.
     19  *
     20  * You should have received a copy of the GNU General Public License
     21  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     22 
     23 #include "rnsf.h"
     24 #include "rnsf_c.h"
     25 #include "rnsf_log.h"
     26 #include "rnsf_phase_fn.h"
     27 
     28 #include <rsys/algorithm.h>
     29 
     30 /*******************************************************************************
     31  * Helper functions
     32  ******************************************************************************/
     33 static INLINE int
     34 cmp_phase(const void* key, const void* item)
     35 {
     36   const double wavelength = (*(const double*)key);
     37   const struct rnsf_phase_fn* phase = item;
     38   if(wavelength < phase->wlen[0]) {
     39     return -1;
     40   } else if(wavelength > phase->wlen[1]) {
     41     return 1;
     42   } else {
     43     return 0;
     44   }
     45 }
     46 
     47 /*******************************************************************************
     48  * Exported functions
     49  ******************************************************************************/
     50 res_T
     51 rnsf_fetch_phase_fn
     52   (const struct rnsf* rnsf,
     53    const double wavelength,
     54    const double sample, /* In [0, 1[ */
     55    size_t* out_iphase_fn)
     56 {
     57   const struct rnsf_phase_fn* phases = NULL;
     58   const struct rnsf_phase_fn* phase = NULL;
     59   size_t nphases = 0;
     60   size_t iphase_fn = 0;
     61   res_T res = RES_OK;
     62 
     63   if(!rnsf
     64   || wavelength < 0
     65   || sample < 0
     66   || sample >= 1
     67   || !out_iphase_fn) {
     68     res = RES_BAD_ARG;
     69     goto error;
     70   }
     71 
     72   phases = darray_phase_fn_cdata_get(&rnsf->phases);
     73   nphases = darray_phase_fn_size_get(&rnsf->phases);
     74 
     75   if(!nphases) {
     76     log_warn(rnsf, "%s: no phase function to fetch.\n", FUNC_NAME);
     77     *out_iphase_fn = 0;
     78     goto exit;
     79   }
     80 
     81   phase = search_lower_bound
     82     (&wavelength, phases, nphases, sizeof(*phases), cmp_phase);
     83 
     84   /* All phase functions are set for a wavelength greater or equal to the
     85    * submitted wavelength */
     86   if(phase == phases) {
     87     ASSERT(phases[0].wlen[0] > wavelength 
     88        || (phases[0].wlen[0] <= wavelength && wavelength <= phases[0].wlen[1]));
     89     iphase_fn = 0;
     90   }
     91 
     92   /* All phase functions are set for a wavelength less than the submitted
     93    * wavelength */
     94   else if(!phase) {
     95     ASSERT(phases[nphases-1].wlen[1] < wavelength);
     96     iphase_fn = nphases - 1;
     97   }
     98 
     99   /* The wavelength range of the found phase function overlaps the submitted
    100    * wavelength */
    101   else if(wavelength >= phase->wlen[0] && wavelength <= phase->wlen[1]) {
    102     iphase_fn = (size_t)(phase - phases);
    103   }
    104 
    105   /* The found phase function is defined for a wavelength greater or equal to
    106    * the submitted wavelength */
    107   else {
    108     /* Compute the linear interpolation parameter that defines the submitted
    109      * wavelength wrt to the wavelenths boundaries of the phase functions
    110      * surrounding it */
    111     const struct rnsf_phase_fn* phase0 = phase-1;
    112     const struct rnsf_phase_fn* phase1 = phase;
    113     const double u =
    114       (phase1->wlen[0] - wavelength)
    115     / (phase1->wlen[0] - phase0->wlen[1]);
    116 
    117     /* Consider the linear interplation parameter previously defined as a
    118      * probability and use the submitted sample to select one of the phase
    119      * function surrounding the submitted wavelength. */
    120     if(sample < u) {
    121       iphase_fn = (size_t)(phase - phases - 1);
    122     } else {
    123       iphase_fn = (size_t)(phase - phases);
    124     }
    125   }
    126 
    127   ASSERT(iphase_fn < nphases);
    128   *out_iphase_fn = iphase_fn;
    129 
    130 exit:
    131   return res;
    132 error:
    133   goto exit;
    134 }