commit 0a4c1facc1a8dd69f1c66892f1706fca85920f09
parent e052aa12853fb92b01ab04a99b5c4bec3334e9d9
Author: Najda Villefranque <najda.villefranque@gmail.com>
Date: Mon, 18 May 2020 14:18:04 +0200
Fix the spectral sampling in the longwave
Diffstat:
6 files changed, 129 insertions(+), 19 deletions(-)
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -714,7 +714,7 @@ brightness_temperature
(struct htrdr* htrdr,
const double lambda_min,
const double lambda_max,
- const double radiance,
+ const double radiance, /* In W/m2/sr/m */
double* temperature)
{
const size_t MAX_ITER = 100;
@@ -761,8 +761,8 @@ exit:
return res;
error:
htrdr_log_err(htrdr,
- "Could not compute the brightness temperature for the radiance %g "
- "estimated over [%g, %g] nanometers.\n",
+ "Could not compute the brightness temperature for the estimated radiance %g "
+ "averaged over [%g, %g] nanometers.\n",
radiance,
lambda_min*1e9,
lambda_max*1e9);
diff --git a/src/htrdr_c.h b/src/htrdr_c.h
@@ -197,7 +197,7 @@ brightness_temperature
(struct htrdr* htrdr,
const double lambda_min, /* In meters */
const double lambda_max, /* In meters */
- /* Integrated over [lambda_min, lambda_max]. In W/m^2/sr/m */
+ /* Averaged over [lambda_min, lambda_max]. In W/m^2/sr/m */
const double radiance,
double* temperature);
diff --git a/src/htrdr_compute_radiance_lw.c b/src/htrdr_compute_radiance_lw.c
@@ -19,6 +19,7 @@
#include "htrdr_interface.h"
#include "htrdr_ground.h"
#include "htrdr_solve.h"
+#include "htrdr_ran_lw.h"
#include <high_tune/htsky.h>
@@ -152,21 +153,23 @@ htrdr_compute_radiance_lw
double range[2];
double pos_next[3];
double dir_next[3];
- double band_bounds[2]; /* In nanometers */
- double band_bounds_m[2]; /* In meters */
+ double wlen_bounds[2]; /* sub-interval bounds in nano meters */
+ double wlen_bounds_m[2]; /* sub-interval bounds in meters */
+ double delta_wlen ; /* sub-interval size in meters */
double temperature;
double g;
double w = 0; /* Weight */
ASSERT(htrdr && rng && pos_in && dir_in && ithread < htrdr->nthreads);
- /* Retrieve the band boundaries */
- htsky_get_spectral_band_bounds(htrdr->sky, iband, band_bounds);
+ /* Retrieve the band boundaries in nano meters*/
+ htrdr_ran_lw_get_wlen_band_bounds(htrdr->ran_lw, wlen, wlen_bounds);
+ /* Convert to meters */
+ wlen_bounds_m[0] = wlen_bounds[0] * 1e-9;
+ wlen_bounds_m[1] = wlen_bounds[1] * 1e-9;
- /* Transform the band boundaries in meters and clamp them to the integration
- * domain */
- band_bounds_m[0] = MMAX(band_bounds[0] * 1e-9, htrdr->wlen_range_m[0]);
- band_bounds_m[1] = MMIN(band_bounds[1] * 1e-9, htrdr->wlen_range_m[1]);
+ /* If monochromatic, set delta_wlen to 1 m, otherwise lmax - lmin in meters */
+ delta_wlen = eq_eps(wlen_bounds[0],wlen_bounds[1],1e-6) ? 1 : wlen_bounds_m[1]-wlen_bounds_m[0] ;
/* Setup the phase function for this spectral band & quadrature point */
CHK(RES_OK == ssf_phase_create
@@ -220,7 +223,8 @@ htrdr_compute_radiance_lw
if(ctx.event_type == EVENT_ABSORPTION) {
ASSERT(!SVX_HIT_NONE(&svx_hit));
temperature = htsky_fetch_temperature(htrdr->sky, pos_next);
- w = planck(band_bounds_m[0], band_bounds_m[1], temperature);
+ /* weight is planck integrated over the spectral sub-interval */
+ w = planck(wlen_bounds_m[0],wlen_bounds_m[1],temperature)*delta_wlen ;
break;
}
@@ -268,7 +272,8 @@ htrdr_compute_radiance_lw
if(temperature <= 0) {
w = 0;
} else {
- w = planck(band_bounds_m[0], band_bounds_m[1], temperature);
+ /* weight is planck integrated over the spectral sub-interval */
+ w = planck(wlen_bounds_m[0],wlen_bounds_m[1],temperature) * delta_wlen ;
}
break;
}
@@ -279,7 +284,6 @@ htrdr_compute_radiance_lw
d3_set(dir, dir_next);
}
SSF(phase_ref_put(phase_hg));
-
return w;
}
diff --git a/src/htrdr_draw_radiance.c b/src/htrdr_draw_radiance.c
@@ -601,6 +601,7 @@ draw_pixel_lw
double ray_dir[3];
double weight;
double r0, r1;
+ double wlen_range_nm[2] ;
double wlen;
size_t iband;
size_t iquad;
@@ -629,12 +630,21 @@ draw_pixel_lw
iquad = htsky_spectral_band_sample_quadrature(htrdr->sky, r1, iband);
/* Retrieve the PDF to sample this sky band */
- band_pdf = htrdr_ran_lw_get_sky_band_pdf(htrdr->ran_lw, iband);
+ /* band_pdf is the pdf in the band of ran_lw->band_len nm */
+ band_pdf = htrdr_ran_lw_get_wlen_band_pdf(htrdr->ran_lw,wlen) ;
- /* Compute the 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 */
weight /= band_pdf;
+
+ htrdr_ran_lw_get_wlen_range(htrdr->ran_lw, wlen_range_nm) ;
+ /* From integrated radiance to average radiance in W/m^2/sr/m */
+ if (!eq_eps(wlen_range_nm[0], wlen_range_nm[1], 1e-6)) {
+ 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;
struct darray_double cdf;
/* Probas to sample a sky band overlapped by the range */
struct darray_double sky_bands_pdf;
@@ -64,9 +65,16 @@ setup_ran_lw_cdf
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 long wave spectral bands -- %s.\n",
+ func_name, res_to_cstr(res));
+ goto error;
+ }
cdf = darray_double_data_get(&ran_lw->cdf);
- pdf = cdf; /* Alias the CDF by the PDF since only one array is necessary */
+ pdf = darray_double_data_get(&ran_lw->pdf); /* Now save pdf to correct weight */
htrdr_log(ran_lw->htrdr,
"Number of bands of the long wave cumulative: %lu\n",
@@ -85,7 +93,7 @@ setup_ran_lw_cdf
lambda_hi *= 1.e-9;
/* Compute the probability of the current band */
- pdf[i] = planck(lambda_lo, lambda_hi, ran_lw->ref_temperature);
+ pdf[i] = blackbody_fraction(lambda_lo, lambda_hi, ran_lw->ref_temperature);
/* Update the norm */
sum += pdf[i];
@@ -111,6 +119,7 @@ exit:
return res;
error:
darray_double_clear(&ran_lw->cdf);
+ darray_double_clear(&ran_lw->pdf);
goto exit;
}
@@ -291,6 +300,7 @@ release_ran_lw(ref_T* ref)
ASSERT(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->sky_bands_pdf);
MEM_RM(ran_lw->htrdr->allocator, ran_lw);
}
@@ -326,6 +336,7 @@ htrdr_ran_lw_create
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);
darray_double_init(htrdr->allocator, &ran_lw->sky_bands_pdf);
ran_lw->range[0] = range[0];
@@ -392,6 +403,70 @@ htrdr_ran_lw_sample
}
}
+res_T
+htrdr_ran_lw_get_wlen_range
+ (const struct htrdr_ran_lw* ran_lw,
+ double * wlens)
+{
+ wlens[0] = ran_lw->range[0] ;
+ wlens[1] = ran_lw->range[1] ;
+ return RES_OK ;
+}
+
+res_T
+htrdr_ran_lw_get_wlen_band_bounds
+ (const struct htrdr_ran_lw* ran_lw,
+ const double wlen,
+ double * wlens)
+{
+ size_t i = htrdr_ran_lw_get_index_wlen(ran_lw,wlen) ;
+ wlens[0] = ran_lw->range[0] + (double) i*ran_lw->band_len ;
+ wlens[1] = ran_lw->range[0] + (double)(i+1)*ran_lw->band_len ;
+ return RES_OK ;
+}
+
+size_t
+htrdr_ran_lw_get_index_wlen
+ (const struct htrdr_ran_lw* ran_lw,
+ const double wlen)
+{
+ /* i==0 corresponds to pdf for band
+ * [range[0], range[0] + band_len]
+ * i==j corresponds to pdf for band
+ * [range[0] + j*band_len , range[0] +(j+1)* band_len]
+ * to find i, floor((wlen - range[0]) / band_len) */
+ if (eq_eps(ran_lw->range[0], ran_lw->range[1], 1.e-6)) {
+ return 0 ;
+ }
+ else {
+ return (size_t) ((wlen - ran_lw->range[0]) / ran_lw->band_len );
+ }
+}
+
+double
+htrdr_ran_lw_get_wlen_band_pdf
+ (const struct htrdr_ran_lw* ran_lw,
+ const double wlen)
+{
+ /* need to retrieve ran_lw->pdf[i]
+ * where i corresponds to the wlen */
+ size_t i ;
+ ASSERT(ran_lw);
+ ASSERT(wlen >= ran_lw->range[0]);
+ ASSERT(wlen <= ran_lw->range[1]);
+
+ i = htrdr_ran_lw_get_index_wlen(ran_lw,wlen) ;
+
+ if (darray_double_size_get(&ran_lw->pdf)==0) {
+ return 1 ;
+ }
+ else {
+ /* Fetch its PDF */
+ ASSERT(i < darray_double_size_get(&ran_lw->pdf));
+ return darray_double_cdata_get(&ran_lw->pdf)[i];
+ }
+}
+
double
htrdr_ran_lw_get_sky_band_pdf
(const struct htrdr_ran_lw* ran_lw,
diff --git a/src/htrdr_ran_lw.h b/src/htrdr_ran_lw.h
@@ -47,6 +47,27 @@ htrdr_ran_lw_sample
(const struct htrdr_ran_lw* ran_lw,
const double r); /* Canonical number in [0, 1[ */
+size_t
+htrdr_ran_lw_get_index_wlen
+ (const struct htrdr_ran_lw* ran_lw,
+ const double wlen);
+
+res_T
+htrdr_ran_lw_get_wlen_range
+ (const struct htrdr_ran_lw* ran_lw,
+ double * wlens) ;
+
+res_T
+htrdr_ran_lw_get_wlen_band_bounds
+ (const struct htrdr_ran_lw* ran_lw,
+ const double wlen,
+ double * wlens) ;
+
+double
+htrdr_ran_lw_get_wlen_band_pdf
+ (const struct htrdr_ran_lw* ran_lw,
+ const double wlen);
+
extern LOCAL_SYM double
htrdr_ran_lw_get_sky_band_pdf
(const struct htrdr_ran_lw* ran_lw,
