commit fb44c28d0176fa571d5c96bd754e4c27903c45a0
parent 69d3ddf421892235c04c476ac25f4dc53fe5e6f4
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 31 Jul 2018 11:19:07 +0200
Refactoring and implement htgop_fetch_<sw|lw>_<ka|ks> functions
Diffstat:
9 files changed, 526 insertions(+), 371 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -41,9 +41,12 @@ set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(HTGOP_FILES_SRC
- htgop.c)
+ htgop.c
+ htgop_fetch_radiative_properties.c
+ htgop_sample_sw_spectral_interval_CIE_1931_XYZ.c)
set(SSOL_FILES_INC
htgop_c.h
+ htgop_fetch_radiative_properties.h
htgop_layer.h
htgop_parse_layers_spectral_intervals_data.h
htgop_reader.h
diff --git a/src/htgop.c b/src/htgop.c
@@ -26,14 +26,6 @@
#include <math.h>
-#define WNUM_TO_WLEN(Num) (1.0e7 / ((double)Num))
-#define WLEN_TO_WNUM(Len) WNUM_TO_WLEN(Len)
-
-#define CIE_XYZ_WAVELENGTH_MIN 380.0 /* In nanometer */
-#define CIE_XYZ_WAVELENGTH_MAX 780.0 /* In nanometer */
-#define CIE_XYZ_WAVENUMBER_MIN WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MAX) /*In cm^-1*/
-#define CIE_XYZ_WAVENUMBER_MAX WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MIN) /*In cm^-1*/
-
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -53,14 +45,6 @@ log_msg
}
static INLINE int
-cmp_dbl(const void* a, const void* b)
-{
- const double d0 = *((const double*)a);
- const double d1 = *((const double*)b);
- return d0 < d1 ? -1 : (d0 > d1 ? 1 : 0);
-}
-
-static INLINE int
cmp_lvl(const void* key, const void* data)
{
const double height = *((const double*)key);
@@ -68,215 +52,6 @@ cmp_lvl(const void* key, const void* data)
return height < lvl->height ? -1 : (height > lvl->height ? 1 : 0);
}
-static FINLINE double
-trapezoidal_integration
- (const double lambda_lo, /* Integral lower bound. In nanometer */
- const double lambda_hi, /* Integral upper bound. In nanometer */
- double (*f_bar)(const double lambda)) /* Function to integrate */
-{
- double dlambda;
- size_t i, n;
- double integral = 0;
- ASSERT(lambda_lo <= lambda_hi);
- ASSERT(lambda_lo > 0);
-
- n = (size_t)(lambda_hi - lambda_lo) + 1;
- dlambda = (lambda_hi - lambda_lo) / (double)n;
-
- FOR_EACH(i, 0, n) {
- const double lambda1 = lambda_lo + dlambda*(double)(i+0);
- const double lambda2 = lambda_hi + dlambda*(double)(i+1);
- const double f1 = f_bar(lambda1);
- const double f2 = f_bar(lambda2);
- integral += (f1 + f2)*dlambda*0.5;
- }
- return integral;
-}
-
-/* The following 3 functions are used to fit the CIE Xbar, Ybar and Zbar curved
- * has defined by the 1931 standard. These analytical fits are propsed by C.
- * Wyman, P. P. Sloan & P. Shirley in "Simple Analytic Approximations to the
- * CIE XYZ Color Matching Functions" - JCGT 2013. */
-static FINLINE double
-fit_x_bar_1931(const double lambda)
-{
- const double a = (lambda - 442.0) * (lambda < 442.0 ? 0.0624 : 0.0374);
- const double b = (lambda - 599.8) * (lambda < 599.8 ? 0.0264 : 0.0323);
- const double c = (lambda - 501.1) * (lambda < 501.1 ? 0.0490 : 0.0382);
- return 0.362*exp(-0.5*a*a) + 1.056*exp(-0.5f*b*b) - 0.065*exp(-0.5*c*c);
-}
-
-static FINLINE double
-fit_y_bar_1931(const double lambda)
-{
- const double a = (lambda - 568.8) * (lambda < 568.8 ? 0.0213 : 0.0247);
- const double b = (lambda - 530.9) * (lambda < 530.9 ? 0.0613 : 0.0322);
- return 0.821*exp(-0.5*a*a) + 0.286*exp(-0.5*b*b);
-}
-
-static FINLINE double
-fit_z_bar_1931(const double lambda)
-{
- const double a = (lambda - 437.0) * (lambda < 437.0 ? 0.0845 : 0.0278);
- const double b = (lambda - 459.0) * (lambda < 459.0 ? 0.0385 : 0.0725);
- return 1.217*exp(-0.5*a*a) + 0.681*exp(-0.5*b*b);
-}
-
-static res_T
-compute_CIE_XYZ_cdf(struct htgop* htgop)
-{
- enum { X, Y, Z }; /* Helper constant */
- size_t iband_min, iband_max;
- const double* wnums = NULL;
- double sum[3] = {0,0,0};
- double* pdf[3] = {NULL, NULL, NULL};
- double* cdf[3] = {NULL, NULL, NULL};
- size_t i, n;
- res_T res = RES_OK;
- ASSERT(htgop);
-
- HTGOP(get_sw_spectral_intervals_wave_numbers(htgop, &wnums));
- HTGOP(get_sw_spectral_intervals_count(htgop, &n));
- ASSERT(n);
-
- /* Allocate and reset the memory space for the tristimulus CDF */
- #define SETUP(Stimulus) { \
- res = darray_double_resize(&htgop->sw_ ## Stimulus ## _cdf, n); \
- if(res != RES_OK) { \
- log_err(htgop, \
- "Could not reserve the memory space for the CDF " \
- "of the "STR(X)" stimulus.\n"); \
- goto error; \
- } \
- memset(darray_double_data_get(&htgop->sw_ ## Stimulus ## _cdf), 0, \
- n*sizeof(double)); \
- cdf[Stimulus] = darray_double_data_get(&htgop->sw_ ## Stimulus ## _cdf); \
- pdf[Stimulus] = cdf[Stimulus]; \
- } (void)0
- SETUP(X);
- SETUP(Y);
- SETUP(Z);
- #undef RESERVE
-
- ASSERT(eq_eps(CIE_XYZ_WAVENUMBER_MAX, WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MIN), 1.e-6));
- ASSERT(eq_eps(CIE_XYZ_WAVENUMBER_MIN, WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MAX), 1.e-6));
- ASSERT(eq_eps(CIE_XYZ_WAVELENGTH_MAX, WNUM_TO_WLEN(CIE_XYZ_WAVENUMBER_MIN), 1.e-6));
- ASSERT(eq_eps(CIE_XYZ_WAVELENGTH_MIN, WNUM_TO_WLEN(CIE_XYZ_WAVENUMBER_MAX), 1.e-6));
-
- /* The sw spectral intervals does not intersect the CIE XYZ interval */
- if(wnums[0] > CIE_XYZ_WAVENUMBER_MAX || wnums[n] < CIE_XYZ_WAVENUMBER_MIN) {
- log_warn(htgop,
- "The loaded SW spectral intervals does not intersect "
- "the CIE XYZ interval.\n");
- goto exit;
- }
-
- /* Search the lower and upper bound of the bands that overlap the XYZ
- * spectral interval. Use a simple linear search since the overall number of
- * bands should be small */
- iband_min = 0, iband_max = n-1;
- FOR_EACH(i, 0, n) {
- if(wnums[i]<=CIE_XYZ_WAVENUMBER_MIN && CIE_XYZ_WAVENUMBER_MIN<=wnums[i+1])
- iband_min = i;
- if(wnums[i]<=CIE_XYZ_WAVENUMBER_MAX && CIE_XYZ_WAVENUMBER_MAX<=wnums[i+1])
- iband_max = i;
- }
- ASSERT(iband_min <= iband_max);
-
- /* Compute the *unormalized* pdf of the tristimulus */
- FOR_EACH(i, iband_min, iband_max+1) {
- const double lambda_lo = WNUM_TO_WLEN(wnums[i+1]);
- const double lambda_hi = WNUM_TO_WLEN(wnums[i+0]);
- ASSERT(lambda_lo <= lambda_hi);
- pdf[X][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_x_bar_1931);
- pdf[Y][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_y_bar_1931);
- pdf[Z][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_z_bar_1931);
- sum[X] += pdf[X][i];
- sum[Y] += pdf[Y][i];
- sum[Z] += pdf[Z][i];
- }
-
- /* Compute the cdf of the tristimulus. Note that the upper bound is 'n' and not
- * 'iband_max+1' in order to ensure that the CDF is strictly increasing. */
- FOR_EACH(i, iband_min, n) {
- /* Normalize the pdf */
- pdf[X][i] /= sum[X];
- pdf[Y][i] /= sum[Y];
- pdf[Z][i] /= sum[Z];
- /* Setup the cumulative */
- if(i == 0) {
- cdf[X][i] = pdf[X][i];
- cdf[Y][i] = pdf[Y][i];
- cdf[Z][i] = pdf[Z][i];
- } else {
- cdf[X][i] = pdf[X][i] + cdf[X][i-1];
- cdf[Y][i] = pdf[Y][i] + cdf[Y][i-1];
- cdf[Z][i] = pdf[Z][i] + cdf[Z][i-1];
- }
- }
- ASSERT(eq_eps(cdf[X][n-1], 1, 1.e-6));
- ASSERT(eq_eps(cdf[Y][n-1], 1, 1.e-6));
- ASSERT(eq_eps(cdf[Z][n-1], 1, 1.e-6));
-
-#ifndef NDEBUG
- FOR_EACH(i, 1, n) {
- ASSERT(cdf[X][i] >= cdf[X][i-1]);
- ASSERT(cdf[Y][i] >= cdf[Y][i-1]);
- ASSERT(cdf[Z][i] >= cdf[Z][i-1]);
- }
-#endif
-
- /* Handle numerical issue */
- cdf[X][n-1] = 1.0;
- cdf[Y][n-1] = 1.0;
- cdf[Z][n-1] = 1.0;
-
-exit:
- return res;
-error:
- darray_double_clear(&htgop->sw_X_cdf);
- darray_double_clear(&htgop->sw_Y_cdf);
- darray_double_clear(&htgop->sw_Z_cdf);
- goto exit;
-}
-
-static INLINE res_T
-sample_sw_spectral_interval_CIE
- (const struct htgop* htgop,
- const char* caller,
- const double* cdf,
- const size_t cdf_length,
- const double r, /* Canonical number in [0, 1[ */
- size_t* id)
-{
- double r_next = nextafter(r, DBL_MAX);
- double* find;
- size_t i;
- ASSERT(htgop && cdf && cdf_length);
-
- if(!id) return RES_BAD_ARG;
- if(r < 0 || r >= 1) {
- log_err(htgop, "%s: invalid canonical random number `%g'.\n", caller, r);
- return RES_BAD_ARG;
- }
-
- /* Use r_next rather than r in order to find the first entry that is not less
- * that *or equal* to r */
- find = search_lower_bound(&r_next, cdf, cdf_length, sizeof(double), cmp_dbl);
- if(!find) { /* <=> CIE XYZ band does not intersect the loaded intervals */
- log_err(htgop,
- "%s: could not sample a spectral interval in the CIE XYZ spectral band.\n",
- caller);
- return RES_BAD_OP;
- }
-
- i = (size_t)(find - cdf);
- ASSERT(i < cdf_length && cdf[i] > r && (!i || cdf[i-1] <= r));
- *id = i;
- return RES_OK;
-}
-
-
static res_T
parse_spectral_intervals
(struct htgop* htgop,
@@ -529,7 +304,7 @@ load_stream(struct htgop* htgop, FILE* stream, const char* stream_name)
#undef CALL
/* Compute the cdf for the SW CIE XYZ spectral interval */
- res = compute_CIE_XYZ_cdf(htgop);
+ res = compute_CIE_1931_XYZ_cdf(htgop);
if(res != RES_OK) goto error;
exit:
@@ -912,9 +687,8 @@ htgop_spectral_interval_sample_quadrature
return RES_BAD_ARG;
}
-
/* Use r_next rather than r in order to find the first entry that is not less
- * that *or equal* to r */
+ * than *or equal* to r */
find = search_lower_bound(&r_next, specint->quadrature_cdf,
specint->quadrature_length, sizeof(double), cmp_dbl);
ASSERT(find);
@@ -927,45 +701,6 @@ htgop_spectral_interval_sample_quadrature
}
res_T
-htgop_sample_sw_spectral_interval_CIE_X
- (const struct htgop* htgop,
- const double r, /* Canonical number in [0, 1[ */
- size_t* ispecint)
-{
- if(!htgop) return RES_BAD_ARG;
- return sample_sw_spectral_interval_CIE(htgop, FUNC_NAME,
- darray_double_cdata_get(&htgop->sw_X_cdf),
- darray_double_size_get(&htgop->sw_X_cdf),
- r, ispecint);
-}
-
-res_T
-htgop_sample_sw_spectral_interval_CIE_Y
- (const struct htgop* htgop,
- const double r, /* Canonical number in [0, 1[ */
- size_t* ispecint)
-{
- if(!htgop) return RES_BAD_ARG;
- return sample_sw_spectral_interval_CIE(htgop, FUNC_NAME,
- darray_double_cdata_get(&htgop->sw_Y_cdf),
- darray_double_size_get(&htgop->sw_Y_cdf),
- r, ispecint);
-}
-
-res_T
-htgop_sample_sw_spectral_interval_CIE_Z
- (const struct htgop* htgop,
- const double r, /* Canonical number in [0, 1[ */
- size_t* ispecint)
-{
- if(!htgop) return RES_BAD_ARG;
- return sample_sw_spectral_interval_CIE(htgop, FUNC_NAME,
- darray_double_cdata_get(&htgop->sw_Z_cdf),
- darray_double_size_get(&htgop->sw_Z_cdf),
- r, ispecint);
-}
-
-res_T
htgop_position_to_layer_id
(const struct htgop* htgop, const double pos, size_t* ilayer)
{
@@ -1001,62 +736,6 @@ htgop_position_to_layer_id
return RES_OK;
}
-res_T
-htgop_layer_fetch_sw_ka
- (const struct htgop_layer* layer,
- const size_t ispecint,
- const size_t iquad_point,
- const double x_h2o, /* H2O mol / mixture mol */
- double* out_ka)
-{
- struct htgop_layer_sw_spectral_interval specint;
- struct htgop_layer_sw_spectral_interval_tab tab;
- double* find;
- double x1, x2, k1, k2;
- double ka;
- res_T res = RES_OK;
-
- if(!layer || !out_ka) return RES_BAD_ARG;
-
- find = search_lower_bound(&x_h2o, layer->x_h2o_tab, layer->tab_length,
- sizeof(double), cmp_dbl);
- if(!find) return RES_BAD_ARG;
-
- res = htgop_layer_get_sw_spectral_interval(layer, ispecint, &specint);
- if(res != RES_OK) return res;
- res = htgop_layer_sw_spectral_interval_get_tab(&specint, iquad_point, &tab);
- if(res != RES_OK) return res;
-
- ASSERT(layer->tab_length);
- ASSERT(tab.tab_length == layer->tab_length);
-
- if(find == layer->x_h2o_tab) {
- ASSERT(layer->x_h2o_tab[0] >= x_h2o);
- x2 = layer->x_h2o_tab[0];
- k2 = tab.ka_tab[0];
- ka = k2 / x2*x_h2o;
- } else {
- const size_t i = (size_t)(find - layer->x_h2o_tab);
- ASSERT(i < layer->tab_length);
- ASSERT(layer->x_h2o_tab[i-0] >= x_h2o);
- ASSERT(layer->x_h2o_tab[i-1] < x_h2o);
- x1 = layer->x_h2o_tab[i-1];
- x2 = layer->x_h2o_tab[i-0];
- k1 = tab.ka_tab[i-1];
- k2 = tab.ka_tab[i-0];
-
- if(!k1 || !k2) {
- ka = k1 + (k2-k1)/(x2-x1) * (x_h2o-x1);
- } else {
- const double alpha = (log(k2) - log(k1))/(log(x1) - log(x2));
- const double beta = log(k1) - alpha * log(x1);
- ka = exp(alpha * log(x_h2o) + beta);
- }
- }
- *out_ka = ka;
- return RES_OK;
-}
-
/*******************************************************************************
* Local functions
******************************************************************************/
diff --git a/src/htgop.h b/src/htgop.h
@@ -51,7 +51,7 @@ struct htgop_ground {
struct htgop_spectral_interval {
double wave_numbers[2]; /* Lower/Upper wave number of the interval in cm^-1 */
- const double* quadrature_pdf; /* PDF of the quadrature */
+ const double* quadrature_pdf; /* Normalized PDF of the quadrature */
const double* quadrature_cdf; /* CDF of the quadrature */
size_t quadrature_length; /* #quadrature points of the interval */
const void* data__; /* Internal data */
@@ -80,7 +80,7 @@ struct htgop_layer_lw_spectral_interval {
struct htgop_layer_lw_spectral_interval_tab {
const double* ka_tab; /* Tabulated absorption coef */
- size_t tab_length;
+ size_t tab_length; /* == lengh of the tabulated xH2O */
};
struct htgop_layer_sw_spectral_interval {
@@ -93,13 +93,13 @@ struct htgop_layer_sw_spectral_interval {
struct htgop_layer_sw_spectral_interval_tab {
const double* ka_tab; /* Tabulated absorption coef */
const double* ks_tab; /* Tabulated scattering coef */
- size_t tab_length;
+ size_t tab_length; /* == lengh of the tabulated xH2O */
};
-BEGIN_DECLS
+BEGIN_DECLS /* HTGOP API */
/*******************************************************************************
- * HTGOP API
+ * HTGOP device management
******************************************************************************/
HTGOP_API res_T
htgop_create
@@ -116,6 +116,9 @@ HTGOP_API res_T
htgop_ref_put
(struct htgop* htgop);
+/*******************************************************************************
+ * Loading functions
+ ******************************************************************************/
HTGOP_API res_T
htgop_load
(struct htgop* htgop,
@@ -126,6 +129,9 @@ htgop_load_stream
(struct htgop* htgop,
FILE* stream);
+/*******************************************************************************
+ * Getters of the the loaded data
+ ******************************************************************************/
HTGOP_API res_T
htgop_get_ground
(const struct htgop* htgop,
@@ -164,24 +170,6 @@ htgop_get_sw_spectral_intervals_count
size_t* nspecints);
HTGOP_API res_T
-htgop_sample_sw_spectral_interval_CIE_X
- (const struct htgop* htgop,
- const double r, /* Canonical random number in [0, 1[ */
- size_t* ispecint); /* Id of the sampled interval */
-
-HTGOP_API res_T
-htgop_sample_sw_spectral_interval_CIE_Y
- (const struct htgop* htgop,
- const double r, /* Canonical random number in [0, 1[ */
- size_t* ispecint); /* Id of the sampled interval */
-
-HTGOP_API res_T
-htgop_sample_sw_spectral_interval_CIE_Z
- (const struct htgop* htgop,
- const double r, /* Canonical random number in [0, 1[ */
- size_t* ispecint); /* Id of the sampled interval */
-
-HTGOP_API res_T
htgop_get_lw_spectral_intervals_wave_numbers
(const struct htgop* htgop,
const double* wave_numbers[]);
@@ -204,12 +192,6 @@ htgop_get_sw_spectral_interval
struct htgop_spectral_interval* interval);
HTGOP_API res_T
-htgop_spectral_interval_sample_quadrature
- (const struct htgop_spectral_interval* interval,
- const double r, /* Canonical random number in [0, 1[ */
- size_t* iquad_point); /* Id of the sample quadrature point */
-
-HTGOP_API res_T
htgop_layer_get_lw_spectral_interval
(const struct htgop_layer* layer,
const size_t ispectral_interval,
@@ -233,13 +215,30 @@ htgop_layer_sw_spectral_interval_get_tab
const size_t iquadrature_point,
struct htgop_layer_sw_spectral_interval_tab* tab);
-HTGOP_API res_T
-htgop_position_to_layer_id
- (const struct htgop* htgop,
- const double pos,
- size_t* ilayer);
-
-/* TODO comment the interpolation procedure */
+/*******************************************************************************
+ * The following functions Retrieve the absorption/scattering coefficient in a
+ * layer, in short/long wave, at a given spectral interval for a specific
+ * quadrature point with respect to a submitted water vapor molar fraction
+ * x_h2o. x_h2o is compared to the values of tabulated water vapor molar
+ * fractions of the layer: x_h2o in [x1, x2], with the corresponding tabulated
+ * values of the required radiative property k(x1) = k1 and k(x2) = k2. Now
+ * there is 3 cases:
+ *
+ * - If x_h2o is smaller than the first tabulated value of the water vapor molar
+ * fraction, x1 and k1 is assumed to be null and k(x_h2o) is computed with a
+ * linear interpolation:
+ * k(x_h2o) = x_h2o * k2/x2
+ *
+ * - If either k1 or k2 is null, k(x_h2o) is also computed with a simple linear
+ * interpolation:
+ * k(x_h2o) = k1 + (k2-k1)/(x2-x1) * (x_h2o-x1)
+ *
+ * - Otherwise, k(x_h2o) is computed as bellow:
+ * k(x_h2o) = exp(alpha * ln(x_h2o) + beta)
+ * with:
+ * alpha = (ln(k2) - ln(k1)) / (ln(x2) - ln(x1))
+ * beta = ln(k1) - alpha * ln(x1)
+ ******************************************************************************/
HTGOP_API res_T
htgop_layer_fetch_lw_ka
(const struct htgop_layer* layer,
@@ -248,16 +247,14 @@ htgop_layer_fetch_lw_ka
const double x_h2o, /* H2O mol / mixture mol */
double* ka);
-/* TODO comment the interpolation procedure */
HTGOP_API res_T
htgop_layer_fetch_sw_ka
(const struct htgop_layer* layer,
const size_t ispecint,
const size_t iquad_point,
const double x_h2o, /* H2O mol / mixture mol */
- double* ka);
+ double* ka);
-/* TODO comment the interpolation procedure */
HTGOP_API res_T
htgop_layer_fetch_sw_ks
(const struct htgop_layer* layer,
@@ -266,6 +263,44 @@ htgop_layer_fetch_sw_ks
const double x_h2o, /* H2O mol / mixture mol */
double* ks);
+/*******************************************************************************
+ * Sample of a short wave spectral interval according to the CIE 1931 XYZ
+ * tristimulus values.
+ ******************************************************************************/
+HTGOP_API res_T
+htgop_sample_sw_spectral_interval_CIE_1931_X
+ (const struct htgop* htgop,
+ const double r, /* Canonical random number in [0, 1[ */
+ size_t* ispecint); /* Id of the sampled interval */
+
+HTGOP_API res_T
+htgop_sample_sw_spectral_interval_CIE_1931_Y
+ (const struct htgop* htgop,
+ const double r, /* Canonical random number in [0, 1[ */
+ size_t* ispecint); /* Id of the sampled interval */
+
+HTGOP_API res_T
+htgop_sample_sw_spectral_interval_CIE_1931_Z
+ (const struct htgop* htgop,
+ const double r, /* Canonical random number in [0, 1[ */
+ size_t* ispecint); /* Id of the sampled interval */
+
+/*******************************************************************************
+ * Miscellaneous functions
+ ******************************************************************************/
+/* Return the layer index into which the subimitted a 1D position lies */
+HTGOP_API res_T
+htgop_position_to_layer_id
+ (const struct htgop* htgop,
+ const double pos,
+ size_t* ilayer);
+
+HTGOP_API res_T
+htgop_spectral_interval_sample_quadrature
+ (const struct htgop_spectral_interval* interval,
+ const double r, /* Canonical random number in [0, 1[ */
+ size_t* iquad_point); /* Id of the sample quadrature point */
+
END_DECLS
#endif /* HTGOP */
diff --git a/src/htgop_c.h b/src/htgop_c.h
@@ -34,7 +34,7 @@ struct htgop {
struct spectral_intervals lw_specints;
struct spectral_intervals sw_specints;
- /* Cumulative distribution function for the Short Wave CIE XYZ tristimulus */
+ /* Cumulative distribution function for the Short Wave CIE XYZ values */
struct darray_double sw_X_cdf;
struct darray_double sw_Y_cdf;
struct darray_double sw_Z_cdf;
@@ -48,6 +48,14 @@ struct htgop {
ref_T ref;
};
+static INLINE int
+cmp_dbl(const void* a, const void* b)
+{
+ const double d0 = *((const double*)a);
+ const double d1 = *((const double*)b);
+ return d0 < d1 ? -1 : (d0 > d1 ? 1 : 0);
+}
+
extern LOCAL_SYM void
log_err
(const struct htgop* htgop,
@@ -68,5 +76,9 @@ log_warn
#endif
;
+extern LOCAL_SYM res_T
+compute_CIE_1931_XYZ_cdf
+ (struct htgop* htgop);
+
#endif /* HTGOP_C_H */
diff --git a/src/htgop_fetch_radiative_properties.c b/src/htgop_fetch_radiative_properties.c
@@ -0,0 +1,38 @@
+/* Copyright (C) |Meso|Star> 2018 (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/>. */
+
+#include "htgop.h"
+#include "htgop_c.h"
+
+#include <rsys/algorithm.h>
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+/* Generate the htgop_layer_fetch_lw_ka function */
+#define LAYER_SPECINT_DOMAIN lw
+#define LAYER_SPECINT_DATA ka
+#include "htgop_fetch_radiative_properties.h"
+
+/* Generate the htgop_layer_fetch_sw_ka function */
+#define LAYER_SPECINT_DOMAIN sw
+#define LAYER_SPECINT_DATA ka
+#include "htgop_fetch_radiative_properties.h"
+
+/* Generate the htgop_layer_fetch_sw_ks function */
+#define LAYER_SPECINT_DOMAIN sw
+#define LAYER_SPECINT_DATA ks
+#include "htgop_fetch_radiative_properties.h"
+
diff --git a/src/htgop_fetch_radiative_properties.h b/src/htgop_fetch_radiative_properties.h
@@ -0,0 +1,101 @@
+/* Copyright (C) |Meso|Star> 2018 (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/>. */
+
+#if !defined(LAYER_SPECINT_DATA) || !defined(LAYER_SPECINT_DOMAIN)
+ #error "Missing the LAYER_SPECINT_<DATA|DOMAIN> macro."
+#endif
+
+/*
+ * Generate the functions that fetch and interpolate a radiative property in a
+ * given spectral domain
+ */
+
+/* Helper macros */
+#define SPECINT \
+ CONCAT(CONCAT(htgop_layer_, LAYER_SPECINT_DOMAIN),_spectral_interval)
+#define SPECINT_TAB \
+ CONCAT(CONCAT(htgop_layer_, LAYER_SPECINT_DOMAIN),_spectral_interval_tab)
+#define GET_SPECINT \
+ CONCAT(CONCAT(htgop_layer_get_, LAYER_SPECINT_DOMAIN), _spectral_interval)
+#define GET_SPECINT_TAB \
+ CONCAT(CONCAT(htgop_layer_, LAYER_SPECINT_DOMAIN), _spectral_interval_get_tab)
+#define K_tab CONCAT(LAYER_SPECINT_DATA, _tab)
+
+res_T
+CONCAT(CONCAT(CONCAT(htgop_layer_fetch_,LAYER_SPECINT_DOMAIN),_),LAYER_SPECINT_DATA)
+ (const struct htgop_layer* layer,
+ const size_t ispecint,
+ const size_t iquad_point,
+ const double x_h2o, /* H2O mol / mixture mol */
+ double* out_k)
+{
+ struct SPECINT specint;
+ struct SPECINT_TAB tab;
+ double* find;
+ double x1, x2, k1, k2;
+ double k;
+ res_T res = RES_OK;
+
+ if(!layer || !out_k) return RES_BAD_ARG;
+
+ find = search_lower_bound(&x_h2o, layer->x_h2o_tab, layer->tab_length,
+ sizeof(double), cmp_dbl);
+ if(!find) return RES_BAD_ARG;
+
+ res = GET_SPECINT(layer, ispecint, &specint);
+ if(res != RES_OK) return res;
+ res = GET_SPECINT_TAB(&specint, iquad_point, &tab);
+ if(res != RES_OK) return res;
+
+ ASSERT(layer->tab_length);
+ ASSERT(tab.tab_length == layer->tab_length);
+
+ if(find == layer->x_h2o_tab) {
+ /* The submitted x_h2o is smaller that the first tabulated value of the
+ * water vapor molar fraction. Use a simple linear interpolation to define
+ * K by assuming that k1 and x1 is null */
+ ASSERT(layer->x_h2o_tab[0] >= x_h2o);
+ x2 = layer->x_h2o_tab[0];
+ k2 = tab.K_tab[0];
+ k = k2 / x2*x_h2o;
+ } else {
+ const size_t i = (size_t)(find - layer->x_h2o_tab);
+ ASSERT(i < layer->tab_length);
+ ASSERT(layer->x_h2o_tab[i-0] >= x_h2o);
+ ASSERT(layer->x_h2o_tab[i-1] < x_h2o);
+ x1 = layer->x_h2o_tab[i-1];
+ x2 = layer->x_h2o_tab[i-0];
+ k1 = tab.K_tab[i-1];
+ k2 = tab.K_tab[i-0];
+
+ if(!k1 || !k2) { /* Linear interpolation */
+ k = k1 + (k2-k1)/(x2-x1) * (x_h2o-x1);
+ } else {
+ const double alpha = (log(k2) - log(k1))/(log(x1) - log(x2));
+ const double beta = log(k1) - alpha * log(x1);
+ k = exp(alpha * log(x_h2o) + beta);
+ }
+ }
+ *out_k = k;
+ return RES_OK;
+}
+
+#undef SPECINT
+#undef SPECINT_TAB
+#undef GET_SPECINT
+#undef GET_SPECINT_TAB
+#undef K_tab
+#undef LAYER_SPECINT_DATA
+#undef LAYER_SPECINT_DOMAIN
diff --git a/src/htgop_parse_layers_spectral_intervals_data.h b/src/htgop_parse_layers_spectral_intervals_data.h
@@ -14,12 +14,12 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
/*
- * Generate the function parsing a per layer and per spectral interval data for
- * a given domain long/short wave.
+ * Generate the function that parses a per layer and per spectral interval data
+ * for a given domain.
*/
#if !defined(LAYER_SPECINT_DATA) || !defined(LAYER_SPECINT_DOMAIN)
- #error "Missing the LAYER_SPECINT_<DATA|XDOMAIN> macro."
+ #error "Missing the LAYER_SPECINT_<DATA|DOMAIN> macro."
#endif
#define DATA CONCAT(CONCAT(LAYER_SPECINT_DATA, _), LAYER_SPECINT_DOMAIN)
diff --git a/src/htgop_sample_sw_spectral_interval_CIE_1931_XYZ.c b/src/htgop_sample_sw_spectral_interval_CIE_1931_XYZ.c
@@ -0,0 +1,285 @@
+/* Copyright (C) |Meso|Star> 2018 (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 support */
+
+#include "htgop.h"
+#include "htgop_c.h"
+
+#include <rsys/algorithm.h>
+
+/* Helper macros to convert from wave number to wavelength and vice versa */
+#define WNUM_TO_WLEN(Num) (1.0e7 / ((double)Num))
+#define WLEN_TO_WNUM(Len) WNUM_TO_WLEN(Len)
+
+/* Boundaries of the CIE XYZ color space in wavelength and wave number */
+#define CIE_XYZ_WAVELENGTH_MIN 380.0 /* In nanometer */
+#define CIE_XYZ_WAVELENGTH_MAX 780.0 /* In nanometer */
+#define CIE_XYZ_WAVENUMBER_MIN WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MAX) /*In cm^-1*/
+#define CIE_XYZ_WAVENUMBER_MAX WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MIN) /*In cm^-1*/
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static FINLINE double
+trapezoidal_integration
+ (const double lambda_lo, /* Integral lower bound. In nanometer */
+ const double lambda_hi, /* Integral upper bound. In nanometer */
+ double (*f_bar)(const double lambda)) /* Function to integrate */
+{
+ double dlambda;
+ size_t i, n;
+ double integral = 0;
+ ASSERT(lambda_lo <= lambda_hi);
+ ASSERT(lambda_lo > 0);
+
+ n = (size_t)(lambda_hi - lambda_lo) + 1;
+ dlambda = (lambda_hi - lambda_lo) / (double)n;
+
+ FOR_EACH(i, 0, n) {
+ const double lambda1 = lambda_lo + dlambda*(double)(i+0);
+ const double lambda2 = lambda_hi + dlambda*(double)(i+1);
+ const double f1 = f_bar(lambda1);
+ const double f2 = f_bar(lambda2);
+ integral += (f1 + f2)*dlambda*0.5;
+ }
+ return integral;
+}
+
+/* The following 3 functions are used to fit the CIE Xbar, Ybar and Zbar curved
+ * has defined by the 1931 standard. These analytical fits are propsed by C.
+ * Wyman, P. P. Sloan & P. Shirley in "Simple Analytic Approximations to the
+ * CIE XYZ Color Matching Functions" - JCGT 2013. */
+static FINLINE double
+fit_x_bar_1931(const double lambda)
+{
+ const double a = (lambda - 442.0) * (lambda < 442.0 ? 0.0624 : 0.0374);
+ const double b = (lambda - 599.8) * (lambda < 599.8 ? 0.0264 : 0.0323);
+ const double c = (lambda - 501.1) * (lambda < 501.1 ? 0.0490 : 0.0382);
+ return 0.362*exp(-0.5*a*a) + 1.056*exp(-0.5f*b*b) - 0.065*exp(-0.5*c*c);
+}
+
+static FINLINE double
+fit_y_bar_1931(const double lambda)
+{
+ const double a = (lambda - 568.8) * (lambda < 568.8 ? 0.0213 : 0.0247);
+ const double b = (lambda - 530.9) * (lambda < 530.9 ? 0.0613 : 0.0322);
+ return 0.821*exp(-0.5*a*a) + 0.286*exp(-0.5*b*b);
+}
+
+static FINLINE double
+fit_z_bar_1931(const double lambda)
+{
+ const double a = (lambda - 437.0) * (lambda < 437.0 ? 0.0845 : 0.0278);
+ const double b = (lambda - 459.0) * (lambda < 459.0 ? 0.0385 : 0.0725);
+ return 1.217*exp(-0.5*a*a) + 0.681*exp(-0.5*b*b);
+}
+
+static INLINE res_T
+sample_sw_spectral_interval_CIE
+ (const struct htgop* htgop,
+ const double* cdf,
+ const size_t cdf_length,
+ const double r, /* Canonical number in [0, 1[ */
+ size_t* id)
+{
+ double r_next = nextafter(r, DBL_MAX);
+ double* find;
+ size_t i;
+ ASSERT(htgop && cdf && cdf_length);
+
+ if(!id) return RES_BAD_ARG;
+ if(r < 0 || r >= 1) {
+ log_err(htgop, "Could nont sample a spectral interval in the CIE XYZ color "
+ "space: invalid canonical random number `%g'.\n", r);
+ return RES_BAD_ARG;
+ }
+
+ /* Use r_next rather than r in order to find the first entry that is not less
+ * than *or equal* to r */
+ find = search_lower_bound(&r_next, cdf, cdf_length, sizeof(double), cmp_dbl);
+ if(!find) { /* <=> CIE XYZ band does not intersect the loaded intervals */
+ log_err(htgop, "Could not sample a spectral interval in the CIE XYZ color "
+ "space : the loaded data does not overlap the CIE XYZ color space.\n");
+ return RES_BAD_OP;
+ }
+
+ i = (size_t)(find - cdf);
+ ASSERT(i < cdf_length && cdf[i] > r && (!i || cdf[i-1] <= r));
+ *id = i;
+ return RES_OK;
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+htgop_sample_sw_spectral_interval_CIE_1931_X
+ (const struct htgop* htgop,
+ const double r, /* Canonical number in [0, 1[ */
+ size_t* ispecint)
+{
+ if(!htgop) return RES_BAD_ARG;
+ return sample_sw_spectral_interval_CIE
+ (htgop, darray_double_cdata_get(&htgop->sw_X_cdf),
+ darray_double_size_get(&htgop->sw_X_cdf), r, ispecint);
+}
+
+res_T
+htgop_sample_sw_spectral_interval_CIE_1931_Y
+ (const struct htgop* htgop,
+ const double r, /* Canonical number in [0, 1[ */
+ size_t* ispecint)
+{
+ if(!htgop) return RES_BAD_ARG;
+ return sample_sw_spectral_interval_CIE
+ (htgop, darray_double_cdata_get(&htgop->sw_Y_cdf),
+ darray_double_size_get(&htgop->sw_Y_cdf), r, ispecint);
+}
+
+res_T
+htgop_sample_sw_spectral_interval_CIE_1931_Z
+ (const struct htgop* htgop,
+ const double r, /* Canonical number in [0, 1[ */
+ size_t* ispecint)
+{
+ if(!htgop) return RES_BAD_ARG;
+ return sample_sw_spectral_interval_CIE
+ (htgop, darray_double_cdata_get(&htgop->sw_Z_cdf),
+ darray_double_size_get(&htgop->sw_Z_cdf), r, ispecint);
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+compute_CIE_1931_XYZ_cdf(struct htgop* htgop)
+{
+ enum { X, Y, Z }; /* Helper constant */
+ size_t iband_min, iband_max;
+ const double* wnums = NULL;
+ double sum[3] = {0,0,0};
+ double* pdf[3] = {NULL, NULL, NULL};
+ double* cdf[3] = {NULL, NULL, NULL};
+ size_t i, n;
+ res_T res = RES_OK;
+ ASSERT(htgop);
+
+ HTGOP(get_sw_spectral_intervals_wave_numbers(htgop, &wnums));
+ HTGOP(get_sw_spectral_intervals_count(htgop, &n));
+ ASSERT(n);
+
+ /* Allocate and reset the memory space for the tristimulus CDF */
+ #define SETUP_STIMULUS(Stimulus) { \
+ res = darray_double_resize(&htgop->sw_ ## Stimulus ## _cdf, n); \
+ if(res != RES_OK) { \
+ log_err(htgop, \
+ "Could not reserve the memory space for the CDF " \
+ "of the "STR(X)" stimulus.\n"); \
+ goto error; \
+ } \
+ memset(darray_double_data_get(&htgop->sw_ ## Stimulus ## _cdf), 0, \
+ n*sizeof(double)); \
+ cdf[Stimulus] = darray_double_data_get(&htgop->sw_ ## Stimulus ## _cdf); \
+ pdf[Stimulus] = cdf[Stimulus]; \
+ } (void)0
+ SETUP_STIMULUS(X);
+ SETUP_STIMULUS(Y);
+ SETUP_STIMULUS(Z);
+ #undef RESERVE
+
+ ASSERT(eq_eps(CIE_XYZ_WAVENUMBER_MAX, WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MIN), 1.e-6));
+ ASSERT(eq_eps(CIE_XYZ_WAVENUMBER_MIN, WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MAX), 1.e-6));
+ ASSERT(eq_eps(CIE_XYZ_WAVELENGTH_MAX, WNUM_TO_WLEN(CIE_XYZ_WAVENUMBER_MIN), 1.e-6));
+ ASSERT(eq_eps(CIE_XYZ_WAVELENGTH_MIN, WNUM_TO_WLEN(CIE_XYZ_WAVENUMBER_MAX), 1.e-6));
+
+ /* The sw spectral intervals does not intersect the CIE XYZ interval */
+ if(wnums[0] > CIE_XYZ_WAVENUMBER_MAX || wnums[n] < CIE_XYZ_WAVENUMBER_MIN) {
+ log_warn(htgop,
+ "The loaded SW spectral intervals does not intersect "
+ "the CIE XYZ interval.\n");
+ goto exit;
+ }
+
+ /* Search the lower and upper bound of the bands that overlap the XYZ
+ * spectral interval. Use a simple linear search since the overall number of
+ * bands should be small */
+ iband_min = 0, iband_max = n-1;
+ FOR_EACH(i, 0, n) {
+ if(wnums[i]<=CIE_XYZ_WAVENUMBER_MIN && CIE_XYZ_WAVENUMBER_MIN<=wnums[i+1])
+ iband_min = i;
+ if(wnums[i]<=CIE_XYZ_WAVENUMBER_MAX && CIE_XYZ_WAVENUMBER_MAX<=wnums[i+1])
+ iband_max = i;
+ }
+ ASSERT(iband_min <= iband_max);
+
+ /* Compute the *unormalized* pdf of the tristimulus */
+ FOR_EACH(i, iband_min, iband_max+1) {
+ const double lambda_lo = WNUM_TO_WLEN(wnums[i+1]);
+ const double lambda_hi = WNUM_TO_WLEN(wnums[i+0]);
+ ASSERT(lambda_lo <= lambda_hi);
+ pdf[X][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_x_bar_1931);
+ pdf[Y][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_y_bar_1931);
+ pdf[Z][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_z_bar_1931);
+ sum[X] += pdf[X][i];
+ sum[Y] += pdf[Y][i];
+ sum[Z] += pdf[Z][i];
+ }
+
+ /* Compute the cdf of the tristimulus. Note that the upper bound is 'n' and not
+ * 'iband_max+1' in order to ensure that the CDF is strictly increasing. */
+ FOR_EACH(i, iband_min, n) {
+ /* Normalize the pdf */
+ pdf[X][i] /= sum[X];
+ pdf[Y][i] /= sum[Y];
+ pdf[Z][i] /= sum[Z];
+ /* Setup the cumulative */
+ if(i == 0) {
+ cdf[X][i] = pdf[X][i];
+ cdf[Y][i] = pdf[Y][i];
+ cdf[Z][i] = pdf[Z][i];
+ } else {
+ cdf[X][i] = pdf[X][i] + cdf[X][i-1];
+ cdf[Y][i] = pdf[Y][i] + cdf[Y][i-1];
+ cdf[Z][i] = pdf[Z][i] + cdf[Z][i-1];
+ }
+ }
+ ASSERT(eq_eps(cdf[X][n-1], 1, 1.e-6));
+ ASSERT(eq_eps(cdf[Y][n-1], 1, 1.e-6));
+ ASSERT(eq_eps(cdf[Z][n-1], 1, 1.e-6));
+
+#ifndef NDEBUG
+ FOR_EACH(i, 1, n) {
+ ASSERT(cdf[X][i] >= cdf[X][i-1]);
+ ASSERT(cdf[Y][i] >= cdf[Y][i-1]);
+ ASSERT(cdf[Z][i] >= cdf[Z][i-1]);
+ }
+#endif
+
+ /* Handle numerical issue */
+ cdf[X][n-1] = 1.0;
+ cdf[Y][n-1] = 1.0;
+ cdf[Z][n-1] = 1.0;
+
+exit:
+ return res;
+error:
+ darray_double_clear(&htgop->sw_X_cdf);
+ darray_double_clear(&htgop->sw_Y_cdf);
+ darray_double_clear(&htgop->sw_Z_cdf);
+ goto exit;
+}
+
+
diff --git a/src/test_htgop_sample.c b/src/test_htgop_sample.c
@@ -104,7 +104,6 @@ check_sample_quadrature
}
}
-
int
main(int argc, char** argv)
{
@@ -135,9 +134,12 @@ main(int argc, char** argv)
CHK(hist_Y = MEM_CALLOC(&allocator, nspecints, sizeof(*hist_Y)));
CHK(hist_Z = MEM_CALLOC(&allocator, nspecints, sizeof(*hist_Z)));
- check_sw_sample_specint(htgop, htgop_sample_sw_spectral_interval_CIE_X, hist_X);
- check_sw_sample_specint(htgop, htgop_sample_sw_spectral_interval_CIE_Y, hist_Y);
- check_sw_sample_specint(htgop, htgop_sample_sw_spectral_interval_CIE_Z, hist_Z);
+ check_sw_sample_specint
+ (htgop, htgop_sample_sw_spectral_interval_CIE_1931_X, hist_X);
+ check_sw_sample_specint
+ (htgop, htgop_sample_sw_spectral_interval_CIE_1931_Y, hist_Y);
+ check_sw_sample_specint
+ (htgop, htgop_sample_sw_spectral_interval_CIE_1931_Z, hist_Z);
FOR_EACH(ispecint, 0, nspecints) {
if(wnums[ispecint+0] > CIE_XYZ_WAVENUMBER_MAX
