Upated setup.py for v1.3 release

AstroAtmosphere/__init__.py

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

AstroAtmosphere/misc.py

@@ -3,12 +3,93 @@ 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
+import numpy as np
 
 '''
 @ 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, 
@@ -23,8 +104,10 @@ def quick_refractive_index(l, conditions='STANDARD', T=288.15, p=101325,
 		l 	: float
 			Wavelength in microns
 		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,23 +126,8 @@ 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()
@@ -84,8 +152,10 @@ def quick_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)
@@ -104,23 +174,8 @@ 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()
@@ -154,8 +209,10 @@ def quick_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15,
 				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,23 +231,8 @@ 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()
@@ -228,8 +270,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 +289,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 +301,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 +355,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 +374,11 @@ def slalib_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15, p=101325,
 		tlr :	float
 			Temperature lapse rate (default=0.0065)
 		eps :	float
-			Precision at which the calculation is aborted (default=1e-10)
+			Precision at which the calculation is aborted (default=1e-8)
 		use_pyslalib : bool
 			Flag that indicates that the user wants to use the pySlalib
 			package. If False, the internal python port of the original
-			FORTRAN code will be used. (Default=False).
+			FORTRAN code will be used. (default=False).
 			
 		Returns
 		-------

AstroAtmosphere/observatories.py

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

AstroAtmosphere/refractivityModels.py

@@ -9,7 +9,7 @@ import numpy as np
 				atmospheric dispersion, given a set of conditions.
 '''
 
-def BarrellAndSears(l, T=288.15, p=101325, RH=0):
+def BarrellAndSears(l, T=288.15, p=101325, RH=0, f=None):
 	'''
 		From Barrell & Sears, 1939
 		Equation 7.7, page 52
@@ -26,15 +26,20 @@ def BarrellAndSears(l, T=288.15, p=101325, RH=0):
 			Atmospheric pressure in Pascal
 		RH 	: 	float
 			Relative humidity (0=<H=<1)
+		f 	: 	float
+			Water vapour pressure in Pascal
 
 		Returns
 		-------
 		n 	: 	float
 			Refractive index of atmospheric air at the given conditions
 	'''
+	if f is None:
+		f = RH * svp(T)
+	f_mmHg = hPA2mmHg(f/100)
 	T = K2C(T)
 	P = hPA2mmHg(p/100)
-	R = (0.378125 + 0.0021414 / l**2 + 0.00001793 / l**4) * (P*(1 + (1.049 - 0.0157*T) * P * 1e-6)/(1 + 0.003661*T)) - (0.0624 - 0.000680 / l**2) * RH / (1 + 0.003661*T)
+	R = (0.378125 + 0.0021414 / l**2 + 0.00001793 / l**4) * (P*(1 + (1.049 - 0.0157*T) * P * 1e-6)/(1 + 0.003661*T)) - (0.0624 - 0.000680 / l**2) * f_mmHg / (1 + 0.003661*T)
 	n = R/1e6 + 1
 	return n
 
@@ -64,17 +69,25 @@ 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
+		f 	: 	float
+			Water vapor pressure in Pascal
 
 		Returns
 		-------
@@ -83,7 +96,20 @@ 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
 
 
@@ -481,6 +507,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
@@ -514,3 +541,23 @@ def hPA2mmHg(hPa):
 	'''
 	return 0.7500616827 * hPa
 
+
+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

@@ -63,6 +63,8 @@ The package contains the following files:
 - `misc.py` 				: Contains several miscellaneous functions, e.g. to quickly calculate the refractive index.
 - `refractionModels.py`		 		: Contains various atmospheric models that can be used to calculate the atmospheric refraction
 - `refractivityModels.py` 		: Contains various functions from literature to calculate the refractive index. Some examples of works included are Ciddor (1996), Barrell & Sears (1939), SLALIB, Owens (1967) and Edlen (1966).
+- `observatories` 				: Contains a dictionary with the average atmospheric conditions for various large observatories around the world.
+
 
 
 ## The Ciddor refractivity model
@@ -186,7 +188,7 @@ 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. |
 | Owens (1967) | `Owens()` | - |
 | Bonsch & Potulski (1998) | `BonschPotulski()` | Function does not include humidity |
@@ -199,9 +201,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 +215,19 @@ 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.                       |
+
 ## 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()`.
+
+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.

setup.py

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