commit 7335d853d797ee252a3c432aceadc95a880dd8a7
parent d375ca7222a77bcf9b34f7e3523f1238a923de44
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 3 Jul 2018 10:21:09 +0200
Pursue the loading of the mie data
Diffstat:
2 files changed, 199 insertions(+), 25 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -39,6 +39,10 @@ include_directories(
${NETCDF_C_INCLUDE_DIRS}
${RSys_INCLUDE_DIR})
+if(CMAKE_COMPILER_IS_GNUCC)
+ set(MATH_LIB m)
+endif()
+
################################################################################
# Configure and define targets
################################################################################
@@ -59,7 +63,7 @@ rcmake_prepend_path(HTMIE_FILES_INC_API ${HTMIE_SOURCE_DIR})
rcmake_prepend_path(HTMIE_FILES_DOC ${PROJECT_SOURCE_DIR}/../../)
add_library(htmie SHARED ${HTMIE_FILES_SRC} ${HTMIE_FILES_INC} ${HTMIE_FILES_INC_API})
-target_link_libraries(htmie RSys ${NETCDF_C_LIBRARIES})
+target_link_libraries(htmie RSys ${NETCDF_C_LIBRARIES} ${MATH_LIB})
set_target_properties(htmie PROPERTIES
DEFINE_SYMBOL HTMIE_SHARED_BUILD
diff --git a/src/htmie.c b/src/htmie.c
@@ -1,6 +1,6 @@
/* Copyright (C) 2018 |Meso|Star> (contact@meso-star.com)
*
-* This program is free software: you can redistribute it and/or modify
+ * 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.
@@ -35,10 +35,12 @@
struct htmie {
struct darray_double wavelengths; /* In nano-meters */
- struct darray_double xsections_absorption;
- struct darray_double xsections_scattering;
+ struct darray_double xsections_absorption; /* In m^2.particle^-1 */
+ struct darray_double xsections_scattering; /* In m^2.particle^-1 */
+ struct darray_double g; /* Asymetric parameter */
- double rmod; /* Modal mean radius in micro-meters */
+ double rmod; /* Modal mean radius of droplet in micro-meters */
+ double smod; /* Modal std dev of droplet size in micro-meters */
int verbose; /* Verbosity level */
struct logger* logger;
@@ -161,7 +163,7 @@ sizeof_nctype(const nc_type type)
}
static res_T
-load_wavelengths(struct htmie* htmie, const int nc)
+load_wavelengths(struct htmie* htmie, const int nc, const double range[2])
{
size_t len;
size_t start;
@@ -172,7 +174,7 @@ load_wavelengths(struct htmie* htmie, const int nc)
int err= NC_NOERR;
int type;
res_T res = RES_OK;
- ASSERT(htmie);
+ ASSERT(htmie && nc != INVALID_NC_ID && range && range[0] <= range[1]);
err = nc_inq_varid(nc, "lambda", &id);
if(err != NC_NOERR) {
@@ -195,7 +197,7 @@ load_wavelengths(struct htmie* htmie, const int nc)
if(type != NC_DOUBLE && type != NC_FLOAT) {
log_err(htmie,
- "Type type of the 'lambda' variable cannot be %s. "
+ "The type of the 'lambda' variable cannot be %s. "
"Expecting floating point data.\n",
nctype_to_str(type));
res = RES_BAD_ARG;
@@ -212,6 +214,18 @@ load_wavelengths(struct htmie* htmie, const int nc)
NC(get_vara_double(nc, id, &start, &len,
darray_double_data_get(&htmie->wavelengths)));
+ /* Check data validity */
+ FOR_EACH(i, 0, len) {
+ if(darray_double_cdata_get(&htmie->wavelengths)[i] < range[0]
+ || darray_double_cdata_get(&htmie->wavelengths)[i] > range[1]) {
+ log_err(htmie, "Invalid wavelength %g. It must be in [%g, %g]\n",
+ darray_double_cdata_get(&htmie->wavelengths)[i],
+ range[0], range[1]);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+
FOR_EACH(i, 1, len) {
if(darray_double_cdata_get(&htmie->wavelengths)[i-1]
> darray_double_cdata_get(&htmie->wavelengths)[i-0]) {
@@ -236,40 +250,46 @@ error:
}
static res_T
-load_modal_mean_radius(struct htmie* htmie, const int nc)
+load_droplet_distribution_variable
+ (struct htmie* htmie,
+ const int nc,
+ const char* var_name,
+ const double range[2],
+ double* value)
{
- size_t len;
size_t start;
+ size_t len;
int id;
int ndims;
int dimid;
int err = NC_NOERR;
int type;
res_T res = RES_OK;
- ASSERT(htmie);
+ ASSERT(htmie && nc != INVALID_NC_ID && range && var_name && value);
+ ASSERT(range[0] <= range[1]);
- err = nc_inq_varid(nc, "rmod", &id);
+ err = nc_inq_varid(nc, var_name, &id);
if(err != NC_NOERR) {
- log_err(htmie, "Could not inquiere the 'rmod' variable -- %s\n",
- ncerr_to_str(err));
+ log_err(htmie, "Could not inquire the '%s' variable -- %s\n",
+ var_name, ncerr_to_str(err));
res = RES_BAD_ARG;
goto error;
}
NC(inq_varndims(nc, id, &ndims));
if(ndims != 1) {
- log_err(htmie, "The dimension of the 'rmod' variable must be 1.\n");
+ log_err(htmie,
+ "The dimension of the '%s' variable must be 1.\n", var_name);
res = RES_BAD_ARG;
goto error;
}
NC(inq_vardimid(nc, id, &dimid));
-
NC(inq_dimlen(nc, dimid, &len));
if(len != 1) {
log_err(htmie,
- "Only 1 distribution is currently supported while the #distributions is "
- "%lu\n", (unsigned long)len);
+ "Only 1 distribution is currently supported while the #distributions of "
+ "the '%s' variable is %lu\n", var_name, (unsigned long)len);
res = RES_BAD_ARG;
goto error;
}
@@ -277,15 +297,24 @@ load_modal_mean_radius(struct htmie* htmie, const int nc)
NC(inq_vartype(nc, id, &type));
if(type != NC_DOUBLE && type != NC_FLOAT) {
log_err(htmie,
- "Type type of the 'rmod' variable cannot be %s. "
+ "The type of the '%s' variable cannot be %s. "
"Expecting floating point data.\n",
- nctype_to_str(type));
+ var_name, nctype_to_str(type));
res = RES_BAD_ARG;
goto error;
}
start = 0;
- NC(get_vara_double(nc, id, &start, &len, &htmie->rmod));
+ NC(get_vara_double(nc, id, &start, &len, value));
+
+ if(*value < range[0] || *value > range[1]) {
+ log_err(htmie,
+ "Invalid droplet distribution variable '%s = %g'. "
+ "It must be in [%g, %g]\n",
+ var_name, var_name, *value, range[0], range[1]);
+ res = RES_BAD_ARG;
+ goto error;
+ }
exit:
return res;
@@ -293,6 +322,132 @@ error:
goto exit;
}
+static res_T
+load_distrib_x_lambda_array
+ (struct htmie* htmie,
+ const int nc,
+ const char* var_name,
+ const double range[2], /* Validity range */
+ struct darray_double* values)
+{
+ size_t start[2];
+ size_t end[2];
+ size_t len;
+ size_t i;
+ int dimids[2];
+ int id;
+ int ndims;
+ int err;
+ int type;
+ res_T res = RES_OK;
+ ASSERT(htmie && nc != INVALID_NC_ID && var_name && values && range);
+ ASSERT(range[0] <= range[1]);
+
+ err = nc_inq_varid(nc, var_name, &id);
+ if(err != NC_NOERR) {
+ log_err(htmie, "Could not inquire the '%s' variable -- %s\n",
+ var_name, ncerr_to_str(err));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ NC(inq_varndims(nc, id, &ndims));
+ if(ndims != 2) {
+ log_err(htmie,
+ "The dimension of the '%s' variable must be 2.\n", var_name);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ NC(inq_vardimid(nc, id, dimids));
+ NC(inq_vartype(nc, id, &type));
+
+ NC(inq_dimlen(nc, dimids[0], &len));
+ if(len != 1) {
+ log_err(htmie,
+ "Only 1 distribution is currently supported while the #distributions of "
+ "the '%s' variable is %lu\n", var_name, (unsigned long)len);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ NC(inq_dimlen(nc, dimids[1], &len));
+ if(len != darray_double_size_get(&htmie->wavelengths)) {
+ log_err(htmie,
+ "Inconsistent wavelengths count. The number of defined wavelengths is %lu "
+ "while the per distribution length of the '%s' variable is %lu.\n",
+ darray_double_size_get(&htmie->wavelengths), var_name, len);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(type != NC_DOUBLE && type != NC_FLOAT) {
+ log_err(htmie,
+ "The type of the '%s' variable cannot be %s. "
+ "Expecting floating point data.\n",
+ var_name, nctype_to_str(type));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = darray_double_resize(values, len);
+ if(res != RES_OK) {
+ log_err(htmie,
+ "Could not allocate the memory space to load the data of the '%s' variable.\n",
+ var_name);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Read the raw data sections */
+ start[0] = 0, start[1] = 0;
+ end[0] = 1, end[1] = len;
+ NC(get_vara_double(nc, id, start, end, darray_double_data_get(values)));
+
+ FOR_EACH(i, 0, darray_double_size_get(values)) {
+ const double* d = darray_double_cdata_get(values);
+ if(d[i] < range[0] || d[i] > range[1]) {
+ log_err(htmie,
+ "Invalid value for the %lu^th data of the '%s' variable : %g. "
+ "The data must be in [%g, %g]\n",
+ (unsigned long)i, var_name, d[i], range[0], range[1]);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+
+exit:
+ return res;
+error:
+ darray_double_clear(values);
+ goto exit;
+}
+
+static INLINE res_T
+massic_to_per_particle_xsections
+ (const struct htmie* htmie,
+ const char* var_name,
+ struct darray_double* xsections)
+{
+ double radius_barre;
+ double avg_volume;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(htmie && var_name && xsections);
+
+ /* Convert data from massic cross sections to per particle cross sections */
+ radius_barre = exp(htmie->rmod); /* in micrometer */
+ avg_volume = /* in micrometer^3 */
+ 4.0/3.0 * PI * radius_barre*radius_barre*radius_barre
+ * exp(4.5*htmie->smod*htmie->smod);
+
+ FOR_EACH(i, 0, darray_double_size_get(xsections)) {
+ double* d = darray_double_data_get(xsections);
+ d[i] = d[i] * avg_volume * 1.e-15; /* in m^2/particle */
+ }
+ return res;
+}
+
static void
release_htmie(ref_T* ref)
{
@@ -302,6 +457,7 @@ release_htmie(ref_T* ref)
darray_double_release(&htmie->wavelengths);
darray_double_release(&htmie->xsections_absorption);
darray_double_release(&htmie->xsections_scattering);
+ darray_double_release(&htmie->g);
MEM_RM(htmie->allocator, htmie);
}
@@ -344,6 +500,7 @@ htmie_create
darray_double_init(htmie->allocator, &htmie->wavelengths);
darray_double_init(htmie->allocator, &htmie->xsections_absorption);
darray_double_init(htmie->allocator, &htmie->xsections_scattering);
+ darray_double_init(htmie->allocator, &htmie->g);
exit:
if(out_htmie) *out_htmie = htmie;
@@ -377,6 +534,7 @@ htmie_load(struct htmie* htmie, const char* path)
{
int err_nc = NC_NOERR;
int nc = INVALID_NC_ID;
+ double range[2];
res_T res = RES_OK;
if(!htmie || !path) {
@@ -391,10 +549,22 @@ htmie_load(struct htmie* htmie, const char* path)
goto error;
}
- res = load_wavelengths(htmie, nc);
- if(res != RES_OK) goto error;
- res = load_modal_mean_radius(htmie, nc);
- if(res != RES_OK) goto error;
+ /* Define the validity range for the loaded data */
+ range[0] = DBL_EPSILON;
+ range[1] = DBL_MAX;
+
+ #define CALL(Func) { res = Func; if(res != RES_OK) goto error; } (void)0
+ CALL(load_wavelengths(htmie, nc, range));
+ CALL(load_droplet_distribution_variable(htmie, nc, "rmod", range, &htmie->rmod));
+ CALL(load_droplet_distribution_variable(htmie, nc, "smod", range, &htmie->smod));
+ CALL(load_distrib_x_lambda_array(htmie, nc, "macs", range, &htmie->xsections_absorption));
+ CALL(load_distrib_x_lambda_array(htmie, nc, "mscs", range, &htmie->xsections_scattering));
+ CALL(load_distrib_x_lambda_array(htmie, nc, "g", range, &htmie->g));
+
+ /* The rmod and smod parameter must be loaded */
+ CALL(massic_to_per_particle_xsections(htmie, "macs", &htmie->xsections_absorption));
+ CALL(massic_to_per_particle_xsections(htmie, "mscs", &htmie->xsections_scattering));
+ #undef CALL
exit:
if(nc != INVALID_NC_ID) NC(close(nc));
