commit f696b9bd29b4d34ca4630b30d6309425cc668b03
parent e5488d977315421230a3b006bc58571ca9684b70
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 29 May 2020 16:57:40 +0200
In longwave strictly use the sky spectral bands data
Diffstat:
5 files changed, 126 insertions(+), 29 deletions(-)
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -499,15 +499,16 @@ htrdr_init
} else { /* Long Wave random variate */
const double Tref = 290; /* In Kelvin */
+ double sky_range[2];
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 */
+ HTSKY(get_spectral_bounds(htrdr->sky, sky_range));
+ htrdr->wlen_range_m[0] = sky_range[0]*1e-9; /* Convert in meters */
+ htrdr->wlen_range_m[1] = sky_range[1]*1e-9; /* Convert in meters */
ASSERT(htrdr->wlen_range_m[0] <= htrdr->wlen_range_m[1]);
n = (size_t)(args->wlen_lw_range[1] - args->wlen_lw_range[0]);
- res = htrdr_ran_lw_create
- (htrdr, args->wlen_lw_range, n, Tref, &htrdr->ran_lw);
+ res = htrdr_ran_lw_create(htrdr, n, Tref, &htrdr->ran_lw);
if(res != RES_OK) goto error;
}
diff --git a/src/htrdr_compute_radiance_lw.c b/src/htrdr_compute_radiance_lw.c
@@ -154,7 +154,9 @@ htrdr_compute_radiance_lw
double pos_next[3];
double dir_next[3];
double temperature;
- double wlen_m = wlen * 1.e-9;
+ double sky_band_bounds[2];
+ double sky_band_bounds_m[2];
+ double sky_band_len_m;
double g;
double w = 0; /* Weight */
@@ -163,10 +165,16 @@ htrdr_compute_radiance_lw
/* Setup the phase function for this spectral band & quadrature point */
CHK(RES_OK == ssf_phase_create
(&htrdr->lifo_allocators[ithread], &ssf_phase_hg, &phase_hg));
- g = htsky_fetch_per_wavelength_particle_phase_function_asymmetry_parameter
- (htrdr->sky, wlen);
+ g = htsky_fetch_particle_phase_function_asymmetry_parameter
+ (htrdr->sky, iband, iquad);
SSF(phase_hg_setup(phase_hg, g));
+ /* Fetch the boundaries of the sampled sky band */
+ HTSKY(get_spectral_band_bounds(htrdr->sky, iband, sky_band_bounds));
+ sky_band_bounds_m[0] = sky_band_bounds[0] * 1.e-9;
+ sky_band_bounds_m[1] = sky_band_bounds[1] * 1.e-9;
+ sky_band_len_m = sky_band_bounds_m[1] - sky_band_bounds_m[0];
+
/* Initialise the random walk */
d3_set(pos, pos_in);
d3_set(dir, dir_in);
@@ -213,7 +221,8 @@ htrdr_compute_radiance_lw
ASSERT(!SVX_HIT_NONE(&svx_hit));
temperature = htsky_fetch_temperature(htrdr->sky, pos_next);
/* weight is planck integrated over the spectral sub-interval */
- w = planck_monochromatic(wlen_m, temperature);
+ w = planck(sky_band_bounds_m[0], sky_band_bounds_m[1], temperature);
+ w = w * sky_band_len_m;
break;
}
@@ -262,7 +271,8 @@ htrdr_compute_radiance_lw
w = 0;
} else {
/* weight is planck integrated over the spectral sub-interval */
- w = planck_monochromatic(wlen_m, temperature);
+ w = planck(sky_band_bounds_m[0], sky_band_bounds_m[1], temperature);
+ w = w * sky_band_len_m;
}
break;
}
diff --git a/src/htrdr_draw_radiance.c b/src/htrdr_draw_radiance.c
@@ -627,24 +627,24 @@ 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_lw_sample(htrdr->ran_lw, r0, r1, NULL);
/* 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 */
+ /* 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);
- /* Importance sampling: correct weight with pdf */
+ /* Importance sampling: correct weight with pdf
+ * W/m^2/sr/m => W/m^2/sr */
+ band_pdf = htrdr_ran_lw_get_sky_band_pdf(htrdr->ran_lw, iband);
weight /= band_pdf;
- /* 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]) ;
- }
+ /* Is not monochromatic */
+ ASSERT(htrdr->wlen_range_m[0] != htrdr->wlen_range_m[1]);
+ weight /= (htrdr->wlen_range_m[1] - htrdr->wlen_range_m[0]) ;
ASSERT(weight >= 0);
/* End the registration of the per realisation time */
diff --git a/src/htrdr_ran_lw.c b/src/htrdr_ran_lw.c
@@ -30,6 +30,7 @@
#include <math.h> /* nextafter */
struct htrdr_ran_lw {
+ struct darray_double pdf_sky_bands;
struct darray_double pdf;
struct darray_double cdf;
double range[2]; /* Boundaries of the spectral integration interval */
@@ -121,6 +122,69 @@ error:
goto exit;
}
+static res_T
+compute_sky_bands_pdf(struct htrdr_ran_lw* ran_lw, const char* func_name)
+{
+ double* pdf = NULL;
+ size_t nbands;
+ res_T res = RES_OK;
+ ASSERT(ran_lw && htsky_is_long_wave(ran_lw->htrdr->sky) && func_name);
+
+ nbands = htsky_get_spectral_bands_count(ran_lw->htrdr->sky);
+
+ res = darray_double_resize(&ran_lw->pdf_sky_bands, nbands);
+ if(res != RES_OK) {
+ htrdr_log_err(ran_lw->htrdr,
+ "%s: error allocating the PDF of the long wave spectral bands "
+ "of the sky -- %s.\n",
+ func_name, res_to_cstr(res));
+ goto error;
+ }
+
+ pdf = darray_double_data_get(&ran_lw->pdf_sky_bands);
+
+ if(eq_eps(ran_lw->range[0], ran_lw->range[1], 1.e-6)) {
+ ASSERT(nbands == 1);
+ pdf[0] = 1;
+ } else {
+ double sum = 0;
+ size_t i = 0;
+
+ /* Compute the *unormalized* probability to sample a long wave band */
+ FOR_EACH(i, 0, nbands) {
+ const size_t iband = htsky_get_spectral_band_id(ran_lw->htrdr->sky, i);
+ double wlens[2];
+
+ HTSKY(get_spectral_band_bounds(ran_lw->htrdr->sky, iband, wlens));
+ /* Ensure that the whole sky band is included into the random variate range */
+ ASSERT(wlens[0] >= ran_lw->range[0]);
+ ASSERT(wlens[1] <= ran_lw->range[1]);
+
+ /* Convert from nanometer to meter */
+ wlens[0] *= 1.e-9;
+ wlens[1] *= 1.e-9;
+
+ /* Compute the probability of the current band */
+ pdf[i] = blackbody_fraction(wlens[0], wlens[1], ran_lw->ref_temperature);
+
+ /* Update the norm */
+ sum += pdf[i];
+ }
+ ASSERT(eq_eps(sum, blackbody_fraction(ran_lw->range[0]*1.e-9,
+ ran_lw->range[1]*1.e-9, ran_lw->ref_temperature), 1.e-6f));
+
+ /* Normalize the probabilities */
+ FOR_EACH(i, 0, nbands) pdf[i] /= sum;
+ }
+
+exit:
+ return res;
+error:
+ darray_double_clear(&ran_lw->pdf_sky_bands);
+ goto exit;
+}
+
+
static double
ran_lw_sample_continue
(const struct htrdr_ran_lw* ran_lw,
@@ -258,6 +322,7 @@ release_ran_lw(ref_T* ref)
ran_lw = CONTAINER_OF(ref, struct htrdr_ran_lw, ref);
darray_double_release(&ran_lw->cdf);
darray_double_release(&ran_lw->pdf);
+ darray_double_release(&ran_lw->pdf_sky_bands);
MEM_RM(ran_lw->htrdr->allocator, ran_lw);
}
@@ -267,19 +332,16 @@ release_ran_lw(ref_T* ref)
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 sky_range[2];
double min_band_len = 0;
res_T res = RES_OK;
- ASSERT(htrdr && range && out_ran_lw && ref_temperature > 0);
+ ASSERT(htrdr && 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) {
@@ -293,9 +355,12 @@ htrdr_ran_lw_create
ran_lw->htrdr = htrdr;
darray_double_init(htrdr->allocator, &ran_lw->cdf);
darray_double_init(htrdr->allocator, &ran_lw->pdf);
+ darray_double_init(htrdr->allocator, &ran_lw->pdf_sky_bands);
+
+ HTSKY(get_spectral_bounds(htrdr->sky, sky_range));
- ran_lw->range[0] = range[0];
- ran_lw->range[1] = range[1];
+ ran_lw->range[0] = sky_range[0];
+ ran_lw->range[1] = sky_range[1];
ran_lw->ref_temperature = ref_temperature;
ran_lw->nbands = nbands;
@@ -304,21 +369,24 @@ htrdr_ran_lw_create
if(nbands == HTRDR_RAN_LW_CONTINUE) {
ran_lw->band_len = 0;
} else {
- ran_lw->band_len = (range[1] - range[0]) / (double)ran_lw->nbands;
+ ran_lw->band_len = (sky_range[1] - sky_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);
+ ran_lw->nbands = (size_t)ceil((sky_range[1] - sky_range[0]) / ran_lw->band_len);
}
res = setup_ran_lw_cdf(ran_lw, FUNC_NAME);
if(res != RES_OK) goto error;
}
+ res = compute_sky_bands_pdf(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]);
+ sky_range[0], sky_range[1]);
exit:
*out_ran_lw = ran_lw;
@@ -353,6 +421,7 @@ htrdr_ran_lw_sample
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)) {
+ FATAL("Unexpected behaviour\n");
if(pdf) *pdf = 1;
return ran_lw->range[0];
} else { /* Continue */
@@ -361,3 +430,16 @@ htrdr_ran_lw_sample
}
}
+double
+htrdr_ran_lw_get_sky_band_pdf
+ (const struct htrdr_ran_lw* ran_lw,
+ const size_t iband)
+{
+ size_t i, n;
+ ASSERT(ran_lw);
+
+ n = htsky_get_spectral_bands_count(ran_lw->htrdr->sky);
+ i = iband - htsky_get_spectral_band_id(ran_lw->htrdr->sky, 0);
+ ASSERT(i < n);
+ return darray_double_cdata_get(&ran_lw->pdf_sky_bands)[i];
+}
diff --git a/src/htrdr_ran_lw.h b/src/htrdr_ran_lw.h
@@ -26,10 +26,9 @@ 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 th CDF */
+ const size_t nbands, /* Hint on #bands used to discretised the CDF */
const double ref_temperature, /* Reference temperature */
struct htrdr_ran_lw** ran_lw);
@@ -49,5 +48,10 @@ htrdr_ran_lw_sample
const double r1, /* Canonical number in [0, 1[ */
double* pdf); /* May be NULL */
+extern LOCAL_SYM double
+htrdr_ran_lw_get_sky_band_pdf
+ (const struct htrdr_ran_lw* ran_lw,
+ const size_t iband);
+
#endif /* HTRDR_RAN_LW_H */
