commit 14f76f7e274d6364db4fb5daeb5c72ff36522f35
parent eb83e35c47e11f837b72594f6033fc208c195882
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 4 Jun 2020 15:21:43 +0200
Add the Planck reference temperature as a command line option
Diffstat:
9 files changed, 255 insertions(+), 185 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -103,6 +103,7 @@ set(HTRDR_FILES_SRC
htrdr_ran_wlen.c
htrdr_rectangle.c
htrdr_slab.c
+ htrdr_spectral.c
htrdr_sun.c)
set(HTRDR_FILES_INC
htrdr.h
@@ -117,6 +118,7 @@ set(HTRDR_FILES_INC
htrdr_ran_wlen.h
htrdr_rectangle.h
htrdr_slab.h
+ htrdr_spectral.h
htrdr_sun.h
htrdr_solve.h)
set(HTRDR_FILES_INC2 # Generated files
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -428,6 +428,7 @@ htrdr_init
htrdr->grid_max_definition[1] = args->grid_max_definition[1];
htrdr->grid_max_definition[2] = args->grid_max_definition[2];
htrdr->spectral_type = args->spectral_type;
+ htrdr->ref_temperature = args->ref_temperature;
res = init_mpi(htrdr);
if(res != RES_OK) goto error;
@@ -524,20 +525,20 @@ htrdr_init
if(res != RES_OK) goto error;
} else {
- double Tref; /* In Kelvin */
size_t n;
- switch(htrdr->spectral_type) {
- case HTRDR_SPECTRAL_LW: Tref = 290; break;
- case HTRDR_SPECTRAL_SW: Tref = 5778; /* Tsun */ break;
- default: FATAL("Unreachable code.\n"); break;
+ if(htrdr->ref_temperature <= 0) {
+ htrdr_log_err(htrdr, "%s: invalid reference temperature %g K.\n",
+ FUNC_NAME, htrdr->ref_temperature);
+ res = RES_BAD_ARG;
+ goto error;
}
ASSERT(htrdr->wlen_range_m[0] <= htrdr->wlen_range_m[1]);
n = (size_t)(spectral_range[1] - spectral_range[0]);
res = htrdr_ran_wlen_create
- (htrdr, spectral_range, n, Tref, &htrdr->ran_wlen);
+ (htrdr, spectral_range, n, htrdr->ref_temperature, &htrdr->ran_wlen);
if(res != RES_OK) goto error;
}
@@ -766,66 +767,6 @@ compute_sky_min_band_len
}
res_T
-brightness_temperature
- (struct htrdr* htrdr,
- const double lambda_min,
- const double lambda_max,
- const double radiance, /* In W/m2/sr/m */
- double* temperature)
-{
- const size_t MAX_ITER = 100;
- const double epsilon_T = 1e-4; /* In K */
- const double epsilon_B = radiance * 1e-8;
- double T, T0, T1, T2;
- double B, B0;
- size_t i;
- res_T res = RES_OK;
- ASSERT(temperature && lambda_min <= lambda_max);
-
- /* Search for a brightness temperature whose radiance is greater than or
- * equal to the estimated radiance */
- T2 = 200;
- FOR_EACH(i, 0, MAX_ITER) {
- const double B2 = planck(lambda_min, lambda_max, T2);
- if(B2 >= radiance) break;
- T2 *= 2;
- }
- if(i >= MAX_ITER) { res = RES_BAD_OP; goto error; }
-
- B0 = T0 = T1 = 0;
- FOR_EACH(i, 0, MAX_ITER) {
- T = (T1+T2)*0.5;
- B = planck(lambda_min, lambda_max, T);
-
- if(B < radiance) {
- T1 = T;
- } else {
- T2 = T;
- }
-
- if(fabs(T-T0) < epsilon_T || fabs(B-B0) < epsilon_B)
- break;
-
- T0 = T;
- B0 = B;
- }
- if(i >= MAX_ITER) { res = RES_BAD_OP; goto error; }
-
- *temperature = T;
-
-exit:
- return res;
-error:
- htrdr_log_err(htrdr,
- "Could not compute the brightness temperature for the estimated radiance %g "
- "averaged over [%g, %g] nanometers.\n",
- radiance,
- lambda_min*1e9,
- lambda_max*1e9);
- goto exit;
-}
-
-res_T
open_output_stream
(struct htrdr* htrdr,
const char* filename,
diff --git a/src/htrdr.h b/src/htrdr.h
@@ -17,6 +17,8 @@
#ifndef HTRDR_H
#define HTRDR_H
+#include "htrdr_spectral.h"
+
#include <rsys/logger.h>
#include <rsys/ref_count.h>
#include <rsys/str.h>
@@ -30,12 +32,6 @@
#define HTRDR(Func) htrdr_ ## Func
#endif
-enum htrdr_spectral_type {
- HTRDR_SPECTRAL_LW, /* Longwave */
- HTRDR_SPECTRAL_SW, /* Shortwave */
- HTRDR_SPECTRAL_SW_CIE_XYZ /* Shortwave wrt the CIE XYZ tristimulus */
-};
-
/* Forward declarations */
struct htsky;
struct htrdr_args;
@@ -64,6 +60,7 @@ struct htrdr {
struct htsky* sky;
enum htrdr_spectral_type spectral_type;
double wlen_range_m[2]; /* Integration range in *meters* */
+ double ref_temperature; /* Reference temperature in Kelvin */
size_t spp; /* #samples per pixel */
size_t width; /* Image width */
diff --git a/src/htrdr_args.c b/src/htrdr_args.c
@@ -323,6 +323,13 @@ parse_spectral_parameter(struct htrdr_args* args, const char* str)
args->spectral_type = HTRDR_SPECTRAL_LW;
res = parse_spectral_range(val, args->wlen_range);
if(res != RES_OK) goto error;
+ } else if(!strcmp(key, "Tref")) {
+ res = cstr_to_double(val, &args->ref_temperature);
+ if(res == RES_OK && args->ref_temperature < 0) res = RES_BAD_ARG;
+ if(res != RES_OK) {
+ fprintf(stderr, "Invalid reference temperature Tref=%s.\n", val);
+ goto error;
+ }
} else {
fprintf(stderr, "Invalid spectral parameter `%s'.\n", key);
res = RES_BAD_ARG;
@@ -534,6 +541,22 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
goto error;
}
+ /* Setup default ref temperature if necessary */
+ if(args->ref_temperature <= 0) {
+ switch(args->spectral_type) {
+ case HTRDR_SPECTRAL_LW:
+ args->ref_temperature = HTRDR_DEFAULT_LW_REF_TEMPERATURE;
+ break;
+ case HTRDR_SPECTRAL_SW:
+ args->ref_temperature = HTRDR_SUN_TEMPERATURE;
+ break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ args->ref_temperature = -1; /* Unused */
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ }
+
exit:
return res;
error:
diff --git a/src/htrdr_args.h.in b/src/htrdr_args.h.in
@@ -16,7 +16,7 @@
#ifndef HTRDR_ARGS_H
#define HTRDR_ARGS_H
-#include "htrdr.h"
+#include "htrdr_spectral.h"
#include "htrdr_cie_xyz.h"
#include <float.h>
@@ -58,6 +58,7 @@ struct htrdr_args {
enum htrdr_spectral_type spectral_type;
double wlen_range[2]; /* Spectral range of integration in nm */
+ double ref_temperature; /* Planck reference temperature in Kelvin */
unsigned nthreads; /* Hint on the number of threads to use */
int force_overwriting;
@@ -94,8 +95,9 @@ struct htrdr_args {
90, /* Sun elevation */ \
@HTRDR_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD@, /* Optical thickness */ \
{UINT_MAX, UINT_MAX, UINT_MAX}, /* Maximum definition of the grid */ \
- HTRDR_SPECTRAL_SW_CIE_XYZ, \
+ HTRDR_SPECTRAL_SW_CIE_XYZ, /* Spectral type */ \
HTRDR_CIE_XYZ_RANGE_DEFAULT__, /* Spectral range */ \
+ -1, /* Reference temperature */ \
(unsigned)~0, /* #threads */ \
0, /* Force overwriting */ \
0, /* dump VTK */ \
diff --git a/src/htrdr_c.h b/src/htrdr_c.h
@@ -18,7 +18,6 @@
#define HTRDR_C_H
#include <rsys/rsys.h>
-#include <rsys/math.h> /* PI support */
#ifndef NDEBUG
#define MPI(Func) ASSERT(MPI_##Func == MPI_SUCCESS)
@@ -36,9 +35,6 @@ enum htrdr_mpi_message {
struct htrdr;
-#define SW_WAVELENGTH_MIN 380 /* In nanometer */
-#define SW_WAVELENGTH_MAX 780 /* In nanometer */
-
/* In nanometer */
static FINLINE double
wavenumber_to_wavelength(const double nu/*In cm^-1*/)
@@ -96,102 +92,6 @@ morton_xyz_decode_u21(const uint64_t code, uint32_t xyz[3])
xyz[2] = (uint32_t)morton3D_decode_u21(code >> 0);
}
-static INLINE double
-wiebelt(const double v)
-{
- int m;
- double w, v2, v4;
- /*.153989717364e+00;*/
- const double fifteen_over_pi_power_4 = 15.0/(PI*PI*PI*PI);
- const double z0 = 1.0/3.0;
- const double z1 = 1.0/8.0;
- const double z2 = 1.0/60.0;
- const double z4 = 1.0/5040.0;
- const double z6 = 1.0/272160.0;
- const double z8 = 1.0/13305600.0;
-
- if(v >= 2.) {
- w = 0;
- for(m=1; m<6 ;m++)
- w+=exp(-m*v)/(m*m*m*m) * (((m*v+3)*m*v+6)*m*v+6);
- w = w * fifteen_over_pi_power_4;
- } else {
- v2 = v*v;
- v4 = v2*v2;
- w = z0 - z1*v + z2*v2 - z4*v2*v2 + z6*v4*v2 - z8*v4*v4;
- w = 1. - fifteen_over_pi_power_4*v2*v*w;
- }
- ASSERT(w >= 0.0 && w <= 1.0);
- return w;
-}
-
-static INLINE double
-blackbody_fraction
- (const double lambda0, /* In meters */
- const double lambda1, /* In meters */
- const double temperature) /* In Kelvin */
-{
- const double C2 = 1.43877735e-2; /* m.K */
- double x0 = C2 / lambda0;
- double x1 = C2 / lambda1;
- double v0 = x0 / temperature;
- double v1 = x1 / temperature;
- double w0 = wiebelt(v0);
- double w1 = wiebelt(v1);
- return w1 - w0;
-}
-
-
-
-/* Return the Planck value in W/m^2/sr/m at a given wavelength */
-static INLINE double
-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;
- const double T4 = T2*T2;
- 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 * (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);
- }
-}
-
/* Return the minimum length in nanometer of the sky spectral bands
* clamped to in [range[0], range[1]]. */
extern LOCAL_SYM double
@@ -199,15 +99,6 @@ compute_sky_min_band_len
(struct htsky* sky,
const double range[2]);
-extern LOCAL_SYM res_T
-brightness_temperature
- (struct htrdr* htrdr,
- const double lambda_min, /* In meters */
- const double lambda_max, /* In meters */
- /* Averaged over [lambda_min, lambda_max]. In W/m^2/sr/m */
- const double radiance,
- double* temperature);
-
extern LOCAL_SYM res_T
open_output_stream
(struct htrdr* htrdr,
diff --git a/src/htrdr_ran_wlen.h b/src/htrdr_ran_wlen.h
@@ -27,8 +27,6 @@ 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 */
diff --git a/src/htrdr_spectral.c b/src/htrdr_spectral.c
@@ -0,0 +1,79 @@
+/* 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/>. */
+
+#include "htrdr.h"
+#include "htrdr_spectral.h"
+
+res_T
+brightness_temperature
+ (struct htrdr* htrdr,
+ const double lambda_min,
+ const double lambda_max,
+ const double radiance, /* In W/m2/sr/m */
+ double* temperature)
+{
+ const size_t MAX_ITER = 100;
+ const double epsilon_T = 1e-4; /* In K */
+ const double epsilon_B = radiance * 1e-8;
+ double T, T0, T1, T2;
+ double B, B0;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(temperature && lambda_min <= lambda_max);
+
+ /* Search for a brightness temperature whose radiance is greater than or
+ * equal to the estimated radiance */
+ T2 = 200;
+ FOR_EACH(i, 0, MAX_ITER) {
+ const double B2 = planck(lambda_min, lambda_max, T2);
+ if(B2 >= radiance) break;
+ T2 *= 2;
+ }
+ if(i >= MAX_ITER) { res = RES_BAD_OP; goto error; }
+
+ B0 = T0 = T1 = 0;
+ FOR_EACH(i, 0, MAX_ITER) {
+ T = (T1+T2)*0.5;
+ B = planck(lambda_min, lambda_max, T);
+
+ if(B < radiance) {
+ T1 = T;
+ } else {
+ T2 = T;
+ }
+
+ if(fabs(T-T0) < epsilon_T || fabs(B-B0) < epsilon_B)
+ break;
+
+ T0 = T;
+ B0 = B;
+ }
+ if(i >= MAX_ITER) { res = RES_BAD_OP; goto error; }
+
+ *temperature = T;
+
+exit:
+ return res;
+error:
+ htrdr_log_err(htrdr,
+ "Could not compute the brightness temperature for the estimated radiance %g "
+ "averaged over [%g, %g] nanometers.\n",
+ radiance,
+ lambda_min*1e9,
+ lambda_max*1e9);
+ goto exit;
+}
+
diff --git a/src/htrdr_spectral.h b/src/htrdr_spectral.h
@@ -0,0 +1,137 @@
+/* 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_SPECTRAL_H
+#define HTRDR_SPECTRAL_H
+
+#include <rsys/rsys.h>
+#include <rsys/math.h> /* PI support */
+
+#define HTRDR_SUN_TEMPERATURE 5778 /* In K */
+#define HTRDR_DEFAULT_LW_REF_TEMPERATURE 290 /* Default longwave temperature in K */
+
+enum htrdr_spectral_type {
+ HTRDR_SPECTRAL_LW, /* Longwave */
+ HTRDR_SPECTRAL_SW, /* Shortwave */
+ HTRDR_SPECTRAL_SW_CIE_XYZ /* Shortwave wrt the CIE XYZ tristimulus */
+};
+
+struct htrdr;
+
+static INLINE double
+wiebelt(const double v)
+{
+ int m;
+ double w, v2, v4;
+ /*.153989717364e+00;*/
+ const double fifteen_over_pi_power_4 = 15.0/(PI*PI*PI*PI);
+ const double z0 = 1.0/3.0;
+ const double z1 = 1.0/8.0;
+ const double z2 = 1.0/60.0;
+ const double z4 = 1.0/5040.0;
+ const double z6 = 1.0/272160.0;
+ const double z8 = 1.0/13305600.0;
+
+ if(v >= 2.) {
+ w = 0;
+ for(m=1; m<6 ;m++)
+ w+=exp(-m*v)/(m*m*m*m) * (((m*v+3)*m*v+6)*m*v+6);
+ w = w * fifteen_over_pi_power_4;
+ } else {
+ v2 = v*v;
+ v4 = v2*v2;
+ w = z0 - z1*v + z2*v2 - z4*v2*v2 + z6*v4*v2 - z8*v4*v4;
+ w = 1. - fifteen_over_pi_power_4*v2*v*w;
+ }
+ ASSERT(w >= 0.0 && w <= 1.0);
+ return w;
+}
+
+static INLINE double
+blackbody_fraction
+ (const double lambda0, /* In meters */
+ const double lambda1, /* In meters */
+ const double temperature) /* In Kelvin */
+{
+ const double C2 = 1.43877735e-2; /* m.K */
+ double x0 = C2 / lambda0;
+ double x1 = C2 / lambda1;
+ double v0 = x0 / temperature;
+ double v1 = x1 / temperature;
+ double w0 = wiebelt(v0);
+ double w1 = wiebelt(v1);
+ return w1 - w0;
+}
+
+/* Return the Planck value in W/m^2/sr/m at a given wavelength */
+static INLINE double
+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;
+ const double T4 = T2*T2;
+ 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 * (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 meters */
+ const double lambda_max, /* In meters */
+ /* Averaged over [lambda_min, lambda_max]. In W/m^2/sr/m */
+ const double radiance,
+ double* temperature);
+
+#endif /* HTRDR_SPECTRAL_H */
