1 files changed, 7 insertions, 31 deletions

diff --git a/src/scattering.c b/src/scattering.c
index c4babae..25b17d4 100644
--- a/src/scattering.c
+++ b/src/scattering.c
@@ -18,8 +18,6 @@
 #include "quantum_efficiency.h"
 #include <gsl/gsl_sf_bessel.h>
 #include <gsl/gsl_integration.h>
-#include <gsl/gsl_interp2d.h>
-#include <gsl/gsl_spline2d.h>
 #include "optics.h"
 #include "sno.h"
 #include <stddef.h> /* for size_t */
@@ -27,8 +25,7 @@
 #include <stdio.h> /* for fprintf() */
 #include <gsl/gsl_errno.h> /* for gsl_strerror() */
 #include "pdg.h"
-#include <gsl/gsl_interp.h>
-#include <gsl/gsl_spline.h>
+#include "misc.h"
 
 static int initialized = 0;
 
@@ -43,10 +40,6 @@ static double *x;
 static double *y;
 static double *z;
 
-static gsl_spline2d *spline;
-static gsl_interp_accel *xacc;
-static gsl_interp_accel *yacc;
-
 static double x2lo = 0.0;
 static double x2hi = 1.0;
 static size_t nx2 = 1000;
@@ -57,9 +50,6 @@ static double *y2;
 /* Quantum efficiency weighted by the Cerenkov spectrum. */
 static double weighted_qe;
 
-static gsl_spline *spline2;
-static gsl_interp_accel *xacc2;
-
 static double prob_scatter(double wavelength, void *params)
 {
     /* Calculate the number of photons emitted per unit solid angle per cm at
@@ -132,10 +122,6 @@ void init_interpolation(void)
     y = malloc(ny*sizeof(double));
     z = malloc(nx*ny*sizeof(double));
 
-    spline = gsl_spline2d_alloc(gsl_interp2d_bilinear, nx, ny);
-    xacc = gsl_interp_accel_alloc();
-    yacc = gsl_interp_accel_alloc();
-
     for (i = 0; i < nx; i++) {
         x[i] = xlo + (xhi-xlo)*i/(nx-1);
     }
@@ -160,18 +146,14 @@ void init_interpolation(void)
                 fprintf(stderr, "error integrating photon angular distribution: %s\n", gsl_strerror(status));
                 exit(1);
             }
-            gsl_spline2d_set(spline, z, i, j, result);
+
+            z[i*ny + j] = result;
         }
     }
 
-    gsl_spline2d_init(spline, x, y, z, nx, ny);
-
     x2 = malloc(nx2*sizeof(double));
     y2 = malloc(nx2*sizeof(double));
 
-    spline2 = gsl_spline_alloc(gsl_interp_linear, nx2);
-    xacc2 = gsl_interp_accel_alloc();
-
     for (i = 0; i < nx2; i++) {
         x2[i] = x2lo + (x2hi-x2lo)*i/(nx2-1);
     }
@@ -192,8 +174,6 @@ void init_interpolation(void)
         y2[i] = result;
     }
 
-    gsl_spline_init(spline2, x2, y2, nx2);
-
     F.function = &gsl_quantum_efficiency;
     F.params = params;
 
@@ -223,12 +203,12 @@ double get_probability(double beta, double cos_theta, double theta0)
      * gsl_spline2d_eval() will quit. */
     if (theta0 > MAX_THETA0) theta0 = MAX_THETA0;
 
-    return gsl_spline2d_eval(spline, beta*cos_theta, beta*sin_theta*theta0, xacc, yacc)/(theta0*sin_theta);
+    return interp2d(beta*cos_theta, beta*sin_theta*theta0, x, y, z, nx, ny)/(theta0*sin_theta);
 }
 
 double get_probability2(double beta)
 {
-    return gsl_spline_eval(spline2, beta, xacc2);
+    return interp1d(beta, x2, y2, nx2);
 }
 
 void free_interpolation(void)
@@ -237,10 +217,6 @@ void free_interpolation(void)
     free(y);
     free(z);
 
-    gsl_spline2d_free(spline);
-    gsl_interp_accel_free(xacc);
-    gsl_interp_accel_free(yacc);
-
-    gsl_spline_free(spline2);
-    gsl_interp_accel_free(xacc2);
+    free(x2);
+    free(y2);
 }