commit ed6cefe7e8b7a97565aba1c8278ea8430b73093b
parent f696b9bd29b4d34ca4630b30d6309425cc668b03
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 29 May 2020 17:01:17 +0200
Revert "In longwave strictly use the sky spectral bands data"
This reverts commit f696b9bd29b4d34ca4630b30d6309425cc668b03.
Diffstat:
5 files changed, 29 insertions(+), 126 deletions(-)
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -499,16 +499,15 @@ htrdr_init
} else { /* Long Wave random variate */
const double Tref = 290; /* In Kelvin */
- double sky_range[2];
size_t n;
- 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 */
+ htrdr->wlen_range_m[0] = args->wlen_lw_range[0]*1e-9; /* Convert in meters */
+ htrdr->wlen_range_m[1] = args->wlen_lw_range[1]*1e-9; /* Convert in meters */
ASSERT(htrdr->wlen_range_m[0] <= htrdr->wlen_range_m[1]);
n = (size_t)(args->wlen_lw_range[1] - args->wlen_lw_range[0]);
- res = htrdr_ran_lw_create(htrdr, n, Tref, &htrdr->ran_lw);
+ res = htrdr_ran_lw_create
+ (htrdr, args->wlen_lw_range, 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,9 +154,7 @@ htrdr_compute_radiance_lw
double pos_next[3];
double dir_next[3];
double temperature;
- double sky_band_bounds[2];
- double sky_band_bounds_m[2];
- double sky_band_len_m;
+ double wlen_m = wlen * 1.e-9;
double g;
double w = 0; /* Weight */
@@ -165,16 +163,10 @@ 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_particle_phase_function_asymmetry_parameter
- (htrdr->sky, iband, iquad);
+ g = htsky_fetch_per_wavelength_particle_phase_function_asymmetry_parameter
+ (htrdr->sky, wlen);
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);
@@ -221,8 +213,7 @@ 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(sky_band_bounds_m[0], sky_band_bounds_m[1], temperature);
- w = w * sky_band_len_m;
+ w = planck_monochromatic(wlen_m, temperature);
break;
}
@@ -271,8 +262,7 @@ htrdr_compute_radiance_lw
w = 0;
} else {
/* weight is planck integrated over the spectral sub-interval */
- w = planck(sky_band_bounds_m[0], sky_band_bounds_m[1], temperature);
- w = w * sky_band_len_m;
+ w = planck_monochromatic(wlen_m, temperature);
}
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, NULL);
+ wlen = htrdr_ran_lw_sample(htrdr->ran_lw, 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 */
+ /* Compute the integrated luminance in W/m^2/sr */
weight = htrdr_compute_radiance_lw(htrdr, ithread, rng, ray_org, ray_dir,
wlen, iband, iquad);
- /* 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);
+ /* Importance sampling: correct weight with pdf */
weight /= band_pdf;
- /* 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]) ;
+ /* 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 */
diff --git a/src/htrdr_ran_lw.c b/src/htrdr_ran_lw.c
@@ -30,7 +30,6 @@
#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 */
@@ -122,69 +121,6 @@ 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,
@@ -322,7 +258,6 @@ 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);
}
@@ -332,16 +267,19 @@ 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 && out_ran_lw && ref_temperature > 0);
+ 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) {
@@ -355,12 +293,9 @@ 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] = sky_range[0];
- ran_lw->range[1] = sky_range[1];
+ ran_lw->range[0] = range[0];
+ ran_lw->range[1] = range[1];
ran_lw->ref_temperature = ref_temperature;
ran_lw->nbands = nbands;
@@ -369,24 +304,21 @@ htrdr_ran_lw_create
if(nbands == HTRDR_RAN_LW_CONTINUE) {
ran_lw->band_len = 0;
} else {
- ran_lw->band_len = (sky_range[1] - sky_range[0]) / (double)ran_lw->nbands;
+ 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((sky_range[1] - sky_range[0]) / ran_lw->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;
}
- 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",
- sky_range[0], sky_range[1]);
+ range[0], range[1]);
exit:
*out_ran_lw = ran_lw;
@@ -421,7 +353,6 @@ 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 */
@@ -430,16 +361,3 @@ 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,9 +26,10 @@ 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 size_t nbands, /* Hint on #bands used to discretised th CDF */
const double ref_temperature, /* Reference temperature */
struct htrdr_ran_lw** ran_lw);
@@ -48,10 +49,5 @@ 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 */
