+
+
+
+
+
+
+
+
+
The dispersion is 31.019 +/- 0.164 milli arc seconds
diff --git a/atmospheric refraction example.ipynb b/atmospheric refraction example.ipynb
index 8da9d5b..59a1783 100644
--- a/atmospheric refraction example.ipynb
+++ b/atmospheric refraction example.ipynb
@@ -31,12 +31,11 @@
"l2 = 1.78 # micron\n",
"\n",
"# Initializing dispersion model\n",
- "at = Ciddor1996(p=P, T=T, H=H, xc=xc)\n",
+ "at = Observatory()\n",
"\n",
"# Calculating indices of refraction for l1 and l2\n",
"n1 = at.n_tph(l=l1, T=T, p=P, RH=H, xc=xc)\n",
- "n2 = at.n_tph(l=l2, T=T, p=P, RH=H, xc=xc)\n",
- "\n"
+ "n2 = at.n_tph(l=l2, T=T, p=P, RH=H, xc=xc)\n"
]
},
{
@@ -56,8 +55,9 @@
"# Density of the atmosphere (following CIPM-81/91 equations)\n",
"rho = at.rho(p=P, T=T, RH=H, xc=xc)\n",
"\n",
- "# Initializing refraction model\n",
- "ref = refraction(lat, h, P, rho)"
+ "# Initializing refraction model and setting the reduced height\n",
+ "ref = refraction(lat, h)\n",
+ "ref.setReducedHeight(P, rho)"
]
},
{
@@ -91,27 +91,74 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "Dispersion is simply the refraction differential between two wavelengths. This is calculated with this module by"
+ "Some error propagation is included in the module as well. It takes is able to calculate the uncertainty in the refractive index and uses that to propagate it into the refraction formule for the cassini model. Note: The Cassini model is the only refraction model for which error propagation is possible, in this module.\n",
+ "\n",
+ "In this example we assume the following uncertainties:\n",
+ "* Wavelength: 1 nm\n",
+ "* Temperature: 0.2 K\n",
+ "* Pressure: 30 Pa\n",
+ "* Relative humidity: 0.03 (3%)\n",
+ "* CO2 density: 25 ppm\n",
+ "* Observed zenith angle: 1 arc second"
]
},
{
"cell_type": "code",
- "execution_count": 4,
+ "execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
- "The dispersion is 31.019 milli arc seconds\n"
+ "The refraction is: 40.697 +/- 0.035 arc seconds\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Calculating the uncertainty in the refractive index\n",
+ "n1 = at.n_tph(l=l1, T=T, p=P, RH=H, xc=xc)\n",
+ "dn1 = at.dn_tph(l1, T, P, H, xc, dl=0.001, dT=0.2, dP=30, dRH=0.03, dCO2=25)\n",
+ "\n",
+ "# Calculating the uncertainty in the atmospheric refraction\n",
+ "refError = ref.cassiniError(n1, dn1, zenith_obs, dz=1/3600)*3.6e3\n",
+ "print ('The refraction is: %.3f +/- %.3f arc seconds' %(refraction_obs, refError))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Dispersion is simply the refraction differential between two wavelengths. This is calculated with a dispersion object. Again, we can also calculate the uncertainty in the dispersion based on the measurement uncertainties. We take the same values as above.\n",
+ "\n",
+ "Note that, since the unertainties in the refractive index at two wavelengths observated at the same time are correlated, we cannot simply use n1, n2, dn1 and dn2 as an input but require the whole set of conditions and uncertainties."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "The dispersion is 31.019 +/- 0.164 milli arc seconds\n"
]
}
],
"source": [
"# Initializing dispersion object\n",
- "disp = dispersion(lat, h, P, rho)\n",
- "dispersion_obs = disp.cassini(n1, n2, zenith=zenith_obs)*3.6e6\n",
- "print ('The dispersion is %.03f milli arc seconds' %(dispersion_obs))"
+ "disp = dispersion(lat, h)\n",
+ "disp.setReducedHeight(P, rho)\n",
+ "\n",
+ "dispersion_obs = disp.cassini(n1, n2, zenith_obs)*3.6e6\n",
+ "\n",
+ "# Calculating the uncertainty in the atmospheric dispersion\n",
+ "dispersion_unc = disp.cassiniError(zenith_obs, l1, l2, T, P, H, xc, \n",
+ " dl1=0.001, dl2=0.001, dT=0.2, dP=30, dRH=0.03, dCO2=25, dz=1/3600) * 3.6e6\n",
+ "\n",
+ "print ('The dispersion is %.03f +/- %.03f milli arc seconds' %(dispersion_obs, dispersion_unc))"
]
}
],
diff --git a/atmosphericModel.py b/atmosphericModel.py
index 3c7cdf5..fcb148b 100644
--- a/atmosphericModel.py
+++ b/atmosphericModel.py
@@ -17,14 +17,13 @@ from scipy import integrate
'''
-class Ciddor1996:
- def __init__(self, p=101325, T=288.15, H=0, xc=450):
- self.p = p # Pressure in Pa
- self.T = T # Temperature in Kelvin
- self.H = H # Relative humidity on a scale from 0 to 1
- self.xc = xc # ppm of CO2
- pass
-
+class Observatory:
+ def __init__(self, dT=0.2, dP=20, dH=0.02, dCO2=20, dl=0.001):
+ self.dT = dT # in K
+ self.dP = dP # in Pa
+ self.dH = dH # relative humidity between 0 and 1
+ self.dCO2 = dCO2 # in ppm
+ self.dl = dl # in micron
def n_as(self, l):
'''
@@ -191,15 +190,216 @@ class Ciddor1996:
return np.sqrt( (1 + 2 * _Ltotal) / (1 - _Ltotal) )
-class error(Ciddor1996):
- def __init__(self, p=101325, T=288.15, H=0, xc=450, dT=0.2, dP=20, dH=0.02, dCO2=20, dl=0.001):
- Ciddor1996.__init__(self, p, T, H)
- self.dT = 0#dT # in K
- self.dP = 0#dP # in Pa
- self.dH = 0#dH # relative humidity between 0 and 1
- self.dCO2 = 0#dCO2 # in ppm
- self.dl = 0#dl # in micron
+ def dn_tph(self, l, T, p, RH, xc, dl=None, dT=None, dP=None, dRH=None, dCO2=None):
+ '''
+ Calculates the uncertainty in the refractive index, following the Ciddor model,
+ by using fully propagated error analysis.
+ This involves almost endless use of the chain rule and evaluation of both the moist
+ air and dry air components. Note: dry air does not contain any water, i.e. H=0, x_w=0, etc.
+
+ Parameters:
+ l - Wavelength in micron (0.4 - 2.5 micron)
+ T - Temperature in Kelvin
+ p - Pressure in Pascals
+ RH - Relative humidity (no units, don't use percentage but value 0