Fixed the setup

.gitignore

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

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

@@ -4,8 +4,9 @@ from .refractionModels import refraction
 from .dispersionModels import dispersion
 from .neoslalib import sla_refro
 from .observatories import observatories
-from .refractivityModels import Edlen1953
+from .refractivityModels import Edlen1953, Mathar
 import numpy as np
+import warnings
 
 '''
 @ Author: 		Joost van den Born
@@ -102,7 +103,9 @@ 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
 			Refers to the standard conditions of the 
 			observatory. See get_conditions() for more info.
@@ -129,11 +132,15 @@ def quick_refractive_index(l, conditions='STANDARD', T=288.15, p=101325,
 	if conditions is not None:
 		T, p, RH, xc, lat, h, _ = get_conditions(conditions).values()
 
-	# Initializing dispersion model
-	at  = Observatory()
+	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)
+		# 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
 
 
@@ -147,7 +154,8 @@ 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.
@@ -180,8 +188,12 @@ def quick_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325,
 	# 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)
+	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)
@@ -202,9 +214,11 @@ 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.
@@ -237,9 +251,21 @@ def quick_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15,
 	# 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)
+	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)
@@ -310,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:
@@ -319,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:

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

@@ -36,9 +36,9 @@ def BarrellAndSears(l, T=288.15, p=101325, RH=0, f=None):
 	'''
 	if f is None:
 		f = RH * svp(T)
-	f_mmHg = hPA2mmHg(f/100)
+	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) * f_mmHg / (1 + 0.003661*T)
 	n = R/1e6 + 1
 	return n
@@ -55,6 +55,8 @@ def Ciddor(l, CO2=450):
 		---------- 
 		l 	: 	float
 			Wavelength in microns
+		CO2 : 	float
+			CO2 density in parts per million
 
 		Returns
 		-------
@@ -86,6 +88,8 @@ def Edlen1953(l, T=288.15, p=101325, RH=0, f=None):
 			Temperature in Kelvin
 		p 	: 	float
 			Atmospheric pressure in Pascal
+		RH 	: 	float
+			Relative humidity (0=<H=<1)
 		f 	: 	float
 			Water vapor pressure in Pascal
 
@@ -98,11 +102,11 @@ def Edlen1953(l, T=288.15, p=101325, RH=0, f=None):
 	R = 6432.8 + 2949810 / (146 - s**2) + 25540/(41-s**2)
 	n_standard = R/1e8 + 1
 
-	pHg = hPA2mmHg(p/100)
+	pHg = hPa2mmHg(p/100)
 	Tc  = K2C(T)
 	if f is None:
 		f = RH * svp(T)
-	f_mmHg = hPA2mmHg(f/100)
+	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
@@ -113,10 +117,10 @@ def Edlen1953(l, T=288.15, p=101325, RH=0, f=None):
 	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
 
@@ -124,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
 		-------
@@ -133,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
-		Note that the pressure P and H are in torricelli (Torr),
-		however, this is almost precisely 1 Torr = 1 mmHg
+		Equations 29 to 31
+		
+		Code heavily based on ScopeSIM (https://scopesim.readthedocs.io/)
+		and http://www.caha.es/newsletter/news03b/pedraz/newslet.html
 
 		Parameters
 		----------
@@ -159,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)
-	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)) - RH * (5.722 - 0.0457*s**2)
+	# Calculation from ScopeSIM.effects.atmospheric_refraction()
+	pres = p / 100 # from mbar to Pa
+	Ps = -10474. + 116.43 * T - 0.43284 * T**2 + 0.0005384 * T**3
+	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
 		---------- 
@@ -180,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
 		-------
@@ -189,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)
-	n = Rtp + 1
+	Rtp = Rx * p / 93214.60 * (1 + 1e-8 * (0.5953 - 0.009876 * Tc)*p) / (1 + 0.0036610 * Tc)
+
+	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
 
@@ -211,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)
-	n = R/1e8 + 1
+	Rs = 8342.54 + 2406147 / (130 - s**2) + 15998 / (38.9-s**2)
+	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
 
 
@@ -348,10 +427,10 @@ def slalib(l, T=288.15, p=101325, RH=0, h=0, lat=-24, tlr=0.0065, eps=1e-10):
 			conditions.
 	'''
 	try:
-	    from pyslalib import slalib
+		from pyslalib import slalib
 	except ImportError as e:
-	    print('pySLALIB package not found! Try installing it first.\nReturning n=1.')
-	    return 1
+		print('pySLALIB package not found! Try installing it first.\nReturning n=1.')
+		return 1
 	
 	zObs = np.deg2rad(45.06563)
 	p = p/100
@@ -394,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
@@ -414,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
@@ -434,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
@@ -453,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
@@ -473,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
@@ -490,7 +569,7 @@ def Mathar(l, T=288.15, p=101325, RH=0):
 		σref = 1e4/20      # cm^−1, reference wavenumber
 		Tref = 273.15+17.5 # K, reference temperature
 		pref = 75000       # Pa, reference pressure
-		Href = 10          #%, reference humdity in percent
+		Href = 10          # %, reference humdity in percent
 		########################################################################
 
 	s = 1e4/l # cm^-1
@@ -525,7 +604,7 @@ def K2C(K):
 	return K - 273.15
 
 
-def hPA2mmHg(hPa):
+def hPa2mmHg(hPa):
 	'''
 		hPa to mmHg conversion
 
@@ -542,6 +621,24 @@ 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.

README.md

@@ -60,10 +60,11 @@ 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
-- `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
+- `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).
-- `observatories` 				: Contains a dictionary with the average atmospheric conditions for various large observatories around the world.
 
 
 
@@ -109,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.
@@ -120,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)
@@ -189,11 +190,12 @@ 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()` | 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. |
-| Owens (1967) | `Owens()` | - |
-| Bonsch & Potulski (1998) | `BonschPotulski()` | Function does not include humidity |
+| Edlen (1966) | `Edlen1966()` | Allows for any temperature, pressure, humidity and CO2 density.|
+| Owens (1967) | `Owens()` | Allows for any temperature, pressure, 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.|
 
@@ -225,9 +227,10 @@ Presently, several observatories are known. These can all be found in the `obser
 | 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
-		return(AstroAtmosphere.slalib(l, T=self.T, p=self.P, RH=self.H)-1)
+		try:
+			from pyslalib import slalib
+			return(AstroAtmosphere.slalib(l, T=self.T, p=self.P, RH=self.H)-1)
+		except ImportError as e:
+			print('Could not import pyslalib.slalib. Using slower Python native function sla_refro().')
+			h=0
+			lat=-24
+			tlr=0.0065
+			eps=1e-10
+			zObs = np.deg2rad(45.07182)
+			p = self.P/100
+			R = []
+			if type(l) == float:
+				l = [l]
+			for wl in l:
+				R.append(AstroAtmosphere.sla_refro(zObs, h, self.T, p, self.H, wl, lat, tlr, eps))
+			if len(l) == 1:
+				return R[0]
+			else:
+				return np.asarray(R)
 
 
 if __name__ == '__main__':

setup.py

@@ -21,14 +21,14 @@ def setup_package():
     # write_version_py()
 
     setup(name='AstroAtmosphere',
-          version='1.3',
+          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.3/atmosphericmodel-v1.3.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'},