1 files changed, 85 insertions, 0 deletions
diff --git a/utils/index.py b/utils/index.py
new file mode 100644
index 0000000..4edcd59
--- /dev/null
+++ b/utils/index.py
@@ -0,0 +1,85 @@
+#!/usr/bin/env python
+"""
+This is a script to calculate the index of refraction of water as a function of
+wavelength, temperature, and pressure and fit it to a linear approximation in a
+wavelength region where the PMTs are sensitive to. The data comes from [1].
+
+[1] "Refractive index of Water and Steam as Function of Wavelength, Temperature, and Density". Schiebener et al. 1989.
+"""
+from __future__ import print_function, division
+import numpy as np
+
+A0 = 0.243905091
+A1 = 9.53518094e-3
+A2 = -3.64358110e-3
+A3 = 2.65666426e-4
+A4 = 1.59189325e-3
+A5 = 2.45733798e-3
+A6 = 0.897478251
+A7 = -1.63066183e-2
+UV = 0.2292020
+IR = 5.432937
+
+def get_index(p, wavelength, T):
+    """
+    Returns the index of pure water for a given density, wavelength, and
+    temperature. The density should be in units of kg/m^3, the wavelength in nm,
+    and the temperature in Celsius.
+    """
+    # normalize the density, temperature, and pressure
+    p = p/1000.0
+    wavelength = wavelength/589.0
+    T = (T+273.15)/273.15
+
+    # first we compute the right hand side of Equation 7
+    c = A0 + A1*p + A2*T + A3*wavelength**2*T + A4/wavelength**2 + A5/(wavelength**2-UV**2) + A6/(wavelength**2-IR**2) + A7*p**2
+    c *= p
+
+    return np.sqrt((2*c+1)/(1-c))
+
+def approximate_index(x, a, b):
+    """
+    Returns the index of refraction using the approximation:
+
+    n = a + b/x**2
+    """
+    return a + b/x**2
+
+def chi2(pars, x, data):
+    """
+    Returns the squares of the difference between the data and the model for
+    the inverse of the index of refraction, i.e.
+
+    1/n = lambda**2/(a*lambda**2+b)
+    """
+    a, b = pars
+    return np.sum((approximate_index(x,a,b)-data)**2)
+
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    import argparse
+    from scipy.optimize import fmin
+
+    parser = argparse.ArgumentParser("fit for an approximation to the inverse of the index of refraction of water")
+    parser.add_argument("--temperature", type=float, default=10.0)
+    parser.add_argument("--density", type=float, default=1000.0)
+    parser.add_argument("--low", type=float, default=300.0)
+    parser.add_argument("--high", type=float, default=600.0)
+    args = parser.parse_args()
+
+    x = np.linspace(args.low,args.high,1000)
+    n = get_index(args.density,x,args.temperature)
+    xmin = fmin(chi2,[1,1],args=(x,1/n))
+
+    print(xmin)
+
+    a, b = xmin
+    y = approximate_index(x,a,b)
+
+    plt.plot(x,1/n,label="data")
+    plt.plot(x,y,label="fit")
+    plt.xlabel("Wavelength (nm)")
+    plt.ylabel("1/n")
+    plt.title("Index of Refraction of Water at 273.15 K")
+    plt.legend()
+    plt.show()