htrdr

commit 07c2335981bd8dc9f57b41053ae2a980b5f829c4
parent 33d12b91fa68c4fb016b68c0b991e3df0a113c67
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Sun, 30 Mar 2025 12:39:59 +0200

core: correction for discrete wavelength distribution

It was failing when the number of wavelengths was 1.

The planets command source test has been enhanced to highlight this
problem, which has been corrected by this commit.

Diffstat:
M
src/core/htrdr_ran_wlen_discrete.c
 | 
83
++++++++++++++++++++++++++++++++++++++++++++++---------------------------------
M
src/planets/test_htrdr_planets_source.c
 | 
46
++++++++++++++++++++++++++++++++++++++++++++++

2 files changed, 94 insertions(+), 35 deletions(-)
diff --git a/src/core/htrdr_ran_wlen_discrete.c b/src/core/htrdr_ran_wlen_discrete.c
@@ -239,8 +239,10 @@ htrdr_ran_wlen_discrete_create
 
   res = setup_per_wlen_radiance(ran, args);
   if(res != RES_OK) goto error;
-  res = setup_distribution(ran, args);
-  if(res != RES_OK) goto error;
+  if(ran->nbands != 0) {
+    res = setup_distribution(ran, args);
+    if(res != RES_OK) goto error;
+  }
 
 exit:
   *out_ran = ran;
@@ -271,42 +273,53 @@ htrdr_ran_wlen_discrete_sample
    const double r1,
    double* pdf) /* In nm⁻¹ */
 {
-  double* find = NULL;
-  const double* proba = NULL;
-  const double* cumul = NULL;
-  const double* wlens = NULL;
-  double w0, w1; /* Wavelengths of the sampled band in nm */
-  double lambda; /* Sampled wavelength in nm */
-  size_t iband = 0;
-  double r0_next = nextafter(r0, DBL_MAX);
+  double lambda = 0;
+
+  ASSERT(ran);
   ASSERT(0 <= r0 && r0 < 1);
   ASSERT(0 <= r1 && r1 < 1);
 
-  cumul = darray_double_cdata_get(&ran->cumul);
-  proba = darray_double_cdata_get(&ran->proba);
-  wlens = darray_double_cdata_get(&ran->wlens);
-
-  /* Sample a band. Use r0_next rather than r0 to find the first entry that is
-   * not less than *or equal* to r0 */
-  find = search_lower_bound
-    (&r0_next, cumul, ran->nbands, sizeof(double), cmp_dbl);
-  ASSERT(find);
-
-  iband = (size_t)(find - cumul);
-  ASSERT(iband < ran->nbands);
-  ASSERT(cumul[iband] > r0 && (!iband || cumul[iband-1] <= r0));
-
-  /* Retrieve the boundaries of the sampled band */
-  w0 = wlens[iband+0];
-  w1 = wlens[iband+1];
-
-  /* Uniformely sample the wavelength in [w0, w1[ */
-  lambda = w0 + r1 * (w1 - w0);
-
-  if(pdf) {
-    const double pdf_wlen = 1.f / (w1-w0);
-    *pdf = proba[iband] * pdf_wlen;
+  if(ran->nbands == 0) {
+    ASSERT(darray_double_size_get(&ran->wlens) == 1);
+    lambda = darray_double_cdata_get(&ran->wlens)[0];
+    if(pdf) *pdf = 1;
+
+  } else {
+    double* find = NULL;
+    const double* proba = NULL;
+    const double* cumul = NULL;
+    const double* wlens = NULL;
+    double w0, w1; /* Wavelengths of the sampled band in nm */
+    size_t iband = 0;
+    double r0_next = nextafter(r0, DBL_MAX);
+    ASSERT(0 <= r0 && r0 < 1);
+    ASSERT(0 <= r1 && r1 < 1);
+
+    cumul = darray_double_cdata_get(&ran->cumul);
+    proba = darray_double_cdata_get(&ran->proba);
+    wlens = darray_double_cdata_get(&ran->wlens);
+
+    /* Sample a band. Use r0_next rather than r0 to find the first entry that is
+     * not less than *or equal* to r0 */
+    find = search_lower_bound
+      (&r0_next, cumul, ran->nbands, sizeof(double), cmp_dbl);
+    ASSERT(find);
+
+    iband = (size_t)(find - cumul);
+    ASSERT(iband < ran->nbands);
+    ASSERT(cumul[iband] > r0 && (!iband || cumul[iband-1] <= r0));
+
+    /* Retrieve the boundaries of the sampled band */
+    w0 = wlens[iband+0];
+    w1 = wlens[iband+1];
+
+    /* Uniformely sample the wavelength in [w0, w1[ */
+    lambda = w0 + r1 * (w1 - w0);
+
+    if(pdf) {
+      const double pdf_wlen = 1.f / (w1-w0);
+      *pdf = proba[iband] * pdf_wlen;
+    }
   }
-
   return lambda;
 }
diff --git a/src/planets/test_htrdr_planets_source.c b/src/planets/test_htrdr_planets_source.c
@@ -177,6 +177,51 @@ test_spectrum(struct htrdr* htrdr)
 }
 
 static void
+test_spectrum_1_entry(struct htrdr* htrdr)
+{
+  struct htrdr_ran_wlen_discrete_create_args distrib_args =
+    HTRDR_RAN_WLEN_DISCRETE_CREATE_ARGS_NULL;
+  struct htrdr_ran_wlen_discrete* distrib = NULL;
+
+  struct htrdr_planets_source_args source_args = HTRDR_PLANETS_SOURCE_ARGS_NULL;
+  struct htrdr_planets_source_spectrum spectrum;
+  struct htrdr_planets_source* source = NULL;
+
+  FILE* fp = NULL;
+  char rnrl_filename[] = "rnrl.bin";
+  double range[2] = {0,0};
+  double lambda = 0;
+  double pdf = 0;
+
+  CHK(fp = fopen(rnrl_filename, "w"));
+  write_per_wlen_radiances(fp, 4096, 1, 16, 16);
+  CHK(fclose(fp) == 0);
+
+  source_args.rnrl_filename = rnrl_filename;
+  source_args.longitude = 0;
+  source_args.latitude = 0;
+  source_args.distance = 0;
+  source_args.radius = 1;
+  source_args.temperature = -1;
+  CHK(htrdr_planets_source_create(htrdr, &source_args, &source) == RES_OK);
+
+  range[0] = range[1] = 0;
+  CHK(htrdr_planets_source_get_spectrum(source, range, &spectrum) == RES_OK);
+
+  distrib_args.get = htrdr_planets_source_spectrum_at;
+  distrib_args.nwavelengths = spectrum.size;
+  distrib_args.context = &spectrum;
+  CHK(htrdr_ran_wlen_discrete_create(htrdr, &distrib_args, &distrib) == RES_OK);
+
+  lambda = htrdr_ran_wlen_discrete_sample(distrib, 0.3, 0.5, &pdf);
+  CHK(lambda == 0);
+  CHK(pdf == 1);
+
+  htrdr_planets_source_ref_put(source);
+  htrdr_ran_wlen_discrete_ref_put(distrib);
+}
+
+static void
 test_spectrum_fail(struct htrdr* htrdr)
 {
   struct htrdr_planets_source_args source_args = HTRDR_PLANETS_SOURCE_ARGS_NULL;
@@ -291,6 +336,7 @@ main(int argc, char** argv)
     test_spectrum_from_files(htrdr, argc, argv);
   } else {
     test_spectrum(htrdr);
+    test_spectrum_1_entry(htrdr);
     test_spectrum_fail(htrdr);
   }