Fixed the setup

.gitignore

@@ -1,3 +1,6 @@
 .DS_Store
 __pychache__/
 .ipynb_checkpoints/
+dist/
+AstroAtmosphere.egg-info/
+

AstroAtmosphere/__init__.py

@@ -5,6 +5,7 @@ from .ciddorModel import *
 from .dispersionModels import *
 from .refractionModels import *
 from .refractionModels import *
-from .examples import *
+from .misc import *
 from .refractivityModels import *
-
+from .neoslalib import *
+from .observatories import *

AstroAtmosphere/dispersionModels.py

@@ -14,6 +14,7 @@ from .ciddorModel import Observatory
 				The actual refraction is calculated by any of the following models:
 				- Plane parallel atmosphere
 				- Cassini spherical homogeneous atmosphere
+				- Mathar's Barometric exponential model
 				- Other methods described in Corbard et al., 2019.
 '''
 
@@ -38,39 +39,45 @@ class dispersion(Observatory):
 
 		Methods
 		-------
-		setReducedHeight(p, rho)
-			Sets the reduced atmospheric height.
-		planeParallel(n1, n2, zenith)
-			Returns the atmospheric dispersion assuming a plane parallel
-			atmosphere
 		cassini(n1, n2, zenith, reduced_height=None)
-			Returns the atmospheric dispersion asusming Cassini's homegeneous
+			Returns the atmospheric dispersion assuming Cassini's homegeneous
-			atmosphere model
+			atmosphere model.
-		tan5(n1, n2, zenith, reduced_height=None)
-			Returns the atmospheric dispersion following Oriani's theorem,
-			including a fifth order term in the expansion.
-		oriani(n1, n2, zenith, reduced_height=None)
-			Returns the atmospheric dispersion following the classic version
-			of Oriani's theorem.
-		corbard(n1, n2, zenith, reduced_height=None)
-			Returns the atmospheric dispersion following the method described 
-			in Corbard et al., 2019.
-		refractionIntegral(n1, n2, z0, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None)
-			Returns the atmospheric dispersion following the refraction integral.
-		psi(x)
-			Submethod for corbard()
-		set_rc()
-			Sets the local radius of curvature of the Earth at the observer.
-		set_H(p, rho)
-			Sets the reduced height for the conditions at the observer, assuming
-			an adiabatic atmosphere.
-		H_isotherm(T)
-			Sets the reduced height for the conditions at the observer, assuming
-			an isothermal atmosphere.
 		cassiniError(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)
 			Returns the uncertainty in the atmospheric dispersion, assuming 
 			Cassini's homegeneous atmosphere model and the given parameters.
+		corbard(n1, n2, zenith, reduced_height=None)
+			Returns the atmospheric dispersion following the method described 
+			in Corbard et al., 2019.
+		gammaCoefficients(n1, reduced_height=None)
+			Helper function for matharExponential(), that calculates the gamma
+			coefficients.
+		H_isotherm(T)
+			Sets the reduced height for the conditions at the observer, assuming
+			an isothermal atmosphere.
+		matharExponential(n1, n2, zenith, reduced_height=None)
+			Returns the atmospheric dispersion assuming	Mathar's barometric 
+			exponential model.	
+		oriani(n1, n2, zenith, reduced_height=None)
+			Returns the atmospheric dispersion following the classic version
+			of Oriani's theorem.
+		planeParallel(n1, n2, zenith)
+			Returns the atmospheric dispersion assuming a plane parallel
+			atmosphere.	
+		psi(x)
+			Submethod for corbard().
+		refractionIntegral(n1, n2, zenith, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None)
+			Returns the atmospheric dispersion following the refraction integral.	
+		setReducedHeight(p, rho)
+			Sets the reduced atmospheric height.
+		set_H(p, rho)
+			Sets the reduced height for the conditions at the observer, assuming
+			an adiabatic atmosphere.
+		set_rc()
+			Sets the local radius of curvature of the Earth at the observer.
+		tan5(n1, n2, zenith, reduced_height=None)
+			Returns the atmospheric dispersion following Oriani's theorem,
+			including a fifth order term in the expansion.		
 		Other Observatory() methods
 	'''
 	def __init__(self, lat, h):
@@ -79,44 +86,6 @@ class dispersion(Observatory):
 		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):
-		'''
-			Sets the reduced atmospheric height.
-
-			Parameters
-			----------
-			p 	: float
-				Atmospheric pressure in Pascal
-			rho : float
-				Density of the atmospheric air at the observer
-		'''
-		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.
-
-			Parameters
-			----------
-			n1 : float
-				Refractive index at the shortest wavelength of interest
-			n2 : float
-				Refractive index at the longest wavelength of interest
-			zenith : float
-				Angle of observation in degrees.
-				zenith = 0 means the observation is directly overhead.
-
-			Returns
-			-------
-			dR : float
-				Atmospheric dispersion in degrees
-		'''
-		_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):
 		'''
@@ -150,297 +119,7 @@ class dispersion(Observatory):
 		return np.degrees(_R1) - np.degrees(_R2)
 
 
-
+	def cassiniError(self, zenith, l1, l2, T, p, RH, xc, dl1=0, dl2=0, dT=0.2, dP=20,
-	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.
-
-			Parameters
-			----------
-			n1 : float
-				Refractive index at the shortest wavelength of interest
-			n2 : float
-				Refractive index at the longest wavelength of interest
-			zenith : float
-				Angle of observation in degrees.
-				zenith = 0 means the observation is directly overhead.
-			reduced_height : float (optional)
-				Reduced height of the atmosphere. Default value is calculated
-				from the object attributes.
-
-			Returns
-			-------
-			dR : float
-				Atmospheric dispersion in degrees
-		'''
-		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.
-
-			Parameters
-			----------
-			n1 : float
-				Refractive index at the shortest wavelength of interest
-			n2 : float
-				Refractive index at the longest wavelength of interest
-			zenith : float
-				Angle of observation in degrees.
-				zenith = 0 means the observation is directly overhead.
-			reduced_height : float (optional)
-				Reduced height of the atmosphere. Default value is calculated
-				from the object attributes.
-
-			Returns
-			-------
-			dR : float
-				Atmospheric dispersion in degrees
-		'''
-		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'.
-
-			Parameters
-			----------
-			n1 : float
-				Refractive index at the shortest wavelength of interest
-			n2 : float
-				Refractive index at the longest wavelength of interest
-			zenith : float
-				Angle of observation in degrees.
-				zenith = 0 means the observation is directly overhead.
-			reduced_height : float (optional)
-				Reduced height of the atmosphere. Default value is calculated
-				from the object attributes.
-
-			Returns
-			-------
-			dR : float
-				Atmospheric dispersion in degrees
-		'''
-		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.
-
-			Parameters
-			----------
-			n1 : float
-				Refractive index at the shortest wavelength of interest
-			n2 : float
-				Refractive index at the longest wavelength of interest
-			z0 : float
-				Angle of observation in degrees.
-				zenith = 0 means the observation is directly overhead.
-			R0 : float (optional)
-				Local radius of curvature of the Earth in km.
-			useStandardAtmosphere : bool
-				Flag to use the included USSA1976 data
-			heightData : 1D Array
-				If "useStandardAtmosphere" is False, the function will ask
-				for alternative atmospheric data to integrate over. This
-				gives the altitudinal points
-			rhoData : 1D Array
-				If useStandardAtmosphere is False, the function will ask
-				for alternative atmospheric data to integrate over. This
-				gives the atmospheric density at the points defined by 
-				"heightData".
-
-			Returns
-			-------
-			dR : float
-				Atmospheric dispersion in degrees
-		'''
-		if not R0:
-			R0 = self.rc
-		if useStandardAtmosphere:
-			_datafileLocation = os.path.join(os.path.dirname(__file__), 'data/1976USSA.txt')
-			_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt(_datafileLocation, 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):
-		'''
-			Submethod for corbard(). Described explicitly in Corbard et al.,2019.
-
-			Parameters
-			----------
-			x 	: 	float
-				Input of the function
-
-			Returns
-			-------
-			y 	: 	float
-				Output of the function
-		'''
-		_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.
-
-			Returns
-			-------
-			rc 	: float
-				Local radius of curvature of the Earth in kilometers.
-		'''
-		_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.
-
-			Parameters
-			----------
-			p 	: float
-				Atmospheric pressure in Pascal
-			rho : float
-				Density of the atmospheric air at the observer
-
-			Returns
-			-------
-			H 	: float
-				Reduced height of the atmosphere in meters
-		'''
-		_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.
-
-			Parameters
-			----------
-			T 	: float
-				Atmospheric temperature in Kelvin
-
-			Returns
-			-------
-			H 	: float
-				Reduced height of the atmosphere in meters
-		'''
-		_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, 
@@ -455,7 +134,7 @@ class dispersion(Observatory):
 
 			Parameters
 			----------
-			z 	: float
+			zenith : float
 				Angle of observation in degrees.
 				zenith = 0 means the observation is directly overhead.
 			l1 	: float
@@ -549,7 +228,7 @@ class dispersion(Observatory):
 		_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)
+		_alpha   = self.rc * np.sin(np.deg2rad(zenith)) / (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)
@@ -696,7 +375,7 @@ class dispersion(Observatory):
 
 
 		# Partial derivatives for derivation to z
-		_zenith = np.deg2rad(z)
+		_zenith = np.deg2rad(zenith)
 		_dz   	= np.deg2rad(dz)
 		_R 		= self.rc
 		_H 		= self.H
@@ -706,3 +385,524 @@ class dispersion(Observatory):
 		_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)
+
+
+	def independentCassiniError(self, n1, n2, zenith, dn1=0, dn2=0, dz=0, 
+									lat=None, h=None):
+		'''
+			Calculates the uncertainty in the atmospheric dispersion, 
+			following the Cassini model, without assumptions about the
+			refractivity model.
+			This involves a much shorter use the chain rule. 
+
+			Not taken into account is the errror propagation in the scale height
+			of the atmosphere, from my results this should give a small or even 
+			negligible difference.
+
+			Parameters
+			----------
+			n1 	: 	float
+				Refractive index at the first wavelength of interest at the 
+				observer
+			n2 	: 	float
+				Refractive index at the second wavelength of interest at the 
+				observer
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			dn1 : float
+				Uncertainty in the provided value of the refractive index
+				at the first wavelength
+			dn2 : float
+				Uncertainty in the provided value of the refractive index
+				at the second wavelength
+			dz 	: float (default: 0)
+				Uncertainty in the zenith angle, in degrees.
+			lat : float (default: None)
+				If entered, overrides the latitude attribute of the class.
+			h 	: float (default=None)
+				If entered, overrides the altitude attribute of the class.
+
+			Returns
+			-------
+			dRR	: float
+				Uncertainty in the atmospheric dispersion in degrees.
+		'''
+		if lat == None:
+			lat = self.lat
+		else:
+			self.lat = lat
+		if h == None:
+			h 	= self.h
+		else:
+			self.h = h
+
+		# Some precalculation and variable setting
+		_zenith = np.deg2rad(zenith)
+		_dz   	= np.deg2rad(dz)
+		_R 		= self.rc
+		_H 		= self.H
+
+		# Partial derivatives for derivation to n1, n2, and z
+		_dddn1  = - _R * np.sin(_zenith) / ( (_R + _H) * np.sqrt(1 - (n1 * _R * np.sin(_zenith) / (_R + _H))**2 ) )
+		_dddn2  =   _R * np.sin(_zenith) / ( (_R + _H) * np.sqrt(1 - (n2 * _R * np.sin(_zenith) / (_R + _H))**2 ) )
+		_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) )		
+		
+		# Combining the partial derivatives
+		_dispersionError = np.sqrt( (_dddn1 * dn1)**2 + (_dddn2 * dn2)**2 + (_dddz * dz)**2)
+
+		return np.degrees(_dispersionError)
+
+
+	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'.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the shortest wavelength of interest
+			n2 : float
+				Refractive index at the longest wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			dR : float
+				Atmospheric dispersion in degrees
+		'''
+		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 gammaCoefficients(self, n1, reduced_height=None):
+		'''
+			Submethod for matharExponential(). Calculates the gamma
+			coefficients.
+
+			Parameters
+			----------
+			n1 	: 	float
+				Refractive index at the wavelength of interest at the observer
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			g 	: 	array
+				List of values, corresponding to the gamma coefficients.
+				g[0] corresponds to g_1, g[1] to g_3, g[2] to g_5, and so on.
+		'''
+		if reduced_height:
+			H = reduced_height
+		else:
+			H = self.H
+
+		_chi0   = n1**2 - 1
+		_Khat 	= H / self.rc
+
+		_g1 = (1 - 5/24 * _Khat * _chi0**2 - 3/4 * _chi0 - 3/8 * _Khat**2 * _chi0 \
+			- _Khat - 6 * _Khat**3 + 5/8 * _chi0**2 + 63/128 * _chi0**4 + 3/8 \
+			* _Khat * _chi0 + 35/256 * _Khat * _chi0**3 + 2 * _Khat**2 + 5/36 \
+			* _Khat**2 * _chi0**2 - 35/64 * _chi0**3 + 9/16 * _Khat**3 * _chi0 \
+			- 120 * _Khat**5 + 24 * _Khat**4 - 231 / 512 * _chi0**5 - 63 / 640 \
+			* _Khat * _chi0**4 - 35 / 512 * _Khat**2 * _chi0**3 - 5 / 36 * _Khat**3\
+			* _chi0**2 - 9 / 8 * _Khat**4 * _chi0 )
+		
+		_g3 = (35/96 * _chi0**3 - 5/12 * _chi0**2 + 333 * _Khat**4 * _chi0 - 2400 \
+			* _Khat**5 + 3605/4608 * _Khat**2 * _chi0**3 + 63/256 * _Khat * _chi0**4\
+			 - 819/16 * _Khat**3 * _chi0 - 35/18 * _Khat**2 * _chi0**2 - 35/96 \
+			 * _Khat * _chi0**3 + 69/8 * _Khat**2 * _chi0 + 5/8 * _Khat * _chi0**2 \
+			 - 3/2 * _Khat * _chi0 + 1355/216 * _Khat**3 * _chi0**2 + 10 * _Khat**2 \
+			 + 336 * _Khat**4 - 2 * _Khat + _chi0/2 - 21/64 * _chi0**4 - 54 * _Khat**3 \
+			 + 77/256 * _chi0**5 ) * 1/2
+		
+		_g5 = (-96 * _Khat**3 + 8 * _Khat**2 - 7/24 * _chi0**3 + 65/4 * _Khat**2 \
+			* _chi0**2 + 35/48 * _Khat * _chi0**3 + 30 * _Khat**2 * _chi0 - 5/3 \
+			* _Khat * _chi0**2 - 3 * _Khat * _chi0 + 20091/8 * _Khat**4 * _chi0 \
+			- 32605/216 * _Khat**3 * _chi0**2 - 35/9 * _Khat**2 * _chi0**3 + 984 \
+			* _Khat**4 + _chi0**2/3 - 10200 * _Khat**5 - 77/320 * _chi0**5 + 21/80 \
+			* _chi0**4 - 1089/4 * _Khat**3 * _chi0 - 7/16 * _Khat * _chi0**4 ) * 3/8
+		
+		_g7 = (-17520 * _Khat**5 - 9/40 * _chi0**4 - 8755/12 * _Khat**3 * _chi0**2 \
+			+ 1056 * _Khat**4 + _chi0**3/4 + 23517/4 * _Khat**4 * _chi0 - 11/3 \
+			* _Khat * _chi0**2 + 21 * _Khat**2 * _chi0 - 7/4 * _Khat * _chi0**3 \
+			+ 505/9 * _Khat**2 * _chi0**2 - 801/2 * _Khat**3 * _chi0 + 33/160 \
+			* _chi0**5 - 48 * _Khat**3 + 63/80 * _Khat * _chi0**4 + 2345/96 \
+			* _Khat**2 * _chi0**3) * 5/16
+		
+		_g9 = (_chi0**4 / 5 + 384 * _Khat**4 - 13440 * _Khat**5 + 5562 * _Khat**4 \
+			* _chi0 + 3125/36 * _Khat**2 * _chi0**3 - 11/60 * _chi0**5 - 25/6 \
+			* _Khat * _chi0**3 - 26950/27 * _Khat**3 * _chi0**2 + 340/9 * _Khat**2 \
+			* _chi0**2 - 180 * _Khat**3 * _chi0 - 9/5 * _Khat * _chi0**4) * 35/128
+
+		g = np.array([_g1, _g3, _g5, _g7, _g9])
+		g*= n1*_chi0/2
+
+		return g
+
+
+	def H_isotherm(self, T):
+		'''
+			Calculates the reduced height of the atmosphere, 
+			assuming an isotherm distribution.
+
+			Parameters
+			----------
+			T 	: float
+				Atmospheric temperature in Kelvin
+
+			Returns
+			-------
+			H 	: float
+				Reduced height of the atmosphere in meters
+		'''
+		_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 matharExponential(self, n1, n2, zenith, reduced_height=None):
+		'''
+			Model on atmospheric refraction described by Richard Mathar in 
+			"A Barometric Exponential Model of the Atmosphere’s Refractive Index:  ZenithAngles and Second Order Aberration in the Entrance Pupil" (2020)
+			ArXiv id: 2004.11808 (https://arxiv.org/abs/2004.11808)
+			
+			The code in this function is adapted from eqs 3, 16, 17 and 18 in this
+			work. It is an adaptation of the RefrExpo.cxx C++ file written by Mathar.
+
+			This model assumes the refractive index, n, is related to the susceptibility, 
+			chi, following: n(r) = sqrt(1+chi(r)) (eq. 9).
+
+			R. Mathar then assumes that the susceptibility is proportional to the
+			air density and reduced height, following eq. 10:
+			chi(r) = chi0 * exp(-(r-rho)/K).
+
+			The refraction is described using a Taylor expansion; R in terms of
+			the odd powers of the tangent of the observed zenith angle (eq. 3), not
+			unlike Oriani's theorem. This requires coefficients g_i.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the wavelength of interest at the observer
+			n2 : float
+				Refractive index at the longest wavelength of interest at the
+				observer
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			dR : float
+				Atmospheric refraction in degrees
+		'''
+		if reduced_height:
+			_H = reduced_height
+		else:
+			_H = self.H
+
+		_zenith	= np.deg2rad(zenith)
+		_rc 	= self.rc # rho
+
+		# Atmospheric refraction for the first wavelength
+		_gamma1 = self.gammaCoefficients(n1, reduced_height=_H)
+		_R1 = sum([_gamma1[i]*np.tan(_zenith)**(2*i+1) for i in range(len(_gamma1))])
+
+		# Atmospheric refraction for the second wavelength
+		_gamma2 = self.gammaCoefficients(n2, reduced_height=_H)
+		_R2 = sum([_gamma2[i]*np.tan(_zenith)**(2*i+1) for i in range(len(_gamma2))])
+
+		# The atmospheric dispersion in radians
+		_dR = _R1 - _R2
+
+		return np.degrees(_dR)
+
+
+	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.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the shortest wavelength of interest
+			n2 : float
+				Refractive index at the longest wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			dR : float
+				Atmospheric dispersion in degrees
+		'''
+		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 planeParallel(self, n1, n2, zenith):
+		'''
+			Refraction based on a flat atmosphere, easily derived from Snell's law.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the shortest wavelength of interest
+			n2 : float
+				Refractive index at the longest wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+
+			Returns
+			-------
+			dR : float
+				Atmospheric dispersion in degrees
+		'''
+		_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 psi(self, x):
+		'''
+			Submethod for corbard(). Described explicitly in Corbard et al.,2019.
+
+			Parameters
+			----------
+			x 	: 	float
+				Input of the function
+
+			Returns
+			-------
+			y 	: 	float
+				Output of the function
+		'''
+		_f   = lambda a: np.exp(-a**2)
+		_int = integrate.quad(_f, x, np.inf)
+		return np.exp(x**2) * _int[0]
+
+
+	def refractionIntegral(self, n1, n2, zenith, 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.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the shortest wavelength of interest
+			n2 : float
+				Refractive index at the longest wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			R0 : float (optional)
+				Local radius of curvature of the Earth in km.
+			useStandardAtmosphere : bool
+				Flag to use the included USSA1976 data
+			heightData : 1D Array
+				If "useStandardAtmosphere" is False, the function will ask
+				for alternative atmospheric data to integrate over. This
+				gives the altitudinal points
+			rhoData : 1D Array
+				If useStandardAtmosphere is False, the function will ask
+				for alternative atmospheric data to integrate over. This
+				gives the atmospheric density at the points defined by 
+				"heightData".
+
+			Returns
+			-------
+			dR : float
+				Atmospheric dispersion in degrees
+		'''
+		if not R0:
+			R0 = self.rc
+		if useStandardAtmosphere:
+			_datafileLocation = os.path.join(os.path.dirname(__file__), 'data/1976USSA.txt')
+			_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt(_datafileLocation, unpack=True, skip_header=6)
+		else:
+			_Dens = rhoData
+			_h 	  = heightData 
+
+		
+		_rList 	= []
+		zenith 	= np.deg2rad(zenith)
+		
+		# 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(zenith) )
+			
+			# 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 setReducedHeight(self, p, rho):
+		'''
+			Sets the reduced atmospheric height.
+
+			Parameters
+			----------
+			p 	: float
+				Atmospheric pressure in Pascal
+			rho : float
+				Density of the atmospheric air at the observer
+		'''
+		self.H  	= self.set_H(p, rho) 		# Reduced height of the atmosphere assuming ideal gas law
+
+
+	def set_H(self, p, rho):
+		'''
+			Calculates the reduced height of the atmosphere, assuming ideal gas law.
+
+			Parameters
+			----------
+			p 	: float
+				Atmospheric pressure in Pascal
+			rho : float
+				Density of the atmospheric air at the observer
+
+			Returns
+			-------
+			H 	: float
+				Reduced height of the atmosphere in meters
+		'''
+		_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 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.
+
+			Returns
+			-------
+			rc 	: float
+				Local radius of curvature of the Earth in kilometers.
+		'''
+		_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 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.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the shortest wavelength of interest
+			n2 : float
+				Refractive index at the longest wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			dR : float
+				Atmospheric dispersion in degrees
+		'''
+		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]

AstroAtmosphere/examples.py

@@ -1,205 +0,0 @@
-#!/usr/bin/env Python3
-from .ciddorModel import Observatory
-from .refractionModels import refraction
-from .dispersionModels import dispersion
-
-'''
-@ Author: 		Joost van den Born
-@ Contact: 		born@astron.nl
-@ Description: 	Some functions for quick testing.
-'''
-
-
-def quick_refractive_index(l, conditions='STANDARD', T=288.15, p=101325, 
-							H=0.0, xc=450, lat=0, h=0):
-	'''
-		Function to quickly calculate the refractive index
-		of atmospheric air at reference conditions, given 
-		a wavelength l.
-
-		Parameters
-		----------
-		l 	: float
-			Wavelength in microns
-		conditions: string
-			'STANDARD' or 'CERRO_ARMAZONES', refers
-			to the standard atmospheric conditions.
-		T 	: float (optional)
-			Temperature in Kelvin
-		p 	: float (optional)
-			Atmospheric pressure in Pa
-		H 	: float (optional)
-			Relative humidity
-		xc 	: float (optional)
-			CO2 density in parts per million
-		lat : float (optional)
-			Latitude of the observer in degrees
-		h 	: float (optional)
-			Altitude of the observer in meters
-	
-		Returns
-		-------
-		n 	: 	float
-			Refractive index for the given conditions.
-	'''
-	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', T=288.15, p=101325, 
-							H=0.0, xc=450, lat=0, h=0):
-	'''
-		Function to quickly calculate the refraction
-		of atmospheric air at reference conditions, 
-		given a wavelength l.
-
-		Parameters
-		----------
-		l  	: 	float 
-			Wavelength in microns
-		z  	: 	float 
-			Zenith angle of the observation
-		conditions 	: string
-			'STANDARD' or 'CERRO_ARMAZONES', refers
-			to the standard atmospheric conditions.
-		T  	: 	float (optional)
-			Temperature in Kelvin
-		p  	: 	float (optional)
-			Atmospheric pressure in Pa
-		H  	: 	float (optional)
-			Relative humidity
-		xc  : 	float (optional)
-			CO2 density in parts per million
-		lat : 	float (optional)
-			Latitude of the observer in degrees
-		h  	: 	float (optional)
-			Altitude of the observer in meters
-	
-		Returns
-		-------
-		R 	: 	float
-			Refraction in degrees
-	'''
-	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', T=288.15, 
-						p=101325, H=0.0, xc=450, lat=0, h=0):
-	'''
-		Function to quickly calculate the dispersion
-		of atmospheric air at reference conditions, 
-		given a wavelength l.
-
-		Parameters
-		----------
-		l1 	: 	float
-			Shortest wavelength of interest in micron.
-		l2 	: 	float
-			Longest wavelength of interest in micron.
-		z 	: 	float
-			Zenith angle of the observation
-		conditions 	: string
-			'STANDARD' or 'CERRO_ARMAZONES', refers
-			to the standard atmospheric conditions.
-		T 	: 	float (optional)
-			Temperature in Kelvin
-		p 	: 	float (optional)
-			Atmospheric pressure in Pa
-		H 	: 	float (optional)
-			Relative humidity
-		xc 	: 	float (optional)
-			CO2 density in parts per million
-		lat : 	float (optional)
-			Latitude of the observer in degrees
-		h 	: 	float (optional)
-			Altitude of the observer in meters
-	
-		Returns
-		-------
-		dR 	: float
-		 	Atmospheric dispersion in degrees
-	'''
-	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)

AstroAtmosphere/misc.py

@@ -0,0 +1,421 @@
+#!/usr/bin/env Python3
+from .ciddorModel import Observatory
+from .refractionModels import refraction
+from .dispersionModels import dispersion
+from .neoslalib import sla_refro
+from .observatories import observatories
+from .refractivityModels import Edlen1953, Mathar
+import numpy as np
+import warnings
+
+'''
+@ Author: 		Joost van den Born
+@ Contact: 		born@astron.nl
+@ Description: 	Some functions for quick testing.
+'''
+def get_conditions(location):
+	'''
+		Returns the mean conditions at various site. Currently
+		supported locations:
+		- 'STANDARD'
+		- 'CERRO_ARMAZONES'
+		- 'CERRO_PARANAL'
+		- 'LA_PALMA'
+		- 'LA_SILLA'
+		- 'LAS_CAMPANAS'
+		- 'MAUNA_KEA'
+
+		Parameters
+		----------
+		location 	: string
+			String denoting the desired location
+
+		Returns
+		-------
+		observatory : dict
+			Dictionary containing the atmospheric parameters
+	'''
+	observatory = observatories[location]
+	return observatory
+
+
+def Filippenko1982(l1, l2, zenith, conditions='STANDARD', T=288.15, p=101325,
+							RH=0.0, xc=450, lat=0, h=0, f=None):
+	'''
+		Calculates the refraction following the 1982 paper by 
+		Alexei V. Filippenko.
+
+		The refraction is based on a flat atmosphere and the refractivity
+		values are calculated from a modified Edlen (1953) equation.
+
+		Parameters
+		----------
+		l1 	: 	float
+			Shortest wavelength of interest in micron.
+		l2 	: 	float
+			Longest wavelength of interest in micron.
+		zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+		conditions 	: string
+			Refers to the standard conditions of the 
+			observatory. See get_conditions() for more info.
+			If set to 'None', the conditions specified 
+			after will be used.
+		T 	: float (optional)
+			Temperature in Kelvin
+		p 	: float (optional)
+			Atmospheric pressure in Pa
+		RH 	: float (optional)
+			Relative humidity
+		xc 	: float (optional)
+			CO2 density in parts per million
+		lat : float (optional)
+			Latitude of the observer in degrees
+		h 	: float (optional)
+			Altitude of the observer in meters
+		f 	: float
+			Water vapour pressure (default=None). If None, it is
+			calculated automatically.
+
+		Returns
+		-------
+		R : float
+			Atmospheric refraction in degrees
+	'''
+
+	if conditions is not None:
+		T, p, RH, xc, lat, h, _ = get_conditions(conditions).values()
+
+	n1 = Edlen1953(l1, T=T, p=p, RH=RH, f=f)
+	n2 = Edlen1953(l2, T=T, p=p, RH=RH, f=f)
+	dR  = (n1 - n2) * np.tan(np.deg2rad(zenith))
+	return np.degrees(dR)
+
+
+def quick_refractive_index(l, conditions='STANDARD', T=288.15, p=101325, 
+							RH=0.0, xc=450, lat=0, h=0):
+	'''
+		Function to quickly calculate the refractive index
+		of atmospheric air at reference conditions, given 
+		a wavelength l.
+
+		Parameters
+		----------
+		l 	: float
+			Wavelength in microns. If the wavelength is above
+			2.5 microns, then the Mathar refractivity model
+			is used.
+		conditions: string
+			Refers to the standard conditions of the 
+			observatory. See get_conditions() for more info.
+			If set to 'None', the conditions specified 
+			after will be used.
+		T 	: float (optional)
+			Temperature in Kelvin
+		p 	: float (optional)
+			Atmospheric pressure in Pa
+		RH 	: float (optional)
+			Relative humidity
+		xc 	: float (optional)
+			CO2 density in parts per million
+		lat : float (optional)
+			Latitude of the observer in degrees
+		h 	: float (optional)
+			Altitude of the observer in meters
+	
+		Returns
+		-------
+		n 	: 	float
+			Refractive index for the given conditions.
+	'''
+	if conditions is not None:
+		T, p, RH, xc, lat, h, _ = get_conditions(conditions).values()
+
+	if l <= 2.5:
+		# Initializing dispersion model
+		at  = Observatory()
+
+		# Calculating indices of refraction for l
+		n 	= at.n_tph(l=l, T=T, p=p, RH=RH, xc=xc)
+	else:
+		# Using the Mathar description for infrared wavelengths
+		n   = Mathar(l=l, T=T, p=p, RH=RH)
+	return n
+
+
+def quick_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325, 
+							RH=0.0, xc=450, lat=0, h=0):
+	'''
+		Function to quickly calculate the refraction
+		of atmospheric air at reference conditions, 
+		given a wavelength l.
+
+		Parameters
+		----------
+		l  	: 	float 
+			Wavelength in microns. If the wavelength is above 2.5 micron
+			then the Mathar refractivity model is used.
+		zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+		conditions 	: string
+			Refers to the standard conditions of the 
+			observatory. See get_conditions() for more info.
+			If set to 'None', the conditions specified 
+			after will be used.
+		T  	: 	float (optional)
+			Temperature in Kelvin
+		p  	: 	float (optional)
+			Atmospheric pressure in Pa
+		RH  	: 	float (optional)
+			Relative humidity
+		xc  : 	float (optional)
+			CO2 density in parts per million
+		lat : 	float (optional)
+			Latitude of the observer in degrees
+		h  	: 	float (optional)
+			Altitude of the observer in meters
+	
+		Returns
+		-------
+		R 	: 	float
+			Refraction in degrees
+	'''
+	if conditions is not None:
+		T, p, RH, xc, lat, h, _ = get_conditions(conditions).values()
+
+	# Initializing dispersion model
+	at  = Observatory()
+
+	if l <= 2.5:
+		# Calculating indices of refraction for l
+		n 	= at.n_tph(l=l, T=T, p=p, RH=RH, xc=xc)
+	else:
+		# Using the Mathar description for infrared wavelengths
+		n   = Mathar(l=l, T=T, p=p, RH=RH)
+
+	# Density of the atmosphere (following CIPM-81/91 equations)
+	rho = at.rho(p=p, T=T, RH=RH, xc=xc)
+
+	# Initializing refraction model and setting the reduced height
+	ref = refraction(lat, h)
+	ref.setReducedHeight(p, rho)
+	return ref.cassini(n, zenith)
+
+
+def quick_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15, 
+						p=101325, RH=0.0, xc=450, lat=0, h=0):
+	'''
+		Function to quickly calculate the dispersion
+		of atmospheric air at reference conditions, 
+		given a wavelength l.
+
+		Parameters
+		----------
+		l1 	: 	float
+			Shortest wavelength of interest in micron. If the wavelength
+			is above 2.5 microns, then the Mathar model is used.
+		l2 	: 	float
+			Longest wavelength of interest in micron. If the wavelength
+			is above 2.5 microns, then the Mathar model is used.
+		zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+		conditions 	: string
+			Refers to the standard conditions of the 
+			observatory. See get_conditions() for more info.
+			If set to 'None', the conditions specified 
+			after will be used.
+		T 	: 	float (optional)
+			Temperature in Kelvin
+		p 	: 	float (optional)
+			Atmospheric pressure in Pa
+		RH 	: 	float (optional)
+			Relative humidity
+		xc 	: 	float (optional)
+			CO2 density in parts per million
+		lat : 	float (optional)
+			Latitude of the observer in degrees
+		h 	: 	float (optional)
+			Altitude of the observer in meters
+	
+		Returns
+		-------
+		dR 	: float
+		 	Atmospheric dispersion in degrees
+	'''
+	if conditions is not None:
+		T, p, RH, xc, lat, h, _ = get_conditions(conditions).values()
+
+	# Initializing dispersion model
+	at  = Observatory()
+
+	if l1 <= 2.5 and l2 <= 2.5:
+		# Calculating indices of refraction for l1 and l2
+		n1 	= at.n_tph(l=l1, T=T, p=p, RH=RH, xc=xc)
+		n2 	= at.n_tph(l=l2, T=T, p=p, RH=RH, xc=xc)
+	elif l1 > 1.4 and l2 > 1.4:
+		# Using the Mathar description for infrared wavelengths
+		n1   = Mathar(l=l1, T=T, p=p, RH=RH)
+		n2   = Mathar(l=l2, T=T, p=p, RH=RH)
+	else:
+		warnings.warn('Using two refractivity models (Ciddor and Mathar). \
+Make sure you trust the results!')
+		l1n  = np.min([l1, l2])
+		l2n  = np.max([l1, l2])
+		n1 	 = at.n_tph(l=l1n, T=T, p=p, RH=RH, xc=xc)
+		n2   = Mathar(l=l2n, T=T, p=p, RH=RH)
+
+	# Density of the atmosphere (following CIPM-81/91 equations)
+	rho = at.rho(p=p, T=T, RH=RH, xc=xc)
+
+	# Initializing refraction model and setting the reduced height
+	disp = dispersion(lat, h)
+	disp.setReducedHeight(p, rho)
+	return disp.cassini(n1, n2, zenith)
+
+
+def slalib_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325, 
+				RH=0.0, xc=450, lat=0, h=0, tlr=0.0065, eps=1e-8, use_pyslalib=False):
+	'''
+		Calculates the atmospheric refraction assuming the SLALIB
+		package. 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.
+		
+		The pySLALIB FORTRAN has been ported, so pyslalib is not
+		needed, but can still be used. There are minor numerical
+		differences, likely due to precision or declaration issues.
+
+		Parameters
+		---------- 
+		l 	:	float
+			Wavelength in microns
+		zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+		conditions 	: string
+			Refers to the standard conditions of the 
+			observatory. See get_conditions() for more info.
+			If set to 'None', the conditions specified 
+			after will be used.
+		T 	: 	float (optional)
+			Temperature in Kelvin
+		p 	: 	float (optional)
+			Atmospheric pressure in Pa
+		RH 	: 	float (optional)
+			Relative humidity
+		xc 	: 	float (optional)
+			CO2 density in parts per million
+		lat : 	float (optional)
+			Latitude of the observer in degrees (default=0)
+		h 	: 	float (optional)
+			Altitude of the observer in meters
+		tlr :	float
+			Temperature lapse rate (default=0.0065)
+		eps :	float
+			Precision at which the calculation is aborted (default=1e-8)
+		use_pyslalib : bool
+			Flag that indicates that the user wants to use the pySlalib
+			package. If False, the internal python port of the original
+			FORTRAN code will be used. (default=False).
+			
+		Returns
+		-------
+		R 	: 	float
+			Refraction in degrees
+	'''
+
+	if conditions is not None:
+		T, p, RH, xc, lat, h, _ = get_conditions(conditions).values()
+
+	# Converting units to the units that sla_refro() requires.
+	_zenith = np.deg2rad(zenith) 	# From degrees to radians
+	_P  	= p / 100 				# From Pa to mbar
+	_phi  	= np.deg2rad(lat) 		# From degrees to radians
+
+	R = []
+	if isinstance(l, float) or isinstance(l, int): 
+		l = [l]
+	if use_pyslalib:
+		try:
+		    from pyslalib import slalib
+		except ImportError as e:
+		    print('pySLALIB package not found! Try installing it first.')
+		for wl in l:
+			# If using original slalib (e.g. PySLALIB)
+			Refrac = slalib.sla_refro(_zenith, h, T, _P, RH, wl, _phi, tlr, eps)
+			R.append(np.degrees(Refrac))
+	
+	# If using internal slalib port
+	else:
+		for wl in l:
+			Refrac = sla_refro(_zenith, h, T, _P, RH, wl, _phi, tlr, eps)
+			R.append(np.degrees(Refrac))
+
+	if len(l) == 1:
+		return R[0]
+	else:
+		return np.asarray(R)
+
+def slalib_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15, p=101325, 
+				RH=0.0, xc=450, lat=0, h=0, tlr=0.0065, eps=1e-8, use_pyslalib=False):
+	'''
+		Calculates the atmospheric refraction assuming the SLALIB
+		package. 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.
+		
+		The pySLALIB FORTRAN has been ported, so pyslalib is not
+		needed, but can still be used. There are minor numerical
+		differences, likely due to precision or declaration issues.
+
+		Parameters
+		---------- 
+		l1 	: 	float
+			Shortest wavelength of interest in micron.
+		l2 	: 	float
+			Longest wavelength of interest in micron.
+		zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+		conditions 	: string
+			Refers to the standard conditions of the 
+			observatory. See get_conditions() for more info.
+			If set to 'None', the conditions specified 
+			after will be used.
+		T 	: 	float (optional)
+			Temperature in Kelvin
+		p 	: 	float (optional)
+			Atmospheric pressure in Pa
+		RH 	: 	float (optional)
+			Relative humidity
+		xc 	: 	float (optional)
+			CO2 density in parts per million
+		lat : 	float (optional)
+			Latitude of the observer in degrees (default=0)
+		h 	: 	float (optional)
+			Altitude of the observer in meters
+		tlr :	float
+			Temperature lapse rate (default=0.0065)
+		eps :	float
+			Precision at which the calculation is aborted (default=1e-8)
+		use_pyslalib : bool
+			Flag that indicates that the user wants to use the pySlalib
+			package. If False, the internal python port of the original
+			FORTRAN code will be used. (default=False).
+			
+		Returns
+		-------
+		dR 	: 	float
+			Refraction in degrees
+		'''
+	R1 = slalib_refraction(l1, zenith, conditions, T, p, RH, xc,\
+	 lat, h, tlr, eps, use_pyslalib)
+
+	R2 = slalib_refraction(l2, zenith, conditions, T, p, RH, xc,\
+	 lat, h, tlr, eps, use_pyslalib)
+
+	dR = R1 - R2
+	return dR

AstroAtmosphere/neoslalib.py

@@ -0,0 +1,754 @@
+#!/usr/bin/env Python3
+import numpy as np
+
+def sla_atms(RT, TT, DNT, GAMAL, R):
+	'''
+		Internal routine used by sla_refro(), ported from the SLALIB 
+		FORTRAN code (Copyright P.T.Wallace.  All rights reserved.)
+
+		The method is also described in "Explanatory Supplement 
+		Astronomical Almanac" by P.K. Seidelmann, 1992, section 3.281.
+
+		Refractive index and derivative with respect to height for the
+		stratosphere.
+
+		Parameters
+		----------
+		RT 		: 	float
+			Height of tropopause from center of the Earth (meters)
+		TT 		: 	float
+			Temperature at the tropopause (K)
+		DNT     : 	float
+			Refractive index at the tropopause
+		GAMAL 	: 	float
+			Constant of the atmospheric model = G*MD/R
+		R 		: 	float
+			Current distance from the centre of the Earth (metre)
+
+		Returns
+		-------
+		DN 		: 	float
+			Refractive index at R
+		RDNDR 	: 	float
+			R * rate the refractive index is changing at R
+	'''
+	_B 		= GAMAL / TT
+	_W 		= (DNT - 1) * np.exp( -_B * (R - RT) )
+	DN  	= 1+_W
+	RDNDR 	= -R * _B * _W
+	return DN, RDNDR
+
+
+def sla_atmt(R0, T0, ALPHA, GAMM2, DELM2, C1, C2, C3, C4, C5, C6, R):
+	'''
+		Internal routine used by sla_refro(), ported from the SLALIB 
+		FORTRAN code (Copyright P.T.Wallace.  All rights reserved.)
+
+		The method is also described in "Explanatory Supplement 
+		Astronomical Almanac" by P.K. Seidelmann, 1992, section 3.281.
+
+		Refractive index and derivative with respect to height for the
+		troposphere.
+
+		Parameters
+		----------
+		R0 		: 	float
+			Height of observer from center of the Earth (meter)
+		T0 		: 	float
+			Temperature at the observer (K)
+		ALPHA 	: 	float
+			Alpha
+		GAMM2 	: 	float 
+			Gamma minus 2 - see HMNAO paper
+		DELM2 	: 	float
+			Delta minus 2
+		C1  	: 	float
+			Useful term
+		C2  	: 	float
+			Useful term
+		C3  	: 	float
+			Useful term
+		C4  	: 	float
+			Useful term
+		C5  	: 	float
+			Useful term
+		C6  	: 	float
+			Useful term
+		R    	:  	float
+			Current distance from the center of the Earth (meters)
+
+		Returns
+		-------
+		T 		: 	float
+			Temperature at R (K)
+		DN 		:   float 
+			Refractive index at R
+		RDNDR 	: 	float
+			R * rate the refractive index is changing at R
+
+		Note that in the optical case C5 and C6 are zero.
+	'''
+	T 		= max(min(T0 - ALPHA * (R - R0), 320), 100)
+	TT0 	= T / T0
+	TT0GM2 	= TT0**GAMM2
+	TT0DM2 	= TT0**DELM2
+	DN 		= 1 + ( C1 * TT0GM2 - (C2 - C5 / T) * TT0DM2) * TT0
+	RDNDR 	= R * (-C3 * TT0GM2 + (C4 - C6 / TT0) * TT0DM2)
+	return T, DN, RDNDR
+
+
+def sla_drange(a):
+	'''
+		Helper function for sla_refro().
+	'''
+	drange = a % (2*np.pi)
+	if drange >= np.pi:
+		drange -= np.abs(2*np.pi)*np.sign(a)
+	return drange
+
+
+def sla_refi(dn, rdndr):
+	'''
+		Helper function for sla_refro(), that calculates
+		the refraction integrand.
+	'''
+	return rdndr/(dn+rdndr)
+
+def sla_refro(ZOBS, HM, TDK, PMB, RH, WL, PHI, TLR, EPS):
+	'''
+		This is a direct port from the original FORTRAN SLA_REFRO()
+		routine.
+
+		Calculates Atmospheric refraction for radio and 
+		optical/IR wavelengths, ported from the original SLALIB
+		FORTRAN code (Copyright P.T.Wallace.  All rights reserved.).
+
+		The method is also described in "Explanatory Supplement 
+		Astronomical Almanac" by P.K. Seidelmann, 1992, section 3.281.
+		
+		Parameters
+		----------
+		ZOBS 	:	float
+			Observed zenith distance of the source (radian)
+		HM 		:	float
+			Height of the observer above sea level (meter)
+		TDK 	:	float
+			Ambient temperature at the observer (K)
+		PMB 	:	float
+			Pressure at the observer (millibar)
+		RH 		:	float
+			Relative humidity at the observer (range 0-1)
+		WL 		:	float
+			Effective wavelength of the source (micrometre)
+		PHI 	:	float
+			Latitude of the observer (radian, astronomical)
+		TLR 	:	float
+			Temperature lapse rate in the troposphere (K/metre)
+		EPS 	:	float
+			Precision required to terminate iteration (radian)
+		
+		Returns
+		-------
+		REF 	:	float
+			Refraction: in vacuo ZD minus observed ZD (radian)
+		
+		Notes
+		-----
+		1  A suggested value for the TLR argument is 0.0065D0.  The
+		   refraction is significantly affected by TLR, and if studies
+		   of the local atmosphere have been carried out a better TLR
+		   value may be available.  The sign of the supplied TLR value
+		   is ignored.
+		2  A suggested value for the EPS argument is 1D-8.  The result is
+		   usually at least two orders of magnitude more computationally
+		   precise than the supplied EPS value.
+		3  The routine computes the refraction for zenith distances up
+		   to and a little beyond 90 deg using the method of Hohenkerk
+		   and Sinclair (NAO Technical Notes 59 and 63, subsequently adopted
+		   in the Explanatory Supplement, 1992 edition - see section 3.281).
+		4  The code is a development of the optical/IR refraction subroutine
+		   AREF of C.Hohenkerk (HMNAO, September 1984), with extensions to
+		   support the radio case.  Apart from merely cosmetic changes, the
+		   following modifications to the original HMNAO optical/IR refraction
+		   code have been made:
+		   .  The angle arguments have been changed to radians.
+		   .  Any value of ZOBS is allowed (see note 6, below).
+		   .  Other argument values have been limited to safe values.
+		   .  Murray's values for the gas constants have been used
+		      (Vectorial Astrometry, Adam Hilger, 1983).
+		   .  The numerical integration phase has been rearranged for
+		      extra clarity.
+		   .  A better model for Ps(T) has been adopted (taken from
+		      Gill, Atmosphere-Ocean Dynamics, Academic Press, 1982).
+		   .  More accurate expressions for Pwo have been adopted
+		      (again from Gill 1982).
+		   .  The formula for the water vapour pressure, given the
+		      saturation pressure and the relative humidity, is from
+		      Crane (1976), expression 2.5.5.
+		   .  Provision for radio wavelengths has been added using
+		      expressions devised by A.T.Sinclair, RGO (private
+		      communication 1989).  The refractivity model currently
+		      used is from J.M.Rueger, "Refractive Index Formulae for
+		      Electronic Distance Measurement with Radio and Millimetre
+		      Waves", in Unisurv Report S-68 (2002), School of Surveying
+		      and Spatial Information Systems, University of New South
+		      Wales, Sydney, Australia.
+		   .  The optical refractivity for dry air is from Resolution 3 of
+		      the International Association of Geodesy adopted at the XXIIth
+		      General Assembly in Birmingham, UK, 1999.
+		   .  Various small changes have been made to gain speed.
+		5  The radio refraction is chosen by specifying WL > 100 micrometres.
+		   Because the algorithm takes no account of the ionosphere, the
+		   accuracy deteriorates at low frequencies, below about 30 MHz.
+		6  Before use, the value of ZOBS is expressed in the range +/- pi.
+		   If this ranged ZOBS is -ve, the result REF is computed from its
+		   absolute value before being made -ve to match.  In addition, if
+		   it has an absolute value greater than 93 deg, a fixed REF value
+		   equal to the result for ZOBS = 93 deg is returned, appropriately
+		   signed.
+		7  As in the original Hohenkerk and Sinclair algorithm, fixed values
+		   of the water vapour polytrope exponent, the height of the
+		   tropopause, and the height at which refraction is negligible are
+		   used.
+		8  The radio refraction has been tested against work done by
+		   Iain Coulson, JACH, (private communication 1995) for the
+		   James Clerk Maxwell Telescope, Mauna Kea.  For typical conditions,
+		   agreement at the 0.1 arcsec level is achieved for moderate ZD,
+		   worsening to perhaps 0.5-1.0 arcsec at ZD 80 deg.  At hot and
+		   humid sea-level sites the accuracy will not be as good.
+		9  It should be noted that the relative humidity RH is formally
+		   defined in terms of "mixing ratio" rather than pressures or
+		   densities as is often stated.  It is the mass of water per unit
+		   mass of dry air divided by that for saturated air at the same
+		   temperature and pressure (see Gill 1982).
+		10 The algorithm is designed for observers in the troposphere.  The
+		   supplied temperature, pressure and lapse rate are assumed to be
+		   for a point in the troposphere and are used to define a model
+		   atmosphere with the tropopause at 11km altitude and a constant
+		   temperature above that.  However, in practice, the refraction
+		   values returned for stratospheric observers, at altitudes up to
+		   25km, are quite usable.
+	'''
+
+	# Fixed parameters
+	D93 	= np.deg2rad(93)# 93 degrees in radians
+	GCR 	= 8314.32 		# Universal gas constant
+	DMD 	= 28.9644 		# Molecular weight of dry air
+	DMW 	= 18.0152 		# Molecular weight of water vapour
+	S 		= 6378120 		# Mean Earth radius (metre)
+	DELTA 	= 18.36 		# Exponent of temperature dependence of water vapour pressure
+	HT 		= 11000			# Height of tropopause (metre)
+	HS 		= 80000			# Upper limit for refractive effects (metre)
+	ISMAX 	= 16384			# Numerical integration: maximum number of strips.
+	
+	# Transform ZOBS into the normal range.
+	ZOBS1 = sla_drange(ZOBS)
+	ZOBS2 = min(abs(ZOBS1),D93)
+
+	# Keep other arguments within safe bounds.
+	HMOK  = min(max(HM,-1E3),HS)
+	TDKOK = min(max(TDK,100),500)
+	PMBOK = min(max(PMB,0),10000)
+	RHOK  = min(max(RH,0),1)
+	WLOK  = max(WL,0.1)
+	ALPHA = min(max(abs(TLR),0.001),0.01)
+
+	# Tolerance for iteration.
+	TOL   = min(max(abs(EPS),1E-12),0.1)/2
+
+	# Decide whether optical/IR or radio case - switch at 100 microns.
+	OPTIC = (WLOK<=100)
+
+	# Set up model atmosphere parameters defined at the observer.
+	WLSQ  = WLOK*WLOK
+
+	# Local gravity, can be made consistent with AstroAtmosphere
+	GB    = 9.784*(1-0.0026*np.cos(PHI+PHI)-0.00000028*HMOK)
+
+	if OPTIC:
+		# Refractivity from Barrell & Sears (1939)!
+		A = (287.6155+(1.62887+0.01360/WLSQ)/WLSQ)*273.15E-6/1013.25
+	else:
+		A = 77.6890E-6
+
+	GAMAL = (GB*DMD)/GCR
+	GAMMA = GAMAL/ALPHA
+	GAMM2 = GAMMA-2
+	DELM2 = DELTA-2
+	TDC   = TDKOK-273.15
+	PSAT  = 10**((0.7859+0.03477*TDC)/(1+0.00412*TDC))\
+				*(1+PMBOK*(4.5E-6+6E-10*TDC*TDC))
+
+	if PMBOK > 0:
+		PWO = RHOK*PSAT/(1-(1-RHOK)*PSAT/PMBOK)
+	else:
+		PWO = 0
+
+	W  = PWO*(1-DMW/DMD)*GAMMA/(DELTA-GAMMA)
+	C1 = A*(PMBOK+W)/TDKOK
+	if OPTIC:
+		C2 = (A*W+11.2684E-6*PWO)/TDKOK
+	else:
+		C2 = (A*W+6.3938E-6*PWO)/TDKOK
+	
+	C3 = (GAMMA-1)*ALPHA*C1/TDKOK
+	C4 = (DELTA-1)*ALPHA*C2/TDKOK
+	if OPTIC:
+		C5 = 0
+		C6 = 0
+	else:
+		C5 = 375463E-6*PWO/TDKOK
+		C6 = C5*DELM2*ALPHA/(TDKOK*TDKOK)
+
+	# Conditions at the observer.
+	R0  = S+HMOK
+	TEMPO, DN0, RDNDR0 = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,R0)
+	SK0 = DN0*R0*np.sin(ZOBS2)
+	F0  = sla_refi(DN0,RDNDR0)
+
+	# Conditions in the troposphere at the tropopause.
+	RT = S+max(HT,HMOK)
+	TT, DNT, RDNDRT = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,RT)
+	SINE = SK0/(RT*DNT)
+	ZT = np.arctan2(SINE,np.sqrt(max(1-SINE*SINE,0)))
+	FT = sla_refi(DNT,RDNDRT)
+
+	# Conditions in the stratosphere at the tropopause.
+	DNTS,RDNDRP = sla_atms(RT,TT,DNT,GAMAL,RT)
+	SINE = SK0/(RT*DNTS)
+	ZTS  = np.arctan2(SINE,np.sqrt(max(1-SINE*SINE,0)))
+	FTS  = sla_refi(DNTS,RDNDRP)
+
+	# Conditions at the stratosphere limit.
+	RS   = S+HS
+	DNS,RDNDRS = sla_atms(RT,TT,DNT,GAMAL,RS)
+	SINE = SK0/(RS*DNS)
+	ZS   = np.arctan2(SINE,np.sqrt(max(1-SINE*SINE,0)))
+	FS   = sla_refi(DNS,RDNDRS)
+
+	# Variable initialization to avoid compiler warning.
+	REFT = 0
+
+	# Integrate the refraction integral in two parts;  first in the
+	# troposphere (K=1), then in the stratosphere (K=2).
+
+	for K in [1,2]:
+		# Initialize previous refraction to ensure at least two iterations.
+		REFOLD = 1
+
+		# Start off with 8 strips.
+		IS = 8
+
+		# Start Z, Z range, and start and end values.
+		if K==1:
+			Z0 = ZOBS2
+			ZRANGE = ZT-Z0
+			FB = F0
+			FF = FT
+		else:
+			Z0 = ZTS
+			ZRANGE = ZS-Z0
+			FB = FTS
+			FF = FS
+
+		# Sums of odd and even values.
+		FO = 0
+		FE = 0
+
+		# First time through the loop we have to do every point.
+		N = 1
+
+		# Start of iteration loop (terminates at specified precision).
+		LOOP = True
+		while LOOP:
+
+			# Strip width.
+			H = ZRANGE/IS
+
+			# Initialize distance from Earth centre for quadrature pass.
+			if K==1:
+				R = R0
+			else:
+				R = RT
+
+			# One pass (no need to compute evens after first time).
+			# DO I=1,IS-1,N
+			for I in np.arange(1, IS-1, N):
+
+				# Sine of observed zenith distance.
+				SZ = np.sin(Z0+H*I)
+
+				# Find R (to the nearest metre, maximum four iterations).
+				if (SZ>1E-20):
+					W  = SK0/SZ
+					RG = R
+					DR = 1E6
+					J  = 0
+					while (abs(DR)>1 and J<4):
+						J+=1
+						if K==1:
+							TG,DN,RDNDR = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,RG)
+						else:
+							DN,RDNDR = sla_atms(RT,TT,DNT,GAMAL,RG)
+						DR = (RG*DN-W)/(DN+RDNDR)
+						RG = RG-DR
+					R = RG
+
+				# Find the refractive index and integrand at R.
+				if K==1:
+					T,DN,RDNDR = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,R)
+				else:
+					DN,RDNDR = sla_atms(RT,TT,DNT,GAMAL,R)
+
+				F = sla_refi(DN,RDNDR)
+
+				# Accumulate odd and (first time only) even values.
+				if N==1 and (I%2)==0:
+					FE = FE+F
+				else:
+					FO = FO+F
+
+
+			# Evaluate the integrand using Simpson's Rule.
+			REFP = H*(FB+4*FO+2*FE+FF)/3
+
+			# Has the required precision been achieved (or can't be)?
+			if (abs(REFP-REFOLD)>TOL and IS<ISMAX):
+				# No: prepare for next iteration.
+
+				# Save current value for convergence test.
+				REFOLD = REFP
+
+				# Double the number of strips.
+				IS = IS+IS
+
+				# Sum of all current values = sum of next pass's even values.
+				FE = FE+FO
+
+				# Prepare for new odd values.
+				FO = 0
+
+				# Skip even values next time.
+				N = 2
+			else:
+				# Yes: save troposphere component and terminate the loop.
+				if (K==1):
+					REFT = REFP
+				LOOP = False
+
+	# Result.
+	REF = REFT+REFP
+	if (ZOBS1<0): 
+		REF = -REF
+	return REF
+
+
+# =============================================================================
+# WORK IN PROGRESS; Implementation of the slalib model for 
+# 					any refractive index model.
+# Main issue is in the calculation of dn/dr for any model
+# of the refractive index. One possibility is following 
+# the Corbard error function model, which assumes the
+# atmospheric density follows the equation rho(r)=exp(-u/a).
+# But then the SLA_REFRO routine would not have to be modified
+# and one would simply use the refraction.corbard() function.
+# =============================================================================
+#
+# def get_rc(lat):
+# 		'''
+# 			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.
+
+# 			Returns
+# 			-------
+# 			rc 	: float
+# 				Local radius of curvature of the Earth in kilometers.
+# 		'''
+# 		_A   = 6378.137 # km
+# 		_B   = 6356.752 # km
+# 		_phi = np.deg2rad(lat)
+		
+# 		_rc0 = (_A*_B)**2 / (_A**2 * np.cos(_phi)**2 + _B**2 * np.sin(_phi)**2)**(3/2)
+# 		return _rc0*1000
+
+
+# def get_local_g(lat, h):
+# 	'''
+# 		Returns the local gravitic acceleration for a given latitude 
+# 		and altitude.
+
+# 		Parameters
+# 		----------
+# 		lat : float
+# 			Latitude in degrees
+# 		h 	: float
+# 			Altitude above sea level in meters
+
+# 		Returns
+# 		-------
+# 		g 	: float
+# 			Local acceleration due to gravity in m/s^2
+# 	'''
+# 	_phi = np.deg2rad(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) * h + 7.2e-13 * h**2
+# 	return _g
+
+# def aa_refro(zenith, l, h, T, p, RH, xc, lat, TLR, EPS):
+# 	'''
+# 		Modified version of sla_refro(), updating some constants
+# 		and changing the units on input and output parameters.
+
+# 		The function calculates the atmospheric refraction
+# 		optical/IR wavelengths (not radio!), ported from the original 
+# 		SLALIB FORTRAN code (Copyright P.T.Wallace.  All rights 
+# 		reserved.).
+
+# 		The method is also described in "Explanatory Supplement 
+# 		Astronomical Almanac" by P.K. Seidelmann, 1992, section 3.281.
+
+# 		See also the description of sla_refro().
+		
+# 		Parameters
+# 		----------
+# 		zenith 	:	float
+# 			Observed zenith distance of the source (degrees)
+# 		l 		:	float
+# 			Effective wavelength of the source (micrometre)
+# 		h 		:	float
+# 			Height of the observer above sea level (meter)
+# 		T 		:	float
+# 			Ambient temperature at the observer (K)
+# 		p 		:	float
+# 			Pressure at the observer (Pascal)
+# 		RH 		:	float
+# 			Relative humidity at the observer (range 0-1)
+# 		xc 		: float
+# 			CO2 density in parts per million
+# 		lat 	:	float
+# 			Latitude of the observer (degrees)
+# 		TLR 	:	float
+# 			Temperature lapse rate in the troposphere (K/metre)
+# 		EPS 	:	float
+# 			Precision required to terminate iteration (radian)
+		
+# 		Returns
+# 		-------
+# 		REF 	:	float
+# 			Refraction: in vacuo ZD minus observed ZD (radian)
+# 	'''
+	
+# 	# Converting some units for use later
+# 	ZOBS = np.deg2rad(zenith)
+# 	HM 	 = h
+# 	TDK  = T
+# 	PMB  = p / 100
+# 	WL   = l
+# 	PHI  = np.deg2rad(lat)
+
+# 	# Fixed parameters
+# 	D93 	= np.deg2rad(93)# 93 degrees in radians
+# 	GCR 	= 8314.510 		# Universal gas constant
+# 	DMD 	= 28.9635 + 12.011e-6 * (xc - 400) # Molar mass in g/mol of dry air containing xc ppm of CO2
+# 	DMW 	= 18.0152 		# Molecular weight of water vapour
+# 	S 		= get_rc(lat)	# Mean Earth radius (meters)
+# 	DELTA 	= 18.36 		# Exponent of temperature dependence of water vapour pressure
+# 	HT 		= 11000			# Height of tropopause (metre)
+# 	HS 		= 80000			# Upper limit for refractive effects (meters)
+# 	ISMAX 	= 16384			# Numerical integration: maximum number of strips.
+	
+# 	# Transform ZOBS into the normal range.
+# 	ZOBS1 = sla_drange(ZOBS)
+# 	ZOBS2 = min(abs(ZOBS1),D93)
+
+# 	# Keep other arguments within safe bounds.
+# 	HMOK  = min(max(HM,-1E3),HS)
+# 	TDKOK = min(max(TDK,100),500)
+# 	PMBOK = min(max(PMB,0),10000)
+# 	RHOK  = min(max(RH,0),1)
+# 	WLOK  = max(WL,0.1)
+# 	ALPHA = min(max(abs(TLR),0.001),0.01)
+
+# 	# Tolerance for iteration.
+# 	TOL   = min(max(abs(EPS),1E-12),0.1)/2
+
+# 	# Set up model atmosphere parameters defined at the observer.
+# 	WLSQ  = WLOK*WLOK
+
+# 	# Local gravity
+# 	GB = get_local_g(lat, h)
+
+# 	# Refractivity of air
+# 	A = (287.6155+(1.62887+0.01360/WLSQ)/WLSQ)*273.15E-6/1013.25
+
+# 	GAMAL = (GB * DMD) / GCR
+# 	GAMMA = GAMAL / ALPHA
+# 	GAMM2 = GAMMA - 2
+# 	DELM2 = DELTA - 2
+# 	TDC   = TDKOK - 273.15
+# 	PSAT  = 10**( (0.7859+0.03477*TDC) / (1+0.00412*TDC) )\
+# 				* (1 + PMBOK * (4.5E-6 + 6E-10 * TDC*TDC) )
+
+# 	if PMBOK > 0:
+# 		PWO = RHOK*PSAT/(1-(1-RHOK)*PSAT/PMBOK)
+# 	else:
+# 		PWO = 0
+
+# 	W  = PWO*(1-DMW/DMD)*GAMMA/(DELTA-GAMMA)
+# 	C1 = A*(PMBOK+W)/TDKOK
+# 	C2 = (A*W+11.2684E-6*PWO)/TDKOK
+	
+# 	C3 = (GAMMA-1)*ALPHA*C1/TDKOK
+# 	C4 = (DELTA-1)*ALPHA*C2/TDKOK
+# 	C5 = 0
+# 	C6 = 0
+
+# 	# Conditions at the observer.
+# 	R0  = S+HMOK
+# 	TEMPO, DN0, RDNDR0 = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,R0)
+# 	SK0 = DN0*R0*np.sin(ZOBS2)
+# 	F0  = sla_refi(DN0,RDNDR0)
+
+# 	# Conditions in the troposphere at the tropopause.
+# 	RT = S+max(HT,HMOK)
+# 	TT, DNT, RDNDRT = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,RT)
+# 	SINE = SK0/(RT*DNT)
+# 	ZT = np.arctan2(SINE,np.sqrt(max(1-SINE*SINE,0)))
+# 	FT = sla_refi(DNT,RDNDRT)
+
+# 	# Conditions in the stratosphere at the tropopause.
+# 	DNTS,RDNDRP = sla_atms(RT,TT,DNT,GAMAL,RT)
+# 	SINE = SK0/(RT*DNTS)
+# 	ZTS  = np.arctan2(SINE,np.sqrt(max(1-SINE*SINE,0)))
+# 	FTS  = sla_refi(DNTS,RDNDRP)
+
+# 	# Conditions at the stratosphere limit.
+# 	RS   = S+HS
+# 	DNS,RDNDRS = sla_atms(RT,TT,DNT,GAMAL,RS)
+# 	SINE = SK0/(RS*DNS)
+# 	ZS   = np.arctan2(SINE,np.sqrt(max(1-SINE*SINE,0)))
+# 	FS   = sla_refi(DNS,RDNDRS)
+
+# 	# Variable initialization to avoid compiler warning.
+# 	REFT = 0
+
+# 	# Integrate the refraction integral in two parts;  first in the
+# 	# troposphere (K=1), then in the stratosphere (K=2).
+# 	for K in [1,2]:
+# 		# Initialize previous refraction to ensure at least two iterations.
+# 		REFOLD = 1
+
+# 		# Start off with 8 strips.
+# 		IS = 8
+
+# 		# Start Z, Z range, and start and end values.
+# 		if K==1:
+# 			Z0 = ZOBS2
+# 			ZRANGE = ZT-Z0
+# 			FB = F0
+# 			FF = FT
+# 		else:
+# 			Z0 = ZTS
+# 			ZRANGE = ZS-Z0
+# 			FB = FTS
+# 			FF = FS
+
+# 		# Sums of odd and even values.
+# 		FO = 0
+# 		FE = 0
+
+# 		# First time through the loop we have to do every point.
+# 		N = 1
+
+# 		# Start of iteration loop (terminates at specified precision).
+# 		LOOP = True
+# 		while LOOP:
+
+# 			# Strip width.
+# 			H = ZRANGE/IS
+
+# 			# Initialize distance from Earth centre for quadrature pass.
+# 			if K==1:
+# 				R = R0
+# 			else:
+# 				R = RT
+
+# 			# One pass (no need to compute evens after first time).
+# 			# DO I=1,IS-1,N
+# 			for I in np.arange(1, IS-1, N):
+
+# 				# Sine of observed zenith distance.
+# 				SZ = np.sin(Z0+H*I)
+
+# 				# Find R (to the nearest metre, maximum four iterations).
+# 				if (SZ>1E-20):
+# 					W  = SK0/SZ
+# 					RG = R
+# 					DR = 1E6
+# 					J  = 0
+# 					while (abs(DR)>1 and J<4):
+# 						J+=1
+# 						if K==1:
+# 							TG,DN,RDNDR = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,RG)
+# 						else:
+# 							DN,RDNDR = sla_atms(RT,TT,DNT,GAMAL,RG)
+# 						DR = (RG*DN-W)/(DN+RDNDR)
+# 						RG = RG-DR
+# 					R = RG
+
+# 				# Find the refractive index and integrand at R.
+# 				if K==1:
+# 					T,DN,RDNDR = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,R)
+# 				else:
+# 					DN,RDNDR = sla_atms(RT,TT,DNT,GAMAL,R)
+
+# 				F = sla_refi(DN,RDNDR)
+
+# 				# Accumulate odd and (first time only) even values.
+# 				if N==1 and (I%2)==0:
+# 					FE = FE+F
+# 				else:
+# 					FO = FO+F
+
+
+# 			# Evaluate the integrand using Simpson's Rule.
+# 			REFP = H*(FB+4*FO+2*FE+FF)/3
+
+# 			# Has the required precision been achieved (or can't be)?
+# 			if (abs(REFP-REFOLD)>TOL and IS<ISMAX):
+# 				# No: prepare for next iteration.
+
+# 				# Save current value for convergence test.
+# 				REFOLD = REFP
+
+# 				# Double the number of strips.
+# 				IS = IS+IS
+
+# 				# Sum of all current values = sum of next pass's even values.
+# 				FE = FE+FO
+
+# 				# Prepare for new odd values.
+# 				FO = 0
+
+# 				# Skip even values next time.
+# 				N = 2
+# 			else:
+# 				# Yes: save troposphere component and terminate the loop.
+# 				if (K==1):
+# 					REFT = REFP
+# 				LOOP = False
+
+# 	# Result.
+# 	REF = REFT+REFP
+# 	if (ZOBS1<0): 
+# 		REF = -REF
+# 	return REF
+
+
+

AstroAtmosphere/observatories.py

@@ -0,0 +1,85 @@
+observatories = { 
+					'STANDARD': {
+						'T': 288.15, 
+						'p': 101325,
+						'RH': 0.0,
+						'xc': 450,
+						'latitude': 0,
+						'altitude': 0,
+						'comment': 'Standard atmospheric conditions'
+						},
+					'CERRO_ARMAZONES': {
+						'T': 279.65, 
+						'p': 71200,
+						'RH': 0.22,
+						'xc': 450,
+						'latitude': -24.5983,
+						'altitude': 3064,
+						'comment': 'Default conditions for Cerro Armazones, site of the ELT.\nValues are taken from the optical model of MICADO.'
+						},
+					'CERRO_PARANAL': {
+						'T': 285.95, 
+						'p': 74300,
+						'RH': 0.15,
+						'xc': 450,
+						'latitude': -24.6230,
+						'altitude': 2635,
+						'comment': 'Average conditions at Cerro Paranal, site of the VLT.\nValues are taken from https://www.eso.org/sci/facilities/paranal/astroclimate/site.html and https://academic.oup.com/mnras/article/399/2/783/1061377'
+						},
+					'LA_SILLA': {
+						'T': 284.65, 
+						'p': 77220,
+						'RH': 0.44,
+						'xc': 450,
+						'latitude': -29.25,
+						'altitude': 2400,
+						'comment': 'Average conditions of the La Silla Observatory.\nValues based on https://www.eso.org/gen-fac/pubs/astclim/lasilla/diffrefr.html'
+						},
+					'LAS_CAMPANAS': {
+						'T': 285.85, 
+						'p': 76400,
+						'RH': 0.0,
+						'xc': 450,
+						'latitude': -29.016,
+						'altitude': 2380,
+						'comment': 'Mean conditions at Las Campanas observatory, site of the GMT.\nValues based on https://www.gmto.org/Resources/GMT-ID-01466-Chapter_5_Site_Evaluation.pdf and http://weather.lco.cl/'
+						},
+					'LA_PALMA': {
+						'T': 281.95, 
+						'p': 77500,
+						'RH': 0.36,
+						'xc': 450,
+						'latitude': 28.764,
+						'altitude': 2396,
+						'comment': 'Average conditions for Roque de los Muchachos Observatory on La Palma, Spain.\nValues taken from https://academic.oup.com/mnras/article/399/2/783/1061377'
+						},
+					'MAUNA_KEA': {
+						'T': 274.45, 
+						'p': 61400,
+						'RH': 0.48,
+						'xc': 450,
+						'latitude': 19.826,
+						'altitude': 4205,
+						'comment': 'Average conditions at the Mauna Kea Observatory.\nValues taken from http://www.soest.hawaii.edu/MET/Faculty/businger/mauna_kea/daSilvaThesisFinal.pdf, Ch. 3. and https://www.eso.org/gen-fac/pubs/astclim/espas/espas_reports/ESPAS-MaunaKea.pdf'
+						}
+					# 'MT_GRAHAM': {
+					# 	'T': 278.15, 
+					# 	'p': 68500,
+					# 	'RH': 0.0,
+					# 	'xc': 450,
+					# 	'latitude': 32.701,
+					# 	'altitude': 3191,
+					# 	'comment': 'Mean conditions at Mount Graham International Observatory,\
+					# 				site of the LBT. Values are estimates.'
+					# 	},
+					# 'FRED_LAWRENCE_WHIPPLE': {
+					# 	'T': 288.15, 
+					# 	'p': 74300,
+					# 	'RH': 0.0,
+					# 	'xc': 450,
+					# 	'latitude': 31.681,
+					# 	'altitude': 2616,
+					# 	'comment': 'Mean conditions at Fred Lawrence Whipple Observatory, site of\
+					# 				the MMT. Values are estimates.'
+					#	}
+					}

AstroAtmosphere/refractionModels.py

@@ -9,11 +9,11 @@ from .ciddorModel import Observatory
 @ 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
+				index of refraction for a wide range of conditions, with corresponding error propagation if needed.
-				error propagation if needed.
 				The actual refraction is calculated by any of the following models:
 				- Plane parallel atmosphere
 				- Cassini spherical homogeneous atmosphere
+				- Mathar's Barometric exponential model
 				- Other methods described in Corbard et al., 2019.
 '''
 
@@ -38,39 +38,45 @@ class refraction(Observatory):
 
 		Methods
 		-------
-		setReducedHeight(p, rho)
-			Sets the reduced atmospheric height.
-		planeParallel(n1, n2, zenith)
-			Returns the atmospheric refraction assuming a plane parallel
-			atmosphere
 		cassini(n1, n2, zenith, reduced_height=None)
-			Returns the atmospheric refraction asusming Cassini's homegeneous
+			Returns the atmospheric refraction assuming Cassini's homegeneous
-			atmosphere model
+			atmosphere model.
-		tan5(n1, n2, zenith, reduced_height=None)
-			Returns the atmospheric refraction following Oriani's theorem,
-			including a fifth order term in the expansion.
-		oriani(n1, n2, zenith, reduced_height=None)
-			Returns the atmospheric refraction following the classic version
-			of Oriani's theorem.
-		corbard(n1, n2, zenith, reduced_height=None)
-			Returns the atmospheric refraction following the method described 
-			in Corbard et al., 2019.
-		refractionIntegral(n1, n2, z0, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None)
-			Returns the atmospheric refraction following the refraction integral.
-		psi(x)
-			Submethod for corbard()
-		set_rc()
-			Sets the local radius of curvature of the Earth at the observer.
-		set_H(p, rho)
-			Sets the reduced height for the conditions at the observer, assuming
-			an adiabatic atmosphere.
-		H_isotherm(T)
-			Sets the reduced height for the conditions at the observer, assuming
-			an isothermal atmosphere.
 		cassiniError(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)
 			Returns the uncertainty in the atmospheric refraction, assuming 
 			Cassini's homegeneous atmosphere model and the given parameters.
+		corbard(n1, n2, zenith, reduced_height=None)
+			Returns the atmospheric refraction following the method described 
+			in Corbard et al., 2019.
+		gammaCoefficients(n1, reduced_height=None)
+			Helper function for matharExponential(), that calculates the gamma
+			coefficients.
+		H_isotherm(T)
+			Sets the reduced height for the conditions at the observer, assuming
+			an isothermal atmosphere.
+		matharExponential(n1, zenith, reduced_height=None)
+			Returns the atmospheric refraction assuming	Mathar's barometric 
+			exponential model.
+		oriani(n1, n2, zenith, reduced_height=None)
+			Returns the atmospheric refraction following the classic version
+			of Oriani's theorem.	
+		planeParallel(n1, n2, zenith)
+			Returns the atmospheric refraction assuming a plane parallel
+			atmosphere.
+		psi(x)
+			Submethod for corbard().
+		refractionIntegral(n1, n2, zenith, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None)
+			Returns the atmospheric refraction following the refraction integral.
+		setReducedHeight(p, rho)
+			Sets the reduced atmospheric height.
+		set_H(p, rho)
+			Sets the reduced height for the conditions at the observer, assuming
+			an adiabatic atmosphere.
+		set_rc()
+			Sets the local radius of curvature of the Earth at the observer.			
+		tan5(n1, n2, zenith, reduced_height=None)
+			Returns the atmospheric refraction following Oriani's theorem,
+			including a fifth order term in the expansion.				
 		Other Observatory() methods
 	'''
 	def __init__(self, lat, h):
@@ -80,41 +86,6 @@ class refraction(Observatory):
 		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):
-		'''
-			Sets the reduced atmospheric height.
-
-			Parameters
-			----------
-			p 	: float
-				Atmospheric pressure in Pascal
-			rho : float
-				Density of the atmospheric air at the observer
-		'''
-		self.H  	= self.set_H(p, rho) 		# Reduced height of the atmosphere assuming ideal gas law
-
-
-	def planeParallel(self, n1, zenith):
-		'''
-			Refraction based on a flat atmosphere, easily derived from Snell's law.
-
-			Parameters
-			----------
-			n1 : float
-				Refractive index at the wavelength of interest
-			zenith : float
-				Angle of observation in degrees.
-				zenith = 0 means the observation is directly overhead.
-
-			Returns
-			-------
-			R : float
-				Atmospheric refraction in degrees
-		'''
-		_R = np.arcsin(n * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
-		return np.degrees(_R)
-
-
 
 	def cassini(self, n1, zenith, reduced_height=None):
 		'''
@@ -145,276 +116,6 @@ class refraction(Observatory):
 		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.
-
-			Parameters
-			----------
-			n1 : float
-				Refractive index at the wavelength of interest
-			zenith : float
-				Angle of observation in degrees.
-				zenith = 0 means the observation is directly overhead.
-			reduced_height : float (optional)
-				Reduced height of the atmosphere. Default value is calculated
-				from the object attributes.
-
-			Returns
-			-------
-			R : float
-				Atmospheric refraction in degrees
-		'''
-		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.
-
-			Parameters
-			----------
-			n1 : float
-				Refractive index at the wavelength of interest
-			zenith : float
-				Angle of observation in degrees.
-				zenith = 0 means the observation is directly overhead.
-			reduced_height : float (optional)
-				Reduced height of the atmosphere. Default value is calculated
-				from the object attributes.
-
-			Returns
-			-------
-			R : float
-				Atmospheric refraction in degrees
-		'''
-		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'.
-
-			Parameters
-			----------
-			n1 : float
-				Refractive index at the wavelength of interest
-			zenith : float
-				Angle of observation in degrees.
-				zenith = 0 means the observation is directly overhead.
-			reduced_height : float (optional)
-				Reduced height of the atmosphere. Default value is calculated
-				from the object attributes.
-
-			Returns
-			-------
-			R : float
-				Atmospheric refraction in degrees
-		'''
-		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, n1, 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.
-
-			Parameters
-			----------
-			n1 : float
-				Refractive index at the wavelength of interest
-			z0 : float
-				Angle of observation in degrees.
-				zenith = 0 means the observation is directly overhead.
-			R0 : float (optional)
-				Local radius of curvature of the Earth in km.
-			useStandardAtmosphere : bool
-				Flag to use the included USSA1976 data
-			heightData : 1D Array
-				If "useStandardAtmosphere" is False, the function will ask
-				for alternative atmospheric data to integrate over. This
-				gives the altitudinal points
-			rhoData : 1D Array
-				If useStandardAtmosphere is False, the function will ask
-				for alternative atmospheric data to integrate over. This
-				gives the atmospheric density at the points defined by 
-				"heightData".
-
-			Returns
-			-------
-			R : float
-				Atmospheric refraction in degrees
-
-		'''
-		if not R0:
-			R0 = self.rc
-		if useStandardAtmosphere:
-			_datafileLocation = os.path.join(os.path.dirname(__file__), 'data/1976USSA.txt')
-			_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt(_datafileLocation, 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):
-		'''
-			Submethod for corbard(). Described explicitly in Corbard et al.,2019.
-
-			Parameters
-			----------
-			x 	: 	float
-				Input of the function
-
-			Returns
-			-------
-			y 	: 	float
-				Output of the function
-		'''
-		_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.
-
-			Returns
-			-------
-			rc 	: float
-				Local radius of curvature of the Earth in kilometers.
-		'''
-		_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.
-
-			Parameters
-			----------
-			p 	: float
-				Atmospheric pressure in Pascal
-			rho : float
-				Density of the atmospheric air at the observer
-
-			Returns
-			-------
-			H 	: float
-				Reduced height of the atmosphere in meters
-		'''
-		_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.
-
-			Parameters
-			----------
-			T 	: float
-				Atmospheric temperature in Kelvin
-
-			Returns
-			-------
-			H 	: float
-				Reduced height of the atmosphere in meters
-		'''
-		_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.
@@ -455,3 +156,426 @@ class refraction(Observatory):
 		_dR 	= np.sqrt( (_dRdn * dn)**2 + (_dRdz * _dz)**2)
 
 		return np.degrees(_dR)
+
+
+	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'.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			R : float
+				Atmospheric refraction in degrees
+		'''
+		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 gammaCoefficients(self, n1, reduced_height=None):
+		'''
+			Submethod for matharExponential(). Calculates the gamma
+			coefficients.
+
+			Parameters
+			----------
+			n1 	: 	float
+				Refractive index at the wavelength of interest at the observer
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			g 	: 	array
+				List of values, corresponding to the gamma coefficients.
+				g[0] corresponds to g_1, g[1] to g_3, g[2] to g_5, and so on.
+		'''
+
+		if reduced_height:
+			H = reduced_height
+		else:
+			H = self.H
+
+		_chi0   = n1**2 - 1
+		_Khat 	= H / self.rc
+
+		_g1 = (1 - 5/24 * _Khat * _chi0**2 - 3/4 * _chi0 - 3/8 * _Khat**2 * _chi0 \
+			- _Khat - 6 * _Khat**3 + 5/8 * _chi0**2 + 63/128 * _chi0**4 + 3/8 \
+			* _Khat * _chi0 + 35/256 * _Khat * _chi0**3 + 2 * _Khat**2 + 5/36 \
+			* _Khat**2 * _chi0**2 - 35/64 * _chi0**3 + 9/16 * _Khat**3 * _chi0 \
+			- 120 * _Khat**5 + 24 * _Khat**4 - 231 / 512 * _chi0**5 - 63 / 640 \
+			* _Khat * _chi0**4 - 35 / 512 * _Khat**2 * _chi0**3 - 5 / 36 * _Khat**3\
+			* _chi0**2 - 9 / 8 * _Khat**4 * _chi0 )
+		
+		_g3 = (35/96 * _chi0**3 - 5/12 * _chi0**2 + 333 * _Khat**4 * _chi0 - 2400 \
+			* _Khat**5 + 3605/4608 * _Khat**2 * _chi0**3 + 63/256 * _Khat * _chi0**4\
+			 - 819/16 * _Khat**3 * _chi0 - 35/18 * _Khat**2 * _chi0**2 - 35/96 \
+			 * _Khat * _chi0**3 + 69/8 * _Khat**2 * _chi0 + 5/8 * _Khat * _chi0**2 \
+			 - 3/2 * _Khat * _chi0 + 1355/216 * _Khat**3 * _chi0**2 + 10 * _Khat**2 \
+			 + 336 * _Khat**4 - 2 * _Khat + _chi0/2 - 21/64 * _chi0**4 - 54 * _Khat**3 \
+			 + 77/256 * _chi0**5 ) * 1/2
+		
+		_g5 = (-96 * _Khat**3 + 8 * _Khat**2 - 7/24 * _chi0**3 + 65/4 * _Khat**2 \
+			* _chi0**2 + 35/48 * _Khat * _chi0**3 + 30 * _Khat**2 * _chi0 - 5/3 \
+			* _Khat * _chi0**2 - 3 * _Khat * _chi0 + 20091/8 * _Khat**4 * _chi0 \
+			- 32605/216 * _Khat**3 * _chi0**2 - 35/9 * _Khat**2 * _chi0**3 + 984 \
+			* _Khat**4 + _chi0**2/3 - 10200 * _Khat**5 - 77/320 * _chi0**5 + 21/80 \
+			* _chi0**4 - 1089/4 * _Khat**3 * _chi0 - 7/16 * _Khat * _chi0**4 ) * 3/8
+		
+		_g7 = (-17520 * _Khat**5 - 9/40 * _chi0**4 - 8755/12 * _Khat**3 * _chi0**2 \
+			+ 1056 * _Khat**4 + _chi0**3/4 + 23517/4 * _Khat**4 * _chi0 - 11/3 \
+			* _Khat * _chi0**2 + 21 * _Khat**2 * _chi0 - 7/4 * _Khat * _chi0**3 \
+			+ 505/9 * _Khat**2 * _chi0**2 - 801/2 * _Khat**3 * _chi0 + 33/160 \
+			* _chi0**5 - 48 * _Khat**3 + 63/80 * _Khat * _chi0**4 + 2345/96 \
+			* _Khat**2 * _chi0**3) * 5/16
+		
+		_g9 = (_chi0**4 / 5 + 384 * _Khat**4 - 13440 * _Khat**5 + 5562 * _Khat**4 \
+			* _chi0 + 3125/36 * _Khat**2 * _chi0**3 - 11/60 * _chi0**5 - 25/6 \
+			* _Khat * _chi0**3 - 26950/27 * _Khat**3 * _chi0**2 + 340/9 * _Khat**2 \
+			* _chi0**2 - 180 * _Khat**3 * _chi0 - 9/5 * _Khat * _chi0**4) * 35/128
+
+		g = np.array([_g1, _g3, _g5, _g7, _g9])
+		g*= n1*_chi0/2
+
+		return g		
+	
+
+	def H_isotherm(self, T):
+		'''
+			Calculates the reduced height of the atmosphere, 
+			assuming an isotherm distribution.
+
+			Parameters
+			----------
+			T 	: float
+				Atmospheric temperature in Kelvin
+
+			Returns
+			-------
+			H 	: float
+				Reduced height of the atmosphere in meters
+		'''
+		_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 matharExponential(self, n1, zenith, reduced_height=None):
+		'''
+			Model on atmospheric refraction described by Richard Mathar in 
+			"A Barometric Exponential Model of the Atmosphere’s Refractive Index:  ZenithAngles and Second Order Aberration in the Entrance Pupil" (2020)
+			ArXiv id: 2004.11808 (https://arxiv.org/abs/2004.11808)
+			
+			The code in this function is adapted from eqs 3, 16, 17 and 18 in this
+			work. It is an adaptation of the RefrExpo.cxx C++ file written by Mathar.
+
+			This model assumes the refractive index, n, is related to the susceptibility, 
+			chi, following: n(r) = sqrt(1+chi(r)) (eq. 9).
+
+			R. Mathar then assumes that the susceptibility is proportional to the
+			air density and reduced height, following eq. 10:
+			chi(r) = chi0 * exp(-(r-rho)/K).
+
+			The refraction is described using a Taylor expansion; R in terms of
+			the odd powers of the tangent of the observed zenith angle (eq. 3), not
+			unlike Oriani's theorem. This requires coefficients g_i.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the wavelength of interest at the observer
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			R : float
+				Atmospheric refraction in degrees
+		'''
+
+		if reduced_height:
+			_H = reduced_height
+		else:
+			_H = self.H
+
+		_zenith	= np.deg2rad(zenith)
+		_rc 	= self.rc # rho
+
+		# Atmospheric refraction in radians
+		_gamma  = self.gammaCoefficients(n1, reduced_height=_H)
+		_R = sum([_gamma[i]*np.tan(_zenith)**(2*i+1) for i in range(len(_gamma))])
+		
+		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.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			R : float
+				Atmospheric refraction in degrees
+		'''
+		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 planeParallel(self, n1, zenith):
+		'''
+			Refraction based on a flat atmosphere, easily derived from Snell's law.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+
+			Returns
+			-------
+			R : float
+				Atmospheric refraction in degrees
+		'''
+		_R = np.arcsin(n1 * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
+		return np.degrees(_R)
+
+
+	def psi(self, x):
+		'''
+			Submethod for corbard(). Described explicitly in Corbard et al.,2019.
+
+			Parameters
+			----------
+			x 	: 	float
+				Input of the function
+
+			Returns
+			-------
+			y 	: 	float
+				Output of the function
+		'''
+		_f   = lambda a: np.exp(-a**2)
+		_int = integrate.quad(_f, x, np.inf)
+		return np.exp(x**2) * _int[0]
+
+
+	def refractionIntegral(self, n1, zenith, 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.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the wavelength of interest at the observer
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			R0 : float (optional)
+				Local radius of curvature of the Earth in km.
+			useStandardAtmosphere : bool
+				Flag to use the included USSA1976 data
+			heightData : 1D Array
+				If "useStandardAtmosphere" is False, the function will ask
+				for alternative atmospheric data to integrate over. This
+				gives the altitudinal points
+			rhoData : 1D Array
+				If useStandardAtmosphere is False, the function will ask
+				for alternative atmospheric data to integrate over. This
+				gives the atmospheric density at the points defined by 
+				"heightData".
+
+			Returns
+			-------
+			R : float
+				Atmospheric refraction in degrees
+
+		'''
+		if not R0:
+			R0 = self.rc
+		if useStandardAtmosphere:
+			_datafileLocation = os.path.join(os.path.dirname(__file__), 'data/1976USSA.txt')
+			_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt(_datafileLocation, unpack=True, skip_header=6)
+		else:
+			_Dens = rhoData
+			_h 	  = heightData
+
+		zenith 		= np.deg2rad(zenith)
+		_Dens 	= np.asarray(_Dens)
+		_h 	 	= np.asarray(_h) 
+		_R 		= R0 + _h * 1000
+		
+		# Gladstone-Dale relation between refractive index and density
+		_n_h = 1 + (n1 - 1) / _Dens[0] * _Dens 
+
+		# Calculation of the zenith angle for given height and refractive index, from the refractive invariant
+		_z 		= np.arcsin( (n1 * R0) / (_n_h * _R) * np.sin(zenith) )
+		
+		# 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 setReducedHeight(self, p, rho):
+		'''
+			Sets the reduced atmospheric height.
+
+			Parameters
+			----------
+			p 	: float
+				Atmospheric pressure in Pascal
+			rho : float
+				Density of the atmospheric air at the observer
+		'''
+		self.H  	= self.set_H(p, rho) 		# Reduced height of the atmosphere assuming ideal gas law
+	
+
+	def set_H(self, p, rho):
+		'''
+			Calculates the reduced height of the atmosphere, assuming ideal gas law.
+
+			Parameters
+			----------
+			p 	: float
+				Atmospheric pressure in Pascal
+			rho : float
+				Density of the atmospheric air at the observer
+
+			Returns
+			-------
+			H 	: float
+				Reduced height of the atmosphere in meters
+		'''
+		_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 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.
+
+			Returns
+			-------
+			rc 	: float
+				Local radius of curvature of the Earth in kilometers.
+		'''
+		_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 tan5(self, n1, zenith, reduced_height=None):
+		'''
+			Same as oriani's theorem, but including a higher order term.
+			Found in Corbard et al., 2019.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			R : float
+				Atmospheric refraction in degrees
+		'''
+		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)
+

AstroAtmosphere/refractivityModels.py

@@ -9,7 +9,7 @@ import numpy as np
 				atmospheric dispersion, given a set of conditions.
 '''
 
-def BarrellAndSears(l, T=288.15, p=101325, H=0):
+def BarrellAndSears(l, T=288.15, p=101325, RH=0, f=None):
 	'''
 		From Barrell & Sears, 1939
 		Equation 7.7, page 52
@@ -24,20 +24,26 @@ def BarrellAndSears(l, T=288.15, p=101325, H=0):
 			Temperature in Kelvin
 		p 	: 	float
 			Atmospheric pressure in Pascal
-		H 	: 	float
+		RH 	: 	float
 			Relative humidity (0=<H=<1)
+		f 	: 	float
+			Water vapour pressure in Pascal
 
 		Returns
 		-------
 		n 	: 	float
 			Refractive index of atmospheric air at the given conditions
 	'''
+	if f is None:
+		f = RH * svp(T)
+	f_mmHg = hPa2mmHg(f/100)
 	T = K2C(T)
-	P = hPA2mmHg(p/100)
+	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)
+	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) * f_mmHg / (1 + 0.003661*T)
 	n = R/1e6 + 1
 	return n
 
+
 def Ciddor(l, CO2=450):
 	'''
 		From Ciddor, 1996
@@ -49,6 +55,8 @@ def Ciddor(l, CO2=450):
 		---------- 
 		l 	: 	float
 			Wavelength in microns
+		CO2 : 	float
+			CO2 density in parts per million
 
 		Returns
 		-------
@@ -62,17 +70,28 @@ def Ciddor(l, CO2=450):
 	n = R/1e8 + 1
 	return n
 
-def Edlen1953(l):
+
+def Edlen1953(l, T=288.15, p=101325, RH=0, f=None):
 	'''
 		From Edlen, 1953
 		Equation 2
 
-		At standard conditions only
+		At standard conditions only.
+		Following Filippenko (1982), eqs. 2 & 3 we also include
+		the	modification proposed by Barrell in 1951.
 
 		Parameters
 		---------- 
 		l 	: 	float
 			Wavelength in microns
+		T 	: 	float
+			Temperature in Kelvin
+		p 	: 	float
+			Atmospheric pressure in Pascal
+		RH 	: 	float
+			Relative humidity (0=<H=<1)
+		f 	: 	float
+			Water vapor pressure in Pascal
 
 		Returns
 		-------
@@ -81,13 +100,27 @@ def Edlen1953(l):
 	'''
 	s = 1 / l
 	R = 6432.8 + 2949810 / (146 - s**2) + 25540/(41-s**2)
-	n = R/1e8 + 1
+	n_standard = R/1e8 + 1
+
+	pHg = hPa2mmHg(p/100)
+	Tc  = K2C(T)
+	if f is None:
+		f = RH * svp(T)
+	f_mmHg = hPa2mmHg(f/100)
+
+	# Changes in temperature and pressure, from Barrell (1951)
+	n_tp = R/1e8 * (pHg * (1 + (1.049 - 0.0157*Tc) * 1e-6 * pHg ) ) / ( 720.883 * (1 + 0.003661 * Tc)) + 1
+
+	# Correction for water vapour pressure
+	n_svp = (n_tp - 1) * 1e6 - ( 0.0624 - 0.000680 * s**2 ) / (1 + 0.003661 * Tc) * f_mmHg
+	n = n_svp/1e6 + 1
 	return n
 
-def Edlen1966(l):
+
+def Edlen1966(l, T=288.15, p=101325, RH=0, CO2=300, f=None):
 	'''
 		From Edlen, 1966
-		Equation 1
+		Equations 1, 12, 17, 22
 
 		At standard conditions only
 
@@ -95,6 +128,16 @@ def Edlen1966(l):
 		---------- 
 		l 	: 	float
 			Wavelength in microns
+		T 	: 	float
+			Temperature in Kelvin
+		p 	: 	float
+			Atmospheric pressure in Pascal
+		RH 	: 	float
+			Relative humidity (0=<H=<1)
+		CO2 : 	float
+			CO2 density in parts per million
+		f 	: 	float
+			Water vapor pressure in Pascal
 
 		Returns
 		-------
@@ -104,14 +147,37 @@ def Edlen1966(l):
 	s = 1 / l
 	R = 8342.13 + 2406030 / (130 - s**2) + 15997 / (38.9 - s**2)
 	n = R/1e8 + 1
+
+	# Influence of temperature and pressure, eq. 12
+	Tc = K2C(T)
+	Ptorr = Pa2Torr(p)
+	R_tp = ( Ptorr * R/1e8 ) / 720.775 * (1 + Ptorr * (0.817 - 0.0133 * Tc) * 1e-6) / ( 1 + 0.0036610 * Tc)
+
+	# Allow for a CO2 fraction, eq. 17, but we include T and P dependence.
+	# This is okay from "Asthe refractivity of CO, is about 50% higher than
+	# that of air this CO,-content increases the refractivity over that of 
+	# CO2-free air by approximately 0.5 x 0.0003 x (n-1)s" and "As the 
+	# density factor of CO, is not much different from that of air this ratio r
+	# can be taken as independent of t and p ."
+	R_x = R_tp * (1 + 0.540 * (CO2/1e6 - 0.0003))
+	n_x = R_x + 1
+
+	# Correction for water vapour, eq. 22
+	if f is None:
+		f = RH * svp(T)
+	ftorr = Pa2Torr(f)
+
+	n = n_x - ftorr * (5.7224 - 0.0457 * s**2) * 1e-8
 	return n
 
-def Owens(l, T=288.15, p=101325, H=0):
+
+def Owens(l, T=288.15, p=101325, RH=0):
 	'''
 		From Owens, 1967
-		Equation 32
+		Equations 29 to 31
-		Note that the pressure P and H are in torricelli (Torr),
+		
-		however, this is almost precisely 1 Torr = 1 mmHg
+		Code heavily based on ScopeSIM (https://scopesim.readthedocs.io/)
+		and http://www.caha.es/newsletter/news03b/pedraz/newslet.html
 
 		Parameters
 		----------
@@ -121,25 +187,43 @@ def Owens(l, T=288.15, p=101325, H=0):
 			Temperature in Kelvin
 		p 	: 	float
 			Atmospheric pressure in Pa
-		H 	: 	float
+		RH 	: 	float
-			Humdity (0<=H<=1)
+			Relative humdity (0<=H<=1)
 
 		Returns
 		-------
 		n 	: 	float
 			Refractive index of atmospheric air at the given conditions
 	'''
-	P = hPA2mmHg(p/100)
+	# Calculation from ScopeSIM.effects.atmospheric_refraction()
-	T = K2C(T)
+	pres = p / 100 # from mbar to Pa
-	s = 1 / l
+	Ps = -10474. + 116.43 * T - 0.43284 * T**2 + 0.0005384 * T**3
-	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)
+	Pw = Ps * RH / 100.
+	Pa = pres - Pw
+
+	# Wavenumber s
+	s = 1. / l
+
+	# Density factor of dry air, eq. 30
+	Da = (Pa / T) * (1. + Pa * (57.9E-8 - 0.0009325 / T + 0.25844 / T**2))
+	
+	# Density factor of moist air, eq. 31
+	Dw = (Pw / T) * (1. + Pw * (1. + 3.7E-4 * Pw) *
+					 (-2.37321E-3 + 2.23366 / T - 710.792 / T**2 +
+					  77514.1 / T**3))
+
+	# Combining the two into eq. 29
+	na = Da * (2371.34 + 683939.7 / (130. - s**2) + 4547.3 / (38.9 - s**2))
+	nw = Dw * (6487.31 + 58.058 * s**2 - 0.7115 * s**4 + 0.08851 * s**6)
+	R = na + nw
 	n = R/1e8 + 1
 	return n
 
-def BonschPotulski(l, T=288.15, p=101325, CO2=450):
+
+def BonschPotulski(l, T=288.15, p=101325, RH=0, CO2=400, f=None):
 	'''
 		From Bonsch & Potulski, 1998
-		Equations 6, 7, 8
+		Equations 6, 7, 8 and 9a
 		
 		Parameters
 		---------- 
@@ -149,8 +233,12 @@ def BonschPotulski(l, T=288.15, p=101325, CO2=450):
 			Temperature in Kelvin
 		p 	: 	float
 			Atmospheric pressure in Pa
+		RH 	: 	float
+			Relative humidity (0<=H<=1)
 		CO2 : 	float
 			Particles per million of CO2
+		f 	: 	float
+			Water vapor pressure in Pascal
 
 		Returns
 		-------
@@ -158,20 +246,26 @@ def BonschPotulski(l, T=288.15, p=101325, CO2=450):
 			Refractive index of atmospheric air at the given conditions
 	'''
 	x = CO2 / 1e6
-	T = K2C(T)
+	Tc = 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)
+	Rtp = Rx * p / 93214.60 * (1 + 1e-8 * (0.5953 - 0.009876 * Tc)*p) / (1 + 0.0036610 * Tc)
-	n = Rtp + 1
+
+	if f is None:
+		f = RH * svp(T)
+	
+	n_tp = Rtp + 1
+	n = n_tp - f * (3.8020 - 0.0384 * s**2) * 1e-10
 	return n
 
-def BirchDowns(l):
+
+def BirchDowns(l, T=288.15, p=101325, RH=0, f=None):
 	'''
 		From Birch & Downs, 1993 and 1994 (correction)
-		Equations 1 and 2
+		Equations 1, 2 and 3
 
 		Only at standard conditions
 
@@ -179,17 +273,35 @@ def BirchDowns(l):
 		----------
 		l 	: 	float
 			wavelength in micron
+		T 	: 	float
+			Temperature in Kelvin
+		p 	: 	float
+			Atmospheric pressure in Pa
+		RH 	: 	float
+			Relative humidity (0<=H<=1)
+		f 	: 	float
+			Water vapor pressure in Pascal
 
 		Returns
 		-------
 		n 	: 	float
 			Refractive index of atmospheric air at the given conditions
 	'''
+	Tc = K2C(T)
 	s = 1 / l
-	R = 8342.54 + 2406147 / (130 - s**2) + 15998 / (38.9-s**2)
+	Rs = 8342.54 + 2406147 / (130 - s**2) + 15998 / (38.9-s**2)
-	n = R/1e8 + 1
+	n = Rs/1e8 + 1
+
+	R_tp = (p * Rs/1e8) / 96095.43 * (1 + 1e-8 * (0.601 - 0.00972 * Tc) * p) / (1 + 0.0036610 * Tc) 
+	n_tp = R_tp + 1
+	
+	if f is None:
+		f = RH * svp(T)
+
+	n = n_tp - f * (3.7345 - 0.0401 * s**2) * 1e-10
 	return n
 
+
 def PeckReeder(l):
 	'''
 		From Peck & Reeder, 1972
@@ -213,8 +325,70 @@ def PeckReeder(l):
 	return n
 
 
-def slalib(l, T=288.15, p=101325, H=0, h=0, lat=-24, tlr=0.0065, eps=1e-10):
+def HohenkerkAndSinclair(l, T=288.15, p=101325, RH=0, h=0, lat=0):
 	'''
+		Calcalation of the refractive index following Hohenkerk & Sinclair
+		(1985), which is the same one as used in SLALIB. 
+
+		This equation is based on the Barrell & Sears equation, but with
+		slightly different values. The refractive index values here seem 
+		to be a bit higher than the original Barrel and Sears values.
+
+		Parameters
+		----------
+		l 	:	float
+			Wavelength in microns
+		T 	:	float
+			Temperature in Kelvin
+		p 	:	float
+			Atmospheric pressure in Pa
+		RH 	:	float
+			fractional humidity (0=<H=<1)
+		h 	:	float
+			Altitude of the observer (default=0 m)
+		lat :	float
+			Latitude of the observer (default=-24)
+		
+		Returns
+		-------
+		n 	:	float
+			Refractive index of atmospheric air at the given conditions
+	'''
+
+	# Some constants
+	alpha = 0.0065
+	delta = 18.36
+	R = 8314.36
+	Md = 28.966
+	Mw = 18.016
+
+	# Pressure in millibars
+	P0 = p/100 
+
+	# Local acceleration due to gravity
+	g_local = 9.784 * (1 - 0.0026 * np.cos(2*np.deg2rad(lat))) - 0.00000028*h
+	gamma = g_local * Md / (R * alpha)
+
+	# Partial pressure of water vapour at observer
+	Pw = RH * (T/247.1)**18.36
+	
+	# Refractivity * 1e6 at reference conditions
+	A = (287.604 + 1.6288*l**-2 + 0.0136*l**-4) * 273.15/1013.25
+	# A = (287.5522 + 1.6288*l**-2 + 0.0136*l**-4) * 273.15/1013.25
+	
+	P = (P0 + (1 - Mw/Md) * gamma/(delta - gamma) * Pw )\
+		- (1 - Mw/Md) * gamma / (delta-gamma) * Pw
+	
+	# Refractive index at given conditions
+	n = 1 + 1e-6 * (A * P - 11.2684*Pw)/T
+	return n
+
+
+def slalib(l, T=288.15, p=101325, RH=0, h=0, lat=-24, tlr=0.0065, eps=1e-10):
+	'''
+		NOTE: This function still works for now, but is considered DEPRECATED.
+			  Use BarrellAndSears() or HohenkerkAndSinclair() instead.
+
 		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.
@@ -228,8 +402,8 @@ def slalib(l, T=288.15, p=101325, H=0, h=0, lat=-24, tlr=0.0065, eps=1e-10):
 			Temperature in Kelvin
 		p 	:	float
 			Atmospheric pressure in Pa
-		H 	:	float
+		RH 	:	float
-			fractional humidity (0=<H=<1)
+			Relative humidity (0=<H=<1)
 		h 	:	float
 			Altitude of the observer (default=0 m)
 		lat :	float
@@ -243,6 +417,14 @@ def slalib(l, T=288.15, p=101325, H=0, h=0, lat=-24, tlr=0.0065, eps=1e-10):
 		-------
 		n 	:	float
 			Refractive index of atmospheric air at the given conditions
+	
+		Notes
+		-----
+		1. 	I don't understand why this function works. It calculates the 
+			refraction and somehow returns the refractive index. Perhaps 
+			this is a very particular set of input parameters that make it 
+			work. It does closely follow BarrellAndSears() for standard
+			conditions.
 	'''
 	try:
 		from pyslalib import slalib
@@ -256,14 +438,14 @@ def slalib(l, T=288.15, p=101325, H=0, h=0, lat=-24, tlr=0.0065, eps=1e-10):
 	if type(l) == float:
 		l = [l]
 	for wl in l:
-		R.append(slalib.sla_refro(zObs, h, T, p, H, wl, lat, tlr, eps))
+		R.append(slalib.sla_refro(zObs, h, T, p, RH, 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):
+def Mathar(l, T=288.15, p=101325, RH=0):
 	'''
 		Refractivity model by Mathar (2007), valid above 1.4 microns
 		Code taken from
@@ -291,7 +473,7 @@ def Mathar(l, T=288.15, p=101325, H=0):
 	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:
+	elif 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
@@ -311,7 +493,7 @@ def Mathar(l, T=288.15, p=101325, H=0):
 		Href = 10          #%, reference humdity in percent
 		########################################################################
 	
-	if l > 2.65 and l <= 4.275:
+	elif 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
@@ -331,7 +513,7 @@ def Mathar(l, T=288.15, p=101325, H=0):
 		Href = 10          #%, reference humdity in percent
 		########################################################################
 
-	if l > 4.275 and l <= 6.35:
+	elif 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
@@ -350,7 +532,7 @@ def Mathar(l, T=288.15, p=101325, H=0):
 		Href = 10          #%, reference humdity in percent
 		########################################################################
 
-	if l > 6.35 and l <= 15.05:
+	elif 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
@@ -370,7 +552,7 @@ def Mathar(l, T=288.15, p=101325, H=0):
 		########################################################################
 
 
-	if l > 15.05:
+	elif 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
@@ -387,12 +569,12 @@ def Mathar(l, T=288.15, p=101325, H=0):
 		σ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
+		Href = 10          # %, reference humdity in percent
 		########################################################################
 
 	s = 1e4/l # cm^-1
 	n = 1
-	H = H*100
+	H = RH*100
 
 	for j in range(0, 6):
 		n += ( cref[j] + cT[j]*(1/T-1/Tref) + cTT[j]*(1/T-1/Tref)**2
@@ -404,6 +586,7 @@ def Mathar(l, T=288.15, p=101325, H=0):
 	return n 
 
 
+# Other necessary conversion or calculations
 def K2C(K):
 	'''
 		Kelvin to Celsius conversion
@@ -420,7 +603,8 @@ def K2C(K):
 	'''
 	return K - 273.15
 
-def hPA2mmHg(hPa):
+
+def hPa2mmHg(hPa):
 	'''
 		hPa to mmHg conversion
 
@@ -436,3 +620,41 @@ def hPA2mmHg(hPa):
 	'''
 	return 0.7500616827 * hPa
 
+
+def Pa2Torr(Pa):
+	'''
+		Pa to Torr conversion
+
+		Parameters
+		----------
+		Pa 	:	float
+			Atmospheric pressure in Pa (Pascal)
+
+		Returns
+		-------
+		Torr 	: 	float
+			Atmospheric pressure in Torr
+	'''
+
+	return Pa * 760/101325
+
+
+def svp(T):
+	'''
+		Saturation vapour pressure in Pa.
+
+		Parameters
+		----------
+		T 	: float
+			Atmospheric temperature in Kelvin
+
+		Returns
+		-------
+		svp	: float
+			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)

MANIFEST.in

@@ -0,0 +1 @@
+include AstroAtmosphere/data/1976USSA.txt

README.md

@@ -2,7 +2,7 @@
 
 **Note:**
 It would be appreciated if a reference to the following work, for which this package was originally build, is included whenever this code is used for a publication: 
-*Quantification of the expected residual dispersion of the MICADO Near-IR imaging instrument*, van den Born & Jellema, 2020, MNRAS. Accepted 2020 June 23. Received 2020 June 23; in original form 2020 January 24.
+*Quantification of the expected residual dispersion of the MICADO Near-IR imaging instrument*, van den Born & Jellema, 2020, MNRAS, DOI: [10.1093/mnras/staa1870](https://doi.org/10.1093/mnras/staa1870).
 
 
 **Author info**: *Joost van den Born*, *born@astron.nl*
@@ -41,7 +41,7 @@ disp = dispersion(lat, h)
 disp.setReducedHeight(P, rho)
 
 # Calculation of the atmopheric dipsersion
-atm_dispersion = ref.cassini(n1, z0) * 3.6e6
+atm_dispersion = disp.cassini(n1, n2, z0) * 3.6e6
 print ('The dispersion is %.03f milli arc seconds' %(atm_dispersion))
 ```
 
@@ -60,11 +60,14 @@ The package contains the following files:
 - `/data/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.
+- `misc.py` 				    : Contains several miscellaneous functions, e.g. to quickly calculate the refractive index.
+- `neoslalib.py`                : Contains a Python native port of the `pyslalib.SLA_REFRO()` routine, translated from the original FORTRAN code.
+- `observatories.py` 			: Contains a dictionary with the average atmospheric conditions for various large observatories around the world.
 - `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).
 
 
+
 ## The Ciddor refractivity model
 
 The refractivity model from Ciddor, 1996, calculates the index of refraction of air using the following parameters:
@@ -85,7 +88,7 @@ Basic functionality of other refractivity models is included as well and is disc
 
 ## The atmospheric refraction and atmospheric dispersion models
 
-If the atmospheric refraction is denoted as R(n), then atmospheric dispersion is defined as R(n1) - R(n2). Thus they are very similar in function and also in use.
+If the atmospheric refraction is denoted as R(n), then atmospheric dispersion is defined as R(n1) - R(n2). Thus they are very similar in function and use.
 
 
 The available models in this package:
@@ -95,6 +98,9 @@ The available models in this package:
 - `oriani()` - Oriani's theorem ( Atan(z) + Btan^3(z) )
 - `tan5()` - Oriani's theorem expanded ( Atan(z) + Btan^3(z) + Ctan^5(z)  )
 - `corbard()` - Error function refraction model
+- `matharExponential()` - Barometric exponential model by R. Mathar.
+
+The package also includes a python native port of the `SLA_REFRO()` FORTRAN routine, ported from pyslalib. See below for more information.
 
 
 Input parameters for the spherical atmosphere models (all but the first):
@@ -104,7 +110,7 @@ Input parameters for the spherical atmosphere models (all but the first):
 - xc - CO2 density              [ppm]
 - lat - Latitude                [deg]
 - h - Height above sea level    [m]
-- l - Wavelength(s)             [um]
+- l - Wavelength(s)             [μm]
 - z - Zenith angle              [deg]
 
 An example - To calculate the atmospheric refraction we set up a `refraction()` object. This requires the latitude and altitude of the observer. The atmospheric dispersion is calculated similarly, but by exchanging the `refraction()` object by a `dispersion()` object and parsing two refractive indices into the `cassini()` function.
@@ -115,7 +121,7 @@ rho = at.rho(p=101325, T=273.15, RH=0.2, xc=300)
 
 # Initializing refraction model and setting the reduced height
 ref = refraction(lat=30, h=2200)
-ref.setReducedHeight(p=101325, rho)
+ref.setReducedHeight(p=101325, rho=rho)
 
 # Calculating the atmospheric refraction
 atm_ref = ref.cassini(n, zenith=30)
@@ -183,21 +189,22 @@ Included are
 | ---------- | -------- | ------- |
 | Barrell & Sears (1939) | `BarrellAndSears()` | Outdated model, especially poor for short wavelengths. |
 | Ciddor (1996) | `Ciddor()` | Only at standard atmosphere. CO2 density can be changed. |
-| Edlen (1953) | `Edlen1953()` | Only at standard atmosphere. Function requires wavelength only. |
+| Edlen (1953) | `Edlen1953()` | The equation by Edlèn, 1953, combined with modification from Barrel (1951), to include temperature, pressure and humidity dependence. |
-| Edlen (1966) | `Edlen1966()` | Only at standard atmosphere. Function requires wavelength only. |
+| Edlen (1966) | `Edlen1966()` | Allows for any temperature, pressure, humidity and CO2 density.|
-| Owens (1967) | `Owens()` | - |
+| Owens (1967) | `Owens()` | Allows for any temperature, pressure, humidity. |
-| Bonsch & Potulski (1998) | `BonschPotulski()` | Function does not include humidity |
+| Bonsch & Potulski (1998) | `BonschPotulski()` | Allows for any temperature, pressure, humidity and CO2 density. |
 | Birch & Downs (1994) | `BirchDowns()` | Only at standard atmosphere. Function requires wavelength only. |
 | Peck & Reeder (1972) | `PeckReeder()` | Only at standard atmosphere. Function requires wavelength only. |
-| Hohenkerk & Sinclair (1985) | `slalib()` | Only at standard atmosphere. Function requires wavelength only. |
+| Mathar (2006) | `Mathar()` | Allows for temperature, pressure and humidity specification. Works best at wavelengths above 1.4 μm. |
+| Hohenkerk & Sinclair (1985) | `HohenkerkAndSinclair()` | Based on Barrell & Sears equation. Generally, gives slightly larger values.|
+| Hohenkerk & Sinclair (1985) | `slalib()` | Uses pySLALIB to retrieve the refractive index at standard atmospheric conditions.|
 
 The last function extracts the refractive index from the `pyslalib.sla_refro()` function and thus requires [pyslalib](https://pypi.org/project/pySLALIB/). This is an optional dependency for the setup of AstroAtmosphere.
 
 ## Quick functions
 
-Some functions that assume default conditions for the atmosphere. Then only the wavelength and zenith angle are required to calculate the refraction and dispersion. 
+Some functions can make use of the `observatories` dictionary as input. These include default average conditions of the atmosphere at these locations. When these are used, only the wavelength and zenith angle are required to calculate the refraction and dispersion. 
 
-Presently, only the `'STANDARD'` and `'CERRO_ARMAZONES'` flags are known. Contact the author to add additional standard conditions.
 
 ```python
  # Quick refractive index
@@ -209,3 +216,21 @@ rq = quick_refraction(l1, zenith_obs, conditions='STANDARD')
 # Quick dispersion
 dq = quick_dispersion(l1, l2, zenith_obs, conditions='STANDARD')
 ```
+
+Presently, several observatories are known. These can all be found in the `observatories` dictionary. Contact the author to add additional standard conditions.
+
+| Observatory                                   | Flag              | Notes                                         | 
+| --------------------------------------------- | ----------------- | --------------------------------------------- |
+| Standard conditions                           | `STANDARD`        |                                               |
+| Cerro Armazones                               | `CERRO_ARMAZONES` | Location of the future ELT                    |
+| Cerro Paranal                                 | `CERRO_PARANAL`   | Location of the VLT, VISTA, VST, NGTS and SSO |
+| Roque de los Muchachos Observatory (La Palma) | `LA_PALMA`        | GTC, INT, WHT, NOT, and others.               |
+| La Silla Observatory                          | `LA_SILLA`        | ESO 3.6m telescope, NTG, BlackGEM and others. | 
+| Las Campanas Observatory                      | `LAS_CAMPANAS`    | Site of the future GMT.                       |
+| Mauna Kea Observatory (Hawaii)                | `MAUNA_KEA`       | Location of Keck, Subaru, Gemini North and others.|
+
+## Other
+
+This package also contains a python native port of the `SLA_REFRO()` routine of the pyslalib package, which was originally written in FORTRAN by P.T. Wallace. It can be called directly using `sla_refro()`. When using the same units as the other refraction functions, you might prefer using `slalib_refraction()` or `slalib_dispersion()`. Do note that this implementation is slower than the (compiled) FORTRAN implementation.
+
+Finally, the dispersion model used by Alexei Filippenko in 1982 is also included. This model is a combination of the (modified) Edlèn equation and a plane-parallel atmospheric geometry. This calculation can be accessed by calling the `Filippenko1982()` function.

scripts/ADC.py

@@ -10,8 +10,10 @@ class ADC:
 		self.p3 = prism('S-FPL51', 15, wedge[2], wedge[1])
 		self.p4 = prism('S-LAH71', 15, wedge[1], wedge[0])
 		self.theta 	= 90
+		self.FR 	= 0
 
 	def setOrientation(self, theta):
+		theta = float(theta)
 		self.p1.setPrismOrientation(theta)
 		self.p2.setPrismOrientation(theta)
 		self.p3.setPrismOrientation(180-1*theta)
@@ -19,10 +21,12 @@ class ADC:
 		self.theta = theta
 
 	def setFieldRotation(self, FR):
+		FR = float(FR)
 		self.p1.setPrismOrientation(self.theta-FR)
 		self.p2.setPrismOrientation(self.theta-FR)
 		self.p3.setPrismOrientation(180-1*self.theta-FR)
 		self.p4.setPrismOrientation(180-1*self.theta-FR)
+		self.FR = FR
 
 	def calculateDispersion(self, l1, l2, hx=0, hy=0, refraction=False):
 		'''
@@ -96,6 +100,7 @@ class prism:
 		self.tilt1y 	= 0
 		self.tilt2x		= tilt2
 		self.tilt2y 	= 0
+		self.theta 		= 90
 
 		# Remember original tilts for reset:
 		self.tilt1x_OG	= tilt1
@@ -116,7 +121,7 @@ class prism:
 	def setPrismOrientation(self, theta):
 		self.tilt1x, self.tilt1y = tiltRotation(self.tilt1x_OG, self.tilt1y_OG, theta)
 		self.tilt2x, self.tilt2y = tiltRotation(self.tilt2x_OG, self.tilt2y_OG, theta)
-	
+		self.theta = theta
 
 
 def tiltRotation(tx,ty, theta):

scripts/compare_refractivity_models.py

@@ -13,7 +13,8 @@ import AstroAtmosphere
 @ Description:	A comparison of various atmospheric refraction models.
 				Also these are used to calculate the difference in 
 				atmospheric dispersion, given a set of conditions.
-				This is a reproduction of Figure 1 of Spàno (2014).
+				This is a reproduction of and expansion on Figure 1 of 
+				Spàno (2014).
 '''
 
 
@@ -54,7 +55,7 @@ class atmosphericModel:
 	def Owens(self, l=None):
 		if not l:
 			l = self.lrange
-		return (AstroAtmosphere.Owens(l, T=self.T, p=self.P, H=self.H)-1)
+		return (AstroAtmosphere.Owens(l, T=self.T, p=self.P, RH=self.H)-1)
 
 
 	def BonschPotulski(self, l=None):
@@ -72,10 +73,34 @@ class atmosphericModel:
 			l = self.lrange
 		return (AstroAtmosphere.PeckReeder(l)-1)
 
+	def HohenkerkAndSinclair(self, l=None):
+		if not l:
+			l = self.lrange
+		return (AstroAtmosphere.HohenkerkAndSinclair(l, T=self.T, p=self.P, RH=self.H)-1)
+
 	def slalib(self, l=None):
 		if not l:
 			l = self.lrange
-		return(AstroAtmosphere.slalib(l, T=self.T, p=self.P, H=self.H)-1)
+		try:
+			from pyslalib import slalib
+			return(AstroAtmosphere.slalib(l, T=self.T, p=self.P, RH=self.H)-1)
+		except ImportError as e:
+			print('Could not import pyslalib.slalib. Using slower Python native function sla_refro().')
+			h=0
+			lat=-24
+			tlr=0.0065
+			eps=1e-10
+			zObs = np.deg2rad(45.07182)
+			p = self.P/100
+			R = []
+			if type(l) == float:
+				l = [l]
+			for wl in l:
+				R.append(AstroAtmosphere.sla_refro(zObs, h, self.T, p, self.H, wl, lat, tlr, eps))
+			if len(l) == 1:
+				return R[0]
+			else:
+				return np.asarray(R)
 
 
 if __name__ == '__main__':
@@ -92,6 +117,7 @@ if __name__ == '__main__':
 	BonschPotulskiData 	= 1e8 * at.BonschPotulski()
 	BirchDownsData 		= 1e8 * at.BirchDowns()
 	PeckReederData 		= 1e8 * at.PeckReeder()
+	HSdata 				= 1e8 * at.HohenkerkAndSinclair()
 	slalibData 			= 1e8 * at.slalib()
 	
 	# Plotting the results
@@ -106,13 +132,14 @@ if __name__ == '__main__':
 	ax2.plot(at.lrange, BonschPotulskiData - CiddorData, 'g-.', label='Bonsch & Potulski (1998)')
 	ax2.plot(at.lrange, BirchDownsData - CiddorData, 'r-.', label='Birch & Downs (1993, 1994)')
 	ax2.plot(at.lrange, PeckReederData - CiddorData, 'b--', label='Peck & Reeder (1972)')
+	ax2.plot(at.lrange, HSdata - CiddorData, 'y--', label='Hohenkerk & Sinclair (1985)')
 	ax2.plot(at.lrange, slalibData - CiddorData, 'y-', label='ZEMAX / SLALIB')
 	box = ax2.get_position()
 	ax2.set_position([box.x0, box.y0, box.width, box.height])
 	ax2.legend(loc='center left', bbox_to_anchor=(1, 0.5))
 	ax2.set_xlabel('Wavelength ($\mu m$)')
 	ax2.set_ylabel('$\Delta$n $\\times$ $10^{8}$')
-	ax2.set_ylim([-3,3])
+	ax2.set_ylim([-3,7])
 	plt.tight_layout()
 	plt.show()
 
@@ -128,6 +155,7 @@ if __name__ == '__main__':
 						at.BonschPotulski,
 						at.BirchDowns,
 						at.PeckReeder,
+						at.HohenkerkAndSinclair,
 						at.slalib
 					]
 	
@@ -166,11 +194,14 @@ if __name__ == '__main__':
 						'Ciddor (1996)',
 						'Edlen (1953)',
 						'Edlen (1966)',
+						'Owens (1967)',
 						'Bonsch & Potulski (1998)',
 						'Birch & Downs (1993, 1994)',
 						'Peck & Reeder (1972)',
+						'Hohenkerk & Sinclair (1985)',
 						'slalib']
 	tableData = []
+
 	for i in range(len(methodStringList)):
 		tableData.append([methodStringList[i], filterDispersion[0][i], filterDispersion[1][i], filterDispersion[2][i], filterDispersion[3][i]])
 	print(tabulate(tableData, headers=['Model', 'I-band\n(0.8 - 0.934)', 'J-band\n(1.15 - 1.345)', 'H-band\n(1.49 - 1.78)', 'K-band\n(1.97 - 2.32)']))

setup.py

@@ -21,14 +21,14 @@ def setup_package():
     # write_version_py()
 
     setup(name='AstroAtmosphere',
-          version='1.0',
+          version='1.6',
           description='Collection of equations for refractive index, refraction and dispersion calculations of atmospheric air.',
           long_description=__readme__,
           long_description_content_type='text/markdown',
           author='Joost van den Born',
           author_email='born@astron.nl',
           url='https://gitlab.astro-wise.org/micado/atmosphericmodel',
-          download_url='https://gitlab.astro-wise.org/micado/atmosphericmodel/-/archive/v_10/atmosphericmodel-v_10.tar.gz',
+          download_url='https://gitlab.astro-wise.org/micado/atmosphericmodel/-/archive/v1.6/atmosphericmodel-v1.6.tar.gz',
           license='GNU General Public License v3',
           data_files = [('',['LICENSE.txt'])],
           package_dir={'AstroAtmosphere': 'AstroAtmosphere'},