Fixed the setup

.gitignore

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

AstroAtmosphere/__init__.py

@@ -8,3 +8,4 @@ from .refractionModels import *
 from .misc import *
 from .refractivityModels import *
 from .neoslalib import *
+from .observatories import *

AstroAtmosphere/dispersionModels.py

@@ -387,6 +387,73 @@ class dispersion(Observatory):
 		return np.degrees(_dispersionError)
 
 
+	def independentCassiniError(self, n1, n2, zenith, dn1=0, dn2=0, dz=0, 
+									lat=None, h=None):
+		'''
+			Calculates the uncertainty in the atmospheric dispersion, 
+			following the Cassini model, without assumptions about the
+			refractivity model.
+			This involves a much shorter use the chain rule. 
+
+			Not taken into account is the errror propagation in the scale height
+			of the atmosphere, from my results this should give a small or even 
+			negligible difference.
+
+			Parameters
+			----------
+			n1 	: 	float
+				Refractive index at the first wavelength of interest at the 
+				observer
+			n2 	: 	float
+				Refractive index at the second wavelength of interest at the 
+				observer
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			dn1 : float
+				Uncertainty in the provided value of the refractive index
+				at the first wavelength
+			dn2 : float
+				Uncertainty in the provided value of the refractive index
+				at the second wavelength
+			dz 	: float (default: 0)
+				Uncertainty in the zenith angle, in degrees.
+			lat : float (default: None)
+				If entered, overrides the latitude attribute of the class.
+			h 	: float (default=None)
+				If entered, overrides the altitude attribute of the class.
+
+			Returns
+			-------
+			dRR	: float
+				Uncertainty in the atmospheric dispersion in degrees.
+		'''
+		if lat == None:
+			lat = self.lat
+		else:
+			self.lat = lat
+		if h == None:
+			h 	= self.h
+		else:
+			self.h = h
+
+		# Some precalculation and variable setting
+		_zenith = np.deg2rad(zenith)
+		_dz   	= np.deg2rad(dz)
+		_R 		= self.rc
+		_H 		= self.H
+
+		# Partial derivatives for derivation to n1, n2, and z
+		_dddn1  = - _R * np.sin(_zenith) / ( (_R + _H) * np.sqrt(1 - (n1 * _R * np.sin(_zenith) / (_R + _H))**2 ) )
+		_dddn2  =   _R * np.sin(_zenith) / ( (_R + _H) * np.sqrt(1 - (n2 * _R * np.sin(_zenith) / (_R + _H))**2 ) )
+		_dddz 	= (_R * np.cos(_zenith)) / (_R + _H) * ( n1 / np.sqrt(1-(n1**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2) - n2 / np.sqrt(1-(n2**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2) )		
+		
+		# Combining the partial derivatives
+		_dispersionError = np.sqrt( (_dddn1 * dn1)**2 + (_dddn2 * dn2)**2 + (_dddz * dz)**2)
+
+		return np.degrees(_dispersionError)
+
+
 	def corbard(self, n1, n2, zenith, reduced_height=None):
 		'''
 			Corbard et al., 2019, mention an additional formula based on

AstroAtmosphere/misc.py

@@ -3,12 +3,94 @@ from .ciddorModel import Observatory
 from .refractionModels import refraction
 from .dispersionModels import dispersion
 from .neoslalib import sla_refro
+from .observatories import observatories
+from .refractivityModels import Edlen1953, Mathar
+import numpy as np
+import warnings
 
 '''
 @ Author: 		Joost van den Born
 @ Contact: 		born@astron.nl
 @ Description: 	Some functions for quick testing.
 '''
+def get_conditions(location):
+	'''
+		Returns the mean conditions at various site. Currently
+		supported locations:
+		- 'STANDARD'
+		- 'CERRO_ARMAZONES'
+		- 'CERRO_PARANAL'
+		- 'LA_PALMA'
+		- 'LA_SILLA'
+		- 'LAS_CAMPANAS'
+		- 'MAUNA_KEA'
+
+		Parameters
+		----------
+		location 	: string
+			String denoting the desired location
+
+		Returns
+		-------
+		observatory : dict
+			Dictionary containing the atmospheric parameters
+	'''
+	observatory = observatories[location]
+	return observatory
+
+
+def Filippenko1982(l1, l2, zenith, conditions='STANDARD', T=288.15, p=101325,
+							RH=0.0, xc=450, lat=0, h=0, f=None):
+	'''
+		Calculates the refraction following the 1982 paper by 
+		Alexei V. Filippenko.
+
+		The refraction is based on a flat atmosphere and the refractivity
+		values are calculated from a modified Edlen (1953) equation.
+
+		Parameters
+		----------
+		l1 	: 	float
+			Shortest wavelength of interest in micron.
+		l2 	: 	float
+			Longest wavelength of interest in micron.
+		zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+		conditions 	: string
+			Refers to the standard conditions of the 
+			observatory. See get_conditions() for more info.
+			If set to 'None', the conditions specified 
+			after will be used.
+		T 	: float (optional)
+			Temperature in Kelvin
+		p 	: float (optional)
+			Atmospheric pressure in Pa
+		RH 	: float (optional)
+			Relative humidity
+		xc 	: float (optional)
+			CO2 density in parts per million
+		lat : float (optional)
+			Latitude of the observer in degrees
+		h 	: float (optional)
+			Altitude of the observer in meters
+		f 	: float
+			Water vapour pressure (default=None). If None, it is
+			calculated automatically.
+
+		Returns
+		-------
+		R : float
+			Atmospheric refraction in degrees
+	'''
+
+	if conditions is not None:
+		T, p, RH, xc, lat, h, _ = get_conditions(conditions).values()
+
+	n1 = Edlen1953(l1, T=T, p=p, RH=RH, f=f)
+	n2 = Edlen1953(l2, T=T, p=p, RH=RH, f=f)
+	dR  = (n1 - n2) * np.tan(np.deg2rad(zenith))
+	return np.degrees(dR)
 
 
 def quick_refractive_index(l, conditions='STANDARD', T=288.15, p=101325, 
@@ -21,10 +103,14 @@ def quick_refractive_index(l, conditions='STANDARD', T=288.15, p=101325,
 		Parameters
 		----------
 		l 	: float
-			Wavelength in microns
+			Wavelength in microns. If the wavelength is above
+			2.5 microns, then the Mathar refractivity model
+			is used.
 		conditions: string
-			'STANDARD' or 'CERRO_ARMAZONES', refers
+			Refers to the standard conditions of the 
-			to the standard atmospheric conditions.
+			observatory. See get_conditions() for more info.
+			If set to 'None', the conditions specified 
+			after will be used.
 		T 	: float (optional)
 			Temperature in Kelvin
 		p 	: float (optional)
@@ -43,29 +129,18 @@ def quick_refractive_index(l, conditions='STANDARD', T=288.15, p=101325,
 		n 	: 	float
 			Refractive index for the given conditions.
 	'''
-	if conditions == 'STANDARD':
+	if conditions is not None:
-		# Parameters at Cerro Armazones
+		T, p, RH, xc, lat, h, _ = get_conditions(conditions).values()
-		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
 
+	if l <= 2.5:
 		# Initializing dispersion model
 		at  = Observatory()
 
 		# Calculating indices of refraction for l
 		n 	= at.n_tph(l=l, T=T, p=p, RH=RH, xc=xc)
+	else:
+		# Using the Mathar description for infrared wavelengths
+		n   = Mathar(l=l, T=T, p=p, RH=RH)
 	return n
 
 
@@ -79,13 +154,16 @@ def quick_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325,
 		Parameters
 		----------
 		l  	: 	float 
-			Wavelength in microns
+			Wavelength in microns. If the wavelength is above 2.5 micron
+			then the Mathar refractivity model is used.
 		zenith : float
 				Angle of observation in degrees.
 				zenith = 0 means the observation is directly overhead.
 		conditions 	: string
-			'STANDARD' or 'CERRO_ARMAZONES', refers
+			Refers to the standard conditions of the 
-			to the standard atmospheric conditions.
+			observatory. See get_conditions() for more info.
+			If set to 'None', the conditions specified 
+			after will be used.
 		T  	: 	float (optional)
 			Temperature in Kelvin
 		p  	: 	float (optional)
@@ -104,29 +182,18 @@ def quick_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325,
 		R 	: 	float
 			Refraction in degrees
 	'''
-	if conditions == 'STANDARD':
+	if conditions is not None:
-		# Parameters at Cerro Armazones
+		T, p, RH, xc, lat, h, _ = get_conditions(conditions).values()
-		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()
 
+	if l <= 2.5:
 		# Calculating indices of refraction for l
 		n 	= at.n_tph(l=l, T=T, p=p, RH=RH, xc=xc)
+	else:
+		# Using the Mathar description for infrared wavelengths
+		n   = Mathar(l=l, T=T, p=p, RH=RH)
 
 	# Density of the atmosphere (following CIPM-81/91 equations)
 	rho = at.rho(p=p, T=T, RH=RH, xc=xc)
@@ -147,15 +214,19 @@ def quick_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15,
 		Parameters
 		----------
 		l1 	: 	float
-			Shortest wavelength of interest in micron.
+			Shortest wavelength of interest in micron. If the wavelength
+			is above 2.5 microns, then the Mathar model is used.
 		l2 	: 	float
-			Longest wavelength of interest in micron.
+			Longest wavelength of interest in micron. If the wavelength
+			is above 2.5 microns, then the Mathar model is used.
 		zenith : float
 				Angle of observation in degrees.
 				zenith = 0 means the observation is directly overhead.
 		conditions 	: string
-			'STANDARD' or 'CERRO_ARMAZONES', refers
+			Refers to the standard conditions of the 
-			to the standard atmospheric conditions.
+			observatory. See get_conditions() for more info.
+			If set to 'None', the conditions specified 
+			after will be used.
 		T 	: 	float (optional)
 			Temperature in Kelvin
 		p 	: 	float (optional)
@@ -174,30 +245,27 @@ def quick_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15,
 		dR 	: float
 		 	Atmospheric dispersion in degrees
 	'''
-	if conditions == 'STANDARD':
+	if conditions is not None:
-		# Parameters at Cerro Armazones
+		T, p, RH, xc, lat, h, _ = get_conditions(conditions).values()
-		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()
 
+	if l1 <= 2.5 and l2 <= 2.5:
 		# Calculating indices of refraction for l1 and l2
 		n1 	= at.n_tph(l=l1, T=T, p=p, RH=RH, xc=xc)
 		n2 	= at.n_tph(l=l2, T=T, p=p, RH=RH, xc=xc)
+	elif l1 > 1.4 and l2 > 1.4:
+		# Using the Mathar description for infrared wavelengths
+		n1   = Mathar(l=l1, T=T, p=p, RH=RH)
+		n2   = Mathar(l=l2, T=T, p=p, RH=RH)
+	else:
+		warnings.warn('Using two refractivity models (Ciddor and Mathar). \
+Make sure you trust the results!')
+		l1n  = np.min([l1, l2])
+		l2n  = np.max([l1, l2])
+		n1 	 = at.n_tph(l=l1n, T=T, p=p, RH=RH, xc=xc)
+		n2   = Mathar(l=l2n, T=T, p=p, RH=RH)
 
 	# Density of the atmosphere (following CIPM-81/91 equations)
 	rho = at.rho(p=p, T=T, RH=RH, xc=xc)
@@ -228,8 +296,10 @@ def slalib_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325,
 				Angle of observation in degrees.
 				zenith = 0 means the observation is directly overhead.
 		conditions 	: string
-			'STANDARD' or 'CERRO_ARMAZONES', refers
+			Refers to the standard conditions of the 
-			to the standard atmospheric conditions.
+			observatory. See get_conditions() for more info.
+			If set to 'None', the conditions specified 
+			after will be used.
 		T 	: 	float (optional)
 			Temperature in Kelvin
 		p 	: 	float (optional)
@@ -245,11 +315,11 @@ def slalib_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325,
 		tlr :	float
 			Temperature lapse rate (default=0.0065)
 		eps :	float
-			Precision at which the calculation is aborted (default=1e-10)
+			Precision at which the calculation is aborted (default=1e-8)
 		use_pyslalib : bool
 			Flag that indicates that the user wants to use the pySlalib
 			package. If False, the internal python port of the original
-			FORTRAN code will be used. (Default=False).
+			FORTRAN code will be used. (default=False).
 			
 		Returns
 		-------
@@ -257,23 +327,8 @@ def slalib_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325,
 			Refraction in degrees
 	'''
 
-	if conditions == 'STANDARD':
+	if conditions is not None:
-		# Parameters at Cerro Armazones
+		T, p, RH, xc, lat, h, _ = get_conditions(conditions).values()
-		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
@@ -281,7 +336,7 @@ def slalib_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325,
 	_phi  	= np.deg2rad(lat) 		# From degrees to radians
 
 	R = []
-	if type(l) == float or type(l) == int:
+	if isinstance(l, float) or isinstance(l, int): 
 		l = [l]
 	if use_pyslalib:
 		try:
@@ -290,13 +345,13 @@ def slalib_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325,
 		    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)
+			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)
+			Refrac = sla_refro(_zenith, h, T, _P, RH, wl, _phi, tlr, eps)
 			R.append(np.degrees(Refrac))
 
 	if len(l) == 1:
@@ -326,8 +381,10 @@ def slalib_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15, p=101325,
 				Angle of observation in degrees.
 				zenith = 0 means the observation is directly overhead.
 		conditions 	: string
-			'STANDARD' or 'CERRO_ARMAZONES', refers
+			Refers to the standard conditions of the 
-			to the standard atmospheric conditions.
+			observatory. See get_conditions() for more info.
+			If set to 'None', the conditions specified 
+			after will be used.
 		T 	: 	float (optional)
 			Temperature in Kelvin
 		p 	: 	float (optional)
@@ -343,11 +400,11 @@ def slalib_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15, p=101325,
 		tlr :	float
 			Temperature lapse rate (default=0.0065)
 		eps :	float
-			Precision at which the calculation is aborted (default=1e-10)
+			Precision at which the calculation is aborted (default=1e-8)
 		use_pyslalib : bool
 			Flag that indicates that the user wants to use the pySlalib
 			package. If False, the internal python port of the original
-			FORTRAN code will be used. (Default=False).
+			FORTRAN code will be used. (default=False).
 			
 		Returns
 		-------

AstroAtmosphere/observatories.py

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

AstroAtmosphere/refractionModels.py

@@ -391,7 +391,7 @@ class refraction(Observatory):
 			R : float
 				Atmospheric refraction in degrees
 		'''
-		_R = np.arcsin(n * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
+		_R = np.arcsin(n1 * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
 		return np.degrees(_R)
 
 
@@ -430,7 +430,7 @@ class refraction(Observatory):
 			Parameters
 			----------
 			n1 : float
-				Refractive index at the wavelength of interest
+				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.
@@ -469,10 +469,10 @@ class refraction(Observatory):
 		_R 		= R0 + _h * 1000
 		
 		# Gladstone-Dale relation between refractive index and density
-		_n_h = 1 + (n0 - 1) / _Dens[0] * _Dens 
+		_n_h = 1 + (n1 - 1) / _Dens[0] * _Dens 
 
 		# Calculation of the zenith angle for given height and refractive index, from the refractive invariant
-		_z 		= np.arcsin( (n0 * R0) / (_n_h * _R) * np.sin(zenith) )
+		_z 		= np.arcsin( (n1 * R0) / (_n_h * _R) * np.sin(zenith) )
 		
 		# Calculation of log space derivative of n and r
 		_lnn = np.log(_n_h)

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, RH=0):
+def BarrellAndSears(l, T=288.15, p=101325, RH=0, f=None):
 	'''
 		From Barrell & Sears, 1939
 		Equation 7.7, page 52
@@ -26,15 +26,20 @@ def BarrellAndSears(l, T=288.15, p=101325, RH=0):
 			Atmospheric pressure in Pascal
 		RH 	: 	float
 			Relative humidity (0=<H=<1)
+		f 	: 	float
+			Water vapour pressure in Pascal
 
 		Returns
 		-------
 		n 	: 	float
 			Refractive index of atmospheric air at the given conditions
 	'''
+	if f is None:
+		f = RH * svp(T)
+	f_mmHg = hPa2mmHg(f/100)
 	T = K2C(T)
-	P = hPA2mmHg(p/100)
+	P = hPa2mmHg(p/100)
-	R = (0.378125 + 0.0021414 / l**2 + 0.00001793 / l**4) * (P*(1 + (1.049 - 0.0157*T) * P * 1e-6)/(1 + 0.003661*T)) - (0.0624 - 0.000680 / l**2) * RH / (1 + 0.003661*T)
+	R = (0.378125 + 0.0021414 / l**2 + 0.00001793 / l**4) * (P*(1 + (1.049 - 0.0157*T) * P * 1e-6)/(1 + 0.003661*T)) - (0.0624 - 0.000680 / l**2) * f_mmHg / (1 + 0.003661*T)
 	n = R/1e6 + 1
 	return n
 
@@ -50,6 +55,8 @@ def Ciddor(l, CO2=450):
 		---------- 
 		l 	: 	float
 			Wavelength in microns
+		CO2 : 	float
+			CO2 density in parts per million
 
 		Returns
 		-------
@@ -64,17 +71,27 @@ def Ciddor(l, CO2=450):
 	return n
 
 
-def Edlen1953(l):
+def Edlen1953(l, T=288.15, p=101325, RH=0, f=None):
 	'''
 		From Edlen, 1953
 		Equation 2
 
-		At standard conditions only
+		At standard conditions only.
+		Following Filippenko (1982), eqs. 2 & 3 we also include
+		the	modification proposed by Barrell in 1951.
 
 		Parameters
 		---------- 
 		l 	: 	float
 			Wavelength in microns
+		T 	: 	float
+			Temperature in Kelvin
+		p 	: 	float
+			Atmospheric pressure in Pascal
+		RH 	: 	float
+			Relative humidity (0=<H=<1)
+		f 	: 	float
+			Water vapor pressure in Pascal
 
 		Returns
 		-------
@@ -83,14 +100,27 @@ def Edlen1953(l):
 	'''
 	s = 1 / l
 	R = 6432.8 + 2949810 / (146 - s**2) + 25540/(41-s**2)
-	n = R/1e8 + 1
+	n_standard = R/1e8 + 1
+
+	pHg = hPa2mmHg(p/100)
+	Tc  = K2C(T)
+	if f is None:
+		f = RH * svp(T)
+	f_mmHg = hPa2mmHg(f/100)
+
+	# Changes in temperature and pressure, from Barrell (1951)
+	n_tp = R/1e8 * (pHg * (1 + (1.049 - 0.0157*Tc) * 1e-6 * pHg ) ) / ( 720.883 * (1 + 0.003661 * Tc)) + 1
+
+	# Correction for water vapour pressure
+	n_svp = (n_tp - 1) * 1e6 - ( 0.0624 - 0.000680 * s**2 ) / (1 + 0.003661 * Tc) * f_mmHg
+	n = n_svp/1e6 + 1
 	return n
 
 
-def Edlen1966(l):
+def Edlen1966(l, T=288.15, p=101325, RH=0, CO2=300, f=None):
 	'''
 		From Edlen, 1966
-		Equation 1
+		Equations 1, 12, 17, 22
 
 		At standard conditions only
 
@@ -98,6 +128,16 @@ def Edlen1966(l):
 		---------- 
 		l 	: 	float
 			Wavelength in microns
+		T 	: 	float
+			Temperature in Kelvin
+		p 	: 	float
+			Atmospheric pressure in Pascal
+		RH 	: 	float
+			Relative humidity (0=<H=<1)
+		CO2 : 	float
+			CO2 density in parts per million
+		f 	: 	float
+			Water vapor pressure in Pascal
 
 		Returns
 		-------
@@ -107,15 +147,37 @@ def Edlen1966(l):
 	s = 1 / l
 	R = 8342.13 + 2406030 / (130 - s**2) + 15997 / (38.9 - s**2)
 	n = R/1e8 + 1
+
+	# Influence of temperature and pressure, eq. 12
+	Tc = K2C(T)
+	Ptorr = Pa2Torr(p)
+	R_tp = ( Ptorr * R/1e8 ) / 720.775 * (1 + Ptorr * (0.817 - 0.0133 * Tc) * 1e-6) / ( 1 + 0.0036610 * Tc)
+
+	# Allow for a CO2 fraction, eq. 17, but we include T and P dependence.
+	# This is okay from "Asthe refractivity of CO, is about 50% higher than
+	# that of air this CO,-content increases the refractivity over that of 
+	# CO2-free air by approximately 0.5 x 0.0003 x (n-1)s" and "As the 
+	# density factor of CO, is not much different from that of air this ratio r
+	# can be taken as independent of t and p ."
+	R_x = R_tp * (1 + 0.540 * (CO2/1e6 - 0.0003))
+	n_x = R_x + 1
+
+	# Correction for water vapour, eq. 22
+	if f is None:
+		f = RH * svp(T)
+	ftorr = Pa2Torr(f)
+
+	n = n_x - ftorr * (5.7224 - 0.0457 * s**2) * 1e-8
 	return n
 
 
 def Owens(l, T=288.15, p=101325, RH=0):
 	'''
 		From Owens, 1967
-		Equation 32
+		Equations 29 to 31
-		Note that the pressure P and H are in torricelli (Torr),
+		
-		however, this is almost precisely 1 Torr = 1 mmHg
+		Code heavily based on ScopeSIM (https://scopesim.readthedocs.io/)
+		and http://www.caha.es/newsletter/news03b/pedraz/newslet.html
 
 		Parameters
 		----------
@@ -133,18 +195,35 @@ def Owens(l, T=288.15, p=101325, RH=0):
 		n 	: 	float
 			Refractive index of atmospheric air at the given conditions
 	'''
-	P = hPA2mmHg(p/100)
+	# Calculation from ScopeSIM.effects.atmospheric_refraction()
-	T = K2C(T)
+	pres = p / 100 # from mbar to Pa
-	s = 1 / l
+	Ps = -10474. + 116.43 * T - 0.43284 * T**2 + 0.0005384 * T**3
-	R = 8342.13 + 2406030 / (130 - s**2) + 15997 / (38.9 - s**2) * (P/720.775) * ((1 + P * (0.817 - 0.0133) * 1e-6)/(1 + 0.0036610*T)) - RH * (5.722 - 0.0457*s**2)
+	Pw = Ps * RH / 100.
+	Pa = pres - Pw
+
+	# Wavenumber s
+	s = 1. / l
+
+	# Density factor of dry air, eq. 30
+	Da = (Pa / T) * (1. + Pa * (57.9E-8 - 0.0009325 / T + 0.25844 / T**2))
+	
+	# Density factor of moist air, eq. 31
+	Dw = (Pw / T) * (1. + Pw * (1. + 3.7E-4 * Pw) *
+					 (-2.37321E-3 + 2.23366 / T - 710.792 / T**2 +
+					  77514.1 / T**3))
+
+	# Combining the two into eq. 29
+	na = Da * (2371.34 + 683939.7 / (130. - s**2) + 4547.3 / (38.9 - s**2))
+	nw = Dw * (6487.31 + 58.058 * s**2 - 0.7115 * s**4 + 0.08851 * s**6)
+	R = na + nw
 	n = R/1e8 + 1
 	return n
 
 
-def BonschPotulski(l, T=288.15, p=101325, CO2=450):
+def BonschPotulski(l, T=288.15, p=101325, RH=0, CO2=400, f=None):
 	'''
 		From Bonsch & Potulski, 1998
-		Equations 6, 7, 8
+		Equations 6, 7, 8 and 9a
 		
 		Parameters
 		---------- 
@@ -154,8 +233,12 @@ def BonschPotulski(l, T=288.15, p=101325, CO2=450):
 			Temperature in Kelvin
 		p 	: 	float
 			Atmospheric pressure in Pa
+		RH 	: 	float
+			Relative humidity (0<=H<=1)
 		CO2 : 	float
 			Particles per million of CO2
+		f 	: 	float
+			Water vapor pressure in Pascal
 
 		Returns
 		-------
@@ -163,21 +246,26 @@ def BonschPotulski(l, T=288.15, p=101325, CO2=450):
 			Refractive index of atmospheric air at the given conditions
 	'''
 	x = CO2 / 1e6
-	T = K2C(T)
+	Tc = K2C(T) 
 	s = 1 / l
-	# Rs = 8092.33 + 2333983 / (130 - s**2) + 15518 / (38.9 - s**2)
+
 	Rs = 8091.37 + 2333983 / (130 - s**2) + 15518 / (38.9 - s**2)
 	Rs = Rs / 1e8
 	Rx = Rs * (1 + 0.5327 * (x - 0.0004))
-	Rtp = Rx * p / 93214.60 * (1 + 1e-8 * (0.5953 - 0.009876 * T)*p) / (1 + 0.0036610 * T)
+	Rtp = Rx * p / 93214.60 * (1 + 1e-8 * (0.5953 - 0.009876 * Tc)*p) / (1 + 0.0036610 * Tc)
-	n = Rtp + 1
+
+	if f is None:
+		f = RH * svp(T)
+	
+	n_tp = Rtp + 1
+	n = n_tp - f * (3.8020 - 0.0384 * s**2) * 1e-10
 	return n
 
 
-def BirchDowns(l):
+def BirchDowns(l, T=288.15, p=101325, RH=0, f=None):
 	'''
 		From Birch & Downs, 1993 and 1994 (correction)
-		Equations 1 and 2
+		Equations 1, 2 and 3
 
 		Only at standard conditions
 
@@ -185,15 +273,32 @@ def BirchDowns(l):
 		----------
 		l 	: 	float
 			wavelength in micron
+		T 	: 	float
+			Temperature in Kelvin
+		p 	: 	float
+			Atmospheric pressure in Pa
+		RH 	: 	float
+			Relative humidity (0<=H<=1)
+		f 	: 	float
+			Water vapor pressure in Pascal
 
 		Returns
 		-------
 		n 	: 	float
 			Refractive index of atmospheric air at the given conditions
 	'''
+	Tc = K2C(T)
 	s = 1 / l
-	R = 8342.54 + 2406147 / (130 - s**2) + 15998 / (38.9-s**2)
+	Rs = 8342.54 + 2406147 / (130 - s**2) + 15998 / (38.9-s**2)
-	n = R/1e8 + 1
+	n = Rs/1e8 + 1
+
+	R_tp = (p * Rs/1e8) / 96095.43 * (1 + 1e-8 * (0.601 - 0.00972 * Tc) * p) / (1 + 0.0036610 * Tc) 
+	n_tp = R_tp + 1
+	
+	if f is None:
+		f = RH * svp(T)
+
+	n = n_tp - f * (3.7345 - 0.0401 * s**2) * 1e-10
 	return n
 
 
@@ -368,7 +473,7 @@ def Mathar(l, T=288.15, p=101325, RH=0):
 	if l < 1.3 or l > 24:
 		print('Using the Mathar model at this wavelength (l={:.3f}) is not recommended!'.format(l))
 
-	if l <= 2.65:
+	elif l <= 2.65:
 		########################################################################
 		# model parameters for wavelengths between 1.3 and 2.5 microns
 		cref = [ 0.200192e-3,  0.113474e-9,  -0.424595e-14,  0.100957e-16, -0.293315e-20,  0.307228e-24] # cm^j
@@ -388,7 +493,7 @@ def Mathar(l, T=288.15, p=101325, RH=0):
 		Href = 10          #%, reference humdity in percent
 		########################################################################
 	
-	if l > 2.65 and l <= 4.275:
+	elif l > 2.65 and l <= 4.275:
 		########################################################################
 		# model parameters for wavelengths between 2.8 and 4.2 microns
 		cref = [ 0.200049e-3,  0.145221e-9,   0.250951e-12, -0.745834e-15, -0.161432e-17,  0.352780e-20] # cm^j
@@ -408,7 +513,7 @@ def Mathar(l, T=288.15, p=101325, RH=0):
 		Href = 10          #%, reference humdity in percent
 		########################################################################
 
-	if l > 4.275 and l <= 6.35:
+	elif l > 4.275 and l <= 6.35:
 		########################################################################
 		# model parameters for wavelengths between 4.35 and 5.2 microns
 		cref = [ 0.200020e-3,  0.275346e-9,   0.325702e-12, -0.693603e-14,  0.285610e-17,  0.338758e-18] # cm^j
@@ -427,7 +532,7 @@ def Mathar(l, T=288.15, p=101325, RH=0):
 		Href = 10          #%, reference humdity in percent
 		########################################################################
 
-	if l > 6.35 and l <= 15.05:
+	elif l > 6.35 and l <= 15.05:
 		########################################################################
 		# model parameters for wavelengths between 7.5 and 14.1 microns
 		cref = [ 0.199885e-3,  0.344739e-9,  -0.273714e-12,  0.393383e-15, -0.569488e-17,  0.164556e-19] # cm^j
@@ -447,7 +552,7 @@ def Mathar(l, T=288.15, p=101325, RH=0):
 		########################################################################
 
 
-	if l > 15.05:
+	elif l > 15.05:
 		########################################################################
 		# model parameters for wavelengths between 16 and 24 microns
 		cref = [ 0.199436e-3,  0.299123e-8,  -0.214862e-10,  0.143338e-12,  0.122398e-14, -0.114628e-16] # cm^j
@@ -481,6 +586,7 @@ def Mathar(l, T=288.15, p=101325, RH=0):
 	return n 
 
 
+# Other necessary conversion or calculations
 def K2C(K):
 	'''
 		Kelvin to Celsius conversion
@@ -498,7 +604,7 @@ def K2C(K):
 	return K - 273.15
 
 
-def hPA2mmHg(hPa):
+def hPa2mmHg(hPa):
 	'''
 		hPa to mmHg conversion
 
@@ -514,3 +620,41 @@ def hPA2mmHg(hPa):
 	'''
 	return 0.7500616827 * hPa
 
+
+def Pa2Torr(Pa):
+	'''
+		Pa to Torr conversion
+
+		Parameters
+		----------
+		Pa 	:	float
+			Atmospheric pressure in Pa (Pascal)
+
+		Returns
+		-------
+		Torr 	: 	float
+			Atmospheric pressure in Torr
+	'''
+
+	return Pa * 760/101325
+
+
+def svp(T):
+	'''
+		Saturation vapour pressure in Pa.
+
+		Parameters
+		----------
+		T 	: float
+			Atmospheric temperature in Kelvin
+
+		Returns
+		-------
+		svp	: float
+			Saturation vapour pressure in Pa
+	'''
+	_A =  1.2378847e-5 # K^-2
+	_B = -1.9121316e-2 # K^-2
+	_C = 33.93711047
+	_D = -6.3431645e3  # K^-2
+	return np.exp(_A*T**2 + _B*T + _C + _D/T)

MANIFEST.in

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

README.md

@@ -61,10 +61,13 @@ The package contains the following files:
 - `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
 - `misc.py` 				    : Contains several miscellaneous functions, e.g. to quickly calculate the refractive index.
+- `neoslalib.py`                : Contains a Python native port of the `pyslalib.SLA_REFRO()` routine, translated from the original FORTRAN code.
+- `observatories.py` 			: Contains a dictionary with the average atmospheric conditions for various large observatories around the world.
 - `refractionModels.py`		 	: Contains various atmospheric models that can be used to calculate the atmospheric refraction
 - `refractivityModels.py` 		: Contains various functions from literature to calculate the refractive index. Some examples of works included are Ciddor (1996), Barrell & Sears (1939), SLALIB, Owens (1967) and Edlen (1966).
 
 
+
 ## The Ciddor refractivity model
 
 The refractivity model from Ciddor, 1996, calculates the index of refraction of air using the following parameters:
@@ -107,7 +110,7 @@ Input parameters for the spherical atmosphere models (all but the first):
 - xc - CO2 density              [ppm]
 - lat - Latitude                [deg]
 - h - Height above sea level    [m]
-- l - Wavelength(s)             [um]
+- l - Wavelength(s)             [μm]
 - z - Zenith angle              [deg]
 
 An example - To calculate the atmospheric refraction we set up a `refraction()` object. This requires the latitude and altitude of the observer. The atmospheric dispersion is calculated similarly, but by exchanging the `refraction()` object by a `dispersion()` object and parsing two refractive indices into the `cassini()` function.
@@ -118,7 +121,7 @@ rho = at.rho(p=101325, T=273.15, RH=0.2, xc=300)
 
 # Initializing refraction model and setting the reduced height
 ref = refraction(lat=30, h=2200)
-ref.setReducedHeight(p=101325, rho)
+ref.setReducedHeight(p=101325, rho=rho)
 
 # Calculating the atmospheric refraction
 atm_ref = ref.cassini(n, zenith=30)
@@ -186,12 +189,13 @@ Included are
 | ---------- | -------- | ------- |
 | Barrell & Sears (1939) | `BarrellAndSears()` | Outdated model, especially poor for short wavelengths. |
 | Ciddor (1996) | `Ciddor()` | Only at standard atmosphere. CO2 density can be changed. |
-| Edlen (1953) | `Edlen1953()` | Only at standard atmosphere. Function requires wavelength only. |
+| Edlen (1953) | `Edlen1953()` | The equation by Edlèn, 1953, combined with modification from Barrel (1951), to include temperature, pressure and humidity dependence. |
-| Edlen (1966) | `Edlen1966()` | Only at standard atmosphere. Function requires wavelength only. |
+| Edlen (1966) | `Edlen1966()` | Allows for any temperature, pressure, humidity and CO2 density.|
-| Owens (1967) | `Owens()` | - |
+| Owens (1967) | `Owens()` | Allows for any temperature, pressure, humidity. |
-| Bonsch & Potulski (1998) | `BonschPotulski()` | Function does not include humidity |
+| Bonsch & Potulski (1998) | `BonschPotulski()` | Allows for any temperature, pressure, humidity and CO2 density. |
 | Birch & Downs (1994) | `BirchDowns()` | Only at standard atmosphere. Function requires wavelength only. |
 | Peck & Reeder (1972) | `PeckReeder()` | Only at standard atmosphere. Function requires wavelength only. |
+| Mathar (2006) | `Mathar()` | Allows for temperature, pressure and humidity specification. Works best at wavelengths above 1.4 μm. |
 | Hohenkerk & Sinclair (1985) | `HohenkerkAndSinclair()` | Based on Barrell & Sears equation. Generally, gives slightly larger values.|
 | Hohenkerk & Sinclair (1985) | `slalib()` | Uses pySLALIB to retrieve the refractive index at standard atmospheric conditions.|
 
@@ -199,9 +203,8 @@ The last function extracts the refractive index from the `pyslalib.sla_refro()`
 
 ## Quick functions
 
-Some functions that assume default conditions for the atmosphere. Then only the wavelength and zenith angle are required to calculate the refraction and dispersion. 
+Some functions can make use of the `observatories` dictionary as input. These include default average conditions of the atmosphere at these locations. When these are used, only the wavelength and zenith angle are required to calculate the refraction and dispersion. 
 
-Presently, only the `'STANDARD'` and `'CERRO_ARMAZONES'` flags are known. Contact the author to add additional standard conditions.
 
 ```python
  # Quick refractive index
@@ -214,6 +217,20 @@ rq = quick_refraction(l1, zenith_obs, conditions='STANDARD')
 dq = quick_dispersion(l1, l2, zenith_obs, conditions='STANDARD')
 ```
 
+Presently, several observatories are known. These can all be found in the `observatories` dictionary. Contact the author to add additional standard conditions.
+
+| Observatory                                   | Flag              | Notes                                         | 
+| --------------------------------------------- | ----------------- | --------------------------------------------- |
+| Standard conditions                           | `STANDARD`        |                                               |
+| Cerro Armazones                               | `CERRO_ARMAZONES` | Location of the future ELT                    |
+| Cerro Paranal                                 | `CERRO_PARANAL`   | Location of the VLT, VISTA, VST, NGTS and SSO |
+| Roque de los Muchachos Observatory (La Palma) | `LA_PALMA`        | GTC, INT, WHT, NOT, and others.               |
+| La Silla Observatory                          | `LA_SILLA`        | ESO 3.6m telescope, NTG, BlackGEM and others. | 
+| Las Campanas Observatory                      | `LAS_CAMPANAS`    | Site of the future GMT.                       |
+| Mauna Kea Observatory (Hawaii)                | `MAUNA_KEA`       | Location of Keck, Subaru, Gemini North and others.|
+
 ## Other
 
-This package also contains a python native port of the `SLA_REFRO()` routine of the pyslalib package, which was originally written in FORTRAN by P.T. Wallace. It can be called directly using `sla_refro()`. When using the same units as the other refraction functions, you might prefer using `slalib_refraction()` or `slalib_dispersion()`.
+This package also contains a python native port of the `SLA_REFRO()` routine of the pyslalib package, which was originally written in FORTRAN by P.T. Wallace. It can be called directly using `sla_refro()`. When using the same units as the other refraction functions, you might prefer using `slalib_refraction()` or `slalib_dispersion()`. Do note that this implementation is slower than the (compiled) FORTRAN implementation.
+
+Finally, the dispersion model used by Alexei Filippenko in 1982 is also included. This model is a combination of the (modified) Edlèn equation and a plane-parallel atmospheric geometry. This calculation can be accessed by calling the `Filippenko1982()` function.

scripts/compare_refractivity_models.py

@@ -81,7 +81,26 @@ class atmosphericModel:
 	def slalib(self, l=None):
 		if not l:
 			l = self.lrange
+		try:
+			from pyslalib import slalib
 			return(AstroAtmosphere.slalib(l, T=self.T, p=self.P, RH=self.H)-1)
+		except ImportError as e:
+			print('Could not import pyslalib.slalib. Using slower Python native function sla_refro().')
+			h=0
+			lat=-24
+			tlr=0.0065
+			eps=1e-10
+			zObs = np.deg2rad(45.07182)
+			p = self.P/100
+			R = []
+			if type(l) == float:
+				l = [l]
+			for wl in l:
+				R.append(AstroAtmosphere.sla_refro(zObs, h, self.T, p, self.H, wl, lat, tlr, eps))
+			if len(l) == 1:
+				return R[0]
+			else:
+				return np.asarray(R)
 
 
 if __name__ == '__main__':

setup.py

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