add a function to compute the angular distribution normalization

This seems to speed things up a little bit.

1 files changed, 67 insertions, 0 deletions

diff --git a/src/test.c b/src/test.c
index 7053a94..5c44243 100644
--- a/src/test.c
+++ b/src/test.c
@@ -1976,6 +1976,72 @@ int test_argmin(char *err)
     return 0;
 }
 
+static double gsl_electron_get_angular_pdf_no_norm(double cos_theta, void *params)
+{
+    double alpha = ((double *) params)[0];
+    double beta = ((double *) params)[1];
+    double mu = ((double *) params)[2];
+    return electron_get_angular_pdf_no_norm(cos_theta,alpha,beta,mu);
+}
+
+static double electron_get_angular_pdf_norm_test(double alpha, double beta, double mu)
+{
+    double params[3];
+    gsl_integration_cquad_workspace *w;
+    gsl_function F;
+    double result, error;
+    int status;
+    size_t nevals;
+
+    w = gsl_integration_cquad_workspace_alloc(100);
+
+    F.function = &gsl_electron_get_angular_pdf_no_norm;
+    F.params = params;
+
+    params[0] = alpha;
+    params[1] = beta;
+    params[2] = mu;
+
+    status = gsl_integration_cquad(&F, -1, 1, 0, 1e-9, w, &result, &error, &nevals);
+
+    if (status) {
+        fprintf(stderr, "error integrating electron angular distribution: %s\n", gsl_strerror(status));
+        exit(1);
+    }
+
+    gsl_integration_cquad_workspace_free(w);
+
+    return result;
+}
+
+int test_electron_get_angular_pdf_norm(char *err)
+{
+    /* Tests that the electron_get_angular_pdf_norm() function returns the correct value. */
+    size_t i;
+    double a, b, mu, result, expected;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        a = fmax(0.1,genrand_real2());
+        b = genrand_real2();
+        mu = genrand_real2()*2 - 1;
+
+        result = electron_get_angular_pdf_norm(a,b,mu);
+        expected = electron_get_angular_pdf_norm_test(a,b,mu);
+
+        if (!isclose(result, expected, 0, 1e-9)) {
+            sprintf(err, "electron_get_angular_pdf_norm() returned %.5g, but expected %.5g", result, expected);
+            goto err;
+        }
+    }
+
+    return 0;
+
+err:
+    return 1;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -2026,6 +2092,7 @@ struct tests {
     {test_get_dir, "test_get_dir"},
     {test_argmax, "test_argmax"},
     {test_argmin, "test_argmin"},
+    {test_electron_get_angular_pdf_norm, "test_electron_get_angular_pdf_norm"},
 };
 
 int main(int argc, char **argv)