diff options
| -rw-r--r-- | src/path.c | 21 | ||||
| -rw-r--r-- | src/path.h | 3 |
2 files changed, 22 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) @@ -2,6 +2,7 @@ #define PATH_H #include <gsl/gsl_spline.h> +#include <stddef.h> /* for size_t */ /* Minimum value for the scattering RMS `theta0`. * @@ -20,6 +21,8 @@ typedef struct path { double dir[3]; double theta0; + size_t n; + double *s; double *x; double *y; |