add absorption lengths for D2O and H2O weighted by the Cerenkov spectrum and the quantum efficiency

1 files changed, 119 insertions, 0 deletions

diff --git a/src/optics.c b/src/optics.c
index af89802..b01ce7d 100644
--- a/src/optics.c
+++ b/src/optics.c
@@ -1,6 +1,10 @@
 #include <math.h>
 #include "optics.h"
 #include "misc.h"
+#include "quantum_efficiency.h"
+#include <stdio.h> /* for fprintf() */
+#include <gsl/gsl_integration.h>
+#include <gsl/gsl_errno.h> /* for gsl_strerror() */
 
 /* Absorption coefficients for H2O and D2O as a function of wavelength from
  * SNOMAN.
@@ -14,6 +18,15 @@ static double absorption_coefficient_h2o[7] = {1e-5, 1e-5, 1.3e-5, 6.0e-5, 2.0e-
 static double absorption_coefficient_d2o_wavelengths[7] = {254.0, 313.0, 366.0, 406.0, 436.0, 548.0, 578.0};
 static double absorption_coefficient_d2o[7] = {1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5};
 
+static int initialized = 0;
+
+/* D2O absorption length weighted by the Cerenkov spectrum and the quantum
+ * efficiency. */
+static double absorption_length_d2o_weighted = 0.0;
+/* H2O absorption length weighted by the Cerenkov spectrum and the quantum
+ * efficiency. */
+static double absorption_length_h2o_weighted = 0.0;
+
 /* From Table 4 in the paper. */
 static double A0 = 0.243905091;
 static double A1 = 9.53518094e-3;
@@ -39,12 +52,118 @@ static double RIND_D2O_C3 = 0.32;
  * From http://pdg.lbl.gov/2018/reviews/rpp2018-rev-phys-constants.pdf. */
 static double HC = 1.239841973976e3;
 
+static double gsl_absorption_length_snoman_d2o(double wavelength, void *params)
+{
+    /* Returns the absorption length in D2O weighted by the quantum efficiency
+     * and the Cerenkov spectrum. */
+    double qe;
+
+    qe = get_quantum_efficiency(wavelength);
+
+    return qe*get_absorption_length_snoman_d2o(wavelength)/pow(wavelength,2);
+}
+
+static double gsl_absorption_length_snoman_h2o(double wavelength, void *params)
+{
+    /* Returns the absorption length in H2O weighted by the quantum efficiency
+     * and the Cerenkov spectrum. */
+    double qe;
+
+    qe = get_quantum_efficiency(wavelength);
+
+    return qe*get_absorption_length_snoman_h2o(wavelength)/pow(wavelength,2);
+}
+
+static double gsl_cerenkov(double wavelength, void *params)
+{
+    /* Returns the quantum efficiency multiplied by the Cerenkov spectrum. */
+    double qe;
+
+    qe = get_quantum_efficiency(wavelength);
+
+    return qe/pow(wavelength,2);
+}
+
+int optics_init(void)
+{
+    double norm;
+    double result, error;
+    size_t nevals;
+    int status;
+    gsl_integration_cquad_workspace *w;
+    gsl_function F;
+
+    w = gsl_integration_cquad_workspace_alloc(100);
+
+    F.function = &gsl_cerenkov;
+    F.params = NULL;
+
+    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &norm, &error, &nevals);
+
+    if (status) {
+        fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+        exit(1);
+    }
+
+    F.function = &gsl_absorption_length_snoman_d2o;
+
+    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+
+    absorption_length_d2o_weighted = result/norm;
+
+    if (status) {
+        fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+        exit(1);
+    }
+
+    F.function = &gsl_absorption_length_snoman_h2o;
+
+    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+
+    absorption_length_h2o_weighted = result/norm;
+
+    if (status) {
+        fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+        exit(1);
+    }
+
+    gsl_integration_cquad_workspace_free(w);
+
+    initialized = 1;
+
+    return 0;
+}
+
+double get_weighted_absorption_length_snoman_h2o(void)
+{
+    /* Returns the weighted absorption length in H2O in cm. */
+
+    if (!initialized) {
+        fprintf(stderr, "weighted absorption length hasn't been initialized!\n");
+        exit(1);
+    }
+
+    return absorption_length_h2o_weighted;
+}
+
 double get_absorption_length_snoman_h2o(double wavelength)
 {
     /* Returns the absorption length in H2O in cm. */
     return 1.0/interp1d(wavelength, absorption_coefficient_h2o_wavelengths, absorption_coefficient_h2o, sizeof(absorption_coefficient_h2o_wavelengths)/sizeof(absorption_coefficient_h2o_wavelengths[0]));
 }
 
+double get_weighted_absorption_length_snoman_d2o(void)
+{
+    /* Returns the weighted absorption length in D2O in cm. */
+
+    if (!initialized) {
+        fprintf(stderr, "weighted absorption length hasn't been initialized!\n");
+        exit(1);
+    }
+
+    return absorption_length_d2o_weighted;
+}
+
 double get_absorption_length_snoman_d2o(double wavelength)
 {
     /* Returns the absorption length in D2O in cm. */