add absorption length for acrylic

1 files changed, 63 insertions, 0 deletions

diff --git a/src/optics.c b/src/optics.c
index 57c9d57..5f01089 100644
--- a/src/optics.c
+++ b/src/optics.c
@@ -24,6 +24,10 @@ 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};
 
+/* From prod/media.dat for Acrylic (good). */
+static double absorption_coefficient_acrylic_wavelengths[10] = {300.0, 310.0, 320.0, 330.0, 340.0, 350.0, 360.0, 380.0, 400.0, 450.0};
+static double absorption_coefficient_acrylic[10] = {0.3641, 0.0966, 0.0667, 0.049, 0.035, 0.0261, 0.0194, 0.0115, 0.0086, 0.0071};
+
 static int initialized = 0;
 
 /* D2O absorption length weighted by the Cerenkov spectrum and the quantum
@@ -32,6 +36,10 @@ 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;
+/* Acrylic absorption length weighted by the Cerenkov spectrum and the quantum
+ * efficiency. */
+static double absorption_length_acrylic_weighted = 0.0;
+
 
 /* From Table 4 in the paper. */
 static double A0 = 0.243905091;
@@ -80,6 +88,17 @@ static double gsl_absorption_length_snoman_h2o(double wavelength, void *params)
     return qe*get_absorption_length_snoman_h2o(wavelength)/pow(wavelength,2);
 }
 
+static double gsl_absorption_length_snoman_acrylic(double wavelength, void *params)
+{
+    /* Returns the absorption length in acrylic weighted by the quantum
+     * efficiency and the Cerenkov spectrum. */
+    double qe;
+
+    qe = get_quantum_efficiency(wavelength);
+
+    return qe*get_absorption_length_snoman_acrylic(wavelength)/pow(wavelength,2);
+}
+
 static double gsl_cerenkov(double wavelength, void *params)
 {
     /* Returns the quantum efficiency multiplied by the Cerenkov spectrum. */
@@ -133,6 +152,17 @@ int optics_init(void)
         exit(1);
     }
 
+    F.function = &gsl_absorption_length_snoman_acrylic;
+
+    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+
+    absorption_length_acrylic_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;
@@ -206,6 +236,39 @@ double get_absorption_length_snoman_d2o(double wavelength)
     return 1.0/gsl_spline_eval(spline, wavelength, acc);
 }
 
+double get_weighted_absorption_length_snoman_acrylic(void)
+{
+    /* Returns the weighted absorption length in acrylic in cm. */
+
+    if (!initialized) {
+        fprintf(stderr, "weighted absorption length hasn't been initialized!\n");
+        exit(1);
+    }
+
+    return absorption_length_acrylic_weighted;
+}
+
+double get_absorption_length_snoman_acrylic(double wavelength)
+{
+    /* Returns the absorption length in acrylic in cm. */
+    static gsl_spline *spline;
+    static gsl_interp_accel *acc;
+
+    if (!spline) {
+        spline = gsl_spline_alloc(gsl_interp_linear, LEN(absorption_coefficient_acrylic_wavelengths));
+        gsl_spline_init(spline, absorption_coefficient_acrylic_wavelengths, absorption_coefficient_acrylic, LEN(absorption_coefficient_acrylic_wavelengths));
+        acc = gsl_interp_accel_alloc();
+    }
+
+    if (wavelength < absorption_coefficient_acrylic_wavelengths[0])
+        wavelength = absorption_coefficient_acrylic_wavelengths[0];
+
+    if (wavelength > absorption_coefficient_acrylic_wavelengths[LEN(absorption_coefficient_acrylic_wavelengths)-1])
+        wavelength = absorption_coefficient_acrylic_wavelengths[LEN(absorption_coefficient_acrylic_wavelengths)-1];
+
+    return 1.0/gsl_spline_eval(spline, wavelength, acc);
+}
+
 double get_index(double p, double wavelength, double T)
 {
     /* Returns the index of refraction of pure water for a given density,