update theta0 calculation

This commit updates path_eval() to calculate theta0 using the residual
scattering RMS for a truncated KL expansion. Since there isn't a nice closed
form solution for this, we instead compute a rough approximation by evaluating
the residual scattering RMS at the center of the track.

1 files changed, 19 insertions, 2 deletions

diff --git a/src/path.c b/src/path.c
index 5a87cf6..673792c 100644
--- a/src/path.c
+++ b/src/path.c
@@ -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)