Fixed Owens (1967) model

Update README.md
Improved quick functions
2026-10-30 05:01:02 +00:00 · 2022-07-06 12:19:31 +02:00 · 2022-03-23 08:55:21 +00:00 · 2022-03-23 08:42:45 +01:00 · 2021-03-10 12:27:07 +01:00 · 2021-03-10 12:04:27 +01:00
10 changed files with 473 additions and 149 deletions
@@ -8,3 +8,4 @@ from .refractionModels import *
 from .misc import *
 from .refractivityModels import *
 from .neoslalib import *
+from .observatories import *
@@ -3,12 +3,94 @@ from .ciddorModel import Observatory
 from .refractionModels import refraction
 from .dispersionModels import dispersion
 from .neoslalib import sla_refro
+from .observatories import observatories
+from .refractivityModels import Edlen1953, Mathar
+import numpy as np
+import warnings

 '''
@ Author: 		Joost van den Born
@ Contact: 		born@astron.nl
@ Description: 	Some functions for quick testing.
 '''
+def get_conditions(location):
+	'''
+		Returns the mean conditions at various site. Currently
+		supported locations:
+		- 'STANDARD'
+		- 'CERRO_ARMAZONES'
+		- 'CERRO_PARANAL'
+		- 'LA_PALMA'
+		- 'LA_SILLA'
+		- 'LAS_CAMPANAS'
+		- 'MAUNA_KEA'
+
+		Parameters
+		----------
+		location 	: string
+			String denoting the desired location
+
+		Returns
+		-------
+		observatory : dict
+			Dictionary containing the atmospheric parameters
+	'''
+	observatory = observatories[location]
+	return observatory
+
+
+def Filippenko1982(l1, l2, zenith, conditions='STANDARD', T=288.15, p=101325,
+							RH=0.0, xc=450, lat=0, h=0, f=None):
+	'''
+		Calculates the refraction following the 1982 paper by 
+		Alexei V. Filippenko.
+
+		The refraction is based on a flat atmosphere and the refractivity
+		values are calculated from a modified Edlen (1953) equation.
+
+		Parameters
+		----------
+		l1 	: 	float
+			Shortest wavelength of interest in micron.
+		l2 	: 	float
+			Longest wavelength of interest in micron.
+		zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+		conditions 	: string
+			Refers to the standard conditions of the 
+			observatory. See get_conditions() for more info.
+			If set to 'None', the conditions specified 
+			after will be used.
+		T 	: float (optional)
+			Temperature in Kelvin
+		p 	: float (optional)
+			Atmospheric pressure in Pa
+		RH 	: float (optional)
+			Relative humidity
+		xc 	: float (optional)
+			CO2 density in parts per million
+		lat : float (optional)
+			Latitude of the observer in degrees
+		h 	: float (optional)
+			Altitude of the observer in meters
+		f 	: float
+			Water vapour pressure (default=None). If None, it is
+			calculated automatically.
+
+		Returns
+		-------
+		R : float
+			Atmospheric refraction in degrees
+	'''
+
+	if conditions is not None:
+		T, p, RH, xc, lat, h, _ = get_conditions(conditions).values()
+
+	n1 = Edlen1953(l1, T=T, p=p, RH=RH, f=f)
+	n2 = Edlen1953(l2, T=T, p=p, RH=RH, f=f)
+	dR  = (n1 - n2) * np.tan(np.deg2rad(zenith))
+	return np.degrees(dR)


 def quick_refractive_index(l, conditions='STANDARD', T=288.15, p=101325, 
@@ -21,10 +103,14 @@ def quick_refractive_index(l, conditions='STANDARD', T=288.15, p=101325,
 		Parameters
 		----------
 		l 	: float
-			Wavelength in microns
+			Wavelength in microns. If the wavelength is above
+			2.5 microns, then the Mathar refractivity model
+			is used.
 		conditions: string
-			'STANDARD' or 'CERRO_ARMAZONES', refers
-			to the standard atmospheric conditions.
+			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)
@@ -43,29 +129,18 @@ def quick_refractive_index(l, conditions='STANDARD', T=288.15, p=101325,
 		n 	: 	float
 			Refractive index for the given conditions.
 	'''
-	if conditions == 'STANDARD':
-		# Parameters at Cerro Armazones
-		T   = 288.15    # k
-		p   = 101325    # Pa
-		RH  = 0.0
-		xc  = 450       # ppm
-		lat = 0  		# degrees
-		h   = 0     	# m
-		
-	elif conditions == 'CERRO_ARMAZONES':
-		# Parameters at Cerro Armazones
-		T   = 279.65    # k
-		p   = 71200     # Pa
-		RH  = 0.22
-		xc  = 450       # ppm
-		lat = -24.5983  # degrees
-		h   = 3064      # m
+	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


@@ -79,13 +154,16 @@ def quick_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325,
 		Parameters
 		----------
 		l  	: 	float 
-			Wavelength in microns
+			Wavelength in microns. If the wavelength is above 2.5 micron
+			then the Mathar refractivity model is used.
 		zenith : float
 				Angle of observation in degrees.
 				zenith = 0 means the observation is directly overhead.
 		conditions 	: string
-			'STANDARD' or 'CERRO_ARMAZONES', refers
-			to the standard atmospheric conditions.
+			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)
@@ -104,29 +182,18 @@ def quick_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325,
 		R 	: 	float
 			Refraction in degrees
 	'''
-	if conditions == 'STANDARD':
-		# Parameters at Cerro Armazones
-		T   = 288.15    # k
-		p   = 101325    # Pa
-		RH  = 0.0
-		xc  = 450       # ppm
-		lat = 0  		# degrees
-		h   = 0     	# m
-		
-	elif conditions == 'CERRO_ARMAZONES':
-		# Parameters at Cerro Armazones
-		T   = 279.65    # k
-		p   = 71200     # Pa
-		RH  = 0.22
-		xc  = 450       # ppm
-		lat = -24.5983  # degrees
-		h   = 3064      # m
+	if conditions is not None:
+		T, p, RH, xc, lat, h, _ = get_conditions(conditions).values()

 	# 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)
@@ -147,15 +214,19 @@ def quick_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15,
 		Parameters
 		----------
 		l1 	: 	float
-			Shortest wavelength of interest in micron.
+			Shortest wavelength of interest in micron. If the wavelength
+			is above 2.5 microns, then the Mathar model is used.
 		l2 	: 	float
-			Longest wavelength of interest in micron.
+			Longest wavelength of interest in micron. If the wavelength
+			is above 2.5 microns, then the Mathar model is used.
 		zenith : float
 				Angle of observation in degrees.
 				zenith = 0 means the observation is directly overhead.
 		conditions 	: string
-			'STANDARD' or 'CERRO_ARMAZONES', refers
-			to the standard atmospheric conditions.
+			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)
@@ -174,30 +245,27 @@ def quick_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15,
 		dR 	: float
 		 	Atmospheric dispersion in degrees
 	'''
-	if conditions == 'STANDARD':
-		# Parameters at Cerro Armazones
-		T   = 288.15    # k
-		p   = 101325    # Pa
-		RH  = 0.0
-		xc  = 450       # ppm
-		lat = 0  		# degrees
-		h   = 0     	# m
-		
-	elif conditions == 'CERRO_ARMAZONES':
-		# Parameters at Cerro Armazones
-		T   = 279.65    # k
-		p   = 71200     # Pa
-		RH  = 0.22
-		xc  = 450       # ppm
-		lat = -24.5983  # degrees
-		h   = 3064      # m
+	if conditions is not None:
+		T, p, RH, xc, lat, h, _ = get_conditions(conditions).values()

 	# 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)
@@ -228,8 +296,10 @@ def slalib_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325,
 				Angle of observation in degrees.
 				zenith = 0 means the observation is directly overhead.
 		conditions 	: string
-			'STANDARD' or 'CERRO_ARMAZONES', refers
-			to the standard atmospheric conditions.
+			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)
@@ -245,11 +315,11 @@ def slalib_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325,
 		tlr :	float
 			Temperature lapse rate (default=0.0065)
 		eps :	float
-			Precision at which the calculation is aborted (default=1e-10)
+			Precision at which the calculation is aborted (default=1e-8)
 		use_pyslalib : bool
 			Flag that indicates that the user wants to use the pySlalib
 			package. If False, the internal python port of the original
-			FORTRAN code will be used. (Default=False).
+			FORTRAN code will be used. (default=False).
 			
 		Returns
 		-------
@@ -257,23 +327,8 @@ def slalib_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325,
 			Refraction in degrees
 	'''

-	if conditions == 'STANDARD':
-		# Parameters at Cerro Armazones
-		T   = 288.15    # k
-		p   = 101325    # Pa
-		RH  = 0.0
-		xc  = 450       # ppm
-		lat = 0  		# degrees
-		h   = 0     	# m
-		
-	elif conditions == 'CERRO_ARMAZONES':
-		# Parameters at Cerro Armazones
-		T   = 279.65    # k
-		p   = 71200     # Pa
-		RH  = 0.22
-		xc  = 450       # ppm
-		lat = -24.5983  # degrees
-		h   = 3064      # m
+	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
@@ -326,8 +381,10 @@ def slalib_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15, p=101325,
 				Angle of observation in degrees.
 				zenith = 0 means the observation is directly overhead.
 		conditions 	: string
-			'STANDARD' or 'CERRO_ARMAZONES', refers
-			to the standard atmospheric conditions.
+			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)
@@ -343,11 +400,11 @@ def slalib_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15, p=101325,
 		tlr :	float
 			Temperature lapse rate (default=0.0065)
 		eps :	float
-			Precision at which the calculation is aborted (default=1e-10)
+			Precision at which the calculation is aborted (default=1e-8)
 		use_pyslalib : bool
 			Flag that indicates that the user wants to use the pySlalib
 			package. If False, the internal python port of the original
-			FORTRAN code will be used. (Default=False).
+			FORTRAN code will be used. (default=False).
 			
 		Returns
 		-------
@@ -33,7 +33,7 @@ def sla_atms(RT, TT, DNT, GAMAL, R):
 			R * rate the refractive index is changing at R
 	'''
 	_B 		= GAMAL / TT
-	_W 		= (DNT - 1) * np.exp( -B * (R - RT) )
+	_W 		= (DNT - 1) * np.exp( -_B * (R - RT) )
 	DN  	= 1+_W
 	RDNDR 	= -R * _B * _W
 	return DN, RDNDR
@@ -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.'
+					#	}
+					}
@@ -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)
@@ -9,7 +9,7 @@ import numpy as np
 				atmospheric dispersion, given a set of conditions.
 '''

-def BarrellAndSears(l, T=288.15, p=101325, RH=0):
+def BarrellAndSears(l, T=288.15, p=101325, RH=0, f=None):
 	'''
 		From Barrell & Sears, 1939
 		Equation 7.7, page 52
@@ -26,15 +26,20 @@ def BarrellAndSears(l, T=288.15, p=101325, RH=0):
 			Atmospheric pressure in Pascal
 		RH 	: 	float
 			Relative humidity (0=<H=<1)
+		f 	: 	float
+			Water vapour pressure in Pascal

 		Returns
 		-------
 		n 	: 	float
 			Refractive index of atmospheric air at the given conditions
 	'''
+	if f is None:
+		f = RH * svp(T)
+	f_mmHg = hPa2mmHg(f/100)
 	T = K2C(T)
-	P = hPA2mmHg(p/100)
-	R = (0.378125 + 0.0021414 / l**2 + 0.00001793 / l**4) * (P*(1 + (1.049 - 0.0157*T) * P * 1e-6)/(1 + 0.003661*T)) - (0.0624 - 0.000680 / l**2) * RH / (1 + 0.003661*T)
+	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

@@ -50,6 +55,8 @@ def Ciddor(l, CO2=450):
 		---------- 
 		l 	: 	float
 			Wavelength in microns
+		CO2 : 	float
+			CO2 density in parts per million

 		Returns
 		-------
@@ -64,17 +71,27 @@ def Ciddor(l, CO2=450):
 	return n


-def Edlen1953(l):
+def Edlen1953(l, T=288.15, p=101325, RH=0, f=None):
 	'''
 		From Edlen, 1953
 		Equation 2

-		At standard conditions only
+		At standard conditions only.
+		Following Filippenko (1982), eqs. 2 & 3 we also include
+		the	modification proposed by Barrell in 1951.

 		Parameters
 		---------- 
 		l 	: 	float
 			Wavelength in microns
+		T 	: 	float
+			Temperature in Kelvin
+		p 	: 	float
+			Atmospheric pressure in Pascal
+		RH 	: 	float
+			Relative humidity (0=<H=<1)
+		f 	: 	float
+			Water vapor pressure in Pascal

 		Returns
 		-------
@@ -83,14 +100,27 @@ def Edlen1953(l):
 	'''
 	s = 1 / l
 	R = 6432.8 + 2949810 / (146 - s**2) + 25540/(41-s**2)
-	n = R/1e8 + 1
+	n_standard = R/1e8 + 1
+
+	pHg = hPa2mmHg(p/100)
+	Tc  = K2C(T)
+	if f is None:
+		f = RH * svp(T)
+	f_mmHg = hPa2mmHg(f/100)
+
+	# Changes in temperature and pressure, from Barrell (1951)
+	n_tp = R/1e8 * (pHg * (1 + (1.049 - 0.0157*Tc) * 1e-6 * pHg ) ) / ( 720.883 * (1 + 0.003661 * Tc)) + 1
+
+	# Correction for water vapour pressure
+	n_svp = (n_tp - 1) * 1e6 - ( 0.0624 - 0.000680 * s**2 ) / (1 + 0.003661 * Tc) * f_mmHg
+	n = n_svp/1e6 + 1
 	return n


-def Edlen1966(l):
+def Edlen1966(l, T=288.15, p=101325, RH=0, CO2=300, f=None):
 	'''
 		From Edlen, 1966
-		Equation 1
+		Equations 1, 12, 17, 22

 		At standard conditions only

@@ -98,6 +128,16 @@ def Edlen1966(l):
 		---------- 
 		l 	: 	float
 			Wavelength in microns
+		T 	: 	float
+			Temperature in Kelvin
+		p 	: 	float
+			Atmospheric pressure in Pascal
+		RH 	: 	float
+			Relative humidity (0=<H=<1)
+		CO2 : 	float
+			CO2 density in parts per million
+		f 	: 	float
+			Water vapor pressure in Pascal

 		Returns
 		-------
@@ -107,15 +147,37 @@ def Edlen1966(l):
 	s = 1 / l
 	R = 8342.13 + 2406030 / (130 - s**2) + 15997 / (38.9 - s**2)
 	n = R/1e8 + 1
+
+	# Influence of temperature and pressure, eq. 12
+	Tc = K2C(T)
+	Ptorr = Pa2Torr(p)
+	R_tp = ( Ptorr * R/1e8 ) / 720.775 * (1 + Ptorr * (0.817 - 0.0133 * Tc) * 1e-6) / ( 1 + 0.0036610 * Tc)
+
+	# Allow for a CO2 fraction, eq. 17, but we include T and P dependence.
+	# This is okay from "Asthe refractivity of CO, is about 50% higher than
+	# that of air this CO,-content increases the refractivity over that of 
+	# CO2-free air by approximately 0.5 x 0.0003 x (n-1)s" and "As the 
+	# density factor of CO, is not much different from that of air this ratio r
+	# can be taken as independent of t and p ."
+	R_x = R_tp * (1 + 0.540 * (CO2/1e6 - 0.0003))
+	n_x = R_x + 1
+
+	# Correction for water vapour, eq. 22
+	if f is None:
+		f = RH * svp(T)
+	ftorr = Pa2Torr(f)
+
+	n = n_x - ftorr * (5.7224 - 0.0457 * s**2) * 1e-8
 	return n


 def Owens(l, T=288.15, p=101325, RH=0):
 	'''
 		From Owens, 1967
-		Equation 32
-		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
 		----------
@@ -133,18 +195,35 @@ def Owens(l, T=288.15, p=101325, RH=0):
 		n 	: 	float
 			Refractive index of atmospheric air at the given conditions
 	'''
-	P = hPA2mmHg(p/100)
-	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
 		---------- 
@@ -154,8 +233,12 @@ def BonschPotulski(l, T=288.15, p=101325, CO2=450):
 			Temperature in Kelvin
 		p 	: 	float
 			Atmospheric pressure in Pa
+		RH 	: 	float
+			Relative humidity (0<=H<=1)
 		CO2 : 	float
 			Particles per million of CO2
+		f 	: 	float
+			Water vapor pressure in Pascal

 		Returns
 		-------
@@ -163,21 +246,26 @@ def BonschPotulski(l, T=288.15, p=101325, CO2=450):
 			Refractive index of atmospheric air at the given conditions
 	'''
 	x = CO2 / 1e6
-	T = K2C(T)
+	Tc = K2C(T) 
 	s = 1 / l
-	# Rs = 8092.33 + 2333983 / (130 - s**2) + 15518 / (38.9 - s**2)
+
 	Rs = 8091.37 + 2333983 / (130 - s**2) + 15518 / (38.9 - s**2)
 	Rs = Rs / 1e8
 	Rx = Rs * (1 + 0.5327 * (x - 0.0004))
-	Rtp = Rx * p / 93214.60 * (1 + 1e-8 * (0.5953 - 0.009876 * T)*p) / (1 + 0.0036610 * T)
-	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

@@ -185,15 +273,32 @@ def BirchDowns(l):
 		----------
 		l 	: 	float
 			wavelength in micron
+		T 	: 	float
+			Temperature in Kelvin
+		p 	: 	float
+			Atmospheric pressure in Pa
+		RH 	: 	float
+			Relative humidity (0<=H<=1)
+		f 	: 	float
+			Water vapor pressure in Pascal

 		Returns
 		-------
 		n 	: 	float
 			Refractive index of atmospheric air at the given conditions
 	'''
+	Tc = K2C(T)
 	s = 1 / l
-	R = 8342.54 + 2406147 / (130 - s**2) + 15998 / (38.9-s**2)
-	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


@@ -322,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
@@ -368,7 +473,7 @@ def Mathar(l, T=288.15, p=101325, RH=0):
 	if l < 1.3 or l > 24:
 		print('Using the Mathar model at this wavelength (l={:.3f}) is not recommended!'.format(l))

-	if l <= 2.65:
+	elif l <= 2.65:
 		########################################################################
 		# model parameters for wavelengths between 1.3 and 2.5 microns
 		cref = [ 0.200192e-3,  0.113474e-9,  -0.424595e-14,  0.100957e-16, -0.293315e-20,  0.307228e-24] # cm^j
@@ -388,7 +493,7 @@ def Mathar(l, T=288.15, p=101325, RH=0):
 		Href = 10          #%, reference humdity in percent
 		########################################################################
 	
-	if l > 2.65 and l <= 4.275:
+	elif l > 2.65 and l <= 4.275:
 		########################################################################
 		# model parameters for wavelengths between 2.8 and 4.2 microns
 		cref = [ 0.200049e-3,  0.145221e-9,   0.250951e-12, -0.745834e-15, -0.161432e-17,  0.352780e-20] # cm^j
@@ -408,7 +513,7 @@ def Mathar(l, T=288.15, p=101325, RH=0):
 		Href = 10          #%, reference humdity in percent
 		########################################################################

-	if l > 4.275 and l <= 6.35:
+	elif l > 4.275 and l <= 6.35:
 		########################################################################
 		# model parameters for wavelengths between 4.35 and 5.2 microns
 		cref = [ 0.200020e-3,  0.275346e-9,   0.325702e-12, -0.693603e-14,  0.285610e-17,  0.338758e-18] # cm^j
@@ -427,7 +532,7 @@ def Mathar(l, T=288.15, p=101325, RH=0):
 		Href = 10          #%, reference humdity in percent
 		########################################################################

-	if l > 6.35 and l <= 15.05:
+	elif l > 6.35 and l <= 15.05:
 		########################################################################
 		# model parameters for wavelengths between 7.5 and 14.1 microns
 		cref = [ 0.199885e-3,  0.344739e-9,  -0.273714e-12,  0.393383e-15, -0.569488e-17,  0.164556e-19] # cm^j
@@ -447,7 +552,7 @@ def Mathar(l, T=288.15, p=101325, RH=0):
 		########################################################################


-	if l > 15.05:
+	elif l > 15.05:
 		########################################################################
 		# model parameters for wavelengths between 16 and 24 microns
 		cref = [ 0.199436e-3,  0.299123e-8,  -0.214862e-10,  0.143338e-12,  0.122398e-14, -0.114628e-16] # cm^j
@@ -464,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
@@ -481,6 +586,7 @@ def Mathar(l, T=288.15, p=101325, RH=0):
 	return n 


+# Other necessary conversion or calculations
 def K2C(K):
 	'''
 		Kelvin to Celsius conversion
@@ -498,7 +604,7 @@ def K2C(K):
 	return K - 273.15


-def hPA2mmHg(hPa):
+def hPa2mmHg(hPa):
 	'''
 		hPa to mmHg conversion

@@ -514,3 +620,41 @@ def hPA2mmHg(hPa):
 	'''
 	return 0.7500616827 * hPa

+
+def Pa2Torr(Pa):
+	'''
+		Pa to Torr conversion
+
+		Parameters
+		----------
+		Pa 	:	float
+			Atmospheric pressure in Pa (Pascal)
+
+		Returns
+		-------
+		Torr 	: 	float
+			Atmospheric pressure in Torr
+	'''
+
+	return Pa * 760/101325
+
+
+def svp(T):
+	'''
+		Saturation vapour pressure in Pa.
+
+		Parameters
+		----------
+		T 	: float
+			Atmospheric temperature in Kelvin
+
+		Returns
+		-------
+		svp	: float
+			Saturation vapour pressure in Pa
+	'''
+	_A =  1.2378847e-5 # K^-2
+	_B = -1.9121316e-2 # K^-2
+	_C = 33.93711047
+	_D = -6.3431645e3  # K^-2
+	return np.exp(_A*T**2 + _B*T + _C + _D/T)
@@ -0,0 +1 @@
+include AstroAtmosphere/data/1976USSA.txt
@@ -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
- `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).


+
 ## The Ciddor refractivity model

 The refractivity model from Ciddor, 1996, calculates the index of refraction of air using the following parameters:
@@ -107,7 +110,7 @@ Input parameters for the spherical atmosphere models (all but the first):
 - xc - CO2 density              [ppm]
 - lat - Latitude                [deg]
 - h - Height above sea level    [m]
- l - Wavelength(s)             [um]
+- l - Wavelength(s)             [μm]
 - z - Zenith angle              [deg]

 An example - To calculate the atmospheric refraction we set up a `refraction()` object. This requires the latitude and altitude of the observer. The atmospheric dispersion is calculated similarly, but by exchanging the `refraction()` object by a `dispersion()` object and parsing two refractive indices into the `cassini()` function.
@@ -118,7 +121,7 @@ rho = at.rho(p=101325, T=273.15, RH=0.2, xc=300)

 # Initializing refraction model and setting the reduced height
 ref = refraction(lat=30, h=2200)
-ref.setReducedHeight(p=101325, rho)
+ref.setReducedHeight(p=101325, rho=rho)

 # Calculating the atmospheric refraction
 atm_ref = ref.cassini(n, zenith=30)
@@ -186,12 +189,13 @@ Included are
 | ---------- | -------- | ------- |
 | Barrell & Sears (1939) | `BarrellAndSears()` | Outdated model, especially poor for short wavelengths. |
 | Ciddor (1996) | `Ciddor()` | Only at standard atmosphere. CO2 density can be changed. |
-| Edlen (1953) | `Edlen1953()` | Only at standard atmosphere. Function requires wavelength only. |
-| Edlen (1966) | `Edlen1966()` | Only at standard atmosphere. Function requires wavelength only. |
-| Owens (1967) | `Owens()` | - |
-| Bonsch & Potulski (1998) | `BonschPotulski()` | Function does not include humidity |
+| 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()` | 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.|

@@ -199,9 +203,8 @@ The last function extracts the refractive index from the `pyslalib.sla_refro()`

 ## Quick functions

-Some functions that assume default conditions for the atmosphere. Then only the wavelength and zenith angle are required to calculate the refraction and dispersion. 
+Some functions can make use of the `observatories` dictionary as input. These include default average conditions of the atmosphere at these locations. When these are used, only the wavelength and zenith angle are required to calculate the refraction and dispersion. 

-Presently, only the `'STANDARD'` and `'CERRO_ARMAZONES'` flags are known. Contact the author to add additional standard conditions.

 ```python
 # Quick refractive index
@@ -214,6 +217,20 @@ rq = quick_refraction(l1, zenith_obs, conditions='STANDARD')
 dq = quick_dispersion(l1, l2, zenith_obs, conditions='STANDARD')
 ```

+Presently, several observatories are known. These can all be found in the `observatories` dictionary. Contact the author to add additional standard conditions.
+
+| Observatory                                   | Flag              | Notes                                         | 
+| --------------------------------------------- | ----------------- | --------------------------------------------- |
+| Standard conditions                           | `STANDARD`        |                                               |
+| Cerro Armazones                               | `CERRO_ARMAZONES` | Location of the future ELT                    |
+| Cerro Paranal                                 | `CERRO_PARANAL`   | Location of the VLT, VISTA, VST, NGTS and SSO |
+| Roque de los Muchachos Observatory (La Palma) | `LA_PALMA`        | GTC, INT, WHT, NOT, and others.               |
+| La Silla Observatory                          | `LA_SILLA`        | ESO 3.6m telescope, NTG, BlackGEM and others. | 
+| Las Campanas Observatory                      | `LAS_CAMPANAS`    | Site of the future GMT.                       |
+| Mauna Kea Observatory (Hawaii)                | `MAUNA_KEA`       | Location of Keck, Subaru, Gemini North and others.|
+
 ## Other

-This package also contains a python native port of the `SLA_REFRO()` routine of the pyslalib package, which was originally written in FORTRAN by P.T. Wallace. It can be called directly using `sla_refro()`. When using the same units as the other refraction functions, you might prefer using `slalib_refraction()` or `slalib_dispersion()`.
+This package also contains a python native port of the `SLA_REFRO()` routine of the pyslalib package, which was originally written in FORTRAN by P.T. Wallace. It can be called directly using `sla_refro()`. When using the same units as the other refraction functions, you might prefer using `slalib_refraction()` or `slalib_dispersion()`. Do note that this implementation is slower than the (compiled) FORTRAN implementation.
+
+Finally, the dispersion model used by Alexei Filippenko in 1982 is also included. This model is a combination of the (modified) Edlèn equation and a plane-parallel atmospheric geometry. This calculation can be accessed by calling the `Filippenko1982()` function.
@@ -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__':
@@ -21,14 +21,14 @@ def setup_package():
    # write_version_py()

    setup(name='AstroAtmosphere',
-          version='1.2',
+          version='1.5',
          description='Collection of equations for refractive index, refraction and dispersion calculations of atmospheric air.',
          long_description=__readme__,
          long_description_content_type='text/markdown',
          author='Joost van den Born',
          author_email='born@astron.nl',
          url='https://gitlab.astro-wise.org/micado/atmosphericmodel',
-          download_url='https://gitlab.astro-wise.org/micado/atmosphericmodel/-/archive/v1.2/atmosphericmodel-v1.2.tar.gz',
+          download_url='https://gitlab.astro-wise.org/micado/atmosphericmodel/-/archive/v1.5/atmosphericmodel-v1.5.tar.gz',
          license='GNU General Public License v3',
          data_files = [('',['LICENSE.txt'])],
          package_dir={'AstroAtmosphere': 'AstroAtmosphere'},
Author	SHA1	Message	Date
yosata	ee7134da7d	Fixed Owens (1967) model	2022-07-06 12:19:31 +02:00
Joost van den Born	27abba5f5d	Update README.md	2022-03-23 08:55:21 +00:00
yosata	03d583a9cd	Improved quick functions	2022-03-23 08:42:45 +01:00
yosata	3cdf4ccd39	Added TPH dependence to BirchDowns()	2021-03-10 12:27:07 +01:00
yosata	2652d75eec	Update README	2021-03-10 12:04:27 +01:00
yosata	f4e77ffaf7	Ready for v1.4 release.	2021-03-10 11:52:14 +01:00
yosata	7771c59268	Added TPH dependency to Edlen and BonschPotulski functions.	2021-03-10 11:50:00 +01:00
yosata	160818cdf7	Merge of local branch with master branch - readying for 1.3.1 release.	2020-11-15 18:05:37 +01:00
yosata	a3521a2426	Small bugfixes in refractionModels.py	2020-11-15 18:04:39 +01:00
Joost van den Born	f4b37bd9bf	Update README.md	2020-09-08 10:38:24 +02:00
Joost van den Born	b804b15405	Update README.md	2020-09-08 10:33:28 +02:00
yosata	7a084d1635	Upated setup.py for v1.3 release	2020-09-07 21:29:56 +02:00
Joost van den Born	f1aa2ee903	Update README.md	2020-09-07 21:24:43 +02:00
yosata	8f9a0bbbea	Added Filippenko model, updated B&S refractivity	2020-09-07 21:20:09 +02:00
yosata	cce60766fc	included data files into installation	2020-08-07 13:32:24 +02:00
yosata	2b1835c814	fixed minor bug	2020-08-05 15:59:25 +02:00