start working on the code to numerically integrate to find the event rate

2 files changed, 79 insertions, 1 deletions

diff --git a/Makefile b/Makefile
index 8f1188a..3b4ff53 100644
--- a/Makefile
+++ b/Makefile
@@ -1,3 +1,3 @@
-LDLIBS=-lm
+LDLIBS=-lm -lgsl -lgslcblas
 
 calculate_limits: calculate_limits.c
diff --git a/calculate_limits.c b/calculate_limits.c
index 47ed1c9..c4919d0 100644
--- a/calculate_limits.c
+++ b/calculate_limits.c
@@ -2,6 +2,12 @@
 #include <gsl/gsl_integration.h>
 #include <math.h> /* For M_PI */
 
+/* Decay length of mediator V (in mm). */
+double decay_length = 1000.0;
+
+/* Mass of dark matter particle (MeV). */
+double mass = 1000.0;
+
 /* From Google maps. Probably not very accurate, but should be good enough for
  * this calculation. */
 double latitude = 46.471857;
@@ -112,6 +118,76 @@ void rotate_earth_to_sno(double *x_earth, double *x_sno)
     rotate(x_sno, x_earth, dir, theta);
 }
 
+/* Integral over phi. */
+double f3(double phi, void *params)
+{
+    double result, error;
+    gsl_function F;
+    double *data = (double *) params;
+    data[2] = phi;
+
+    return 1;
+}
+
+/* Integral over theta. */
+double f2(double theta, void *params)
+{
+    double result, error;
+    gsl_function F;
+    double *data = (double *) params;
+    data[1] = theta;
+
+    gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
+
+    F.function = &f3;
+    F.params = params;
+
+    gsl_integration_qags(&F, 0, 2*M_PI, 0, 1e-7, 1000, w, &result, &error);
+
+    gsl_integration_workspace_free(w);
+
+    return result;
+}
+
+/* Integral over r. */
+double f1(double r, void *params)
+{
+    double result, error;
+    gsl_function F;
+    double *data = (double *) params;
+    data[0] = r;
+
+    gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
+
+    F.function = &f2;
+    F.params = params;
+
+    gsl_integration_qags(&F, 0, M_PI, 0, 1e-7, 1000, w, &result, &error);
+
+    gsl_integration_workspace_free(w);
+
+    return result;
+}
+
+/* Integral over r. */
+double get_event_rate()
+{
+    double result, error;
+    double params[3];
+    gsl_function F;
+
+    gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
+
+    F.function = &f2;
+    F.params = params;
+
+    gsl_integration_qags(&F, 0, radius_fiducial, 0, 1e-7, 1000, w, &result, &error);
+
+    gsl_integration_workspace_free(w);
+
+    return result;
+}
+
 int main(int argc, char **argv)
 {
     /* Spherical angles for the SNO detector in the earth frame which has z
@@ -122,5 +198,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());
+
     return 0;
 }