commit 54f3f18b8ee8ca8e2fe130102bbb86e8d3a71888
parent 98ff91202ebc9388e8c5dd4481fe9c1cf3326da2
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 8 Aug 2018 12:38:26 +0200
Full rewrite of the fetch/get radiative properties API
Diffstat:
9 files changed, 316 insertions(+), 310 deletions(-)
diff --git a/src/htgop.c b/src/htgop.c
@@ -537,6 +537,7 @@ htgop_layer_get_lw_spectral_interval
specint->quadrature_length = darray_double_size_get(&lspecint->ka_nominal);
specint->htgop = layer->htgop;
specint->data__ = lspecint;
+ specint->layer__ = l;
return RES_OK;
}
@@ -561,6 +562,7 @@ htgop_layer_get_sw_spectral_interval
specint->quadrature_length = darray_double_size_get(&lspecint->ka_nominal);
specint->htgop = layer->htgop;
specint->data__ = lspecint;
+ specint->layer__ = l;
return RES_OK;
}
@@ -571,6 +573,7 @@ htgop_layer_lw_spectral_interval_get_tab
struct htgop_layer_lw_spectral_interval_tab* tab)
{
const struct layer_lw_spectral_interval* lspecint;
+ const struct layer* l;
const struct darray_double* t;
if(!specint || !tab) return RES_BAD_ARG;
if(iquad >= specint->quadrature_length) {
@@ -579,7 +582,9 @@ htgop_layer_lw_spectral_interval_get_tab
return RES_BAD_ARG;
}
lspecint = specint->data__;
+ l = specint->layer__;
t = darray_dbllst_cdata_get(&lspecint->ka_tab) + iquad;
+ tab->x_h2o_tab = darray_double_cdata_get(&l->x_h2o_tab);
tab->ka_tab = darray_double_cdata_get(t);
tab->tab_length = darray_double_size_get(t);
tab->htgop = specint->htgop;
@@ -593,6 +598,7 @@ htgop_layer_sw_spectral_interval_get_tab
struct htgop_layer_sw_spectral_interval_tab* tab)
{
const struct layer_sw_spectral_interval* lspecint;
+ const struct layer* l;
const struct darray_double* ka_tab;
const struct darray_double* ks_tab;
if(!specint || !tab) return RES_BAD_ARG;
@@ -602,8 +608,11 @@ htgop_layer_sw_spectral_interval_get_tab
return RES_BAD_ARG;
}
lspecint = specint->data__;
+ l = specint->layer__;
+
ka_tab = darray_dbllst_cdata_get(&lspecint->ka_tab) + iquad;
ks_tab = darray_dbllst_cdata_get(&lspecint->ks_tab) + iquad;
+ tab->x_h2o_tab = darray_double_cdata_get(&l->x_h2o_tab);
tab->ka_tab = darray_double_cdata_get(ka_tab);
tab->ks_tab = darray_double_cdata_get(ks_tab);
tab->tab_length = darray_double_size_get(ka_tab);
@@ -785,6 +794,12 @@ htgop_position_to_layer_id
#define DATA ks
#include "htgop_fetch_radiative_properties.h"
+/* Generate the htgop_layer_get_sw_kext function */
+#define GET_DATA(Tab, Id) ((Tab)->ka_tab[Id] + (Tab)->ks_tab[Id])
+#define DOMAIN sw
+#define DATA kext
+#include "htgop_fetch_radiative_properties.h"
+
/* Generate the functions that get the boundaries of ka in LW */
#define DOMAIN lw
#define DATA ka
diff --git a/src/htgop.h b/src/htgop.h
@@ -77,10 +77,13 @@ struct htgop_layer_lw_spectral_interval {
const double* ka_nominal; /* Per quadrature point absorption coef */
size_t quadrature_length;
const struct htgop* htgop;
- const void* data__; /* Internal data */
+ /* Internal data */
+ const void* data__;
+ const void* layer__;
};
struct htgop_layer_lw_spectral_interval_tab {
+ const double* x_h2o_tab; /* <=> layer.x_h2o_tab */
const double* ka_tab; /* Tabulated absorption coef */
size_t tab_length; /* == lengh of the tabulated xH2O */
const struct htgop* htgop;
@@ -91,10 +94,13 @@ struct htgop_layer_sw_spectral_interval {
const double* ks_nominal; /* Per quadrature point scattering coef */
size_t quadrature_length;
const struct htgop* htgop;
- const void* data__; /* Internal data */
+ /* Internal data */
+ const void* data__;
+ const void* layer__;
};
struct htgop_layer_sw_spectral_interval_tab {
+ const double* x_h2o_tab; /* <=> layer.x_h2o_tab */
const double* ka_tab; /* Tabulated absorption coef */
const double* ks_tab; /* Tabulated scattering coef */
size_t tab_length; /* == lengh of the tabulated xH2O */
@@ -245,29 +251,29 @@ htgop_layer_sw_spectral_interval_get_tab
* beta = ln(k1) - alpha * ln(x1)
******************************************************************************/
HTGOP_API res_T
-htgop_layer_fetch_lw_ka
- (const struct htgop_layer* layer,
- const size_t ispecint,
- const size_t iquad_point,
+htgop_layer_lw_spectral_interval_tab_fetch_ka
+ (const struct htgop_layer_lw_spectral_interval_tab* tab,
const double x_h2o, /* H2O mol / mixture mol */
double* ka);
HTGOP_API res_T
-htgop_layer_fetch_sw_ka
- (const struct htgop_layer* layer,
- const size_t ispecint,
- const size_t iquad_point,
+htgop_layer_sw_spectral_interval_tab_fetch_ka
+ (const struct htgop_layer_sw_spectral_interval_tab* tab,
const double x_h2o, /* H2O mol / mixture mol */
double* ka);
HTGOP_API res_T
-htgop_layer_fetch_sw_ks
- (const struct htgop_layer* layer,
- const size_t ispecint,
- const size_t iquad_point,
+htgop_layer_sw_spectral_interval_tab_fetch_ks
+ (const struct htgop_layer_sw_spectral_interval_tab* tab,
const double x_h2o, /* H2O mol / mixture mol */
double* ks);
+HTGOP_API res_T
+htgop_layer_sw_spectral_interval_tab_fetch_kext
+ (const struct htgop_layer_sw_spectral_interval_tab* tab,
+ const double x_h2o, /* H2O mol / mixture mol */
+ double* kext);
+
/*******************************************************************************
* Sample of a short wave spectral interval according to the CIE 1931 XYZ
* tristimulus values.
@@ -294,75 +300,55 @@ htgop_sample_sw_spectral_interval_CIE_1931_Z
* Retrieve the boundaries of the radiative properties
******************************************************************************/
HTGOP_API res_T
-htgop_layers_get_lw_ka_bounds
- (const struct htgop* htgop,
- const size_t ilayer_range[2], /* Range of layers to handle */
- const size_t ispecint_range[2], /* Range of spectral intervals to handle */
- double bounds[2]);
-
-HTGOP_API res_T
-htgop_layers_get_sw_ka_bounds
- (const struct htgop* htgop,
- const size_t ilayer_range[2], /* Range of layers to handle */
- const size_t ispecint_range[2], /* Range of spectral intervals to handle */
- double bounds[2]);
-
-HTGOP_API res_T
-htgop_layers_get_sw_ks_bounds
- (const struct htgop* htgop,
- const size_t ilayer_range[2], /* Range of layers to handle */
- const size_t ispecint_range[2], /* Range of spectral intervals to handle */
- double bounds[2]);
-
-HTGOP_API res_T
-htgop_layers_get_sw_kext_bounds
- (const struct htgop* htgop,
- const size_t ilayer_range[2], /* Range of layers to handle */
- const size_t ispecint_range[2], /* Range of spectral intervals to handle */
- double bounds[2]);
-
-HTGOP_API res_T
htgop_layer_lw_spectral_interval_quadpoints_get_ka_bounds
(const struct htgop_layer_lw_spectral_interval* specint,
const size_t iquad_range[2], /* Range of quadrature points to handle */
+ const double x_h2o_range[2],
double bounds[2]);
HTGOP_API res_T
htgop_layer_sw_spectral_interval_quadpoints_get_ka_bounds
(const struct htgop_layer_sw_spectral_interval* specint,
const size_t iquad_range[2], /* Range of quadrature points to handle */
+ const double x_h2o_range[2],
double bounds[2]);
HTGOP_API res_T
htgop_layer_sw_spectral_interval_quadpoints_get_ks_bounds
(const struct htgop_layer_sw_spectral_interval* specint,
const size_t iquad_range[2], /* Range of quadrature points to handle */
+ const double x_h2o_range[2],
double bounds[2]);
HTGOP_API res_T
htgop_layer_sw_spectral_interval_quadpoints_get_kext_bounds
(const struct htgop_layer_sw_spectral_interval* specint,
const size_t iquad_range[2], /* Range of quadrature points to handle */
+ const double x_h2o_range[2],
double bounds[2]);
HTGOP_API res_T
htgop_layer_lw_spectral_interval_tab_get_ka_bounds
(const struct htgop_layer_lw_spectral_interval_tab* tab,
+ const double x_h2o_range[2],
double bounds[2]);
HTGOP_API res_T
htgop_layer_sw_spectral_interval_tab_get_ka_bounds
(const struct htgop_layer_sw_spectral_interval_tab* tab,
+ const double x_h2o_range[2],
double bounds[2]);
HTGOP_API res_T
htgop_layer_sw_spectral_interval_tab_get_ks_bounds
(const struct htgop_layer_sw_spectral_interval_tab* tab,
+ const double x_h2o_range[2],
double bounds[2]);
HTGOP_API res_T
htgop_layer_sw_spectral_interval_tab_get_kext_bounds
(const struct htgop_layer_sw_spectral_interval_tab* tab,
+ const double x_h2o_range[2],
double bounds[2]);
/*******************************************************************************
@@ -385,7 +371,7 @@ htgop_spectral_interval_sample_quadrature
* color space */
HTGOP_API res_T
htgop_get_sw_spectral_intervals_CIE_XYZ
- (const struct htgop* htgop,
+ (const struct htgop* htgop,
size_t specint_raneg[2]);
END_DECLS
diff --git a/src/htgop_c.h b/src/htgop_c.h
@@ -80,5 +80,29 @@ extern LOCAL_SYM res_T
compute_CIE_1931_XYZ_cdf
(struct htgop* htgop);
+extern LOCAL_SYM double
+layer_lw_spectral_interval_tab_interpolate_ka
+ (const struct htgop_layer_lw_spectral_interval_tab* tab,
+ const double x_h2o,
+ const double* x_h2o_upp); /* Pointer toward the 1st tabbed xH2O >= x_h2o */
+
+extern LOCAL_SYM double
+layer_sw_spectral_interval_tab_interpolate_ka
+ (const struct htgop_layer_sw_spectral_interval_tab* tab,
+ const double x_h2o,
+ const double* x_h2o_upp); /* Pointer toward the 1st tabbed xH2O >= x_h2o */
+
+extern LOCAL_SYM double
+layer_sw_spectral_interval_tab_interpolate_ks
+ (const struct htgop_layer_sw_spectral_interval_tab* tab,
+ const double x_h2o,
+ const double* x_h2o_upp); /* Pointer toward the 1st tabbed xH2O >= x_h2o */
+
+extern LOCAL_SYM double
+layer_sw_spectral_interval_tab_interpolate_kext
+ (const struct htgop_layer_sw_spectral_interval_tab* tab,
+ const double x_h2o,
+ const double* x_h2o_upp); /* Pointer toward the 1st tabbed xH2O >= x_h2o */
+
#endif /* HTGOP_C_H */
diff --git a/src/htgop_fetch_radiative_properties.h b/src/htgop_fetch_radiative_properties.h
@@ -27,59 +27,68 @@
* given spectral domain
*/
-/* Helper macros */
-#define SPECINT CONCAT(CONCAT(htgop_layer_, DOMAIN),_spectral_interval)
-#define SPECINT_TAB CONCAT(CONCAT(htgop_layer_, DOMAIN),_spectral_interval_tab)
-#define GET_SPECINT CONCAT(CONCAT(htgop_layer_get_, DOMAIN), _spectral_interval)
-#define GET_SPECINT_TAB CONCAT(CONCAT(htgop_layer_, DOMAIN), _spectral_interval_get_tab)
-#define K_tab CONCAT(DATA, _tab)
+#ifndef GET_DATA
+ #define GET_DATA(Tab, Id) ((Tab)->CONCAT(DATA, _tab)[Id])
+#endif
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
res_T
-CONCAT(CONCAT(CONCAT(htgop_layer_fetch_,DOMAIN),_),DATA)
- (const struct htgop_layer* layer,
- const size_t ispecint,
- const size_t iquad_point,
+CONCAT(CONCAT(CONCAT(htgop_layer_,DOMAIN),_spectral_interval_tab_fetch_),DATA)
+ (const struct CONCAT(CONCAT(htgop_layer_,DOMAIN),_spectral_interval_tab)* tab,
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;
+ if(!tab || !out_k) return RES_BAD_ARG;
- find = search_lower_bound(&x_h2o, layer->x_h2o_tab, layer->tab_length,
+ find = search_lower_bound(&x_h2o, tab->x_h2o_tab, tab->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;
+ *out_k = CONCAT(CONCAT(CONCAT(
+ layer_,DOMAIN),_spectral_interval_tab_interpolate_),DATA)(tab, x_h2o, find);
+ return RES_OK;
+}
- ASSERT(layer->tab_length);
- ASSERT(tab.tab_length == layer->tab_length);
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+double
+CONCAT(CONCAT(CONCAT(layer_,DOMAIN),_spectral_interval_tab_interpolate_),DATA)
+ (const struct CONCAT(CONCAT(htgop_layer_,DOMAIN),_spectral_interval_tab)* tab,
+ const double x_h2o,
+ const double* x_h2o_upp) /* Pointer toward the upper bound used to interpolate */
+{
+ double x1, x2, k1, k2;
+ double k;
+
+ ASSERT(tab && x_h2o_upp);
- if(find == layer->x_h2o_tab) {
+ ASSERT(x_h2o_upp >= tab->x_h2o_tab);
+ ASSERT(x_h2o_upp < tab->x_h2o_tab + tab->tab_length);
+ ASSERT(*x_h2o_upp >= x_h2o);
+ ASSERT(x_h2o_upp == tab->x_h2o_tab || x_h2o_upp[-1] < x_h2o);
+
+ if(x_h2o_upp == tab->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];
+ ASSERT(tab->x_h2o_tab[0] >= x_h2o);
+ x2 = tab->x_h2o_tab[0];
+ k2 = GET_DATA(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];
+ const size_t i = (size_t)(x_h2o_upp - tab->x_h2o_tab);
+ ASSERT(i < tab->tab_length);
+ ASSERT(tab->x_h2o_tab[i-0] >= x_h2o);
+ ASSERT(tab->x_h2o_tab[i-1] < x_h2o);
+ x1 = tab->x_h2o_tab[i-1];
+ x2 = tab->x_h2o_tab[i-0];
+ k1 = GET_DATA(tab, i-1);
+ k2 = GET_DATA(tab, i-0);
if(!k1 || !k2) { /* Linear interpolation */
k = k1 + (k2-k1)/(x2-x1) * (x_h2o-x1);
@@ -89,14 +98,9 @@ CONCAT(CONCAT(CONCAT(htgop_layer_fetch_,DOMAIN),_),DATA)
k = exp(alpha * log(x_h2o) + beta);
}
}
- *out_k = k;
- return RES_OK;
+ return k;
}
-#undef SPECINT
-#undef SPECINT_TAB
-#undef GET_SPECINT
-#undef GET_SPECINT_TAB
-#undef K_tab
+#undef GET_DATA
#undef DATA
#undef DOMAIN
diff --git a/src/htgop_get_radiative_properties_bounds.h b/src/htgop_get_radiative_properties_bounds.h
@@ -28,100 +28,27 @@
* Generate functions that retrieve the boundaries of the radiative properties
*/
-res_T
-CONCAT(CONCAT(CONCAT(CONCAT(htgop_layers_get_,DOMAIN),_),DATA),_bounds)
- (const struct htgop* htgop,
- const size_t ilayer_range[2],
- const size_t ispecint_range[2],
- double bounds[2])
-{
- size_t ilayer;
- size_t nlayers;
- size_t nspecints;
- res_T res = RES_OK;
-
- if(!htgop || !bounds || !ilayer_range || !ispecint_range || !bounds) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- HTGOP(get_layers_count(htgop, &nlayers));
- if(ilayer_range[0] > ilayer_range[1]
- || ilayer_range[0] >= nlayers
- || ilayer_range[1] >= nlayers) {
- log_err(htgop, "%s: invalid range of layers [%lu, %lu].\n", FUNC_NAME,
- (unsigned long)ilayer_range[0], (unsigned long)ilayer_range[1]);
- res = RES_BAD_ARG;
- goto error;
- }
-
- HTGOP(CONCAT(CONCAT(get_,DOMAIN),_spectral_intervals_count)(htgop, &nspecints));
- if(ispecint_range[0] > ispecint_range[1]
- || ispecint_range[0] >= nspecints
- || ispecint_range[1] >= nspecints) {
- log_err(htgop, "%s: invalid range of spectral intervals [%lu, %lu].\n",
- FUNC_NAME,
- (unsigned long)ispecint_range[0],
- (unsigned long)ispecint_range[1]);
- res = RES_BAD_ARG;
- goto error;
- }
-
- bounds[0] = DBL_MAX;
- bounds[1] =-DBL_MAX;
- FOR_EACH(ilayer, ilayer_range[0], ilayer_range[1]+1) {
- struct htgop_layer layer;
- size_t ispecint;
- HTGOP(get_layer(htgop, ilayer, &layer));
-
- FOR_EACH(ispecint, ispecint_range[0], ispecint_range[1]+1) {
- struct CONCAT(CONCAT(htgop_layer_,DOMAIN),_spectral_interval) specint;
- size_t iquad_range[2];
- double quad_bounds[2];
-
- HTGOP(CONCAT(CONCAT(layer_get_,DOMAIN),_spectral_interval)
- (&layer, ispecint, &specint));
-
- iquad_range[0] = 0;
- iquad_range[1] = specint.quadrature_length - 1;
- HTGOP(CONCAT(CONCAT(CONCAT(CONCAT(
- layer_,DOMAIN),_spectral_interval_quadpoints_get_),DATA),_bounds)
- (&specint, iquad_range, quad_bounds));
-
- bounds[0] = MMIN(bounds[0], quad_bounds[0]);
- bounds[1] = MMAX(bounds[1], quad_bounds[1]);
- }
- }
-exit:
- return res;
-error:
- if(bounds) {
- bounds[0] = DBL_MAX;
- bounds[1] =-DBL_MAX;
- }
- goto exit;
-}
-
res_T CONCAT(CONCAT(CONCAT(CONCAT(
htgop_layer_,DOMAIN),_spectral_interval_quadpoints_get_),DATA),_bounds)
(const struct CONCAT(CONCAT(htgop_layer_,DOMAIN),_spectral_interval)* specint,
- const size_t quad_range[2],
+ const size_t iquad_range[2],
+ const double x_h2o_range[2],
double bounds[2])
{
size_t iquad;
res_T res = RES_OK;
- if(!specint || !bounds || !quad_range) {
+ if(!specint || !bounds || !iquad_range) {
res = RES_BAD_ARG;
goto error;
}
- if(quad_range[0] > quad_range[1]
- || quad_range[0] >= specint->quadrature_length
- || quad_range[1] >= specint->quadrature_length) {
+ if(iquad_range[0] > iquad_range[1]
+ || iquad_range[0] >= specint->quadrature_length
+ || iquad_range[1] >= specint->quadrature_length) {
log_err(specint->htgop,
"%s: invalid range of quadrature points [%lu, %lu].\n", FUNC_NAME,
- (unsigned long)quad_range[0], (unsigned long)quad_range[1]);
+ (unsigned long)iquad_range[0], (unsigned long)iquad_range[1]);
res = RES_BAD_ARG;
goto error;
}
@@ -129,15 +56,15 @@ htgop_layer_,DOMAIN),_spectral_interval_quadpoints_get_),DATA),_bounds)
bounds[0] = DBL_MAX;
bounds[1] =-DBL_MAX;
- FOR_EACH(iquad, quad_range[0], quad_range[1]+1) {
+ FOR_EACH(iquad, iquad_range[0], iquad_range[1]+1) {
struct CONCAT(CONCAT(htgop_layer_,DOMAIN),_spectral_interval_tab) tab;
double tab_bounds[2];
HTGOP(CONCAT(CONCAT(layer_,DOMAIN),_spectral_interval_get_tab)
(specint, iquad, &tab));
- HTGOP(CONCAT(CONCAT(CONCAT(CONCAT(
- layer_,DOMAIN),_spectral_interval_tab_get_),DATA),_bounds)
- (&tab, tab_bounds));
+ res = CONCAT(CONCAT(CONCAT(CONCAT(
+ htgop_layer_,DOMAIN),_spectral_interval_tab_get_),DATA),_bounds)
+ (&tab, x_h2o_range, tab_bounds);
bounds[0] = MMIN(bounds[0], tab_bounds[0]);
bounds[1] = MMAX(bounds[1], tab_bounds[1]);
@@ -156,25 +83,68 @@ error:
res_T CONCAT(CONCAT(CONCAT(CONCAT(
htgop_layer_,DOMAIN),_spectral_interval_tab_get_),DATA),_bounds)
(const struct CONCAT(CONCAT(htgop_layer_,DOMAIN),_spectral_interval_tab)* tab,
+ const double x_h2o_range[2],
double bounds[2])
{
+ size_t itab_range[2];
+ double* x_h2o_low;
+ double* x_h2o_upp;
+ double k1, k2;
size_t i;
res_T res = RES_OK;
- if(!tab || !bounds) {
+ if(!tab || !bounds || !x_h2o_range) {
res = RES_BAD_ARG;
goto error;
}
- bounds[0] = DBL_MAX;
- bounds[1] =-DBL_MAX;
- FOR_EACH(i, 0, tab->tab_length) {
+ if(x_h2o_range[0] > x_h2o_range[1]) {
+ log_err(tab->htgop,
+ "%s: invalid water vapor molar fraction range [%g, %g].\n",
+ FUNC_NAME, SPLIT2(x_h2o_range));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Find the pointers toward the first x_h2o greater than or equal to the
+ * x_h2o boundaries */
+ x_h2o_low = search_lower_bound(&x_h2o_range[0], tab->x_h2o_tab,
+ tab->tab_length, sizeof(double), cmp_dbl);
+ x_h2o_upp = search_lower_bound(&x_h2o_range[1], tab->x_h2o_tab,
+ tab->tab_length, sizeof(double), cmp_dbl);
+ if(!x_h2o_low || !x_h2o_upp || (x_h2o_low > x_h2o_upp)) {
+ log_err(tab->htgop, "%s: invalid x_h2o range [%g, %g].\n",
+ FUNC_NAME, SPLIT2(x_h2o_range));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Compute the K at the boundaries of the x_h2o_range */
+ k1 = CONCAT(CONCAT(CONCAT(
+ layer_,DOMAIN),_spectral_interval_tab_interpolate_),DATA)
+ (tab, x_h2o_range[0], x_h2o_low);
+
+ if(x_h2o_range[0] == x_h2o_range[1]) {
+ k2 = k1;
+ } else {
+ k2 = CONCAT(CONCAT(CONCAT(
+ layer_,DOMAIN),_spectral_interval_tab_interpolate_),DATA)
+ (tab, x_h2o_range[1], x_h2o_upp);
+ }
+ bounds[0] = MMIN(k1, k2);
+ bounds[1] = MMAX(k1, k2);
+
+ /* Define the remaining tabulation point to handle */
+ itab_range[0] = (size_t)(x_h2o_low - tab->x_h2o_tab);
+ itab_range[1] = (size_t)(x_h2o_upp - tab->x_h2o_tab);
+
+ /* Iterate on the remaining tabulation point */
+ FOR_EACH(i, itab_range[0], itab_range[1]) {
const double k = GET_DATA(tab, i);
bounds[0] = MMIN(bounds[0], k);
bounds[1] = MMAX(bounds[1], k);
}
- return RES_OK;
exit:
return res;
error:
diff --git a/src/test_htgop_fetch_radiative_properties.c b/src/test_htgop_fetch_radiative_properties.c
@@ -30,6 +30,11 @@
#define DOMAIN sw
#define DATA ks
#include "test_htgop_fetch_radiative_properties.h"
+/* Generate the check_layer_fetch_sw_ks function */
+#define GET_K(Tab, Id) ((Tab)->ka_tab[Id] + (Tab)->ks_tab[Id])
+#define DOMAIN sw
+#define DATA kext
+#include "test_htgop_fetch_radiative_properties.h"
int
main(int argc, char** argv)
@@ -56,6 +61,7 @@ main(int argc, char** argv)
check_layer_fetch_lw_ka(htgop, &layer);
check_layer_fetch_sw_ka(htgop, &layer);
check_layer_fetch_sw_ks(htgop, &layer);
+ check_layer_fetch_sw_kext(htgop, &layer);
}
CHK(htgop_ref_put(htgop) == RES_OK);
diff --git a/src/test_htgop_fetch_radiative_properties.h b/src/test_htgop_fetch_radiative_properties.h
@@ -34,8 +34,11 @@
#define GET_LAY_SPECINT_TAB \
CONCAT(CONCAT(htgop_layer_, DOMAIN), _spectral_interval_get_tab)
#define FETCH_K \
- CONCAT(CONCAT(CONCAT(htgop_layer_fetch_, DOMAIN),_), DATA)
-#define K_TAB CONCAT(DATA, _tab)
+ CONCAT(CONCAT(CONCAT(htgop_layer_, DOMAIN),_spectral_interval_tab_fetch_), DATA)
+
+#ifndef GET_K
+ #define GET_K(Tab, Id) ((Tab)->CONCAT(DATA, _tab)[Id])
+#endif
static void
CONCAT(CONCAT(CONCAT(check_layer_fetch_, DOMAIN),_), DATA)
@@ -45,8 +48,6 @@ CONCAT(CONCAT(CONCAT(check_layer_fetch_, DOMAIN),_), DATA)
struct htgop_spectral_interval specint;
size_t ispecint;
size_t nspecints;
- size_t quad_len;
- double x_h2o;
double k;
CHK(htgop && layer);
@@ -57,15 +58,6 @@ CONCAT(CONCAT(CONCAT(check_layer_fetch_, DOMAIN),_), DATA)
CHK(specint.quadrature_length);
CHK(layer->tab_length);
- quad_len = specint.quadrature_length;
-
- x_h2o = layer->x_h2o_tab[0];
- CHK(FETCH_K(NULL, 0, 0, x_h2o, &k) == RES_BAD_ARG);
- CHK(FETCH_K(layer, nspecints, 0, x_h2o, &k) == RES_BAD_ARG);
- CHK(FETCH_K(layer, 0, quad_len, x_h2o, &k) == RES_BAD_ARG);
- CHK(FETCH_K(layer, 0, 0, 1.0000001, &k) == RES_BAD_ARG);
- CHK(FETCH_K(layer, 0, 0, x_h2o, NULL) == RES_BAD_ARG);
-
FOR_EACH(ispecint, 0, nspecints) {
struct LAY_SPECINT lay_specint;
struct LAY_SPECINT_TAB lay_tab;
@@ -81,21 +73,22 @@ CONCAT(CONCAT(CONCAT(check_layer_fetch_, DOMAIN),_), DATA)
CHK(GET_LAY_SPECINT_TAB(&lay_specint, iquad, &lay_tab) == RES_OK);
CHK(lay_tab.tab_length == layer->tab_length);
- /* Check boundaries */
FOR_EACH(itab, 0, layer->tab_length) {
+ double x_h2o;
const size_t N = 10;
size_t i;
+ /* Check boundaries */
x_h2o = layer->x_h2o_tab[itab];
- CHK(FETCH_K(layer, ispecint, iquad, x_h2o, &k) == RES_OK);
- CHK(eq_eps(k, lay_tab.K_TAB[itab], 1.e-6));
+ CHK(FETCH_K(&lay_tab, x_h2o, &k) == RES_OK);
+ CHK(eq_eps(k, GET_K(&lay_tab,itab), 1.e-6));
FOR_EACH(i, 0, N) {
const double r = rand_canonic();
const double x1 = itab ? layer->x_h2o_tab[itab-1] : 0;
- const double k1 = itab ? lay_tab.K_TAB[itab-1] : 0;
+ const double k1 = itab ? GET_K(&lay_tab, itab-1) : 0;
const double x2 = layer->x_h2o_tab[itab];
- const double k2 = lay_tab.K_TAB[itab];
+ const double k2 = GET_K(&lay_tab, itab);
double ref;
x_h2o = x1 + (x2 - x1) * r;
@@ -107,7 +100,7 @@ CONCAT(CONCAT(CONCAT(check_layer_fetch_, DOMAIN),_), DATA)
ref = exp(alpha*log(x_h2o) + beta);
}
- CHK(FETCH_K(layer, ispecint, iquad, x_h2o, &k) == RES_OK);
+ CHK(FETCH_K(&lay_tab, x_h2o, &k) == RES_OK);
CHK(eq_eps(k, ref, 1.e-6));
}
}
@@ -121,8 +114,7 @@ CONCAT(CONCAT(CONCAT(check_layer_fetch_, DOMAIN),_), DATA)
#undef LAY_SPECINT_TAB
#undef GET_LAY_SPECINT
#undef GET_LAY_SPECINT_TAB
-#undef FETCH_K
-#undef K_TAB
+#undef GET_K
#undef DOMAIN
#undef DATA
diff --git a/src/test_htgop_get_radiative_properties_bounds.h b/src/test_htgop_get_radiative_properties_bounds.h
@@ -25,14 +25,7 @@
CONCAT(CONCAT(htgop_get_, DOMAIN), _spectral_intervals_count)
#define GET_SPECINT \
CONCAT(CONCAT(htgop_get_, DOMAIN), _spectral_interval)
-#define LAY_GET_BOUNDS \
- CONCAT(CONCAT(CONCAT(CONCAT(htgop_layers_get_,DOMAIN),_),DATA),_bounds)
-#define LAY_SPECINT_QUADS_GET_BOUNDS \
- CONCAT(CONCAT(CONCAT(CONCAT( \
- htgop_layer_,DOMAIN),_spectral_interval_quadpoints_get_),DATA),_bounds)
-#define LAY_SPECINT_TAB_GET_BOUNDS \
- CONCAT(CONCAT(CONCAT(CONCAT( \
- htgop_layer_,DOMAIN),_spectral_interval_tab_get_),DATA),_bounds)
+
#define LAY_SPECINT \
CONCAT(CONCAT(htgop_layer_, DOMAIN),_spectral_interval)
#define LAY_SPECINT_TAB \
@@ -42,6 +35,13 @@
#define GET_LAY_SPECINT_TAB \
CONCAT(CONCAT(htgop_layer_, DOMAIN), _spectral_interval_get_tab)
+#define LAY_SPECINT_QUADS_GET_BOUNDS \
+ CONCAT(CONCAT(CONCAT(CONCAT( \
+ htgop_layer_,DOMAIN),_spectral_interval_quadpoints_get_),DATA),_bounds)
+#define LAY_SPECINT_TAB_GET_BOUNDS \
+ CONCAT(CONCAT(CONCAT(CONCAT( \
+ htgop_layer_,DOMAIN),_spectral_interval_tab_get_),DATA),_bounds)
+
#ifndef GET_K
#define GET_K(Tab, Id) ((Tab)->CONCAT(DATA,_tab)[Id])
#endif
@@ -51,150 +51,156 @@ check_layer_,DOMAIN),_),DATA),_bounds)(struct htgop* htgop)
{
struct htgop_layer lay;
struct LAY_SPECINT band;
- size_t nlays, nbands;
- size_t ilay, iband;
- size_t band_range[2];
- size_t lay_range[2];
- size_t quad_range[2];
+ struct LAY_SPECINT_TAB tab;
+ double x_h2o_range[2];
double bounds[2];
-
+ size_t quad_range[2];
+ size_t ilay, iband, iquad, itab;
+ size_t nlays, nbands;
+
CHK(GET_NSPECINTS(htgop, &nbands) == RES_OK);
CHK(htgop_get_layers_count(htgop, &nlays) == RES_OK);
CHK(nlays && nbands);
- lay_range[0] = lay_range[1] = 0;
- band_range[0] = band_range[1] = 0;
- CHK(LAY_GET_BOUNDS(NULL, lay_range, band_range, bounds) == RES_BAD_ARG);
- CHK(LAY_GET_BOUNDS(htgop, NULL, band_range, bounds) == RES_BAD_ARG);
- CHK(LAY_GET_BOUNDS(htgop, lay_range, NULL, bounds) == RES_BAD_ARG);
- CHK(LAY_GET_BOUNDS(htgop, lay_range, band_range, NULL) == RES_BAD_ARG);
-
- lay_range[1] = nlays;
- CHK(LAY_GET_BOUNDS(htgop, lay_range, band_range, bounds) == RES_BAD_ARG);
- lay_range[1] = 0;
- band_range[1] = nlays;
- CHK(LAY_GET_BOUNDS(htgop, lay_range, band_range, bounds) == RES_BAD_ARG);
- band_range[1] = 0;
-
- if(nlays > 1) {
- lay_range[0] = 1;
- CHK(LAY_GET_BOUNDS(htgop, lay_range, band_range, bounds) == RES_BAD_ARG);
- lay_range[0] = 0;
- }
- if(nbands > 1) {
- band_range[0] = 1;
- CHK(LAY_GET_BOUNDS(htgop, lay_range, band_range, bounds) == RES_BAD_ARG);
- band_range[0] = 0;
- }
-
CHK(htgop_get_layer(htgop, 0, &lay) == RES_OK);
CHK(GET_LAY_SPECINT(&lay, 0, &band) == RES_OK);
-
CHK(band.quadrature_length);
- quad_range[0] = quad_range[1] = 0;
- CHK(LAY_SPECINT_QUADS_GET_BOUNDS(NULL, quad_range, bounds) == RES_BAD_ARG);
- CHK(LAY_SPECINT_QUADS_GET_BOUNDS(&band, NULL, bounds) == RES_BAD_ARG);
- CHK(LAY_SPECINT_QUADS_GET_BOUNDS(&band, quad_range, NULL) == RES_BAD_ARG);
+ CHK(GET_LAY_SPECINT_TAB(&band, 0, &tab) == RES_OK);
- quad_range[1] = band.quadrature_length;
- CHK(LAY_SPECINT_QUADS_GET_BOUNDS(&band, quad_range, bounds) == RES_BAD_ARG);
+ quad_range[0] = 0;
quad_range[1] = 0;
- if(band.quadrature_length) {
- quad_range[0] = 1;
- CHK(LAY_SPECINT_QUADS_GET_BOUNDS(&band, quad_range, bounds) == RES_BAD_ARG);
- quad_range[0] = 0;
- }
-
- FOR_EACH(iband, 0, nbands) {
- double band_bounds[2] = {DBL_MAX, -DBL_MAX};
-
- FOR_EACH(ilay, 0, nlays) {
- bounds[0] = DBL_MAX;
- bounds[1] =-DBL_MAX;
-
- band_range[0] = iband;
- band_range[1] = iband;
- lay_range[0] = ilay;
- lay_range[1] = ilay;
- CHK(LAY_GET_BOUNDS(htgop, lay_range, band_range, bounds) == RES_OK);
-
- band_bounds[0] = MMIN(band_bounds[0], bounds[0]);
- band_bounds[1] = MMAX(band_bounds[1], bounds[1]);
-
- lay_range[0] = 0;
- lay_range[1] = ilay;
- CHK(LAY_GET_BOUNDS(htgop, lay_range, band_range, bounds) == RES_OK);
-
- CHK(bounds[0] == band_bounds[0]);
- CHK(bounds[1] == band_bounds[1]);
- }
- }
+ x_h2o_range[0] = 0;
+ x_h2o_range[1] = 0;
+ CHK(LAY_SPECINT_QUADS_GET_BOUNDS
+ (NULL, quad_range, x_h2o_range, bounds) == RES_BAD_ARG);
+ CHK(LAY_SPECINT_QUADS_GET_BOUNDS
+ (&band, NULL, x_h2o_range, bounds) == RES_BAD_ARG);
+ CHK(LAY_SPECINT_QUADS_GET_BOUNDS
+ (&band, quad_range, NULL, bounds) == RES_BAD_ARG);
+ CHK(LAY_SPECINT_QUADS_GET_BOUNDS
+ (&band, quad_range, x_h2o_range, NULL) == RES_BAD_ARG);
+ x_h2o_range[0] = 1;
+ x_h2o_range[1] = 0;
+ CHK(LAY_SPECINT_QUADS_GET_BOUNDS
+ (&band, quad_range, x_h2o_range, bounds) == RES_BAD_ARG);
+ x_h2o_range[0] = 0;
+ x_h2o_range[1] = 1;
+ quad_range[0] = 0;
+ quad_range[1] = band.quadrature_length;
+ CHK(LAY_SPECINT_QUADS_GET_BOUNDS
+ (&band, quad_range, x_h2o_range, bounds) == RES_BAD_ARG);
+
+ x_h2o_range[0] = 0;
+ x_h2o_range[1] = 0;
+ CHK(LAY_SPECINT_TAB_GET_BOUNDS(NULL, x_h2o_range, bounds) == RES_BAD_ARG);
+ CHK(LAY_SPECINT_TAB_GET_BOUNDS(&tab, NULL, bounds) == RES_BAD_ARG);
+ CHK(LAY_SPECINT_TAB_GET_BOUNDS(&tab, x_h2o_range, NULL) == RES_BAD_ARG);
+ x_h2o_range[0] = 1;
+ x_h2o_range[1] = 0;
+ CHK(LAY_SPECINT_TAB_GET_BOUNDS(&tab, x_h2o_range, bounds) == RES_BAD_ARG);
FOR_EACH(ilay, 0, nlays) {
- double lay_bounds[2] = {DBL_MAX, -DBL_MAX};
-
CHK(htgop_get_layer(htgop, ilay, &lay) == RES_OK);
FOR_EACH(iband, 0, nbands) {
- double band_bounds[2] = {DBL_MAX, -DBL_MAX};
- double band_bounds_tmp[2] = {DBL_MAX, -DBL_MAX};
- size_t iquad;
-
+ double quad_bounds2[2] = {DBL_MAX, -DBL_MAX};
CHK(GET_LAY_SPECINT(&lay, iband, &band) == RES_OK);
FOR_EACH(iquad, 0, band.quadrature_length) {
- struct LAY_SPECINT_TAB tab;
double quad_bounds[2] = {DBL_MAX, -DBL_MAX};
- double quad_bounds_tmp[2] = {DBL_MAX, -DBL_MAX};
- size_t itab;
+ double tab_bounds2[2] = {DBL_MAX, -DBL_MAX};
+ size_t itest;
CHK(GET_LAY_SPECINT_TAB(&band, iquad, &tab) == RES_OK);
- CHK(LAY_SPECINT_TAB_GET_BOUNDS(&tab, quad_bounds) == RES_OK);
-
FOR_EACH(itab, 0, tab.tab_length) {
const double k = GET_K(&tab, itab);
- quad_bounds_tmp[0] = MMIN(quad_bounds_tmp[0], k);
- quad_bounds_tmp[1] = MMAX(quad_bounds_tmp[1], k);
+ double tab_bounds[2] = {DBL_MAX, -DBL_MAX};
+ tab_bounds2[0] = MMIN(tab_bounds2[0], k);
+ tab_bounds2[1] = MMAX(tab_bounds2[1], k);
+
+ x_h2o_range[0] = tab.x_h2o_tab[0];
+ x_h2o_range[1] = tab.x_h2o_tab[itab];
+ CHK(LAY_SPECINT_TAB_GET_BOUNDS
+ (&tab, x_h2o_range, tab_bounds) == RES_OK);
+ CHK(eq_eps(tab_bounds2[0], tab_bounds[0], 1.e-6));
+ CHK(eq_eps(tab_bounds2[1], tab_bounds[1], 1.e-6));
}
- CHK(quad_bounds[0] == quad_bounds_tmp[0]);
- CHK(quad_bounds[1] == quad_bounds_tmp[1]);
- band_bounds_tmp[0] = MMIN(band_bounds_tmp[0], quad_bounds[0]);
- band_bounds_tmp[1] = MMAX(band_bounds_tmp[1], quad_bounds[1]);
- quad_range[0] = 0;
+ FOR_EACH(itest, 0, 10) {
+ const double bias = tab.x_h2o_tab[0];
+ const double scale = tab.x_h2o_tab[tab.tab_length-1] - bias;
+ const double r0 = rand_canonic();
+ const double r1 = rand_canonic();
+ double tab_bounds[2] = {DBL_MAX, -DBL_MAX};
+ double tab_bounds3[2] = {DBL_MAX, -DBL_MAX};
+ double k0, k1;
+ size_t tab_range[2];
+
+ x_h2o_range[0] = r0*scale + bias;
+ x_h2o_range[1] = r1*scale + bias;
+ if(x_h2o_range[0] > x_h2o_range[1]) {
+ SWAP(double, x_h2o_range[0], x_h2o_range[1]);
+ }
+ HTGOP(CONCAT(CONCAT(CONCAT(
+ layer_,DOMAIN),_spectral_interval_tab_fetch_),DATA)
+ (&tab, x_h2o_range[0], &k0));
+ HTGOP(CONCAT(CONCAT(CONCAT(
+ layer_,DOMAIN),_spectral_interval_tab_fetch_),DATA)
+ (&tab, x_h2o_range[1], &k1));
+
+ tab_bounds[0] = MMIN(k0, k1);
+ tab_bounds[1] = MMAX(k0, k1);
+
+ tab_range[0] = SIZE_MAX;
+ tab_range[1] = SIZE_MAX;
+ FOR_EACH(itab, 0, tab.tab_length) {
+ if(tab_range[0]==SIZE_MAX && tab.x_h2o_tab[itab] > x_h2o_range[0])
+ tab_range[0] = itab;
+ if(tab_range[1]==SIZE_MAX && tab.x_h2o_tab[itab] >= x_h2o_range[1])
+ tab_range[1] = itab;
+ if(tab_range[0] != SIZE_MAX && tab_range[1] != SIZE_MAX)
+ break;
+ }
+ CHK(tab_range[0] != SIZE_MAX && tab_range[1] != SIZE_MAX);
+
+ FOR_EACH(itab, tab_range[0], tab_range[1]) {
+ tab_bounds[0] = MMIN(GET_K(&tab, itab), tab_bounds[0]);
+ tab_bounds[1] = MMAX(GET_K(&tab, itab), tab_bounds[1]);
+ }
+
+ CHK(LAY_SPECINT_TAB_GET_BOUNDS
+ (&tab, x_h2o_range, tab_bounds3) == RES_OK);
+ CHK(eq_eps(tab_bounds[0], tab_bounds3[0], 1.e-6));
+ CHK(eq_eps(tab_bounds[1], tab_bounds3[1], 1.e-6));
+ }
+
+ quad_range[0] = iquad;
quad_range[1] = iquad;
+ x_h2o_range[0] = lay.x_h2o_tab[0];
+ x_h2o_range[1] = lay.x_h2o_tab[lay.tab_length-1];
CHK(LAY_SPECINT_QUADS_GET_BOUNDS
- (&band, quad_range, quad_bounds_tmp) == RES_OK);
+ (&band, quad_range, x_h2o_range, quad_bounds) == RES_OK);
+ CHK(eq_eps(quad_bounds[0], tab_bounds2[0], 1.e-6));
+ CHK(eq_eps(quad_bounds[1], tab_bounds2[1], 1.e-6));
- CHK(quad_bounds_tmp[0] == band_bounds_tmp[0]);
- CHK(quad_bounds_tmp[1] == band_bounds_tmp[1]);
- }
+ quad_bounds2[0] = MMIN(quad_bounds[0], quad_bounds2[0]);
+ quad_bounds2[1] = MMAX(quad_bounds[1], quad_bounds2[1]);
- band_range[0] = iband;
- band_range[1] = iband;
- lay_range[0] = ilay;
- lay_range[1] = ilay;
- CHK(LAY_GET_BOUNDS(htgop, lay_range, band_range, band_bounds) == RES_OK);
- CHK(band_bounds[0] == band_bounds_tmp[0]);
- CHK(band_bounds[1] == band_bounds_tmp[1]);
-
- lay_bounds[0] = MMIN(lay_bounds[0], band_bounds[0]);
- lay_bounds[1] = MMAX(lay_bounds[1], band_bounds[1]);
-
- band_range[0] = 0;
- band_range[1] = iband;
- CHK(LAY_GET_BOUNDS(htgop, lay_range, band_range, band_bounds) == RES_OK);
- CHK(lay_bounds[0] == band_bounds[0]);
- CHK(lay_bounds[1] == band_bounds[1]);
+ quad_range[0] = 0;
+ quad_range[1] = iquad;
+ CHK(LAY_SPECINT_QUADS_GET_BOUNDS
+ (&band, quad_range, x_h2o_range, quad_bounds) == RES_OK);
+ CHK(eq_eps(quad_bounds2[0], quad_bounds[0], 1.e-6));
+ CHK(eq_eps(quad_bounds2[1], quad_bounds[1], 1.e-6));
+ }
}
}
}
#undef GET_NSPECINTS
#undef GET_SPECINT
-#undef LAY_GET_BOUNDS
#undef LAY_SPECINT_QUADS_GET_BOUNDS
#undef LAY_SPECINT_TAB_GET_BOUNDS
#undef LAY_SPECINT
diff --git a/src/test_htgop_load.c b/src/test_htgop_load.c
@@ -352,6 +352,7 @@ main(int argc, char** argv)
CHK(lw_tab.tab_length == tab_len);
CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK);
CHK(lw_tab.ka_tab[itab] == dbl);
+ CHK(lw_tab.x_h2o_tab[itab] == lay.x_h2o_tab[itab]);
}
}
}
@@ -376,6 +377,7 @@ main(int argc, char** argv)
CHK(sw_tab.tab_length == tab_len);
CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK);
CHK(sw_tab.ka_tab[itab] == dbl);
+ CHK(sw_tab.x_h2o_tab[itab] == lay.x_h2o_tab[itab]);
}
}
}
@@ -400,6 +402,7 @@ main(int argc, char** argv)
CHK(sw_tab.tab_length == tab_len);
CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK);
CHK(sw_tab.ks_tab[itab] == dbl);
+ CHK(sw_tab.x_h2o_tab[itab] == lay.x_h2o_tab[itab]);
}
}
}
