fix a bug in the absorption and scattering probabilities

Previously I was computing the fraction of light absorbed and scattered by
calculating an average absorption and scattering length weighted by the
Cerenkov spectrum and the PMT quantum efficiency, which isn't correct since we
should be averaging the absorption and scattering probabilities, not the
absorption and scattering lengths.

This commit fixes this by instead computing the average probability that a
photon is absorbed or scattered as a function of the distance travelled by
integrating the absorption and scattering probabilities over all wavelengths
weighted by the PMT quantum efficiency and the Cerenkov spectrum.

1 files changed, 360 insertions, 240 deletions

diff --git a/src/optics.c b/src/optics.c
index d626bf7..830ada9 100644
--- a/src/optics.c
+++ b/src/optics.c
@@ -25,11 +25,12 @@
 #include "misc.h"
 #include "db.h"
 #include "dict.h"
+#include "sno.h"
 
 char optics_err[256];
 
-/* Absorption coefficients for H2O and D2O as a function of wavelength from
- * SNOMAN.
+/* Absorption coefficients for H2O, D2O, and acrylic as a function of
+ * wavelength from SNOMAN.
  *
  * Note: These numbers come from mcprod/media_qoca_d2o_20060110.cmd. */
 static double absorption_coefficient_h2o_wavelengths[6] = {337.0, 365.0, 386.0, 420.0, 500.0, 620.0};
@@ -47,21 +48,32 @@ static double absorption_coefficient_acrylic[6] = {5.610e-02, 2.279e-02, 1.204e-
 
 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;
-/* Acrylic absorption length weighted by the Cerenkov spectrum and the quantum
- * efficiency. */
-static double absorption_length_acrylic_weighted = 0.0;
-/* D2O Rayleigh scattering length weighted by the Cerenkov spectrum and the
- * quantum efficiency. */
-static double rayleigh_scattering_length_d2o_weighted = 0.0;
-/* H2O Rayleigh scattering length weighted by the Cerenkov spectrum and the
- * quantum efficiency. */
-static double rayleigh_scattering_length_h2o_weighted = 0.0;
+/* Number of points to sample when computing the average absorption and
+ * scattering probabilities.
+ *
+ * Note that we use this array to sample for very small distances when
+ * computing the absorption in the acrylic so (xhi-xlo)/N should be less than
+ * 10 cm. */
+#define N 10000
+
+/* Lower and upper bounds for the distances used to calculate the average
+ * absorption and scattering probabilities. */
+static const double xlo = 0.0;
+static const double xhi = 2*PSUP_RADIUS;
+
+/* Array of lengths used in the lookup tables for the average absorption and
+ * scattering probabilities. */
+static double fabs_x[N];
+
+/* Average fraction of light which is not absorbed by D2O as a function of
+ * length weighted by the Cerenkov spectrum and the quantum efficiency.
+ *
+ * We calculate the fraction not absorbed instead of absorbed */
+static double fabs_d2o[N];
+static double fabs_h2o[N];
+static double fabs_acrylic[N];
+static double fsct_d2o[N];
+static double fsct_h2o[N];
 
 /* From Table 4 in the paper. */
 static double A0 = 0.243905091;
@@ -88,59 +100,329 @@ static double RIND_D2O_C3 = 0.32;
  * From http://pdg.lbl.gov/2018/reviews/rpp2018-rev-phys-constants.pdf. */
 static double HC = 1.239841973976e3;
 
+double get_index(double p, double wavelength, double T)
+{
+    /* Returns the index of refraction of pure water for a given density,
+     * wavelength, and temperature. The density should be in units of g/cm^3,
+     * the wavelength in nm, and the temperature in Celsius.
+     *
+     * See "Refractive Index of Water and Steam as a function of Wavelength,
+     * Temperature, and Density" by Schiebener et al. 1990. */
+    /* normalize the temperature and pressure */
+    wavelength = wavelength/589.0;
+    T = (T+273.15)/273.15;
+
+    /* first we compute the right hand side of Equation 7 */
+    double c = A0 + A1*p + A2*T + A3*pow(wavelength,2)*T + A4/pow(wavelength,2) + A5/(pow(wavelength,2)-pow(UV,2)) + A6/(pow(wavelength,2)-pow(IR,2)) + A7*pow(p,2);
+    c *= p;
+
+    return sqrt((2*c+1)/(1-c));
+}
+
+double get_index_snoman_h2o(double wavelength)
+{
+    /* Returns the index of refraction of light water that SNOMAN uses for a
+     * given wavelength.  The wavelength should be in units of nm.
+     *
+     * The coefficients come from media.dat in SNOMAN version 5.0294 and the
+     * formula used to compute the index of refraction comes from the SNOMAN
+     * companion page for titles MEDA. */
+    double E;
+
+    /* Calculate the energy of the photon in eV. */
+    E = HC/wavelength;
+
+    return RIND_H2O_C1 + RIND_H2O_C2*exp(RIND_H2O_C3*E);
+}
+
+double get_index_snoman_d2o(double wavelength)
+{
+    /* Returns the index of refraction of heavy water that SNOMAN uses for a
+     * given wavelength.  The wavelength should be in units of nm.
+     *
+     * The coefficients come from media.dat in SNOMAN version 5.0294 and the
+     * formula used to compute the index of refraction comes from the SNOMAN
+     * companion page for titles MEDA. */
+    double E;
+
+    /* Calculate the energy of the photon in eV. */
+    E = HC/wavelength;
+
+    return RIND_D2O_C1 + RIND_D2O_C2*exp(RIND_D2O_C3*E);
+}
+
+/* Returns the average probability that a photon is not absorbed after a
+ * distance `x` in cm in D2O.
+ *
+ * This average probability is found by averaging
+ *
+ *     e^(-x/l(wavelength))
+ *
+ * over the wavelength range 200-800 nm weighted by the Cerenkov wavelenght
+ * distribution and the quantum efficiency. */
+double get_fabs_d2o(double x)
+{
+    if (!initialized) {
+        fprintf(stderr, "weighted absorption length hasn't been initialized!\n");
+        exit(1);
+    }
+
+    return interp1d(x,fabs_x,fabs_d2o,LEN(fabs_x));
+}
+
+/* Returns the average probability that a photon is not absorbed after a
+ * distance `x` in cm in H2O. See get_fabs_d2o() for how this is calculated. */
+double get_fabs_h2o(double x)
+{
+    if (!initialized) {
+        fprintf(stderr, "average probability hasn't been initialized!\n");
+        exit(1);
+    }
+
+    return interp1d(x,fabs_x,fabs_h2o,LEN(fabs_x));
+}
+
+/* Returns the average probability that a photon is not absorbed after a
+ * distance `x` in cm in the SNO acrylic. See get_fabs_d2o() for how this is
+ * calculated.
+ *
+ * FIXME: Should include a term for the "dark" acrylic. */
+double get_fabs_acrylic(double x)
+{
+    if (!initialized) {
+        fprintf(stderr, "weighted absorption length hasn't been initialized!\n");
+        exit(1);
+    }
+
+    return interp1d(x,fabs_x,fabs_acrylic,LEN(fabs_x));
+}
+
+/* Returns the average probability that a photon is not scattered after a
+ * distance `x` in cm in D2O. See get_fabs_d2o() for how this is calculated. */
+double get_fsct_d2o(double x)
+{
+    if (!initialized) {
+        fprintf(stderr, "weighted rayleigh scattering length hasn't been initialized!\n");
+        exit(1);
+    }
+
+    return interp1d(x,fabs_x,fsct_d2o,LEN(fabs_x));
+}
+
+/* Returns the average probability that a photon is not scattered after a
+ * distance `x` in cm in H2O. See get_fabs_d2o() for how this is calculated. */
+double get_fsct_h2o(double x)
+{
+    if (!initialized) {
+        fprintf(stderr, "weighted rayleigh scattering length hasn't been initialized!\n");
+        exit(1);
+    }
+
+    return interp1d(x,fabs_x,fsct_h2o,LEN(fabs_x));
+}
+
+/* Returns the absorption length as a function of wavelength in D2O.
+ *
+ * The wavelength should be in nm and the absorption length is in cm. */
+double get_absorption_length_snoman_d2o(double wavelength)
+{
+    /* Note: We use GSL splines here instead of interp1d() since we aren't
+     * guaranteed that the values in the lookup table are evenly spaced. */
+    static gsl_spline *spline;
+    static gsl_interp_accel *acc;
+
+    if (!spline) {
+        spline = gsl_spline_alloc(gsl_interp_linear, LEN(absorption_coefficient_d2o_wavelengths));
+        gsl_spline_init(spline, absorption_coefficient_d2o_wavelengths, absorption_coefficient_d2o, LEN(absorption_coefficient_d2o_wavelengths));
+        acc = gsl_interp_accel_alloc();
+    }
+
+    /* If the wavelength is outside the range of the lookup table we just
+     * return the value at the beginning or end. We do this because by default
+     * gsl_spline_eval will exit if the wavelength is out of bounds. */
+    if (wavelength < absorption_coefficient_d2o_wavelengths[0])
+        wavelength = absorption_coefficient_d2o_wavelengths[0];
+
+    if (wavelength > absorption_coefficient_d2o_wavelengths[LEN(absorption_coefficient_d2o_wavelengths)-1])
+        wavelength = absorption_coefficient_d2o_wavelengths[LEN(absorption_coefficient_d2o_wavelengths)-1];
+
+    return 1.0/gsl_spline_eval(spline, wavelength, acc);
+}
+
+/* Returns the absorption length as a function of wavelength in H2O.
+ *
+ * The wavelength should be in nm and the absorption length is in cm. */
+double get_absorption_length_snoman_h2o(double wavelength)
+{
+    /* Note: We use GSL splines here instead of interp1d() since we aren't
+     * guaranteed that the values in the lookup table are evenly spaced. */
+    static gsl_spline *spline;
+    static gsl_interp_accel *acc;
+
+    if (!spline) {
+        spline = gsl_spline_alloc(gsl_interp_linear, LEN(absorption_coefficient_h2o_wavelengths));
+        gsl_spline_init(spline, absorption_coefficient_h2o_wavelengths, absorption_coefficient_h2o, LEN(absorption_coefficient_h2o_wavelengths));
+        acc = gsl_interp_accel_alloc();
+    }
+
+    /* If the wavelength is outside the range of the lookup table we just
+     * return the value at the beginning or end. We do this because by default
+     * gsl_spline_eval will exit if the wavelength is out of bounds. */
+    if (wavelength < absorption_coefficient_h2o_wavelengths[0])
+        wavelength = absorption_coefficient_h2o_wavelengths[0];
+
+    if (wavelength > absorption_coefficient_h2o_wavelengths[LEN(absorption_coefficient_h2o_wavelengths)-1])
+        wavelength = absorption_coefficient_h2o_wavelengths[LEN(absorption_coefficient_h2o_wavelengths)-1];
+
+    return 1.0/gsl_spline_eval(spline, wavelength, acc);
+}
+
+/* Returns the absorption length as a function of wavelength in acrylic.
+ *
+ * The wavelength should be in nm and the absorption length is in cm. */
+double get_absorption_length_snoman_acrylic(double wavelength)
+{
+    /* Note: We use GSL splines here instead of interp1d() since we aren't
+     * guaranteed that the values in the lookup table are evenly spaced. */
+    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 the wavelength is outside the range of the lookup table we just
+     * return the value at the beginning or end. We do this because by default
+     * gsl_spline_eval will exit if the wavelength is out of bounds. */
+    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);
+}
+
+/* Returns the Rayleigh scattering length at a particular wavelength in D2O.
+ *
+ * The wavelength should be in nm and the returned scattering length is in cm. */
+double get_rayleigh_scattering_length_snoman_d2o(double wavelength)
+{
+    /* Note: We use GSL splines here instead of interp1d() since we aren't
+     * guaranteed that the values in the lookup table are evenly spaced. */
+    static gsl_spline *spline;
+    static gsl_interp_accel *acc;
+
+    if (!spline) {
+        spline = gsl_spline_alloc(gsl_interp_linear, LEN(rayleigh_coefficient_wavelengths));
+        gsl_spline_init(spline, rayleigh_coefficient_wavelengths, rayleigh_coefficient_d2o, LEN(rayleigh_coefficient_wavelengths));
+        acc = gsl_interp_accel_alloc();
+    }
+
+    /* If the wavelength is outside the range of the lookup table we just
+     * return the value at the beginning or end. We do this because by default
+     * gsl_spline_eval will exit if the wavelength is out of bounds. */
+    if (wavelength < rayleigh_coefficient_wavelengths[0])
+        wavelength = rayleigh_coefficient_wavelengths[0];
+
+    if (wavelength > rayleigh_coefficient_wavelengths[LEN(rayleigh_coefficient_wavelengths)-1])
+        wavelength = rayleigh_coefficient_wavelengths[LEN(rayleigh_coefficient_wavelengths)-1];
+
+    return 1.0/gsl_spline_eval(spline, wavelength, acc);
+}
+
+/* Returns the Rayleigh scattering length at a particular wavelength in H2O.
+ *
+ * The wavelength should be in nm and the returned scattering length is in cm. */
+double get_rayleigh_scattering_length_snoman_h2o(double wavelength)
+{
+    /* Returns the Rayleigh scattering length in H2O in cm. */
+    static gsl_spline *spline;
+    static gsl_interp_accel *acc;
+
+    if (!spline) {
+        spline = gsl_spline_alloc(gsl_interp_linear, LEN(rayleigh_coefficient_wavelengths));
+        gsl_spline_init(spline, rayleigh_coefficient_wavelengths, rayleigh_coefficient_h2o, LEN(rayleigh_coefficient_wavelengths));
+        acc = gsl_interp_accel_alloc();
+    }
+
+    /* If the wavelength is outside the range of the lookup table we just
+     * return the value at the beginning or end. We do this because by default
+     * gsl_spline_eval will exit if the wavelength is out of bounds. */
+    if (wavelength < rayleigh_coefficient_wavelengths[0])
+        wavelength = rayleigh_coefficient_wavelengths[0];
+
+    if (wavelength > rayleigh_coefficient_wavelengths[LEN(rayleigh_coefficient_wavelengths)-1])
+        wavelength = rayleigh_coefficient_wavelengths[LEN(rayleigh_coefficient_wavelengths)-1];
+
+    return 1.0/gsl_spline_eval(spline, wavelength, acc);
+}
+
 static double gsl_rayleigh_scattering_length_snoman_d2o(double wavelength, void *params)
 {
     /* Returns the Rayleigh scattering length in D2O weighted by the quantum efficiency
      * and the Cerenkov spectrum. */
-    double qe;
+    double qe, x;
+
+    x = *((double *) params);
 
     qe = get_quantum_efficiency(wavelength);
 
-    return qe*get_rayleigh_scattering_length_snoman_d2o(wavelength)/pow(wavelength,2);
+    return qe*exp(-x/get_rayleigh_scattering_length_snoman_d2o(wavelength))/pow(wavelength,2);
 }
 
 static double gsl_rayleigh_scattering_length_snoman_h2o(double wavelength, void *params)
 {
     /* Returns the Rayleigh scattering length in H2O weighted by the quantum efficiency
      * and the Cerenkov spectrum. */
-    double qe;
+    double qe, x;
+
+    x = *((double *) params);
 
     qe = get_quantum_efficiency(wavelength);
 
-    return qe*get_rayleigh_scattering_length_snoman_h2o(wavelength)/pow(wavelength,2);
+    return qe*exp(-x/get_rayleigh_scattering_length_snoman_h2o(wavelength))/pow(wavelength,2);
 }
 
 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;
+    double qe, x;
+
+    x = *((double *) params);
 
     qe = get_quantum_efficiency(wavelength);
 
-    return qe*get_absorption_length_snoman_d2o(wavelength)/pow(wavelength,2);
+    return qe*exp(-x/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;
+    double qe, x;
+
+    x = *((double *) params);
 
     qe = get_quantum_efficiency(wavelength);
 
-    return qe*get_absorption_length_snoman_h2o(wavelength)/pow(wavelength,2);
+    return qe*exp(-x/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;
+    double qe, x;
+
+    x = *((double *) params);
 
     qe = get_quantum_efficiency(wavelength);
 
-    return qe*get_absorption_length_snoman_acrylic(wavelength)/pow(wavelength,2);
+    return qe*exp(-x/get_absorption_length_snoman_acrylic(wavelength))/pow(wavelength,2);
 }
 
 static double gsl_cerenkov(double wavelength, void *params)
@@ -172,41 +454,57 @@ int optics_init(dict *db)
     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);
+        sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+        return -1;
+    }
+
+    for (i = 0; i < LEN(fabs_x); i++) {
+        fabs_x[i] = xlo + (xhi-xlo)*i/(LEN(fabs_x)-1);
     }
 
     F.function = &gsl_absorption_length_snoman_d2o;
 
-    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+    for (i = 0; i < LEN(fabs_x); i++) {
+        F.params = fabs_x+i;
 
-    absorption_length_d2o_weighted = result/norm;
+        status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
 
-    if (status) {
-        fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
-        exit(1);
+        if (status) {
+            sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+            return -1;
+        }
+
+        fabs_d2o[i] = result/norm;
     }
 
     F.function = &gsl_absorption_length_snoman_h2o;
 
-    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+    for (i = 0; i < LEN(fabs_x); i++) {
+        F.params = fabs_x+i;
 
-    absorption_length_h2o_weighted = result/norm;
+        status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
 
-    if (status) {
-        fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
-        exit(1);
+        if (status) {
+            sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+            return -1;
+        }
+
+        fabs_h2o[i] = result/norm;
     }
 
     F.function = &gsl_absorption_length_snoman_acrylic;
 
-    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+    for (i = 0; i < LEN(fabs_x); i++) {
+        F.params = fabs_x+i;
 
-    absorption_length_acrylic_weighted = result/norm;
+        status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
 
-    if (status) {
-        fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
-        exit(1);
+        if (status) {
+            sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+            return -1;
+        }
+
+        fabs_acrylic[i] = result/norm;
     }
 
     rspr = get_bank(db, "RSPR", 1);
@@ -224,24 +522,32 @@ int optics_init(dict *db)
 
     F.function = &gsl_rayleigh_scattering_length_snoman_d2o;
 
-    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+    for (i = 0; i < LEN(fabs_x); i++) {
+        F.params = fabs_x+i;
 
-    rayleigh_scattering_length_d2o_weighted = result/norm;
+        status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
 
-    if (status) {
-        fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
-        exit(1);
+        if (status) {
+            sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+            return -1;
+        }
+
+        fsct_d2o[i] = result/norm;
     }
 
     F.function = &gsl_rayleigh_scattering_length_snoman_h2o;
 
-    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+    for (i = 0; i < LEN(fabs_x); i++) {
+        F.params = fabs_x+i;
 
-    rayleigh_scattering_length_h2o_weighted = result/norm;
+        status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
 
-    if (status) {
-        fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
-        exit(1);
+        if (status) {
+            sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+            return -1;
+        }
+
+        fsct_h2o[i] = result/norm;
     }
 
     gsl_integration_cquad_workspace_free(w);
@@ -250,189 +556,3 @@ int optics_init(dict *db)
 
     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_rayleigh_scattering_length_snoman_d2o(double wavelength)
-{
-    /* Returns the Rayleigh scattering length in D2O in cm. */
-    return 1.0/interp1d(wavelength,rayleigh_coefficient_wavelengths, rayleigh_coefficient_d2o,LEN(rayleigh_coefficient_wavelengths));
-}
-
-double get_rayleigh_scattering_length_snoman_h2o(double wavelength)
-{
-    /* Returns the Rayleigh scattering length in H2O in cm. */
-    return 1.0/interp1d(wavelength,rayleigh_coefficient_wavelengths, rayleigh_coefficient_h2o,LEN(rayleigh_coefficient_wavelengths));
-}
-
-double get_weighted_rayleigh_scattering_length_snoman_d2o(void)
-{
-    /* Returns the weighted Rayleigh scattering length in D2O in cm. */
-
-    if (!initialized) {
-        fprintf(stderr, "weighted rayleigh scattering length hasn't been initialized!\n");
-        exit(1);
-    }
-
-    return rayleigh_scattering_length_d2o_weighted;
-}
-
-double get_weighted_rayleigh_scattering_length_snoman_h2o(void)
-{
-    /* Returns the weighted Rayleigh scattering length in H2O in cm. */
-
-    if (!initialized) {
-        fprintf(stderr, "weighted rayleigh scattering length hasn't been initialized!\n");
-        exit(1);
-    }
-
-    return rayleigh_scattering_length_h2o_weighted;
-}
-
-double get_absorption_length_snoman_h2o(double wavelength)
-{
-    /* Returns the absorption length in H2O in cm. */
-    static gsl_spline *spline;
-    static gsl_interp_accel *acc;
-
-    if (!spline) {
-        spline = gsl_spline_alloc(gsl_interp_linear, LEN(absorption_coefficient_h2o_wavelengths));
-        gsl_spline_init(spline, absorption_coefficient_h2o_wavelengths, absorption_coefficient_h2o, LEN(absorption_coefficient_h2o_wavelengths));
-        acc = gsl_interp_accel_alloc();
-    }
-
-    if (wavelength < absorption_coefficient_h2o_wavelengths[0])
-        wavelength = absorption_coefficient_h2o_wavelengths[0];
-
-    if (wavelength > absorption_coefficient_h2o_wavelengths[LEN(absorption_coefficient_h2o_wavelengths)-1])
-        wavelength = absorption_coefficient_h2o_wavelengths[LEN(absorption_coefficient_h2o_wavelengths)-1];
-
-    return 1.0/gsl_spline_eval(spline, wavelength, acc);
-}
-
-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. */
-    static gsl_spline *spline;
-    static gsl_interp_accel *acc;
-
-    if (!spline) {
-        spline = gsl_spline_alloc(gsl_interp_linear, LEN(absorption_coefficient_d2o_wavelengths));
-        gsl_spline_init(spline, absorption_coefficient_d2o_wavelengths, absorption_coefficient_d2o, LEN(absorption_coefficient_d2o_wavelengths));
-        acc = gsl_interp_accel_alloc();
-    }
-
-    if (wavelength < absorption_coefficient_d2o_wavelengths[0])
-        wavelength = absorption_coefficient_d2o_wavelengths[0];
-
-    if (wavelength > absorption_coefficient_d2o_wavelengths[LEN(absorption_coefficient_d2o_wavelengths)-1])
-        wavelength = absorption_coefficient_d2o_wavelengths[LEN(absorption_coefficient_d2o_wavelengths)-1];
-
-    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,
-     * wavelength, and temperature. The density should be in units of g/cm^3,
-     * the wavelength in nm, and the temperature in Celsius.
-     *
-     * See "Refractive Index of Water and Steam as a function of Wavelength,
-     * Temperature, and Density" by Schiebener et al. 1990. */
-    /* normalize the temperature and pressure */
-    wavelength = wavelength/589.0;
-    T = (T+273.15)/273.15;
-
-    /* first we compute the right hand side of Equation 7 */
-    double c = A0 + A1*p + A2*T + A3*pow(wavelength,2)*T + A4/pow(wavelength,2) + A5/(pow(wavelength,2)-pow(UV,2)) + A6/(pow(wavelength,2)-pow(IR,2)) + A7*pow(p,2);
-    c *= p;
-
-    return sqrt((2*c+1)/(1-c));
-}
-
-double get_index_snoman_h2o(double wavelength)
-{
-    /* Returns the index of refraction of light water that SNOMAN uses for a
-     * given wavelength.  The wavelength should be in units of nm.
-     *
-     * The coefficients come from media.dat in SNOMAN version 5.0294 and the
-     * formula used to compute the index of refraction comes from the SNOMAN
-     * companion page for titles MEDA. */
-    double E;
-
-    /* Calculate the energy of the photon in eV. */
-    E = HC/wavelength;
-
-    return RIND_H2O_C1 + RIND_H2O_C2*exp(RIND_H2O_C3*E);
-}
-
-double get_index_snoman_d2o(double wavelength)
-{
-    /* Returns the index of refraction of heavy water that SNOMAN uses for a
-     * given wavelength.  The wavelength should be in units of nm.
-     *
-     * The coefficients come from media.dat in SNOMAN version 5.0294 and the
-     * formula used to compute the index of refraction comes from the SNOMAN
-     * companion page for titles MEDA. */
-    double E;
-
-    /* Calculate the energy of the photon in eV. */
-    E = HC/wavelength;
-
-    return RIND_D2O_C1 + RIND_D2O_C2*exp(RIND_D2O_C3*E);
-}