add refractive index for heavy and light water from snoman

3 files changed, 52 insertions, 4 deletions

diff --git a/src/muon.c b/src/muon.c
index 52c73ed..0279984 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -240,16 +240,17 @@ double get_expected_charge(double x, double T, double *pos, double *dir, double
 
     R = NORM(pos);
 
+    /* FIXME: I just calculate delta assuming 400 nm light. */
+    wavelength0 = 400.0;
+
     if (R <= AV_RADIUS) {
         rho = HEAVY_WATER_DENSITY;
+        n = get_index_snoman_d2o(wavelength0);
     } else {
         rho = WATER_DENSITY;
+        n = get_index_snoman_h2o(wavelength0);
     }
 
-    /* FIXME: I just calculate delta assuming 400 nm light. */
-    wavelength0 = 400.0;
-    n = get_index(rho, wavelength0, 10.0);
-
     if (beta < 1/n) return 0;
 
     /* FIXME: is this formula valid for muons? */
diff --git a/src/optics.c b/src/optics.c
index 23e043d..dd9b82d 100644
--- a/src/optics.c
+++ b/src/optics.c
@@ -13,6 +13,19 @@ static double A7 = -1.63066183e-2;
 static double UV = 0.2292020;
 static double IR = 5.432937;
 
+static double RIND_H2O_C1 = 1.302;
+static double RIND_H2O_C2 = 0.01562;
+static double RIND_H2O_C3 = 0.32;
+
+static double RIND_D2O_C1 = 1.302;
+static double RIND_D2O_C2 = 0.01333;
+static double RIND_D2O_C3 = 0.32;
+
+/* Wavelength of 1 eV particle in nm.
+ *
+ * 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,
@@ -31,3 +44,35 @@ double get_index(double p, double wavelength, double T)
 
     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);
+}
diff --git a/src/optics.h b/src/optics.h
index 5e39ea1..e2df659 100644
--- a/src/optics.h
+++ b/src/optics.h
@@ -2,5 +2,7 @@
 #define OPTICS_H
 
 double get_index(double p, double wavelength, double T);
+double get_index_snoman_h2o(double wavelength);
+double get_index_snoman_d2o(double wavelength);
 
 #endif