updated README

AstroAtmosphere/__init__.py

@@ -5,6 +5,6 @@ from .ciddorModel import *
 from .dispersionModels import *
 from .refractionModels import *
 from .refractionModels import *
-from .examples import *
+from .misc import *
 from .refractivityModels import *
-
+from .neoslalib 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
-			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.
+			Returns the atmospheric dispersion assuming Cassini's homegeneous
+			atmosphere model.
 		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 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,
+	def cassiniError(self, zenith, 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,457 @@ 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 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,364 @@
+#!/usr/bin/env Python3
+from .ciddorModel import Observatory
+from .refractionModels import refraction
+from .dispersionModels import dispersion
+from .neoslalib import sla_refro
+
+'''
+@ 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, 
+							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
+		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
+		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 == 'STANDARD':
+		# Parameters at Cerro Armazones
+		T   = 288.15    # k
+		p   = 101325    # Pa
+		RH  = 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
+		RH  = 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=RH, xc=xc)
+	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
+		zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+		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
+		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 == 'STANDARD':
+		# Parameters at Cerro Armazones
+		T   = 288.15    # k
+		p   = 101325    # Pa
+		RH  = 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
+		RH  = 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=RH, xc=xc)
+
+	# 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.
+		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
+			'STANDARD' or 'CERRO_ARMAZONES', refers
+			to the standard atmospheric conditions.
+		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 == 'STANDARD':
+		# Parameters at Cerro Armazones
+		T   = 288.15    # k
+		p   = 101325    # Pa
+		RH  = 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
+		RH  = 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=RH, xc=xc)
+	n2 	= at.n_tph(l=l2, T=T, p=p, RH=RH, xc=xc)
+
+	# 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
+			'STANDARD' or 'CERRO_ARMAZONES', refers
+			to the standard atmospheric conditions.
+		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-10)
+		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 == 'STANDARD':
+		# Parameters at Cerro Armazones
+		T   = 288.15    # k
+		p   = 101325    # Pa
+		RH  = 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
+		RH  = 0.22
+		xc  = 450       # ppm
+		lat = -24.5983  # degrees
+		h   = 3064      # m
+
+	# 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 type(l) == float or type(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
+			'STANDARD' or 'CERRO_ARMAZONES', refers
+			to the standard atmospheric conditions.
+		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-10)
+		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/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
-				error propagation if needed.
+				index of refraction for a wide range of conditions, with corresponding error propagation if needed.
 				The actual refraction is calculated by any of the following models:
 				- Plane parallel atmosphere
 				- Cassini spherical homogeneous atmosphere
+				- 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
-			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.
+			Returns the atmospheric refraction assuming Cassini's homegeneous
+			atmosphere model.
 		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,72 +116,46 @@ class refraction(Observatory):
 		return np.degrees(_R1)
 
 
-
-	def tan5(self, n1, zenith, reduced_height=None):
+	def cassiniError(self, n, dn, zenith, dz=0, reduced_height=None):
 		'''
-			Same as oriani's theorem, but including a higher order term.
-			Found in Corbard et al., 2019.
+		Calculates the uncertainty in the refraction angle in degrees.
+		See Observatory().n_tph() and Observatory().dn_tph() on how to
+		calculate the refractive index and it's uncertainty.
 
-			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.
+		Parameters
+		----------
+		n 	: 	float
+			Refractive index of air
+		dn 	: 	float
+			Uncertainty in the refractive index of air
+		zenith : float
+			Angle of observation in degrees.
+			zenith = 0 means the observation is directly overhead.
+		dz  : 	float (optional)
+			Uncertainty in the zenith angle, in degrees.
+		reduced_height : float (optional)
+			Reduced height of the atmosphere. Default value is calculated
+			from the object attributes.
 
-			Returns
-			-------
-			R : float
-				Atmospheric refraction in degrees
+		Returns
+		-------
+		dR 	: float
+			Uncertainty in the atmospheric refraction in degrees.
 		'''
 		if reduced_height:
-			H = reduced_height
+			_H = reduced_height
 		else:
-			H = self.H
+			_H = self.H
+		_zenith = np.deg2rad(zenith)
+		_dz 	= np.deg2rad(dz)
+		_R 		= self.rc
+		_dRdn 	= (_R * np.sin(_zenith)) / ((_R + _H) * np.sqrt(1 - ((n**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2)))
 
-		_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)
+		_dRdz 	= (_R * np.cos(_zenith)) / (_R + _H) * ( n / np.sqrt(1-(n**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2) - 1 / np.sqrt(1-(_R**2 * np.sin(_zenith)**2)/(_R + _H)**2) )
 
+		_dR 	= np.sqrt( (_dRdn * dn)**2 + (_dRdz * _dz)**2)
 
-
-	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)
-
+		return np.degrees(_dR)
 
 
 	def corbard(self, n1, zenith, reduced_height=None):
@@ -246,8 +191,231 @@ class refraction(Observatory):
 		return np.degrees(_R)
 
 
+	def gammaCoefficients(self, n1, reduced_height=None):
+		'''
+			Submethod for matharExponential(). Calculates the gamma
+			coefficients.
 
-	def refractionIntegral(self, n1, z0, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None):
+			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(n * 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
@@ -263,7 +431,7 @@ class refraction(Observatory):
 			----------
 			n1 : float
 				Refractive index at the wavelength of interest
-			z0 : float
+			zenith : float
 				Angle of observation in degrees.
 				zenith = 0 means the observation is directly overhead.
 			R0 : float (optional)
@@ -295,7 +463,7 @@ class refraction(Observatory):
 			_Dens = rhoData
 			_h 	  = heightData
 
-		z0 						= np.deg2rad(z0)
+		zenith 		= np.deg2rad(zenith)
 		_Dens 	= np.asarray(_Dens)
 		_h 	 	= np.asarray(_h) 
 		_R 		= R0 + _h * 1000
@@ -304,7 +472,7 @@ class refraction(Observatory):
 		_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) )
+		_z 		= np.arcsin( (n0 * R0) / (_n_h * _R) * np.sin(zenith) )
 		
 		# Calculation of log space derivative of n and r
 		_lnn = np.log(_n_h)
@@ -320,45 +488,18 @@ class refraction(Observatory):
 		return np.degrees(np.sum(_Integrand))
 
 
-
-	def psi(self, x):
+	def setReducedHeight(self, p, rho):
 		'''
-			Submethod for corbard(). Described explicitly in Corbard et al.,2019.
+			Sets the reduced atmospheric height.
 
 			Parameters
 			----------
-			x 	: 	float
-				Input of the function
-
-			Returns
-			-------
-			y 	: 	float
-				Output of the function
+			p 	: float
+				Atmospheric pressure in Pascal
+			rho : float
+				Density of the atmospheric air at the observer
 		'''
-		_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
-
+		self.H  	= self.set_H(p, rho) 		# Reduced height of the atmosphere assuming ideal gas law
 	
 
 	def set_H(self, p, rho):
@@ -387,71 +528,54 @@ class refraction(Observatory):
 		_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
+	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
 			-------
-			H 	: float
-				Reduced height of the atmosphere in meters
+			rc 	: float
+				Local radius of curvature of the Earth in kilometers.
 		'''
-		_kb = 1.38064852e-23 # kg m2 s-2 K-1
-		_m  = 4.809651698e-26 # kg (weight per molecule of air, assuming dry air)
-
+		_A   = 6378.137 # km
+		_B   = 6356.752 # km
 		_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)
+		_rc0 = (_A*_B)**2 / (_A**2 * np.cos(_phi)**2 + _B**2 * np.sin(_phi)**2)**(3/2)
+		return _rc0 + self.h/1000
 
-	def cassiniError(self, n, dn, zenith, dz=0, reduced_height=None):
+
+	def tan5(self, n1, zenith, reduced_height=None):
 		'''
-		Calculates the uncertainty in the refraction angle in degrees.
-		See Observatory().n_tph() and Observatory().dn_tph() on how to
-		calculate the refractive index and it's uncertainty.
+			Same as oriani's theorem, but including a higher order term.
+			Found in Corbard et al., 2019.
 
-		Parameters
-		----------
-		n 	: 	float
-			Refractive index of air
-		dn 	: 	float
-			Uncertainty in the refractive index of air
-		zenith : float
-			Angle of observation in degrees.
-			zenith = 0 means the observation is directly overhead.
-		dz  : 	float (optional)
-			Uncertainty in the zenith angle, in degrees.
-		reduced_height : float (optional)
-			Reduced height of the atmosphere. Default value is calculated
-			from the object attributes.
+			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
-		-------
-		dR 	: float
-			Uncertainty in the atmospheric refraction in degrees.
+			Returns
+			-------
+			R : float
+				Atmospheric refraction in degrees
 		'''
 		if reduced_height:
-			_H = reduced_height
+			H = reduced_height
 		else:
-			_H = self.H
-		_zenith = np.deg2rad(zenith)
-		_dz = np.deg2rad(dz)
-		_R = self.rc
-		_dRdn = (_R * np.sin(_zenith)) / ((_R + _H) * np.sqrt(1 - ((n**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2)))
+			H = self.H
 		
-		_dRdz = (_R * np.cos(_zenith)) / (_R + _H) * ( n / np.sqrt(1-(n**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2) - 1 / np.sqrt(1-(_R**2 * np.sin(_zenith)**2)/(_R + _H)**2) )
+		_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)
 
-		_dR = np.sqrt( (_dRdn * dn)**2 + (_dRdz * _dz)**2)
-
-		return np.degrees(_dR)

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):
 	'''
 		From Barrell & Sears, 1939
 		Equation 7.7, page 52
@@ -24,7 +24,7 @@ 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)
 
 		Returns
@@ -34,10 +34,11 @@ def BarrellAndSears(l, T=288.15, p=101325, H=0):
 	'''
 	T = K2C(T)
 	P = hPA2mmHg(p/100)
-	R = (0.378125 + 0.0021414 / l**2 + 0.00001793 / l**4) * (P*(1 + (1.049 - 0.0157*T) * P * 1e-6)/(1 + 0.003661*T)) - (0.0624 - 0.000680 / l**2) * H / (1 + 0.003661*T)
+	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) * RH / (1 + 0.003661*T)
 	n = R/1e6 + 1
 	return n
 
+
 def Ciddor(l, CO2=450):
 	'''
 		From Ciddor, 1996
@@ -62,6 +63,7 @@ def Ciddor(l, CO2=450):
 	n = R/1e8 + 1
 	return n
 
+
 def Edlen1953(l):
 	'''
 		From Edlen, 1953
@@ -84,6 +86,7 @@ def Edlen1953(l):
 	n = R/1e8 + 1
 	return n
 
+
 def Edlen1966(l):
 	'''
 		From Edlen, 1966
@@ -106,7 +109,8 @@ def Edlen1966(l):
 	n = R/1e8 + 1
 	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
@@ -121,8 +125,8 @@ def Owens(l, T=288.15, p=101325, H=0):
 			Temperature in Kelvin
 		p 	: 	float
 			Atmospheric pressure in Pa
-		H 	: 	float
-			Humdity (0<=H<=1)
+		RH 	: 	float
+			Relative humdity (0<=H<=1)
 
 		Returns
 		-------
@@ -132,10 +136,11 @@ def Owens(l, T=288.15, p=101325, H=0):
 	P = hPA2mmHg(p/100)
 	T = K2C(T)
 	s = 1 / l
-	R = 8342.13 + 2406030 / (130 - s**2) + 15997 / (38.9 - s**2) * (P/720.775) * ((1 + P * (0.817 - 0.0133) * 1e-6)/(1 + 0.0036610*T)) - H * (5.722 - 0.0457*s**2)
+	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)) - RH * (5.722 - 0.0457*s**2)
 	n = R/1e8 + 1
 	return n
 
+
 def BonschPotulski(l, T=288.15, p=101325, CO2=450):
 	'''
 		From Bonsch & Potulski, 1998
@@ -168,6 +173,7 @@ def BonschPotulski(l, T=288.15, p=101325, CO2=450):
 	n = Rtp + 1
 	return n
 
+
 def BirchDowns(l):
 	'''
 		From Birch & Downs, 1993 and 1994 (correction)
@@ -190,6 +196,7 @@ def BirchDowns(l):
 	n = R/1e8 + 1
 	return n
 
+
 def PeckReeder(l):
 	'''
 		From Peck & Reeder, 1972
@@ -213,8 +220,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 +297,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
-			fractional humidity (0=<H=<1)
+		RH 	:	float
+			Relative humidity (0=<H=<1)
 		h 	:	float
 			Altitude of the observer (default=0 m)
 		lat :	float
@@ -243,6 +312,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 +333,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
@@ -392,7 +469,7 @@ def Mathar(l, T=288.15, p=101325, H=0):
 
 	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
@@ -420,6 +497,7 @@ def K2C(K):
 	'''
 	return K - 273.15
 
+
 def hPA2mmHg(hPa):
 	'''
 		hPa to mmHg conversion

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,7 +60,7 @@ 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.
 - `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).
 
@@ -85,7 +85,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 +95,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):
@@ -189,7 +192,8 @@ Included are
 | Bonsch & Potulski (1998) | `BonschPotulski()` | Function does not include humidity |
 | 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. |
+| 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.
 
@@ -209,3 +213,7 @@ rq = quick_refraction(l1, zenith_obs, conditions='STANDARD')
 # Quick dispersion
 dq = quick_dispersion(l1, l2, zenith_obs, conditions='STANDARD')
 ```
+
+## 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()`.

scripts/ADC.py

@@ -9,9 +9,11 @@ class ADC:
 		self.p2 = prism('S-FPL51', 15, wedge[1], wedge[2])
 		self.p3 = prism('S-FPL51', 15, wedge[2], wedge[1])
 		self.p4 = prism('S-LAH71', 15, wedge[1], wedge[0])
-		self.theta = 90
+		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,15 @@ 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)
+		return(AstroAtmosphere.slalib(l, T=self.T, p=self.P, RH=self.H)-1)
 
 
 if __name__ == '__main__':
@@ -92,6 +98,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 +113,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 +136,7 @@ if __name__ == '__main__':
 						at.BonschPotulski,
 						at.BirchDowns,
 						at.PeckReeder,
+						at.HohenkerkAndSinclair,
 						at.slalib
 					]
 	
@@ -166,11 +175,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.2',
           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.2/atmosphericmodel-v1.2.tar.gz',
           license='GNU General Public License v3',
           data_files = [('',['LICENSE.txt'])],
           package_dir={'AstroAtmosphere': 'AstroAtmosphere'},