commit 09d151c4f0ad685300f406442f1873685f187ce4
parent 3ef6e51e5d4f766d38ec245c8cdc74ae2d82e9b9
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 22 Jun 2022 11:57:19 +0200
Implement the rnsf_fetch_phase_fn function
Diffstat:
3 files changed, 164 insertions(+), 0 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -52,6 +52,7 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(RNSF_FILES_SRC
rnsf.c
+ rnsf_fetch_phase_fn.c
rnsf_log.c)
set(RNSF_FILES_INC
rnsf_c.h
diff --git a/src/rnsf.h b/src/rnsf.h
@@ -148,6 +148,32 @@ rnsf_phase_fn_get_discrete
(const struct rnsf_phase_fn* phase_fn,
struct rnsf_phase_fn_discrete* discrete);
+/* Fetch a phase function with respect to the submitted wavelength and the
+ * given sample. The phase function is chosen according to the alpha
+ * coefficient between [0, 1] defined as below:
+ *
+ * alpha = (lambda1 - wavelength) / (lambda1 - lambda0)
+ *
+ * with lambda0 < lambda1, the phase function wavelengths surrounding the
+ * submitted wavelength. The index returned refers to the phase function
+ * corresponding to lambda0 or lambda1 depending on whether the given sample is
+ * less than or greater than alpha, respectively.
+ *
+ * If the phase functions are stored by spectral band, lambda0 and lambda1 are
+ * the upper boundary of the lower band and the lower boundary of the upper
+ * band, respectively, which surround the wavelength. If the wavelength is in a
+ * specific band, the index for that band is simply returned.
+ *
+ * If the wavelength is outside the domain, the index returned refers to either
+ * the first or the last phase function depending on whether the wavelength is
+ * below or above the domain, respectively. */
+RNSF_API res_T
+rnsf_fetch_phase_fn
+ (const struct rnsf* rnsf,
+ const double wavelength,
+ const double sample, /* In [0, 1[ */
+ size_t* iphase_fn);
+
END_DECLS
#endif /* RNSF_H */
diff --git a/src/rnsf_fetch_phase_fn.c b/src/rnsf_fetch_phase_fn.c
@@ -0,0 +1,137 @@
+/* Copyright (C) 2022 GSMA - Université de Reims Champgne-Ardenne, CNRS
+ * Copyright (C) 2022 IPGP, Université Paris Cité, CNRS
+ * Copyright (C) 2022 LAPLACE - Université de Toulouse, CNRS, INPT, UPS
+ * Copyright (C) 2022 LATMOS/IPSL - UVSQ, Université Paris-Saclay,
+ * Sorbonne Université, CNRS
+ * Copyright (C) 2022 LESIA - Observatoire de Paris, Université PSL,
+ * Sorbonne Université, Université Paris Cité
+ * Copyright (C) 2022 LMD/IPSL - Sorbonne Université, Université PSL,
+ * Ecole Polytechnique, Institut Polytechnique de Paris,
+ * CNRS
+ * Copyright (C) 2022 |Meso|Star> (contact@meso-star.com)
+ *
+ * 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 "rnsf.h"
+#include "rnsf_c.h"
+#include "rnsf_log.h"
+#include "rnsf_phase_fn.h"
+
+#include <rsys/algorithm.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE int
+cmp_phase(const void* key, const void* item)
+{
+ const double wavelength = (*(const double*)key);
+ const struct rnsf_phase_fn* phase = item;
+ if(wavelength < phase->wlen[0]) {
+ return -1;
+ } else if(wavelength > phase->wlen[1]) {
+ return 1;
+ } else {
+ return 0;
+ }
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+rnsf_fetch_phase_fn
+ (const struct rnsf* rnsf,
+ const double wavelength,
+ const double sample, /* In [0, 1[ */
+ size_t* out_iphase_fn)
+{
+ const struct rnsf_phase_fn* phases = NULL;
+ const struct rnsf_phase_fn* phase = NULL;
+ size_t nphases = 0;
+ size_t iphase_fn = 0;
+ res_T res = RES_OK;
+
+ if(!rnsf
+ || wavelength < 0
+ || sample < 0
+ || sample >= 1
+ || !out_iphase_fn) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ phases = darray_phase_fn_cdata_get(&rnsf->phases);
+ nphases = darray_phase_fn_size_get(&rnsf->phases);
+
+ if(!nphases) {
+ log_warn(rnsf, "%s: no phase function to fetch.\n", FUNC_NAME);
+ *out_iphase_fn = 0;
+ goto exit;
+ }
+
+ phase = search_lower_bound
+ (&wavelength, phases, nphases, sizeof(*phases), cmp_phase);
+
+ /* All phase functions are set for a wavelength greater or equal to the
+ * subimitted wavelength */
+ if(phase == phases) {
+ ASSERT(phases[0].wlen[0] > wavelength
+ || (phases[0].wlen[0] <= wavelength && wavelength <= phases[0].wlen[1]));
+ iphase_fn = 0;
+ }
+
+ /* All phase functions are set for a wavelength less than the submitted
+ * wavelength */
+ else if(!phase) {
+ ASSERT(phases[nphases-1].wlen[1] < wavelength);
+ iphase_fn = nphases - 1;
+ }
+
+ /* The wavelength range of the found phase function overlaps the submitted
+ * wavelength */
+ else if(wavelength >= phase->wlen[0] && wavelength <= phase->wlen[1]) {
+ iphase_fn = (size_t)(phase - phases);
+ }
+
+ /* The found phase function is defined for a wavelength greater or equal to
+ * the submitted wavelength */
+ else {
+ /* Compute the linear interpolation parameter that defines the submitted
+ * wavelength wrt to the wavelenths boundaries of the phase functions
+ * surrounding it */
+ const struct rnsf_phase_fn* phase0 = phase-1;
+ const struct rnsf_phase_fn* phase1 = phase;
+ const double u =
+ (phase1->wlen[0] - wavelength)
+ / (phase1->wlen[0] - phase0->wlen[1]);
+
+ /* Consider the linear interplation parameter previously defined as a
+ * probability and use the submitted sample to select one of the phase
+ * function surrounding the submitted wavelength. */
+ if(sample < u) {
+ iphase_fn = (size_t)(phase - phases - 1);
+ } else {
+ iphase_fn = (size_t)(phase - phases);
+ }
+ }
+
+ ASSERT(iphase_fn < nphases);
+ *out_iphase_fn = iphase_fn;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
