htmie

commit 830e199c75c36875cf2ada239757f972b2bd00ca
parent 23a872c03a59bc04f20c51deb4a6a144036c70ec
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 17 Jul 2018 14:06:10 +0200

Fix the computation of the average cross sections

Diffstat:
M
src/htmie.c
 | 
43
++++++++++++++++++++++++++++++++-----------
M
src/test_htmie_load.c
 | 
20
+++++++++++++-------

2 files changed, 45 insertions(+), 18 deletions(-)
diff --git a/src/htmie.c b/src/htmie.c
@@ -34,6 +34,16 @@
     }                                                                          \
   } (void)0
 
+/* Context used in the sum_xsections functor */
+struct sum_xsections_context {
+  double prev_data;
+  double prev_wavelength;
+  double sum;
+  int first_entry;
+};
+static const struct sum_xsections_context SUM_XSECTIONS_CONTEXT_NULL =
+  {DBL_MAX, DBL_MAX, 0, 1};
+
 struct htmie {
   struct darray_double wavelengths; /* In nano-meters */
   struct darray_double xsections_absorption; /* In m^2.kg^-1 */
@@ -527,12 +537,21 @@ min_max(const double wavelength, const double data, void* ctx)
 }
 
 static FINLINE void
-add(const double wavelength, const double data, void* ctx)
+sum_xsections(const double wavelength, const double data, void* context)
 {
-  double* accum = ctx;
-  ASSERT(ctx);
-  (void)wavelength;
-  *accum += data;
+  struct sum_xsections_context* ctx = context;
+  ASSERT(context);
+
+  if(ctx->first_entry) {
+    ASSERT(ctx->sum == 0);
+    ctx->first_entry = 0;
+  } else {
+    ASSERT(wavelength > ctx->prev_wavelength);
+    ctx->sum +=
+      0.5 * (ctx->prev_data + data)*(wavelength - ctx->prev_wavelength);
+  }
+  ctx->prev_wavelength = wavelength;
+  ctx->prev_data = data;
 }
 
 static INLINE void
@@ -818,11 +837,12 @@ htmie_compute_xsection_absorption_average
    const enum htmie_filter_type type)
 {
   const double* xsections;
-  double sum = 0;
+  struct sum_xsections_context ctx = SUM_XSECTIONS_CONTEXT_NULL;
   ASSERT(htmie && band[0] < band[1]);
   xsections = darray_double_cdata_get(&htmie->xsections_absorption);
-  for_each_wavelength_in_spectral_band(htmie, band, type, xsections, add, &sum);
-  return sum / (band[1] - band[0]);
+  for_each_wavelength_in_spectral_band
+    (htmie, band, type, xsections, sum_xsections, &ctx);
+  return ctx.sum / (band[1] - band[0]);
 }
 
 double
@@ -832,11 +852,12 @@ htmie_compute_xsection_scattering_average
    const enum htmie_filter_type type)
 {
   const double* xsections;
-  double sum = 0;
+  struct sum_xsections_context ctx = SUM_XSECTIONS_CONTEXT_NULL;
   ASSERT(htmie && band[0] < band[1]);
   xsections = darray_double_cdata_get(&htmie->xsections_scattering);
-  for_each_wavelength_in_spectral_band(htmie, band, type, xsections, add, &sum);
-  return sum / (band[1] - band[0]);
+  for_each_wavelength_in_spectral_band
+    (htmie, band, type, xsections, sum_xsections, &ctx);
+  return ctx.sum / (band[1] - band[0]);
 }
 
 double
diff --git a/src/test_htmie_load.c b/src/test_htmie_load.c
@@ -187,8 +187,6 @@ static void
 test_avg(struct htmie* htmie)
 {
   double band[2];
-  double avg;
-  double ref;
   const double* wlens;
   const size_t ntests = 128;
   size_t itest;
@@ -214,14 +212,22 @@ test_avg(struct htmie* htmie)
     band[1] = wlens[iupp];
 
     #define TEST(Name, Filter) {                                               \
-      ref  = htmie_fetch_xsection_ ## Name(htmie, band[0], Filter);            \
-      ref += htmie_fetch_xsection_ ## Name(htmie, band[1], Filter);            \
+      double avg = 0;                                                          \
+      double sum = 0;                                                          \
+      double data = 0;                                                         \
+      double prev_wlen = band[0];                                              \
+      double prev_data = htmie_fetch_xsection_ ## Name(htmie, band[0], Filter);\
       FOR_EACH(i, ilow+1, iupp) {                                              \
-        ref += htmie_fetch_xsection_ ## Name(htmie, wlens[i], Filter);         \
+        data = htmie_fetch_xsection_ ## Name(htmie, wlens[i], Filter);         \
+        sum += 0.5*(prev_data + data) * (wlens[i] - prev_wlen);                \
+        prev_wlen = wlens[i];                                                  \
+        prev_data = data;                                                      \
       }                                                                        \
-      ref /= (band[1] - band[0]);                                              \
+      data = htmie_fetch_xsection_ ## Name(htmie, band[1], Filter);            \
+      sum += 0.5*(prev_data + data) * (band[1] - prev_wlen);                   \
+      sum /= (band[1] - band[0]);                                              \
       avg = htmie_compute_xsection_## Name ## _average(htmie, band, Filter);   \
-      CHK(eq_eps(avg, ref, 1.e-6));                                            \
+      CHK(eq_eps(avg, sum, 1.e-6));                                            \
     } (void)0
 
     TEST(absorption, NEAREST);