Restructured the code

.gitignore

@@ -1,2 +1,3 @@
 .DS_Store
+__pychache__/
 

README.md

@@ -10,16 +10,21 @@ Author info: *Joost van den Born*, *born@astron.nl*
 *Below you can find a set of models describing the atmospheric refraction and dispersion.*
 
 **Contents:**
-- atmosphericModel.py           : Python model that can be used to calculate the atmospheric refraction and dispersion
+- atmosphericModel          	: Python model that can be used to calculate the atmospheric refraction and dispersion, the main content of this repository.
-- refractionModel.py            : Python script that compares the output of the different rerfraction models. 
+- refractionModel.py            : Python script that compares the output of the different refraction models. 
                                   Can be used as an example on how to use atmosphericModel.py
-- 1976USSA.txt                  : Table with temperature, pressure, density and more as a function of altitude.
-                                  Values derived from the 1976 U.S. Standard Atmosphere.
-- manual.pdf                    : Manual for using atmosphericModel.py
 - atmospheric refraction example: jupyter notebook example
-- [Tech note]...models.pdf      : Tech note describing the models.
 
 
+**On the contents of the atmosphericModel package**
+The package contains the following files:
+- 1976USSA.txt 					: Data file of the 1976 US Standard Atmosphere.
+- ciddorModel.py 				: Contains the `Observatory()` class to calculate the refractive index of atmospheric air, using the Ciddor (1996) model.
+- dispersionModels.py 			: Contains various atmospheric models that can be used to calculate the atmospheric dispersion
+- examples.py 					: Contains several quick functions, e.g. to quickly calculate the refractive index.
+- refractionModels.py 			: Contains various atmospheric models that can be used to calculate the atmospheric refraction
+- refractivityModels.py 		: Contains various functions from literature to calculate the refractive index. Some examples of works included are Ciddor (1996), Barrell & Sears (1939), SLALIB, Owens (1967) and Edlen (1966).
+
 
 **On the dispersion model:**
 

atmospheric refraction example.html

@@ -13286,6 +13286,57 @@ div#notebook {
 </div>
 </div>
 
+</div>
+<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
+</div><div class="inner_cell">
+<div class="text_cell_render border-box-sizing rendered_html">
+<p>There are also some quick calculation functions defined in the <code>examples.py</code> file, if you just want to know the refractive index, the refraction or the dispersion at standard atmospheric conditions (<code>conditions='STANDARD'</code>) or at Cerro Armazones (<code>conditions='CERRO_ARMAZONES'</code>).</p>
+
+</div>
+</div>
+</div>
+<div class="cell border-box-sizing code_cell rendered">
+<div class="input">
+<div class="prompt input_prompt">In&nbsp;[6]:</div>
+<div class="inner_cell">
+    <div class="input_area">
+<div class=" highlight hl-ipython3"><pre><span></span><span class="c1"># Quick refractive index</span>
+<span class="n">nq</span> <span class="o">=</span> <span class="n">quick_refractive_index</span><span class="p">(</span><span class="n">l</span><span class="o">=</span><span class="mf">1.49</span><span class="p">,</span> <span class="n">conditions</span><span class="o">=</span><span class="s1">&#39;STANDARD&#39;</span><span class="p">)</span>
+<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Quick refractive index: </span><span class="si">{}</span><span class="s1">&#39;</span><span class="o">.</span><span class="n">format</span><span class="p">(</span><span class="n">nq</span><span class="p">))</span>
+
+<span class="c1"># Quick refraction</span>
+<span class="n">rq</span> <span class="o">=</span> <span class="n">quick_refraction</span><span class="p">(</span><span class="n">l1</span><span class="p">,</span> <span class="n">zenith_obs</span><span class="p">,</span> <span class="n">conditions</span><span class="o">=</span><span class="s1">&#39;CERRO_ARMAZONES&#39;</span><span class="p">)</span> <span class="o">*</span> <span class="mf">3.6e3</span>
+<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Quick refraction: </span><span class="si">{:.3f}</span><span class="s1"> arcsec&#39;</span><span class="o">.</span><span class="n">format</span><span class="p">(</span><span class="n">rq</span><span class="p">))</span>
+
+<span class="c1"># Quick dispersion</span>
+<span class="n">dq</span> <span class="o">=</span> <span class="n">quick_dispersion</span><span class="p">(</span><span class="n">l1</span><span class="p">,</span> <span class="n">l2</span><span class="p">,</span> <span class="n">zenith_obs</span><span class="p">,</span> <span class="n">conditions</span><span class="o">=</span><span class="s1">&#39;CERRO_ARMAZONES&#39;</span><span class="p">)</span> <span class="o">*</span> <span class="mf">3.6e6</span>
+<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Quick dispersion: </span><span class="si">{:.3f}</span><span class="s1"> mas&#39;</span><span class="o">.</span><span class="n">format</span><span class="p">(</span><span class="n">dq</span><span class="p">))</span>
+</pre></div>
+
+    </div>
+</div>
+</div>
+
+<div class="output_wrapper">
+<div class="output">
+
+
+<div class="output_area">
+
+    <div class="prompt"></div>
+
+
+<div class="output_subarea output_stream output_stdout output_text">
+<pre>Quick refractive index: 1.0002733131243147
+Quick refraction: 40.697 arcsec
+Quick dispersion: 31.019 mas
+</pre>
+</div>
+</div>
+
+</div>
+</div>
+
 </div>
     </div>
   </div>

atmospheric refraction example.ipynb

@@ -159,6 +159,42 @@
     "\n",
     "print ('The dispersion is %.03f +/- %.03f milli arc seconds' %(dispersion_obs, dispersion_unc))"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "There are also some quick calculation functions defined in the `examples.py` file, if you just want to know the refractive index, the refraction or the dispersion at standard atmospheric conditions (`conditions='STANDARD'`) or at Cerro Armazones (`conditions='CERRO_ARMAZONES'`)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Quick refractive index: 1.0002733131243147\n",
+      "Quick refraction: 40.697 arcsec\n",
+      "Quick dispersion: 31.019 mas\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Quick refractive index\n",
+    "nq = quick_refractive_index(l=1.49, conditions='STANDARD')\n",
+    "print('Quick refractive index: {}'.format(nq))\n",
+    "\n",
+    "# Quick refraction\n",
+    "rq = quick_refraction(l1, zenith_obs, conditions='CERRO_ARMAZONES') * 3.6e3\n",
+    "print('Quick refraction: {:.3f} arcsec'.format(rq))\n",
+    "\n",
+    "# Quick dispersion\n",
+    "dq = quick_dispersion(l1, l2, zenith_obs, conditions='CERRO_ARMAZONES') * 3.6e6\n",
+    "print('Quick dispersion: {:.3f} mas'.format(dq))\n"
+   ]
   }
  ],
  "metadata": {

atmosphericModel/__init__.py

@@ -0,0 +1,10 @@
+#!/usr/bin/env Python3
+import os
+
+from .ciddorModel import *
+from .dispersionModels import *
+from .refractionModels import *
+from .refractionModels import *
+from .examples import *
+from .refractivityModels import *
+

atmosphericModel/ciddorModel.py

@@ -0,0 +1,385 @@
+#!/usr/bin/env Python3
+import numpy as np
+from scipy import integrate
+
+'''
+@ Author:		Joost van den Born 
+@ Contact:		born@astron.nl
+@ Description:	This script offers a complete model of the atmospheric refraction
+				and dispersion. It uses the Ciddor1996 paper to calculate the 
+				index of refraction for a wide range of conditions, with corresponding
+				error propagation if needed.
+				The actual refraction is calculated by any of the following models:
+				- Plane parallel atmosphere
+				- Cassini spherical homogeneous atmosphere
+				- Other methods described in Corbard et al., 2019.
+'''
+
+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):
+		'''
+			Formula that calculate the refractive index for standard atmospheric
+			conditions, i.e. 15 C, 101,325 Pa, 0% RH, 450 ppm CO2.
+			It returns n_as, not (n_as - 1) * 1e8!
+		'''
+		_s = 1 / l
+		_k0 = 238.0185 	# um^-2
+		_k1 = 5792105  	# um^-2
+		_k2 = 57.362		# um^-2
+		_k3 = 167917		# um^-2
+		return (_k1 / (_k0 - _s**2) + _k3 / (_k2 - _s**2)) / 1e8 + 1
+
+
+	def n_axs(self, l, xc):
+		'''
+			Formula that calculate the refractive index for standard atmospheric
+			conditions, i.e. 15 C, 101,325 Pa, 0% RH, but with xc ppm of CO2.
+			It returns n_axs, not (n_axs - 1)
+		'''
+		_n_as = self.n_as(l)
+		return ( (_n_as - 1) * (1 + 0.534e-6 * (xc - 450) ) ) + 1
+
+
+	def n_ws(self, l):
+		'''
+			Formula that calculates the refractivity of air for standard moist conditions
+			as defined by B&S (1939), i.e. 20 C, 1333 Pa
+			Returns n_ws, not (n_ws - 1) * 1e8
+		'''
+		_s  = 1 / l
+		_cf = 1.022
+		_w0 = 295.235 	# um^-2
+		_w1 = 2.6422	# um^-2
+		_w2 = -0.032380	# um^-4
+		_w3 = 0.004028	# um^-6
+		return (_cf * (_w0 + _w1 * _s**2 + _w2 * _s**4 + _w3 * _s**6) ) / 1e8 + 1
+
+
+	def f(self, p, T):
+		'''
+			Enhancement factor of water vapour in air
+			t in degrees celsius
+			p in Pascal (Pa)
+		'''
+		_t = T - 273.15
+		_a = 1.00062
+		_b = 3.14e-8 	# Pa^-1
+		_c = 5.6e-7		# deg C^-2
+		return _a + _b * p + _c * _t**2
+
+
+	def svp(self, T):
+		'''
+			Saturation vapour pressure in Pa.
+		'''
+		_A =  1.2378847e-5 # K^-2
+		_B = -1.9121316e-2 # K^-2
+		_C = 33.93711047
+		_D = -6.3431645e3  # K^-2
+		return np.exp(_A*T**2 + _B*T + _C + _D/T)
+
+
+	def svp2(self, T):
+		'''
+			Calculation of svp from IAPWS formula, based on documentation by NIST toolbox.
+		'''
+		K1 =  1.16705214528e+03
+		K2 = -7.24213167032e+05
+		K3 = -1.70738469401e+01
+		K4 =  1.20208247025e+04
+		K5 = -3.23255503223e+06
+		K6 =  1.49151086135e+01
+		K7 = -4.82326573616e+03
+		K8 =  4.05113405421e+05
+		K9 = -2.38555575678e-01
+		K10 = 6.50175348448e+02
+		_W = T + K9 / (T - K10)
+		_A = _W**2 + K1*_W + K2
+		_B = K3*_W**2 + K4*_W + K5
+		_C = K6*_W**2 + K7*_W + K8
+		_X = -_B + np.sqrt(_B**2 - 4*_A*_C)
+		return (2*_C / _X)**4 * 1e6
+
+
+	def xw(self, p, T, RH):
+		'''
+			Molar fraction of water vapour in moist air.
+		'''
+		_f 		= self.f(p, T)
+		_svp 	= self.svp(T)
+		return _f * RH * _svp / p
+
+
+	def Z(self, p, T, RH, xw=''):
+		'''
+			Compressibility of moist air
+		'''
+		_a0 =  1.58123e-6 # K Pa^-1
+		_a1 = -2.9331e-8  # Pa^-1
+		_a2 =  1.1043e-10 # K^-1 Pa^-1
+		_b0 =  5.707e-6   # K Pa^-1
+		_b1 = -2.051e-8   # Pa^-1
+		_c0 =  1.9898e-4  # K Pa^-1
+		_c1 = -2.376e-6   # Pa^-1
+		_d  =  1.83e-11   # K^2 Pa^-2
+		_e  = -0.765e-8   # K^2 Pa^-2
+		_t  = T - 273.15
+		if xw == '':
+			_xw = self.xw(p, T, RH)
+		else:
+			_xw = xw
+		Z = 1 - (p/T) * (_a0 + _a1*_t + _a2*_t**2 + (_b0 + _b1*_t)*_xw + (_c0 + _c1*_t)*_xw**2) + (p/T)**2 * (_d + _e*_xw**2)
+		return Z
+
+
+	def rho(self, p, T, RH, xc, xw=''):
+		'''
+			Calculates the density of air for the given parameters.
+			T in Kelvin
+			p in Pascal
+			RH between 0 and 1
+			xc in ppm
+		'''
+		_Ma = 1e-3 * (28.9635 + 12.011e-6 * (xc - 400) )	# Molar mass in kg/mol of dry air containing xc ppm of CO2
+		_Mw = 0.018015 										# kg/mol
+		_R  = 8.314510										# J mol^-1 K^-1, gas constant
+		if xw == '':
+			_xw = self.xw(p, T, RH)
+		else: 
+			_xw = xw
+		_Z  = self.Z(p, T, RH, xw=_xw)
+		return ( ( p * _Ma ) / (_Z * _R * T) ) * ( 1 - _xw * (1 - _Mw/_Ma)) 
+
+
+	def n_tph(self, l, T, p, RH, xc):
+		'''
+			Calculates the refraction of the atmosphere for any conditions.
+		'''
+		_n_a 	 = self.n_axs(l, xc)
+		_n_w 	 = self.n_ws(l)
+		_La 	 = (_n_a**2 - 1) / (_n_a**2 + 2)
+		_Lw 	 = (_n_w**2 - 1) / (_n_w**2 + 2)
+		
+		# Calculate the density of the dry air components, using the partial pressure of the dry air, i.e. p_dry = p_atm - p_water
+		_rho_a 	 = self.rho(p=p-RH*self.svp(T), T=T, RH=0, xc=xc)				# Dry air component for the actual conditions.
+		_rho_axs = self.rho(p=101325, T=288.15, RH=0, xc=xc) 					# Dry air component at standard conditions.
+		
+		if RH == 0:
+			# nan is returned for _rho_w if we put in a (partial) pressure of 0
+			_rho_w = 0
+		else:
+			# Partial pressure of water vapour is p_water = h*svp
+			_rho_w 	 = self.rho(p=RH*self.svp(T), T=T, RH=RH, xc=xc)			# Moist air component for the actual conditions.
+		_rho_ws	 = self.rho(p=1333, T=293.15, RH=1333/self.svp(293.15), xc=xc) 	# Moist air component at standard conditions.
+		
+		_Ltotal  = ( _rho_a / _rho_axs ) * _La + ( _rho_w / _rho_ws ) * _Lw
+
+		# print 'Dry air component = ', (_rho_a/_rho_axs)*(_n_a-1)*1e8
+		# print 'Wet component = ', ( _rho_w / _rho_ws ) * (_n_w - 1)*1e8
+		# print 'Sum of components = ', (_rho_a/_rho_axs)*(_n_a-1)*1e8+( _rho_w / _rho_ws ) * (_n_w - 1)*1e8
+		return np.sqrt( (1 + 2 * _Ltotal) / (1 - _Ltotal) )
+
+
+	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<x<1)
+			xc 	- CO2 density in ppm
+			dl - Uncertainty in the blue wavelength
+			dT 	- Uncertainty in the temperature, default: 0.2 K
+			dP 	- Uncertainty in the pressure, default: 20 Pa
+			dRH - Uncertainty in the relative humidity, default: 0.02
+			dCO2- Uncertainty in the CO2 density, default: 20
+		'''
+
+		if dP == None:
+			dP = self.dP
+		if dT == None:
+			dT = self.dT
+		if dRH == None:
+			dRH = self.dH
+		if dCO2 == None:
+			dCO2 = self.dCO2
+
+		# constants
+		_cf = 1.022
+		_w0 = 295.235 	# um^-2
+		_w1 = 2.6422	# um^-2
+		_w2 = -0.032380	# um^-4
+		_w3 = 0.004028	# um^-6
+
+		_k0 = 238.0185 		# um^-2
+		_k1 = 5792105  		# um^-2
+		_k2 = 57.362		# um^-2
+		_k3 = 167917		# um^-2	
+
+		_A =  1.2378847e-5 # K^-2
+		_B = -1.9121316e-2 # K^-2
+		_C = 33.93711047
+		_D = -6.3431645e3  # K^-2
+
+		_Mw  = 0.018015   # kg/mol
+		_R   = 8.314510	  # J mol^-1 K^-1, gas constant
+
+		_a0 =  1.58123e-6 # K Pa^-1
+		_a1 = -2.9331e-8  # Pa^-1
+		_a2 =  1.1043e-10 # K^-1 Pa^-1
+		_b0 =  5.707e-6   # K Pa^-1
+		_b1 = -2.051e-8   # Pa^-1
+		_c0 =  1.9898e-4  # K Pa^-1
+		_c1 = -2.376e-6   # Pa^-1
+		_d  =  1.83e-11   # K^2 Pa^-2
+		_e  = -0.765e-8   # K^2 Pa^-2
+		_t  = T - 273.15
+
+		_aa = 1.00062
+		_bb = 3.14e-8 	# Pa^-1
+		_cc = 5.6e-7	# deg C^-2
+
+		_RH_a = 0
+		_RH_w = RH
+
+		# Calculated values
+		# _rho_obs = self.rho(p, T, RH, xc)
+		# _ref	 = refraction(lat, h, p, _rho_obs)
+		# _alpha   = _ref.rc * np.sin(np.deg2rad(z)) / (_ref.rc + _ref.H)
+		_p_a 	 = p - RH*self.svp(T) 
+		_p_w     = RH*self.svp(T)
+		_x_w_a	 = self.xw(_p_a, T, _RH_a)
+		_x_w_w 	 = self.f(_p_w, T)
+		_Z_a 	 = self.Z(p=_p_a, T=T, RH=_RH_a)
+		_Z_w 	 = self.Z(p=_p_w, T=T, RH=_RH_w, xw=_x_w_w)
+		_Ma  	 = 1e-3 * (28.9635 + 12.011e-6 * (xc - 400) )	# Molar mass in kg/mol of dry air containing xc ppm of CO2
+		_p_sv 	 = self.svp(T)
+		_rho_a 	 = self.rho(p=p-RH*self.svp(T), T=T, RH=0, xc=xc)
+		_rho_axs = self.rho(p=101325, T=288.15, RH=0, xc=xc)
+		if RH == 0.0:
+			# nan is returned for _rho_w if we put in a (partial) pressure of 0
+			_rho_w = 0
+		else:
+			_rho_w = self.rho(p=RH*self.svp(T), T=T, RH=RH, xc=xc) 
+		_rho_ws  = self.rho(p=1333, T=293.15, RH=1333/self.svp(293.15), xc=xc)
+		_n_a1 = self.n_axs(l, xc)
+		_n_w1 = self.n_ws(l)
+
+		_L1_a = (_n_a1**2 - 1) / (_n_a1**2 + 2)  # Dry air component of L
+		_L1_w = (_n_w1**2 - 1) / (_n_w1**2 + 2)
+		_L1   = ( _rho_a / _rho_axs ) * _L1_a + ( _rho_w / _rho_ws ) * _L1_w # Overall L
+		_n1   = np.sqrt( (1 + 2 * _L1) / (1 - _L1) )
+
+		# Partial derivatives for derivation to l1
+		_dn_wdl  = _cf * 1e-8 * (-1 * _w1 * 2/l**3 - _w2 * 4/l**5 - _w3 * 6/l**7)
+		_dLdn_w1 = _rho_w / _rho_ws * 6 * _n_w1 / (_n_w1**2 + 2)**2
+		_dn_adl  = 1e-8 * (1 + 0.534e-6 * (xc - 450)) * ( - (2 * _k1 * l ) / ((_k0 * l**2 - 1)**2) - (2 * _k3 * l) / ((_k2 * l**2 - 1)**2) )
+		_dLdn_a1 = _rho_a / _rho_axs * 6 * _n_a1 / (_n_a1**2 + 2)**2
+		_dndL1   = ( ( (1 + 2*_L1) / (1 - _L1)**2 ) + (2 / (1-_L1)) ) / (2 * np.sqrt( (1 + 2*_L1) / (1 - _L1) ) ) 
+		_dndl   = _dndL1 * ( _dLdn_a1 * _dn_adl + _dLdn_w1 * _dn_wdl)
+
+		# # Partial derivatives for derivation to l2
+
+
+		# Partial derivatives for derivation to T
+		# Dry air parts:
+		_drho_adp_a  = _Ma / (_Z_a * _R * T) * (1 - _x_w_a * (1 - _Mw / _Ma))
+		_dp_adp_sv 	 = - RH
+		_dp_svdT 	 = (2*_A*T + _B - _D/(T**2)) * np.exp(_A * T**2 + _B * T + _C + _D/T )
+		_drho_adTZ	 = - _p_a * _Ma / (_R * ( _Z_a * T)**2) * (1 - _x_w_a * (1 - _Mw / _Ma))
+		_dTZdT_a	 = _p_a / T * (_a2 * _t**2 + _a1 * _t + _a0 + (_b1*_t + _b0) * _x_w_a + (_c0 + _c1*_t) * _x_w_a**2 ) - _p_a * (_a1 + 2*_a2*_t + _b1*_x_w_a + _c1*_x_w_a**2) - 2*(_p_a/T)**2 * (_d + _e * _x_w_a**2) + _Z_a
+		_dTZdxw_a	 = 2*_e*_p_a**2 / T * _x_w_a - _p_a * (_b1*_t + _b0 + 2 * (_c0 + _c1*_t) * _x_w_a)
+		_dxw_adf	 = _RH_a * _p_sv / _p_a
+		_dfdT 		 = 2*_cc*_t
+		_dxw_adp_sv	 = self.f(_p_a, T) * _RH_a / _p_a
+		_drho_adxw_a = 0 # as rho_a has not water content, i.e. rho = p*M/ZRT
+
+
+		# Moist air parts:
+		_drho_wdp_w  = _Ma / (_Z_w * _R * T) * (1 - _x_w_w * (1 - _Mw / _Ma))
+		_dp_wdp_sv 	 = RH
+		_dp_svdT 	 = (2*_A*T + _B - _D/(T**2)) * np.exp(_A * T**2 + _B * T + _C + _D/T )
+		_drho_wdTZ	 = - _p_w * _Ma / (_R * (_Z_w * T)**2) * (1 - _x_w_w * (1 - _Mw / _Ma))
+		_dTZdT_w	 = _p_w / T * (_a2 * _t**2 + _a1 * _t + _a0 + (_b1*_t + _b0) * _x_w_w + (_c0 + _c1*_t) * _x_w_w**2 ) - _p_w * (_a1 + 2*_a2*_t + _b1*_x_w_w + _c1*_x_w_w**2) - 2*(_p_w/T)**2 * (_d + _e * _x_w_w**2) + _Z_w
+		_dTZdxw_w	 = 2*_e*_p_w**2 / T * _x_w_w - _p_w * (_b1*_t + _b0 + 2 * (_c0 + _c1*_t) * _x_w_w)
+		_dxw_wdf	 = 1
+		_dfdT 		 = 2*_cc*_t
+		_dxw_wdp_sv	 = self.f(_p_w, T) / _p_sv
+		_drho_wdxw_w = _p_w * _Ma / (_Z_w * _R * T) * (_Mw / _Ma - 1)
+
+		_tA =  _drho_adp_a * _dp_adp_sv * _dp_svdT +  _drho_adTZ * (_dTZdT_a + _dTZdxw_a * (_dxw_adf * _dfdT + _dxw_adp_sv * _dp_svdT) ) + _drho_adxw_a * ( _dxw_adf * _dfdT + _dxw_adp_sv * _dp_svdT )
+		_tB =  _drho_wdp_w * _dp_wdp_sv * _dp_svdT +  _drho_wdTZ * (_dTZdT_w + _dTZdxw_w * (_dxw_wdf * _dfdT + _dxw_wdp_sv * _dp_svdT) ) + _drho_wdxw_w * ( _dxw_wdf * _dfdT + _dxw_wdp_sv * _dp_svdT )
+		_dL1drho_a 	 = _L1_a / _rho_axs
+		_dL1drho_w 	 = _L1_w / _rho_ws
+		
+		_dndT = _tA * (_dndL1 * _dL1drho_a) + _tB * (_dndL1 * _dL1drho_w)
+
+		# Partial derivatives for derivation to p
+		# Dry air parts:
+		_drho_adZ 	= - _p_a * _Ma / ( _Z_a**2 * _R * T) * (1 - _x_w_a * (1 - _Mw / _Ma))
+		_dZdp_a 	= - 1 / T * (_a0 + _a1 * _t + _a2 * _t**2 + (_b0 + _b1*_t)*_x_w_a + (_c0 + _c1*_t)*_x_w_a**2 ) + 2*_p_a/(T**2) * (_d + _e * _x_w_a**2)
+		_dfdp 		= _bb
+		_dxw_adp_a 	= - self.f(_p_a, T) * _RH_a * _p_sv / (_p_a)**2 
+
+		# Moist air parts
+		_dZdp_w 	= - 1 / T * (_a0 + _a1 * _t + _a2 * _t**2 + (_b0 + _b1*_t)*_x_w_w + (_c0 + _c1*_t)*_x_w_w**2 ) + 2*_p_w/(T**2) * (_d + _e * _x_w_w**2)
+		# Not necessary, since dpwdp=0
+
+		_pA = _drho_adp_a + _drho_adZ * _dZdp_a + _drho_adxw_a * (_dxw_adf * _dfdp + _dxw_adp_a)
+		_dndP = _pA * (_dndL1 * _dL1drho_a)
+
+		
+		# Partial derivatives for derivation to H
+		# Dry air parts:
+		_dp_adH 	= - _p_sv
+		_dZdxw_a 	= (_p_a / T)**2 * 2 * _e * _x_w_a - (_p_a / T) * ( (_b0 + _b1 * _t) + 2 * (_c0 + _c1 * _t) * _x_w_a )
+		_dxw_adH 	= self.f(_p_a, T) * _p_sv / _p_a
+
+		# Moist air parts:
+		_dp_wdH 	= _p_sv
+		_dZdxw_w 	= (_p_w / T)**2 * 2 * _e * _x_w_w - (_p_w / T) * ( (_b0 + _b1 * _t) + 2 * (_c0 + _c1 * _t) * _x_w_w )
+		_dxw_wdH 	= 0
+		_dxw_wdp_w 	= 0		
+		_drho_wdZ 	= - _p_w * _Ma / ( _Z_w**2 * _R * T) * (1 - _x_w_w * (1 - _Mw / _Ma))
+
+		_hA = _drho_adp_a * _dp_adH + _drho_adZ * ( _dZdp_a * _dp_adH + _dZdxw_a * (_dxw_adH + _dxw_adp_a * _dp_adH + _dxw_adf * _dfdp * _dp_adH) ) + _drho_adxw_a * (_dxw_adH + _dxw_adp_a * _dp_adH + _dxw_adf * _dfdp * _dp_adH)
+		_hB = _drho_wdp_w * _dp_wdH + _drho_wdZ * ( _dZdp_w * _dp_wdH + _dZdxw_w * ( _bb * _p_sv ) ) + _drho_wdxw_w * ( _bb * _p_sv )
+		_dndRH = _hA * (_dndL1 * _dL1drho_a) + _hB * (_dndL1 * _dL1drho_w)
+
+
+		# Partial derivatives for derivation to xc
+		# Dry air parts:
+		_drho_adMa   = _p_a * (1 - _x_w_a) / (_Z_a * _R * T)
+		_dL1drho_axs = - (_rho_a * _L1_a)/(_rho_axs**2)
+		_drho_axsdMa = _drho_adMa
+		_dL1dL_a	 = _rho_a / _rho_axs
+		_dL_adn_a1	 = 6 * _n_a1 / (_n_a1**2 + 2)**2
+
+		# Moist air parts:
+		_drho_wdMa   = _p_w * (1 - _x_w_w) / (_Z_w * _R * T)
+		_dL1drho_ws = - (_rho_w * _L1_w)/(_rho_ws**2)
+		_drho_wsdMa = _drho_wdMa
+
+		# Other parts:
+		_dMadx_c   = 12.011e-9
+		_dn_adx_c1 = (self.n_as(l) - 1) * 0.534e-6
+
+		_dL1dx_c = _dL1drho_a * _drho_adMa * _dMadx_c + _dL1drho_w * _drho_wdMa * _dMadx_c + _dL1drho_axs * _drho_axsdMa * _dMadx_c + _dL1drho_ws * _drho_wsdMa * _dMadx_c + _dL1dL_a * _dL_adn_a1 * _dn_adx_c1
+		_dndCO2 = _dndL1 * _dL1dx_c
+
+		# Final error:
+		_nError = np.sqrt( (_dndl * dl)**2 + (_dndT * dT)**2  + (_dndP * dP)**2 + (_dndRH * dRH)**2 + (_dndCO2 * dCO2)**2 )
+
+		return _nError	

atmosphericModel/dispersionModels.py

@@ -0,0 +1,448 @@
+#!/usr/bin/env Python3
+import numpy as np
+from scipy import integrate
+from .ciddorModel import Observatory
+'''
+@ Author:		Joost van den Born 
+@ Contact:		born@astron.nl
+@ Description:	This script offers a complete model of the atmospheric refraction
+				and dispersion. It uses the Ciddor1996 paper to calculate the 
+				index of refraction for a wide range of conditions, with corresponding
+				error propagation if needed.
+				The actual refraction is calculated by any of the following models:
+				- Plane parallel atmosphere
+				- Cassini spherical homogeneous atmosphere
+				- Other methods described in Corbard et al., 2019.
+'''
+
+class dispersion(Observatory):
+	def __init__(self, lat, h):
+		Observatory.__init__(self)
+		self.h  	= h 						# height above sea level in meters
+		self.lat 	= lat						# Latitude of the observer in degrees
+		self.rc 	= self.set_rc() * 1000		# Radius of curvature of the earth at the observer
+
+	def setReducedHeight(self, P, rho):
+		self.H  	= self.set_H(P, rho) 		# Reduced height of the atmosphere assuming ideal gas law
+
+
+	def planeParallel(self, n1, n2, zenith):
+		'''
+			Refraction based on a flat atmosphere, easily derived from Snell's law.
+		'''
+		_R1 = np.arcsin(n1 * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
+		_R2 = np.arcsin(n2 * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
+		return np.degrees(_R1) - np.degrees(_R2)
+
+
+
+	def cassini(self, n1, n2, zenith, reduced_height=None):
+		'''
+			Refraction of spherical atmosphere, derived from geometry using Sine law
+		'''
+		if reduced_height:
+			H = reduced_height
+		else:
+			H = self.H
+		
+		_R1 = np.arcsin(n1 * self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H)) - np.arcsin(self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H) )
+		_R2 = np.arcsin(n2 * self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H)) - np.arcsin(self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H) )	
+		return np.degrees(_R1) - np.degrees(_R2)
+
+
+
+	def tan5(self, n1, n2, zenith, reduced_height=None):
+		'''
+			Same as oriani's theorem, but including a higher order term.
+			Found in Corbard et al., 2019.
+		'''
+		if reduced_height:
+			H = reduced_height
+		else:
+			H = self.H
+		
+		_rList = []
+		for _n in [n1, n2]:
+			_a = _n - 1
+			_b = H / self.rc
+			_R = _a * (1 - _b) * np.tan( np.deg2rad(zenith) )  - _a * (_b - _a/2) * np.tan( np.deg2rad(zenith) )**3 + 3 * _a * (_b - _a/2)**2 * np.tan( np.deg2rad(zenith) )**5
+			_rList.append(np.degrees(_R))
+		return _rList[0] - _rList[1]
+
+
+
+	def oriani(self, n1, n2, zenith, reduced_height=None):
+		'''
+			Classic refraction formula with two tan(z) terms. This does not assume any 
+			information about the structure of the atmosphere. Found with Oriani's theorem.
+		'''
+		if reduced_height:
+			H = reduced_height
+		else:
+			H = self.H
+
+		_rList = []
+		for _n in [n1, n2]:
+			_a = _n - 1
+			_b = H / self.rc
+			_R = _a * (1 - _b) * np.tan( np.deg2rad(zenith) ) - _a * (_b - _a/2) * np.tan( np.deg2rad(zenith) )**3
+			_rList.append(np.degrees(_R))
+		return _rList[0] - _rList[1]
+
+
+
+	def corbard(self, n1, n2, zenith, reduced_height=None):
+		'''
+			Corbard et al., 2019, mention an additional formula based on the error function.
+			Oriani's theorem can be derived from this equation by 'keeping only the three first
+			terms of its asymptotic expansion'.
+		'''
+		if reduced_height:
+			H = reduced_height
+		else:
+			H = self.H
+
+		_rList = []
+		for _n in [n1, n2]:
+			_a = _n - 1
+			_b = H / self.rc
+			_R = _a * ( (2 - _a) / (np.sqrt( 2*_b - _a )) ) * np.sin( np.deg2rad(zenith) ) * self.psi( (np.cos(np.deg2rad(zenith))) / np.sqrt(2*_b - _a) )
+			_rList.append(np.degrees(_R))
+		return _rList[0] - _rList[1]
+
+
+
+	def refractionIntegral(self, n1, n2, z0, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None):
+		'''
+			Using data that gives temperature, pressure and density as a function of height, 
+			the gladstone-dale relation is invoked to calculate the path of a light ray 
+			along the atmosphere. This is allows us to also include atmospheric activity in 
+			different atmosphere layers.
+			Recommended use with the US Standard Atmosphere from 1976.
+			USSA1976 data generated from code supplied by http://www.pdas.com/atmos.html
+			The method described in eq. 3 in Auer & Standish (2000) is used for the calculation.
+		'''
+		if not R0:
+			R0 = self.rc
+		if useStandardAtmosphere:
+			_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('1976USSA.txt', unpack=True, skip_header=6)
+		else:
+			_Dens = rhoData
+			_h 	  = heightData 
+
+		
+		_rList 					= []
+		z0 						= np.deg2rad(z0)
+		
+		# Integration of the refraction integral over the atmosphere layers.
+		for _n in [n1, n2]:
+			_Dens 	= np.asarray(_Dens)
+			_h 	 	= np.asarray(_h) 
+			_R 		= R0 + _h * 1000
+			
+			# Gladstone-Dale relation between refractive index and density
+			_n_h = 1 + (_n - 1) / _Dens[0] * _Dens 
+
+			# Calculation of the zenith angle for given height and refractive index, from the refractive invariant
+			_z 		= np.arcsin( (_n * R0) / (_n_h * _R) * np.sin(z0) )
+			
+			# Calculation of log space derivative of n and r
+			_lnn = np.log(_n_h)
+			_lnr = np.log(_R)
+			_dlnn = np.asarray([_lnn[i+1] - _lnn[i] for i in range(len(_lnn)-1)])
+			_dlnr = np.asarray([_lnr[i+1] - _lnr[i] for i in range(len(_lnr)-1)])
+			_dz   = np.asarray([(_z[i] - _z[i+1]) for i in range(len(_z)-1)])
+
+			# Calculation of the integrand in eq. 3 in Auer & Standish
+			_Integrand = -1*(_dlnn/_dlnr) / (1 + _dlnn/_dlnr) * _dz
+			_rList.append(np.degrees(np.sum(_Integrand)))
+		return _rList[0] - _rList[1]
+
+
+
+	def psi(self, x):
+		_f   = lambda a: np.exp(-a**2)
+		_int = integrate.quad(_f, x, np.inf)
+		return np.exp(x**2) * _int[0]
+
+
+	def set_rc(self):
+		'''
+			Returns the radius of the curvature at the observer in km,
+			taking into account the latitude of the observer and the corresponding
+			ellipsoidality of the earth.
+		'''
+		_A   = 6378.137 # km
+		_B   = 6356.752 # km
+		_phi = np.deg2rad(self.lat)
+		
+		_rc0 = (_A*_B)**2 / (_A**2 * np.cos(_phi)**2 + _B**2 * np.sin(_phi)**2)**(3/2)
+		return _rc0 + self.h/1000
+
+
+	def set_H(self, P, rho):
+		'''
+			Calculates the reduced height of the atmosphere, assuming ideal gas law.
+		'''
+		_phi = np.deg2rad(self.lat)
+		_c1 = 5.2790414e-3
+		_c2 = 2.32718e-5
+		_c3 = 1.262e-7
+		_c4 = 7e-10
+		_g0 = 9.780327 # m/s^2
+		_g0_local = _g0 * (1 + _c1 * np.sin(_phi)**2 + _c2 * np.sin(_phi)**4 + _c3 * np.sin(_phi)**6 + _c4 * np.sin(_phi)**8)
+		_g = _g0_local - (3.0877e-6 - 4.3e-9 * np.sin(_phi)**2) * self.h + 7.2e-13 * self.h**2
+		return P / (rho * _g)
+
+	def H_isotherm(self, T):
+		'''
+			Calculates the reduced height of the atmosphere, assuming an isotherm distribution.
+		'''
+		_kb = 1.38064852e-23 # kg m2 s-2 K-1
+		_m  = 4.809651698e-26 # kg (weight per molecule of air, assuming dry air)
+
+		_phi = np.deg2rad(self.lat)
+		_c1 = 5.2790414e-3
+		_c2 = 2.32718e-5
+		_c3 = 1.262e-7
+		_c4 = 7e-10
+		_g0 = 9.780327 # m/s^2
+		_g0_local = _g0 * (1 + _c1 * np.sin(_phi)**2 + _c2 * np.sin(_phi)**4 + _c3 * np.sin(_phi)**6 + _c4 * np.sin(_phi)**8)
+
+		return (_kb*T) / (_m * _g0_local)
+
+	def cassiniError(self, z, l1, l2, T, p, RH, xc, dl1=0, dl2=0, dT=0.2, dP=20, dRH=0.02, dCO2=20, dz=0, lat=None, h=None):
+		'''
+			Calculates the uncertainty in the atmospheric dispersion, following the Cassini 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.
+
+			Not taken into account is the errror propagation in the scale height of the atmosphere or 
+			the zenith error, from my results this should give a small or even negligible difference.
+
+			Parameters:
+			z 	- Zenith angle in degrees
+			l1 	- Bluest wavelength in micron (0.4 - 2.5 micron)
+			l2 	- Reddest 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<x<1)
+			xc 	- CO2 density in ppm
+			dl1 - Uncertainty in the blue wavelength
+			dl2 - Uncertainty in the red wavelength
+			dT 	- Uncertainty in the temperature, default: 0.2 K
+			dP 	- Uncertainty in the pressure, default: 20 Pa
+			dRH - Uncertainty in the relative humidity, default: 0.02
+			dCO2- Uncertainty in the CO2 density, default: 20
+			lat - [Optional] Latitude of the observatory in degrees
+			h 	- [Optional] Height in m of the observatory
+		'''
+		if lat == None:
+			lat = self.lat
+		else:
+			self.lat = lat
+		if h == None:
+			h 	= self.h
+		else:
+			self.h = h
+
+		# constants
+		_cf = 1.022
+		_w0 = 295.235 	# um^-2
+		_w1 = 2.6422	# um^-2
+		_w2 = -0.032380	# um^-4
+		_w3 = 0.004028	# um^-6
+
+		_k0 = 238.0185 		# um^-2
+		_k1 = 5792105  		# um^-2
+		_k2 = 57.362		# um^-2
+		_k3 = 167917		# um^-2	
+
+		_A =  1.2378847e-5 # K^-2
+		_B = -1.9121316e-2 # K^-2
+		_C = 33.93711047
+		_D = -6.3431645e3  # K^-2
+
+		_Mw  = 0.018015   # kg/mol
+		_R   = 8.314510	  # J mol^-1 K^-1, gas constant
+
+		_a0 =  1.58123e-6 # K Pa^-1
+		_a1 = -2.9331e-8  # Pa^-1
+		_a2 =  1.1043e-10 # K^-1 Pa^-1
+		_b0 =  5.707e-6   # K Pa^-1
+		_b1 = -2.051e-8   # Pa^-1
+		_c0 =  1.9898e-4  # K Pa^-1
+		_c1 = -2.376e-6   # Pa^-1
+		_d  =  1.83e-11   # K^2 Pa^-2
+		_e  = -0.765e-8   # K^2 Pa^-2
+		_t  = T - 273.15
+
+		_aa = 1.00062
+		_bb = 3.14e-8 	# Pa^-1
+		_cc = 5.6e-7	# deg C^-2
+
+		_RH_a = 0
+		_RH_w = RH
+
+		# Calculated values
+		_rho_obs = self.rho(p, T, RH, xc)
+		self.rc  = self.set_rc() * 1000
+		self.H   = self.set_H(p, _rho_obs) 
+		_alpha   = self.rc * np.sin(np.deg2rad(z)) / (self.rc + self.H)
+		_p_a 	 = p - RH*self.svp(T) 
+		_p_w     = RH*self.svp(T)
+		_x_w_a	 = self.xw(_p_a, T, _RH_a)
+		_x_w_w 	 = self.f(_p_w, T)
+		_Z_a 	 = self.Z(p=_p_a, T=T, RH=_RH_a)
+		_Z_w 	 = self.Z(p=_p_w, T=T, RH=_RH_w, xw=_x_w_w)
+		_Ma  	 = 1e-3 * (28.9635 + 12.011e-6 * (xc - 400) )	# Molar mass in kg/mol of dry air containing xc ppm of CO2
+		_p_sv 	 = self.svp(T)
+		_rho_a 	 = self.rho(p=p-RH*self.svp(T), T=T, RH=0, xc=xc)
+		_rho_axs = self.rho(p=101325, T=288.15, RH=0, xc=xc)
+		if RH == 0.0:
+			# nan is returned for _rho_w if we put in a (partial) pressure of 0
+			_rho_w = 0
+		else:
+			_rho_w = self.rho(p=RH*self.svp(T), T=T, RH=RH, xc=xc) 
+		_rho_ws  = self.rho(p=1333, T=293.15, RH=1333/self.svp(293.15), xc=xc)
+		_n_a1 	 = self.n_axs(l1, xc)
+		_n_w1 	 = self.n_ws(l1)
+		_n_a2 	 = self.n_axs(l2, xc)
+		_n_w2 	 = self.n_ws(l2)
+
+		_L1_a 	 = (_n_a1**2 - 1) / (_n_a1**2 + 2)  # Dry air component of L
+		_L1_w 	 = (_n_w1**2 - 1) / (_n_w1**2 + 2)
+		_L1   	 = ( _rho_a / _rho_axs ) * _L1_a + ( _rho_w / _rho_ws ) * _L1_w # Overall L
+		_L2_a 	 = (_n_a2**2 - 1) / (_n_a2**2 + 2) # Dry air component of L
+		_L2_w 	 = (_n_w2**2 - 1) / (_n_w2**2 + 2)
+		_L2    	 = ( _rho_a / _rho_axs ) * _L2_a + ( _rho_w / _rho_ws ) * _L2_w # Overall L
+
+		_n1   	 = np.sqrt( (1 + 2 * _L1) / (1 - _L1) )
+		_n2   	 = np.sqrt( (1 + 2 * _L2) / (1 - _L2) )
+
+		# Partial derivatives for derivation to l1
+		_dn_wdl1 = _cf * 1e-8 * (-1 * _w1 * 2/l1**3 - _w2 * 4/l1**5 - _w3 * 6/l1**7)
+		_dLdn_w1 = _rho_w / _rho_ws * 6 * _n_w1 / (_n_w1**2 + 2)**2
+		_dn_adl1 = 1e-8 * (1 + 0.534e-6 * (xc - 450)) * ( - (2 * _k1 * l1 ) / ((_k0 * l1**2 - 1)**2) - (2 * _k3 * l1) / ((_k2 * l1**2 - 1)**2) )
+		_dLdn_a1 = _rho_a / _rho_axs * 6 * _n_a1 / (_n_a1**2 + 2)**2
+		_dndL1   = ( ( (1 + 2*_L1) / (1 - _L1)**2 ) + (2 / (1-_L1)) ) / (2 * np.sqrt( (1 + 2*_L1) / (1 - _L1) ) ) 
+		_dddn1   = _alpha / np.sqrt( 1 - _alpha**2 * _n1**2 )
+		_dddl1   = _dddn1 * _dndL1 * ( _dLdn_a1 * _dn_adl1 + _dLdn_w1 * _dn_wdl1)
+
+		# Partial derivatives for derivation to l2
+		_dn_wdl2 = _cf * 1e-8 * (-1 * _w1 * 2/l2**3 - _w2 * 4/l2**5 - _w3 * 6/l2**7)
+		_dLdn_w2 = _rho_w / _rho_ws * 6 * _n_w2 / (_n_w2**2 + 2)**2
+		_dn_adl2 = 1e-8 * (1 + 0.534e-6 * (xc - 450)) * ( - (2 * _k1 * l2 ) / ((_k0 * l2**2 - 1)**2) - (2 * _k3 * l2) / ((_k2 * l2**2 - 1)**2) )
+		_dLdn_a2 = _rho_a / _rho_axs * 6 * _n_a2 / (_n_a2**2 + 2)**2
+		_dndL2   = ( ( (1 + 2*_L2) / (1 - _L2)**2 ) + (2 / (1-_L2)) ) / (2 * np.sqrt( (1 + 2*_L2) / (1 - _L2) ) ) 
+		_dddn2   = _alpha / np.sqrt( 1 - _alpha**2 * _n2**2 )
+		_dddl2   = _dddn2 * _dndL2 * ( _dLdn_a2 * _dn_adl2 + _dLdn_w2 * _dn_wdl2)
+
+
+		# Partial derivatives for derivation to T
+		# Dry air parts:
+		_drho_adp_a  = _Ma / (_Z_a * _R * T) * (1 - _x_w_a * (1 - _Mw / _Ma))
+		_dp_adp_sv 	 = - RH
+		_dp_svdT 	 = (2*_A*T + _B - _D/(T**2)) * np.exp(_A * T**2 + _B * T + _C + _D/T )
+		_drho_adTZ	 = - _p_a * _Ma / (_R * ( _Z_a * T)**2) * (1 - _x_w_a * (1 - _Mw / _Ma))
+		_dTZdT_a	 = _p_a / T * (_a2 * _t**2 + _a1 * _t + _a0 + (_b1*_t + _b0) * _x_w_a + (_c0 + _c1*_t) * _x_w_a**2 ) - _p_a * (_a1 + 2*_a2*_t + _b1*_x_w_a + _c1*_x_w_a**2) - 2*(_p_a/T)**2 * (_d + _e * _x_w_a**2) + _Z_a
+		_dTZdxw_a	 = 2*_e*_p_a**2 / T * _x_w_a - _p_a * (_b1*_t + _b0 + 2 * (_c0 + _c1*_t) * _x_w_a)
+		_dxw_adf	 = _RH_a * _p_sv / _p_a
+		_dfdT 		 = 2*_cc*_t
+		_dxw_adp_sv	 = self.f(_p_a, T) * _RH_a / _p_a
+		_drho_adxw_a = 0 # as rho_a has not water content, i.e. rho = p*M/ZRT
+
+
+		# Moist air parts:
+		_drho_wdp_w  = _Ma / (_Z_w * _R * T) * (1 - _x_w_w * (1 - _Mw / _Ma))
+		_dp_wdp_sv 	 = RH
+		_dp_svdT 	 = (2*_A*T + _B - _D/(T**2)) * np.exp(_A * T**2 + _B * T + _C + _D/T )
+		_drho_wdTZ	 = - _p_w * _Ma / (_R * (_Z_w * T)**2) * (1 - _x_w_w * (1 - _Mw / _Ma))
+		_dTZdT_w	 = _p_w / T * (_a2 * _t**2 + _a1 * _t + _a0 + (_b1*_t + _b0) * _x_w_w + (_c0 + _c1*_t) * _x_w_w**2 ) - _p_w * (_a1 + 2*_a2*_t + _b1*_x_w_w + _c1*_x_w_w**2) - 2*(_p_w/T)**2 * (_d + _e * _x_w_w**2) + _Z_w
+		_dTZdxw_w	 = 2*_e*_p_w**2 / T * _x_w_w - _p_w * (_b1*_t + _b0 + 2 * (_c0 + _c1*_t) * _x_w_w)
+		_dxw_wdf	 = 1 # _RH_w * _p_sv / _p_w
+		_dfdT 		 = 2*_cc*_t
+		_dxw_wdp_sv	 = self.f(_p_w, T) / _p_sv # * _RH_w / _p_w
+		_drho_wdxw_w = _p_w * _Ma / (_Z_w * _R * T) * (_Mw / _Ma - 1)
+
+		_tA =  _drho_adp_a * _dp_adp_sv * _dp_svdT +  _drho_adTZ * (_dTZdT_a + _dTZdxw_a * (_dxw_adf * _dfdT + _dxw_adp_sv * _dp_svdT) ) + _drho_adxw_a * ( _dxw_adf * _dfdT + _dxw_adp_sv * _dp_svdT )
+		_tB =  _drho_wdp_w * _dp_wdp_sv * _dp_svdT +  _drho_wdTZ * (_dTZdT_w + _dTZdxw_w * (_dxw_wdf * _dfdT + _dxw_wdp_sv * _dp_svdT) ) + _drho_wdxw_w * ( _dxw_wdf * _dfdT + _dxw_wdp_sv * _dp_svdT )
+		_dL1drho_a 	 = _L1_a / _rho_axs
+		_dL2drho_a 	 = _L2_a / _rho_axs
+		_dL1drho_w 	 = _L1_w / _rho_ws
+		_dL2drho_w 	 = _L2_w / _rho_ws
+		
+		_dddT = _tA * (_dddn1 * _dndL1 * _dL1drho_a - _dddn2 * _dndL2 * _dL2drho_a) + _tB * (_dddn1 * _dndL1 * _dL1drho_w - _dddn2 * _dndL2 * _dL2drho_w)
+
+		# Partial derivatives for derivation to p
+		# Dry air parts:
+		_drho_adZ 	= - _p_a * _Ma / ( _Z_a**2 * _R * T) * (1 - _x_w_a * (1 - _Mw / _Ma))
+		_dZdp_a 	= - 1 / T * (_a0 + _a1 * _t + _a2 * _t**2 + (_b0 + _b1*_t)*_x_w_a + (_c0 + _c1*_t)*_x_w_a**2 ) + 2*_p_a/(T**2) * (_d + _e * _x_w_a**2)
+		_dfdp 		= _bb
+		_dxw_adp_a 	= - self.f(_p_a, T) * _RH_a * _p_sv / (_p_a)**2 
+
+		# Moist air parts
+		_dZdp_w 	= - 1 / T * (_a0 + _a1 * _t + _a2 * _t**2 + (_b0 + _b1*_t)*_x_w_w + (_c0 + _c1*_t)*_x_w_w**2 ) + 2*_p_w/(T**2) * (_d + _e * _x_w_w**2)
+		# Not necessary, since dpwdp=0
+
+		_pA = _drho_adp_a + _drho_adZ * _dZdp_a + _drho_adxw_a * (_dxw_adf * _dfdp + _dxw_adp_a)
+		_dddP = _pA * (_dddn1 * _dndL1 * _dL1drho_a - _dddn2 * _dndL2 * _dL2drho_a)
+
+		
+		# Partial derivatives for derivation to H
+		# Dry air parts:
+		_dp_adH 	= - _p_sv
+		_dZdxw_a 	= (_p_a / T)**2 * 2 * _e * _x_w_a - (_p_a / T) * ( (_b0 + _b1 * _t) + 2 * (_c0 + _c1 * _t) * _x_w_a )
+		_dxw_adH 	= self.f(_p_a, T) * _p_sv / _p_a
+
+		# Moist air parts:
+		_dp_wdH 	= _p_sv
+		_dZdxw_w 	= (_p_w / T)**2 * 2 * _e * _x_w_w - (_p_w / T) * ( (_b0 + _b1 * _t) + 2 * (_c0 + _c1 * _t) * _x_w_w )
+		_dxw_wdH 	= 0
+		_dxw_wdp_w 	= 0		
+		_drho_wdZ 	= - _p_w * _Ma / ( _Z_w**2 * _R * T) * (1 - _x_w_w * (1 - _Mw / _Ma))
+
+		_hA = _drho_adp_a * _dp_adH + _drho_adZ * ( _dZdp_a * _dp_adH + _dZdxw_a * (_dxw_adH + _dxw_adp_a * _dp_adH + _dxw_adf * _dfdp * _dp_adH) ) + _drho_adxw_a * (_dxw_adH + _dxw_adp_a * _dp_adH + _dxw_adf * _dfdp * _dp_adH)
+		_hB = _drho_wdp_w * _dp_wdH + _drho_wdZ * ( _dZdp_w * _dp_wdH + _dZdxw_w * ( _bb * _p_sv ) ) + _drho_wdxw_w * ( _bb * _p_sv )
+		_dddRH = _hA * (_dddn1 * _dndL1 * _dL1drho_a - _dddn2 * _dndL2 * _dL2drho_a) + _hB * (_dddn1 * _dndL1 * _dL1drho_w - _dddn2 * _dndL2 * _dL2drho_w)
+
+
+		# Partial derivatives for derivation to xc
+		# Dry air parts:
+		_drho_adMa   = _p_a * (1 - _x_w_a) / (_Z_a * _R * T)
+		_dL1drho_axs = - (_rho_a * _L1_a)/(_rho_axs**2)
+		_dL2drho_axs = - (_rho_a * _L2_a)/(_rho_axs**2)
+		_drho_axsdMa = _drho_adMa
+		_dL1dL_a	 = _rho_a / _rho_axs
+		_dL_adn_a1	 = 6 * _n_a1 / (_n_a1**2 + 2)**2
+		_dL2dL_a	 = _dL1dL_a
+		_dL_adn_a2	 = 6 * _n_a2 / (_n_a2**2 + 2)**2
+
+		# Moist air parts:
+		_drho_wdMa   = _p_w * (1 - _x_w_w) / (_Z_w * _R * T)
+		_dL1drho_ws = - (_rho_w * _L1_w)/(_rho_ws**2)
+		_dL2drho_ws = - (_rho_w * _L2_w)/(_rho_ws**2)
+		_drho_wsdMa = _drho_wdMa
+
+		# Other parts:
+		_dMadx_c   = 12.011e-9
+		_dn_adx_c1 = (self.n_as(l1) - 1) * 0.534e-6
+		_dn_adx_c2 = (self.n_as(l2) - 1) * 0.534e-6
+
+		_dL1dx_c = _dL1drho_a * _drho_adMa * _dMadx_c + _dL1drho_w * _drho_wdMa * _dMadx_c + _dL1drho_axs * _drho_axsdMa * _dMadx_c + _dL1drho_ws * _drho_wsdMa * _dMadx_c + _dL1dL_a * _dL_adn_a1 * _dn_adx_c1
+		_dL2dx_c = _dL2drho_a * _drho_adMa * _dMadx_c + _dL2drho_w * _drho_wdMa * _dMadx_c + _dL2drho_axs * _drho_axsdMa * _dMadx_c + _dL2drho_ws * _drho_wsdMa * _dMadx_c + _dL2dL_a * _dL_adn_a2 * _dn_adx_c2
+		_dddCO2 = _dddn1 * _dndL1 * _dL1dx_c - _dddn2 * _dndL2 * _dL2dx_c
+
+
+		# Partial derivatives for derivation to z
+		_zenith = np.deg2rad(z)
+		_dz   	= np.deg2rad(dz)
+		_R 		= self.rc
+		_H 		= self.H
+		_dddz 	= (_R * np.cos(_zenith)) / (_R + _H) * ( _n1 / np.sqrt(1-(_n1**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2) - _n2 / np.sqrt(1-(_n2**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2) )		
+
+		# Final error:
+		_dispersionError = np.sqrt( (_dddl1 * dl1)**2 + (_dddl2 * dl2)**2 + (_dddT * dT)**2  + (_dddP * dP)**2 + (_dddRH * dRH)**2 + (_dddCO2 * dCO2)**2  + (_dddz * _dz)**2)
+
+		return np.degrees(_dispersionError)

atmosphericModel/examples.py

@@ -0,0 +1,106 @@
+#!/usr/bin/env Python3
+from .ciddorModel import Observatory
+from .refractionModels import refraction
+from .dispersionModels import dispersion
+'''
+	Author: Joost van den Born
+	Some functions for quick testing.
+'''
+
+
+def quick_refractive_index(l, conditions='STANDARD'):
+	if conditions == 'STANDARD':
+		# Parameters at Cerro Armazones
+		T   = 288.15    # k
+		P   = 101325    # Pa
+		H   = 0.0
+		xc  = 450       # ppm
+		lat = 0  		# degrees
+		h   = 0     	# m
+		
+
+	elif conditions == 'CERRO_ARMAZONES':
+		# Parameters at Cerro Armazones
+		T   = 279.65    # k
+		P   = 71200     # Pa
+		H   = 0.22
+		xc  = 450       # ppm
+		lat = -24.5983  # degrees
+		h   = 3064      # m
+		
+
+	# Initializing dispersion model
+	at  = Observatory()
+
+	# Calculating indices of refraction for l
+	n 	= at.n_tph(l=l, T=T, p=P, RH=H, xc=xc)
+	return n
+
+def quick_refraction(l, z, conditions='STANDARD'):
+	if conditions == 'STANDARD':
+		# Parameters at Cerro Armazones
+		T   = 288.15    # k
+		P   = 101325    # Pa
+		H   = 0.0
+		xc  = 450       # ppm
+		lat = 0  		# degrees
+		h   = 0     	# m
+		
+
+	elif conditions == 'CERRO_ARMAZONES':
+		# Parameters at Cerro Armazones
+		T   = 279.65    # k
+		P   = 71200     # Pa
+		H   = 0.22
+		xc  = 450       # ppm
+		lat = -24.5983  # degrees
+		h   = 3064      # m
+
+	# Initializing dispersion model
+	at  = Observatory()
+
+	# Calculating indices of refraction for l
+	n 	= at.n_tph(l=l, T=T, p=P, RH=H, xc=xc)
+
+	# Density of the atmosphere (following CIPM-81/91 equations)
+	rho = at.rho(p=P, T=T, RH=H, xc=xc)
+
+	# Initializing refraction model and setting the reduced height
+	ref = refraction(lat, h)
+	ref.setReducedHeight(P, rho)
+	return ref.cassini(n, z)
+
+def quick_dispersion(l1, l2, z, conditions='STANDARD'):
+	if conditions == 'STANDARD':
+		# Parameters at Cerro Armazones
+		T   = 288.15    # k
+		P   = 101325    # Pa
+		H   = 0.0
+		xc  = 450       # ppm
+		lat = 0  		# degrees
+		h   = 0     	# m
+		
+
+	elif conditions == 'CERRO_ARMAZONES':
+		# Parameters at Cerro Armazones
+		T   = 279.65    # k
+		P   = 71200     # Pa
+		H   = 0.22
+		xc  = 450       # ppm
+		lat = -24.5983  # degrees
+		h   = 3064      # m
+
+	# Initializing dispersion model
+	at  = Observatory()
+
+	# Calculating indices of refraction for l1 and l2
+	n1 	= at.n_tph(l=l1, T=T, p=P, RH=H, xc=xc)
+	n2 	= at.n_tph(l=l2, T=T, p=P, RH=H, xc=xc)
+
+	# Density of the atmosphere (following CIPM-81/91 equations)
+	rho = at.rho(p=P, T=T, RH=H, xc=xc)
+
+	# Initializing refraction model and setting the reduced height
+	disp = dispersion(lat, h)
+	disp.setReducedHeight(P, rho)
+	return disp.cassini(n1, n2, z)

atmosphericModel/refractionModels.py

@@ -0,0 +1,218 @@
+#!/usr/bin/env Python3
+import numpy as np
+from scipy import integrate
+from .ciddorModel import Observatory
+'''
+@ Author:		Joost van den Born 
+@ Contact:		born@astron.nl
+@ Description:	This script offers a complete model of the atmospheric refraction
+				and dispersion. It uses the Ciddor1996 paper to calculate the 
+				index of refraction for a wide range of conditions, with corresponding
+				error propagation if needed.
+				The actual refraction is calculated by any of the following models:
+				- Plane parallel atmosphere
+				- Cassini spherical homogeneous atmosphere
+				- Other methods described in Corbard et al., 2019.
+'''
+
+class refraction(Observatory):
+	def __init__(self, lat, h):
+		Observatory.__init__(self)
+		self.h  	= h 					# height above sea level in meters
+		self.lat 	= lat					# Latitude of the observer in degrees
+		self.rc 	= self.set_rc() * 1000	# Radius of curvature of the earth at the observer
+		# self.H2 	= self.H_isotherm(T) 	# Reduced height of isotherm atmosphere
+
+	def setReducedHeight(self, P, rho):
+		self.H  	= self.set_H(P, rho) 		# Reduced height of the atmosphere assuming ideal gas law
+
+
+	def planeParallel(self, n, zenith):
+		'''
+			Refraction based on a flat atmosphere, easily derived from Snell's law.
+		'''
+		_R = np.arcsin(n * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
+		return np.degrees(_R)
+
+
+
+	def cassini(self, n1, zenith, reduced_height=None):
+		'''
+			Refraction of spherical atmosphere, derived from geometry using Sine law
+		'''
+		if reduced_height:
+			H = reduced_height
+		else:
+			H = self.H
+		
+		_R1 = np.arcsin(n1 * self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H)) - np.arcsin(self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H) )
+		return np.degrees(_R1)
+
+
+
+	def tan5(self, n1, zenith, reduced_height=None):
+		'''
+			Same as oriani's theorem, but including a higher order term.
+			Found in Corbard et al., 2019.
+		'''
+		if reduced_height:
+			H = reduced_height
+		else:
+			H = self.H
+		
+		_a = n1 - 1
+		_b = H / self.rc
+		_R = _a * (1 - _b) * np.tan( np.deg2rad(zenith) )  - _a * (_b - _a/2) * np.tan( np.deg2rad(zenith) )**3 + 3 * _a * (_b - _a/2)**2 * np.tan( np.deg2rad(zenith) )**5
+		return np.degrees(_R)
+
+
+
+	def oriani(self, n1, zenith, reduced_height=None):
+		'''
+			Classic refraction formula with two tan(z) terms. This does not assume any 
+			information about the structure of the atmosphere. Found with Oriani's theorem.
+		'''
+		if reduced_height:
+			H = reduced_height
+		else:
+			H = self.H
+
+		_a = n1 - 1
+		_b = H / self.rc
+		_R = _a * (1 - _b) * np.tan( np.deg2rad(zenith) ) - _a * (_b - _a/2) * np.tan( np.deg2rad(zenith) )**3
+		return np.degrees(_R)
+
+
+
+	def corbard(self, n1, zenith, reduced_height=None):
+		'''
+			Corbard et al., 2019, mention an additional formula based on the error function.
+			Oriani's theorem can be derived from this equation by 'keeping only the three first
+			terms of its asymptotic expansion'.
+		'''
+		if reduced_height:
+			H = reduced_height
+		else:
+			H = self.H
+
+		_a = n1 - 1
+		_b = H / self.rc
+		_R = _a * ( (2 - _a) / (np.sqrt( 2*_b - _a )) ) * np.sin( np.deg2rad(zenith) ) * self.psi( (np.cos(np.deg2rad(zenith))) / np.sqrt(2*_b - _a) )
+		return np.degrees(_R)
+
+
+
+	def refractionIntegral(self, n0, z0, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None):
+		'''
+			Using data that gives temperature, pressure and density as a function of height, 
+			the gladstone-dale relation is invoked to calculate the path of a light ray 
+			along the atmosphere. This is allows us to also include atmospheric activity in 
+			different atmosphere layers.
+			Recommended use with the US Standard Atmosphere from 1976.
+			USSA1976 data generated from code supplied by http://www.pdas.com/atmos.html
+		'''
+		if not R0:
+			R0 = self.rc
+		if useStandardAtmosphere:
+			_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('1976USSA.txt', unpack=True, skip_header=6)
+		else:
+			_Dens = rhoData
+			_h 	  = heightData
+
+		z0 						= np.deg2rad(z0)
+		_Dens 	= np.asarray(_Dens)
+		_h 	 	= np.asarray(_h) 
+		_R 		= R0 + _h * 1000
+		
+		# Gladstone-Dale relation between refractive index and density
+		_n_h = 1 + (n0 - 1) / _Dens[0] * _Dens 
+
+		# Calculation of the zenith angle for given height and refractive index, from the refractive invariant
+		_z 		= np.arcsin( (n0 * R0) / (_n_h * _R) * np.sin(z0) )
+		
+		# Calculation of log space derivative of n and r
+		_lnn = np.log(_n_h)
+		_lnr = np.log(_R)
+		_dlnn = np.asarray([_lnn[i+1] - _lnn[i] for i in range(len(_lnn)-1)])
+		_dlnr = np.asarray([_lnr[i+1] - _lnr[i] for i in range(len(_lnr)-1)])
+		_dz   = np.asarray([(_z[i] - _z[i+1]) for i in range(len(_z)-1)])
+
+		# Calculation of the integrand in eq. 3 in Auer & Standish
+		_Integrand = -1*(_dlnn/_dlnr) / (1 + _dlnn/_dlnr) * _dz
+
+
+		return np.degrees(np.sum(_Integrand))
+
+
+
+	def psi(self, x):
+		_f   = lambda a: np.exp(-a**2)
+		_int = integrate.quad(_f, x, np.inf)
+		return np.exp(x**2) * _int[0]
+
+
+
+	def set_rc(self):
+		'''
+			Returns the radius of the curvature at the observer in km,
+			taking into account the latitude of the observer and the corresponding
+			ellipsoidality of the earth.
+		'''
+		_A   = 6378.137 # km
+		_B   = 6356.752 # km
+		_phi = np.deg2rad(self.lat)
+		
+		_rc0 = (_A*_B)**2 / (_A**2 * np.cos(_phi)**2 + _B**2 * np.sin(_phi)**2)**(3/2)
+		return _rc0 + self.h/1000
+
+
+
+	def set_H(self, P, rho):
+		'''
+			Calculates the reduced height of the atmosphere, assuming ideal gas law.
+		'''
+		_phi = np.deg2rad(self.lat)
+		_c1 = 5.2790414e-3
+		_c2 = 2.32718e-5
+		_c3 = 1.262e-7
+		_c4 = 7e-10
+		_g0 = 9.780327 # m/s^2
+		_g0_local = _g0 * (1 + _c1 * np.sin(_phi)**2 + _c2 * np.sin(_phi)**4 + _c3 * np.sin(_phi)**6 + _c4 * np.sin(_phi)**8)
+		_g = _g0_local - (3.0877e-6 - 4.3e-9 * np.sin(_phi)**2) * self.h + 7.2e-13 * self.h**2
+		return P / (rho * _g)
+
+	def H_isotherm(self, T):
+		'''
+			Calculates the reduced height of the atmosphere, assuming an isotherm distribution.
+		'''
+		_kb = 1.38064852e-23 # kg m2 s-2 K-1
+		_m  = 4.809651698e-26 # kg (weight per molecule of air, assuming dry air)
+
+		_phi = np.deg2rad(self.lat)
+		_c1 = 5.2790414e-3
+		_c2 = 2.32718e-5
+		_c3 = 1.262e-7
+		_c4 = 7e-10
+		_g0 = 9.780327 # m/s^2
+		_g0_local = _g0 * (1 + _c1 * np.sin(_phi)**2 + _c2 * np.sin(_phi)**4 + _c3 * np.sin(_phi)**6 + _c4 * np.sin(_phi)**8)
+
+		return (_kb*T) / (_m * _g0_local)
+
+	def cassiniError(self, n, dn, zenith, dz=0, reduced_height=None):
+		'''
+		Calculates the uncertainty in the refraction angle in degrees.
+		'''
+		if reduced_height:
+			_H = reduced_height
+		else:
+			_H = self.H
+		_zenith = np.deg2rad(zenith)
+		_dz = np.deg2rad(dz)
+		_R = self.rc
+		_dRdn = (_R * np.sin(_zenith)) / ((_R + _H) * np.sqrt(1 - ((n**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2)))
+
+		_dRdz = (_R * np.cos(_zenith)) / (_R + _H) * ( n / np.sqrt(1-(n**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2) - 1 / np.sqrt(1-(_R**2 * np.sin(_zenith)**2)/(_R + _H)**2) )
+
+		_dR = np.sqrt( (_dRdn * dn)**2 + (_dRdz * _dz)**2)
+
+		return np.degrees(_dR)

atmosphericModel/refractivityModels.py

@@ -0,0 +1,349 @@
+#!/usr/bin/env Python3
+import numpy as np
+
+
+'''
+@ Author:		Joost van den Born 
+@ Contact:		born@astron.nl
+@ Description:	A comparison of various atmospheric refraction models.
+				Also these are used to calculate the difference in 
+				atmospheric dispersion, given a set of conditions.
+'''
+
+def BarrellAndSears(l, T=288.15, p=101325, H=0):
+	'''
+		From Barrell & Sears, 1939
+		Equation 7.7, page 52
+		For atmospheric air (i.e. CO2 = 300 ppm)
+		This equation is applicable for the  ranges T=10-30C and p=720-800 mmHg.
+
+		Input: 
+		l - Wavelength in microns
+		T - Temperature in Kelvin
+		p - Atmospheric pressure in Pa
+		H - fractional humidity (0=<H=<1)
+
+		Returns:
+		n - Refractive index of atmospheric air at the given conditions
+	'''
+	T = K2C(T)
+	P = hPA2mmHg(p/100)
+	R = (0.378125 + 0.0021414 / l**2 + 0.00001793 / l**4) * (P*(1 + (1.049 - 0.0157*T) * P * 1e-6)/(1 + 0.003661*T)) - (0.0624 - 0.000680 / l**2) * H / (1 + 0.003661*T)
+	n = R/1e6 + 1
+	return n
+
+def Ciddor(l):
+	'''
+		From Ciddor, 1996
+		Equations 1
+		At default conditions only. For more flexibility use
+		Observatory.n_tph()
+
+		Input: 
+		l - Wavelength in microns
+
+		Returns:
+		n - Refractive index of atmospheric air at the given conditions
+	'''
+	s = 1 / l
+	R = 5792105 / (238.0185 - s**2) + 167917 / (57.362 - s**2)
+	n = R/1e8 + 1
+	return n
+
+def Edlen1953(l):
+	'''
+		From Edlen, 1953
+		Equation 2
+
+		At standard conditions only
+
+		Input: 
+		l - Wavelength in microns
+
+		Returns:
+		n - Refractive index of atmospheric air at the given conditions
+	'''
+	s = 1 / l
+	R = 6432.8 + 2949810 / (146 - s**2) + 25540/(41-s**2)
+	n = R/1e8 + 1
+	return n
+
+def Edlen1966(l):
+	'''
+		From Edlen, 1966
+		Equation 1
+
+		At standard conditions only
+
+		Input: 
+		l - Wavelength in microns
+
+		Returns:
+		n - Refractive index of atmospheric air at the given conditions
+	'''
+	s = 1 / l
+	R = 8342.13 + 2406030 / (130 - s**2) + 15997 / (38.9 - s**2)
+	n = R/1e8 + 1
+	return n
+
+def Owens(l, T=288.15, p=101325, H=0):
+	'''
+		From Owens, 1967
+		Equation 32
+		Note that the pressure P and H are in torricelli (Torr),
+		however, this is almost precisely 1 Torr = 1 mmHg
+
+		Input:
+		l - wavelength in micron
+		T - Temperature in Kelvin
+		p - Atmospheric pressure in Pa
+		H - Humdity (0<=H<=1)
+
+		Returns:
+		n - Refractive index of atmospheric air at the given conditions
+	'''
+	P = hPA2mmHg(p/100)
+	T = K2C(T)
+	s = 1 / l
+	R = 8342.13 + 2406030 / (130 - s**2) + 15997 / (38.9 - s**2) * (P/720.775) * ((1 + P * (0.817 - 0.0133) * 1e-6)/(1 + 0.0036610*T)) - H * (5.722 - 0.0457*s**2)
+	n = R/1e8 + 1
+	return n
+
+def BonschPotulski(l, T=288.15, p=101325, CO2=450):
+	'''
+		From Bonsch & Potulski, 1998
+		Equations 6, 7, 8
+		
+		Input: 
+		l - Wavelength in microns
+		T - Temperature in Kelvin
+		p - Atmospheric pressure in Pa
+		CO2 - Particles per million of CO2
+
+		Returns:
+		n - Refractive index of atmospheric air at the given conditions
+	'''
+	x = CO2 / 1e6
+	T = K2C(T)
+	s = 1 / l
+	# Rs = 8092.33 + 2333983 / (130 - s**2) + 15518 / (38.9 - s**2)
+	Rs = 8091.37 + 2333983 / (130 - s**2) + 15518 / (38.9 - s**2)
+	Rs = Rs / 1e8
+	Rx = Rs * (1 + 0.5327 * (x - 0.0004))
+	Rtp = Rx * p / 93214.60 * (1 + 1e-8 * (0.5953 - 0.009876 * T)*p) / (1 + 0.0036610 * T)
+	n = Rtp + 1
+	return n
+
+def BirchDowns(l):
+	'''
+		From Birch & Downs, 1993 and 1994 (correction)
+		Equations 1 and 2
+
+		Only at standard conditions
+
+		Input:
+		l - wavelength in micron
+
+		Returns:
+		n - Refractive index of atmospheric air at the given conditions
+	'''
+	s = 1 / l
+	R = 8342.54 + 2406147 / (130 - s**2) + 15998 / (38.9-s**2)
+	n = R/1e8 + 1
+	return n
+
+def PeckReeder(l):
+	'''
+		From Peck & Reeder, 1972
+		Equation 3
+
+		Only at standard conditions
+
+		Input:
+		l - wavelength in micron
+
+		Returns:
+		n - Refractive index of atmospheric air at the given conditions
+	'''
+	s = 1 / l
+	R = 8060.51 + 2480990 / (132.274 - s**2) + 17455.7 / (39.32957 - s**2)
+	n = R/1e8 + 1
+	return n
+
+
+def slalib(l, T=288.15, p=101325, H=0, h=0, lat=-24, tlr=0.0065, eps=1e-10):
+	'''
+		SLALIB uses the same atmospheric model as ZEMAX
+		based on the works of Hohenkerk & Sinclair and Seidelmann.
+		They take the B&S refraction formula as a basis.
+		Requires the pySLALIB package to work.
+
+		Input: 
+		l - Wavelength in microns
+		T - Temperature in Kelvin
+		p - Atmospheric pressure in Pa
+		H - fractional humidity (0=<H=<1)
+		h - Altitude of the observer (default=0 m)
+		lat - Latitude of the observer (default=-24)
+		tlr - temperature lapse rate (default=0.0065)
+		eps - precision at which the calculation is aborted (default=1e-10)
+
+		Returns:
+		n - Refractive index of atmospheric air at the given conditions
+	'''
+	from pyslalib import slalib
+	zObs = np.deg2rad(45.06563)
+	p = p/100
+	R = []
+	if type(l) == float:
+		l = [l]
+	for wl in l:
+		R.append(slalib.sla_refro(zObs, h, T, p, H, wl, lat, tlr, eps))
+	if len(l) == 1:
+		return R[0]+1
+	else:
+		return np.asarray(R)+1
+
+
+def Mathar(l, T=288.15, p=101325, H=0):
+	'''
+		Refractivity model by Mathar (2007), valid above 1.4 microns
+		Code taken and modified from: 
+		https://github.com/polyanskiy/refractiveindex.info-scripts/blob/master/scripts/Mathar%202007%20-%20Air%201.3-2.5um.py
+
+		Input: 
+		l - Wavelength in microns
+		T - Temperature in Kelvin
+		p - Atmospheric pressure in Pa
+		H - fractional humidity (0=<H=<1)
+
+		Returns:
+		n - Refractive index of atmospheric air at the given conditions
+	'''
+	if l < 1.3 or l > 24:
+		print('Using the Mathar model at this wavelength (l={:.3f}) is not recommended!'.format(l))
+
+	if l <= 2.65:
+		########################################################################
+		# model parameters for wavelengths between 1.3 and 2.5 microns
+		cref = [ 0.200192e-3,  0.113474e-9,  -0.424595e-14,  0.100957e-16, -0.293315e-20,  0.307228e-24] # cm^j
+		cT   = [ 0.588625e-1, -0.385766e-7,   0.888019e-10, -0.567650e-13,  0.166615e-16, -0.174845e-20] # cm^j · K
+		cTT  = [-3.01513,      0.406167e-3,  -0.514544e-6,   0.343161e-9,  -0.101189e-12,  0.106749e-16] # cm^j · K^2
+		cH   = [-0.103945e-7,  0.136858e-11, -0.171039e-14,  0.112908e-17, -0.329925e-21,  0.344747e-25] # cm^j · %^-1
+		cHH  = [ 0.573256e-12, 0.186367e-16, -0.228150e-19,  0.150947e-22, -0.441214e-26,  0.461209e-30] # cm^j · %^-2
+		cp   = [ 0.267085e-8,  0.135941e-14,  0.135295e-18,  0.818218e-23, -0.222957e-26,  0.249964e-30] # cm^j · Pa^-1
+		cpp  = [ 0.609186e-17, 0.519024e-23, -0.419477e-27,  0.434120e-30, -0.122445e-33,  0.134816e-37] # cm^j · Pa^-2
+		cTH  = [ 0.497859e-4, -0.661752e-8,   0.832034e-11, -0.551793e-14,  0.161899e-17, -0.169901e-21] # cm^j · K · %^-1
+		cTp  = [ 0.779176e-6,  0.396499e-12,  0.395114e-16,  0.233587e-20, -0.636441e-24,  0.716868e-28] # cm^j · K · Pa^-1
+		cHp  = [-0.206567e-15, 0.106141e-20, -0.149982e-23,  0.984046e-27, -0.288266e-30,  0.299105e-34] # cm^j · %^-1 · Pa^-1
+
+		sref = 1e4/2.25    # cm^−1, reference wavenumber
+		Tref = 273.15+17.5 # K, reference temperature
+		pref = 75000       # Pa, reference pressure
+		Href = 10          #%, reference humdity in percent
+		########################################################################
+	
+	if l > 2.65 and l <= 4.275:
+		########################################################################
+		# model parameters for wavelengths between 2.8 and 4.2 microns
+		cref = [ 0.200049e-3,  0.145221e-9,   0.250951e-12, -0.745834e-15, -0.161432e-17,  0.352780e-20] # cm^j
+		cT   = [ 0.588432e-1, -0.825182e-7,   0.137982e-9,   0.352420e-13, -0.730651e-15, -0.167911e-18] # cm^j · K
+		cTT  = [-3.13579,      0.694124e-3,  -0.500604e-6,  -0.116668e-8,   0.209644e-11,  0.591037e-14] # cm^j · K^2
+		cH   = [-0.108142e-7,  0.230102e-11, -0.154652e-14, -0.323014e-17,  0.630616e-20,  0.173880e-22] # cm^j · %^-1
+		cHH  = [ 0.586812e-12, 0.312198e-16, -0.197792e-19, -0.461945e-22,  0.788398e-25,  0.245580e-27] # cm^j · %^-2
+		cp   = [ 0.266900e-8,  0.168162e-14,  0.353075e-17, -0.963455e-20, -0.223079e-22,  0.453166e-25] # cm^j · Pa^-1
+		cpp  = [ 0.608860e-17, 0.461560e-22,  0.184282e-24, -0.524471e-27, -0.121299e-29,  0.246512e-32] # cm^j · Pa^-2
+		cTH  = [ 0.517962e-4, -0.112149e-7,   0.776507e-11,  0.172569e-13, -0.320582e-16, -0.899435e-19] # cm^j · K · %^-1
+		cTp  = [ 0.778638e-6,  0.446396e-12,  0.784600e-15, -0.195151e-17, -0.542083e-20,  0.103530e-22] # cm^j · K · Pa^-1
+		cHp  = [-0.217243e-15, 0.104747e-20, -0.523689e-23,  0.817386e-26,  0.309913e-28, -0.363491e-31] # cm^j · %^-1 · Pa^-1
+
+		sref = 1e4/3.4     # cm^−1, reference wavenumber
+		Tref = 273.15+17.5 # K, reference temperature
+		pref = 75000       # Pa, reference pressure
+		Href = 10          #%, reference humdity in percent
+		########################################################################
+
+	if l > 4.275 and l <= 6.35:
+		########################################################################
+		# model parameters for wavelengths between 4.35 and 5.2 microns
+		cref = [ 0.200020e-3,  0.275346e-9,   0.325702e-12, -0.693603e-14,  0.285610e-17,  0.338758e-18] # cm^j
+		cT   = [ 0.590035e-1, -0.375764e-6,   0.134585e-9,   0.124316e-11,  0.508510e-13, -0.189245e-15] # cm^j · K
+		cTT  = [-4.09830,      0.250037e-2,   0.275187e-6,  -0.653398e-8,  -0.310589e-9,   0.127747e-11] # cm^j · K^2
+		cH   = [-0.140463e-7,  0.839350e-11, -0.190929e-14, -0.121399e-16, -0.898863e-18,  0.364662e-20] # cm^j · %^-1
+		cHH  = [ 0.543605e-12, 0.112802e-15, -0.229979e-19, -0.191450e-21, -0.120352e-22,  0.500955e-25] # cm^j · %^-2
+		cp   = [ 0.266898e-8,  0.273629e-14,  0.463466e-17, -0.916894e-23,  0.136685e-21,  0.413687e-23] # cm^j · Pa^-1
+		cpp  = [ 0.610706e-17, 0.116620e-21,  0.244736e-24, -0.497682e-26,  0.742024e-29,  0.224625e-30] # cm^j · Pa^-2
+		cTH  = [ 0.674488e-4, -0.406775e-7,   0.289063e-11,  0.819898e-13,  0.468386e-14, -0.191182e-16] # cm^j · K · %^-1
+		cTp  = [ 0.778627e-6,  0.593296e-12,  0.145042e-14,  0.489815e-17,  0.327941e-19,  0.128020e-21] # cm^j · K · Pa^-1
+		cHp  = [-0.211676e-15, 0.487921e-20, -0.682545e-23,  0.942802e-25, -0.946422e-27, -0.153682e-29] # cm^j · %^-1 · Pa^-1
+		σref = 1e4/4.8     # cm^−1, reference wavenumber
+		Tref = 273.15+17.5 # K, reference temperature
+		pref = 75000       # Pa, reference pressure
+		Href = 10          #%, reference humdity in percent
+		########################################################################
+
+	if l > 6.35 and l <= 15.05:
+		########################################################################
+		# model parameters for wavelengths between 7.5 and 14.1 microns
+		cref = [ 0.199885e-3,  0.344739e-9,  -0.273714e-12,  0.393383e-15, -0.569488e-17,  0.164556e-19] # cm^j
+		cT   = [ 0.593900e-1, -0.172226e-5,   0.237654e-8,  -0.381812e-11,  0.305050e-14, -0.157464e-16] # cm^j · K
+		cTT  = [-6.50355,      0.103830e-1,  -0.139464e-4,   0.220077e-7,  -0.272412e-10,  0.126364e-12] # cm^j · K^2
+		cH   = [-0.221938e-7,  0.347377e-10, -0.465991e-13,  0.735848e-16, -0.897119e-19,  0.380817e-21] # cm^j · %^-1
+		cHH  = [ 0.393524e-12, 0.464083e-15, -0.621764e-18,  0.981126e-21, -0.121384e-23,  0.515111e-26] # cm^j · %^-2
+		cp   = [ 0.266809e-8,  0.695247e-15,  0.159070e-17, -0.303451e-20, -0.661489e-22,  0.178226e-24] # cm^j · Pa^-1
+		cpp  = [ 0.610508e-17, 0.227694e-22,  0.786323e-25, -0.174448e-27, -0.359791e-29,  0.978307e-32] # cm^j · Pa^-2
+		cTH  = [ 0.106776e-3, -0.168516e-6,   0.226201e-9,  -0.356457e-12,  0.437980e-15, -0.194545e-17] # cm^j · K · %^-1
+		cTp  = [ 0.77368e-6,   0.216404e-12,  0.581805e-15, -0.189618e-17, -0.198869e-19,  0.589381e-22] # cm^j · K · Pa^-1
+		cHp  = [-0.206365e-15, 0.300234e-19, -0.426519e-22,  0.684306e-25, -0.467320e-29,  0.126117e-30] # cm^j · %^-1 · Pa^-1
+		σref = 1e4/10.1    # cm^−1, reference wavenumber
+		Tref = 273.15+17.5 # K, reference temperature
+		pref = 75000       # Pa, reference pressure
+		Href = 10          #%, reference humdity in percent
+		########################################################################
+
+
+	if l > 15.05:
+		########################################################################
+		# model parameters for wavelengths between 16 and 24 microns
+		cref = [ 0.199436e-3,  0.299123e-8,  -0.214862e-10,  0.143338e-12,  0.122398e-14, -0.114628e-16] # cm^j
+		# something seems to be wrong with cT...
+		cT   = [ 0.621723e-1, -0.177074e-4,   0.152213e-6,  -0.954584-9,   -0.996706e-11,  0.921476e-13] # cm^j · K
+		cTT  = [-23.2409,      0.108557,     -0.102439e-2,   0.634072e-5,   0.762517e-7,  -0.675587e-9 ] # cm^j · K^2
+		cH   = [-0.772707e-7,  0.347237e-9,  -0.272675e-11,  0.170858e-13,  0.156889e-15, -0.150004e-17] # cm^j · %^-1
+		cHH  = [-0.326604e-12, 0.463606e-14, -0.364272e-16,  0.228756e-18,  0.209502e-20, -0.200547e-22] # cm^j · %^-2
+		cp   = [ 0.266827e-8,  0.120788e-14,  0.522646e-17,  0.783027e-19,  0.753235e-21, -0.228819e-24] # cm^j · Pa^-1
+		cpp  = [ 0.613675e-17, 0.585494e-22,  0.286055e-24,  0.425193e-26,  0.413455e-28, -0.812941e-32] # cm^j · Pa^-2
+		cTH  = [ 0.375974e-3, -0.171849e-5,   0.146704e-7,  -0.917231e-10, -0.955922e-12,  0.880502e-14] # cm^j · K · %^-1
+		cTp  = [ 0.778436e-6,  0.461840e-12,  0.306229e-14, -0.623183e-16, -0.161119e-18,  0.800756e-20] # cm^j · K · Pa^-1
+		cHp  = [-0.272614e-15, 0.304662e-18, -0.239590e-20,  0.149285e-22,  0.136086e-24, -0.130999e-26] # cm^j · %^-1 · Pa^-1
+		σref = 1e4/20      # cm^−1, reference wavenumber
+		Tref = 273.15+17.5 # K, reference temperature
+		pref = 75000       # Pa, reference pressure
+		Href = 10          #%, reference humdity in percent
+		########################################################################
+
+	s = 1e4/l # cm^-1
+	n = 1
+	H = H*100
+
+	for j in range(0, 6):
+		n += ( cref[j] + cT[j]*(1/T-1/Tref) + cTT[j]*(1/T-1/Tref)**2
+			+ cH[j]*(H-Href) + cHH[j]*(H-Href)**2
+			+ cp[j]*(p-pref) + cpp[j]*(p-pref)**2
+			+ cTH[j]*(1/T-1/Tref)*(H-Href)
+			+ cTp[j]*(1/T-1/Tref)*(p-pref)
+			+ cHp[j]*(H-Href)*(p-pref) ) * (s-sref)**j   
+	return n 
+
+
+def K2C(K):
+	'''
+		Kelvin to Celsius conversion
+	'''
+	return K - 273.15
+
+def hPA2mmHg(hPa):
+	'''
+		hPa to mmHg conversion
+	'''
+	return 0.7500616827 * hPa