commit f2b34e491acf232348062eb57fe6e71063963630
parent d28e18ce65bb4dee17c87953acaf4dbbdb0ab298
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 12 Mar 2020 14:24:33 +0100
Add and test the htgop_get_lw_intervals function
Diffstat:
5 files changed, 216 insertions(+), 37 deletions(-)
diff --git a/src/htgop.c b/src/htgop.c
@@ -26,6 +26,12 @@
#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_WAVENUMBER_MIN WLEN_TO_WNUM(LW_WAVELENGTH_MAX) /* In cm^-1 */
+#define LW_WAVENUMBER_MAX WLEN_TO_WNUM(LW_WAVELENGTH_MIN) /* In cm^-1 */
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -779,6 +785,46 @@ htgop_position_to_layer_id
return RES_OK;
}
+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 +892,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
@@ -372,7 +372,13 @@ htgop_spectral_interval_sample_quadrature
HTGOP_API res_T
htgop_get_sw_spectral_intervals_CIE_XYZ
(const struct htgop* htgop,
- size_t specint_raneg[2]);
+ size_t specint_range[2]);
+
+HTGOP_API res_T
+htgop_get_lw_spectral_intervals
+ (const struct htgop* htgop,
+ const double wnum_range[2],
+ size_t specint_range[2]);
END_DECLS
diff --git a/src/htgop_c.h b/src/htgop_c.h
@@ -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
@@ -104,5 +108,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_sw_spectral_interval_CIE_XYZ.c b/src/htgop_sw_spectral_interval_CIE_XYZ.c
@@ -20,10 +20,6 @@
#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 */
@@ -167,48 +163,27 @@ htgop_get_sw_spectral_intervals_CIE_XYZ
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;
+ double wnum_range[2] = {0, 0};
+ size_t nwnums = 0;
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;
+ if(!htgop || !specint_range) {
+ res = RES_BAD_ARG;
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;
+ 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);
- 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);
+ res = get_spectral_intervals
+ (htgop, FUNC_NAME, wnum_range, wnums, nwnums, specint_range);
+ if(res != RES_OK) goto error;
exit:
return res;
error:
- if(specint_range) {
- specint_range[0] = SIZE_MAX;
- specint_range[1] = 0;
- }
goto exit;
}
diff --git a/src/test_htgop_load.c b/src/test_htgop_load.c
@@ -111,6 +111,84 @@ check_sw_specints_CIE_XYZ(const struct htgop* htgop)
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]);
+
+ 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]);
+ }
+}
+
int
main(int argc, char** argv)
{
@@ -410,6 +488,10 @@ main(int argc, char** argv)
check_position_to_layer_id(htgop);
check_sw_specints_CIE_XYZ(htgop);
+ 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);
