use interp1d() to interpolate path to speed things up

1 files changed, 27 insertions, 39 deletions

diff --git a/src/path.c b/src/path.c
index 1e53d79..a15cba6 100644
--- a/src/path.c
+++ b/src/path.c
@@ -6,6 +6,7 @@
 #include "pdg.h"
 #include "vector.h"
 #include <stdlib.h> /* for malloc(), calloc(), etc. */
+#include "misc.h"
 
 #define N 100
 
@@ -144,50 +145,34 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
         dz[i] = tmp2[2];
     }
 
-    for (i = 0; i < 8; i++) {
-        p->acc[i] = gsl_interp_accel_alloc();
-        p->spline[i] = gsl_spline_alloc(gsl_interp_linear,N);
-    }
-
-    gsl_spline_init(p->spline[0],s,x,N);
-    gsl_spline_init(p->spline[1],s,y,N);
-    gsl_spline_init(p->spline[2],s,z,N);
-    gsl_spline_init(p->spline[3],s,T,N);
-    gsl_spline_init(p->spline[4],s,t,N);
-    gsl_spline_init(p->spline[5],s,dx,N);
-    gsl_spline_init(p->spline[6],s,dy,N);
-    gsl_spline_init(p->spline[7],s,dz,N);
+    p->s = s;
+    p->x = x;
+    p->y = y;
+    p->z = z;
+    p->T = T;
+    p->t = t;
+    p->dx = dx;
+    p->dy = dy;
+    p->dz = dz;
 
-    free(s);
     free(theta1);
     free(theta2);
-    free(x);
-    free(y);
-    free(z);
-    free(T);
-    free(t);
-    free(dx);
-    free(dy);
-    free(dz);
 
     return p;
 }
 
 void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t, double *theta0)
 {
-    if (s > p->spline[0]->x[p->spline[0]->size-1])
-        s = p->spline[0]->x[p->spline[0]->size-1];
-
-    pos[0] = gsl_spline_eval(p->spline[0],s,p->acc[0]);
-    pos[1] = gsl_spline_eval(p->spline[1],s,p->acc[1]);
-    pos[2] = gsl_spline_eval(p->spline[2],s,p->acc[2]);
+    pos[0] = interp1d(s,p->s,p->x,N);
+    pos[1] = interp1d(s,p->s,p->y,N);
+    pos[2] = interp1d(s,p->s,p->z,N);
 
-    *T = gsl_spline_eval(p->spline[3],s,p->acc[3]);
-    *t = gsl_spline_eval(p->spline[4],s,p->acc[4]);
+    *T = interp1d(s,p->s,p->T,N);
+    *t = interp1d(s,p->s,p->t,N);
 
-    dir[0] = gsl_spline_eval(p->spline[5],s,p->acc[5]);
-    dir[1] = gsl_spline_eval(p->spline[6],s,p->acc[6]);
-    dir[2] = gsl_spline_eval(p->spline[7],s,p->acc[7]);
+    dir[0] = interp1d(s,p->s,p->dx,N);
+    dir[1] = interp1d(s,p->s,p->dy,N);
+    dir[2] = interp1d(s,p->s,p->dz,N);
 
     /* FIXME: This should be the *residual* scattering RMS. */
     *theta0 = p->theta0*sqrt(s);
@@ -195,12 +180,15 @@ void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t
 
 void path_free(path *p)
 {
-    size_t i;
-
-    for (i = 0; i < 5; i++) {
-        gsl_interp_accel_free(p->acc[i]);
-        gsl_spline_free(p->spline[i]);
-    }
+    free(p->s);
+    free(p->x);
+    free(p->y);
+    free(p->z);
+    free(p->T);
+    free(p->t);
+    free(p->dx);
+    free(p->dy);
+    free(p->dz);
 
     free(p);
 }