htmie

commit dc8518ff8d2f4ae34af2ab71cff43fe0293c2813
parent 7335d853d797ee252a3c432aceadc95a880dd8a7
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue,  3 Jul 2018 12:20:28 +0200

Implement accessors toward the loaded data

Diffstat:
M
src/htmie.c
 | 
135
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-----
M
src/htmie.h
 | 
44
++++++++++++++++++++++++++++++++++++++++++++

2 files changed, 171 insertions(+), 8 deletions(-)
diff --git a/src/htmie.c b/src/htmie.c
@@ -15,6 +15,7 @@
 
 #include "htmie.h"
 
+#include <rsys/algorithm.h>
 #include <rsys/dynamic_array_double.h>
 #include <rsys/logger.h>
 #include <rsys/mem_allocator.h>
@@ -426,14 +427,13 @@ error:
 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);
+  ASSERT(htmie && xsections);
 
   /* Convert data from massic cross sections to per particle cross sections */
   radius_barre = exp(htmie->rmod); /* in micrometer */
@@ -448,6 +448,52 @@ massic_to_per_particle_xsections
   return res;
 }
 
+static FINLINE double
+fetch_data
+  (const struct htmie* htmie,
+   const double wavelength,
+   const enum htmie_filter_type type,
+   const double* data) /* Data to interpolate */
+{
+  size_t nwlens;
+  const double* wlens;
+  const double* wlen_upp;
+  double val;
+  ASSERT(htmie && (unsigned)type < HTMIE_FILTER_TYPES_COUNT__ && data);
+
+  wlens = htmie_get_wavelengths(htmie);
+  nwlens = htmie_get_wavelengths_count(htmie);
+  ASSERT(nwlens);
+
+  wlen_upp = search_lower_bound
+    (&wavelength, wlens, nwlens, sizeof(double), cmp_dbl);
+
+  if(!wlen_upp) { /* Clamp to upper */
+    val = data[nwlens-1];
+  } else if(wlen_upp == wlens) { /* Clamp to lower */
+    val = data[0];
+  } else {
+    const size_t i = (size_t)(wlen_upp - wlens);
+    const double a = data[i-1];
+    const double b = data[i-0];
+    const double u = (wavelength - wlens[i-1]) / (wlens[i] - wlens[i-1]);
+    ASSERT(i < nwlens && (wlens[i] >= wavelength || wlens[i-1] <= wavelength));
+    ASSERT(eq_eps(u, 1, 1.e-6));
+
+    switch(type) {
+      case HTMIE_FILTER_NEAREST:
+        val = u < 0.5 ? a : b;
+        break;
+      case HTMIE_FILTER_LINEAR:
+        CLAMP(u, 0, 1); /* Handle numerical inaccuracy */
+        val = u*b + (1-u)*a;
+        break;
+      default: FATAL("Unreachable code.\n"); break;
+    }
+  }
+  return val;
+}
+
 static void
 release_htmie(ref_T* ref)
 {
@@ -555,15 +601,19 @@ htmie_load(struct htmie* htmie, const char* path)
 
   #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_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));
+  CALL(massic_to_per_particle_xsections(htmie, &htmie->xsections_absorption));
+  CALL(massic_to_per_particle_xsections(htmie, &htmie->xsections_scattering));
   #undef CALL
 
 exit:
@@ -572,3 +622,72 @@ exit:
 error:
   goto exit;
 }
+
+const double*
+htmie_get_wavelengths(const struct htmie* htmie)
+{
+  ASSERT(htmie);
+  return darray_double_cdata_get(&htmie->wavelengths);
+}
+
+const double*
+htmie_get_xsections_absorption(const struct htmie* htmie)
+{
+  ASSERT(htmie);
+  return darray_double_cdata_get(&htmie->xsections_absorption);
+}
+
+const double*
+htmie_get_xsections_scattering(const struct htmie* htmie)
+{
+  ASSERT(htmie);
+  return darray_double_cdata_get(&htmie->xsections_scattering);
+}
+
+const double*
+htmie_get_asymmetric_parameter(const struct htmie* htmie)
+{
+  ASSERT(htmie);
+  return darray_double_cdata_get(&htmie->g);
+}
+
+size_t
+htmie_get_wavelengths_count(const struct htmie* htmie)
+{
+  ASSERT(htmie);
+  return darray_double_size_get(&htmie->wavelengths);
+}
+
+double
+htmie_fetch_xsection_absorption
+  (const struct htmie* htmie,
+   const double wavelength,
+   const enum htmie_filter_type type)
+{
+  ASSERT(htmie);
+  return fetch_data
+    (htmie, wavelength, type, htmie_get_xsections_absorption(htmie));
+}
+
+double
+htmie_fetch_xsection_scattering
+  (const struct htmie* htmie,
+   const double wavelength,
+   const enum htmie_filter_type type)
+{
+  ASSERT(htmie);
+  return fetch_data
+    (htmie, wavelength, type, htmie_get_xsections_scattering(htmie));
+}
+
+double
+htmie_fetch_asymmetric_parameter
+  (const struct htmie* htmie,
+   const double wavelength,
+   const enum htmie_filter_type type)
+{
+  ASSERT(htmie);
+  return fetch_data
+    (htmie, wavelength, type, htmie_get_asymmetric_parameter(htmie));
+}
+
diff --git a/src/htmie.h b/src/htmie.h
@@ -40,6 +40,12 @@
 struct logger;
 struct mem_allocator;
 
+enum htmie_filter_type {
+  HTMIE_FILTER_LINEAR,
+  HTMIE_FILTER_NEAREST,
+  HTMIE_FILTER_TYPES_COUNT__
+};
+
 /* Forward declaration of opaque data types */
 struct htmie;
 
@@ -68,6 +74,44 @@ htmie_load
   (struct htmie *htmie,
    const char* path);
 
+HTMIE_API const double*
+htmie_get_wavelengths
+  (const struct htmie* htmie);
+
+HTMIE_API const double*
+htmie_get_xsections_absorption
+  (const struct htmie* htmie);
+
+HTMIE_API const double*
+htmie_get_xsections_scattering
+  (const struct htmie* htmie);
+
+HTMIE_API const double*
+htmie_get_asymmetric_parameter
+  (const struct htmie* htmie);
+
+HTMIE_API size_t
+htmie_get_wavelengths_count
+  (const struct htmie* htmie);
+
+HTMIE_API double
+htmie_fetch_xsection_absorption
+  (const struct htmie* htmie,
+   const double wavelength,
+   const enum htmie_filter_type type);
+
+HTMIE_API double
+htmie_fetch_xsection_scattering
+  (const struct htmie* htmie,
+   const double wavelength,
+   const enum htmie_filter_type type);
+
+HTMIE_API double
+htmie_fetch_asymmetric_parameter
+  (const struct htmie* htmie,
+   const double wavelength,
+   const enum htmie_filter_type type);
+
 END_DECLS
 
 #endif /* HTMIE_H */