htrdr

commit 59a22ee1441fd8435237d3b15badd9e507906389
parent 10e16f52ce11eb43cc4c7a81c2042aa328c37f3c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 21 Apr 2020 14:26:04 +0200

Add support of a monochromatic LW

Diffstat:
M
src/htrdr.c
 | 
8
++++----
M
src/htrdr_c.h
 | 
55
++++++++++++++++++++++++++++++++++++++++++++++---------

2 files changed, 50 insertions(+), 13 deletions(-)
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -489,10 +489,6 @@ htrdr_init
   res = htsky_create(&htrdr->logger, htrdr->allocator, &htsky_args, &htrdr->sky);
   if(res != RES_OK) goto error;
 
-  htrdr->wlen_range_m[0] = args->wlen_lw_range[0]*1e-9; /* Convert in meters */
-  htrdr->wlen_range_m[1] = args->wlen_lw_range[1]*1e-9; /* Convert in meters */
-  ASSERT(htrdr->wlen_range_m[0] <= htrdr->wlen_range_m[1]);
-
   if(!htsky_is_long_wave(htrdr->sky)) { /* Short wave random variate */
     const double* range = HTRDR_CIE_XYZ_RANGE_DEFAULT;
     size_t n;
@@ -505,6 +501,10 @@ htrdr_init
     const double Tref = 290; /* In Kelvin */
     size_t n;
 
+    htrdr->wlen_range_m[0] = args->wlen_lw_range[0]*1e-9; /* Convert in meters */
+    htrdr->wlen_range_m[1] = args->wlen_lw_range[1]*1e-9; /* Convert in meters */
+    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_ran_lw_create
       (htrdr, args->wlen_lw_range, n, Tref, &htrdr->ran_lw);
diff --git a/src/htrdr_c.h b/src/htrdr_c.h
@@ -127,8 +127,8 @@ wiebelt(const double v)
 
 static INLINE double
 blackbody_fraction
-  (const double lambda0, /* In meter */
-   const double lambda1, /* In meter */
+  (const double lambda0, /* In meters */
+   const double lambda1, /* In meters */
    const double temperature) /* In Kelvin */
 {
   const double C2 = 1.43877735e-2; /* m.K */
@@ -141,10 +141,31 @@ blackbody_fraction
   return w1 - w0;
 }
 
+
+
+/* Return the Planck value in W/m^2/sr/m at a given wavelength */
 static INLINE double
-planck
-  (const double lambda_min, /* In meter */
-   const double lambda_max, /* In meter */
+planck_monochromatic
+  (const double lambda, /* In meters */
+   const double temperature) /* In Kelvin */
+{
+  const double c = 299792458; /* m/s */
+  const double h = 6.62607015e-34; /* J.s */
+  const double k = 1.380649e-23; /* J/K */
+  const double lambda2 = lambda*lambda;
+  const double lambda5 = lambda2*lambda2*lambda;
+  const double B = ((2.0 * h * c*c) / lambda5) /* W/m^2/sr/m */
+                 / (exp(h*c/(lambda*k*temperature))-1.0);
+  ASSERT(temperature > 0);
+  return B;
+}
+
+/* Return the average Planck value in W/m^2/sr/m over the
+ * [lambda_min, lambda_max] interval. */
+static INLINE double
+planck_interval
+  (const double lambda_min, /* In meters */
+   const double lambda_max, /* In meters */
    const double temperature) /* In Kelvin  */
 {
   const double T2 = temperature*temperature;
@@ -152,15 +173,31 @@ planck
   const double BOLTZMANN_CONSTANT = 5.6696e-8; /* W/m^2/K^4 */
   ASSERT(lambda_min < lambda_max && temperature > 0);
   return blackbody_fraction(lambda_min, lambda_max, temperature)
-       * BOLTZMANN_CONSTANT * T4 / PI; /* In W/m^2/sr */
+       * BOLTZMANN_CONSTANT * T4 / (PI * (lambda_max-lambda_min)); /* In W/m^2/sr/m */
+}
+
+/* Invoke planck_monochromatic or planck_interval whether the submitted
+ * interval is null or not, respectively. The returned value is in W/m^2/sr/m */
+static FINLINE double
+planck
+  (const double lambda_min, /* In meters */
+   const double lambda_max, /* In meters */
+   const double temperature) /* In Kelvin  */
+{
+  ASSERT(lambda_min <= lambda_max && temperature > 0);
+  if(lambda_min == lambda_max) {
+    return planck_monochromatic(lambda_min, temperature);
+  } else {
+    return planck_interval(lambda_min, lambda_max, temperature);
+  }
 }
 
 extern LOCAL_SYM res_T
 brightness_temperature
   (struct htrdr* htrdr,
-   const double lambda_min, /* In meter */
-   const double lambda_max, /* In meter */
-   /* Integrated over [lambda_min, lambda_max]. In W/m^2/sr */
+   const double lambda_min, /* In meters */
+   const double lambda_max, /* In meters */
+   /* Integrated over [lambda_min, lambda_max]. In W/m^2/sr/m */
    const double radiance,
    double* temperature);