fix a bug when computing the absorption length in H2O and D2O

Previously I was interpolating the absorption lengths using interp1d() but that
only works when the x array is uniform. Since the wavelengths are not spaced
uniformly, we have to use the GSL interpolation routines.

1 files changed, 38 insertions, 2 deletions

diff --git a/src/optics.c b/src/optics.c
index b01ce7d..57c9d57 100644
--- a/src/optics.c
+++ b/src/optics.c
@@ -5,6 +5,12 @@
 #include <stdio.h> /* for fprintf() */
 #include <gsl/gsl_integration.h>
 #include <gsl/gsl_errno.h> /* for gsl_strerror() */
+#include <gsl/gsl_spline.h>
+
+/* Macro to compute the size of a static C array.
+ *
+ * See https://stackoverflow.com/questions/1598773. */
+#define LEN(x) ((sizeof(x)/sizeof(0[x]))/((size_t)(!(sizeof(x) % sizeof(0[x])))))
 
 /* Absorption coefficients for H2O and D2O as a function of wavelength from
  * SNOMAN.
@@ -149,7 +155,22 @@ double get_weighted_absorption_length_snoman_h2o(void)
 double get_absorption_length_snoman_h2o(double wavelength)
 {
     /* Returns the absorption length in H2O in cm. */
-    return 1.0/interp1d(wavelength, absorption_coefficient_h2o_wavelengths, absorption_coefficient_h2o, sizeof(absorption_coefficient_h2o_wavelengths)/sizeof(absorption_coefficient_h2o_wavelengths[0]));
+    static gsl_spline *spline;
+    static gsl_interp_accel *acc;
+
+    if (!spline) {
+        spline = gsl_spline_alloc(gsl_interp_linear, LEN(absorption_coefficient_h2o_wavelengths));
+        gsl_spline_init(spline, absorption_coefficient_h2o_wavelengths, absorption_coefficient_h2o, LEN(absorption_coefficient_h2o_wavelengths));
+        acc = gsl_interp_accel_alloc();
+    }
+
+    if (wavelength < absorption_coefficient_h2o_wavelengths[0])
+        wavelength = absorption_coefficient_h2o_wavelengths[0];
+
+    if (wavelength > absorption_coefficient_h2o_wavelengths[LEN(absorption_coefficient_h2o_wavelengths)-1])
+        wavelength = absorption_coefficient_h2o_wavelengths[LEN(absorption_coefficient_h2o_wavelengths)-1];
+
+    return 1.0/gsl_spline_eval(spline, wavelength, acc);
 }
 
 double get_weighted_absorption_length_snoman_d2o(void)
@@ -167,7 +188,22 @@ double get_weighted_absorption_length_snoman_d2o(void)
 double get_absorption_length_snoman_d2o(double wavelength)
 {
     /* Returns the absorption length in D2O in cm. */
-    return 1.0/interp1d(wavelength, absorption_coefficient_d2o_wavelengths, absorption_coefficient_d2o, sizeof(absorption_coefficient_d2o_wavelengths)/sizeof(absorption_coefficient_d2o_wavelengths[0]));
+    static gsl_spline *spline;
+    static gsl_interp_accel *acc;
+
+    if (!spline) {
+        spline = gsl_spline_alloc(gsl_interp_linear, LEN(absorption_coefficient_d2o_wavelengths));
+        gsl_spline_init(spline, absorption_coefficient_d2o_wavelengths, absorption_coefficient_d2o, LEN(absorption_coefficient_d2o_wavelengths));
+        acc = gsl_interp_accel_alloc();
+    }
+
+    if (wavelength < absorption_coefficient_d2o_wavelengths[0])
+        wavelength = absorption_coefficient_d2o_wavelengths[0];
+
+    if (wavelength > absorption_coefficient_d2o_wavelengths[LEN(absorption_coefficient_d2o_wavelengths)-1])
+        wavelength = absorption_coefficient_d2o_wavelengths[LEN(absorption_coefficient_d2o_wavelengths)-1];
+
+    return 1.0/gsl_spline_eval(spline, wavelength, acc);
 }
 
 double get_index(double p, double wavelength, double T)