From 643204e807d5e78f883fc30dc7383a209e86dbc5 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Wed, 17 Oct 2018 09:33:29 -0500 Subject: 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 --- src/path.c | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) (limited to 'src/path.c') 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 /* for malloc(), calloc(), etc. */ #include "misc.h" #include /* 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) -- cgit