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, 62 insertions, 0 deletions

diff --git a/src/test.c b/src/test.c
index e46963b..1a30976 100644
--- a/src/test.c
+++ b/src/test.c
@@ -13,6 +13,7 @@
 #include "pdg.h"
 #include <gsl/gsl_spline.h>
 #include "vector.h"
+#include <gsl/gsl_statistics.h>
 
 typedef int testFunction(char *err);
 
@@ -657,6 +658,65 @@ err:
     return 1;
 }
 
+int test_mean(char *err)
+{
+    /* Tests the mean() function. */
+    size_t i, j;
+    double x[100];
+    double mu, expected;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        for (j = 0; j < sizeof(x)/sizeof(x[0]); j++) {
+            x[j] = genrand_real2();
+        }
+
+        mu = mean(x,sizeof(x)/sizeof(x[0]));
+        expected = gsl_stats_mean(x,1,sizeof(x)/sizeof(x[0]));
+
+        if (!isclose(mu, expected, 0, 1e-9)) {
+            sprintf(err, "mean() returned %.5g, but expected %.5g", mu, expected);
+            goto err;
+        }
+    }
+
+    return 0;
+
+err:
+    return 1;
+}
+
+int test_std(char *err)
+{
+    /* Tests the mean() function. */
+    size_t i, j;
+    double x[100];
+    double sigma, expected, mu;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        for (j = 0; j < sizeof(x)/sizeof(x[0]); j++) {
+            x[j] = genrand_real2();
+        }
+
+        sigma = std(x,sizeof(x)/sizeof(x[0]));
+        mu = gsl_stats_mean(x,1,sizeof(x)/sizeof(x[0]));
+        expected = gsl_stats_sd_with_fixed_mean(x,1,sizeof(x)/sizeof(x[0]),mu);
+
+        if (!isclose(sigma, expected, 0, 1e-9)) {
+            sprintf(err, "std() returned %.5g, but expected %.5g", sigma, expected);
+            goto err;
+        }
+    }
+
+    return 0;
+
+err:
+    return 1;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -675,6 +735,8 @@ struct tests {
     {test_lnfact, "test_lnfact"},
     {test_ln, "test_ln"},
     {test_get_path_length, "test_get_path_length"},
+    {test_mean, "test_mean"},
+    {test_std, "test_std"},
 };
 
 int main(int argc, char **argv)