commit ac5a4071d4ab1e6bbf07ef5144d7ca4875708980
parent 8c91240f90f27048e16f570cc0b44991ee2d8c95
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 30 Jul 2018 15:46:15 +0200
Implement and test htgop_position_to_layer_id
Diffstat:
4 files changed, 112 insertions(+), 2 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -24,7 +24,7 @@ option(NO_TEST "Do not build tests" OFF)
# Check dependencies
################################################################################
find_package(RCMake 0.3 REQUIRED)
-find_package(RSys 0.6 REQUIRED)
+find_package(RSys 0.7 REQUIRED)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
include(rcmake)
diff --git a/src/htgop.c b/src/htgop.c
@@ -60,6 +60,14 @@ cmp_dbl(const void* a, const void* b)
return d0 < d1 ? -1 : (d0 > d1 ? 1 : 0);
}
+static INLINE int
+cmp_lvl(const void* key, const void* data)
+{
+ const double height = *((const double*)key);
+ const struct htgop_level* lvl = (const struct htgop_level*)data;
+ return height < lvl->height ? -1 : (height > lvl->height ? 1 : 0);
+}
+
static FINLINE double
trapezoidal_integration
(const double lambda_lo, /* Integral lower bound. In nanometer */
@@ -411,6 +419,7 @@ load_stream(struct htgop* htgop, FILE* stream, const char* stream_name)
struct layer* layers = NULL;
unsigned long nlvls, nlays, tab_len, nspecints_lw, nspecints_sw;
unsigned long ilvl, ilay, itab;
+ unsigned long iline_saved;
res_T res = RES_OK;
ASSERT(htgop && stream && stream_name);
reader_init(&rdr, stream, stream_name);
@@ -424,8 +433,21 @@ load_stream(struct htgop* htgop, FILE* stream, const char* stream_name)
/* Parse the number of levels/layers */
CALL(cstr_to_ulong(read_line(&rdr), &nlvls));
+ iline_saved = rdr.iline;
CALL(cstr_to_ulong(read_line(&rdr), &nlays));
- if(nlvls != nlays + 1) { res = RES_BAD_ARG; goto error; }
+ if(!nlays) {
+ log_err(htgop, "%s:%lu: The number of layers cannot be null.\n",
+ rdr.name, rdr.iline);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ if(nlvls != nlays + 1) {
+ log_err(htgop,
+ "%s:%lu: Invalid number of levels `%lu' (#layers = `%lu').\n",
+ rdr.name, iline_saved, nlvls, nlays);
+ res = RES_BAD_ARG;
+ goto error;
+ }
/* Parse the ground temperature */
CALL(cstr_to_double(read_line(&rdr), &htgop->ground.temperature));
@@ -935,6 +957,42 @@ htgop_sample_sw_spectral_interval_CIE_Z
r, ispecint);
}
+res_T
+htgop_position_to_layer_id
+ (const struct htgop* htgop, const double pos, size_t* ilayer)
+{
+ const struct htgop_level* lvls;
+ size_t nlvls;
+
+ if(!htgop || !ilayer) return RES_BAD_ARG;
+
+ lvls = darray_level_cdata_get(&htgop->levels);
+ nlvls = darray_level_size_get(&htgop->levels);
+ ASSERT(nlvls);
+
+ if(pos == lvls[0].height) {
+ *ilayer = 0;
+ } else if(pos == lvls[nlvls-1].height) {
+ *ilayer = nlvls - 2;
+ } else {
+ const struct htgop_level* find;
+ size_t i;
+
+ find = search_lower_bound(&pos, lvls, nlvls, sizeof(*lvls), cmp_lvl);
+ if(!find || find == lvls) {
+ log_err(htgop, "%s: the position `%g' is outside the atmospheric slices.\n",
+ FUNC_NAME, pos);
+ return RES_BAD_ARG;
+ }
+
+ i = (size_t)(find - lvls);
+ ASSERT(i < nlvls && i > 0);
+ ASSERT(lvls[i].height > pos && (!i || lvls[i-1].height <= pos));
+ *ilayer = i - 1;
+ }
+ return RES_OK;
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
diff --git a/src/htgop.h b/src/htgop.h
@@ -233,6 +233,12 @@ htgop_layer_sw_spectral_interval_get_tab
const size_t iquadrature_point,
struct htgop_layer_sw_spectral_interval_tab* tab);
+HTGOP_API res_T
+htgop_position_to_layer_id
+ (const struct htgop* htgop,
+ const double pos,
+ size_t* ilayer);
+
END_DECLS
#endif /* HTGOP */
diff --git a/src/test_htgop_load.c b/src/test_htgop_load.c
@@ -13,6 +13,8 @@
* 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_reader.h"
@@ -20,6 +22,48 @@
#include <rsys/cstr.h>
+static void
+check_position_to_layer_id(const struct htgop* htgop)
+{
+ struct htgop_level lvl0, lvl1;
+ const size_t N = 10000;
+ double org, sz;
+ size_t ilay;
+ size_t i, n;
+
+ CHK(htgop_get_levels_count(htgop, &n) == RES_OK);
+ CHK(n > 1);
+
+ /* Test the htgop_position_to_layer_id function */
+ CHK(htgop_get_level(htgop, 0, &lvl0) == RES_OK);
+ CHK(htgop_get_level(htgop, n-1, &lvl1) == RES_OK);
+ org = lvl0.height;
+ sz = lvl1.height - lvl0.height;
+
+ CHK(htgop_position_to_layer_id(NULL, lvl0.height, &ilay) == RES_BAD_ARG);
+ CHK(htgop_position_to_layer_id
+ (htgop, nextafter(lvl0.height,-DBL_MAX), &ilay) == RES_BAD_ARG);
+ CHK(htgop_position_to_layer_id
+ (htgop, nextafter(lvl1.height, DBL_MAX), &ilay) == RES_BAD_ARG);
+ CHK(htgop_position_to_layer_id(htgop, lvl0.height, NULL) == RES_BAD_ARG);
+
+ CHK(htgop_position_to_layer_id(htgop, lvl0.height, &ilay) == RES_OK);
+ CHK(ilay == 0);
+ CHK(htgop_position_to_layer_id(htgop, lvl1.height, &ilay) == RES_OK);
+ CHK(ilay == n - 2);
+
+ FOR_EACH(i, 0, N) {
+ const double r = rand_canonic();
+ const double pos = org + r * sz;
+
+ CHK(htgop_position_to_layer_id(htgop, pos, &ilay) == RES_OK);
+ CHK(htgop_get_level(htgop, ilay+0, &lvl0) == RES_OK);
+ CHK(htgop_get_level(htgop, ilay+1, &lvl1) == RES_OK);
+ CHK(pos >= lvl0.height);
+ CHK(pos <= lvl1.height);
+ }
+}
+
int
main(int argc, char** argv)
{
@@ -310,6 +354,8 @@ main(int argc, char** argv)
}
}
+ check_position_to_layer_id(htgop);
+
CHK(fclose(fp) == 0);
reader_release(&rdr);
CHK(htgop_ref_put(htgop) == RES_OK);
