commit a45589a28cbd3d45263cda6700b8ec3d9d1df99c
parent 7cf3fb9c94f52f07efed7808b0b4c12c3fd827ad
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 15 Apr 2020 15:55:00 +0200
Rm the functions explicitly relying on the CIE XYZ colorspace
Diffstat:
8 files changed, 176 insertions(+), 536 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -41,8 +41,7 @@ set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(HTGOP_FILES_SRC
- htgop.c
- htgop_sw_spectral_interval_CIE_XYZ.c)
+ htgop.c)
set(HTGOP_FILES_INC
htgop_c.h
htgop_fetch_radiative_properties.h
diff --git a/src/htgop.c b/src/htgop.c
@@ -27,11 +27,17 @@
#include <math.h>
/* Boundaries of the long wave spectral data */
-#define LW_WAVELENGTH_MIN 1000 /* In nanometer */
-#define LW_WAVELENGTH_MAX 100000 /* In nanometer */
+#define LW_WAVELENGTH_MIN 1000.0 /* In nanometer */
+#define LW_WAVELENGTH_MAX 100000.0 /* In nanometer */
#define LW_WAVENUMBER_MIN WLEN_TO_WNUM(LW_WAVELENGTH_MAX) /* In cm^-1 */
#define LW_WAVENUMBER_MAX WLEN_TO_WNUM(LW_WAVELENGTH_MIN) /* In cm^-1 */
+/* Boundaries of theshort wave spectral data */
+#define SW_WAVELENGTH_MIN 380.0 /* In nanometer */
+#define SW_WAVELENGTH_MAX 780.0 /* In nanometer */
+#define SW_WAVENUMBER_MIN WLEN_TO_WNUM(SW_WAVELENGTH_MAX) /* In cm^-1 */
+#define SW_WAVENUMBER_MAX WLEN_TO_WNUM(SW_WAVELENGTH_MIN) /* In cm^-1 */
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -322,10 +328,6 @@ 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_1931_XYZ_cdf(htgop);
- if(res != RES_OK) goto error;
-
exit:
reader_release(&rdr);
return res;
@@ -851,6 +853,46 @@ htgop_position_to_layer_id
}
res_T
+htgop_get_sw_spectral_intervals
+ (const struct htgop* htgop,
+ const double wnum_range[2],
+ size_t specint_range[2])
+{
+ const double* wnums = NULL;
+ size_t nwnums = 0;
+ res_T res = RES_OK;
+
+ if(!htgop || !wnum_range || wnum_range[0]>wnum_range[1] || !specint_range) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(wnum_range[0] < SW_WAVENUMBER_MIN
+ || wnum_range[1] > SW_WAVENUMBER_MAX) {
+ log_err(htgop,
+ "%s: the range [%g, %g] is not strictly included "
+ "in the long wave interval [%g, %g].\n",
+ FUNC_NAME,
+ wnum_range[0], wnum_range[1],
+ SW_WAVENUMBER_MIN, SW_WAVENUMBER_MAX);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ wnums = darray_double_cdata_get(&htgop->sw_specints.wave_numbers);
+ nwnums = darray_double_size_get(&htgop->sw_specints.wave_numbers);
+
+ res = get_spectral_intervals
+ (htgop, FUNC_NAME, wnum_range, wnums, nwnums, specint_range);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
htgop_get_lw_spectral_intervals
(const struct htgop* htgop,
const double wnum_range[2],
diff --git a/src/htgop.h b/src/htgop.h
@@ -289,28 +289,6 @@ htgop_layer_sw_spectral_interval_tab_fetch_kext
double* kext);
/*******************************************************************************
- * 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 */
-
-/*******************************************************************************
* Retrieve the boundaries of the radiative properties
******************************************************************************/
HTGOP_API res_T
@@ -381,11 +359,10 @@ htgop_spectral_interval_sample_quadrature
const double r, /* Canonical random number in [0, 1[ */
size_t* iquad_point); /* Id of the sample quadrature point */
-/* Return the inclusive range of spectral interval that overlaps the CIE XYZ
- * color space */
HTGOP_API res_T
-htgop_get_sw_spectral_intervals_CIE_XYZ
+htgop_get_sw_spectral_intervals
(const struct htgop* htgop,
+ const double wnum_range[2],
size_t specint_range[2]);
HTGOP_API res_T
diff --git a/src/htgop_c.h b/src/htgop_c.h
@@ -80,10 +80,6 @@ log_warn
#endif
;
-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,
diff --git a/src/htgop_sw_spectral_interval_CIE_XYZ.c b/src/htgop_sw_spectral_interval_CIE_XYZ.c
@@ -1,311 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (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>
-
-/* 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 do 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);
-}
-
-res_T
-htgop_get_sw_spectral_intervals_CIE_XYZ
- (const struct htgop* htgop,
- size_t specint_range[2])
-{
- const double* wnums = NULL;
- double wnum_range[2] = {0, 0};
- size_t nwnums = 0;
- res_T res = RES_OK;
-
- if(!htgop || !specint_range) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- wnum_range[0] = CIE_XYZ_WAVENUMBER_MIN;
- wnum_range[1] = CIE_XYZ_WAVENUMBER_MAX;
- wnums = darray_double_cdata_get(&htgop->sw_specints.wave_numbers);
- nwnums = darray_double_size_get(&htgop->sw_specints.wave_numbers);
-
- res = get_spectral_intervals
- (htgop, FUNC_NAME, wnum_range, wnums, nwnums, specint_range);
- if(res != RES_OK) goto error;
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-/*******************************************************************************
- * 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_check_specints.h b/src/test_htgop_check_specints.h
@@ -0,0 +1,116 @@
+/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (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"
+
+#if !defined(DOMAIN) \
+ || !defined(DOMAIN_WAVELENGTH_MIN) \
+ || !defined(DOMAIN_WAVELENGTH_MAX)
+ #error "Missing the <DATA|DOMAIN> macro."
+#endif
+
+/* Helper macros */
+#define GET_SPECINTS \
+ CONCAT(CONCAT(htgop_get_, DOMAIN), _spectral_intervals)
+#define GET_NSPECINTS \
+ CONCAT(CONCAT(htgop_get_, DOMAIN), _spectral_intervals_count)
+#define GET_SPECINT \
+ CONCAT(CONCAT(htgop_get_, DOMAIN), _spectral_interval)
+
+static void
+CONCAT(CONCAT(check_, DOMAIN), _specints)
+ (const struct htgop* htgop,
+ const double low, /* In wavenumber */
+ const double upp) /* In wavenumber */
+{
+ struct htgop_spectral_interval specint_low;
+ struct htgop_spectral_interval specint_upp;
+ double wnums[2];
+ size_t range[2];
+ size_t range2[2];
+ size_t i;
+ size_t nspecints;
+
+ CHK(low <= upp);
+
+ wnums[0] = low;
+ wnums[1] = upp;
+ CHK(GET_SPECINTS(NULL, wnums, range) == RES_BAD_ARG);
+ CHK(GET_SPECINTS(htgop, NULL, range) == RES_BAD_ARG);
+ CHK(GET_SPECINTS(htgop, wnums, NULL) == RES_BAD_ARG);
+ CHK(GET_SPECINTS(htgop, wnums, range) == RES_OK);
+ CHK(range[0] <= range[1]);
+
+ if(upp != low) {
+ wnums[0] = upp;
+ wnums[1] = low;
+ CHK(GET_SPECINTS(htgop, wnums, range) == RES_BAD_ARG);
+ }
+ wnums[0] = low;
+ wnums[1] = 1.e7/(DOMAIN_WAVELENGTH_MIN-1.0);
+ CHK(GET_SPECINTS(htgop, wnums, range) == RES_BAD_ARG);
+ wnums[0] = 1.e7/(DOMAIN_WAVELENGTH_MAX+1.0);
+ wnums[1] = upp;
+ CHK(GET_SPECINTS(htgop, wnums, range) == RES_BAD_ARG);
+ wnums[0] = low;
+ wnums[1] = upp;
+ CHK(GET_SPECINTS(htgop, wnums, range) == RES_OK);
+
+ CHK(GET_SPECINT(htgop, range[0], &specint_low) == RES_OK);
+ CHK(GET_SPECINT(htgop, range[1], &specint_upp) == RES_OK);
+ CHK(GET_NSPECINTS(htgop, &nspecints) == RES_OK);
+
+ CHK(specint_low.wave_numbers[0] < specint_upp.wave_numbers[1]);
+ CHK(specint_low.wave_numbers[0] <= wnums[0] || range[0] == 0);
+ CHK(specint_upp.wave_numbers[1] >= wnums[1] || range[1] == nspecints-1);
+
+ range2[0] = SIZE_MAX;
+ range2[1] = SIZE_MAX;
+ FOR_EACH(i, 0, nspecints) {
+ struct htgop_spectral_interval specint;
+ CHK(GET_SPECINT(htgop, i, &specint) == RES_OK);
+
+ if(specint.wave_numbers[0]<=wnums[0] && specint.wave_numbers[1]>wnums[0]) {
+ range2[0] = i;
+ }
+ if(specint.wave_numbers[0]<wnums[1] && specint.wave_numbers[1]>=wnums[1]) {
+ range2[1] = i;
+ }
+ }
+
+ CHK(range2[0] <= range2[1]);
+
+ if(range2[0] == SIZE_MAX) {
+ /* The loaded data does not strictly include the submitted range */
+ CHK(range[0] == 0);
+ } else {
+ CHK(range2[0] == range[0]);
+ }
+
+ if(range2[1] == SIZE_MAX) {
+ /* The loaded data does not strictly include the submitted range */
+ CHK(range[1] == nspecints-1);
+ } else {
+ CHK(range2[1] == range[1]);
+ }
+}
+
+#undef GET_SPECINT
+#undef GET_SPECINTS
+#undef GET_NSPECINTS
+#undef DOMAIN
+#undef DOMAIN_WAVELENGTH_MIN
+#undef DOMAIN_WAVELENGTH_MAX
+
diff --git a/src/test_htgop_load.c b/src/test_htgop_load.c
@@ -64,130 +64,15 @@ check_position_to_layer_id(const struct htgop* htgop)
}
}
-static void
-check_sw_specints_CIE_XYZ(const struct htgop* htgop)
-{
- struct htgop_spectral_interval specint_low;
- struct htgop_spectral_interval specint_upp;
- const double wnum_min_CIE_XYZ = 1.e7/780.0;
- const double wnum_max_CIE_XYZ = 1.e7/380.0;
- size_t range[2];
- size_t range2[2];
- size_t i;
- size_t nspecints;
-
- CHK(htgop_get_sw_spectral_intervals_CIE_XYZ(NULL, range) == RES_BAD_ARG);
- CHK(htgop_get_sw_spectral_intervals_CIE_XYZ(htgop, NULL) == RES_BAD_ARG);
- CHK(htgop_get_sw_spectral_intervals_CIE_XYZ(htgop, range) == RES_OK);
- CHK(range[0] <= range[1]);
-
- CHK(htgop_get_sw_spectral_interval(htgop, range[0], &specint_low) == RES_OK);
- CHK(htgop_get_sw_spectral_interval(htgop, range[1], &specint_upp) == RES_OK);
-
- CHK(specint_low.wave_numbers[0] < specint_upp.wave_numbers[1]);
- CHK(specint_low.wave_numbers[0] <= wnum_min_CIE_XYZ);
- CHK(specint_upp.wave_numbers[1] >= wnum_max_CIE_XYZ);
-
- CHK(htgop_get_sw_spectral_intervals_count(htgop, &nspecints) == RES_OK);
- range2[0] = SIZE_MAX;
- range2[1] = SIZE_MAX;
- FOR_EACH(i, 0, nspecints) {
- struct htgop_spectral_interval specint;
- CHK(htgop_get_sw_spectral_interval(htgop, i, &specint) == RES_OK);
-
- if(specint.wave_numbers[0] <= wnum_min_CIE_XYZ
- && specint.wave_numbers[1] > wnum_min_CIE_XYZ) {
- range2[0] = i;
- }
- if(specint.wave_numbers[0] < wnum_max_CIE_XYZ
- && specint.wave_numbers[1] >= wnum_max_CIE_XYZ) {
- range2[1] = i;
- }
- }
- CHK(range2[0] != SIZE_MAX);
- CHK(range2[1] != SIZE_MAX);
- CHK(range2[0] <= range2[1]);
- CHK(range2[0] == range[0]);
- CHK(range2[1] == range[1]);
-}
-
-static void
-check_lw_specints
- (const struct htgop* htgop,
- const double low, /* In wavenumber */
- const double upp) /* In wavenumber */
-{
- struct htgop_spectral_interval specint_low;
- struct htgop_spectral_interval specint_upp;
- double wnums[2];
- size_t range[2];
- size_t range2[2];
- size_t i;
- size_t nspecints;
-
- CHK(low <= upp);
-
- wnums[0] = low;
- wnums[1] = upp;
- CHK(htgop_get_lw_spectral_intervals(NULL, wnums, range) == RES_BAD_ARG);
- CHK(htgop_get_lw_spectral_intervals(htgop, NULL, range) == RES_BAD_ARG);
- CHK(htgop_get_lw_spectral_intervals(htgop, wnums, NULL) == RES_BAD_ARG);
- CHK(htgop_get_lw_spectral_intervals(htgop, wnums, range) == RES_OK);
- CHK(range[0] <= range[1]);
-
- if(upp != low) {
- wnums[0] = upp;
- wnums[1] = low;
- CHK(htgop_get_lw_spectral_intervals(htgop, wnums, range) == RES_BAD_ARG);
- }
- wnums[0] = low;
- wnums[1] = 1.e7/999;
- CHK(htgop_get_lw_spectral_intervals(htgop, wnums, range) == RES_BAD_ARG);
- wnums[0] = 1.e7/100001;
- wnums[1] = upp;
- CHK(htgop_get_lw_spectral_intervals(htgop, wnums, range) == RES_BAD_ARG);
- wnums[0] = low;
- wnums[1] = upp;
- CHK(htgop_get_lw_spectral_intervals(htgop, wnums, range) == RES_OK);
-
- CHK(htgop_get_lw_spectral_interval(htgop, range[0], &specint_low) == RES_OK);
- CHK(htgop_get_lw_spectral_interval(htgop, range[1], &specint_upp) == RES_OK);
- CHK(htgop_get_lw_spectral_intervals_count(htgop, &nspecints) == RES_OK);
-
- CHK(specint_low.wave_numbers[0] < specint_upp.wave_numbers[1]);
- CHK(specint_low.wave_numbers[0] <= wnums[0] || range[0] == 0);
- CHK(specint_upp.wave_numbers[1] >= wnums[1] || range[1] == nspecints-1);
-
- range2[0] = SIZE_MAX;
- range2[1] = SIZE_MAX;
- FOR_EACH(i, 0, nspecints) {
- struct htgop_spectral_interval specint;
- CHK(htgop_get_lw_spectral_interval(htgop, i, &specint) == RES_OK);
-
- if(specint.wave_numbers[0]<=wnums[0] && specint.wave_numbers[1]>wnums[0]) {
- range2[0] = i;
- }
- if(specint.wave_numbers[0]<wnums[1] && specint.wave_numbers[1]>=wnums[1]) {
- range2[1] = i;
- }
- }
-
- CHK(range2[0] <= range2[1]);
+#define DOMAIN sw
+#define DOMAIN_WAVELENGTH_MIN 100.0
+#define DOMAIN_WAVELENGTH_MAX 1200.0
+#include "test_htgop_check_specints.h"
- if(range2[0] == SIZE_MAX) {
- /* The loaded data does not strictly include the submitted range */
- CHK(range[0] == 0);
- } else {
- CHK(range2[0] == range[0]);
- }
-
- if(range2[1] == SIZE_MAX) {
- /* The loaded data does not strictly include the submitted range */
- CHK(range[1] == nspecints-1);
- } else {
- CHK(range2[1] == range[1]);
- }
-}
+#define DOMAIN lw
+#define DOMAIN_WAVELENGTH_MIN 1000.0
+#define DOMAIN_WAVELENGTH_MAX 100000.0
+#include "test_htgop_check_specints.h"
int
main(int argc, char** argv)
@@ -521,7 +406,7 @@ main(int argc, char** argv)
}
check_position_to_layer_id(htgop);
- check_sw_specints_CIE_XYZ(htgop);
+ check_sw_specints(htgop, 1.e7/780.0, 1.e7/380.0);
check_lw_specints(htgop, 1.e7/100000, 1.e7/1000);
check_lw_specints(htgop, 1.e7/5600, 1.e7/3000);
check_lw_specints(htgop, 1.e7/9500, 1.e7/9300);
diff --git a/src/test_htgop_sample.c b/src/test_htgop_sample.c
@@ -20,38 +20,6 @@
#include <string.h>
#define N 100000
-#define CIE_XYZ_WAVENUMBER_MIN (1.0e7/780.0)
-#define CIE_XYZ_WAVENUMBER_MAX (1.0e7/380.0)
-
-static void
-check_sw_sample_specint
- (struct htgop* htgop,
- res_T (*sample_func)(const struct htgop*, const double, size_t*),
- int* hist)
-{
- const double* wnums;
- size_t nspecints;
- size_t ispecint;
- size_t i;
-
- CHK(sample_func && htgop && hist);
- CHK(htgop_get_sw_spectral_intervals_count(htgop, &nspecints) == RES_OK);
- CHK(htgop_get_sw_spectral_intervals_wave_numbers(htgop, &wnums) == RES_OK);
-
- CHK(sample_func(NULL, 0, &ispecint) == RES_BAD_ARG);
- CHK(sample_func(htgop, 1, &ispecint) == RES_BAD_ARG);
- CHK(sample_func(htgop, 0, NULL) == RES_BAD_ARG);
-
- memset(hist, 0, sizeof(*hist)*nspecints);
- FOR_EACH(i, 0, N) {
- const double r = rand_canonic();
- CHK(sample_func(htgop, r, &ispecint) == RES_OK);
- CHK(ispecint < nspecints);
- CHK(wnums[ispecint+0] <= CIE_XYZ_WAVENUMBER_MAX);
- CHK(wnums[ispecint+1] >= CIE_XYZ_WAVENUMBER_MIN);
- hist[ispecint] += 1;
- }
-}
static void
check_sample_quadrature
@@ -110,11 +78,7 @@ main(int argc, char** argv)
struct mem_allocator allocator;
struct htgop* htgop;
const double* wnums;
- int* hist_X;
- int* hist_Y;
- int* hist_Z;
size_t nspecints;
- size_t ispecint;
if(argc < 2) {
fprintf(stderr, "Usage: %s FILENAME\n", argv[0]);
@@ -130,34 +94,6 @@ main(int argc, char** argv)
CHK(htgop_get_sw_spectral_intervals_wave_numbers(htgop, &wnums) == RES_OK);
- CHK(hist_X = MEM_CALLOC(&allocator, nspecints, sizeof(*hist_X)));
- 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_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
- || wnums[ispecint+1] < CIE_XYZ_WAVENUMBER_MIN) {
- CHK(hist_X[ispecint] == 0);
- CHK(hist_Y[ispecint] == 0);
- CHK(hist_Z[ispecint] == 0);
- } else {
- CHK(!hist_X[ispecint] || hist_X[ispecint] != hist_Y[ispecint]);
- CHK(!hist_X[ispecint] || hist_X[ispecint] != hist_Z[ispecint]);
- CHK(!hist_Y[ispecint] || hist_Y[ispecint] != hist_Z[ispecint]);
- }
- }
-
- MEM_RM(&allocator, hist_X);
- MEM_RM(&allocator, hist_Y);
- MEM_RM(&allocator, hist_Z);
-
check_sample_quadrature(htgop, htgop_get_sw_spectral_intervals_count,
htgop_get_sw_spectral_interval);
check_sample_quadrature(htgop, htgop_get_lw_spectral_intervals_count,
