star-blackbody

commit df2e9538aa4221652604881eb58ab02bb388ad8c
parent 640332c1eaced62e981063affc9edcaa84ecd84c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 29 Nov 2023 14:15:47 +0100

Fix the sbb_ran_planck_sample function

It would not sample anything if the spectral interval was less than one
micron. This problem is inherited from the old unit used to define
wavelengths: we then used nanometers, whereas the function now expects
meters. So, while we were comparing the limits of the interval with an
epsilon of 1.e-6 nm, the same epsilon was used, but this time defined in
meters, which was far too large. The epsilon is now set at 1 picometre.

Note that the Planck sampling test has been updated to show this
problem, which this commit now corrects.

Diffstat:
M
src/sbb_ran_planck.c
 | 
2
+-
M
src/test_sbb_ran_planck.c
 | 
23
++++++++++++++++++-----

2 files changed, 19 insertions(+), 6 deletions(-)
diff --git a/src/sbb_ran_planck.c b/src/sbb_ran_planck.c
@@ -372,7 +372,7 @@ sbb_ran_planck_sample
 {
   ASSERT(planck && r0 >= 0 && r1 >= 0 && r0 < 1 && r1 < 1);
 
-  if(eq_eps(planck->range[0], planck->range[1], 1.e-6)) {
+  if(eq_eps(planck->range[0], planck->range[1], 1.e-12 /* <=> 1 pm */)) {
     if(pdf) *pdf = 1;
     return planck->range[0];
 
diff --git a/src/test_sbb_ran_planck.c b/src/test_sbb_ran_planck.c
@@ -107,8 +107,10 @@ planck_integration
 
   ref = sbb_planck(lambda_lo, lambda_hi, T0) * delta_lambda;
 
-  printf("planck(%g m, %g m, %g K) = %g ~ %g +/- %g [W/m^2/sr] (Tref = %g K)\n",
-    lambda_lo, lambda_hi, T0, ref, E, SE, Tref);
+  printf(
+    "planck(%g m, %g m, %g K) = %g ~ %g +/- %g [W/m^2/sr] "
+    "(Tref = %g K; nbands = %lu)\n",
+    lambda_lo, lambda_hi, T0, ref, E, SE, Tref, (unsigned long)nbands);
   CHK(eq_eps(ref, E, 3*SE));
 
   CHK(sbb_ran_planck_ref_put(planck) == RES_OK);
@@ -118,14 +120,17 @@ int
 main(int argc, char** argv)
 {
   /* Input parameters */
-  const double lambda_lo = 3.6e-6; /* [m] */
-  const double lambda_hi = 12.3e-6; /* [m] */
-  const double delta_lambda = lambda_hi - lambda_lo;
+  double lambda_lo = 0;
+  double lambda_hi = 0;
+  double delta_lambda = 0;
   size_t nbands = 0;
   (void)argc, (void)argv;
 
   check_planck();
 
+  lambda_lo = 3.6e-6; /* [m] */
+  lambda_hi = 12.3e-6; /* [m] */
+  delta_lambda = lambda_hi - lambda_lo;
   nbands = (size_t)(delta_lambda * 1e9);
 
   planck_integration(lambda_lo, lambda_hi, 500, 550, 0);
@@ -133,5 +138,13 @@ main(int argc, char** argv)
   planck_integration(lambda_lo, lambda_hi, 300, 550, 0);
   planck_integration(lambda_lo, lambda_hi, 300, 550, nbands);
   planck_integration(lambda_lo, lambda_hi, 300, 550, nbands/4);
+
+  lambda_lo = 3.7e-6; /* [m] */
+  lambda_hi = 4.5e-6; /* [m] */
+  delta_lambda = lambda_hi - lambda_lo;
+  nbands = (size_t)(delta_lambda * 1e9);
+
+  planck_integration(lambda_lo, lambda_hi, 300, 550,  0);
+  planck_integration(lambda_lo, lambda_hi, 300, 550,  nbands);
   return 0;
 }