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

1 files changed, 16 insertions, 0 deletions

diff --git a/src/misc.c b/src/misc.c
index 6a67d58..2f02595 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -2,6 +2,22 @@
 #include <math.h>
 #include <stdlib.h> /* for size_t */
 
+double interp1d(double x, double *xp, double *yp, size_t n)
+{
+    /* A fast interpolation routine which assumes that the values in `xp` are
+     * evenly spaced.
+     *
+     * If x < xp[0] returns yp[0] and if x > xp[n-1] returns yp[n-1].  */
+    size_t i;
+
+    if (x < xp[0]) return yp[0];
+    if (x > xp[n-1]) return yp[n-1];
+
+    i = (x-xp[0])/(xp[1]-xp[0]);
+
+    return yp[i] + (yp[i+1]-yp[i])*(x-xp[i])/(xp[i+1]-xp[i]);
+}
+
 int isclose(double a, double b, double rel_tol, double abs_tol)
 {
     /* Returns 1 if a and b are "close". This algorithm is taken from Python's