commit 7a48312dc934278dab83684b1c237e0a0232ddac
parent d5a0b5dd94ab262f2af18291f11b0c60b9eefc6d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 4 Jun 2020 17:49:28 +0200
Merge branch 'feature_shortwave' into develop
Diffstat:
21 files changed, 1075 insertions(+), 800 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -25,7 +25,7 @@ option(NO_TEST "Do not build the tests" OFF)
# Check dependencies
################################################################################
find_package(AW 2.0 REQUIRED)
-find_package(HTSky 0.1 REQUIRED)
+find_package(HTSky 0.2 REQUIRED)
find_package(MruMtl 0.0 REQUIRED)
find_package(RCMake 0.3 REQUIRED)
find_package(RSys 0.7 REQUIRED)
@@ -100,9 +100,10 @@ set(HTRDR_FILES_SRC
htrdr_interface.c
htrdr_main.c
htrdr_mtl.c
- htrdr_ran_lw.c
+ htrdr_ran_wlen.c
htrdr_rectangle.c
htrdr_slab.c
+ htrdr_spectral.c
htrdr_sun.c)
set(HTRDR_FILES_INC
htrdr.h
@@ -114,9 +115,10 @@ set(HTRDR_FILES_INC
htrdr_ground.h
htrdr_interface.h
htrdr_mtl.h
- htrdr_ran_lw.h
+ 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/doc/htrdr-image.5.txt b/doc/htrdr-image.5.txt
@@ -30,25 +30,33 @@ thus ignored as well as empty lines. The first valid line stores 2 unsigned
integers that represent the image definition, i.e. the number of pixels per
line and per column. Then each line stores 8 pixel components.
-If the image is a regular rendering in the visible part of the spectrum, the
-pixel components are actually 4 pairs of floating points data representing the
-pixel color encoded in the CIE 1931 XYZ color space and the per radiative path
-computation time. The first, second and third pair encode the estimated
-integrated radiance in W.sr^-1.m^-2 of the X, Y and Z pixel components,
-respectively. The first value of each pair is the expected value of the
-estimated radiance while the second one is its associated standard deviation.
-The fourth pair saves the estimate in microseconds of the per radiative path
-computation time and its standard error.
-
-If the image is an infrared rendering, the first and second pixel components
-store the expected value and the standard error, respectively, of the
-estimated brightness temperature in Kelvin. The third and fourth component
-save the expected value and standard deviation of the pixel radiance in
-W.sr^-1.m^-2.nm^-1. The fifth and sixth pixel components are unused. Finally
-the last 2 pixel components save, as for a regular rendering, the estimate in
+If the image is a regular rendering in the visible part of the spectrum
+(*-s* _cie_xyz_ option in *htrdr*(1)), the pixel components are actually 4
+pairs of floating points data representing the pixel color encoded in the CIE
+1931 XYZ color space and the per radiative path computation time. The first,
+second and third pairs encode the estimated integrated radiance in W/sr/m^2 of
+the X, Y and Z pixel components, respectively. The first value of each pair is
+the expected value of the estimated radiance while the second one is its
+associated standard deviation. The fourth pair saves the estimate in
microseconds of the per radiative path computation time and its standard
error.
+If the image is an infrared rendering (*-s* *lw*=_wlen-min_,_wlen_max_ option
+in *htrdr*(1)), the first and second pixel components store the expected value
+and the standard error, respectively, of the estimated brightness temperature
+in Kelvin. The third and fourth components save the expected value and standard
+deviation of the pixel radiance that is either an integrated radiance in
+W/sr/m^2 or a spectral radiance in W/sr/m^2/nm whether this radiance is
+computed for a spectral range or for one wavelength. The fifth and sixth pixel
+components are unused. Finally the last 2 pixel components save, as for a
+regular rendering, the estimate in microseconds of the per radiative path
+computation time and its standard error.
+
+If it was generating from a shortwave rendering (*-s*
+*sw*=_wlen-min_,wlen-max_ option in *htrdr*(1)) the image is formatted as in
+longwave rendering mode exepted that the first and second pixel components are
+unused since no brightness temperature was evaluated in shortwave.
+
Pixels are sorted line by line, with the origin defined at the top left corner
of the image. With an image definition of N by M pixels, with N the number of
pixels per line and M the overall number of lines in the image, the first N
diff --git a/doc/htrdr.1.txt.in b/doc/htrdr.1.txt.in
@@ -30,7 +30,13 @@ SYNOPSIS
DESCRIPTION
-----------
*htrdr* is an image renderer of scenes composed of an atmospheric gas mixture,
-clouds, and a ground. Images can be rendered in the visible or the infrared
+clouds, and a ground. It evaluates the intensity incoming on each pixel of the
+sensor array. The underlying algorithm is based on a Monte-Carlo method: it
+consists in simulating a given number of optical paths originating from the
+camera, directed into the atmosphere, taking into account light absorption and
+scattering phenomena.
+
+Images can be rendered in the visible or the infrared
part of the spectrum. It uses spectral data that should be provided for the
pressure and temperature atmospheric vertical profile [1] (*-a* _atmosphere_),
the liquid water content in suspension within the clouds stored in a *htcp*(5)
@@ -44,27 +50,28 @@ whose materials are listed in the *htrdr-material*(5) file provided through the
*-M* option. Both, the clouds and the ground, can be infinitely repeated along
the X and Y axis by setting the *-r* and the *-R* options, respectively.
-*htrdr* evaluates the intensity incoming on each pixel of the sensor array. The
-underlying algorithm is based on a Monte-Carlo method: it consists in
-simulating a given number of optical paths originating from the camera,
-directed into the atmosphere, taking into account light absorption and
-scattering phenomena. When rendering in the visible part of the spectrum, the
-computation is performed over the visible part of the spectrum in [380, 780]
+Spectral dimension can be integrated in many ways (*-s* option). By default,
+the computation is performed for the visible part of the spectrum in [380, 780]
nanometers, for the three components of the CIE 1931 XYZ colorimetric space
that are subsequently recombined in order to obtain the final color for each
pixel, and finally the whole image of the scene as seen from the set
-observation position. In longwave, the rendering is performed for the range
-of wavelengths defined by the *-l* option. The estimated radiance per pixel is
-then converted to its brightness temperature and both are saved in the output
-image.
-
-In *htrdr* the spatial unit 1.0 corresponds to one meter, the estimated
-radiance is given in W.sr^-1.m^-2 and the temperatures are expressed in Kelvin.
-The results are written to the output file if the *-o* option is defined and
-the standard output otherwise. The output image is a list of raw ASCII data
-formatted with respect to the *htrdr-image*(5) file format. Since *htrdr*
-relies on the Monte-Carlo method, each estimation is given with its numerical
-accuracy.
+observation position. The two other ways consist in explicitly defining the
+longwave or shortwave spectral range to handle and continuously sampling a
+wavelength in this range according to the Planck function for a reference
+temperature. Actually longwave and shortwave are keywords that mean that the
+source of radiation is whether external or internal to the medium,
+respectively. In shortwave, only the pixel radiance is evaluated and stored in
+the output image. In longwave this estimated radiance is then converted to its
+brightness temperature and both are saved in the image.
+
+In *htrdr* the spatial unit 1.0 corresponds to one meter and the temperatures
+are expressed in Kelvin. The estimated radiances are given in W/sr/m^2 excepted
+for monochromatic computations where the computed spectral radiance is defined
+in W/sr/m^2/nm. The results are written to the output file if the *-o* option
+is defined and the standard output otherwise. The output image is a list of raw
+ASCII data formatted with respect to the *htrdr-image*(5) file format. Since
+*htrdr* relies on the Monte-Carlo method, each estimation is given with its
+numerical accuracy.
During the simulation, *htrdr* dynamically loads/unloads cloud properties to
handle clouds whose data that do not feat in main memory. *htrdr* also supports
@@ -146,10 +153,6 @@ OPTIONS
brightness temperature for the submitted range of wavelengths. By default,
*spp* is set to @HTRDR_ARGS_DEFAULT_IMG_SPP@.
-*-l* __wlen-min__,__wlen-max__::
- Switch in infrared rendering for the longwave interval in [__wlen-min__,
- __wlen-max__] nanometers.
-
*-R*::
Infinitely repeat the _ground_ along the X and Y axis.
@@ -180,6 +183,44 @@ OPTIONS
File where *htrdr* writes its _output_ data. If not defined, write results to
standard output.
+*-s* <__spectral-parameter__:...>::
+ define the type and the range of the spectral integration. Available spectral
+ parameters are:
+
+ **cie_xyz**;;
+ the radiance is computed for the visible part of the spectrum in [380, 780]
+ nanometers with respect to the XYZ CIE 1931 tristimulus values. This is the
+ default comportment of *htrdr*.
+
+ **lw**=*_wlen-min_*,*_wlen-max_*;;
+ perform the spectral sampling continuously in the [_wlen-min_, _wlen-max_]
+ wavelength range (wavelength must be provided in nanometers) according to
+ the Planck function for a reference temperature. If _wlen-min_ and
+ _wlen-max_ are equals the computation is monochromatic. *lw* means for
+ longwave but is here a code word that really means "computation of radiance
+ using the internal source of radiation": in other words, radiation is
+ emitted by the medium and its boundaries. Because the application is for
+ the terrestrial atmosphere, internal radiation is emitted in the thermal
+ longwave part of the electromagnetic spectrum. Therefore the default value
+ of the reference temperature used in the spectral sampling is fixed at 290
+ K.
+
+ **sw**=*_wlen-min_*,*_wlen-max_*;;
+ perform the spectral sampling continuously in the [_wlen-min_, _wlen-max_]
+ wavelength range (wavelength must be provided in nanometers) according to
+ the Planck function for a reference temperature. If _wlen-min_ and
+ _wlen-max_ are equals the computation is monochromatic. In the present
+ case, *sw* means that the source of radiation is external to the medium:
+ because the application is for the terrestrial atmosphere, the value of the
+ reference temperature is by default fixed at 5778 K, i.e. the blackbody
+ temperature of the Sun.
+
+ **Tref**=**_temperature_**;;
+ reference temperature of the Planck function used to continuously sample the
+ longwave/shortwave spectral range. In longwave, it is set to 290 K by
+ default while in shortwave its default value is the blackbody temperature
+ of the sun (i.e. 5778 K).
+
*-T* _threshold_::
Optical thickness used as threshold criteria to partition the properties of
the clouds. By default its value is @HTRDR_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD@.
@@ -237,13 +278,14 @@ PPM image [4]:
$ htpp -o image.ppm output
Render the previous scene in infrared for the wavelengths in [9200, 10000]
-nanometers:
+nanometers with a reference temperature of 300 Kelvin:
- $ htrdr -l 9200,10000 -a gas.txt -m Mie.nc -g mountains.obj -R \
+ $ htrdr -a gas.txt -m Mie.nc -g mountains.obj -R \
-M materials.mtl \
-c clouds.htcp \
-C pos=0,0,400:tgt=0,1,0:up=0,0,1 \
-i def=800x600:spp=64 \
+ -s lw=9200,1000:Tref=300
-f -o output
Move the sun by setting its azimuthal and elevation angles to *120* and *40*
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_ran_lw.h"
+#include "htrdr_ran_wlen.h"
#include "htrdr_sun.h"
#include "htrdr_solve.h"
@@ -102,6 +102,23 @@ log_msg
}
}
+static enum htsky_spectral_type
+htrdr_to_sky_spectral_type(const enum htrdr_spectral_type type)
+{
+ enum htsky_spectral_type spectype;
+ switch(type) {
+ case HTRDR_SPECTRAL_LW:
+ spectype = HTSKY_SPECTRAL_LW;
+ break;
+ case HTRDR_SPECTRAL_SW:
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ spectype = HTSKY_SPECTRAL_SW;
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return spectype;
+}
+
static res_T
dump_buffer
(struct htrdr* htrdr,
@@ -117,13 +134,8 @@ dump_buffer
ASSERT(htrdr && buf && stream_name && stream);
(void)stream_name;
- if(htsky_is_long_wave(htrdr->sky)) {
- pixsz = sizeof(struct htrdr_pixel_lw);
- pixal = ALIGNOF(struct htrdr_pixel_lw);
- } else {
- pixsz = sizeof(struct htrdr_pixel_sw);
- pixal = ALIGNOF(struct htrdr_pixel_sw);
- }
+ pixsz = htrdr_spectral_type_get_pixsz(htrdr->spectral_type);
+ pixal = htrdr_spectral_type_get_pixal(htrdr->spectral_type);
htrdr_buffer_get_layout(buf, &layout);
if(layout.elmt_size != pixsz || layout.alignment != pixal) {
@@ -140,8 +152,8 @@ dump_buffer
struct htrdr_estimate pix_time = HTRDR_ESTIMATE_NULL;
const struct htrdr_accum* pix_time_acc = NULL;
- if(htsky_is_long_wave(htrdr->sky)) {
- const struct htrdr_pixel_lw* pix = htrdr_buffer_at(buf, x, y);
+ if(htrdr->spectral_type != HTRDR_SPECTRAL_SW_CIE_XYZ){
+ const struct htrdr_pixel_xwave* pix = htrdr_buffer_at(buf, x, y);
fprintf(stream, "%g %g ",
pix->radiance_temperature.E, pix->radiance_temperature.SE);
fprintf(stream, "%g %g ", pix->radiance.E, pix->radiance.SE);
@@ -149,7 +161,7 @@ dump_buffer
pix_time_acc = &pix->time;
} else {
- const struct htrdr_pixel_sw* pix = htrdr_buffer_at(buf, x, y);
+ const struct htrdr_pixel_image* pix = htrdr_buffer_at(buf, x, y);
fprintf(stream, "%g %g ", pix->X.E, pix->X.SE);
fprintf(stream, "%g %g ", pix->Y.E, pix->Y.SE);
fprintf(stream, "%g %g ", pix->Z.E, pix->Z.SE);
@@ -389,6 +401,7 @@ htrdr_init
struct htsky_args htsky_args = HTSKY_ARGS_DEFAULT;
double proj_ratio;
double sun_dir[3];
+ double spectral_range[2];
const char* output_name = NULL;
size_t ithread;
int nthreads_max;
@@ -414,6 +427,8 @@ htrdr_init
htrdr->grid_max_definition[0] = args->grid_max_definition[0];
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;
@@ -484,32 +499,48 @@ htrdr_init
htsky_args.nthreads = htrdr->nthreads;
htsky_args.repeat_clouds = args->repeat_clouds;
htsky_args.verbose = htrdr->mpi_rank == 0 ? args->verbose : 0;
- htsky_args.wlen_lw_range[0] = args->wlen_lw_range[0];
- htsky_args.wlen_lw_range[1] = args->wlen_lw_range[1];
+ htsky_args.spectral_type = htrdr_to_sky_spectral_type(args->spectral_type);
+ htsky_args.wlen_range[0] = args->wlen_range[0];
+ htsky_args.wlen_range[1] = args->wlen_range[1];
res = htsky_create(&htrdr->logger, htrdr->allocator, &htsky_args, &htrdr->sky);
if(res != RES_OK) goto error;
- if(!htsky_is_long_wave(htrdr->sky)) { /* Short wave random variate */
- const double* range = HTRDR_CIE_XYZ_RANGE_DEFAULT;
- size_t n;
+ HTSKY(get_raw_spectral_bounds(htrdr->sky, spectral_range));
+
+ spectral_range[0] = MMAX(args->wlen_range[0], spectral_range[0]);
+ spectral_range[1] = MMIN(args->wlen_range[1], spectral_range[1]);
+ if(spectral_range[0] != args->wlen_range[0]
+ || spectral_range[1] != args->wlen_range[1]) {
+ htrdr_log_warn(htrdr,
+ "%s: the submitted spectral range overflowed the spectral data.\n", FUNC_NAME);
+ }
+
+ htrdr->wlen_range_m[0] = spectral_range[0]*1e-9; /* Convert in meters */
+ htrdr->wlen_range_m[1] = spectral_range[1]*1e-9; /* Convert in meters */
- n = (size_t)(range[1] - range[0]);
- res = htrdr_cie_xyz_create(htrdr, range, n, &htrdr->cie);
+ if(htrdr->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ) {
+ size_t n;
+ n = (size_t)(spectral_range[1] - spectral_range[0]);
+ res = htrdr_cie_xyz_create(htrdr, spectral_range, n, &htrdr->cie);
if(res != RES_OK) goto error;
- } else { /* Long Wave random variate */
- const double Tref = 290; /* In Kelvin */
+ } else {
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 */
+ 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]);
- 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);
+ res = htrdr_ran_wlen_create
+ (htrdr, spectral_range, n, htrdr->ref_temperature, &htrdr->ran_wlen);
if(res != RES_OK) goto error;
- }
+ }
htrdr->lifo_allocators = MEM_CALLOC
(htrdr->allocator, htrdr->nthreads, sizeof(*htrdr->lifo_allocators));
@@ -535,16 +566,9 @@ htrdr_init
/* Create the image buffer only on the master process; the image parts
* rendered by the processes are gathered onto the master process. */
if(!htrdr->dump_vtk && htrdr->mpi_rank == 0) {
- size_t pixsz = 0; /* sizeof(pixel) */
- size_t pixal = 0; /* alignof(pixel) */
+ const size_t pixsz = htrdr_spectral_type_get_pixsz(htrdr->spectral_type);
+ const size_t pixal = htrdr_spectral_type_get_pixal(htrdr->spectral_type);
- if(htsky_is_long_wave(htrdr->sky)) {
- pixsz = sizeof(struct htrdr_pixel_lw);
- pixal = ALIGNOF(struct htrdr_pixel_lw);
- } else {
- pixsz = sizeof(struct htrdr_pixel_sw);
- pixal = ALIGNOF(struct htrdr_pixel_sw);
- }
res = htrdr_buffer_create(htrdr,
args->image.definition[0], /* Width */
args->image.definition[1], /* Height */
@@ -575,7 +599,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->ran_lw) htrdr_ran_lw_ref_put(htrdr->ran_lw);
+ 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;
@@ -743,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>
@@ -50,13 +52,15 @@ struct htrdr {
struct htrdr_mtl* mtl;
struct htrdr_sun* sun;
struct htrdr_cie_xyz* cie;
- struct htrdr_ran_lw* ran_lw;
+ struct htrdr_ran_wlen* ran_wlen;
struct htrdr_camera* cam;
struct htrdr_buffer* buf;
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
@@ -61,16 +61,6 @@ print_help(const char* cmd)
" -i <image> define the image to compute. Refer to the htrdr man\n"
" page for the list of image options\n");
printf(
-" -l WLEN_MIN,WLEN_MAX\n"
-" switch in infrared rendering for longwave in\n"
-" [WLEN_MIN, WLEN_MAX], in nanometers. By default, the\n"
-" rendering is performed for the visible part of the\n"
-" spectrum in [380, 780] nanometers.\n");
- printf(
-" -R infinitely repeat the ground along the X and Y axis.\n");
- printf(
-" -r infinitely repeat the clouds along the X and Y axis.\n");
- printf(
" -M MATERIALS file listing the ground materials.\n");
printf(
" -m MIE file of Mie's data.\n");
@@ -81,6 +71,14 @@ print_help(const char* cmd)
" -o OUTPUT file where data are written. If not defined, data are\n"
" written to standard output.\n");
printf(
+" -R infinitely repeat the ground along the X and Y axis.\n");
+ printf(
+" -r infinitely repeat the clouds along the X and Y axis.\n");
+ printf(
+" -s <spectral> define the type and range of the spectral\n"
+" integration. Refer to the htrdr man page for the list\n"
+" of spectral options\n");
+ printf(
" -T THRESHOLD optical thickness used as threshold during the octree\n"
" building. By default its value is `%g'.\n",
HTRDR_ARGS_DEFAULT.optical_thickness);
@@ -140,7 +138,7 @@ parse_fov(const char* str, double* out_fov)
fprintf(stderr, "Invalid field of view `%s'.\n", str);
return RES_BAD_ARG;
}
- if(fov < 30 || fov > 120) {
+ if(fov <= 0 || fov >= 180) {
fprintf(stderr, "The field of view %g is not in [30, 120].\n", fov);
return RES_BAD_ARG;
}
@@ -217,10 +215,11 @@ parse_camera_parameter(struct htrdr_args* args, const char* str)
char* val;
char* ctx;
res_T res = RES_OK;
+ ASSERT(args);
if(strlen(str) >= sizeof(buf) -1/*NULL char*/) {
fprintf(stderr,
- "Could not duplicate the rectangle option string `%s'.\n", str);
+ "Could not duplicate the camera option string `%s'.\n", str);
res = RES_MEM_ERR;
goto error;
}
@@ -262,6 +261,89 @@ error:
}
static res_T
+parse_spectral_range(const char* str, double wlen_range[2])
+{
+ double range[2];
+ size_t len;
+ res_T res = RES_OK;
+ ASSERT(wlen_range && str);
+
+ res = cstr_to_list_double(str, ',', range, &len, 2);
+ if(res == RES_OK && len != 2) res = RES_BAD_ARG;
+ if(res == RES_OK && range[0] > range[1]) res = RES_BAD_ARG;
+ if(res != RES_OK) {
+ fprintf(stderr, "Invalid spectral range `%s'.\n", str);
+ goto error;
+ }
+ wlen_range[0] = range[0];
+ wlen_range[1] = range[1];
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_spectral_parameter(struct htrdr_args* args, const char* str)
+{
+ char buf[128];
+ char* key;
+ char* val;
+ char* ctx;
+ res_T res = RES_OK;
+ ASSERT(args);
+
+ if(strlen(str) >= sizeof(buf) -1/*NULL char*/) {
+ fprintf(stderr,
+ "Could not duplicate the spectral option string `%s'.\n", str);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ strncpy(buf, str, sizeof(buf));
+
+ key = strtok_r(buf, "=", &ctx);
+ val = strtok_r(NULL, "", &ctx);
+
+ if(!strcmp(key, "cie_xyz")) {
+ args->spectral_type = HTRDR_SPECTRAL_SW_CIE_XYZ;
+ args->wlen_range[0] = HTRDR_CIE_XYZ_RANGE_DEFAULT[0];
+ args->wlen_range[1] = HTRDR_CIE_XYZ_RANGE_DEFAULT[1];
+ } else {
+ if(!val) {
+ fprintf(stderr, "Missing value to the spectral option `%s'.\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ if(!strcmp(key, "sw")) {
+ args->spectral_type = HTRDR_SPECTRAL_SW;
+ res = parse_spectral_range(val, args->wlen_range);
+ if(res != RES_OK) goto error;
+ } else if(!strcmp(key, "lw")) {
+ 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;
+ goto error;
+ }
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
parse_multiple_parameters
(struct htrdr_args* args,
const char* str,
@@ -365,31 +447,6 @@ error:
goto exit;
}
-static res_T
-parse_lw_range(struct htrdr_args* args, const char* str)
-{
- double range[2];
- size_t len;
- res_T res = RES_OK;
- ASSERT(args && str);
-
- res = cstr_to_list_double(str, ',', range, &len, 2);
- if(res == RES_OK && len != 2) res = RES_BAD_ARG;
- if(res == RES_OK && range[0] > range[1]) res = RES_BAD_ARG;
- if(res != RES_OK) {
- fprintf(stderr, "Invalid longwave range `%s'.\n", str);
- goto error;
- }
-
- args->wlen_lw_range[0] = range[0];
- args->wlen_lw_range[1] = range[1];
-
-exit:
- return res;
-error:
- goto exit;
-}
-
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -414,7 +471,7 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
}
}
- while((opt = getopt(argc, argv, "a:C:c:D:dfg:hi:l:M:m:O:o:RrT:t:V:v")) != -1) {
+ while((opt = getopt(argc, argv, "a:C:c:D:dfg:hi:M:m:O:o:Rrs:T:t:V:v")) != -1) {
switch(opt) {
case 'a': args->filename_gas = optarg; break;
case 'C':
@@ -435,15 +492,16 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
res = parse_multiple_parameters
(args, optarg, parse_image_parameter);
break;
- case 'l':
- res = parse_lw_range(args, optarg);
- break;
case 'M': args->filename_mtl = optarg; break;
case 'm': args->filename_mie = optarg; break;
case 'O': args->cache = optarg; break;
case 'o': args->output = optarg; break;
case 'r': args->repeat_clouds = 1; break;
case 'R': args->repeat_ground = 1; break;
+ case 's':
+ res = parse_multiple_parameters
+ (args, optarg, parse_spectral_parameter);
+ break;
case 'T':
res = cstr_to_double(optarg, &args->optical_thickness);
if(res == RES_OK && args->optical_thickness < 0) res = RES_BAD_ARG;
@@ -483,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,8 @@
#ifndef HTRDR_ARGS_H
#define HTRDR_ARGS_H
-#include "htrdr_ground.h"
+#include "htrdr_spectral.h"
+#include "htrdr_cie_xyz.h"
#include <float.h>
#include <limits.h>
@@ -55,7 +56,9 @@ struct htrdr_args {
double optical_thickness; /* Threshold used during octree building */
unsigned grid_max_definition[3]; /* Maximum definition of the grid */
- double wlen_lw_range[2]; /* Long Wave range to handle in nm */
+ 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;
@@ -92,7 +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 */ \
- {DBL_MAX, -DBL_MAX}, /* Long wave range. Degenerated <=> short wave */ \
+ 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_camera.c b/src/htrdr_camera.c
@@ -73,7 +73,7 @@ htrdr_camera_create
ref_init(&cam->ref);
cam->htrdr = htrdr;
- if(fov <= 0) {
+ if(fov <= 0 || fov >= PI) {
htrdr_log_err(htrdr, "invalid horizontal camera field of view `%g'\n", fov);
res = RES_BAD_ARG;
goto error;
diff --git a/src/htrdr_cie_xyz.c b/src/htrdr_cie_xyz.c
@@ -325,7 +325,7 @@ htrdr_cie_xyz_create
res = setup_cie_xyz(cie, FUNC_NAME, nbands);
if(res != RES_OK) goto error;
- htrdr_log(htrdr, "Short wave interval defined on [%g, %g] nanometers.\n",
+ htrdr_log(htrdr, "CIE XYZ spectral interval defined on [%g, %g] nanometers.\n",
range[0], range[1]);
exit:
diff --git a/src/htrdr_cie_xyz.h b/src/htrdr_cie_xyz.h
@@ -22,7 +22,9 @@ struct htrdr;
struct htrdr_cie_xyz;
/* Wavelength boundaries of the CIE XYZ color space in nanometers */
-static const double HTRDR_CIE_XYZ_RANGE_DEFAULT[2] = {380, 780};
+#define HTRDR_CIE_XYZ_RANGE_DEFAULT__ {380, 780}
+static const double HTRDR_CIE_XYZ_RANGE_DEFAULT[2] =
+ HTRDR_CIE_XYZ_RANGE_DEFAULT__;
extern LOCAL_SYM res_T
htrdr_cie_xyz_create
diff --git a/src/htrdr_compute_radiance_lw.c b/src/htrdr_compute_radiance_lw.c
@@ -1,4 +1,4 @@
-/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
+/*Spectralt (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
@@ -19,7 +19,6 @@
#include "htrdr_interface.h"
#include "htrdr_ground.h"
#include "htrdr_solve.h"
-#include "htrdr_ran_lw.h"
#include <high_tune/htsky.h>
diff --git a/src/htrdr_compute_radiance_sw.c b/src/htrdr_compute_radiance_sw.c
@@ -279,7 +279,8 @@ htrdr_compute_radiance_sw
double g = 0; /* Asymmetry parameter of the HG phase function */
ASSERT(htrdr && rng && pos_in && dir_in && ithread < htrdr->nthreads);
- ASSERT(!htsky_is_long_wave(htrdr->sky));
+ ASSERT(htrdr->spectral_type == HTRDR_SPECTRAL_SW
+ || htrdr->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ);
CHK(RES_OK == ssf_phase_create
(&htrdr->lifo_allocators[ithread], &ssf_phase_hg, &phase_hg));
@@ -298,7 +299,6 @@ htrdr_compute_radiance_sw
ASSERT(band_bounds[0] <= wlen && wlen <= band_bounds[1]);
sun_solid_angle = htrdr_sun_get_solid_angle(htrdr->sun);
L_sun = htrdr_sun_get_radiance(htrdr->sun, wlen);
-
d3_set(pos, pos_in);
d3_set(dir, dir_in);
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_ran_lw.h"
+#include "htrdr_ran_wlen.h"
#include "htrdr_solve.h"
#include <high_tune/htsky.h>
@@ -46,13 +46,13 @@
STATIC_ASSERT(IS_POW2(TILE_SIZE), TILE_SIZE_must_be_a_power_of_2);
enum pixel_format {
- PIXEL_LW,
- PIXEL_SW
+ PIXEL_XWAVE,
+ PIXEL_IMAGE
};
union pixel {
- struct htrdr_pixel_lw lw;
- struct htrdr_pixel_sw sw;
+ struct htrdr_pixel_xwave xwave;
+ struct htrdr_pixel_image image;
};
/* Tile of row ordered image pixels */
@@ -101,6 +101,19 @@ morton2D_encode(const uint16_t u16)
return u32;
}
+static INLINE enum pixel_format
+spectral_type_to_pixfmt(const enum htrdr_spectral_type spectral_type)
+{
+ enum pixel_format pixfmt;
+ switch(spectral_type) {
+ case HTRDR_SPECTRAL_LW: pixfmt = PIXEL_XWAVE; break;
+ case HTRDR_SPECTRAL_SW: pixfmt = PIXEL_XWAVE; break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ: pixfmt = PIXEL_IMAGE; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return pixfmt;
+}
+
static FINLINE struct tile*
tile_create(struct mem_allocator* allocator, const enum pixel_format fmt)
{
@@ -155,9 +168,11 @@ tile_at
ASSERT(tile && x < TILE_SIZE && y < TILE_SIZE);
return tile->data.pixels + (y*TILE_SIZE + x);
}
+
static void
write_tile_data
- (struct htrdr_buffer* buf,
+ (struct htrdr* htrdr,
+ struct htrdr_buffer* buf,
const struct tile_data* tile_data)
{
struct htrdr_buffer_layout layout = HTRDR_BUFFER_LAYOUT_NULL;
@@ -165,13 +180,12 @@ write_tile_data
size_t irow_tile;
size_t ncols_tile, nrows_tile;
char* buf_mem;
- ASSERT(buf && tile_data);
+ ASSERT(htrdr && buf && tile_data);
+ (void)htrdr;
htrdr_buffer_get_layout(buf, &layout);
buf_mem = htrdr_buffer_get_data(buf);
- ASSERT(layout.elmt_size == (tile_data->format == PIXEL_LW
- ? sizeof(struct htrdr_pixel_lw)
- : sizeof(struct htrdr_pixel_sw)));
+ ASSERT(layout.elmt_size==htrdr_spectral_type_get_pixsz(htrdr->spectral_type));
/* Compute the row/column of the tile origin into the buffer */
icol = tile_data->x * (size_t)TILE_SIZE;
@@ -190,11 +204,11 @@ write_tile_data
FOR_EACH(x, 0, ncols_tile) {
switch(tile_data->format) {
- case PIXEL_LW:
- ((struct htrdr_pixel_lw*)buf_col)[x] = tile_row[x].lw;
+ case PIXEL_XWAVE:
+ ((struct htrdr_pixel_xwave*)buf_col)[x] = tile_row[x].xwave;
break;
- case PIXEL_SW:
- ((struct htrdr_pixel_sw*)buf_col)[x] = tile_row[x].sw;
+ case PIXEL_IMAGE:
+ ((struct htrdr_pixel_image*)buf_col)[x] = tile_row[x].image;
break;
default: FATAL("Unreachable code.\n"); break;
}
@@ -445,17 +459,18 @@ mpi_gather_tiles
LIST_FOR_EACH(node, tiles) {
struct tile* t = CONTAINER_OF(node, struct tile, node);
- write_tile_data(buf, &t->data);
+ write_tile_data(htrdr, buf, &t->data);
++itile;
}
if(itile != ntiles) {
+ enum pixel_format pixfmt;
ASSERT(htrdr->mpi_nprocs > 1);
/* Create a temporary tile to receive the tile data computed by the
* concurrent MPI processes */
- tile = tile_create
- (htrdr->allocator, htsky_is_long_wave(htrdr->sky) ? PIXEL_LW : PIXEL_SW);
+ pixfmt = spectral_type_to_pixfmt(htrdr->spectral_type);
+ tile = tile_create(htrdr->allocator, pixfmt);
if(!tile) {
res = RES_MEM_ERR;
htrdr_log_err(htrdr,
@@ -468,7 +483,7 @@ mpi_gather_tiles
FOR_EACH(itile, itile, ntiles) {
MPI(Recv(&tile->data, (int)msg_sz, MPI_CHAR, MPI_ANY_SOURCE,
HTRDR_MPI_TILE_DATA, MPI_COMM_WORLD, MPI_STATUS_IGNORE));
- write_tile_data(buf, &tile->data);
+ write_tile_data(htrdr, buf, &tile->data);
}
}
}
@@ -481,7 +496,7 @@ error:
}
static void
-draw_pixel_sw
+draw_pixel_image
(struct htrdr* htrdr,
const size_t ithread,
const size_t ipix[2],
@@ -489,13 +504,13 @@ draw_pixel_sw
const struct htrdr_camera* cam,
const size_t spp,
struct ssp_rng* rng,
- struct htrdr_pixel_sw* pixel)
+ struct htrdr_pixel_image* pixel)
{
struct htrdr_accum XYZ[3]; /* X, Y, and Z */
struct htrdr_accum time;
size_t ichannel;
ASSERT(ipix && ipix && pix_sz && cam && rng && pixel);
- ASSERT(!htsky_is_long_wave(htrdr->sky));
+ ASSERT(htrdr->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ);
/* Reset accumulators */
XYZ[0] = HTRDR_ACCUM_NULL;
@@ -550,7 +565,7 @@ draw_pixel_sw
(htrdr, ithread, rng, ray_org, ray_dir, wlen, iband, iquad);
ASSERT(weight >= 0);
- pdf *= 1.e9; /* Transform the pdf fro nm^-1 to m^-1 */
+ pdf *= 1.e9; /* Transform the pdf from nm^-1 to m^-1 */
weight /= pdf; /* In W/m^2/sr */
/* End the registration of the per realisation time */
@@ -576,8 +591,39 @@ draw_pixel_sw
pixel->time = time;
}
+
+static INLINE double
+radiance_temperature
+ (struct htrdr* htrdr,
+ const double radiance) /* In W/m^2/sr */
+{
+ double temperature = 0;
+ double radiance_avg = radiance;
+ res_T res = RES_OK;
+ ASSERT(htrdr && radiance >= 0);
+
+ /* From integrated radiance to average radiance in W/m^2/sr/m */
+ if(htrdr->wlen_range_m[0] != htrdr->wlen_range_m[1]) { /* !monochromatic */
+ radiance_avg /= (htrdr->wlen_range_m[1] - htrdr->wlen_range_m[0]);
+ }
+
+ res = brightness_temperature
+ (htrdr,
+ htrdr->wlen_range_m[0],
+ htrdr->wlen_range_m[1],
+ radiance_avg,
+ &temperature);
+ if(res != RES_OK) {
+ htrdr_log_warn(htrdr,
+ "Could not compute the brightness temperature for the radiance %g.\n",
+ radiance_avg);
+ temperature = 0;
+ }
+ return temperature;
+}
+
static void
-draw_pixel_lw
+draw_pixel_xwave
(struct htrdr* htrdr,
const size_t ithread,
const size_t ipix[2],
@@ -585,14 +631,15 @@ draw_pixel_lw
const struct htrdr_camera* cam,
const size_t spp,
struct ssp_rng* rng,
- struct htrdr_pixel_lw* pixel)
+ struct htrdr_pixel_xwave* pixel)
{
struct htrdr_accum radiance;
struct htrdr_accum time;
size_t isamp;
double temp_min, temp_max;
ASSERT(ipix && ipix && pix_sz && cam && rng && pixel);
- ASSERT(htsky_is_long_wave(htrdr->sky));
+ ASSERT(htrdr->spectral_type == HTRDR_SPECTRAL_LW
+ || htrdr->spectral_type == HTRDR_SPECTRAL_SW);
/* Reset the pixel accumulators */
radiance = HTRDR_ACCUM_NULL;
@@ -627,26 +674,29 @@ draw_pixel_lw
r2 = ssp_rng_canonical(rng);
/* Sample a wavelength */
- wlen = htrdr_ran_lw_sample(htrdr->ran_lw, 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);
iquad = htsky_spectral_band_sample_quadrature(htrdr->sky, r2, iband);
- /* Compute the integrated luminance in W/m^2/sr/m */
- weight = htrdr_compute_radiance_lw(htrdr, ithread, rng, ray_org, ray_dir,
- wlen, iband, iquad);
+ /* Compute the spectral radiance in W/m^2/sr/m */
+ switch(htrdr->spectral_type) {
+ case HTRDR_SPECTRAL_LW:
+ weight = htrdr_compute_radiance_lw(htrdr, ithread, rng, ray_org,
+ ray_dir, wlen, iband, iquad);
+ break;
+ case HTRDR_SPECTRAL_SW:
+ weight = htrdr_compute_radiance_sw(htrdr, ithread, rng, ray_org,
+ ray_dir, wlen, iband, iquad);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ ASSERT(weight >= 0);
/* Importance sampling: correct weight with pdf */
weight /= band_pdf; /* In W/m^2/sr */
- /* From integrated radiance to average radiance in W/m^2/sr/m */
- if(htrdr->wlen_range_m[0] != htrdr->wlen_range_m[1]) {
- /* Is not monochromatic */
- weight /= (htrdr->wlen_range_m[1] - htrdr->wlen_range_m[0]) ;
- }
- ASSERT(weight >= 0);
-
/* End the registration of the per realisation time */
time_sub(&t0, time_current(&t1), &t0);
usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
@@ -669,30 +719,13 @@ draw_pixel_lw
pixel->time = time;
/* Compute the brightness_temperature of the pixel and estimate its standard
- * error */
- #define BRIGHTNESS_TEMPERATURE(Radiance, Temperature) { \
- res_T res = brightness_temperature \
- (htrdr, \
- htrdr->wlen_range_m[0], \
- htrdr->wlen_range_m[1], \
- (Radiance), \
- &(Temperature)); \
- if(res != RES_OK) { \
- htrdr_log_warn(htrdr, \
- "Could not compute the brightness temperature for the radiance %g.\n", \
- (Radiance)); \
- (Temperature) = 0; \
- } \
- } (void)0
- BRIGHTNESS_TEMPERATURE(pixel->radiance.E, pixel->radiance_temperature.E);
- BRIGHTNESS_TEMPERATURE(pixel->radiance.E - pixel->radiance.SE, temp_min);
- BRIGHTNESS_TEMPERATURE(pixel->radiance.E + pixel->radiance.SE, temp_max);
- pixel->radiance_temperature.SE = temp_max - temp_min;
- #undef BRIGHTNESS_TEMPERATURE
-
- /* Transform the pixel radiance from W/m^2/sr/m to W/m^/sr/nm */
- pixel->radiance.E *= 1.0e-9;
- pixel->radiance.SE *= 1.0e-9;
+ * error if the sources were in the medium (<=> longwave) */
+ if(htrdr->spectral_type == HTRDR_SPECTRAL_LW) {
+ pixel->radiance_temperature.E = radiance_temperature(htrdr, pixel->radiance.E);
+ temp_min = radiance_temperature(htrdr, pixel->radiance.E - pixel->radiance.SE);
+ temp_max = radiance_temperature(htrdr, pixel->radiance.E + pixel->radiance.SE);
+ pixel->radiance_temperature.SE = temp_max - temp_min;
+ }
}
static res_T
@@ -734,10 +767,15 @@ draw_tile
ipix[1] = tile_org[1] + ipix_tile[1];
/* Draw the pixel */
- if(htsky_is_long_wave(htrdr->sky)) {
- draw_pixel_lw(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->lw);
- } else {
- draw_pixel_sw(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->sw);
+ switch(htrdr->spectral_type) {
+ case HTRDR_SPECTRAL_LW:
+ case HTRDR_SPECTRAL_SW:
+ draw_pixel_xwave(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->xwave);
+ break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ draw_pixel_image(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->image);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
}
}
return RES_OK;
@@ -783,7 +821,7 @@ draw_image
htrdr->mpi_working_procs[htrdr->mpi_rank] = 0;
--htrdr->mpi_nworking_procs;
- pixfmt = htsky_is_long_wave(htrdr->sky) ? PIXEL_LW : PIXEL_SW;
+ pixfmt = spectral_type_to_pixfmt(htrdr->spectral_type);
omp_set_num_threads((int)nthreads);
#pragma omp parallel
@@ -937,19 +975,13 @@ htrdr_draw_radiance
proc_work_init(htrdr->allocator, &work);
if(htrdr->mpi_rank == 0) {
- size_t pixsz, pixal;
+ const size_t pixsz = htrdr_spectral_type_get_pixsz(htrdr->spectral_type);
+ const size_t pixal = htrdr_spectral_type_get_pixal(htrdr->spectral_type);
+
htrdr_buffer_get_layout(buf, &layout);
ASSERT(layout.width || layout.height || layout.elmt_size);
ASSERT(layout.width == width && layout.height == height);
- if(htsky_is_long_wave(htrdr->sky)) {
- pixsz = sizeof(struct htrdr_pixel_lw);
- pixal = ALIGNOF(struct htrdr_pixel_lw);
- } else {
- pixsz = sizeof(struct htrdr_pixel_sw);
- pixal = ALIGNOF(struct htrdr_pixel_sw);
- }
-
if(layout.elmt_size != pixsz || layout.alignment < pixal) {
htrdr_log_err(htrdr, "%s: invalid buffer layout.\n", FUNC_NAME);
res = RES_BAD_ARG;
diff --git a/src/htrdr_ran_lw.c b/src/htrdr_ran_lw.c
@@ -1,363 +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_ran_lw.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_lw {
- 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_ran_lw_cdf
- (struct htrdr_ran_lw* ran_lw,
- const char* func_name)
-{
- double* pdf = NULL;
- double* cdf = NULL;
- double sum = 0;
- size_t i;
- res_T res = RES_OK;
- ASSERT(ran_lw && func_name && ran_lw->nbands != HTRDR_RAN_LW_CONTINUE);
-
- res = darray_double_resize(&ran_lw->cdf, ran_lw->nbands);
- if(res != RES_OK) {
- htrdr_log_err(ran_lw->htrdr,
- "%s: Error allocating the CDF of the longwave spectral bands -- %s.\n",
- func_name, res_to_cstr(res));
- goto error;
- }
- res = darray_double_resize(&ran_lw->pdf, ran_lw->nbands);
- if(res != RES_OK) {
- htrdr_log_err(ran_lw->htrdr,
- "%s: Error allocating the pdf of the longwave spectral bands -- %s.\n",
- func_name, res_to_cstr(res));
- goto error;
- }
-
- cdf = darray_double_data_get(&ran_lw->cdf);
- pdf = darray_double_data_get(&ran_lw->pdf); /* Now save pdf to correct weight */
-
- htrdr_log(ran_lw->htrdr,
- "Number of bands of the longwave cumulative: %lu\n",
- (unsigned long)ran_lw->nbands);
-
- /* Compute the *unormalized* probability to sample a longwave band */
- FOR_EACH(i, 0, ran_lw->nbands) {
- double lambda_lo = ran_lw->range[0] + (double)i * ran_lw->band_len;
- double lambda_hi = MMIN(lambda_lo + ran_lw->band_len, ran_lw->range[1]);
- ASSERT(lambda_lo<= lambda_hi);
- ASSERT(lambda_lo > ran_lw->range[0] || eq_eps(lambda_lo, ran_lw->range[0], 1.e-6));
- ASSERT(lambda_lo < ran_lw->range[1] || eq_eps(lambda_lo, ran_lw->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, ran_lw->ref_temperature);
-
- /* Update the norm */
- sum += pdf[i];
- }
-
- /* Compute the cumulative of the previously computed probabilities */
- FOR_EACH(i, 0, ran_lw->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[ran_lw->nbands-1], 1, 1.e-6));
- cdf[ran_lw->nbands - 1] = 1.0; /* Handle numerical issue */
-
-exit:
- return res;
-error:
- darray_double_clear(&ran_lw->cdf);
- darray_double_clear(&ran_lw->pdf);
- goto exit;
-}
-
-static double
-ran_lw_sample_continue
- (const struct htrdr_ran_lw* ran_lw,
- 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(ran_lw && 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], ran_lw->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, ran_lw->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(ran_lw->htrdr,
- "%s: could not sample a wavelength in the range [%g, %g] nanometers "
- "for the reference temperature %g Kelvin.\n",
- func_name, SPLIT2(range), ran_lw->ref_temperature);
- }
-
- if(pdf) {
- const double Tref = ran_lw->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
-ran_lw_sample_discrete
- (const struct htrdr_ran_lw* ran_lw,
- 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(ran_lw && ran_lw->nbands != HTRDR_RAN_LW_CONTINUE);
- ASSERT(0 <= r0 && r0 < 1);
- ASSERT(0 <= r1 && r1 < 1);
- (void)func_name;
- (void)pdf_band;
-
- cdf = darray_double_cdata_get(&ran_lw->cdf);
- cdf_length = darray_double_size_get(&ran_lw->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] = ran_lw->range[0] + (double)i*ran_lw->band_len;
- band_range[1] = band_range[0] + ran_lw->band_len;
-
- /* Fetch the pdf of the sampled band */
- pdf_band = darray_double_cdata_get(&ran_lw->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_ran_lw(ref_T* ref)
-{
- struct htrdr_ran_lw* ran_lw = NULL;
- ASSERT(ref);
- ran_lw = CONTAINER_OF(ref, struct htrdr_ran_lw, ref);
- darray_double_release(&ran_lw->cdf);
- darray_double_release(&ran_lw->pdf);
- MEM_RM(ran_lw->htrdr->allocator, ran_lw);
-}
-
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-res_T
-htrdr_ran_lw_create
- (struct htrdr* htrdr,
- const double range[2], /* Must be included in [1000, 100000] nanometers */
- const size_t nbands, /* # bands used to discretized CDF */
- const double ref_temperature,
- struct htrdr_ran_lw** out_ran_lw)
-{
- struct htrdr_ran_lw* ran_lw = NULL;
- double min_band_len = 0;
- res_T res = RES_OK;
- ASSERT(htrdr && range && out_ran_lw && ref_temperature > 0);
- ASSERT(ref_temperature > 0);
- ASSERT(range[0] >= 1000);
- ASSERT(range[1] <= 100000);
- ASSERT(range[0] <= range[1]);
-
- ran_lw = MEM_CALLOC(htrdr->allocator, 1, sizeof(*ran_lw));
- if(!ran_lw) {
- 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(&ran_lw->ref);
- ran_lw->htrdr = htrdr;
- darray_double_init(htrdr->allocator, &ran_lw->cdf);
- darray_double_init(htrdr->allocator, &ran_lw->pdf);
-
- ran_lw->range[0] = range[0];
- ran_lw->range[1] = range[1];
- ran_lw->ref_temperature = ref_temperature;
- ran_lw->nbands = nbands;
-
- min_band_len = compute_sky_min_band_len(ran_lw->htrdr->sky, ran_lw->range);
-
- if(nbands == HTRDR_RAN_LW_CONTINUE) {
- ran_lw->band_len = 0;
- } else {
- ran_lw->band_len = (range[1] - range[0]) / (double)ran_lw->nbands;
-
- /* Adjust the band length to ensure that each sky spectral interval is
- * overlapped by at least one band */
- if(ran_lw->band_len > min_band_len) {
- ran_lw->band_len = min_band_len;
- ran_lw->nbands = (size_t)ceil((range[1] - range[0]) / ran_lw->band_len);
- }
-
- res = setup_ran_lw_cdf(ran_lw, FUNC_NAME);
- if(res != RES_OK) goto error;
- }
-
- htrdr_log(htrdr, "Long wave interval defined on [%g, %g] nanometers.\n",
- range[0], range[1]);
-
-exit:
- *out_ran_lw = ran_lw;
- return res;
-error:
- if(ran_lw) htrdr_ran_lw_ref_put(ran_lw);
- goto exit;
-}
-
-void
-htrdr_ran_lw_ref_get(struct htrdr_ran_lw* ran_lw)
-{
- ASSERT(ran_lw);
- ref_get(&ran_lw->ref);
-}
-
-void
-htrdr_ran_lw_ref_put(struct htrdr_ran_lw* ran_lw)
-{
- ASSERT(ran_lw);
- ref_put(&ran_lw->ref, release_ran_lw);
-}
-
-double
-htrdr_ran_lw_sample
- (const struct htrdr_ran_lw* ran_lw,
- const double r0,
- const double r1,
- double* pdf)
-{
- ASSERT(ran_lw);
- if(ran_lw->nbands != HTRDR_RAN_LW_CONTINUE) { /* Discrete */
- return ran_lw_sample_discrete(ran_lw, r0, r1, FUNC_NAME, pdf);
- } else if(eq_eps(ran_lw->range[0], ran_lw->range[1], 1.e-6)) {
- if(pdf) *pdf = 1;
- return ran_lw->range[0];
- } else { /* Continue */
- return ran_lw_sample_continue
- (ran_lw, r0, ran_lw->range, FUNC_NAME, pdf);
- }
-}
-
diff --git a/src/htrdr_ran_lw.h b/src/htrdr_ran_lw.h
@@ -1,53 +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_RAN_LW_H
-#define HTRDR_RAN_LW_H
-
-#include <rsys/rsys.h>
-
-#define HTRDR_RAN_LW_CONTINUE 0
-
-struct htrdr;
-struct htrdr_ran_lw;
-
-extern LOCAL_SYM res_T
-htrdr_ran_lw_create
- (struct htrdr* htrdr,
- const double range[2], /* Must be included in [1000, 100000] nanometers */
- /* # bands used to discretisze the LW domain. HTRDR_RAN_LW_CONTINUE <=> no
- * discretisation */
- const size_t nbands, /* Hint on #bands used to discretised the CDF */
- const double ref_temperature, /* Reference temperature */
- struct htrdr_ran_lw** ran_lw);
-
-extern LOCAL_SYM void
-htrdr_ran_lw_ref_get
- (struct htrdr_ran_lw* ran_lw);
-
-extern LOCAL_SYM void
-htrdr_ran_lw_ref_put
- (struct htrdr_ran_lw* ran_lw);
-
-/* Return a wavelength in nanometer */
-extern LOCAL_SYM double
-htrdr_ran_lw_sample
- (const struct htrdr_ran_lw* ran_lw,
- const double r0, /* Canonical number in [0, 1[ */
- const double r1, /* Canonical number in [0, 1[ */
- double* pdf); /* May be NULL */
-
-#endif /* HTRDR_RAN_LW_H */
-
diff --git a/src/htrdr_ran_wlen.c b/src/htrdr_ran_wlen.c
@@ -0,0 +1,364 @@
+/* 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] <= 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,54 @@
+/* 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,
+ 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_solve.h b/src/htrdr_solve.h
@@ -35,20 +35,20 @@ struct htrdr_estimate {
};
static const struct htrdr_estimate HTRDR_ESTIMATE_NULL;
-struct htrdr_pixel_sw {
- struct htrdr_estimate X; /* In W/m^2/sr */
- struct htrdr_estimate Y; /* In W/m^2/sr */
- struct htrdr_estimate Z; /* In W/m^2/sr */
+struct htrdr_pixel_xwave {
+ struct htrdr_estimate radiance; /* In W/m^2/sr */
+ struct htrdr_estimate radiance_temperature; /* In K */
struct htrdr_accum time; /* In microseconds */
};
-static const struct htrdr_pixel_sw HTRDR_PIXEL_SW_NULL;
+static const struct htrdr_pixel_xwave HTRDR_PIXEL_XWAVE_NULL;
-struct htrdr_pixel_lw {
- struct htrdr_estimate radiance; /* In W/m^2/sr */
- struct htrdr_estimate radiance_temperature; /* In K */
+struct htrdr_pixel_image {
+ struct htrdr_estimate X; /* In W/m^2/sr */
+ struct htrdr_estimate Y; /* In W/m^2/sr */
+ struct htrdr_estimate Z; /* In W/m^2/sr */
struct htrdr_accum time; /* In microseconds */
};
-static const struct htrdr_pixel_lw HTRDR_PIXEL_LW_NULL;
+static const struct htrdr_pixel_image HTRDR_PIXEL_IMAGE_NULL;
/* Forward declarations */
struct htrdr;
@@ -120,4 +120,38 @@ htrdr_accum_get_estimation
}
}
+static INLINE size_t
+htrdr_spectral_type_get_pixsz(const enum htrdr_spectral_type type)
+{
+ size_t sz = 0;
+ switch(type) {
+ case HTRDR_SPECTRAL_LW:
+ case HTRDR_SPECTRAL_SW:
+ sz = sizeof(struct htrdr_pixel_xwave);
+ break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ sz = sizeof(struct htrdr_pixel_image);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return sz;
+}
+
+static INLINE size_t
+htrdr_spectral_type_get_pixal(const enum htrdr_spectral_type type)
+{
+ size_t al = 0;
+ switch(type) {
+ case HTRDR_SPECTRAL_LW:
+ case HTRDR_SPECTRAL_SW:
+ al = ALIGNOF(struct htrdr_pixel_xwave);
+ break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ al = ALIGNOF(struct htrdr_pixel_image);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return al;
+}
+
#endif /* HTRDR_SOLVE_H */
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 */
