commit f187fa13f89458d34d7675d91250c9c04daf6bcf
parent 374d459992b2797e70ff1d1a58c1c8c067fad93f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 30 Jul 2018 11:21:05 +0200
Add the quadrature CDF
Diffstat:
5 files changed, 77 insertions(+), 27 deletions(-)
diff --git a/src/htgop.c b/src/htgop.c
@@ -274,7 +274,8 @@ parse_spectral_intervals
{
unsigned long ispecint;
double* wave_numbers = NULL;
- struct darray_double* quadratures = NULL;
+ struct darray_double* pdfs = NULL;
+ struct darray_double* cdfs = NULL;
res_T res = RES_OK;
(void)htgop;
ASSERT(htgop && rdr && nspecints && specints);
@@ -283,14 +284,18 @@ parse_spectral_intervals
/* Allocate the memory space for the wave numbers and the quadratures */
CALL(darray_double_resize(&specints->wave_numbers, nspecints + 1));
- CALL(darray_dbllst_resize(&specints->quadratures, nspecints));
+ CALL(darray_dbllst_resize(&specints->quadrature_pdfs, nspecints));
+ CALL(darray_dbllst_resize(&specints->quadrature_cdfs, nspecints));
wave_numbers = darray_double_data_get(&specints->wave_numbers);
- quadratures = darray_dbllst_data_get(&specints->quadratures);
+ pdfs = darray_dbllst_data_get(&specints->quadrature_pdfs);
+ cdfs = darray_dbllst_data_get(&specints->quadrature_cdfs);
FOR_EACH(ispecint, 0, nspecints) {
- double* quadrature = NULL;
+ double* pdf = NULL;
+ double* cdf = NULL;
double wave_number_low;
double wave_number_upp;
+ double sum = 0;
unsigned long quad_len;
unsigned long iquad;
@@ -315,13 +320,21 @@ parse_spectral_intervals
CALL(cstr_to_ulong(read_line(rdr), &quad_len));
/* Allocate the points of the quadrature */
- CALL(darray_double_resize(&quadratures[ispecint], quad_len));
- quadrature = darray_double_data_get(&quadratures[ispecint]);
+ CALL(darray_double_resize(&pdfs[ispecint], quad_len));
+ CALL(darray_double_resize(&cdfs[ispecint], quad_len));
+ pdf = darray_double_data_get(&pdfs[ispecint]);
+ cdf = darray_double_data_get(&cdfs[ispecint]);
/* Read the weight of the quadrature points */
+ sum = 0;
FOR_EACH(iquad, 0, quad_len) {
- CALL(cstr_to_double(read_line(rdr), &quadrature[iquad]));
+ CALL(cstr_to_double(read_line(rdr), &pdf[iquad]));
+ cdf[iquad] = iquad == 0 ? pdf[iquad] : pdf[iquad]+cdf[iquad-1];
+ sum += pdf[iquad];
}
+ ASSERT(eq_eps(sum, 1.0, 1.e6));
+ ASSERT(eq_eps(cdf[quad_len-1], 1.0, 1.e6));
+ cdf[quad_len-1] = 1.0; /* Handle numerical imprecision */
}
#undef CALL
@@ -330,7 +343,8 @@ exit:
return res;
error:
darray_double_clear(&specints->wave_numbers);
- darray_dbllst_clear(&specints->quadratures);
+ darray_dbllst_clear(&specints->quadrature_pdfs);
+ darray_dbllst_clear(&specints->quadrature_cdfs);
goto exit;
}
@@ -744,7 +758,7 @@ res_T
htgop_get_lw_spectral_intervals_count(const struct htgop* htgop, size_t* nspecints)
{
if(!htgop || !nspecints) return RES_BAD_ARG;
- *nspecints = darray_dbllst_size_get(&htgop->lw_specints.quadratures);
+ *nspecints = darray_dbllst_size_get(&htgop->lw_specints.quadrature_pdfs);
return RES_OK;
}
@@ -752,7 +766,7 @@ res_T
htgop_get_sw_spectral_intervals_count(const struct htgop* htgop, size_t* nspecints)
{
if(!htgop || !nspecints) return RES_BAD_ARG;
- *nspecints = darray_dbllst_size_get(&htgop->sw_specints.quadratures);
+ *nspecints = darray_dbllst_size_get(&htgop->sw_specints.quadrature_pdfs);
return RES_OK;
}
@@ -781,7 +795,8 @@ htgop_get_lw_spectral_interval
struct htgop_spectral_interval* interval)
{
const double* wave_numbers;
- const struct darray_double* quadratures;
+ const struct darray_double* quad_cdfs;
+ const struct darray_double* quad_pdfs;
size_t n;
if(!htgop || !interval) return RES_BAD_ARG;
HTGOP(get_lw_spectral_intervals_count(htgop, &n));
@@ -791,11 +806,13 @@ htgop_get_lw_spectral_interval
return RES_BAD_ARG;
}
wave_numbers = darray_double_cdata_get(&htgop->lw_specints.wave_numbers);
- quadratures = darray_dbllst_cdata_get(&htgop->lw_specints.quadratures);
+ quad_pdfs = darray_dbllst_cdata_get(&htgop->lw_specints.quadrature_pdfs);
+ quad_cdfs = darray_dbllst_cdata_get(&htgop->lw_specints.quadrature_cdfs);
interval->wave_numbers[0] = wave_numbers[ispecint+0];
interval->wave_numbers[1] = wave_numbers[ispecint+1];
- interval->quadrature = darray_double_cdata_get(&quadratures[ispecint]);
- interval->quadrature_length = darray_double_size_get(&quadratures[ispecint]);
+ interval->quadrature_pdf = darray_double_cdata_get(&quad_pdfs[ispecint]);
+ interval->quadrature_cdf = darray_double_cdata_get(&quad_cdfs[ispecint]);
+ interval->quadrature_length = darray_double_size_get(&quad_pdfs[ispecint]);
return RES_OK;
}
@@ -806,7 +823,8 @@ htgop_get_sw_spectral_interval
struct htgop_spectral_interval* interval)
{
const double* wave_numbers;
- const struct darray_double* quadratures;
+ const struct darray_double* quad_pdfs;
+ const struct darray_double* quad_cdfs;
size_t n;
if(!htgop || !interval) return RES_BAD_ARG;
HTGOP(get_sw_spectral_intervals_count(htgop, &n));
@@ -816,11 +834,13 @@ htgop_get_sw_spectral_interval
return RES_BAD_ARG;
}
wave_numbers = darray_double_cdata_get(&htgop->sw_specints.wave_numbers);
- quadratures = darray_dbllst_cdata_get(&htgop->sw_specints.quadratures);
+ quad_pdfs = darray_dbllst_cdata_get(&htgop->sw_specints.quadrature_pdfs);
+ quad_cdfs = darray_dbllst_cdata_get(&htgop->sw_specints.quadrature_cdfs);
interval->wave_numbers[0] = wave_numbers[ispecint+0];
interval->wave_numbers[1] = wave_numbers[ispecint+1];
- interval->quadrature = darray_double_cdata_get(&quadratures[ispecint]);
- interval->quadrature_length = darray_double_size_get(&quadratures[ispecint]);
+ interval->quadrature_pdf = darray_double_cdata_get(&quad_pdfs[ispecint]);
+ interval->quadrature_cdf = darray_double_cdata_get(&quad_cdfs[ispecint]);
+ interval->quadrature_length = darray_double_size_get(&quad_pdfs[ispecint]);
return RES_OK;
}
diff --git a/src/htgop.h b/src/htgop.h
@@ -51,7 +51,8 @@ struct htgop_ground {
struct htgop_spectral_interval {
double wave_numbers[2]; /* Lower/Upper wave number of the interval in cm^-1 */
- const double* quadrature; /* Weights of the quadrature points */
+ const double* quadrature_pdf; /* PDF of the quadrature */
+ const double* quadrature_cdf; /* CDF of the quadrature */
size_t quadrature_length; /* #quadrature points of the interval */
};
@@ -202,6 +203,13 @@ htgop_get_sw_spectral_interval
struct htgop_spectral_interval* interval);
HTGOP_API res_T
+htgop_sw_spectral_interval_sample_quadrature
+ (const struct htgop* htgop,
+ const size_t ispectral_interval,
+ const double r, /* Canonical random number in [0, 1[ */
+ size_t* iquad); /* Id of the sampled interval */
+
+HTGOP_API res_T
htgop_layer_get_lw_spectral_interval
(const struct htgop_layer* layer,
const size_t ispectral_interval,
diff --git a/src/htgop_parse_layers_spectral_intervals_data.h b/src/htgop_parse_layers_spectral_intervals_data.h
@@ -54,7 +54,7 @@ CONCAT(parse_layers_spectral_intervals_, DATA)
struct darray_double* XDATA(tab);
quad = darray_dbllst_data_get
- (&htgop->XDOMAIN(specints).quadratures) + ispecint;
+ (&htgop->XDOMAIN(specints).quadrature_pdfs) + ispecint;
quad_len = darray_double_size_get(quad);
/* Allocate the per layer data */
diff --git a/src/htgop_spectral_intervals.h b/src/htgop_spectral_intervals.h
@@ -22,7 +22,8 @@ struct spectral_intervals {
/* List of wave numbers, in cm^-1, sorted in ascending order.
* #wave numbers == #intervals + 1 */
struct darray_double wave_numbers;
- struct darray_dbllst quadratures; /* Per spectral interval quadrature weight */
+ struct darray_dbllst quadrature_pdfs; /* Per spectral band quadrature pdf */
+ struct darray_dbllst quadrature_cdfs; /* Per spectral band quadrature cdf */
};
static INLINE void
@@ -31,7 +32,9 @@ spectral_intervals_init
{
ASSERT(sinters);
darray_double_init(allocator, &sinters->wave_numbers);
- darray_dbllst_init(allocator, &sinters->quadratures);
+ darray_dbllst_init(allocator, &sinters->quadrature_pdfs);
+ darray_dbllst_init(allocator, &sinters->quadrature_cdfs);
+
}
static INLINE void
@@ -39,7 +42,8 @@ spectral_intervals_release(struct spectral_intervals* sinters)
{
ASSERT(sinters);
darray_double_release(&sinters->wave_numbers);
- darray_dbllst_release(&sinters->quadratures);
+ darray_dbllst_release(&sinters->quadrature_pdfs);
+ darray_dbllst_release(&sinters->quadrature_cdfs);
}
static INLINE res_T
@@ -50,7 +54,8 @@ spectral_intervals_copy
ASSERT(dst && src);
#define CALL(Func) { if(RES_OK != (res = Func)) return res; } (void)0
CALL(darray_double_copy(&dst->wave_numbers, &src->wave_numbers));
- CALL(darray_dbllst_copy(&dst->quadratures, &src->quadratures));
+ CALL(darray_dbllst_copy(&dst->quadrature_pdfs, &src->quadrature_pdfs));
+ CALL(darray_dbllst_copy(&dst->quadrature_cdfs, &src->quadrature_cdfs));
#undef CALL
return RES_OK;
}
@@ -63,7 +68,8 @@ spectral_intervals_copy_and_release
ASSERT(dst && src);
#define CALL(Func) { if(RES_OK != (res = Func)) return res; } (void)0
CALL(darray_double_copy_and_release(&dst->wave_numbers, &src->wave_numbers));
- CALL(darray_dbllst_copy_and_release(&dst->quadratures, &src->quadratures));
+ CALL(darray_dbllst_copy_and_release(&dst->quadrature_pdfs, &src->quadrature_pdfs));
+ CALL(darray_dbllst_copy_and_release(&dst->quadrature_pdfs, &src->quadrature_cdfs));
#undef CALL
spectral_intervals_release(src);
return RES_OK;
diff --git a/src/test_htgop_load.c b/src/test_htgop_load.c
@@ -182,8 +182,16 @@ main(int argc, char** argv)
FOR_EACH(iquad, 0, specint.quadrature_length) {
CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK);
- CHK(specint.quadrature[iquad] == dbl);
+ CHK(specint.quadrature_pdf[iquad] == dbl);
+ if(iquad == 0) {
+ CHK(specint.quadrature_cdf[iquad] == dbl);
+ } else {
+ const double pdf =
+ specint.quadrature_cdf[iquad] - specint.quadrature_cdf[iquad-1];
+ CHK(eq_eps(pdf, dbl, 1.e-6));
+ }
}
+ CHK(specint.quadrature_cdf[specint.quadrature_length-1] == 1.0);
}
/* Per SW spectral interval data */
@@ -212,7 +220,15 @@ main(int argc, char** argv)
FOR_EACH(iquad, 0, specint.quadrature_length) {
CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK);
- CHK(specint.quadrature[iquad] == dbl);
+ CHK(specint.quadrature_pdf[iquad] == dbl);
+ if(iquad == 0) {
+ CHK(specint.quadrature_cdf[iquad] == dbl);
+ } else {
+ const double pdf =
+ specint.quadrature_cdf[iquad] - specint.quadrature_cdf[iquad-1];
+ CHK(eq_eps(pdf, dbl, 1.e-6));
+ }
+ CHK(specint.quadrature_cdf[specint.quadrature_length-1] == 1.0);
}
}
