From c79b3cf71e2e22f704e724e0c5ecb156550ac32e Mon Sep 17 00:00:00 2001
From: tlatorre <tlatorre@uchicago.edu>
Date: Tue, 29 Jan 2019 13:19:12 -0600
Subject: add a function to compute the angular distribution normalization

This seems to speed things up a little bit.
---
 src/electron.c | 38 +++-----------------------------------
 1 file changed, 3 insertions(+), 35 deletions(-)

(limited to 'src/electron.c')

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)
-- 
cgit