commit c6b3b5ace0ebfdf16ec4a1da54a222344d1c2360
parent d5a0b5dd94ab262f2af18291f11b0c60b9eefc6d
Author: Najda Villefranque <najda.villefranque@gmail.com>
Date: Wed, 27 May 2020 18:59:32 +0200
First shot at SW integration in non-image mode
Argument -s WLEN_MIN,WLEN_MAX same as -l but for solar source
Spectral sampling in SW with planck function at Tref=Tsun=5778 K
(*ran_lw* replaced by *wlen_ran*)
Additional member in struct htrdr : "is_image"
pixel type is chosen accordingly
pixel_image or pixel_integ
instead of pixel_sw and pixel_lw
if is_image, draw_pixel_image (compute_radiance_sw only)
else, draw_pixel_integ :
if is_longwave : compute_radiance_lw
else : compute_radiance_sw
Tbrightness computed only if longwave
L_sun is planck_monochromatic at wlen, Tsun in compute_radiance_sw
Removed asserts that controled integration range
Works with modified version of htsky
Diffstat:
13 files changed, 599 insertions(+), 518 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -100,7 +100,7 @@ set(HTRDR_FILES_SRC
htrdr_interface.c
htrdr_main.c
htrdr_mtl.c
- htrdr_ran_lw.c
+ htrdr_wlen_ran.c
htrdr_rectangle.c
htrdr_slab.c
htrdr_sun.c)
@@ -114,7 +114,7 @@ set(HTRDR_FILES_INC
htrdr_ground.h
htrdr_interface.h
htrdr_mtl.h
- htrdr_ran_lw.h
+ htrdr_wlen_ran.h
htrdr_rectangle.h
htrdr_slab.h
htrdr_sun.h
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_wlen_ran.h"
#include "htrdr_sun.h"
#include "htrdr_solve.h"
@@ -117,12 +117,12 @@ 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);
+ if(!htrdr->is_image) {
+ pixsz = sizeof(struct htrdr_pixel_integ);
+ pixal = ALIGNOF(struct htrdr_pixel_integ);
} else {
- pixsz = sizeof(struct htrdr_pixel_sw);
- pixal = ALIGNOF(struct htrdr_pixel_sw);
+ pixsz = sizeof(struct htrdr_pixel_image);
+ pixal = ALIGNOF(struct htrdr_pixel_image);
}
htrdr_buffer_get_layout(buf, &layout);
@@ -140,8 +140,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->is_image){
+ const struct htrdr_pixel_integ* 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 +149,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);
@@ -484,12 +484,28 @@ 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];
+ /* should the sky load short or long wave data ? */
+ /* if longwave is degenerated => sw ; else : lw */
+ if(args->wlen_lw_range[0] > args->wlen_lw_range[1]) {
+ htsky_args.is_long_wave = 0 ;
+ htsky_args.wlen_range[0] = args->wlen_sw_range[0];
+ htsky_args.wlen_range[1] = args->wlen_sw_range[1];
+ htrdr->is_image=0;
+ } else {
+ htsky_args.is_long_wave = 1 ;
+ htsky_args.wlen_range[0] = args->wlen_lw_range[0];
+ htsky_args.wlen_range[1] = args->wlen_lw_range[1];
+ if(args->wlen_sw_range[0] > args->wlen_sw_range[1]) { /* image */
+ htrdr->is_image = 1 ;
+ } else {
+ htrdr->is_image = 0 ;
+ }
+ }
+
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 */
+ if(htrdr->is_image) {
const double* range = HTRDR_CIE_XYZ_RANGE_DEFAULT;
size_t n;
@@ -497,19 +513,34 @@ htrdr_init
res = htrdr_cie_xyz_create(htrdr, range, n, &htrdr->cie);
if(res != RES_OK) goto error;
- } else { /* Long Wave random variate */
- const double Tref = 290; /* In Kelvin */
- size_t n;
+ } else {
+ if(htsky_is_long_wave(htrdr->sky)) { /* Long wave random variate */
+ const double Tref=290 ; /* In Kelvin */
+ size_t n;
- htrdr->wlen_range_m[0] = args->wlen_lw_range[0]*1e-9; /* Convert in meters */
- htrdr->wlen_range_m[1] = args->wlen_lw_range[1]*1e-9; /* Convert in meters */
- ASSERT(htrdr->wlen_range_m[0] <= htrdr->wlen_range_m[1]);
+ htrdr->wlen_range_m[0] = args->wlen_lw_range[0]*1e-9; /* Convert in meters */
+ htrdr->wlen_range_m[1] = args->wlen_lw_range[1]*1e-9; /* Convert in meters */
+ ASSERT(htrdr->wlen_range_m[0] <= htrdr->wlen_range_m[1]);
+ n = (size_t)(args->wlen_lw_range[1] - args->wlen_lw_range[0]);
- 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);
- if(res != RES_OK) goto error;
- }
+ res = htrdr_wlen_ran_create
+ (htrdr, args->wlen_lw_range, n, Tref, &htrdr->wlen_ran);
+ if(res != RES_OK) goto error;
+
+ } else {
+ const double Tref=5778 ; /* Tsun In Kelvin */
+ size_t n;
+
+ htrdr->wlen_range_m[0] = args->wlen_sw_range[0]*1e-9; /* Convert in meters */
+ htrdr->wlen_range_m[1] = args->wlen_sw_range[1]*1e-9; /* Convert in meters */
+ ASSERT(htrdr->wlen_range_m[0] <= htrdr->wlen_range_m[1]);
+ n = (size_t)(args->wlen_sw_range[1] - args->wlen_sw_range[0]);
+
+ res = htrdr_wlen_ran_create
+ (htrdr, args->wlen_sw_range, n, Tref, &htrdr->wlen_ran);
+ if(res != RES_OK) goto error;
+ }
+ }
htrdr->lifo_allocators = MEM_CALLOC
(htrdr->allocator, htrdr->nthreads, sizeof(*htrdr->lifo_allocators));
@@ -538,12 +569,12 @@ htrdr_init
size_t pixsz = 0; /* sizeof(pixel) */
size_t pixal = 0; /* alignof(pixel) */
- if(htsky_is_long_wave(htrdr->sky)) {
- pixsz = sizeof(struct htrdr_pixel_lw);
- pixal = ALIGNOF(struct htrdr_pixel_lw);
+ if(!htrdr->is_image) {
+ pixsz = sizeof(struct htrdr_pixel_integ);
+ pixal = ALIGNOF(struct htrdr_pixel_integ);
} else {
- pixsz = sizeof(struct htrdr_pixel_sw);
- pixal = ALIGNOF(struct htrdr_pixel_sw);
+ pixsz = sizeof(struct htrdr_pixel_image);
+ pixal = ALIGNOF(struct htrdr_pixel_image);
}
res = htrdr_buffer_create(htrdr,
args->image.definition[0], /* Width */
@@ -575,7 +606,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->wlen_ran) htrdr_wlen_ran_ref_put(htrdr->wlen_ran);
if(htrdr->output && htrdr->output != stdout) fclose(htrdr->output);
if(htrdr->lifo_allocators) {
size_t i;
diff --git a/src/htrdr.h b/src/htrdr.h
@@ -50,13 +50,14 @@ struct htrdr {
struct htrdr_mtl* mtl;
struct htrdr_sun* sun;
struct htrdr_cie_xyz* cie;
- struct htrdr_ran_lw* ran_lw;
+ struct htrdr_wlen_ran* wlen_ran;
struct htrdr_camera* cam;
struct htrdr_buffer* buf;
struct htsky* sky;
double wlen_range_m[2]; /* Integration range in *meters* */
+ int is_image ; /* XYZ Image or spectral integration? */
size_t spp; /* #samples per pixel */
size_t width; /* Image width */
diff --git a/src/htrdr_args.c b/src/htrdr_args.c
@@ -67,10 +67,6 @@ print_help(const char* cmd)
" 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 +77,16 @@ 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 WLEN_MIN,WLEN_MAX\n"
+" switch in solar rendering for the short waves 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 with CIE.\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);
@@ -366,23 +372,23 @@ error:
}
static res_T
-parse_lw_range(struct htrdr_args* args, const char* str)
+parse_spectral_range(const char* str, double wlen_range[2])
{
double range[2];
size_t len;
res_T res = RES_OK;
- ASSERT(args && str);
+ 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 longwave range `%s'.\n", str);
+ fprintf(stderr, "Invalid spectral range `%s'.\n", str);
goto error;
}
- args->wlen_lw_range[0] = range[0];
- args->wlen_lw_range[1] = range[1];
+ wlen_range[0] = range[0];
+ wlen_range[1] = range[1];
exit:
return res;
@@ -390,6 +396,28 @@ error:
goto exit;
}
+static res_T
+parse_lw_range(struct htrdr_args* args, const char* str)
+{
+ res_T res = RES_OK;
+ ASSERT(args && str);
+
+ res = parse_spectral_range(str, args->wlen_lw_range) ;
+
+ return res ;
+}
+
+static res_T
+parse_sw_range(struct htrdr_args* args, const char* str)
+{
+ res_T res = RES_OK;
+ ASSERT(args && str);
+
+ res = parse_spectral_range(str, args->wlen_sw_range) ;
+
+ return res ;
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -414,7 +442,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:l:M:m:O:o:Rrs:T:t:V:v")) != -1) {
switch(opt) {
case 'a': args->filename_gas = optarg; break;
case 'C':
@@ -444,6 +472,9 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
case 'o': args->output = optarg; break;
case 'r': args->repeat_clouds = 1; break;
case 'R': args->repeat_ground = 1; break;
+ case 's':
+ res = parse_sw_range(args, optarg);
+ break;
case 'T':
res = cstr_to_double(optarg, &args->optical_thickness);
if(res == RES_OK && args->optical_thickness < 0) res = RES_BAD_ARG;
diff --git a/src/htrdr_args.h.in b/src/htrdr_args.h.in
@@ -55,7 +55,8 @@ 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 */
+ double wlen_sw_range[2]; /* Spectral range of integration in nm (solar) */
+ double wlen_lw_range[2]; /* Spectral range of integration in nm (thermal) */
unsigned nthreads; /* Hint on the number of threads to use */
int force_overwriting;
@@ -92,7 +93,8 @@ 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 */ \
+ {DBL_MAX, -DBL_MAX}, /* Shortwave range. Degenerated <=> CIE sw */ \
+ {DBL_MAX, -DBL_MAX}, /* Long wave range. Degenerated <=> CIE sw */ \
(unsigned)~0, /* #threads */ \
0, /* Force overwriting */ \
0, /* dump VTK */ \
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
@@ -275,6 +275,8 @@ htrdr_compute_radiance_sw
double L_sun; /* Sun radiance in W.m^-2.sr^-1 */
double sun_dir[3];
double ksi = 1; /* Throughput */
+ double temperature=5778;
+ double wlen_m = wlen * 1.e-9;
double w = 0; /* MC weight */
double g = 0; /* Asymmetry parameter of the HG phase function */
@@ -297,8 +299,8 @@ htrdr_compute_radiance_sw
htsky_get_spectral_band_bounds(htrdr->sky, iband, band_bounds);
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);
-
+ /* L_sun = htrdr_sun_get_radiance(htrdr->sun, wlen);*/
+ L_sun = planck_monochromatic(wlen_m, temperature);
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_wlen_ran.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_INTEG,
+ PIXEL_IMAGE
};
union pixel {
- struct htrdr_pixel_lw lw;
- struct htrdr_pixel_sw sw;
+ struct htrdr_pixel_integ integ;
+ struct htrdr_pixel_image image;
};
/* Tile of row ordered image pixels */
@@ -169,9 +169,9 @@ write_tile_data
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 == (tile_data->format == PIXEL_INTEG
+ ? sizeof(struct htrdr_pixel_integ)
+ : sizeof(struct htrdr_pixel_image)));
/* Compute the row/column of the tile origin into the buffer */
icol = tile_data->x * (size_t)TILE_SIZE;
@@ -190,11 +190,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_INTEG:
+ ((struct htrdr_pixel_integ*)buf_col)[x] = tile_row[x].integ;
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;
}
@@ -455,7 +455,7 @@ mpi_gather_tiles
/* 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);
+ (htrdr->allocator, htrdr->is_image ? PIXEL_IMAGE : PIXEL_INTEG) ;
if(!tile) {
res = RES_MEM_ERR;
htrdr_log_err(htrdr,
@@ -481,7 +481,7 @@ error:
}
static void
-draw_pixel_sw
+draw_pixel_image
(struct htrdr* htrdr,
const size_t ithread,
const size_t ipix[2],
@@ -489,7 +489,7 @@ 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;
@@ -577,7 +577,7 @@ draw_pixel_sw
}
static void
-draw_pixel_lw
+draw_pixel_integ
(struct htrdr* htrdr,
const size_t ithread,
const size_t ipix[2],
@@ -585,14 +585,14 @@ draw_pixel_lw
const struct htrdr_camera* cam,
const size_t spp,
struct ssp_rng* rng,
- struct htrdr_pixel_lw* pixel)
+ struct htrdr_pixel_integ* 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->is_image));
/* Reset the pixel accumulators */
radiance = HTRDR_ACCUM_NULL;
@@ -627,15 +627,20 @@ 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_wlen_ran_sample(htrdr->wlen_ran, 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);
+ if(htsky_is_long_wave(htrdr->sky)) {
+ weight = htrdr_compute_radiance_lw(htrdr, ithread, rng, ray_org, ray_dir,
+ wlen, iband, iquad);
+ } else {
+ weight = htrdr_compute_radiance_sw(htrdr, ithread, rng, ray_org, ray_dir,
+ wlen, iband, iquad);
+ }
/* Importance sampling: correct weight with pdf */
weight /= band_pdf; /* In W/m^2/sr */
@@ -669,26 +674,28 @@ 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
+ * error if the sources were in the medium (i.e., is_long_wave) */
+ if(htsky_is_long_wave(htrdr->sky)) {
+ #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;
@@ -734,10 +741,10 @@ 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);
+ if(!(htrdr->is_image)) {
+ draw_pixel_integ(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->integ);
} else {
- draw_pixel_sw(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->sw);
+ draw_pixel_image(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->image);
}
}
return RES_OK;
@@ -783,7 +790,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 = htrdr->is_image ? PIXEL_IMAGE : PIXEL_INTEG ;
omp_set_num_threads((int)nthreads);
#pragma omp parallel
@@ -942,12 +949,12 @@ htrdr_draw_radiance
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);
+ if(!(htrdr->is_image)) {
+ pixsz = sizeof(struct htrdr_pixel_integ);
+ pixal = ALIGNOF(struct htrdr_pixel_integ);
} else {
- pixsz = sizeof(struct htrdr_pixel_sw);
- pixal = ALIGNOF(struct htrdr_pixel_sw);
+ pixsz = sizeof(struct htrdr_pixel_image);
+ pixal = ALIGNOF(struct htrdr_pixel_image);
}
if(layout.elmt_size != pixsz || layout.alignment < pixal) {
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_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_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_sw HTRDR_PIXEL_SW_NULL;
+static const struct htrdr_pixel_image HTRDR_PIXEL_IMAGE_NULL;
-struct htrdr_pixel_lw {
+struct htrdr_pixel_integ {
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_lw HTRDR_PIXEL_LW_NULL;
+static const struct htrdr_pixel_integ HTRDR_PIXEL_INTEG_NULL;
/* Forward declarations */
struct htrdr;
diff --git a/src/htrdr_wlen_ran.c b/src/htrdr_wlen_ran.c
@@ -0,0 +1,365 @@
+/* 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_wlen_ran.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_wlen_ran {
+ 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_wlen_ran* 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_wlen_ran* 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_wlen_ran* 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_wlen_ran* wlen_ran = NULL;
+ ASSERT(ref);
+ wlen_ran = CONTAINER_OF(ref, struct htrdr_wlen_ran, 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_wlen_ran_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_wlen_ran** out_wlen_ran)
+{
+ struct htrdr_wlen_ran* 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] >= 200);
+ ASSERT(range[1] <= 100000);
+ 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_wlen_ran_ref_put(wlen_ran);
+ goto exit;
+}
+
+void
+htrdr_wlen_ran_ref_get(struct htrdr_wlen_ran* wlen_ran)
+{
+ ASSERT(wlen_ran);
+ ref_get(&wlen_ran->ref);
+}
+
+void
+htrdr_wlen_ran_ref_put(struct htrdr_wlen_ran* wlen_ran)
+{
+ ASSERT(wlen_ran);
+ ref_put(&wlen_ran->ref, release_wlen_ran);
+}
+
+double
+htrdr_wlen_ran_sample
+ (const struct htrdr_wlen_ran* 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_wlen_ran.h b/src/htrdr_wlen_ran.h
@@ -0,0 +1,59 @@
+/* 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_WLEN_RAN_H
+#define HTRDR_WLEN_RAN_H
+
+#include <rsys/rsys.h>
+
+#define HTRDR_WLEN_RAN_CONTINUE 0
+#define HTRDR_WLEN_RAN_SOLAR_WVN_MIN 820 # 12195 nm
+#define HTRDR_WLEN_RAN_SOLAR_WVN_MAX 50000 # 200 nm
+#define HTRDR_WLEN_RAN_THERMAL_WVN_MIN 10 # 1000000 nm
+#define HTRDR_WLEN_RAN_THERMAL_WVN_MAX 3250 # 3077 nm
+
+struct htrdr;
+struct htrdr_wlen_ran;
+
+extern LOCAL_SYM res_T
+htrdr_wlen_ran_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 */
+ const size_t nbands, /* Hint on #bands used to discretised th CDF */
+ const double ref_temperature, /* Reference temperature */
+ struct htrdr_wlen_ran** wlen_ran);
+
+extern LOCAL_SYM void
+htrdr_wlen_ran_ref_get
+ (struct htrdr_wlen_ran* wlen_ran);
+
+extern LOCAL_SYM void
+htrdr_wlen_ran_ref_put
+ (struct htrdr_wlen_ran* wlen_ran);
+
+/* Return a wavelength in nanometer */
+extern LOCAL_SYM double
+htrdr_wlen_ran_sample
+ (const struct htrdr_wlen_ran* 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_WLEN_RAN_H */
+
