1 files changed, 3 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)