1 files changed, 48 insertions, 0 deletions
diff --git a/src/optics.c b/src/optics.c
index 830ada9..6b6d761 100644
--- a/src/optics.c
+++ b/src/optics.c
@@ -29,6 +29,8 @@
 
 char optics_err[256];
 
+double avg_index_h2o, avg_index_d2o;
+
 /* Absorption coefficients for H2O, D2O, and acrylic as a function of
  * wavelength from SNOMAN.
  *
@@ -425,6 +427,28 @@ static double gsl_absorption_length_snoman_acrylic(double wavelength, void *para
     return qe*exp(-x/get_absorption_length_snoman_acrylic(wavelength))/pow(wavelength,2);
 }
 
+/* Returns the product of the index of refraction in D2O multiplied by the
+ * quantum efficiency and the Cerenkov spectrum. */
+static double gsl_index_d2o(double wavelength, void *params)
+{
+    double qe;
+
+    qe = get_quantum_efficiency(wavelength);
+
+    return qe*get_index_snoman_d2o(wavelength)/pow(wavelength,2);
+}
+
+/* Returns the product of the index of refraction in H2O multiplied by the
+ * quantum efficiency and the Cerenkov spectrum. */
+static double gsl_index_h2o(double wavelength, void *params)
+{
+    double qe;
+
+    qe = get_quantum_efficiency(wavelength);
+
+    return qe*get_index_snoman_h2o(wavelength)/pow(wavelength,2);
+}
+
 static double gsl_cerenkov(double wavelength, void *params)
 {
     /* Returns the quantum efficiency multiplied by the Cerenkov spectrum. */
@@ -458,6 +482,30 @@ int optics_init(dict *db)
         return -1;
     }
 
+    F.function = &gsl_index_d2o;
+    F.params = NULL;
+
+    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+
+    if (status) {
+        sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+        return -1;
+    }
+
+    avg_index_d2o = result/norm;
+
+    F.function = &gsl_index_h2o;
+    F.params = NULL;
+
+    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+
+    if (status) {
+        sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+        return -1;
+    }
+
+    avg_index_h2o = result/norm;
+
     for (i = 0; i < LEN(fabs_x); i++) {
         fabs_x[i] = xlo + (xhi-xlo)*i/(LEN(fabs_x)-1);
     }