diff options
Diffstat (limited to 'src/path.c')
| -rw-r--r-- | src/path.c | 21 |
1 files changed, 19 insertions, 2 deletions
@@ -7,6 +7,7 @@ #include "vector.h" #include <stdlib.h> /* for malloc(), calloc(), etc. */ #include "misc.h" +#include <gsl/gsl_sf_psi.h> /* for gsl_sf_psi_n() */ #define N 1000 @@ -48,6 +49,8 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0 p->dir[1] = dir[1]; p->dir[2] = dir[2]; + p->n = n; + double *s = malloc(sizeof(double)*N); double *theta1 = calloc(N,sizeof(double)); double *theta2 = calloc(N,sizeof(double)); @@ -163,6 +166,8 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0 void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t, double *theta0) { + double range; + 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); @@ -174,8 +179,20 @@ void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t 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 = fmax(MIN_THETA0, p->theta0*sqrt(s)); + range = p->s[N-1]; + + /* Calculate the approximate residual scattering RMS. + * + * For a given truncated KL expansion the residual error in the + * approximation is given by: + * + * \sum_{n+1}^\infty sqrt(lambda_k)*sqrt(2/range)*sin(pi*s*(k-0.5)/range) + * + * Since there is no nice closed form solution for this, we estimate it by + * evaluating the result at the middle of the track. */ + *theta0 = 4*range*pow(p->theta0,2)*gsl_sf_psi_n(1,p->n-0.5)/pow(M_PI,2); + + *theta0 = fmax(MIN_THETA0, *theta0); } void path_free(path *p) |