make epsabs and epsrel global variables

1 files changed, 13 insertions, 10 deletions

diff --git a/calculate_limits.c b/calculate_limits.c
index 07fe052..eefb6b8 100644
--- a/calculate_limits.c
+++ b/calculate_limits.c
@@ -3,7 +3,7 @@
 #include <math.h> /* For M_PI */
 
 /* Decay length of mediator V (in mm). */
-double decay_length = 1000.0;
+double decay_length = 1000e9;
 
 /* Mass of dark matter particle (MeV). */
 double mass = 1000.0;
@@ -27,6 +27,9 @@ double radius_fiducial = 5000;
  * they are used in some functions. */
 double x_sno[3];
 
+double epsabs = 1e-2;
+double epsrel = 1e-2;
+
 double deg2rad(double deg)
 {
     return deg*M_PI/180.0;
@@ -160,7 +163,7 @@ double f2(double theta, void *params)
     F.function = &f3;
     F.params = params;
 
-    gsl_integration_qags(&F, 0, 2*M_PI, 0, 1e-7, 1000, w, &result, &error);
+    gsl_integration_qags(&F, 0, 2*M_PI, epsabs, epsrel, 1000, w, &result, &error);
 
     gsl_integration_workspace_free(w);
 
@@ -180,7 +183,7 @@ double f1(double r, void *params)
     F.function = &f2;
     F.params = params;
 
-    gsl_integration_qags(&F, 0, M_PI, 0, 1e-7, 1000, w, &result, &error);
+    gsl_integration_qags(&F, 0, M_PI, epsabs, epsrel, 1000, w, &result, &error);
 
     gsl_integration_workspace_free(w);
 
@@ -201,10 +204,10 @@ double f4_earth(double phi_earth, void *params)
 
     gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
 
-    F.function = &f2;
+    F.function = &f1;
     F.params = params;
 
-    gsl_integration_qags(&F, 0, radius_fiducial, 0, 1e-7, 1000, w, &result, &error);
+    gsl_integration_qags(&F, 0, radius_fiducial, epsabs, epsrel, 1000, w, &result, &error);
 
     gsl_integration_workspace_free(w);
 
@@ -223,7 +226,7 @@ double f3_earth(double theta_earth, void *params)
     F.function = &f4_earth;
     F.params = params;
 
-    gsl_integration_qags(&F, 0, 2*M_PI, 0, 1e-7, 1000, w, &result, &error);
+    gsl_integration_qags(&F, 0, 2*M_PI, epsabs, epsrel, 1000, w, &result, &error);
 
     gsl_integration_workspace_free(w);
 
@@ -239,10 +242,10 @@ double f2_earth(double r_earth, void *params)
 
     gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
 
-    F.function = &f2;
+    F.function = &f3_earth;
     F.params = (void *) data;
 
-    gsl_integration_qags(&F, 0, M_PI, 0, 1e-7, 1000, w, &result, &error);
+    gsl_integration_qags(&F, 0, M_PI, epsabs, epsrel, 1000, w, &result, &error);
 
     gsl_integration_workspace_free(w);
 
@@ -259,7 +262,7 @@ double get_event_rate()
     F.function = &f2_earth;
     F.params = NULL;
 
-    gsl_integration_qags(&F, 0, radius_earth, 0, 1e-7, 1000, w, &result, &error);
+    gsl_integration_qags(&F, 0, radius_earth, epsabs, epsrel, 1000, w, &result, &error);
 
     gsl_integration_workspace_free(w);
 
@@ -276,7 +279,7 @@ int main(int argc, char **argv)
 
     sphere2cartesian(radius_earth - sno_depth, sno_theta, sno_phi, x_sno, x_sno+1, x_sno+2);
 
-    printf("volume = %.2g\n", get_event_rate());
+    printf("volume = %.18e\n", get_event_rate());
 
     return 0;
 }