diff options
-rw-r--r-- | src/path.c | 66 | ||||
-rw-r--r-- | src/path.h | 12 |
2 files changed, 37 insertions, 41 deletions
@@ -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); } @@ -9,8 +9,16 @@ typedef struct path { double pos[3]; double dir[3]; double theta0; - gsl_interp_accel *acc[8]; - gsl_spline *spline[8]; + + double *s; + double *x; + double *y; + double *z; + double *T; + double *t; + double *dx; + double *dy; + double *dz; } path; double path_get_coefficient(unsigned int k, double *s, double *x, double theta0, size_t n); |