add interp1d function to do fast interpolation when the x values are evenly spaced

1 files changed, 46 insertions, 1 deletions

diff --git a/src/test.c b/src/test.c
index cc054c0..2ea73d1 100644
--- a/src/test.c
+++ b/src/test.c
@@ -11,6 +11,7 @@
 #include "path.h"
 #include "random.h"
 #include "pdg.h"
+#include <gsl/gsl_spline.h>
 
 typedef int testFunction(char *err);
 
@@ -396,6 +397,49 @@ int test_path(char *err)
     return 0;
 }
 
+int test_interp1d(char *err)
+{
+    /* Tests the interp1d() function. */
+    size_t i;
+    double xp[100], yp[100], range, x, y, expected;
+    gsl_interp_accel *acc;
+    gsl_spline *spline;
+
+    range = 1.0;
+
+    init_genrand(0);
+
+    for (i = 0; i < sizeof(xp)/sizeof(xp[0]); i++) {
+        xp[i] = range*i/(100-1);
+        yp[i] = genrand_real2();
+    }
+
+    acc = gsl_interp_accel_alloc();
+    spline = gsl_spline_alloc(gsl_interp_linear,100);
+    gsl_spline_init(spline,xp,yp,100);
+
+    for (i = 0; i < 1000; i++) {
+        x = genrand_real2()*range;
+        expected = gsl_spline_eval(spline,x,acc);
+        y = interp1d(x,xp,yp,100);
+        if (!isclose(y, expected, 1e-9, 1e-9)) {
+            sprintf(err, "interp1d(%.2g) returned %.5g, but expected %.5g", x, y, expected);
+            goto err;
+        }
+    }
+
+    gsl_interp_accel_free(acc);
+    gsl_spline_free(spline);
+
+    return 0;
+
+err:
+    gsl_interp_accel_free(acc);
+    gsl_spline_free(spline);
+
+    return 1;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -408,7 +452,8 @@ struct tests {
     {test_muon_get_T, "test_muon_get_T"},
     {test_logsumexp, "test_logsumexp"},
     {test_sno_charge_integral, "test_sno_charge_integral"},
-    {test_path, "test_path"}
+    {test_path, "test_path"},
+    {test_interp1d, "test_interp1d"}
 };
 
 int main(int argc, char **argv)