add interp2d() for fast bilinear 2D interpolation

1 files changed, 62 insertions, 0 deletions

diff --git a/src/test.c b/src/test.c
index 242e728..f8cedb9 100644
--- a/src/test.c
+++ b/src/test.c
@@ -18,6 +18,7 @@
 #include "proton.h"
 #include "likelihood.h"
 #include "id_particles.h"
+#include <gsl/gsl_spline2d.h>
 
 typedef int testFunction(char *err);
 
@@ -909,6 +910,66 @@ err:
     return 1;
 }
 
+int test_interp2d(char *err)
+{
+    /* Tests the interp2d() function. */
+    size_t i, j;
+    double xp[100], yp[200], zp[20000], zp2[20000], range, x, y, z, expected;
+    gsl_interp_accel *xacc;
+    gsl_interp_accel *yacc;
+    gsl_spline2d *spline;
+
+    range = 1.0;
+
+    init_genrand(0);
+
+    for (i = 0; i < LEN(xp); i++) {
+        xp[i] = range*i/(LEN(xp)-1);
+    }
+
+    for (i = 0; i < LEN(yp); i++) {
+        yp[i] = range*i/(LEN(yp)-1);
+    }
+
+    spline = gsl_spline2d_alloc(gsl_interp2d_bilinear,LEN(xp),LEN(yp));
+
+    for (i = 0; i < LEN(xp); i++) {
+        for (j = 0; j < LEN(yp); j++) {
+            zp[i*LEN(yp)+j] = genrand_real2();
+            gsl_spline2d_set(spline, zp2, i, j, zp[i*LEN(yp)+j]);
+        }
+    }
+
+    xacc = gsl_interp_accel_alloc();
+    yacc = gsl_interp_accel_alloc();
+
+    gsl_spline2d_init(spline,xp,yp,zp2,LEN(xp),LEN(yp));
+
+    for (i = 0; i < 1000; i++) {
+        x = genrand_real2()*range;
+        y = genrand_real2()*range;
+        expected = gsl_spline2d_eval(spline,x,y,xacc,yacc);
+        z = interp2d(x,y,xp,yp,zp,LEN(xp),LEN(yp));
+        if (!isclose(z, expected, 1e-9, 1e-9)) {
+            sprintf(err, "interp2d() returned %.5g, but expected %.5g", z, expected);
+            goto err;
+        }
+    }
+
+    gsl_interp_accel_free(xacc);
+    gsl_interp_accel_free(yacc);
+    gsl_spline2d_free(spline);
+
+    return 0;
+
+err:
+    gsl_interp_accel_free(xacc);
+    gsl_interp_accel_free(yacc);
+    gsl_spline2d_free(spline);
+
+    return 1;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -937,6 +998,7 @@ struct tests {
     {test_proton_get_energy, "test_proton_get_energy"},
     {test_log_norm, "test_log_norm"},
     {test_trapz, "test_trapz"},
+    {test_interp2d, "test_interp2d"},
 };
 
 int main(int argc, char **argv)