1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
|
#!/usr/bin/env python
# Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago>
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option)
# any later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
# more details.
#
# You should have received a copy of the GNU General Public License along with
# this program. If not, see <https://www.gnu.org/licenses/>.
"""
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()
|
