htgop

Optical properties of a gas mixture
git clone git://git.meso-star.fr/htgop.git
Log | Files | Refs | README | LICENSE

commit 3cbb00a2bd7ad12bbb0fcf04c9d4ba9f1bfb52ef
parent e6d9162c3b790252d5ad46cfa61ca634abcbf43b
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Mon,  4 May 2020 10:51:01 +0200

Merge branch 'release_0.1'

Diffstat:
MREADME.md | 18+++++++++++++++---
Mcmake/CMakeLists.txt | 9++++-----
Msrc/htgop.c | 218+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--
Msrc/htgop.h | 51++++++++++++++++++++++++---------------------------
Msrc/htgop_c.h | 21++++++++++++++++-----
Msrc/htgop_dbllst.h | 2+-
Msrc/htgop_fetch_radiative_properties.h | 2+-
Msrc/htgop_get_radiative_properties_bounds.h | 2+-
Msrc/htgop_layer.h | 2+-
Msrc/htgop_parse_layers_spectral_intervals_data.h | 2+-
Msrc/htgop_reader.h | 2+-
Msrc/htgop_spectral_intervals.h | 2+-
Dsrc/htgop_sw_spectral_interval_CIE_XYZ.c | 336-------------------------------------------------------------------------------
Msrc/test_htgop.c | 2+-
Asrc/test_htgop_check_specints.h | 116+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/test_htgop_fetch_radiative_properties.c | 2+-
Msrc/test_htgop_fetch_radiative_properties.h | 2+-
Msrc/test_htgop_get_radiative_properties_bounds.c | 2+-
Msrc/test_htgop_get_radiative_properties_bounds.h | 2+-
Msrc/test_htgop_load.c | 97++++++++++++++++++++++++++++++++++++++++---------------------------------------
Msrc/test_htgop_sample.c | 66+-----------------------------------------------------------------
Msrc/test_htgop_utils.h | 2+-
22 files changed, 451 insertions(+), 507 deletions(-)

diff --git a/README.md b/README.md @@ -20,13 +20,25 @@ informations on CMake. ## Release notes +### Version 0.1 + +- Add the `htgop_get_<lw|sw>_spectral_intervals` functions: they return the + indices of the lower and upper spectral intervals that include a given range + of long/short waves. +- Add the `htgop_find_<lw|sw>_spectral_interval_id` functions: they return the + index of the spectral interval that includes the submitted short/long wave. +- Remove the functions explicitly relying onto the CIE 1931 XYZ color space, + i.e. `htgop_sample_sw_spectral_interval_CIE_1931_<X|Y|Z>` and + `htgop_get_sw_spectral_intervals_CIE_XYZ`. + ### Version 0.0.2 - Fix an issue when the parsed line is greater than 128 characters. ## Licenses -Copyright (C) 2018, 2019 |Meso|Star>. htgop is free software released under the -GPL v3+ license: GNU GPL version 3 or later. You are welcome to redistribute it -under certain conditions; refer to the COPYING file for details. +Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com). htgop is +free software released under the GPL v3+ license: GNU GPL version 3 or later. +You are welcome to redistribute it under certain conditions; refer to the +COPYING file for details. diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -1,4 +1,4 @@ -# Copyright (C) 2018, 2019 |Meso|Star> +# 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 @@ -36,13 +36,12 @@ include_directories(${RSys_INCLUDE_DIR}) # Configure and define targets ################################################################################ set(VERSION_MAJOR 0) -set(VERSION_MINOR 0) -set(VERSION_PATCH 2) +set(VERSION_MINOR 1) +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 @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 @@ -26,6 +26,18 @@ #include <math.h> +/* Boundaries of the long wave spectral data */ +#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 100.0 /* In nanometer */ +#define SW_WAVELENGTH_MAX 1200.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 ******************************************************************************/ @@ -316,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; @@ -713,6 +721,71 @@ htgop_get_sw_spectral_interval } res_T +htgop_find_lw_spectral_interval_id + (const struct htgop* htgop, + const double wnum, + size_t* ispecint) +{ + const double* find = NULL; + const double* wnums = NULL; + size_t nwnums = 0; + if(!htgop || !ispecint) return RES_BAD_ARG; + + wnums = darray_double_cdata_get(&htgop->lw_specints.wave_numbers); + nwnums = darray_double_size_get(&htgop->lw_specints.wave_numbers); + + find = search_lower_bound(&wnum, wnums, nwnums, sizeof(*wnums), cmp_dbl); + if(!find || find == wnums) { + log_warn(htgop, "%s: the wavenumber `%g' is not included in any band.\n", + FUNC_NAME, wnum); + *ispecint = SIZE_MAX; + } else { + *ispecint = (size_t)(find - wnums - 1); +#ifndef NDEBUG + { + size_t n; + HTGOP(get_lw_spectral_intervals_count(htgop, &n)); + ASSERT(*ispecint < n); + } +#endif + } + return RES_OK; +} + +res_T +htgop_find_sw_spectral_interval_id + (const struct htgop* htgop, + const double wnum, + size_t* ispecint) +{ + const double* find = NULL; + const double* wnums = NULL; + size_t nwnums = 0; + if(!htgop || !ispecint) return RES_BAD_ARG; + + wnums = darray_double_cdata_get(&htgop->sw_specints.wave_numbers); + nwnums = darray_double_size_get(&htgop->sw_specints.wave_numbers); + + find = search_lower_bound(&wnum, wnums, nwnums, sizeof(*wnums), cmp_dbl); + if(!find || find == wnums) { + log_warn(htgop, "%s: the wavenumber `%g' is not included in any band.\n", + FUNC_NAME, wnum); + *ispecint = SIZE_MAX; + } else { + *ispecint = (size_t)(find - wnums - 1); +#ifndef NDEBUG + { + size_t n; + HTGOP(get_sw_spectral_intervals_count(htgop, &n)); + ASSERT(*ispecint < n); + } +#endif + + } + return RES_OK; +} + +res_T htgop_spectral_interval_sample_quadrature (const struct htgop_spectral_interval* specint, const double r, /* Canonical number in [0, 1[ */ @@ -779,6 +852,86 @@ htgop_position_to_layer_id return RES_OK; } +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], + 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] < LW_WAVENUMBER_MIN + || wnum_range[1] > LW_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], + LW_WAVENUMBER_MIN, LW_WAVENUMBER_MAX); + res = RES_BAD_ARG; + goto error; + } + + wnums = darray_double_cdata_get(&htgop->lw_specints.wave_numbers); + nwnums = darray_double_size_get(&htgop->lw_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; +} + /* Generate the htgop_layer_fetch_lw_ka function */ #define DOMAIN lw #define DATA ka @@ -846,3 +999,58 @@ log_warn(const struct htgop* htgop, const char* msg, ...) va_end(vargs_list); } +res_T +get_spectral_intervals + (const struct htgop* htgop, + const char* func_name, + const double wnum_range[2], + const double* wnums, + const size_t nwnums, + size_t specint_range[2]) +{ + double wnum_min = 0; + double wnum_max = 0; + double* low = NULL; + double* upp = NULL; + res_T res = RES_OK; + ASSERT(htgop && wnum_range && wnums && nwnums && specint_range); + ASSERT(wnum_range[0] <= wnum_range[1]); + + wnum_min = nextafter(wnum_range[0], DBL_MAX); + wnum_max = wnum_range[1]; + low = search_lower_bound(&wnum_min, wnums, nwnums, sizeof(double), cmp_dbl); + upp = search_lower_bound(&wnum_max, wnums, nwnums, sizeof(double), cmp_dbl); + + if(!low || upp == wnums) { + log_err(htgop, + "%s: the loaded data do not overlap the wave numbers in [%g, %g].\n", + func_name, wnum_range[0], wnum_range[1]); + res = RES_BAD_OP; + goto error; + } + ASSERT(low[0] >= wnum_min && (low == wnums || low[-1] < wnum_min)); + ASSERT(!upp || (upp[0] >= wnum_max && upp[-1] < wnum_max)); + + specint_range[0] = low == wnums ? 0 : (size_t)(low - wnums) - 1; + specint_range[1] = !upp ? nwnums - 1 : (size_t)(upp - wnums); + /* Transform the upper wnum bound in spectral interval id */ + specint_range[1] = specint_range[1] - 1; + + ASSERT(specint_range[0] <= specint_range[1]); + ASSERT(wnums[specint_range[0]+1] > wnum_range[0]); + ASSERT(wnums[specint_range[1]+0] < wnum_range[1]); + + ASSERT(wnums[specint_range[0]+0] <= wnum_range[0] + || specint_range[0]==0); /* The data do not include `wnum_range' */ + ASSERT(wnums[specint_range[1]+1] >= wnum_range[1] + || specint_range[1] + 1 == nwnums-1);/* The data do not include `wnum_range' */ + +exit: + return res; +error: + if(specint_range) { + specint_range[0] = SIZE_MAX; + specint_range[1] = 0; + } + goto exit; +} diff --git a/src/htgop.h b/src/htgop.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 @@ -202,6 +202,20 @@ htgop_get_sw_spectral_interval const size_t ispectral_interval, struct htgop_spectral_interval* interval); +/* Find the id of the long wave interval that includes `wave_number'. */ +HTGOP_API res_T +htgop_find_lw_spectral_interval_id + (const struct htgop* htgop, + const double wave_number, /* In cm^-1 */ + size_t* ispectral_interval); /* SIZE_MAX <=> not found */ + +/* Find the id of the short wave interval that includes `wave_number' */ +HTGOP_API res_T +htgop_find_sw_spectral_interval_id + (const struct htgop* htgop, + const double wave_number, /* In cm^-1 */ + size_t* ispectral_interval); /* SIZE_MAX <=> not found */ + HTGOP_API res_T htgop_layer_get_lw_spectral_interval (const struct htgop_layer* layer, @@ -275,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 @@ -367,12 +359,17 @@ 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 +htgop_get_lw_spectral_intervals (const struct htgop* htgop, - size_t specint_raneg[2]); + const double wnum_range[2], + size_t specint_range[2]); END_DECLS diff --git a/src/htgop_c.h b/src/htgop_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 @@ -22,6 +22,10 @@ #include <rsys/dynamic_array.h> #include <rsys/ref_count.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) + /* Generate the dynamic array of level */ #define DARRAY_NAME level #define DARRAY_DATA struct htgop_level @@ -76,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, @@ -104,5 +104,16 @@ layer_sw_spectral_interval_tab_interpolate_kext const double x_h2o, const double* x_h2o_upp); /* Pointer toward the 1st tabbed xH2O >= x_h2o */ +/* Return the *inclusive* lower/upper index of the spectral intervals that + * overlaps the submitted wave number range */ +extern LOCAL_SYM res_T +get_spectral_intervals + (const struct htgop* htgop, + const char* func_name, + const double wnum_range[2], + const double* wnums, + const size_t nwnums, + size_t specint_range[2]); + #endif /* HTGOP_C_H */ diff --git a/src/htgop_dbllst.h b/src/htgop_dbllst.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 diff --git a/src/htgop_fetch_radiative_properties.h b/src/htgop_fetch_radiative_properties.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 diff --git a/src/htgop_get_radiative_properties_bounds.h b/src/htgop_get_radiative_properties_bounds.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 diff --git a/src/htgop_layer.h b/src/htgop_layer.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 diff --git a/src/htgop_parse_layers_spectral_intervals_data.h b/src/htgop_parse_layers_spectral_intervals_data.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 diff --git a/src/htgop_reader.h b/src/htgop_reader.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 diff --git a/src/htgop_spectral_intervals.h b/src/htgop_spectral_intervals.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 diff --git a/src/htgop_sw_spectral_interval_CIE_XYZ.c b/src/htgop_sw_spectral_interval_CIE_XYZ.c @@ -1,336 +0,0 @@ -/* Copyright (C) 2018, 2019 |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> - -/* 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 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; - const double wnum_min = nextafter(CIE_XYZ_WAVENUMBER_MIN, DBL_MAX); - const double wnum_max = CIE_XYZ_WAVENUMBER_MAX; - size_t nwnums; - res_T res = RES_OK; - double* low = NULL; - double* upp = NULL; - if(!htgop || !specint_range) return RES_BAD_ARG; - - wnums = darray_double_cdata_get(&htgop->sw_specints.wave_numbers); - nwnums = darray_double_size_get(&htgop->sw_specints.wave_numbers); - - low = search_lower_bound(&wnum_min, wnums, nwnums, sizeof(double), cmp_dbl); - upp = search_lower_bound(&wnum_max, wnums, nwnums, sizeof(double), cmp_dbl); - - if(!low || upp == wnums) { - log_err(htgop, - "%s: the loaded data do not overlap the CIE XYZ color space.\n", - FUNC_NAME); - res = RES_BAD_OP; - goto error; - } - ASSERT(low[0] >= wnum_min && (low == wnums || low[-1] < wnum_min)); - ASSERT(!upp || (upp[0] >= wnum_max && upp[-1] < wnum_max)); - - specint_range[0] = low == wnums ? 0 : (size_t)(low - wnums) - 1; - specint_range[1] = !upp ? nwnums - 1 : (size_t)(upp - wnums); - /* Transform the upper wnum bound in spectral interval id */ - specint_range[1] = specint_range[1] - 1; - - ASSERT(specint_range[0] < specint_range[1]); - ASSERT(wnums[specint_range[0]+0] <= CIE_XYZ_WAVENUMBER_MIN); - ASSERT(wnums[specint_range[0]+1] > CIE_XYZ_WAVENUMBER_MIN); - ASSERT(wnums[specint_range[1]+0] < CIE_XYZ_WAVENUMBER_MAX); - ASSERT(wnums[specint_range[1]+1] >= CIE_XYZ_WAVENUMBER_MAX); - -exit: - return res; -error: - if(specint_range) { - specint_range[0] = SIZE_MAX; - specint_range[1] = 0; - } - 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.c b/src/test_htgop.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 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_fetch_radiative_properties.c b/src/test_htgop_fetch_radiative_properties.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 diff --git a/src/test_htgop_fetch_radiative_properties.h b/src/test_htgop_fetch_radiative_properties.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 diff --git a/src/test_htgop_get_radiative_properties_bounds.c b/src/test_htgop_get_radiative_properties_bounds.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 diff --git a/src/test_htgop_get_radiative_properties_bounds.h b/src/test_htgop_get_radiative_properties_bounds.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 diff --git a/src/test_htgop_load.c b/src/test_htgop_load.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 @@ -64,52 +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]); -} +#define DOMAIN sw +#define DOMAIN_WAVELENGTH_MIN 100.0 +#define DOMAIN_WAVELENGTH_MAX 1200.0 +#include "test_htgop_check_specints.h" + +#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) @@ -253,6 +216,18 @@ main(int argc, char** argv) CHK(cstr_to_ulong(read_line(&rdr), &ul) == RES_OK); CHK(ul == nspecints_sw); + /* Search for a LW spectral interval */ + CHK(htgop_find_lw_spectral_interval_id(NULL, 0, &ispecint) == RES_BAD_ARG); + CHK(htgop_find_lw_spectral_interval_id(htgop, 0, NULL) == RES_BAD_ARG); + CHK(htgop_find_lw_spectral_interval_id(htgop, 0, &ispecint) == RES_OK); + CHK(ispecint == SIZE_MAX); + + /* Search for a SW spectral interval */ + CHK(htgop_find_sw_spectral_interval_id(NULL, 0, &ispecint) == RES_BAD_ARG); + CHK(htgop_find_sw_spectral_interval_id(htgop, 0, NULL) == RES_BAD_ARG); + CHK(htgop_find_sw_spectral_interval_id(htgop, 0, &ispecint) == RES_OK); + CHK(ispecint == SIZE_MAX); + /* Per LW spectral interval data */ CHK(htgop_get_lw_spectral_intervals_wave_numbers(NULL, NULL) == RES_BAD_ARG); CHK(htgop_get_lw_spectral_intervals_wave_numbers(htgop, NULL) == RES_BAD_ARG); @@ -266,10 +241,21 @@ main(int argc, char** argv) CHK(htgop_get_lw_spectral_interval(htgop, nspecints_lw, &specint) == RES_BAD_ARG); CHK(htgop_get_lw_spectral_interval(NULL, 0, &specint) == RES_BAD_ARG); FOR_EACH(ispecint, 0, nspecints_lw) { + size_t i; CHK(htgop_get_lw_spectral_interval(htgop, ispecint, &specint) == RES_OK); CHK(specint.wave_numbers[0] == wnums[ispecint+0]); CHK(specint.wave_numbers[1] == wnums[ispecint+1]); + /* Check the finding of a LW spectral interval */ + FOR_EACH(i, 0, 10) { + double wnum; + size_t id; + wnum = wnums[ispecint+0] + + rand_canonic() * (wnums[ispecint+1] - wnums[ispecint+0]); + CHK(htgop_find_lw_spectral_interval_id(htgop, wnum, &id) == RES_OK); + CHK(id == ispecint); + } + CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); CHK(specint.wave_numbers[0] == dbl); CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); @@ -304,10 +290,21 @@ main(int argc, char** argv) CHK(htgop_get_sw_spectral_interval(htgop, nspecints_sw, &specint) == RES_BAD_ARG); CHK(htgop_get_sw_spectral_interval(NULL, 0, &specint) == RES_BAD_ARG); FOR_EACH(ispecint, 0, nspecints_sw) { + int i; CHK(htgop_get_sw_spectral_interval(htgop, ispecint, &specint) == RES_OK); CHK(specint.wave_numbers[0] == wnums[ispecint+0]); CHK(specint.wave_numbers[1] == wnums[ispecint+1]); + /* Check the finding of a SW spectral interval */ + FOR_EACH(i, 0, 10) { + double wnum; + size_t id; + wnum = wnums[ispecint+0] + + rand_canonic() * (wnums[ispecint+1] - wnums[ispecint+0]); + CHK(htgop_find_sw_spectral_interval_id(htgop, wnum, &id) == RES_OK); + CHK(id == ispecint); + } + CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); CHK(specint.wave_numbers[0] == dbl); CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); @@ -409,7 +406,11 @@ 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); + check_lw_specints(htgop, 1.e7/5000, 1.e7/5000); CHK(fclose(fp) == 0); reader_release(&rdr); diff --git a/src/test_htgop_sample.c b/src/test_htgop_sample.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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 @@ -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, diff --git a/src/test_htgop_utils.h b/src/test_htgop_utils.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018, 2019 |Meso|Star> (contact@meso-star.com) +/* 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