commit 585067cb887fb9a0ae78ae6acb8b8bea0b96afb7
parent f23ffee33ceae969dc30ed615942dbd219b1bb6a
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 2 Jun 2020 11:02:12 +0200
Rename htrdr_wlen_ran in htrdr_ran_wlen
Diffstat:
8 files changed, 433 insertions(+), 435 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -100,7 +100,7 @@ set(HTRDR_FILES_SRC
htrdr_interface.c
htrdr_main.c
htrdr_mtl.c
- htrdr_wlen_ran.c
+ htrdr_ran_wlen.c
htrdr_rectangle.c
htrdr_slab.c
htrdr_sun.c)
@@ -114,7 +114,7 @@ set(HTRDR_FILES_INC
htrdr_ground.h
htrdr_interface.h
htrdr_mtl.h
- htrdr_wlen_ran.h
+ htrdr_ran_wlen.h
htrdr_rectangle.h
htrdr_slab.h
htrdr_sun.h
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -24,7 +24,7 @@
#include "htrdr_camera.h"
#include "htrdr_ground.h"
#include "htrdr_mtl.h"
-#include "htrdr_wlen_ran.h"
+#include "htrdr_ran_wlen.h"
#include "htrdr_sun.h"
#include "htrdr_solve.h"
@@ -523,8 +523,8 @@ htrdr_init
ASSERT(htrdr->wlen_range_m[0] <= htrdr->wlen_range_m[1]);
n = (size_t)(args->wlen_lw_range[1] - args->wlen_lw_range[0]);
- res = htrdr_wlen_ran_create
- (htrdr, args->wlen_lw_range, n, Tref, &htrdr->wlen_ran);
+ res = htrdr_ran_wlen_create
+ (htrdr, args->wlen_lw_range, n, Tref, &htrdr->ran_wlen);
if(res != RES_OK) goto error;
} else {
@@ -536,8 +536,8 @@ htrdr_init
ASSERT(htrdr->wlen_range_m[0] <= htrdr->wlen_range_m[1]);
n = (size_t)(args->wlen_sw_range[1] - args->wlen_sw_range[0]);
- res = htrdr_wlen_ran_create
- (htrdr, args->wlen_sw_range, n, Tref, &htrdr->wlen_ran);
+ res = htrdr_ran_wlen_create
+ (htrdr, args->wlen_sw_range, n, Tref, &htrdr->ran_wlen);
if(res != RES_OK) goto error;
}
}
@@ -606,7 +606,7 @@ htrdr_release(struct htrdr* htrdr)
if(htrdr->buf) htrdr_buffer_ref_put(htrdr->buf);
if(htrdr->mtl) htrdr_mtl_ref_put(htrdr->mtl);
if(htrdr->cie) htrdr_cie_xyz_ref_put(htrdr->cie);
- if(htrdr->wlen_ran) htrdr_wlen_ran_ref_put(htrdr->wlen_ran);
+ if(htrdr->ran_wlen) htrdr_ran_wlen_ref_put(htrdr->ran_wlen);
if(htrdr->output && htrdr->output != stdout) fclose(htrdr->output);
if(htrdr->lifo_allocators) {
size_t i;
diff --git a/src/htrdr.h b/src/htrdr.h
@@ -50,7 +50,7 @@ struct htrdr {
struct htrdr_mtl* mtl;
struct htrdr_sun* sun;
struct htrdr_cie_xyz* cie;
- struct htrdr_wlen_ran* wlen_ran;
+ struct htrdr_ran_wlen* ran_wlen;
struct htrdr_camera* cam;
struct htrdr_buffer* buf;
diff --git a/src/htrdr_draw_radiance.c b/src/htrdr_draw_radiance.c
@@ -21,7 +21,7 @@
#include "htrdr_buffer.h"
#include "htrdr_camera.h"
#include "htrdr_cie_xyz.h"
-#include "htrdr_wlen_ran.h"
+#include "htrdr_ran_wlen.h"
#include "htrdr_solve.h"
#include <high_tune/htsky.h>
@@ -627,7 +627,7 @@ draw_pixel_integ
r2 = ssp_rng_canonical(rng);
/* Sample a wavelength */
- wlen = htrdr_wlen_ran_sample(htrdr->wlen_ran, r0, r1, &band_pdf);
+ wlen = htrdr_ran_wlen_sample(htrdr->ran_wlen, r0, r1, &band_pdf);
/* Select the associated band and sample a quadrature point */
iband = htsky_find_spectral_band(htrdr->sky, wlen);
diff --git a/src/htrdr_ran_wlen.c b/src/htrdr_ran_wlen.c
@@ -0,0 +1,366 @@
+/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
+ *
+ * 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/>. */
+
+#define _POSIX_C_SOURCE 200112L /* nextafter */
+
+#include "htrdr.h"
+#include "htrdr_c.h"
+#include "htrdr_ran_wlen.h"
+
+#include <high_tune/htsky.h>
+
+#include <rsys/algorithm.h>
+#include <rsys/cstr.h>
+#include <rsys/dynamic_array_double.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+#include <math.h> /* nextafter */
+
+struct htrdr_ran_wlen {
+ struct darray_double pdf;
+ struct darray_double cdf;
+ double range[2]; /* Boundaries of the spectral integration interval */
+ double band_len; /* Length in nanometers of a band */
+ double ref_temperature; /* In Kelvin */
+ size_t nbands; /* # bands */
+
+ struct htrdr* htrdr;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+setup_wlen_ran_cdf
+ (struct htrdr_ran_wlen* wlen_ran,
+ const char* func_name)
+{
+ double* pdf = NULL;
+ double* cdf = NULL;
+ double sum = 0;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(wlen_ran && func_name && wlen_ran->nbands != HTRDR_WLEN_RAN_CONTINUE);
+
+ res = darray_double_resize(&wlen_ran->cdf, wlen_ran->nbands);
+ if(res != RES_OK) {
+ htrdr_log_err(wlen_ran->htrdr,
+ "%s: Error allocating the CDF of the spectral bands -- %s.\n",
+ func_name, res_to_cstr(res));
+ goto error;
+ }
+ res = darray_double_resize(&wlen_ran->pdf, wlen_ran->nbands);
+ if(res != RES_OK) {
+ htrdr_log_err(wlen_ran->htrdr,
+ "%s: Error allocating the pdf of the spectral bands -- %s.\n",
+ func_name, res_to_cstr(res));
+ goto error;
+ }
+
+ cdf = darray_double_data_get(&wlen_ran->cdf);
+ pdf = darray_double_data_get(&wlen_ran->pdf); /* Now save pdf to correct weight */
+
+ htrdr_log(wlen_ran->htrdr,
+ "Number of bands of the spectrum cumulative: %lu\n",
+ (unsigned long)wlen_ran->nbands);
+
+ /* Compute the *unnormalized* probability to sample a small band */
+ FOR_EACH(i, 0, wlen_ran->nbands) {
+ double lambda_lo = wlen_ran->range[0] + (double)i * wlen_ran->band_len;
+ double lambda_hi = MMIN(lambda_lo + wlen_ran->band_len, wlen_ran->range[1]);
+ ASSERT(lambda_lo<= lambda_hi);
+ ASSERT(lambda_lo > wlen_ran->range[0] || eq_eps(lambda_lo, wlen_ran->range[0], 1.e-6));
+ ASSERT(lambda_lo < wlen_ran->range[1] || eq_eps(lambda_lo, wlen_ran->range[1], 1.e-6));
+
+ /* Convert from nanometer to meter */
+ lambda_lo *= 1.e-9;
+ lambda_hi *= 1.e-9;
+
+ /* Compute the probability of the current band */
+ pdf[i] = blackbody_fraction(lambda_lo, lambda_hi, wlen_ran->ref_temperature);
+
+ /* Update the norm */
+ sum += pdf[i];
+ }
+
+ /* Compute the cumulative of the previously computed probabilities */
+ FOR_EACH(i, 0, wlen_ran->nbands) {
+ /* Normalize the probability */
+ pdf[i] /= sum;
+
+ /* Setup the cumulative */
+ if(i == 0) {
+ cdf[i] = pdf[i];
+ } else {
+ cdf[i] = pdf[i] + cdf[i-1];
+ ASSERT(cdf[i] >= cdf[i-1]);
+ }
+ }
+ ASSERT(eq_eps(cdf[wlen_ran->nbands-1], 1, 1.e-6));
+ cdf[wlen_ran->nbands - 1] = 1.0; /* Handle numerical issue */
+
+exit:
+ return res;
+error:
+ darray_double_clear(&wlen_ran->cdf);
+ darray_double_clear(&wlen_ran->pdf);
+ goto exit;
+}
+
+static double
+wlen_ran_sample_continue
+ (const struct htrdr_ran_wlen* wlen_ran,
+ const double r,
+ const double range[2], /* In nanometer */
+ const char* func_name,
+ double* pdf)
+{
+ /* Numerical parameters of the dichotomy search */
+ const size_t MAX_ITER = 100;
+ const double EPSILON_LAMBDA_M = 1e-15;
+ const double EPSILON_BF = 1e-6;
+
+ /* Local variables */
+ double bf = 0;
+ double bf_prev = 0;
+ double bf_min_max = 0;
+ double lambda_m = 0;
+ double lambda_m_prev = 0;
+ double lambda_m_min = 0;
+ double lambda_m_max = 0;
+ double range_m[2] = {0, 0};
+ size_t i;
+
+ /* Check precondition */
+ ASSERT(wlen_ran && func_name);
+ ASSERT(range && range[0] < range[1]);
+ ASSERT(0 <= r && r < 1);
+
+ /* Convert the wavelength range in meters */
+ range_m[0] = range[0] * 1.0e-9;
+ range_m[1] = range[1] * 1.0e-9;
+
+ /* Setup the dichotomy search */
+ lambda_m_min = range_m[0];
+ lambda_m_max = range_m[1];
+ bf_min_max = blackbody_fraction
+ (range_m[0], range_m[1], wlen_ran->ref_temperature);
+
+ /* Numerically search the lambda corresponding to the submitted canonical
+ * number */
+ FOR_EACH(i, 0, MAX_ITER) {
+ double r_test;
+ lambda_m = (lambda_m_min + lambda_m_max) * 0.5;
+ bf = blackbody_fraction(range_m[0], lambda_m, wlen_ran->ref_temperature);
+
+ r_test = bf / bf_min_max;
+ if(r_test < r) {
+ lambda_m_min = lambda_m;
+ } else {
+ lambda_m_max = lambda_m;
+ }
+
+ if(fabs(lambda_m_prev - lambda_m) < EPSILON_LAMBDA_M
+ || fabs(bf_prev - bf) < EPSILON_BF)
+ break;
+
+ lambda_m_prev = lambda_m;
+ bf_prev = bf;
+ }
+ if(i >= MAX_ITER) {
+ htrdr_log_warn(wlen_ran->htrdr,
+ "%s: could not sample a wavelength in the range [%g, %g] nanometers "
+ "for the reference temperature %g Kelvin.\n",
+ func_name, SPLIT2(range), wlen_ran->ref_temperature);
+ }
+
+ if(pdf) {
+ const double Tref = wlen_ran->ref_temperature;
+ const double B_lambda = planck(lambda_m, lambda_m, Tref);
+ const double B_mean = planck(range_m[0], range_m[1], Tref);
+ *pdf = B_lambda / (B_mean * (range_m[1]-range_m[0]));
+ }
+
+ return lambda_m * 1.0e9; /* Convert in nanometers */
+}
+
+static double
+wlen_ran_sample_discrete
+ (const struct htrdr_ran_wlen* wlen_ran,
+ const double r0,
+ const double r1,
+ const char* func_name,
+ double* pdf)
+{
+ const double* cdf = NULL;
+ const double* find = NULL;
+ double r0_next = nextafter(r0, DBL_MAX);
+ double band_range[2];
+ double lambda = 0;
+ double pdf_continue = 0;
+ double pdf_band = 0;
+ size_t cdf_length = 0;
+ size_t i;
+ ASSERT(wlen_ran && wlen_ran->nbands != HTRDR_WLEN_RAN_CONTINUE);
+ ASSERT(0 <= r0 && r0 < 1);
+ ASSERT(0 <= r1 && r1 < 1);
+ (void)func_name;
+ (void)pdf_band;
+
+ cdf = darray_double_cdata_get(&wlen_ran->cdf);
+ cdf_length = darray_double_size_get(&wlen_ran->cdf);
+ ASSERT(cdf_length > 0);
+
+ /* Use r_next rather than r0 in order to find the first entry that is not less
+ * than *or equal* to r0 */
+ find = search_lower_bound(&r0_next, cdf, cdf_length, sizeof(double), cmp_dbl);
+ ASSERT(find);
+
+ i = (size_t)(find - cdf);
+ ASSERT(i < cdf_length && cdf[i] > r0 && (!i || cdf[i-1] <= r0));
+
+ band_range[0] = wlen_ran->range[0] + (double)i*wlen_ran->band_len;
+ band_range[1] = band_range[0] + wlen_ran->band_len;
+
+ /* Fetch the pdf of the sampled band */
+ pdf_band = darray_double_cdata_get(&wlen_ran->pdf)[i];
+
+ /* Uniformly sample a wavelength in the sampled band */
+ lambda = band_range[0] + (band_range[1] - band_range[0]) * r1;
+ pdf_continue = 1.0 / ((band_range[1] - band_range[0])*1.e-9);
+
+ if(pdf) {
+ *pdf = pdf_band * pdf_continue;
+ }
+
+ return lambda;
+}
+
+static void
+release_wlen_ran(ref_T* ref)
+{
+ struct htrdr_ran_wlen* wlen_ran = NULL;
+ ASSERT(ref);
+ wlen_ran = CONTAINER_OF(ref, struct htrdr_ran_wlen, ref);
+ darray_double_release(&wlen_ran->cdf);
+ darray_double_release(&wlen_ran->pdf);
+ MEM_RM(wlen_ran->htrdr->allocator, wlen_ran);
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+htrdr_ran_wlen_create
+ (struct htrdr* htrdr,
+ /* range must be included in [200,1000] nm for solar or in [1000, 100000]
+ * nanometers for longwave (thermal)*/
+ const double range[2],
+ const size_t nbands, /* # bands used to discretized CDF */
+ const double ref_temperature,
+ struct htrdr_ran_wlen** out_wlen_ran)
+{
+ struct htrdr_ran_wlen* wlen_ran = NULL;
+ double min_band_len = 0;
+ res_T res = RES_OK;
+ ASSERT(htrdr && range && out_wlen_ran && ref_temperature > 0);
+ ASSERT(ref_temperature > 0);
+ ASSERT(range[0] >= 200);
+ ASSERT(range[1] <= 100000);
+ ASSERT(range[0] <= range[1]);
+
+ wlen_ran = MEM_CALLOC(htrdr->allocator, 1, sizeof(*wlen_ran));
+ if(!wlen_ran) {
+ res = RES_MEM_ERR;
+ htrdr_log_err(htrdr,
+ "%s: could not allocate longwave random variate data structure -- %s.\n",
+ FUNC_NAME, res_to_cstr(res));
+ goto error;
+ }
+ ref_init(&wlen_ran->ref);
+ wlen_ran->htrdr = htrdr;
+ darray_double_init(htrdr->allocator, &wlen_ran->cdf);
+ darray_double_init(htrdr->allocator, &wlen_ran->pdf);
+
+ wlen_ran->range[0] = range[0];
+ wlen_ran->range[1] = range[1];
+ wlen_ran->ref_temperature = ref_temperature;
+ wlen_ran->nbands = nbands;
+
+ min_band_len = compute_sky_min_band_len(wlen_ran->htrdr->sky, wlen_ran->range);
+
+ if(nbands == HTRDR_WLEN_RAN_CONTINUE) {
+ wlen_ran->band_len = 0;
+ } else {
+ wlen_ran->band_len = (range[1] - range[0]) / (double)wlen_ran->nbands;
+
+ /* Adjust the band length to ensure that each sky spectral interval is
+ * overlapped by at least one band */
+ if(wlen_ran->band_len > min_band_len) {
+ wlen_ran->band_len = min_band_len;
+ wlen_ran->nbands = (size_t)ceil((range[1] - range[0]) / wlen_ran->band_len);
+ }
+
+ res = setup_wlen_ran_cdf(wlen_ran, FUNC_NAME);
+ if(res != RES_OK) goto error;
+ }
+
+ htrdr_log(htrdr, "Spectral interval defined on [%g, %g] nanometers.\n",
+ range[0], range[1]);
+
+exit:
+ *out_wlen_ran = wlen_ran;
+ return res;
+error:
+ if(wlen_ran) htrdr_ran_wlen_ref_put(wlen_ran);
+ goto exit;
+}
+
+void
+htrdr_ran_wlen_ref_get(struct htrdr_ran_wlen* wlen_ran)
+{
+ ASSERT(wlen_ran);
+ ref_get(&wlen_ran->ref);
+}
+
+void
+htrdr_ran_wlen_ref_put(struct htrdr_ran_wlen* wlen_ran)
+{
+ ASSERT(wlen_ran);
+ ref_put(&wlen_ran->ref, release_wlen_ran);
+}
+
+double
+htrdr_ran_wlen_sample
+ (const struct htrdr_ran_wlen* wlen_ran,
+ const double r0,
+ const double r1,
+ double* pdf)
+{
+ ASSERT(wlen_ran);
+ if(wlen_ran->nbands != HTRDR_WLEN_RAN_CONTINUE) { /* Discrete */
+ return wlen_ran_sample_discrete(wlen_ran, r0, r1, FUNC_NAME, pdf);
+ } else if(eq_eps(wlen_ran->range[0], wlen_ran->range[1], 1.e-6)) {
+ if(pdf) *pdf = 1;
+ return wlen_ran->range[0];
+ } else { /* Continue */
+ return wlen_ran_sample_continue
+ (wlen_ran, r0, wlen_ran->range, FUNC_NAME, pdf);
+ }
+}
+
diff --git a/src/htrdr_ran_wlen.h b/src/htrdr_ran_wlen.h
@@ -0,0 +1,56 @@
+/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
+ *
+ * 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/>. */
+
+#ifndef HTRDR_RAN_WLEN_H
+#define HTRDR_RAN_WLEN_H
+
+#include <rsys/rsys.h>
+
+#define HTRDR_WLEN_RAN_CONTINUE 0
+
+struct htrdr;
+struct htrdr_ran_wlen;
+
+extern LOCAL_SYM res_T
+htrdr_ran_wlen_create
+ (struct htrdr* htrdr,
+ /* range must be included in [200,1000] nm for solar or in [1000, 100000]
+ * nanometers for longwave (thermal)*/
+ const double range[2],
+ /* # bands used to discretisze the spectral domain. HTRDR_WLEN_RAN_CONTINUE
+ * <=> no discretisation */
+ const size_t nbands, /* Hint on #bands used to discretised th CDF */
+ const double ref_temperature, /* Reference temperature */
+ struct htrdr_ran_wlen** wlen_ran);
+
+extern LOCAL_SYM void
+htrdr_ran_wlen_ref_get
+ (struct htrdr_ran_wlen* wlen_ran);
+
+extern LOCAL_SYM void
+htrdr_ran_wlen_ref_put
+ (struct htrdr_ran_wlen* wlen_ran);
+
+/* Return a wavelength in nanometer */
+extern LOCAL_SYM double
+htrdr_ran_wlen_sample
+ (const struct htrdr_ran_wlen* wlen_ran,
+ const double r0, /* Canonical number in [0, 1[ */
+ const double r1, /* Canonical number in [0, 1[ */
+ double* pdf); /* May be NULL */
+
+#endif /* HTRDR_RAN_WLEN_H */
+
diff --git a/src/htrdr_wlen_ran.c b/src/htrdr_wlen_ran.c
@@ -1,365 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020 |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/>. */
-
-#define _POSIX_C_SOURCE 200112L /* nextafter */
-
-#include "htrdr.h"
-#include "htrdr_c.h"
-#include "htrdr_wlen_ran.h"
-
-#include <high_tune/htsky.h>
-
-#include <rsys/algorithm.h>
-#include <rsys/cstr.h>
-#include <rsys/dynamic_array_double.h>
-#include <rsys/mem_allocator.h>
-#include <rsys/ref_count.h>
-
-#include <math.h> /* nextafter */
-
-struct htrdr_wlen_ran {
- struct darray_double pdf;
- struct darray_double cdf;
- double range[2]; /* Boundaries of the spectral integration interval */
- double band_len; /* Length in nanometers of a band */
- double ref_temperature; /* In Kelvin */
- size_t nbands; /* # bands */
-
- struct htrdr* htrdr;
- ref_T ref;
-};
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static res_T
-setup_wlen_ran_cdf
- (struct htrdr_wlen_ran* wlen_ran,
- const char* func_name)
-{
- double* pdf = NULL;
- double* cdf = NULL;
- double sum = 0;
- size_t i;
- res_T res = RES_OK;
- ASSERT(wlen_ran && func_name && wlen_ran->nbands != HTRDR_WLEN_RAN_CONTINUE);
-
- res = darray_double_resize(&wlen_ran->cdf, wlen_ran->nbands);
- if(res != RES_OK) {
- htrdr_log_err(wlen_ran->htrdr,
- "%s: Error allocating the CDF of the spectral bands -- %s.\n",
- func_name, res_to_cstr(res));
- goto error;
- }
- res = darray_double_resize(&wlen_ran->pdf, wlen_ran->nbands);
- if(res != RES_OK) {
- htrdr_log_err(wlen_ran->htrdr,
- "%s: Error allocating the pdf of the spectral bands -- %s.\n",
- func_name, res_to_cstr(res));
- goto error;
- }
-
- cdf = darray_double_data_get(&wlen_ran->cdf);
- pdf = darray_double_data_get(&wlen_ran->pdf); /* Now save pdf to correct weight */
-
- htrdr_log(wlen_ran->htrdr,
- "Number of bands of the spectrum cumulative: %lu\n",
- (unsigned long)wlen_ran->nbands);
-
- /* Compute the *unnormalized* probability to sample a small band */
- FOR_EACH(i, 0, wlen_ran->nbands) {
- double lambda_lo = wlen_ran->range[0] + (double)i * wlen_ran->band_len;
- double lambda_hi = MMIN(lambda_lo + wlen_ran->band_len, wlen_ran->range[1]);
- ASSERT(lambda_lo<= lambda_hi);
- ASSERT(lambda_lo > wlen_ran->range[0] || eq_eps(lambda_lo, wlen_ran->range[0], 1.e-6));
- ASSERT(lambda_lo < wlen_ran->range[1] || eq_eps(lambda_lo, wlen_ran->range[1], 1.e-6));
-
- /* Convert from nanometer to meter */
- lambda_lo *= 1.e-9;
- lambda_hi *= 1.e-9;
-
- /* Compute the probability of the current band */
- pdf[i] = blackbody_fraction(lambda_lo, lambda_hi, wlen_ran->ref_temperature);
-
- /* Update the norm */
- sum += pdf[i];
- }
-
- /* Compute the cumulative of the previously computed probabilities */
- FOR_EACH(i, 0, wlen_ran->nbands) {
- /* Normalize the probability */
- pdf[i] /= sum;
-
- /* Setup the cumulative */
- if(i == 0) {
- cdf[i] = pdf[i];
- } else {
- cdf[i] = pdf[i] + cdf[i-1];
- ASSERT(cdf[i] >= cdf[i-1]);
- }
- }
- ASSERT(eq_eps(cdf[wlen_ran->nbands-1], 1, 1.e-6));
- cdf[wlen_ran->nbands - 1] = 1.0; /* Handle numerical issue */
-
-exit:
- return res;
-error:
- darray_double_clear(&wlen_ran->cdf);
- darray_double_clear(&wlen_ran->pdf);
- goto exit;
-}
-
-static double
-wlen_ran_sample_continue
- (const struct htrdr_wlen_ran* wlen_ran,
- const double r,
- const double range[2], /* In nanometer */
- const char* func_name,
- double* pdf)
-{
- /* Numerical parameters of the dichotomy search */
- const size_t MAX_ITER = 100;
- const double EPSILON_LAMBDA_M = 1e-15;
- const double EPSILON_BF = 1e-6;
-
- /* Local variables */
- double bf = 0;
- double bf_prev = 0;
- double bf_min_max = 0;
- double lambda_m = 0;
- double lambda_m_prev = 0;
- double lambda_m_min = 0;
- double lambda_m_max = 0;
- double range_m[2] = {0, 0};
- size_t i;
-
- /* Check precondition */
- ASSERT(wlen_ran && func_name);
- ASSERT(range && range[0] < range[1]);
- ASSERT(0 <= r && r < 1);
-
- /* Convert the wavelength range in meters */
- range_m[0] = range[0] * 1.0e-9;
- range_m[1] = range[1] * 1.0e-9;
-
- /* Setup the dichotomy search */
- lambda_m_min = range_m[0];
- lambda_m_max = range_m[1];
- bf_min_max = blackbody_fraction
- (range_m[0], range_m[1], wlen_ran->ref_temperature);
-
- /* Numerically search the lambda corresponding to the submitted canonical
- * number */
- FOR_EACH(i, 0, MAX_ITER) {
- double r_test;
- lambda_m = (lambda_m_min + lambda_m_max) * 0.5;
- bf = blackbody_fraction(range_m[0], lambda_m, wlen_ran->ref_temperature);
-
- r_test = bf / bf_min_max;
- if(r_test < r) {
- lambda_m_min = lambda_m;
- } else {
- lambda_m_max = lambda_m;
- }
-
- if(fabs(lambda_m_prev - lambda_m) < EPSILON_LAMBDA_M
- || fabs(bf_prev - bf) < EPSILON_BF)
- break;
-
- lambda_m_prev = lambda_m;
- bf_prev = bf;
- }
- if(i >= MAX_ITER) {
- htrdr_log_warn(wlen_ran->htrdr,
- "%s: could not sample a wavelength in the range [%g, %g] nanometers "
- "for the reference temperature %g Kelvin.\n",
- func_name, SPLIT2(range), wlen_ran->ref_temperature);
- }
-
- if(pdf) {
- const double Tref = wlen_ran->ref_temperature;
- const double B_lambda = planck(lambda_m, lambda_m, Tref);
- const double B_mean = planck(range_m[0], range_m[1], Tref);
- *pdf = B_lambda / (B_mean * (range_m[1]-range_m[0]));
- }
-
- return lambda_m * 1.0e9; /* Convert in nanometers */
-}
-
-static double
-wlen_ran_sample_discrete
- (const struct htrdr_wlen_ran* wlen_ran,
- const double r0,
- const double r1,
- const char* func_name,
- double* pdf)
-{
- const double* cdf = NULL;
- const double* find = NULL;
- double r0_next = nextafter(r0, DBL_MAX);
- double band_range[2];
- double lambda = 0;
- double pdf_continue = 0;
- double pdf_band = 0;
- size_t cdf_length = 0;
- size_t i;
- ASSERT(wlen_ran && wlen_ran->nbands != HTRDR_WLEN_RAN_CONTINUE);
- ASSERT(0 <= r0 && r0 < 1);
- ASSERT(0 <= r1 && r1 < 1);
- (void)func_name;
- (void)pdf_band;
-
- cdf = darray_double_cdata_get(&wlen_ran->cdf);
- cdf_length = darray_double_size_get(&wlen_ran->cdf);
- ASSERT(cdf_length > 0);
-
- /* Use r_next rather than r0 in order to find the first entry that is not less
- * than *or equal* to r0 */
- find = search_lower_bound(&r0_next, cdf, cdf_length, sizeof(double), cmp_dbl);
- ASSERT(find);
-
- i = (size_t)(find - cdf);
- ASSERT(i < cdf_length && cdf[i] > r0 && (!i || cdf[i-1] <= r0));
-
- band_range[0] = wlen_ran->range[0] + (double)i*wlen_ran->band_len;
- band_range[1] = band_range[0] + wlen_ran->band_len;
-
- /* Fetch the pdf of the sampled band */
- pdf_band = darray_double_cdata_get(&wlen_ran->pdf)[i];
-
- /* Uniformly sample a wavelength in the sampled band */
- lambda = band_range[0] + (band_range[1] - band_range[0]) * r1;
- pdf_continue = 1.0 / ((band_range[1] - band_range[0])*1.e-9);
-
- if(pdf) {
- *pdf = pdf_band * pdf_continue;
- }
-
- return lambda;
-}
-
-static void
-release_wlen_ran(ref_T* ref)
-{
- struct htrdr_wlen_ran* wlen_ran = NULL;
- ASSERT(ref);
- wlen_ran = CONTAINER_OF(ref, struct htrdr_wlen_ran, ref);
- darray_double_release(&wlen_ran->cdf);
- darray_double_release(&wlen_ran->pdf);
- MEM_RM(wlen_ran->htrdr->allocator, wlen_ran);
-}
-
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-res_T
-htrdr_wlen_ran_create
- (struct htrdr* htrdr,
- /* range must be included in [200,1000] nm for solar or in [1000, 100000]
- * nanometers for longwave (thermal)*/
- const double range[2],
- const size_t nbands, /* # bands used to discretized CDF */
- const double ref_temperature,
- struct htrdr_wlen_ran** out_wlen_ran)
-{
- struct htrdr_wlen_ran* wlen_ran = NULL;
- double min_band_len = 0;
- res_T res = RES_OK;
- ASSERT(htrdr && range && out_wlen_ran && ref_temperature > 0);
- ASSERT(ref_temperature > 0);
- ASSERT(range[0] >= 200);
- ASSERT(range[1] <= 100000);
- ASSERT(range[0] <= range[1]);
-
- wlen_ran = MEM_CALLOC(htrdr->allocator, 1, sizeof(*wlen_ran));
- if(!wlen_ran) {
- res = RES_MEM_ERR;
- htrdr_log_err(htrdr,
- "%s: could not allocate longwave random variate data structure -- %s.\n",
- FUNC_NAME, res_to_cstr(res));
- goto error;
- }
- ref_init(&wlen_ran->ref);
- wlen_ran->htrdr = htrdr;
- darray_double_init(htrdr->allocator, &wlen_ran->cdf);
- darray_double_init(htrdr->allocator, &wlen_ran->pdf);
-
- wlen_ran->range[0] = range[0];
- wlen_ran->range[1] = range[1];
- wlen_ran->ref_temperature = ref_temperature;
- wlen_ran->nbands = nbands;
-
- min_band_len = compute_sky_min_band_len(wlen_ran->htrdr->sky, wlen_ran->range);
-
- if(nbands == HTRDR_WLEN_RAN_CONTINUE) {
- wlen_ran->band_len = 0;
- } else {
- wlen_ran->band_len = (range[1] - range[0]) / (double)wlen_ran->nbands;
-
- /* Adjust the band length to ensure that each sky spectral interval is
- * overlapped by at least one band */
- if(wlen_ran->band_len > min_band_len) {
- wlen_ran->band_len = min_band_len;
- wlen_ran->nbands = (size_t)ceil((range[1] - range[0]) / wlen_ran->band_len);
- }
-
- res = setup_wlen_ran_cdf(wlen_ran, FUNC_NAME);
- if(res != RES_OK) goto error;
- }
-
- htrdr_log(htrdr, "Spectral interval defined on [%g, %g] nanometers.\n",
- range[0], range[1]);
-
-exit:
- *out_wlen_ran = wlen_ran;
- return res;
-error:
- if(wlen_ran) htrdr_wlen_ran_ref_put(wlen_ran);
- goto exit;
-}
-
-void
-htrdr_wlen_ran_ref_get(struct htrdr_wlen_ran* wlen_ran)
-{
- ASSERT(wlen_ran);
- ref_get(&wlen_ran->ref);
-}
-
-void
-htrdr_wlen_ran_ref_put(struct htrdr_wlen_ran* wlen_ran)
-{
- ASSERT(wlen_ran);
- ref_put(&wlen_ran->ref, release_wlen_ran);
-}
-
-double
-htrdr_wlen_ran_sample
- (const struct htrdr_wlen_ran* wlen_ran,
- const double r0,
- const double r1,
- double* pdf)
-{
- ASSERT(wlen_ran);
- if(wlen_ran->nbands != HTRDR_WLEN_RAN_CONTINUE) { /* Discrete */
- return wlen_ran_sample_discrete(wlen_ran, r0, r1, FUNC_NAME, pdf);
- } else if(eq_eps(wlen_ran->range[0], wlen_ran->range[1], 1.e-6)) {
- if(pdf) *pdf = 1;
- return wlen_ran->range[0];
- } else { /* Continue */
- return wlen_ran_sample_continue
- (wlen_ran, r0, wlen_ran->range, FUNC_NAME, pdf);
- }
-}
-
diff --git a/src/htrdr_wlen_ran.h b/src/htrdr_wlen_ran.h
@@ -1,59 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020 |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/>. */
-
-#ifndef HTRDR_WLEN_RAN_H
-#define HTRDR_WLEN_RAN_H
-
-#include <rsys/rsys.h>
-
-#define HTRDR_WLEN_RAN_CONTINUE 0
-#define HTRDR_WLEN_RAN_SOLAR_WVN_MIN 820 # 12195 nm
-#define HTRDR_WLEN_RAN_SOLAR_WVN_MAX 50000 # 200 nm
-#define HTRDR_WLEN_RAN_THERMAL_WVN_MIN 10 # 1000000 nm
-#define HTRDR_WLEN_RAN_THERMAL_WVN_MAX 3250 # 3077 nm
-
-struct htrdr;
-struct htrdr_wlen_ran;
-
-extern LOCAL_SYM res_T
-htrdr_wlen_ran_create
- (struct htrdr* htrdr,
- /* range must be included in [200,1000] nm for solar or in [1000, 100000]
- * nanometers for longwave (thermal)*/
- const double range[2],
- /* # bands used to discretisze the spectral domain. HTRDR_WLEN_RAN_CONTINUE
- * <=> no discretisation */
- const size_t nbands, /* Hint on #bands used to discretised th CDF */
- const double ref_temperature, /* Reference temperature */
- struct htrdr_wlen_ran** wlen_ran);
-
-extern LOCAL_SYM void
-htrdr_wlen_ran_ref_get
- (struct htrdr_wlen_ran* wlen_ran);
-
-extern LOCAL_SYM void
-htrdr_wlen_ran_ref_put
- (struct htrdr_wlen_ran* wlen_ran);
-
-/* Return a wavelength in nanometer */
-extern LOCAL_SYM double
-htrdr_wlen_ran_sample
- (const struct htrdr_wlen_ran* wlen_ran,
- const double r0, /* Canonical number in [0, 1[ */
- const double r1, /* Canonical number in [0, 1[ */
- double* pdf); /* May be NULL */
-
-#endif /* HTRDR_WLEN_RAN_H */
-
