add a function to compute the angular distribution normalization

This seems to speed things up a little bit.

3 files changed, 72 insertions, 35 deletions

diff --git a/src/electron.c b/src/electron.c
index c6b5ed8..73fe24e 100644
--- a/src/electron.c
+++ b/src/electron.c
@@ -17,6 +17,7 @@
 #include "misc.h"
 #include <gsl/gsl_integration.h>
 #include <math.h> /* for fmax() */
+#include <gsl/gsl_sf_gamma.h>
 
 static int initialized = 0;
 
@@ -123,42 +124,9 @@ double electron_get_angular_pdf_no_norm(double cos_theta, double alpha, double b
     return exp(-pow(fabs(cos_theta-mu),alpha)/beta);
 }
 
-static double gsl_electron_get_angular_pdf(double cos_theta, void *params)
+double electron_get_angular_pdf_norm(double alpha, double beta, double mu)
 {
-    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(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;
-    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;
+    return pow(beta,1/alpha)*gsl_sf_gamma(1/alpha)*(gsl_sf_gamma_inc_P(1/alpha,pow(1-mu,alpha)/beta)+gsl_sf_gamma_inc_P(1/alpha,pow(1+mu,alpha)/beta))/alpha;
 }
 
 double electron_get_angular_pdf_delta_ray(double cos_theta, double alpha, double beta, double mu)
diff --git a/src/electron.h b/src/electron.h
index 20e5ea2..61de2a3 100644
--- a/src/electron.h
+++ b/src/electron.h
@@ -13,6 +13,8 @@
 void electron_get_position_distribution_parameters(double T0, double *a, double *b);
 double electron_get_angular_distribution_alpha(double T0);
 double electron_get_angular_distribution_beta(double T0);
+double electron_get_angular_pdf_no_norm(double cos_theta, double alpha, double beta, double mu);
+double electron_get_angular_pdf_norm(double alpha, double beta, double mu);
 double electron_get_angular_pdf_delta_ray(double cos_theta, double alpha, double beta, double mu);
 double electron_get_angular_pdf(double cos_theta, double alpha, double beta, double mu);
 double electron_get_range(double T, double rho);
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)