fix a bug in the theta0 calculation for a path

This commit fixes a bug in the calculation of the average rms width of the
angular distribution for a path with a KL expansion. I also made a lot of
updates to the test-path program:

- plot the distribution of the KL expansion coefficients
- plot the standard deviation of the angular distribution as a function of
  distance along with the prediction
- plot the simulated and reconstructed path in 3D

1 files changed, 7 insertions, 3 deletions

diff --git a/src/path.c b/src/path.c
index d26e461..0aeb20b 100644
--- a/src/path.c
+++ b/src/path.c
@@ -8,6 +8,7 @@
 #include <stdlib.h> /* for malloc(), calloc(), etc. */
 #include "misc.h"
 #include <gsl/gsl_sf_psi.h> /* for gsl_sf_psi_n() */
+#include "scattering.h"
 
 #define N 1000
 
@@ -156,8 +157,9 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
      * Since there is no nice closed form solution for this, we estimate it by
      * evaluating the result at the middle of the track. */
     if (p->n > 0) {
-        p->theta0 = 4*range*pow(theta0,2)*gsl_sf_psi_n(1,p->n-0.5)/pow(M_PI,2);
+        p->theta0 = theta0*sqrt(range)*sqrt(gsl_sf_psi_n(1,p->n+0.5))/M_PI;
         p->theta0 = fmax(MIN_THETA0, p->theta0);
+        p->theta0 = fmin(MAX_THETA0, p->theta0);
     } else {
         p->theta0 = theta0;
     }
@@ -198,10 +200,12 @@ void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t
 
     normalize(dir);
 
-    if (p->n > 0)
+    if (p->n > 0) {
         *theta0 = p->theta0;
-    else
+    } else {
         *theta0 = fmax(MIN_THETA0,p->theta0*sqrt(s));
+        *theta0 = fmin(MAX_THETA0,*theta0);
+    }
 }
 
 void path_free(path *p)