Fixed the setup

.gitignore

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

AstroAtmosphere/__init__.py

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

AstroAtmosphere/dispersionModels.py

@@ -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/examples.py

@@ -1,207 +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, zenith, 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
-		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
-		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, zenith)
-
-
-def quick_dispersion(l1, l2, zenith, 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.
-		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
-		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, zenith)

AstroAtmosphere/misc.py

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

AstroAtmosphere/neoslalib.py

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

AstroAtmosphere/observatories.py

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

AstroAtmosphere/refractionModels.py

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

MANIFEST.in

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

README.md

@@ -41,7 +41,7 @@ disp = dispersion(lat, h)
 disp.setReducedHeight(P, rho)
 
 # Calculation of the atmopheric dipsersion
-atm_dispersion = disp.cassini(n1, z0) * 3.6e6
+atm_dispersion = disp.cassini(n1, n2, z0) * 3.6e6
 print ('The dispersion is %.03f milli arc seconds' %(atm_dispersion))
 ```
 
@@ -60,11 +60,14 @@ The package contains the following files:
 - `/data/1976USSA.txt` 			: Data file of the 1976 US Standard Atmosphere.
 - `ciddorModel.py` 				: Contains the `Observatory()` class to calculate the refractive index of atmospheric air, using the Ciddor (1996) model.
 - `dispersionModels.py` 		: Contains various atmospheric models that can be used to calculate the atmospheric dispersion
-- `examples.py` 				: Contains several quick functions, e.g. to quickly calculate the refractive index.
+- `misc.py` 				    : Contains several miscellaneous functions, e.g. to quickly calculate the refractive index.
+- `neoslalib.py`                : Contains a Python native port of the `pyslalib.SLA_REFRO()` routine, translated from the original FORTRAN code.
+- `observatories.py` 			: Contains a dictionary with the average atmospheric conditions for various large observatories around the world.
 - `refractionModels.py`		 	: Contains various atmospheric models that can be used to calculate the atmospheric refraction
 - `refractivityModels.py` 		: Contains various functions from literature to calculate the refractive index. Some examples of works included are Ciddor (1996), Barrell & Sears (1939), SLALIB, Owens (1967) and Edlen (1966).
 
 
+
 ## The Ciddor refractivity model
 
 The refractivity model from Ciddor, 1996, calculates the index of refraction of air using the following parameters:
@@ -85,7 +88,7 @@ Basic functionality of other refractivity models is included as well and is disc
 
 ## The atmospheric refraction and atmospheric dispersion models
 
-If the atmospheric refraction is denoted as R(n), then atmospheric dispersion is defined as R(n1) - R(n2). Thus they are very similar in function and also in use.
+If the atmospheric refraction is denoted as R(n), then atmospheric dispersion is defined as R(n1) - R(n2). Thus they are very similar in function and use.
 
 
 The available models in this package:
@@ -97,6 +100,8 @@ The available models in this package:
 - `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):
 - T - Temperature               [K]
@@ -105,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.
@@ -116,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)
@@ -184,21 +189,22 @@ Included are
 | ---------- | -------- | ------- |
 | Barrell & Sears (1939) | `BarrellAndSears()` | Outdated model, especially poor for short wavelengths. |
 | Ciddor (1996) | `Ciddor()` | Only at standard atmosphere. CO2 density can be changed. |
-| Edlen (1953) | `Edlen1953()` | Only at standard atmosphere. Function requires wavelength only. |
+| Edlen (1953) | `Edlen1953()` | The equation by Edlèn, 1953, combined with modification from Barrel (1951), to include temperature, pressure and humidity dependence. |
-| Edlen (1966) | `Edlen1966()` | Only at standard atmosphere. Function requires wavelength only. |
+| Edlen (1966) | `Edlen1966()` | Allows for any temperature, pressure, humidity and CO2 density.|
-| Owens (1967) | `Owens()` | - |
+| Owens (1967) | `Owens()` | Allows for any temperature, pressure, humidity. |
-| Bonsch & Potulski (1998) | `BonschPotulski()` | Function does not include humidity |
+| Bonsch & Potulski (1998) | `BonschPotulski()` | Allows for any temperature, pressure, humidity and CO2 density. |
 | Birch & Downs (1994) | `BirchDowns()` | Only at standard atmosphere. Function requires wavelength only. |
 | Peck & Reeder (1972) | `PeckReeder()` | Only at standard atmosphere. Function requires wavelength only. |
-| Hohenkerk & Sinclair (1985) | `slalib()` | Only at standard atmosphere. Function requires wavelength only. |
+| Mathar (2006) | `Mathar()` | Allows for temperature, pressure and humidity specification. Works best at wavelengths above 1.4 μm. |
+| Hohenkerk & Sinclair (1985) | `HohenkerkAndSinclair()` | Based on Barrell & Sears equation. Generally, gives slightly larger values.|
+| Hohenkerk & Sinclair (1985) | `slalib()` | Uses pySLALIB to retrieve the refractive index at standard atmospheric conditions.|
 
 The last function extracts the refractive index from the `pyslalib.sla_refro()` function and thus requires [pyslalib](https://pypi.org/project/pySLALIB/). This is an optional dependency for the setup of AstroAtmosphere.
 
 ## Quick functions
 
-Some functions that assume default conditions for the atmosphere. Then only the wavelength and zenith angle are required to calculate the refraction and dispersion. 
+Some functions can make use of the `observatories` dictionary as input. These include default average conditions of the atmosphere at these locations. When these are used, only the wavelength and zenith angle are required to calculate the refraction and dispersion. 
 
-Presently, only the `'STANDARD'` and `'CERRO_ARMAZONES'` flags are known. Contact the author to add additional standard conditions.
 
 ```python
  # Quick refractive index
@@ -210,3 +216,21 @@ rq = quick_refraction(l1, zenith_obs, conditions='STANDARD')
 # Quick dispersion
 dq = quick_dispersion(l1, l2, zenith_obs, conditions='STANDARD')
 ```
+
+Presently, several observatories are known. These can all be found in the `observatories` dictionary. Contact the author to add additional standard conditions.
+
+| Observatory                                   | Flag              | Notes                                         | 
+| --------------------------------------------- | ----------------- | --------------------------------------------- |
+| Standard conditions                           | `STANDARD`        |                                               |
+| Cerro Armazones                               | `CERRO_ARMAZONES` | Location of the future ELT                    |
+| Cerro Paranal                                 | `CERRO_PARANAL`   | Location of the VLT, VISTA, VST, NGTS and SSO |
+| Roque de los Muchachos Observatory (La Palma) | `LA_PALMA`        | GTC, INT, WHT, NOT, and others.               |
+| La Silla Observatory                          | `LA_SILLA`        | ESO 3.6m telescope, NTG, BlackGEM and others. | 
+| Las Campanas Observatory                      | `LAS_CAMPANAS`    | Site of the future GMT.                       |
+| Mauna Kea Observatory (Hawaii)                | `MAUNA_KEA`       | Location of Keck, Subaru, Gemini North and others.|
+
+## Other
+
+This package also contains a python native port of the `SLA_REFRO()` routine of the pyslalib package, which was originally written in FORTRAN by P.T. Wallace. It can be called directly using `sla_refro()`. When using the same units as the other refraction functions, you might prefer using `slalib_refraction()` or `slalib_dispersion()`. Do note that this implementation is slower than the (compiled) FORTRAN implementation.
+
+Finally, the dispersion model used by Alexei Filippenko in 1982 is also included. This model is a combination of the (modified) Edlèn equation and a plane-parallel atmospheric geometry. This calculation can be accessed by calling the `Filippenko1982()` function.

scripts/compare_refractivity_models.py

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

setup.py

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