Updated documentation, descriptions for v1.0 release.

AstroAtmosphere/ciddorModel.py

@@ -16,6 +16,53 @@ from scipy import integrate
 '''
 
 class Observatory:
+	'''
+		Class used for all function to calculate the refractive index
+		and the corresponding uncertainties.
+
+		...
+
+		Attributes
+		----------
+		dT 	: 	float (optional)
+			Uncertainty in the atmospheric temperature in Kelvin
+		dP 	: 	float (optional)
+			Uncertainty in the atmospheric pressure in Pascal
+		dH 	: 	float (optional)
+			Uncertainty in the atmospheric relative humidity
+		dCO2: 	float (optional)
+			Uncertianty in the atmospheric CO2 density in parts
+			per million.
+		dl  : 	float (optional)
+			Uncertainty in the wavelengths in micron
+
+		Methods
+		-------
+		n_as(l)
+			Returns the refractive index for standard atmospheric conditions.
+		n_axs(l, xc)
+			Same as n_as(l), but with a different CO2 density given by xc
+		n_ws(l)
+			Returns the refractive index of air for standard moist contions.
+		f(p,T)
+			Returns the enchancement factor of water vapour in air.
+		svp(T)
+			Returns the saturation vapour pressure in Pa.
+		svp2(T)
+			Returns the saturation vapour pressure in Pa. 
+			From IAPWS formula, based on documentation by NIST toolbox.
+		xw(p, T, RH)
+			Returns the molar fraction of water vapour in moist air.
+		Z(p, T, RH, xw='')
+			Returns the compressibility of moist air
+		rho(p, T, RH, xc, xw='')
+			Returns the density of atmospheric air at the given parameters.
+		n_tph(l, T, p, RH, xc)
+			Returns the refractive index for any atmospheric conditions.
+		dn_tph(l, T, p, RH, xc, dl=None, dT=None, dP=None, dRH=None, dCO2=None)
+			Returns the uncertainty in the refractive index for any atmospheric
+			conditions, given the parameters.
+	'''
 	def __init__(self, dT=0.2, dP=20, dH=0.02, dCO2=20, dl=0.001):
 		self.dT 	= dT 	# in K
 		self.dP 	= dP 	# in Pa
@@ -28,6 +75,17 @@ class Observatory:
 			Formula that calculate the refractive index for standard atmospheric
 			conditions, i.e. 15 C, 101,325 Pa, 0% RH, 450 ppm CO2.
 			It returns n_as, not (n_as - 1) * 1e8!
+
+			Parameters
+			----------
+			l 	: float
+				Wavelength of light in microns
+
+			Returns
+			-------
+			n 	: float
+				Refractive index of air at standard atmospheric conditions.
+
 		'''
 		_s = 1 / l
 		_k0 = 238.0185 	# um^-2
@@ -42,6 +100,18 @@ class Observatory:
 			Formula that calculate the refractive index for standard atmospheric
 			conditions, i.e. 15 C, 101,325 Pa, 0% RH, but with xc ppm of CO2.
 			It returns n_axs, not (n_axs - 1)
+
+			Parameters
+			----------
+			l 	: float
+				Wavelength of light in microns
+			xc 	: float
+				CO2 density in parts per million
+
+			Returns
+			-------
+			n 	: float
+				Refractive index of air
 		'''
 		_n_as = self.n_as(l)
 		return ( (_n_as - 1) * (1 + 0.534e-6 * (xc - 450) ) ) + 1
@@ -52,6 +122,16 @@ class Observatory:
 			Formula that calculates the refractivity of air for standard moist conditions
 			as defined by B&S (1939), i.e. 20 C, 1333 Pa
 			Returns n_ws, not (n_ws - 1) * 1e8
+
+			Parameters
+			----------
+			l 	: float
+				Wavelength of light in microns
+
+			Returns
+			-------
+			n 	: float
+				Refractive index of standard moist air.
 		'''
 		_s  = 1 / l
 		_cf = 1.022
@@ -65,8 +145,18 @@ class Observatory:
 	def f(self, p, T):
 		'''
 			Enhancement factor of water vapour in air
-			t in degrees celsius
-			p in Pascal (Pa)
+
+			Parameters
+			----------
+			p 	: float
+				Atmospheric pressure in Pascal
+			T 	: float
+				Atmospheric temperature in Kelvin
+
+			Returns
+			-------
+			f 	: float
+				Enhancement factor of water vapour in air
 		'''
 		_t = T - 273.15
 		_a = 1.00062
@@ -78,6 +168,16 @@ class Observatory:
 	def svp(self, T):
 		'''
 			Saturation vapour pressure in Pa.
+
+			Parameters
+			----------
+			T 	: float
+				Atmospheric temperature in Kelvin
+
+			Returns
+			-------
+			svp	: float
+				Saturation vapour pressure
 		'''
 		_A =  1.2378847e-5 # K^-2
 		_B = -1.9121316e-2 # K^-2
@@ -88,7 +188,18 @@ class Observatory:
 
 	def svp2(self, T):
 		'''
-			Calculation of svp from IAPWS formula, based on documentation by NIST toolbox.
+			Calculation of svp from IAPWS formula, based on documentation 
+			by NIST toolbox.
+
+			Parameters
+			----------
+			T 	: float
+				Atmospheric temperature in Kelvin
+
+			Returns
+			-------
+			svp	: float
+				Saturation vapour pressure
 		'''
 		K1 =  1.16705214528e+03
 		K2 = -7.24213167032e+05
@@ -111,6 +222,20 @@ class Observatory:
 	def xw(self, p, T, RH):
 		'''
 			Molar fraction of water vapour in moist air.
+
+			Parameters
+			----------
+			p 	: float
+				Atmospheric pressure in Pascal
+			T 	: float
+				Atmospheric temperature in Kelvin
+			RH 	: float
+				Atmospheric relative humidity (0<=RH<=1)
+
+			Returns
+			-------
+			xw	: float
+				Molar fraction of water vapour in moist air
 		'''
 		_f 		= self.f(p, T)
 		_svp 	= self.svp(T)
@@ -120,6 +245,20 @@ class Observatory:
 	def Z(self, p, T, RH, xw=''):
 		'''
 			Compressibility of moist air
+
+			Parameters
+			----------
+			p 	: float
+				Atmospheric pressure in Pascal
+			T 	: float
+				Atmospheric temperature in Kelvin
+			RH 	: float
+				Atmospheric relative humidity (0<=RH<=1)
+
+			Returns
+			-------
+			Z 	: float
+				Compressibility of moist air
 		'''
 		_a0 =  1.58123e-6 # K Pa^-1
 		_a1 = -2.9331e-8  # Pa^-1
@@ -142,10 +281,25 @@ class Observatory:
 	def rho(self, p, T, RH, xc, xw=''):
 		'''
 			Calculates the density of air for the given parameters.
-			T in Kelvin
-			p in Pascal
-			RH between 0 and 1
-			xc in ppm
+			
+			Parameters
+			----------
+			p 	: float
+				Atmospheric pressure in Pascal
+			T 	: float
+				Atmospheric temperature in Kelvin
+			RH 	: float
+				Atmospheric relative humidity (0<=RH<=1)
+			xc 	: float
+				CO2 density in parts per million
+			xw 	: str or float
+				Placeholder for molar fraction of water vapour in moist air.
+				Used to pass the function along where it should go.
+
+			Returns
+			-------
+			rho : float
+				Density of air
 		'''
 		_Ma = 1e-3 * (28.9635 + 12.011e-6 * (xc - 400) )	# Molar mass in kg/mol of dry air containing xc ppm of CO2
 		_Mw = 0.018015 										# kg/mol
@@ -161,6 +315,25 @@ class Observatory:
 	def n_tph(self, l, T, p, RH, xc):
 		'''
 			Calculates the refraction of the atmosphere for any conditions.
+			Follows the equations by Ciddor (1996).
+
+			Parameters
+			----------
+			l 	: float
+				Wavelength of light in microns
+			p 	: float
+				Atmospheric pressure in Pascal
+			T 	: float
+				Atmospheric temperature in Kelvin
+			RH 	: float
+				Atmospheric relative humidity (0<=RH<=1)
+			xc 	: float
+				CO2 density in parts per million
+
+			Returns
+			-------
+			n 	: float
+				Refractive index of air
 		'''
 		_n_a 	 = self.n_axs(l, xc)
 		_n_w 	 = self.n_ws(l)
@@ -194,17 +367,33 @@ class Observatory:
 			This involves almost endless use of the chain rule and evaluation of both the moist
 			air and dry air components. Note: dry air does not contain any water, i.e. H=0, x_w=0, etc.
 
-			Parameters:
-			l 	- Wavelength in micron (0.4 - 2.5 micron)
-			T 	- Temperature in Kelvin
-			p 	- Pressure in Pascals
-			RH 	- Relative humidity (no units, don't use percentage but value 0<x<1)
-			xc 	- CO2 density in ppm
-			dl - Uncertainty in the blue wavelength
-			dT 	- Uncertainty in the temperature, default: 0.2 K
-			dP 	- Uncertainty in the pressure, default: 20 Pa
-			dRH - Uncertainty in the relative humidity, default: 0.02
-			dCO2- Uncertainty in the CO2 density, default: 20
+			Parameters
+			----------
+			l 	: float
+				Wavelength of light in microns
+			p 	: float
+				Atmospheric pressure in Pascal
+			T 	: float
+				Atmospheric temperature in Kelvin
+			RH 	: float
+				Atmospheric relative humidity (0<=RH<=1)
+			xc 	: float
+				CO2 density in parts per million
+			dl 	: float 
+				Uncertainty in the wavelength in microns
+			dT 	: float (default: 0.2)
+				Uncertainty in the atmospheric temperature in Kelvin
+			dP 	: float (default: 20)
+				Uncertainty in the atmospheric pressure in Pascal
+			dRH : float (default: 0.02)
+				Uncertainty in the atmospheric relative humidity
+			dCO2: float (default: 20)
+				Uncertainty in the atmospheric CO2 density in parts per million
+
+			Returns
+			-------
+			dn 	: float
+				Uncertainty in the refractive index of air
 		'''
 
 		if dP == None:

AstroAtmosphere/dispersionModels.py

@@ -3,6 +3,7 @@ import numpy as np
 from scipy import integrate
 import os
 from .ciddorModel import Observatory
+
 '''
 @ Author:		Joost van den Born 
 @ Contact:		born@astron.nl
@@ -17,19 +18,99 @@ from .ciddorModel import Observatory
 '''
 
 class dispersion(Observatory):
+	'''
+		Class that sets the observer location.
+		These location attributes are then used to calculate the 
+		appropriate atmospheric dispersion.
+
+		See also:
+		Observatory()
+		refraction()
+		...
+
+		Attributes
+		-----------
+		lat : 	float
+			Latitude of the observer in degrees
+		h 	: 	float
+			Altitude above mean sea level of the observer in meters
+		Other Observatory() attributes
+
+		Methods
+		-------
+		setReducedHeight(p, rho)
+			Sets the reduced atmospheric height.
+		planeParallel(n1, n2, zenith)
+			Returns the atmospheric dispersion assuming a plane parallel
+			atmosphere
+		cassini(n1, n2, zenith, reduced_height=None)
+			Returns the atmospheric dispersion asusming Cassini's homegeneous
+			atmosphere model
+		tan5(n1, n2, zenith, reduced_height=None)
+			Returns the atmospheric dispersion following Oriani's theorem,
+			including a fifth order term in the expansion.
+		oriani(n1, n2, zenith, reduced_height=None)
+			Returns the atmospheric dispersion following the classic version
+			of Oriani's theorem.
+		corbard(n1, n2, zenith, reduced_height=None)
+			Returns the atmospheric dispersion following the method described 
+			in Corbard et al., 2019.
+		refractionIntegral(n1, n2, z0, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None)
+			Returns the atmospheric dispersion following the refraction integral.
+		psi(x)
+			Submethod for corbard()
+		set_rc()
+			Sets the local radius of curvature of the Earth at the observer.
+		set_H(p, rho)
+			Sets the reduced height for the conditions at the observer, assuming
+			an adiabatic atmosphere.
+		H_isotherm(T)
+			Sets the reduced height for the conditions at the observer, assuming
+			an isothermal atmosphere.
+		cassiniError(z, l1, l2, T, p, RH, xc, dl1=0, dl2=0, dT=0.2, dP=20, 
+					dRH=0.02, dCO2=20, dz=0, lat=None, h=None)
+			Returns the uncertainty in the atmospheric dispersion, assuming 
+			Cassini's homegeneous atmosphere model and the given parameters.		
+		Other Observatory() methods
+	'''
 	def __init__(self, lat, h):
 		Observatory.__init__(self)
 		self.h  	= h 						# height above sea level in meters
 		self.lat 	= lat						# Latitude of the observer in degrees
 		self.rc 	= self.set_rc() * 1000		# Radius of curvature of the earth at the observer
 
-	def setReducedHeight(self, P, rho):
-		self.H  	= self.set_H(P, rho) 		# Reduced height of the atmosphere assuming ideal gas law
+	def setReducedHeight(self, p, rho):
+		'''
+			Sets the reduced atmospheric height.
+
+			Parameters
+			----------
+			p 	: float
+				Atmospheric pressure in Pascal
+			rho : float
+				Density of the atmospheric air at the observer
+		'''
+		self.H  	= self.set_H(p, rho) 		# Reduced height of the atmosphere assuming ideal gas law
 
 
 	def planeParallel(self, n1, n2, zenith):
 		'''
 			Refraction based on a flat atmosphere, easily derived from Snell's law.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the shortest wavelength of interest
+			n2 : float
+				Refractive index at the longest wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+
+			Returns
+			-------
+			dR : float
+				Atmospheric dispersion in degrees
 		'''
 		_R1 = np.arcsin(n1 * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
 		_R2 = np.arcsin(n2 * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
@@ -40,6 +121,24 @@ class dispersion(Observatory):
 	def cassini(self, n1, n2, zenith, reduced_height=None):
 		'''
 			Refraction of spherical atmosphere, derived from geometry using Sine law
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the shortest wavelength of interest
+			n2 : float
+				Refractive index at the longest wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			dR : float
+				Atmospheric dispersion in degrees
 		'''
 		if reduced_height:
 			H = reduced_height
@@ -56,6 +155,24 @@ class dispersion(Observatory):
 		'''
 			Same as oriani's theorem, but including a higher order term.
 			Found in Corbard et al., 2019.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the shortest wavelength of interest
+			n2 : float
+				Refractive index at the longest wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			dR : float
+				Atmospheric dispersion in degrees
 		'''
 		if reduced_height:
 			H = reduced_height
@@ -74,8 +191,27 @@ class dispersion(Observatory):
 
 	def oriani(self, n1, n2, zenith, reduced_height=None):
 		'''
-			Classic refraction formula with two tan(z) terms. This does not assume any 
-			information about the structure of the atmosphere. Found with Oriani's theorem.
+			Classic refraction formula with two tan(z) terms. 
+			This does not assume any information about the structure
+			of the atmosphere. Found with Oriani's theorem.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the shortest wavelength of interest
+			n2 : float
+				Refractive index at the longest wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			dR : float
+				Atmospheric dispersion in degrees
 		'''
 		if reduced_height:
 			H = reduced_height
@@ -94,9 +230,28 @@ class dispersion(Observatory):
 
 	def corbard(self, n1, n2, zenith, reduced_height=None):
 		'''
-			Corbard et al., 2019, mention an additional formula based on the error function.
-			Oriani's theorem can be derived from this equation by 'keeping only the three first
-			terms of its asymptotic expansion'.
+			Corbard et al., 2019, mention an additional formula based on
+			the error function. Oriani's theorem can be derived from this
+			equation by 'keeping only the three first terms of its 
+			asymptotic expansion'.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the shortest wavelength of interest
+			n2 : float
+				Refractive index at the longest wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			dR : float
+				Atmospheric dispersion in degrees
 		'''
 		if reduced_height:
 			H = reduced_height
@@ -113,15 +268,46 @@ class dispersion(Observatory):
 
 
 
-	def refractionIntegral(self, n1, n2, z0, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None):
+	def refractionIntegral(self, n1, n2, z0, R0=None, useStandardAtmosphere=True, 
+				heightData=None, rhoData=None):
 		'''
-			Using data that gives temperature, pressure and density as a function of height, 
-			the gladstone-dale relation is invoked to calculate the path of a light ray 
-			along the atmosphere. This is allows us to also include atmospheric activity in 
+			Using data that gives temperature, pressure and density as a 
+			function of height, the gladstone-dale relation is invoked to
+			calculate the path of a light ray along the atmosphere. 
+			This is allows us to also include atmospheric activity in 
 			different atmosphere layers.
 			Recommended use with the US Standard Atmosphere from 1976.
 			USSA1976 data generated from code supplied by http://www.pdas.com/atmos.html
-			The method described in eq. 3 in Auer & Standish (2000) is used for the calculation.
+			The method described in eq. 3 in Auer & Standish (2000) is used for the 
+			calculation.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the shortest wavelength of interest
+			n2 : float
+				Refractive index at the longest wavelength of interest
+			z0 : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			R0 : float (optional)
+				Local radius of curvature of the Earth in km.
+			useStandardAtmosphere : bool
+				Flag to use the included USSA1976 data
+			heightData : 1D Array
+				If "useStandardAtmosphere" is False, the function will ask
+				for alternative atmospheric data to integrate over. This
+				gives the altitudinal points
+			rhoData : 1D Array
+				If useStandardAtmosphere is False, the function will ask
+				for alternative atmospheric data to integrate over. This
+				gives the atmospheric density at the points defined by 
+				"heightData".
+
+			Returns
+			-------
+			dR : float
+				Atmospheric dispersion in degrees
 		'''
 		if not R0:
 			R0 = self.rc
@@ -163,6 +349,19 @@ class dispersion(Observatory):
 
 
 	def psi(self, x):
+		'''
+			Submethod for corbard(). Described explicitly in Corbard et al.,2019.
+
+			Parameters
+			----------
+			x 	: 	float
+				Input of the function
+
+			Returns
+			-------
+			y 	: 	float
+				Output of the function
+		'''
 		_f   = lambda a: np.exp(-a**2)
 		_int = integrate.quad(_f, x, np.inf)
 		return np.exp(x**2) * _int[0]
@@ -173,6 +372,11 @@ class dispersion(Observatory):
 			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
@@ -182,9 +386,21 @@ class dispersion(Observatory):
 		return _rc0 + self.h/1000
 
 
-	def set_H(self, P, rho):
+	def set_H(self, p, rho):
 		'''
 			Calculates the reduced height of the atmosphere, assuming ideal gas law.
+
+			Parameters
+			----------
+			p 	: float
+				Atmospheric pressure in Pascal
+			rho : float
+				Density of the atmospheric air at the observer
+
+			Returns
+			-------
+			H 	: float
+				Reduced height of the atmosphere in meters
 		'''
 		_phi = np.deg2rad(self.lat)
 		_c1 = 5.2790414e-3
@@ -194,11 +410,22 @@ class dispersion(Observatory):
 		_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) * self.h + 7.2e-13 * self.h**2
-		return P / (rho * _g)
+		return p / (rho * _g)
 
 	def H_isotherm(self, T):
 		'''
-			Calculates the reduced height of the atmosphere, assuming an isotherm distribution.
+			Calculates the reduced height of the atmosphere, 
+			assuming an isotherm distribution.
+
+			Parameters
+			----------
+			T 	: float
+				Atmospheric temperature in Kelvin
+
+			Returns
+			-------
+			H 	: float
+				Reduced height of the atmosphere in meters
 		'''
 		_kb = 1.38064852e-23 # kg m2 s-2 K-1
 		_m  = 4.809651698e-26 # kg (weight per molecule of air, assuming dry air)
@@ -213,32 +440,63 @@ class dispersion(Observatory):
 
 		return (_kb*T) / (_m * _g0_local)
 
-	def cassiniError(self, z, l1, l2, T, p, RH, xc, dl1=0, dl2=0, dT=0.2, dP=20, dRH=0.02, dCO2=20, dz=0, lat=None, h=None):
+	def cassiniError(self, z, l1, l2, T, p, RH, xc, dl1=0, dl2=0, dT=0.2, dP=20,
+						dRH=0.02, dCO2=20, dz=0, lat=None, h=None):
 		'''
-			Calculates the uncertainty in the atmospheric dispersion, following the Cassini model,
-			by using fully propagated error analysis.
-			This involves almost endless use of the chain rule and evaluation of both the moist
-			air and dry air components. Note: dry air does not contain any water, i.e. H=0, x_w=0, etc.
+			Calculates the uncertainty in the atmospheric dispersion, 
+			following the Cassini model, by using fully propagated error analysis.
+			This involves almost endless use of the chain rule and evaluation 
+			of both the moist air and dry air components. 
+			Note: dry air does not contain any water, i.e. H=0, x_w=0, etc.
 
-			Not taken into account is the errror propagation in the scale height of the atmosphere or 
-			the zenith error, from my results this should give a small or even negligible difference.
+			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:
-			z 	- Zenith angle in degrees
-			l1 	- Bluest wavelength in micron (0.4 - 2.5 micron)
-			l2 	- Reddest wavelength in micron (0.4 - 2.5 micron)
-			T 	- Temperature in Kelvin
-			p 	- Pressure in Pascals
-			RH 	- Relative humidity (no units, don't use percentage but value 0<x<1)
-			xc 	- CO2 density in ppm
-			dl1 - Uncertainty in the blue wavelength
-			dl2 - Uncertainty in the red wavelength
-			dT 	- Uncertainty in the temperature, default: 0.2 K
-			dP 	- Uncertainty in the pressure, default: 20 Pa
-			dRH - Uncertainty in the relative humidity, default: 0.02
-			dCO2- Uncertainty in the CO2 density, default: 20
-			lat - [Optional] Latitude of the observatory in degrees
-			h 	- [Optional] Height in m of the observatory
+			Parameters
+			----------
+			z 	: float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			l1 	: float
+				Wavelength of light in microns, at the shortest wavelength 
+				of interest
+			l2 	: float
+				Wavelength of light in microns, at the longest wavelength 
+				of interest
+			T 	: float
+				Atmospheric temperature in Kelvin
+			p 	: float
+				Atmospheric pressure in Pascal
+			RH 	: float
+				Atmospheric relative humidity (0<=RH<=1)
+			xc 	: float
+				CO2 density in parts per million
+			dl1	: float 
+				Uncertainty in the wavelength in microns, at the shortest
+				wavelength of interest
+			dl2	: float 
+				Uncertainty in the wavelength in microns, at the longest
+				wavelength of interest
+			dT 	: float (default: 0.2)
+				Uncertainty in the atmospheric temperature in Kelvin
+			dP 	: float (default: 20)
+				Uncertainty in the atmospheric pressure in Pascal
+			dRH : float (default: 0.02)
+				Uncertainty in the atmospheric relative humidity
+			dCO2: float (default: 20)
+				Uncertainty in the atmospheric CO2 density in parts per million
+			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

AstroAtmosphere/examples.py

@@ -2,55 +2,119 @@
 from .ciddorModel import Observatory
 from .refractionModels import refraction
 from .dispersionModels import dispersion
+
 '''
-	Author: Joost van den Born
-	Some functions for quick testing.
+@ Author: 		Joost van den Born
+@ Contact: 		born@astron.nl
+@ Description: 	Some functions for quick testing.
 '''
 
 
-def quick_refractive_index(l, conditions='STANDARD'):
+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
+		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
+		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)
+	n 	= at.n_tph(l=l, T=T, p=p, RH=H, xc=xc)
 	return n
 
-def quick_refraction(l, z, conditions='STANDARD'):
+
+def quick_refraction(l, z, 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
+		z  	: 	float 
+			Zenith angle of the observation
+		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
+		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
+		p   = 71200     # Pa
 		H   = 0.22
 		xc  = 450       # ppm
 		lat = -24.5983  # degrees
@@ -60,31 +124,66 @@ def quick_refraction(l, z, conditions='STANDARD'):
 	at  = Observatory()
 
 	# Calculating indices of refraction for l
-	n 	= at.n_tph(l=l, T=T, p=P, RH=H, xc=xc)
+	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)
+	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)
+	ref.setReducedHeight(p, rho)
 	return ref.cassini(n, z)
 
-def quick_dispersion(l1, l2, z, conditions='STANDARD'):
+
+def quick_dispersion(l1, l2, z, 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.
+		z 	: 	float
+			Zenith angle of the observation
+		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
+		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
+		p   = 71200     # Pa
 		H   = 0.22
 		xc  = 450       # ppm
 		lat = -24.5983  # degrees
@@ -94,13 +193,13 @@ def quick_dispersion(l1, l2, z, conditions='STANDARD'):
 	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)
+	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)
+	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)
+	disp.setReducedHeight(p, rho)
 	return disp.cassini(n1, n2, z)

AstroAtmosphere/refractionModels.py

@@ -3,6 +3,7 @@ import numpy as np
 from scipy import integrate
 import os
 from .ciddorModel import Observatory
+
 '''
 @ Author:		Joost van den Born 
 @ Contact:		born@astron.nl
@@ -17,6 +18,61 @@ from .ciddorModel import Observatory
 '''
 
 class refraction(Observatory):
+	'''
+		Class that sets the observer location.
+		These location attributes are then used to calculate the 
+		appropriate atmospheric refraction.
+
+		See also:
+		Observatory()
+		dispersion()
+		...
+
+		Attributes
+		-----------
+		lat : 	float
+			Latitude of the observer in degrees
+		h 	: 	float
+			Altitude above mean sea level of the observer in meters
+		Other Observatory() attributes
+
+		Methods
+		-------
+		setReducedHeight(p, rho)
+			Sets the reduced atmospheric height.
+		planeParallel(n1, n2, zenith)
+			Returns the atmospheric refraction assuming a plane parallel
+			atmosphere
+		cassini(n1, n2, zenith, reduced_height=None)
+			Returns the atmospheric refraction asusming Cassini's homegeneous
+			atmosphere model
+		tan5(n1, n2, zenith, reduced_height=None)
+			Returns the atmospheric refraction following Oriani's theorem,
+			including a fifth order term in the expansion.
+		oriani(n1, n2, zenith, reduced_height=None)
+			Returns the atmospheric refraction following the classic version
+			of Oriani's theorem.
+		corbard(n1, n2, zenith, reduced_height=None)
+			Returns the atmospheric refraction following the method described 
+			in Corbard et al., 2019.
+		refractionIntegral(n1, n2, z0, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None)
+			Returns the atmospheric refraction following the refraction integral.
+		psi(x)
+			Submethod for corbard()
+		set_rc()
+			Sets the local radius of curvature of the Earth at the observer.
+		set_H(p, rho)
+			Sets the reduced height for the conditions at the observer, assuming
+			an adiabatic atmosphere.
+		H_isotherm(T)
+			Sets the reduced height for the conditions at the observer, assuming
+			an isothermal atmosphere.
+		cassiniError(z, l1, l2, T, p, RH, xc, dl1=0, dl2=0, dT=0.2, dP=20, 
+					dRH=0.02, dCO2=20, dz=0, lat=None, h=None)
+			Returns the uncertainty in the atmospheric refraction, assuming 
+			Cassini's homegeneous atmosphere model and the given parameters.		
+		Other Observatory() methods
+	'''
 	def __init__(self, lat, h):
 		Observatory.__init__(self)
 		self.h  	= h 					# height above sea level in meters
@@ -24,13 +80,36 @@ class refraction(Observatory):
 		self.rc 	= self.set_rc() * 1000	# Radius of curvature of the earth at the observer
 		# self.H2 	= self.H_isotherm(T) 	# Reduced height of isotherm atmosphere
 
-	def setReducedHeight(self, P, rho):
-		self.H  	= self.set_H(P, rho) 		# Reduced height of the atmosphere assuming ideal gas law
+	def setReducedHeight(self, p, rho):
+		'''
+			Sets the reduced atmospheric height.
+
+			Parameters
+			----------
+			p 	: float
+				Atmospheric pressure in Pascal
+			rho : float
+				Density of the atmospheric air at the observer
+		'''
+		self.H  	= self.set_H(p, rho) 		# Reduced height of the atmosphere assuming ideal gas law
 
 
-	def planeParallel(self, n, zenith):
+	def planeParallel(self, n1, zenith):
 		'''
 			Refraction based on a flat atmosphere, easily derived from Snell's law.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+
+			Returns
+			-------
+			R : float
+				Atmospheric refraction in degrees
 		'''
 		_R = np.arcsin(n * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
 		return np.degrees(_R)
@@ -40,6 +119,22 @@ class refraction(Observatory):
 	def cassini(self, n1, zenith, reduced_height=None):
 		'''
 			Refraction of spherical atmosphere, derived from geometry using Sine law
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			R : float
+				Atmospheric refraction in degrees
 		'''
 		if reduced_height:
 			H = reduced_height
@@ -55,6 +150,22 @@ class refraction(Observatory):
 		'''
 			Same as oriani's theorem, but including a higher order term.
 			Found in Corbard et al., 2019.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			R : float
+				Atmospheric refraction in degrees
 		'''
 		if reduced_height:
 			H = reduced_height
@@ -70,8 +181,25 @@ class refraction(Observatory):
 
 	def oriani(self, n1, zenith, reduced_height=None):
 		'''
-			Classic refraction formula with two tan(z) terms. This does not assume any 
-			information about the structure of the atmosphere. Found with Oriani's theorem.
+			Classic refraction formula with two tan(z) terms. This does not 
+			assume any information about the structure of the atmosphere. 
+			Found with Oriani's theorem.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			R : float
+				Atmospheric refraction in degrees
 		'''
 		if reduced_height:
 			H = reduced_height
@@ -90,6 +218,22 @@ class refraction(Observatory):
 			Corbard et al., 2019, mention an additional formula based on the error function.
 			Oriani's theorem can be derived from this equation by 'keeping only the three first
 			terms of its asymptotic expansion'.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the wavelength of interest
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			reduced_height : float (optional)
+				Reduced height of the atmosphere. Default value is calculated
+				from the object attributes.
+
+			Returns
+			-------
+			R : float
+				Atmospheric refraction in degrees
 		'''
 		if reduced_height:
 			H = reduced_height
@@ -103,14 +247,44 @@ class refraction(Observatory):
 
 
 
-	def refractionIntegral(self, n0, z0, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None):
+	def refractionIntegral(self, n1, z0, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None):
 		'''
-			Using data that gives temperature, pressure and density as a function of height, 
-			the gladstone-dale relation is invoked to calculate the path of a light ray 
-			along the atmosphere. This is allows us to also include atmospheric activity in 
+			Using data that gives temperature, pressure and density as a 
+			function of height, the gladstone-dale relation is invoked to
+			calculate the path of a light ray along the atmosphere. 
+			This is allows us to also include atmospheric activity in 
 			different atmosphere layers.
 			Recommended use with the US Standard Atmosphere from 1976.
 			USSA1976 data generated from code supplied by http://www.pdas.com/atmos.html
+			The method described in eq. 3 in Auer & Standish (2000) is used for the 
+			calculation.
+
+			Parameters
+			----------
+			n1 : float
+				Refractive index at the wavelength of interest
+			z0 : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			R0 : float (optional)
+				Local radius of curvature of the Earth in km.
+			useStandardAtmosphere : bool
+				Flag to use the included USSA1976 data
+			heightData : 1D Array
+				If "useStandardAtmosphere" is False, the function will ask
+				for alternative atmospheric data to integrate over. This
+				gives the altitudinal points
+			rhoData : 1D Array
+				If useStandardAtmosphere is False, the function will ask
+				for alternative atmospheric data to integrate over. This
+				gives the atmospheric density at the points defined by 
+				"heightData".
+
+			Returns
+			-------
+			R : float
+				Atmospheric refraction in degrees
+
 		'''
 		if not R0:
 			R0 = self.rc
@@ -148,6 +322,19 @@ class refraction(Observatory):
 
 
 	def psi(self, x):
+		'''
+			Submethod for corbard(). Described explicitly in Corbard et al.,2019.
+
+			Parameters
+			----------
+			x 	: 	float
+				Input of the function
+
+			Returns
+			-------
+			y 	: 	float
+				Output of the function
+		'''
 		_f   = lambda a: np.exp(-a**2)
 		_int = integrate.quad(_f, x, np.inf)
 		return np.exp(x**2) * _int[0]
@@ -159,6 +346,11 @@ class refraction(Observatory):
 			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
@@ -169,9 +361,21 @@ class refraction(Observatory):
 
 
 
-	def set_H(self, P, rho):
+	def set_H(self, p, rho):
 		'''
 			Calculates the reduced height of the atmosphere, assuming ideal gas law.
+
+			Parameters
+			----------
+			p 	: float
+				Atmospheric pressure in Pascal
+			rho : float
+				Density of the atmospheric air at the observer
+
+			Returns
+			-------
+			H 	: float
+				Reduced height of the atmosphere in meters
 		'''
 		_phi = np.deg2rad(self.lat)
 		_c1 = 5.2790414e-3
@@ -181,11 +385,22 @@ class refraction(Observatory):
 		_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) * self.h + 7.2e-13 * self.h**2
-		return P / (rho * _g)
+		return p / (rho * _g)
 
 	def H_isotherm(self, T):
 		'''
-			Calculates the reduced height of the atmosphere, assuming an isotherm distribution.
+			Calculates the reduced height of the atmosphere, 
+			assuming an isotherm distribution.
+
+			Parameters
+			----------
+			T 	: float
+				Atmospheric temperature in Kelvin
+
+			Returns
+			-------
+			H 	: float
+				Reduced height of the atmosphere in meters
 		'''
 		_kb = 1.38064852e-23 # kg m2 s-2 K-1
 		_m  = 4.809651698e-26 # kg (weight per molecule of air, assuming dry air)
@@ -203,6 +418,28 @@ class refraction(Observatory):
 	def cassiniError(self, n, dn, zenith, dz=0, reduced_height=None):
 		'''
 		Calculates the uncertainty in the refraction angle in degrees.
+		See Observatory().n_tph() and Observatory().dn_tph() on how to
+		calculate the refractive index and it's uncertainty.
+
+		Parameters
+		----------
+		n 	: 	float
+			Refractive index of air
+		dn 	: 	float
+			Uncertainty in the refractive index of air
+		zenith : float
+			Angle of observation in degrees.
+			zenith = 0 means the observation is directly overhead.
+		dz  : 	float (optional)
+			Uncertainty in the zenith angle, in degrees.
+		reduced_height : float (optional)
+			Reduced height of the atmosphere. Default value is calculated
+			from the object attributes.
+
+		Returns
+		-------
+		dR 	: float
+			Uncertainty in the atmospheric refraction in degrees.
 		'''
 		if reduced_height:
 			_H = reduced_height

AstroAtmosphere/refractivityModels.py

@@ -1,7 +1,6 @@
 #!/usr/bin/env Python3
 import numpy as np
 
-
 '''
 @ Author:		Joost van den Born 
 @ Contact:		born@astron.nl
@@ -17,14 +16,21 @@ def BarrellAndSears(l, T=288.15, p=101325, H=0):
 		For atmospheric air (i.e. CO2 = 300 ppm)
 		This equation is applicable for the  ranges T=10-30C and p=720-800 mmHg.
 
-		Input: 
-		l - Wavelength in microns
-		T - Temperature in Kelvin
-		p - Atmospheric pressure in Pa
-		H - fractional humidity (0=<H=<1)
+		Parameters
+		---------- 
+		l 	:	float
+			Wavelength in microns
+		T 	: 	float
+			Temperature in Kelvin
+		p 	: 	float
+			Atmospheric pressure in Pascal
+		H 	: 	float
+			Relative humidity (0=<H=<1)
 
-		Returns:
-		n - Refractive index of atmospheric air at the given conditions
+		Returns
+		-------
+		n 	: 	float
+			Refractive index of atmospheric air at the given conditions
 	'''
 	T = K2C(T)
 	P = hPA2mmHg(p/100)
@@ -39,11 +45,15 @@ def Ciddor(l, CO2=450):
 		At default conditions only. For more flexibility use
 		Observatory.n_tph()
 
-		Input: 
-		l - Wavelength in microns
+		Parameters
+		---------- 
+		l 	: 	float
+			Wavelength in microns
 
-		Returns:
-		n - Refractive index of atmospheric air at the given conditions
+		Returns
+		-------
+		n 	: 	float
+			Refractive index of atmospheric air at the given conditions
 	'''
 	s = 1 / l
 	R = 5792105 / (238.0185 - s**2) + 167917 / (57.362 - s**2)
@@ -59,11 +69,15 @@ def Edlen1953(l):
 
 		At standard conditions only
 
-		Input: 
-		l - Wavelength in microns
+		Parameters
+		---------- 
+		l 	: 	float
+			Wavelength in microns
 
-		Returns:
-		n - Refractive index of atmospheric air at the given conditions
+		Returns
+		-------
+		n 	: 	float
+			Refractive index of atmospheric air at the given conditions
 	'''
 	s = 1 / l
 	R = 6432.8 + 2949810 / (146 - s**2) + 25540/(41-s**2)
@@ -77,11 +91,15 @@ def Edlen1966(l):
 
 		At standard conditions only
 
-		Input: 
-		l - Wavelength in microns
+		Parameters
+		---------- 
+		l 	: 	float
+			Wavelength in microns
 
-		Returns:
-		n - Refractive index of atmospheric air at the given conditions
+		Returns
+		-------
+		n 	: 	float
+			Refractive index of atmospheric air at the given conditions
 	'''
 	s = 1 / l
 	R = 8342.13 + 2406030 / (130 - s**2) + 15997 / (38.9 - s**2)
@@ -95,14 +113,21 @@ def Owens(l, T=288.15, p=101325, H=0):
 		Note that the pressure P and H are in torricelli (Torr),
 		however, this is almost precisely 1 Torr = 1 mmHg
 
-		Input:
-		l - wavelength in micron
-		T - Temperature in Kelvin
-		p - Atmospheric pressure in Pa
-		H - Humdity (0<=H<=1)
+		Parameters
+		----------
+		l 	: 	float
+			wavelength in micron
+		T 	: 	float
+			Temperature in Kelvin
+		p 	: 	float
+			Atmospheric pressure in Pa
+		H 	: 	float
+			Humdity (0<=H<=1)
 
-		Returns:
-		n - Refractive index of atmospheric air at the given conditions
+		Returns
+		-------
+		n 	: 	float
+			Refractive index of atmospheric air at the given conditions
 	'''
 	P = hPA2mmHg(p/100)
 	T = K2C(T)
@@ -116,14 +141,21 @@ def BonschPotulski(l, T=288.15, p=101325, CO2=450):
 		From Bonsch & Potulski, 1998
 		Equations 6, 7, 8
 		
-		Input: 
-		l - Wavelength in microns
-		T - Temperature in Kelvin
-		p - Atmospheric pressure in Pa
-		CO2 - Particles per million of CO2
+		Parameters
+		---------- 
+		l 	: 	float
+			Wavelength in microns
+		T 	: 	float
+			Temperature in Kelvin
+		p 	: 	float
+			Atmospheric pressure in Pa
+		CO2 : 	float
+			Particles per million of CO2
 
-		Returns:
-		n - Refractive index of atmospheric air at the given conditions
+		Returns
+		-------
+		n 	: 	float
+			Refractive index of atmospheric air at the given conditions
 	'''
 	x = CO2 / 1e6
 	T = K2C(T)
@@ -143,11 +175,15 @@ def BirchDowns(l):
 
 		Only at standard conditions
 
-		Input:
-		l - wavelength in micron
+		Parameters
+		----------
+		l 	: 	float
+			wavelength in micron
 
-		Returns:
-		n - Refractive index of atmospheric air at the given conditions
+		Returns
+		-------
+		n 	: 	float
+			Refractive index of atmospheric air at the given conditions
 	'''
 	s = 1 / l
 	R = 8342.54 + 2406147 / (130 - s**2) + 15998 / (38.9-s**2)
@@ -161,11 +197,15 @@ def PeckReeder(l):
 
 		Only at standard conditions
 
-		Input:
-		l - wavelength in micron
+		Parameters
+		----------
+		l 	: 	float
+			wavelength in micron
 
-		Returns:
-		n - Refractive index of atmospheric air at the given conditions
+		Returns
+		-------
+		n 	: 	float
+			Refractive index of atmospheric air at the given conditions
 	'''
 	s = 1 / l
 	R = 8060.51 + 2480990 / (132.274 - s**2) + 17455.7 / (39.32957 - s**2)
@@ -180,18 +220,29 @@ def slalib(l, T=288.15, p=101325, H=0, h=0, lat=-24, tlr=0.0065, eps=1e-10):
 		They take the B&S refraction formula as a basis.
 		Requires the pySLALIB package to work.
 
-		Input: 
-		l - Wavelength in microns
-		T - Temperature in Kelvin
-		p - Atmospheric pressure in Pa
-		H - fractional humidity (0=<H=<1)
-		h - Altitude of the observer (default=0 m)
-		lat - Latitude of the observer (default=-24)
-		tlr - temperature lapse rate (default=0.0065)
-		eps - precision at which the calculation is aborted (default=1e-10)
+		Parameters
+		---------- 
+		l 	:	float
+			Wavelength in microns
+		T 	:	float
+			Temperature in Kelvin
+		p 	:	float
+			Atmospheric pressure in Pa
+		H 	:	float
+			fractional humidity (0=<H=<1)
+		h 	:	float
+			Altitude of the observer (default=0 m)
+		lat :	float
+			Latitude of the observer (default=-24)
+		tlr :	float
+			temperature lapse rate (default=0.0065)
+		eps :	float
+			precision at which the calculation is aborted (default=1e-10)
 
-		Returns:
-		n - Refractive index of atmospheric air at the given conditions
+		Returns
+		-------
+		n 	:	float
+			Refractive index of atmospheric air at the given conditions
 	'''
 	try:
 	    from pyslalib import slalib
@@ -215,17 +266,27 @@ def slalib(l, T=288.15, p=101325, H=0, h=0, lat=-24, tlr=0.0065, eps=1e-10):
 def Mathar(l, T=288.15, p=101325, H=0):
 	'''
 		Refractivity model by Mathar (2007), valid above 1.4 microns
-		Code taken and modified from: 
+		Code taken from
 		https://github.com/polyanskiy/refractiveindex.info-scripts/blob/master/scripts/Mathar%202007%20-%20Air%201.3-2.5um.py
+		(better known from refractiveindex.info). The code was subsequently 
+		minimally modified to make it suitable for this package.
 
-		Input: 
-		l - Wavelength in microns
-		T - Temperature in Kelvin
-		p - Atmospheric pressure in Pa
-		H - fractional humidity (0=<H=<1)
 
-		Returns:
-		n - Refractive index of atmospheric air at the given conditions
+		Parameters
+		---------- 
+		l 	: 	float
+			Wavelength in microns
+		T 	: 	float
+			Temperature in Kelvin
+		p 	: 	float
+			Atmospheric pressure in Pa
+		H 	: 	float
+			fractional humidity (0=<H=<1)
+
+		Returns
+		-------
+		n 	: 	float
+			Refractive index of atmospheric air at the given conditions
 	'''
 	if l < 1.3 or l > 24:
 		print('Using the Mathar model at this wavelength (l={:.3f}) is not recommended!'.format(l))
@@ -346,11 +407,32 @@ def Mathar(l, T=288.15, p=101325, H=0):
 def K2C(K):
 	'''
 		Kelvin to Celsius conversion
+
+		Parameters
+		----------
+		K : float
+			Temperature in Kelvin
+
+		Returns
+		-------
+		C : float
+			Temperature in Celsius
 	'''
 	return K - 273.15
 
 def hPA2mmHg(hPa):
 	'''
 		hPa to mmHg conversion
+
+		Parameters
+		----------
+		hPa 	:	float
+			Atmospheric pressure in hPa (hecto Pascal)
+
+		Returns
+		-------
+		mmHg 	: 	float
+			Atmospheric pressure in mmHg
 	'''
 	return 0.7500616827 * hPa
+

README.md

@@ -1,66 +1,128 @@
+# AstroAtmosphere
+
 **Note:**
-A manuscript of this atmospheric code is in preparation and will be submitted to the Monthly Notices of the Royal Astronomical Society (MNRAS) journal before the end of the year (2019). 
-If you use this code for your project, please acknowledge this by making a reference to this paper in due course. You can contact me, the author, for the submission and review status.
+It would be appreciated if a reference to the following work, for which this package was originally build, is included whenever this code is used for a publication: 
+*Quantification of the expected residual dispersion of the MICADO Near-IR imaging instrument*, van den Born & Jellema, 2020, MNRAS. Accepted 2020 June 23. Received 2020 June 23; in original form 2020 January 24.
 
 
-Author info: *Joost van den Born*, *born@astron.nl*
+**Author info**: *Joost van den Born*, *born@astron.nl*
 
 
+## Example
 
-*Below you can find a set of models describing the atmospheric refraction and dispersion.*
+A minimal working example to calculate the atmospheric dispersion at Cerro Armazones between 1.49 and 1.78 micron.
 
-**Contents:**
-- AstroAtmosphere         	    : Python package that can be used to calculate the atmospheric refraction and dispersion, the main content of this repository.
-- atmospheric refraction example: jupyter notebook example
-- `setup.py`                      : Install script to install the AstroAtmosphere package locally.
+```python
+from AstroAtmosphere import *
+
+# Parameters at Cerro Armazones
+T   = 279.65    # k
+P   = 71200     # Pa
+H   = 0.22
+xc  = 450       # ppm
+lat = -24.5983  # degrees
+h   = 3064      # m
+l1  = 1.49      # micron
+l2  = 1.78      # micron
+z0  = 30
+
+# 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)
+
+# Calculation of the atmopheric dipsersion
+atm_dispersion = ref.cassini(n1, z0) * 3.6e6
+print ('The dispersion is %.03f milli arc seconds' %(atm_dispersion))
+```
+
+## Installation
+
+The package is available through Git-Lab and through PyPI. The package may be installed using the `setup.py` file, but easier is to install using `pip`.
+
+```bash
+pip install AstroAtmosphere
+```
 
 
-**On the contents of the AstroAtmosphere package**
+##Contents of the AstroAtmosphere package
+
 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.
-- `refractionModels.py` 		: Contains various atmospheric models that can be used to calculate the atmospheric refraction
+- `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).
 
 
-**On the dispersion model:**
+##The Ciddor refractivity model
 
-Dispersion model from Ciddor, 1996, calculates the index of refraction of air using the following parameters:
-- Temperature
-- Pressure
-- Relative humidity
-- CO2 density
+The refractivity model from Ciddor, 1996, calculates the index of refraction of air using the following parameters:
+- T - Temperature in Kelvin
+- p - Pressure in Pascal
+- RH - Relative humidity in fractional units
+- xc - CO2 density in parts per million
+
+It is invoked by first initializing an `Observatory()` object. The atmospheric conditions can be parsed into the `n_tph()` function.
+
+```python
+at  = Observatory()
+n   = at.n_tph(l=0.55, T=273.15, p=101325, RH=0.2, xc=300)
+```
+
+Basic functionality of other refractivity models is included as well and is discussed below.
 
 
+##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.
 
 
-**On the atmospheric refraction / atmospheric dispersion models:**
-
-Available models:
-- Full integration using the refractive integral (requires detailed information about the atmosphere, e.g. `AstroAtmosphere/data/1976USSA.txt`)
-- Plane-parallel atmosphere refraction model
-- Cassini's homogeneous atmosphere refraction model
-- Oriani's theorem ( Atan(z) + Btan^3(z) )
-- Oriani's theorem expanded ( Atan(z) + Btan^3(z) + Ctan^5(z)  )
-- Error function refraction model
-
-
+The available models in this package:
+- `refractionIntegral()` - Full integration using the refractive integral (requires detailed information about the atmosphere, e.g. `AstroAtmosphere/data/1976USSA.txt`)
+- `planeParallel()` - Plane-parallel atmosphere refraction model
+- `cassini()` - Cassini's homogeneous atmosphere refraction model
+- `oriani()` - Oriani's theorem ( Atan(z) + Btan^3(z) )
+- `tan5()` - Oriani's theorem expanded ( Atan(z) + Btan^3(z) + Ctan^5(z)  )
+- `corbard()` - Error function refraction model
 
 
 Input parameters for the spherical atmosphere models (all but the first):
-- Temperature               [K]
-- Pressure                  [Pa]
-- Relative humidity
-- CO2 density               [ppm]
-- Latitude                  [deg]
-- Height above sea level    [m]
-- Wavelength(s)             [um]
-- Zenith angle              [deg]
+- T - Temperature               [K]
+- p - Pressure                  [Pa]
+- RH - Relative humidity
+- xc - CO2 density              [ppm]
+- lat - Latitude                [deg]
+- h - Height above sea level    [m]
+- l - Wavelength(s)             [um]
+- 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.
+
+```python
+# Density of the atmosphere (following CIPM-81/91 equations)
+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)
+
+# Calculating the atmospheric refraction
+atm_ref = ref.cassini(n, zenith=30)
+```
 
 
-**Analytical error propagation**
+##Analytical error propagation
 
 Full analytical error propagation is available for the Ciddor dispersion model (`Observatory.dn_tph()`), for the Cassini refraction model (`refraction.cassiniError()`) and for the Cassini dispersion model (`dispersion.cassiniError()`).
 The results agree well with a Monte Carlo simulation, suggesting that if additional errors have not been considered, they are of neglible importance.
@@ -97,3 +159,52 @@ Cassini dispersion model:
 | dCO2 | CO2 density | ppm |
 | dz | Zenith angle | deg |
  
+An example of its use is given in the code snippet here:
+```python
+ # Initializing dispersion object
+disp = dispersion(lat, h)
+disp.setReducedHeight(P, rho)
+
+# Atmospheric dispersion
+dispersion_obs = disp.cassini(n1, n2, zenith_obs)
+
+# Uncertainty in the atmospheric dispersion
+dispersion_unc = disp.cassiniError(zenith_obs, l1, l2, T, P, H, xc, 
+    dl1=0.001, dl2=0.001, dT=0.2, dP=30, dRH=0.03, dCO2=25, dz=1/3600) * 3.6e6
+```
+
+##Other refractivity models
+
+The package also includes limited functionality of other refractivity models. These are not as extensive, work only under standard atmospheric conditions, or over a limited wavelength range, compared to the Ciddor model.
+
+Included are
+| Model name | Function | Remarks |
+| ---------- | -------- | ------- |
+| 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|
+|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.|
+
+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. 
+
+Presently, only the `'STANDARD'` and `'CERRO_ARMAZONES'` flags are known. Contact the author to add additional standard conditions.
+
+```python
+ # Quick refractive index
+nq = quick_refractive_index(l=1.49, conditions='STANDARD')
+
+# Quick refraction
+rq = quick_refraction(l1, zenith_obs, conditions='STANDARD')
+
+# Quick dispersion
+dq = quick_dispersion(l1, l2, zenith_obs, conditions='STANDARD')
+```

scripts/ADC.py

@@ -0,0 +1,193 @@
+#!/usr/bin/env Python3
+import numpy as np
+from scipy.optimize import fsolve
+
+
+class ADC:
+	def __init__(self, wedge=(1.6554, 6.28927, -1.48306)):
+		self.p1 = prism('S-LAH71', 15, wedge[0], wedge[1])
+		self.p2 = prism('S-FPL51', 15, wedge[1], wedge[2])
+		self.p3 = prism('S-FPL51', 15, wedge[2], wedge[1])
+		self.p4 = prism('S-LAH71', 15, wedge[1], wedge[0])
+		self.theta = 90
+
+	def setOrientation(self, theta):
+		self.p1.setPrismOrientation(theta)
+		self.p2.setPrismOrientation(theta)
+		self.p3.setPrismOrientation(180-1*theta)
+		self.p4.setPrismOrientation(180-1*theta)
+		self.theta = theta
+
+	def setFieldRotation(self, FR):
+		self.p1.setPrismOrientation(self.theta-FR)
+		self.p2.setPrismOrientation(self.theta-FR)
+		self.p3.setPrismOrientation(180-1*self.theta-FR)
+		self.p4.setPrismOrientation(180-1*self.theta-FR)
+
+	def calculateDispersion(self, l1, l2, hx=0, hy=0, refraction=False):
+		'''
+			Calculates the dispersion for the ADC using the current position, 
+			for a normalized on sky field angle (hx, hy). Default field position
+			is in the center of the FoV.
+		'''
+		_prismList = [self.p1, self.p2, self.p3, self.p4]
+		_d1x_list, _d1y_list = [], []
+		_pupilRatio = 0.07943/38.542
+
+		hx = hx * 7.5e-3
+		hy = hy * 7.5e-3
+		ix = hx * 1/_pupilRatio
+		iy = hy * 1/_pupilRatio
+
+		for l in [l1,l2]:
+			_deltaListX = np.zeros(len(_prismList))
+			_deltaListY = np.zeros(len(_prismList))
+			for p in range(4):
+				_n_g = _prismList[p].indexOfRefraction(l)
+				_b1x = _prismList[p].tilt1x
+				_b1y = _prismList[p].tilt1y
+
+				_b2x = _prismList[p].tilt2x
+				_b2y = _prismList[p].tilt2y
+
+				if p==0:
+				    _t_iz_n_x = ix # Zero before first prism
+				    _t_iz_n_y = iy # Zero before first prism
+
+				# X direction    
+				_t_i1_n_x = _t_iz_n_x
+				_d1x = deviationCalculation3(_t_i1_n_x, _b1x, _b2x, _n_g)
+				_t_tz_n_x = _t_iz_n_x - _d1x
+				_t_iz_n_x = _t_tz_n_x         # Set output to new input (for the next prism)
+
+				# Y direction
+				_t_i1_n_y = _t_iz_n_y
+				_d1y = deviationCalculation3(_t_i1_n_y, _b1y, _b2y, _n_g)
+				_t_tz_n_y = _t_iz_n_y - _d1y
+				_t_iz_n_y = _t_tz_n_y         # Set output to new input (for the next prism)
+
+				# Adding the deviations to the list with delta terms
+				_deltaListX[p] = _d1x 
+				_deltaListY[p] = _d1y
+			_d1x_list.append(_deltaListX.sum())
+			_d1y_list.append(_deltaListY.sum())
+
+		self.dispX = ( _d1x_list[1] - _d1x_list[0] ) * 3600 * 1000 * _pupilRatio
+		self.dispY = ( _d1y_list[1] - _d1y_list[0] ) * 3600 * 1000 * _pupilRatio
+		
+		if refraction:
+			rx = ( _d1x_list[0] ) * 3600 * 1000 * _pupilRatio
+			ry = ( _d1y_list[0] ) * 3600 * 1000 * _pupilRatio
+			return rx, ry
+		return self.dispX, self.dispY
+
+	def optimize(self, inputDispersion, l1, l2, x0=50):
+		_solver = fsolve(optimizeFunction, x0=x0, args=(self, inputDispersion, l1, l2))
+		return _solver[0]
+
+
+
+
+class prism:
+	def __init__(self, glass, thickness, tilt1, tilt2=None):
+		self.glass 		= glass
+		self.thickness 	= thickness
+		self.tilt1x		= tilt1
+		self.tilt1y 	= 0
+		self.tilt2x		= tilt2
+		self.tilt2y 	= 0
+
+		# Remember original tilts for reset:
+		self.tilt1x_OG	= tilt1
+		self.tilt1y_OG 	= 0
+		self.tilt2x_OG	= tilt2
+		self.tilt2y_OG 	= 0
+
+
+	def indexOfRefraction(self, l):
+		if self.glass == 'S-FPL51':
+			a = [1.170105050E+00, 4.757107830E-02, 7.638324450E-01]
+			b = [6.162039240E-03, 2.633728760E-02, 1.418826420E+02]
+		elif self.glass == 'S-LAH71':
+			a = [1.982800310E+00, 3.167584500E-01, 2.444726460E+00]
+			b = [1.189874560E-02, 5.271560010E-02, 2.132206970E+02]
+		return dispersionFormula(l, a, b) 
+
+	def setPrismOrientation(self, theta):
+		self.tilt1x, self.tilt1y = tiltRotation(self.tilt1x_OG, self.tilt1y_OG, theta)
+		self.tilt2x, self.tilt2y = tiltRotation(self.tilt2x_OG, self.tilt2y_OG, theta)
+	
+
+
+def tiltRotation(tx,ty, theta):
+	theta = np.deg2rad(theta-90)
+	tx_new = tx * np.cos(theta) - ty * np.sin(theta)
+	ty_new = tx * np.sin(theta) + ty * np.cos(theta)
+	return tx_new, ty_new
+
+
+def dispersionFormula(l, a, b):
+	n = np.sqrt( (a[0]* l**2)/(l**2 - b[0]) + (a[1] * l**2) / (l**2 - b[1]) + (a[2] * l**2) / (l**2 - b[2]) + 1)
+	return n
+
+
+# def deviationCalculation(t_i1, b1, b2, n):
+    # Returns deviation of the ray after the prism with apex angle alpha
+    # Units: degrees
+    # DO NOT USE!
+    # alpha = b2 - b1
+    # delta = t_i1 + np.arcsin( np.sin(np.deg2rad(alpha)) * np.sqrt(n**2-(np.sin(np.deg2rad(t_i1)))**2) - np.sin(np.deg2rad(t_i1)) * np.cos(np.deg2rad(alpha)) )*180/np.pi - alpha
+    # return delta
+
+
+def deviationCalculation2(t_i1, b1, b2, n):
+	# Returns deviation of the ray after the prism with apex angle alpha
+	# only now using Snell's law only
+    # Units: degrees
+    t_i1  = np.deg2rad(t_i1)
+    b1    = np.deg2rad(b1)
+    b2    = np.deg2rad(b2)
+
+    t_i   = b1 - t_i1
+    t_f   = np.arcsin( 1 / n * np.sin( t_i ) )
+    t_f_z = b1 - t_f
+
+    t_i   = b2 - t_f_z
+    t_f   = np.arcsin( n / 1 * np.sin( t_i ) )
+    t_f_z = b2 - t_f
+
+    delta = t_i1 - t_f_z
+    delta = np.degrees(delta)
+    return delta
+
+def deviationCalculation3(t_i1, b1, b2, n):
+	# Using the formula from Hagen & Tkazcyk (2011)
+	# Gives exactly the same results as deviationCaclulation2(), above
+	t_i1 	= np.deg2rad(t_i1)
+	b1 		= np.deg2rad(b1)
+	b2 		= np.deg2rad(b2)
+	alpha 	= b2 - b1
+	delta 	= t_i1 - b2 - np.arcsin( n * np.sin( np.arcsin(1/n * np.sin(t_i1 - b1) ) - alpha ) )
+	return np.degrees(delta)
+
+def optimizeFunction(theta, adcObject, inputDispersion, l1, l2, FR=None):
+	if FR:
+		FRdeg = np.deg2rad(FR) # 30 degrees field rotation
+		inputDispersionX = np.cos(FRdeg) * inputDispersion
+		inputDispersionY = np.sin(FRdeg) * inputDispersion
+		adcObject.setOrientation(theta)
+		ADCdisp = adcObject.calculateDispersion(l1,l2)
+		return np.sqrt((ADCdisp[1] - inputDispersionY)**2 + (ADCdisp[0] - inputDispersionX)**2)
+	else:	
+		adcObject.setOrientation(theta)
+		ADCdisp = adcObject.calculateDispersion(l1,l2)
+		ADCdisp = ADCdisp[1]
+		return ADCdisp - inputDispersion
+
+def angleFunc(x, a):
+	'''
+		Function to fit the ideal angle on
+	'''
+	return np.arccos(a * np.tan(np.deg2rad(x))) * 180/np.pi
+
+

scripts/ADC_instruction.ipynb

@@ -0,0 +1,156 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Atmospheric Dispersion Corrector model\n",
+    "\n",
+    "This notebook explains how you can use the ADC.py analytic ADC model for some quick simulations. For the model to be representative, the following assumptions must apply:\n",
+    "* The ADC consists of two amici prisms, totaling 4 glass prisms in total\n",
+    "* The ADC is located in the pupil plane of the instrument (or very close to it)\n",
+    "* The rotation of the prisms is defined as follows: for $\\theta=0$ the ADC is in it's maximum dispersive configuration. For any other rotation angle, the top prism set rotates by $\\theta$ degrees and the bottom prism set also rotates by $\\theta$ degrees in opposite direction. The total relative angle is thus $2\\theta$. From this also follows that the minimum dispersion magnitude is reached at $\\theta=90$ degrees.\n",
+    "\n",
+    "The ADC.py model is really not made for very general use, so it will likely require some editing in the source code by the user to make it work for a different instrument than MICADO.\n",
+    "\n",
+    "\n",
+    "We first start by importing the necessary files. Here we'll use the `AstroAtmosphere` package to calculate the atmospheric dispersion, but a manual input is also possible."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import ADC\n",
+    "import AstroAtmosphere"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Then we calculate the atmospheric dispersion for a given set of conditions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Atmospheric conditions\n",
+    "T   = 288.15 \t# k\n",
+    "P   = 101325  \t# Pa\n",
+    "H   = 0.3\n",
+    "xc  = 400   \t# ppm\n",
+    "lat = -25 \t\t# degrees\n",
+    "h   = 2500  \t# m\n",
+    "l1  = 1.49 \t\t# micron, start H-band\n",
+    "l2  = 1.78 \t\t# micron, end H-band\n",
+    "zenith = 30 \t# degrees\n",
+    "\n",
+    "# Initializing dispersion model\n",
+    "at  = AstroAtmosphere.Observatory()\n",
+    "rho = at.rho(p=P, T=T, RH=H, xc=xc)\n",
+    "disp = AstroAtmosphere.dispersion(lat, h)\n",
+    "disp.setReducedHeight(P, rho)\n",
+    "\n",
+    "# Calculation of indices of refraction at l1 and l2\n",
+    "n1 = at.n_tph(l=l1, T=T, p=P, RH=H, xc=xc)\n",
+    "n2 = at.n_tph(l=l2, T=T, p=P, RH=H, xc=xc)\n",
+    "\n",
+    "# Calculation of atmospheric dispersion in milli arc sec\n",
+    "dispersion = disp.cassini(n1, n2, zenith)*3.6e6"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now, we initialize the ADC class. We can calculate the dispersion in sky angles, i.e. in terms of the Field of View of the telescope, for a given rotation of the ADC prisms. In this case, if we set $\\theta=30$ degrees, we find that the ADC produces 58 milli arc seconds of dispersion in H-band."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "The dispersion of the ADC for 30 degrees of rotation: 58.011 mas\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Initialize ADC object\n",
+    "adc = ADC.ADC()\n",
+    "\n",
+    "# Calculate the dispersion for a given rotation\n",
+    "theta = 30 # degrees\n",
+    "adc.setOrientation(theta)\n",
+    "adcDispersion = adc.calculateDispersion(l1,l2)[1]\n",
+    "print(\n",
+    "    'The dispersion of the ADC for {:.0f} degrees of rotation: {:.3f} mas'\\\n",
+    "    .format(theta, adcDispersion))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Of course, we're often much more interested in the optimum positioning of the ADC given an amount of atmospheric dispersion. This is done as follows."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "The atmospheric dispersion: 24.771 mas\n",
+      "The optimum ADC angle in this case: 68.1 degrees\n",
+      "The residual dispersion: 0.000 mas\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Optimize for atmospheric dispersion, given wavelength range {l1,l2}\n",
+    "adcOptimumAngle = adc.optimize(dispersion, l1, l2) \n",
+    "adcDispersionOptimized = adc.calculateDispersion(l1,l2)[1]\n",
+    "\n",
+    "print('The atmospheric dispersion: {:.3f} mas'.format(dispersion) )\n",
+    "print('The optimum ADC angle in this case: {:.1f} degrees'.format(adcOptimumAngle))\n",
+    "print('The residual dispersion: {:.3f} mas'.format(dispersion - adcDispersionOptimized))\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

scripts/ADCsteps.py

@@ -0,0 +1,87 @@
+#!/usr/bin/env Python3
+import numpy as np
+from AstroAtmosphere import *
+import matplotlib.pyplot as plt
+
+'''
+@ Author:		Joost van den Born 
+@ Contact:		born@astron.nl
+@ Description:	This script uses the equations of the paper to calculate the 
+				residual dispersion due to the finite step size of the ADC.
+'''
+
+
+def angleFunc(x, a):
+	'''
+		Function to fit the ideal angle on
+	'''
+	return np.arccos(a * np.tan(np.deg2rad(x))) * 180/np.pi
+
+
+def angleFuncDeriv(x, a):
+	return (a * (1+np.tan(np.deg2rad(x))**2)) / np.sqrt(1 - a**2 * np.tan(np.deg2rad(x))**2)
+
+
+def cassiniDispersionDeriv(n1, n2, zenith, alpha):
+	a1 = n1 * alpha * np.cos(np.deg2rad(zenith)) / np.sqrt(1 - n1**2 * alpha**2 * np.sin(np.deg2rad(zenith))**2)
+	a2 = n2 * alpha * np.cos(np.deg2rad(zenith)) / np.sqrt(1 - n2**2 * alpha**2 * np.sin(np.deg2rad(zenith))**2)
+	return a1 - a2
+
+
+# Parameters at Cerro Armazones
+T   = 279.65 	# k
+P   = 71200  	# Pa
+H   = 0.22
+xc  = 450   	# ppm
+lat = -24.5983 	# degrees
+h   = 3064  	# m
+l1  = 1.49 		# micron
+l2  = 1.78 		# micron
+zenith = 30 	# degrees
+
+
+
+at = dispersion(lat, h)
+at.setReducedHeight(P, at.rho(P, T, H, xc))
+
+n1 = at.n_tph(l1, T, P, H, xc)
+n2 = at.n_tph(l2, T, P, H, xc)
+
+alpha = at.rc / (at.rc+at.h)
+# c_fH = 0.499413 # From Zemax
+c_fH = 0.46816061517198737 # Completely analytical cf
+# f_fH = 0.46462883308144837 # Best fit to analytical ADC model
+
+
+dispersion 		= at.cassini(n1, n2, zenith)*3.6e6
+dispersionDeriv = cassiniDispersionDeriv(n1, n2, zenith, alpha)*3.6e6
+optAngle 		= angleFunc(zenith, c_fH)
+optAngleDeriv	= angleFuncDeriv(zenith, c_fH)
+maxDev 			= optAngleDeriv / dispersionDeriv
+print('Zenith = %s degrees' %zenith)
+print('Dispersion in mas: ', dispersion)
+print('Dispersion/degree zenith at z=%s: %s mas/degree z' %(zenith, dispersionDeriv))
+print('Optimum angle ADC: %s' %(optAngle))
+print('Optimum angle ADC/degree zenith at z=%s: %s' %(zenith, optAngleDeriv))
+print('Max deviation per mas dispersion: %s degrees' %(maxDev))
+print('At designed stepsize at most %s mas of residual dispersion' %(1/maxDev*(360/1262.5)))
+
+
+# Plotting as a function of zenith angle
+dataList = []
+for zenith in np.linspace(5,65,100):
+	dispersion 		= at.cassini(n1, n2, zenith)*3.6e6
+	dispersionDeriv = cassiniDispersionDeriv(n1, n2, zenith, alpha)*3.6e6
+	optAngle 		= angleFunc(zenith, c_fH)
+	optAngleDeriv	= angleFuncDeriv(zenith, c_fH)
+	maxDev 			= optAngleDeriv / dispersionDeriv
+	dataList.append(1/maxDev*(360/1262.5))
+
+
+fig = plt.figure()
+ax = fig.add_subplot(111)
+ax.plot(np.linspace(5,65,100), dataList)
+ax.grid(True, alpha=0.2)
+ax.set_xlabel('Zenith (degrees)')
+ax.set_ylabel('Residual dispersion (mas)')
+plt.show()

scripts/NCEP NARR Data analysis/ReadMe.txt

@@ -0,0 +1,21 @@
+NARR data provided by the NOAA/OAR/ESRL PSL, Boulder, Colorado, USA, from their Web
+site at https://psl.noaa.gov/
+
+The original "air.201806.nc" datafile was downloaded from 
+https://psl.noaa.gov/data/gridded/data.narr.pressure.html
+
+Specific download link (~650 MB):
+ftp://ftp.cdc.noaa.gov/Datasets/NARR/pressure/air.201806.nc [Last accessed on 
+2020-06-17]
+
+
+The following files are a direct slice of this data:
+- export_dist.txt
+- export_heights.txt
+- export_press_atmosphere.txt
+- export_temp_atmosphere.txt
+
+
+This data file was not included here, because it has a very niche use for the
+AstroAtmosphere package and is a large file. The tradeoff seems to me to not very
+favourable, and therefore it is omitted.

scripts/NCEP NARR Data analysis/air_201806_description.txt

@@ -0,0 +1,97 @@
+NCEP North American Regional Reanalysis (NARR) gridded weather data.
+
+Data from https://www.esrl.noaa.gov/psd/data/gridded/data.narr.pressure.html
+
+
+
+###############
+# Header info # 
+###############
+$ ncdump -h air.201806.nc 
+netcdf air.201806 {
+dimensions:
+	time = UNLIMITED ; // (240 currently)
+	level = 29 ;
+	y = 277 ;
+	x = 349 ;
+variables:
+	double time(time) ;
+		time:axis = "T" ;
+		time:coordinate_defines = "point" ;
+		time:delta_t = "0000-00-00 03:00:00" ;
+		time:long_name = "Time" ;
+		time:standard_name = "time" ;
+		time:units = "hours since 1800-1-1 00:00:0.0" ;
+		time:actual_range = 1914576., 1915293. ;
+	float level(level) ;
+		level:GRIB_id = 100s ;
+		level:GRIB_name = "hPa" ;
+		level:actual_range = 1000.f, 100.f ;
+		level:axis = "Z" ;
+		level:coordinate_defines = "point" ;
+		level:long_name = "Level" ;
+		level:positive = "down" ;
+		level:standard_name = "level" ;
+		level:units = "millibar" ;
+	float lat(y, x) ;
+		lat:axis = "Y" ;
+		lat:coordinate_defines = "point" ;
+		lat:long_name = "Latitude" ;
+		lat:standard_name = "latitude" ;
+		lat:units = "degrees_north" ;
+	float lon(y, x) ;
+		lon:axis = "X" ;
+		lon:coordinate_defines = "point" ;
+		lon:long_name = "Longitude" ;
+		lon:standard_name = "longitude" ;
+		lon:units = "degrees_east" ;
+	float y(y) ;
+		y:long_name = "northward distance from southwest corner of domain in projection coordinates" ;
+		y:standard_name = "projection_y_coordinate" ;
+		y:units = "m" ;
+	float x(x) ;
+		x:long_name = "eastward distance from southwest corner of domain in projection coordinates" ;
+		x:standard_name = "projection_x_coordinate" ;
+		x:units = "m" ;
+	int Lambert_Conformal ;
+		Lambert_Conformal:false_easting = 5632642.22547 ;
+		Lambert_Conformal:false_northing = 4612545.65137 ;
+		Lambert_Conformal:grid_mapping_name = "lambert_conformal_conic" ;
+		Lambert_Conformal:latitude_of_projection_origin = 50. ;
+		Lambert_Conformal:longitude_of_central_meridian = -107. ;
+		Lambert_Conformal:standard_parallel = 50., 50. ;
+	float air(time, level, y, x) ;
+		air:GRIB_id = 11s ;
+		air:GRIB_name = "TMP" ;
+		air:_FillValue = 9.96921e+36f ;
+		air:coordinates = "lat lon" ;
+		air:dataset = "NARR 3-hourly" ;
+		air:grid_mapping = "Lambert_Conformal" ;
+		air:level_desc = "Pressure Levels" ;
+		air:long_name = "3-hourly Air Temperature on Pressure Levels" ;
+		air:parent_stat = "Other" ;
+		air:standard_name = "air_temperature" ;
+		air:statistic = "Individual Obs" ;
+		air:units = "degK" ;
+		air:valid_range = 137.5f, 362.5f ;
+		air:var_desc = "Air temperature" ;
+		air:missing_value = -9.96921e+36f ;
+		air:actual_range = 187.1279f, 320.1487f ;
+
+// global attributes:
+		:Conventions = "CF-1.2" ;
+		:centerlat = 50. ;
+		:centerlon = -107. ;
+		:comments = "" ;
+		:institution = "National Centers for Environmental Prediction" ;
+		:latcorners = 1.000001f, 0.897945f, 46.3544f, 46.63433f ;
+		:loncorners = -145.5f, -68.32005f, -2.569891f, 148.6418f ;
+		:platform = "Model" ;
+		:standardpar1 = 50. ;
+		:standardpar2 = 50.000001 ;
+		:title = "8x Daily NARR" ;
+		:history = "created Mon Jul 18 21:01:07 MDT 2016 by NOAA/ESRL/PSD" ;
+		:dataset_title = "NCEP North American Regional Reanalysis (NARR)" ;
+		:references = "https://www.esrl.noaa.gov/psd/data/gridded/data.narr.html" ;
+		:source = "http://www.emc.ncep.noaa.gov/mmb/rreanl/index.html" ;
+}

scripts/NCEP NARR Data analysis/export_dist.txt

@@ -0,0 +1 @@
+[0.0, 0.23179527183433432, 0.46359054366866864, 0.6953858155030029, 0.9271810873373373, 1.1589763591716715, 1.3907716310060059, 1.6225669028403402, 1.8543621746746746, 2.086157446509009, 2.317952718343343, 2.5497479901776776, 2.7815432620120117, 3.0133385338463463, 3.2451338056806804, 3.476929077515015, 3.708724349349349, 3.9405196211836833, 4.172314893018018, 4.404110164852352, 4.635905436686686, 4.867700708521021, 5.099495980355355, 5.33129125218969, 5.5630865240240235, 5.794881795858358, 6.026677067692693, 6.258472339527026, 6.490267611361361, 6.7220628831956954, 6.95385815503003, 7.185653426864364, 7.417448698698698, 7.649243970533033, 7.8810392423673665, 8.112834514201701, 8.344629786036036, 8.57642505787037, 8.808220329704705, 9.04001560153904, 9.271810873373372, 9.503606145207707, 9.735401417042041, 9.967196688876376, 10.19899196071071, 10.430787232545045, 10.66258250437938, 10.894377776213712, 11.126173048048047, 11.357968319882382, 11.589763591716716, 11.82155886355105, 12.053354135385385, 12.28514940721972, 12.516944679054053, 12.748739950888387, 12.980535222722722, 13.212330494557056, 13.444125766391391, 13.675921038225725, 13.90771631006006, 14.139511581894393, 14.371306853728727, 14.603102125563062, 14.834897397397397, 15.066692669231731, 15.298487941066066, 15.5302832129004, 15.762078484734733, 15.993873756569068, 16.225669028403402, 16.45746430023774, 16.68925957207207, 16.921054843906404, 17.15285011574074, 17.384645387575073, 17.61644065940941, 17.848235931243742, 18.08003120307808, 18.31182647491241, 18.543621746746744, 18.77541701858108, 19.007212290415413, 19.23900756224975, 19.470802834084083, 19.70259810591842, 19.93439337775275, 20.166188649587085, 20.39798392142142, 20.629779193255754, 20.86157446509009, 21.093369736924423, 21.32516500875876, 21.556960280593092, 21.788755552427425, 22.02055082426176, 22.252346096096094, 22.48414136793043, 22.715936639764763, 22.9477319115991, 23.179527183433432, 23.411322455267765, 23.6431177271021, 23.874912998936434, 24.10670827077077, 24.338503542605103, 24.57029881443944, 24.802094086273772, 25.033889358108105, 25.26568462994244, 25.497479901776774, 25.72927517361111, 25.961070445445444, 26.19286571727978, 26.424660989114113, 26.656456260948445, 26.888251532782782, 27.120046804617115, 27.35184207645145, 27.583637348285784, 27.81543262012012, 28.047227891954453, 28.279023163788786, 28.510818435623122, 28.742613707457455, 28.97440897929179, 29.206204251126124, 29.43799952296046, 29.669794794794793, 29.901590066629126, 30.133385338463462, 30.365180610297795, 30.59697588213213, 30.828771153966464, 31.0605664258008, 31.292361697635133, 31.524156969469466, 31.755952241303802, 31.987747513138135, 32.21954278497247, 32.451338056806804, 32.68313332864114, 32.91492860047548, 33.146723872309806, 33.37851914414414, 33.61031441597848, 33.84210968781281, 34.073904959647145, 34.30570023148148, 34.53749550331582, 34.76929077515015, 35.00108604698448, 35.23288131881882, 35.46467659065315, 35.696471862487485, 35.92826713432182, 36.16006240615616, 36.39185767799049, 36.62365294982482, 36.85544822165916, 37.08724349349349, 37.319038765327825, 37.55083403716216, 37.7826293089965, 38.01442458083083, 38.24621985266516, 38.4780151244995, 38.70981039633383, 38.941605668168165, 39.1734009400025, 39.40519621183684, 39.63699148367117, 39.8687867555055, 40.10058202733984, 40.33237729917417, 40.564172571008505, 40.79596784284284, 41.02776311467718, 41.25955838651151, 41.491353658345844, 41.72314893018018, 41.95494420201451, 42.186739473848846, 42.41853474568318, 42.65033001751752, 42.88212528935185, 43.113920561186184, 43.34571583302052, 43.57751110485485, 43.809306376689186, 44.04110164852352, 44.27289692035786, 44.50469219219219, 44.736487464026524, 44.96828273586086, 45.20007800769519, 45.431873279529526, 45.66366855136386, 45.8954638231982, 46.12725909503253, 46.359054366866864, 46.5908496387012, 46.82264491053553, 47.054440182369866, 47.2862354542042, 47.51803072603854, 47.74982599787287, 47.981621269707205, 48.21341654154154, 48.44521181337587, 48.67700708521021, 48.90880235704454, 49.14059762887888, 49.37239290071321, 49.604188172547545, 49.83598344438188, 50.06777871621621, 50.29957398805055, 50.53136925988488, 50.76316453171922, 50.99495980355355, 51.226755075387885, 51.45855034722222, 51.69034561905655, 51.92214089089089, 52.15393616272522, 52.38573143455956, 52.61752670639389, 52.849321978228225, 53.08111725006256, 53.31291252189689, 53.54470779373123, 53.776503065565564, 54.0082983373999, 54.24009360923423, 54.471888881068566, 54.7036841529029, 54.93547942473723, 55.16727469657157, 55.399069968405904, 55.63086524024024, 55.86266051207457, 56.094455783908906, 56.32625105574324, 56.55804632757757, 56.78984159941191, 57.021636871246244, 57.25343214308058, 57.48522741491491, 57.717022686749246, 57.94881795858358, 58.18061323041791, 58.41240850225225, 58.644203774086584, 58.87599904592092, 59.10779431775525, 59.339589589589586, 59.57138486142392, 59.80318013325825, 60Line truncated

scripts/NCEP NARR Data analysis/export_heights.txt

@@ -0,0 +1 @@
+[0.0, 0.21553884711779447, 0.43107769423558895, 0.6466165413533834, 0.8621553884711779, 1.0776942355889725, 1.2932330827067668, 1.5087719298245612, 1.7243107769423558, 1.9398496240601504, 2.155388471177945, 2.370927318295739, 2.5864661654135337, 2.8020050125313283, 3.0175438596491224, 3.233082706766917, 3.4486215538847116, 3.664160401002506, 3.8796992481203008, 4.095238095238095, 4.31077694235589, 4.526315789473684, 4.741854636591478, 4.957393483709273, 5.172932330827067, 5.3884711779448615, 5.6040100250626566, 5.819548872180451, 6.035087719298245, 6.25062656641604, 6.466165413533834, 6.681704260651629, 6.897243107769423, 7.112781954887217, 7.328320802005012, 7.5438596491228065, 7.7593984962406015, 7.974937343358396, 8.19047619047619, 8.406015037593985, 8.62155388471178, 8.837092731829573, 9.052631578947368, 9.268170426065163, 9.483709273182956, 9.699248120300751, 9.914786967418546, 10.13032581453634, 10.345864661654135, 10.56140350877193, 10.776942355889723, 10.992481203007518, 11.208020050125313, 11.423558897243106, 11.639097744360901, 11.854636591478696, 12.07017543859649, 12.285714285714285, 12.50125313283208, 12.716791979949875, 12.932330827067668, 13.147869674185463, 13.363408521303258, 13.578947368421051, 13.794486215538846, 14.010025062656641, 14.225563909774435, 14.44110275689223, 14.656641604010025, 14.872180451127818, 15.087719298245613, 15.303258145363408, 15.518796992481203, 15.734335839598996, 15.949874686716791, 16.165413533834585, 16.38095238095238, 16.596491228070175, 16.81203007518797, 17.027568922305765, 17.24310776942356, 17.45864661654135, 17.674185463659146, 17.88972431077694, 18.105263157894736, 18.32080200501253, 18.536340852130326, 18.751879699248118, 18.967418546365913, 19.182957393483708, 19.398496240601503, 19.614035087719298, 19.829573934837093, 20.045112781954884, 20.26065162907268, 20.476190476190474, 20.69172932330827, 20.907268170426065, 21.12280701754386, 21.338345864661655, 21.553884711779446, 21.76942355889724, 21.984962406015036, 22.20050125313283, 22.416040100250626, 22.63157894736842, 22.847117794486213, 23.062656641604008, 23.278195488721803, 23.493734335839598, 23.709273182957393, 23.924812030075188, 24.14035087719298, 24.355889724310774, 24.57142857142857, 24.786967418546364, 25.00250626566416, 25.218045112781954, 25.43358395989975, 25.64912280701754, 25.864661654135336, 26.08020050125313, 26.295739348370926, 26.51127819548872, 26.726817042606516, 26.942355889724308, 27.157894736842103, 27.373433583959898, 27.588972431077693, 27.804511278195488, 28.020050125313283, 28.235588972431078, 28.45112781954887, 28.666666666666664, 28.88220551378446, 29.097744360902254, 29.31328320802005, 29.528822055137844, 29.744360902255636, 29.95989974937343, 30.175438596491226, 30.39097744360902, 30.606516290726816, 30.82205513784461, 31.037593984962406, 31.253132832080198, 31.468671679197993, 31.684210526315788, 31.899749373433583, 32.11528822055138, 32.33082706766917, 32.54636591478697, 32.76190476190476, 32.97744360902256, 33.19298245614035, 33.40852130325814, 33.62406015037594, 33.83959899749373, 34.05513784461153, 34.27067669172932, 34.48621553884712, 34.70175438596491, 34.9172932330827, 35.1328320802005, 35.34837092731829, 35.56390977443609, 35.77944862155388, 35.994987468671674, 36.21052631578947, 36.426065162907264, 36.64160401002506, 36.857142857142854, 37.07268170426065, 37.288220551378444, 37.503759398496236, 37.719298245614034, 37.934837092731826, 38.150375939849624, 38.365914786967416, 38.581453634085214, 38.796992481203006, 39.0125313283208, 39.228070175438596, 39.44360902255639, 39.659147869674186, 39.87468671679198, 40.09022556390977, 40.30576441102757, 40.52130325814536, 40.73684210526316, 40.95238095238095, 41.16791979949875, 41.38345864661654, 41.59899749373433, 41.81453634085213, 42.03007518796992, 42.24561403508772, 42.46115288220551, 42.67669172932331, 42.8922305764411, 43.10776942355889, 43.32330827067669, 43.53884711779448, 43.75438596491228, 43.96992481203007, 44.185463659147864, 44.40100250626566, 44.616541353383454, 44.83208020050125, 45.047619047619044, 45.26315789473684, 45.478696741854634, 45.694235588972425, 45.909774436090224, 46.125313283208015, 46.340852130325814, 46.556390977443606, 46.771929824561404, 46.987468671679196, 47.20300751879699, 47.418546365914786, 47.63408521303258, 47.849624060150376, 48.06516290726817, 48.28070175438596, 48.49624060150376, 48.71177944862155, 48.92731829573935, 49.14285714285714, 49.35839598997494, 49.57393483709273, 49.78947368421052, 50.00501253132832, 50.22055137844611, 50.43609022556391, 50.6516290726817, 50.8671679197995, 51.08270676691729, 51.29824561403508, 51.51378446115288, 51.72932330827067, 51.94486215538847, 52.16040100250626, 52.37593984962406, 52.59147869674185, 52.807017543859644, 53.02255639097744, 53.238095238095234, 53.45363408521303, 53.669172932330824, 53.884711779448615, 54.100250626566414, 54.315789473684205, 54.531328320802004, 54.746867167919795, 54.962406015037594, 55.177944862155385, 55.39348370Line truncated

scripts/NCEP NARR Data analysis/export_press_atmosphere.txt

@@ -0,0 +1 @@
+[[100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.00000000001, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 99999.99999999999, 99999.99999999999, 100000.0, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.00000000001, 99999.99999999999, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.00000000001, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.00000000001, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.00000000001, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.00000000001, 100000.00000000001, 100000.00000000001, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.00000000001, 100000.0, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 99999.99999999999, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100000.0, 100Line truncated

scripts/NCEP NARR Data analysis/export_temp_atmosphere.txt

@@ -0,0 +1 @@
+[[284.050537109375, 284.05081584207153, 284.0510945747681, 284.0513733074647, 284.0516520401612, 284.0519307728577, 284.0522095055543, 284.05248823825076, 284.05276697094735, 284.05304570364393, 284.05332443634046, 284.053603169037, 284.0538819017336, 284.05416063443016, 284.0544393671267, 284.05471809982316, 284.05499683251975, 284.0552755652163, 284.0555542979128, 284.05583303060934, 284.0561117633059, 284.05639049600245, 284.05666922869904, 284.0569479613956, 284.05722669409215, 284.0575054267887, 284.0577841594852, 284.05806289218174, 284.05834162487827, 284.0586203575748, 284.0588990902713, 284.05917782296797, 284.0594565556645, 284.059735288361, 284.06001402105755, 284.06029275375414, 284.06057148645067, 284.06085021914726, 284.0611289518438, 284.0614076845403, 284.06168641723684, 284.06196514993337, 284.06224388262996, 284.0625226153265, 284.062801348023, 284.0630800807196, 284.06335881341613, 284.0636375461127, 284.06391627880924, 284.0641950115058, 284.06447374420236, 284.06475247689883, 284.0650312095954, 284.06530994229195, 284.06558867498853, 284.06586740768506, 284.0661461403816, 284.0664248730782, 284.0667036057747, 284.06698233847123, 284.0672610711678, 284.06753980386435, 284.0678185365609, 284.0680972692574, 284.06837600195394, 284.0686547346505, 284.06893346734705, 284.06921220004364, 284.06949093274017, 284.0697696654367, 284.0700483981333, 284.0703271308298, 284.0706058635264, 284.07088459622287, 284.0711633289194, 284.071442061616, 284.0717207943125, 284.0719995270091, 284.0722782597056, 284.07255699240216, 284.0728357250987, 284.07311445779527, 284.07339319049186, 284.0736719231884, 284.0739506558849, 284.07422938858144, 284.07450812127803, 284.0747868539745, 284.0750655866711, 284.0753443193677, 284.0756230520642, 284.07590178476073, 284.07618051745726, 284.0764592501538, 284.0767379828504, 284.0770167155469, 284.0772954482435, 284.07757418094, 284.07785291363655, 284.0781316463331, 284.07841037902966, 284.0786891117262, 284.0789678444228, 284.0792465771193, 284.07952530981584, 284.07980404251236, 284.08008277520895, 284.0803615079055, 284.080640240602, 284.0809189732986, 284.0811977059951, 284.08147643869165, 284.0817551713882, 284.08203390408477, 284.0823126367813, 284.0825913694779, 284.08287010217435, 284.08314883487094, 284.08342756756747, 284.08370630026405, 284.0839850329606, 284.0842637656571, 284.0845424983537, 284.0848212310502, 284.08509996374676, 284.08537869644334, 284.0856574291398, 284.0859361618364, 284.08621489453293, 284.0864936272295, 284.08677235992604, 284.0870510926226, 284.08732982531916, 284.0876085580157, 284.0878872907122, 284.0881660234088, 284.08844475610533, 284.08872348880186, 284.0890022214984, 284.089280954195, 284.0895596868915, 284.08983841958803, 284.0901171522846, 284.09039588498115, 284.0906746176777, 284.09095335037426, 284.0912320830708, 284.0915108157673, 284.09178954846385, 284.09206828116044, 284.09234701385697, 284.0926257465535, 284.0929044792501, 284.0931832119466, 284.09346194464314, 284.0937406773397, 284.09401941003625, 284.09429814273284, 284.0945768754293, 284.0948556081259, 284.0951343408224, 284.095413073519, 284.09569180621554, 284.09597053891207, 284.0962492716086, 284.0965280043052, 284.0968067370017, 284.0970854696983, 284.09736420239483, 284.106526646489, 284.1175328703515, 284.1285390942141, 284.13954531807667, 284.15055154193925, 284.1615577658018, 284.17256398966435, 284.1835702135269, 284.19457643738946, 284.205582661252, 284.21658888511456, 284.22759510897714, 284.23860133283966, 284.24960755670224, 284.2606137805648, 284.2716200044274, 284.2826262282899, 284.2936324521525, 284.3046386760151, 284.31564489987767, 284.3266511237402, 284.33765734760277, 284.3486635714653, 284.35966979532793, 284.37067601919046, 284.38168224305304, 284.39268846691556, 284.40369469077814, 284.4147009146407, 284.4257071385033, 284.4367133623658, 284.4477195862284, 284.458725810091, 284.4697320339535, 284.4807382578161, 284.49174448167867, 284.5027507055412, 284.51375692940377, 284.52476315326635, 284.53576937712893, 284.54677560099145, 284.55778182485403, 284.5687880487166, 284.5797942725792, 284.5908004964417, 284.6018067203043, 284.6128129441669, 284.6238191680294, 284.634825391892, 284.64583161575456, 284.6568378396171, 284.66784406347966, 284.67885028734224, 284.68985651120477, 284.70086273506735, 284.7118689589299, 284.7228751827925, 284.73388140665503, 284.7448876305176, 284.7558938543802, 284.76690007824277, 284.7779063021053, 284.7889125259679, 284.7999187498304, 284.810924973693, 284.82193119755556, 284.8329374214181, 284.84394364528066, 284.8549498691432, 284.86595609300576, 284.87696231686834, 284.8879685407309, 284.89897476459345, 284.909980988456, 284.9209872123186, 284.9319934361812, 284.9429996600437, 284.9540058839063, 284.96501210776887, 284.9760183316314, 284.987024555494, 284.99803077935655, 285.0090370032191, 285.0200432270817, 285.03104945094424, 285.04205567480676, 285.05306189866934, 285.0640681225319, 285.07507434639444, Line truncated

scripts/NCEP NARR Data analysis/ncep_atmosphere_analysis.py

@@ -0,0 +1,168 @@
+#!/usr/bin/env Python3
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.lines import Line2D
+from matplotlib import cm
+from matplotlib.colors import LogNorm
+from AstroAtmosphere import *
+from scipy import interpolate
+import json
+
+
+'''
+@ Author:		Joost van den Born 
+@ Contact:		born@astron.nl
+@ Description:	This script uses the NCEP NARR data to model a real structured atmosphere 
+				from data with a resolution of 2.5 km per data point. We interpolate this
+				data to a higher resolution.
+				Then a full integration over this atmosphere is done and compared to the 
+				values for the refraction and dispersion if we had only taken the atmospheric
+				conditions directly at the observatory.
+'''
+
+
+# Load a variable from a json text file.
+def loadVariableFromTXT(filename):
+	variableFile = open(filename, 'r')
+	variable = json.load(variableFile)
+	variableFile.close()
+	return np.asarray(variable)
+
+
+plt.rc('font', family='serif')
+plt.rc('axes', labelsize='large')
+
+# h, _, _, _, Temp, Press, Dens, vSound, _,_ = np.genfromtxt('si1py_modified.txt', unpack=True, skip_header=6)
+
+dist 	= loadVariableFromTXT('export_dist.txt')
+heights = loadVariableFromTXT('export_heights.txt')
+press 	= loadVariableFromTXT('export_press_atmosphere.txt')
+temp 	= loadVariableFromTXT('export_temp_atmosphere.txt')
+print("Loading complete!")
+
+# Conditions at Cerro Armazones
+T   = temp[0,0]	# k
+P   = press[0,0]# Pa
+H   = 0.0
+xc  = 314   	# ppm
+lat = -24.5983 	# degrees
+h0   = 0  		# m
+l1  = 1.49 		# micron
+l2  = 1.78 		# micron
+zenith = 60 	# degrees
+
+# Calculate the index of refraction
+at  = Observatory()
+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)
+rho = at.rho(p=P, T=T, RH=H, xc=xc)
+
+RHlist = np.asarray([H]*len(press[:,0]))
+xclist = np.asarray([xc]*len(press[:,0]))
+Dens = at.rho(p=press[:,0], T=temp[:,0], RH=RHlist, xc=xclist)
+
+refr = refraction(lat, h0) 
+refr.setReducedHeight(P, rho)
+disp = dispersion(lat, h0)
+disp.setReducedHeight(P, rho)
+
+
+# We must create a cube of weather data:
+# We assume a flat atmosphere for now, so that h does not need to be adjusted
+# x_max = 200 km, gives zenith_max = ~66.7 degrees
+# Every 20 km in horizontal direction a new perturbation and interpolation
+# for every 500 m.
+x_max = dist[-1] # km
+
+Temp_cube_n = []
+Press_cube_n = []
+
+for x in range(len(dist)):
+	Temp_cube_n.append(temp[:,0])
+	Press_cube_n.append(press[:,0])
+
+Temp_cube_n = np.asarray(Temp_cube_n).T
+Press_cube_n = np.asarray(Press_cube_n).T
+
+
+
+fig = plt.figure()
+ax1 = fig.add_subplot(211)
+im = ax1.pcolormesh(dist,heights, temp, cmap=cm.coolwarm)
+fig.colorbar(im, ax=ax1, label='Temperature (K)')
+ax1.set_aspect(1)
+ax1.set_ylabel('Height (km)')
+
+ax2 = fig.add_subplot(212)
+im2 = ax2.pcolormesh(dist, heights, press, cmap=cm.coolwarm, norm=LogNorm())
+ax2.set_aspect(1)
+ax2.set_xlabel('Horizontal distance (km)')
+ax2.set_ylabel('Height (km)')
+fig.colorbar(im2, ax=ax2, label='Pressure (Pa)')
+plt.tight_layout()
+plt.show(block=False)
+
+# Now for a given zenith angle we need to find the temp and press along the line of sight
+f_Temp_2d  		= interpolate.interp2d(dist, heights, temp, kind='linear')
+f_Press_2d 		= interpolate.interp2d(dist, heights, press, kind='linear')
+f_Temp_2d_n  	= interpolate.interp2d(dist, heights, Temp_cube_n, kind='linear')
+f_Press_2d_n 	= interpolate.interp2d(dist, heights, Press_cube_n, kind='linear')
+print('Interpolation complete!')
+
+
+
+zenithRange = np.linspace(30,65,10)
+refrDiff, dispDiff = [], []
+
+for zenith in zenithRange:
+	z = np.deg2rad(zenith)
+	y_los = heights
+	x_los = np.linspace(0, heights.max()*np.tan(z), len(heights))
+
+	Temp_los, Press_los = [], []
+	Temp_los_n, Press_los_n = [], []
+	for i in range(len(heights)):
+		 Temp_los.append(f_Temp_2d(x_los[i], y_los[i])[0])
+		 Press_los.append(f_Press_2d(x_los[i], y_los[i])[0])
+		 Temp_los_n.append(f_Temp_2d_n(x_los[i], y_los[i])[0])
+		 Press_los_n.append(f_Press_2d_n(x_los[i], y_los[i])[0])
+
+	Temp_los  	= np.asarray(Temp_los)
+	Press_los 	= np.asarray(Press_los)
+	Temp_los_n 	= np.asarray(Temp_los_n)
+	Press_los_n = np.asarray(Press_los_n)
+	rho_los   	= at.rho(p=Press_los, T=Temp_los, RH=H, xc=xc)
+	rho_los_n 	= at.rho(p=Press_los_n, T=Temp_los_n, RH=H, xc=xc)
+
+	refrInt = refr.refractionIntegral(n1, z, useStandardAtmosphere=False, heightData=heights, rhoData=rho_los_n)*3600*1000
+	refrInt_pert = refr.refractionIntegral(n1, z, useStandardAtmosphere=False, heightData=heights, rhoData=rho_los)*3600*1000
+	refrDiff.append(refrInt_pert - refrInt)
+
+	dispInt = disp.refractionIntegral(n1, n2, z, useStandardAtmosphere=False, heightData=heights, rhoData=rho_los_n)*3.6e6
+	dispInt_pert = disp.refractionIntegral(n1, n2, z, useStandardAtmosphere=False, heightData=heights, rhoData=rho_los)*3.6e6
+	dispDiff.append(dispInt_pert - dispInt)
+
+
+
+
+
+fig3 = plt.figure(figsize=(9,7))
+ax1 = fig3.add_subplot(211)
+ax2 = fig3.add_subplot(212)
+ax1.tick_params(axis ='both', direction='in', which='both', labelsize=16) 
+ax2.tick_params(axis ='both', direction='in', which='both', labelsize=16) 
+ax1.grid(True, alpha=0.3)#, which='major')
+ax2.grid(True, alpha=0.3)#, which='major')
+
+
+ax1.plot(zenithRange, np.log10(np.abs(refrDiff)), 'k-', linewidth=2)
+ax2.plot(zenithRange, np.log10(np.abs(dispDiff)), 'k-', linewidth=2)
+ax1.set_xlabel('Zenith (degrees)', fontsize=18)
+ax2.set_xlabel('Zenith (degrees)', fontsize=18)
+ax1.set_ylabel('$\log_{10}($|$R - R_{o}$|$)$ (mas)', fontsize=18)
+ax2.set_ylabel('$\log_{10}($|$\Delta R - \Delta R_{o}$|$)$ (mas)', fontsize=18)
+plt.tight_layout()
+plt.show()
+
+
+	

scripts/NCEP NARR Data analysis/ncep_data_analysis.py

@@ -0,0 +1,251 @@
+#!/usr/bin/env Python3
+import numpy as np
+import sys
+from AstroAtmosphere import *
+import matplotlib.pyplot as plt
+from netCDF4 import Dataset
+from pdb import set_trace
+from matplotlib import cm
+from scipy.interpolate import interp1d, interp2d
+from matplotlib.colors import LogNorm
+import json
+
+# Saves a variable, preferably a list or string to json text file.
+# Useful for long simulations, where you're making a plot
+def saveVariableToTXT(variable, filename):
+	variableFile = open(filename, 'w')
+
+	if type(variable) == np.ndarray:
+		json.dump(variable.tolist(), variableFile)
+	elif type(variable) in [list, str, float]:
+		json.dump(variable, variableFile)
+	else:
+		print('Not a list, array or string. Check variable type.')
+	
+	variableFile.close()
+	return
+
+# Load a variable from a json text file.
+def loadVariableFromTXT(filename):
+	variableFile = open(filename, 'r')
+	variable = json.load(variableFile)
+	variableFile.close()
+	return variable
+
+
+def calculateDistanceResolution(lat, lon):
+	# Calculate local radius of the earth
+	_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)
+
+	# Calculate circumference of the earth at constant latitude.
+	latR = 2 * np.pi * np.cos(_phi) * _rc0
+	distRes = latR/360 * lon
+	return distRes
+
+def H_isotherm(T, lat):
+	'''
+		Calculates the reduced height of the atmosphere, assuming an isotherm distribution.
+	'''
+	_kb = 1.38064852e-23 # kg m2 s-2 K-1
+	_m  = 4.809651698e-26 # kg (weight per molecule of air, assuming dry air)
+
+	_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)
+	return (_kb*T) / (_m * _g0_local)
+
+
+def calculatePressure(Th, P0, h, lat):
+	redH = H_isotherm(Th, lat)
+	Ph = P0 * np.exp(-h/redH)
+	return Ph	
+
+plt.rc('font', family='serif')
+plt.rc('axes', labelsize='large')
+
+# Import the real data
+nc = Dataset('air.201806.nc', 'r')
+latitudes  = np.asarray(nc.variables['lat'][:])
+longitudes = np.asarray(nc.variables['lon'][:])
+air_temp   = nc.variables['air']
+air_press  = np.asarray(nc.variables['level'])
+
+
+
+
+at = Observatory()
+
+# Parameters at Cerro Armazones
+T   = air_temp[0,:,0,0][0]    # k
+P   = air_press[0]     # Pa
+H   = 0
+xc  = 450       # ppm
+lat = 24.44768
+h   = 0      	# m
+l1  = 1.49      # micron
+l2  = 1.78      # micron
+
+rho = at.rho(p=P, T=T, RH=H, xc=xc)
+disp = dispersion(lat, h)
+disp.setReducedHeight(P, rho)
+
+air_dens = at.rho(p=air_press, T=air_temp[0,:,0,0], RH=np.asarray([0]*len(air_press)), xc=np.asarray([0]*len(air_press)))
+
+heights = disp.calcHeight(air_press, air_temp[0,:,0,0], RH=np.asarray([0]*len(air_press)), xc=np.asarray([0]*len(air_press)), P0=air_press[0]) / 1e3
+distance_lat = calculateDistanceResolution(latitudes[100], longitudes[100,:]) - calculateDistanceResolution(latitudes[100,0], longitudes[100,0])
+temp_selection_ground = air_temp[0,0,100,:]
+temp_selection_upper = air_temp[0,-1,100,:]
+print(distance_lat[1])
+
+
+#########################################################
+# Extrapolating using the U.S. Standard Atmosphere 1976.
+# The end point of the real data is used to shift the USSA data accordingly (linearly)
+h_ussa1976, _, _, _, air_temp_ussa1976, air_press_ussa1976,_,_,_,_ = np.genfromtxt('../US1976 Standard Atmosphere/si1py_modified.txt', unpack=True, skip_header=6)
+temp_ip 	= interp1d(h_ussa1976, air_temp_ussa1976)
+press_ip 	= interp1d(h_ussa1976, air_press_ussa1976)
+
+ussa1976_start 		= temp_ip(heights[-1])
+diff_temp_list 		= temp_selection_upper - ussa1976_start
+start_i 			= np.argmax(h_ussa1976 > heights[-1])
+ussa1976_temp_array = np.tile(air_temp_ussa1976[start_i:], (len(temp_selection_upper),1))
+ussa1976_temp_array = np.asarray([ussa1976_temp_array[i] + diff_temp_list[i] for i in range(len(diff_temp_list))]).T
+
+# Combining the (masked) arrays
+new_heights 	= np.append(heights, h_ussa1976[start_i:])
+tempData 		= np.append(air_temp[0,:,100,:], ussa1976_temp_array, axis=0)
+tempData_mask 	= np.append(air_temp[0,:,100,:].mask, np.zeros(ussa1976_temp_array.shape), axis=0)
+tempData 		= np.ma.masked_array(tempData, mask=tempData_mask)
+
+# Now the same for the pressure data
+# Combining the (masked) arrays
+start_i 		= np.argmax(h_ussa1976 > heights[-1])
+new_heights 	= np.append(heights, h_ussa1976[start_i:])
+pressData = []
+for i in range(len(temp_selection_upper)):
+	calc_press_data = calculatePressure(Th=ussa1976_temp_array[:,i], P0=air_press[0], h=h_ussa1976[start_i:]*1000, lat=lat)
+	pressData_i		= np.append(air_press*100, calc_press_data*100, axis=0)
+	pressData.append(pressData_i)
+pressData = np.asarray(pressData).T
+#########################################################
+
+# We'd now like to export a subset of the atmospheric data for a distance of ~200 km.
+# Let's use interpolation to create a uniform gridded dataset
+export_heights 			= np.linspace(0,new_heights.max(),400)
+export_dist 			= np.linspace(0,distance_lat[136]-distance_lat[130],1000)
+atmosphere_temp_func 	= interp2d(distance_lat[130:136]-distance_lat[130], new_heights, tempData[:,130:136], kind='linear')
+export_temp_atmosphere 	= atmosphere_temp_func(export_dist, export_heights)
+
+atmosphere_press_func 	= interp2d(distance_lat[130:136]-distance_lat[130], new_heights, pressData[:,130:136], kind='linear')
+export_press_atmosphere = atmosphere_press_func(export_dist, export_heights)
+
+# saveVariableToTXT(export_heights, 'export_heights.txt')
+# saveVariableToTXT(export_dist, 'export_dist.txt')
+# saveVariableToTXT(export_temp_atmosphere, 'export_temp_atmosphere.txt')
+# saveVariableToTXT(export_press_atmosphere, 'export_press_atmosphere.txt')
+
+
+# Plot the real data + extrapolated standard atmosphere.
+# fig = plt.figure()
+# ax  = fig.add_subplot(211)
+# ax2 = fig.add_subplot(212)
+# im1 = ax.pcolormesh(export_dist, export_heights, export_temp_atmosphere, cmap=cm.coolwarm)
+# im2 = ax2.pcolormesh(export_dist, export_heights, export_press_atmosphere, cmap=cm.coolwarm, norm=LogNorm())
+# ax2.set_xlabel('Lateral distance (km)')
+# ax.set_ylabel('Height (km)')
+# ax2.set_ylabel('Height (km)')
+# cbar = fig.colorbar(im1, ax=ax, label='Temperature (K)')
+# lbls = cbar.get_ticks()
+# cbar.ax.invert_yaxis() 
+# cbar.set_ticks(lbls[::-1])
+# cbar2 = fig.colorbar(im2, ax=ax2, label='Pressure (Pa)')
+# lbls = cbar2.get_ticks()
+# cbar2.ax.invert_yaxis() 
+# cbar2.set_ticks(lbls[::-1])
+# plt.tight_layout()
+# plt.show(block=False)
+# set_trace()
+
+
+
+
+# # Plot the real data + extrapolated standard atmosphere.
+# fig = plt.figure()
+# ax  = fig.add_subplot(211)
+# ax2 = fig.add_subplot(212)
+# im1 = ax.imshow(tempData[:,130:137], cmap=cm.coolwarm, interpolation= 'bicubic', aspect='auto', origin='lower')
+# # im1 = ax.pcolormesh(distance_lat[130:137], new_heights, tempData[:,130:137], cmap=cm.coolwarm)
+# im2 = ax2.pcolormesh(distance_lat, new_heights, pressData, cmap=cm.coolwarm, norm=LogNorm())
+# ax2.set_xlabel('Lateral distance (km)')
+# ax.set_ylabel('Height (km)')
+# ax2.set_ylabel('Height (km)')
+# cbar = fig.colorbar(im1, ax=ax, label='Temperature (K)')
+# lbls = cbar.get_ticks()
+# cbar.ax.invert_yaxis() 
+# cbar.set_ticks(lbls[::-1])
+# cbar2 = fig.colorbar(im2, ax=ax2, label='Pressure (Pa)')
+# lbls = cbar2.get_ticks()
+# cbar2.ax.invert_yaxis() 
+# cbar2.set_ticks(lbls[::-1])
+# plt.tight_layout()
+# plt.show(block=False)
+
+
+# Plot the real data and properties
+fig = plt.figure(figsize=(9,4))
+fig2 = plt.figure()
+ax1 = fig2.add_subplot(221)
+ax2 = fig2.add_subplot(222)
+ax3 = fig2.add_subplot(223)
+ax4 = fig.add_subplot(121)
+ax5 = fig2.add_subplot(224)
+ax6 = fig.add_subplot(122)
+
+ax1.grid(alpha=0.3)
+ax2.grid(alpha=0.3)
+ax3.grid(alpha=0.3)
+ax5.grid(alpha=0.3)
+
+ax1.plot(air_temp[0,:,0,0], heights)
+ax1.set_xlabel('Temperature (K)')
+ax1.set_ylabel('Height (km)')
+
+ax2.plot(air_press, heights)
+ax2.set_xlabel('Pressure (Pa)')
+ax2.set_ylabel('Height (km)')
+
+ax3.plot(air_dens, heights)
+ax3.set_xlabel('Density (g cm$^{-3}$)')
+ax3.set_ylabel('Height (km)')
+
+im = ax4.imshow(np.rot90(air_temp[0,0,:,:].T), aspect=1, cmap='jet')
+ax4.set_yticks([])
+ax4.set_xticks([])
+# fig.colorbar(im, ax=ax4, label='Temperature (K)')
+
+ax5.plot(distance_lat, temp_selection_ground, 'b-', label='ground')
+ax5.plot(distance_lat, temp_selection_upper, 'r-', label='upper atmosphere')
+ax5.legend()
+ax5.set_ylabel('Temperature (K)')
+ax5.set_xlabel('Distance (km)')
+
+im2 = ax6.imshow(air_temp[0,:,100,:], aspect='auto', origin='lower', interpolation='bilinear', extent=[0,distance_lat[-1],0,heights[-1] ], cmap='jet')
+ax6.set_xlabel('Distance (km)', fontsize=18)
+ax6.set_ylabel('Height (km)', fontsize=18)
+cbar = fig.colorbar(im2, ax=ax6)#, label='Temperature (K)')
+cbar.ax.tick_params(labelsize=16)
+cbar.set_label('Temperature (K)', size=18)
+lbls = cbar.get_ticks()
+cbar.ax.invert_yaxis() 
+cbar.set_ticks(lbls[::-1])
+ax6.tick_params(axis ='both', direction='in', labelsize=16) 
+
+plt.tight_layout()
+plt.show()

scripts/USSA 1976/create_USSA_tables.py

@@ -0,0 +1,232 @@
+#!/usr/bin/env python
+
+""" Tables.py - Make tables  of atmospheric properties   (Python 3)
+
+Adapted by
+    Richard J. Kwan, Lightsaber Computing
+from original programs by
+    Ralph L. Carmichael, Public Domain Aeronautical Software
+
+Revision History
+Date         Vers Person Statement of Changes
+2004 Oct 04  1.0  RJK    Initial program
+2017 Jun 04  1.1  RLC    All indents are with spaces; all prints in ( )
+                         New version for Python 3 that does integer div with //
+"""
+
+import sys, math
+
+version = "1.1 (2017 Jun 04)"
+greeting = "Tables - A Python program to compute atmosphere tables"
+author = "Ralph L. Carmichael, Public Domain Aeronautical Software"
+modifier = "Richard J. Kwan, Lightsaber Computing"
+farewell = "Four files added to your directory."
+finalmess = "Normal termination of tables."
+
+#   P H Y S I C A L   C O N S T A N T S
+
+FT2METERS = 0.3048      # mult. ft. to get meters (exact)
+KELVIN2RANKINE = 1.8    # mult deg K to get deg R
+PSF2NSM = 47.880258     # mult lb/sq.ft to get sq.m
+SCF2KCM = 515.379       # mult slugs/cu.ft to get kg/cu.m
+TZERO   = 288.15        # sea-level temperature, kelvins
+PZERO   = 101325.0      # sea-level pressure, N/sq.m
+RHOZERO = 1.225         # sea-level density, kg/cu.m
+AZERO   = 340.294       # speed of sound at S.L.  m/sec
+BETAVISC = 1.458E-6     # viscosity constant
+SUTH    = 110.4         # Sutherland's constant, kelvins
+
+def LongUSTable():
+    Itxt = open('us1py.prt', 'w')
+    Itxt.write(' alt   sigma     delta    theta ')
+    Itxt.write(' temp   press     dens       a    visc  k.visc\n')
+    Itxt.write(' Kft                            ')
+    Itxt.write(' degR lb/sq.ft  s/cu.ft     fps s/ft-s sq.ft/s\n')
+
+    for i in range(-1, 57):
+        altKm=5*i*FT2METERS
+        (sigma, delta, theta) = Atmosphere(altKm)
+        Itxt.write("%4i %9.3E %9.3E %6.4f " % 
+            (5*i, sigma, delta, theta))
+        temp=(KELVIN2RANKINE*TZERO)*theta
+        pressure=(PZERO/PSF2NSM)*delta
+        density=(RHOZERO/SCF2KCM)*sigma
+        asound=(AZERO/FT2METERS)*math.sqrt(theta)
+        Itxt.write("%5.1f %8.3E %9.3E %6.1f" %
+            (temp, pressure, density, asound))
+        viscosity=(1.0/PSF2NSM)*MetricViscosity(theta)
+        kinematicViscosity=viscosity/density
+        Itxt.write("%6.3f %8.2E\n" %
+            (1.0E6*viscosity, kinematicViscosity))
+
+    Itxt.close()
+
+def ShortUSTable():
+    Itxt = open('us2py.prt', 'w')
+    Itxt.write(' alt  sigma  delta  theta ')
+    Itxt.write(' temp  press    dens     a     visc   k.visc ratio\n')
+    Itxt.write(' Kft                      ')
+    Itxt.write(' degR   psf   s/cu.ft   fps s/ft-sec sq.ft/s   1/ft\n')
+
+    for i in range(-1, 66):
+        altKm=i*FT2METERS
+        (sigma, delta, theta) = SimpleAtmosphere(altKm)
+        Itxt.write("%4i %6.4f %6.4f %6.4f" % 
+            (i, sigma, delta, theta))
+        temp=KELVIN2RANKINE*TZERO*theta
+        pressure=PZERO*delta/47.88
+        density=RHOZERO*sigma/515.379
+        asound=(AZERO/FT2METERS)*math.sqrt(theta);
+        Itxt.write("%6.1f %6.1f %9.7f %6.1f" %
+            (temp, pressure, density, asound))
+        viscosity=(1.0/PSF2NSM)*MetricViscosity(theta)
+        kinematicViscosity=viscosity/density
+        vratio=asound/kinematicViscosity
+        Itxt.write("%6.3f %8.2E %4.2f\n"%
+            (1.0E6*viscosity, kinematicViscosity, 1.0E-6*vratio))
+            
+    Itxt.close()
+
+def LongSITable():
+    Itxt = open('si1py.prt', 'w')
+    Itxt.write(' alt    sigma      delta    theta ')
+    Itxt.write(' temp   press      dens     a   visc   k.visc\n')
+    Itxt.write('  Km                              ')
+    Itxt.write('   K   N/sq.m    kg/cu.m  m/sec kg/m-s sq.m/s\n')
+
+    for i in range(-4,173):
+        altKm=0.5*i
+        (sigma, delta, theta) = Atmosphere(altKm)
+        Itxt.write("%4.1f %8.4E %8.4E %6.4f" %
+            (0.5*i, sigma,delta,theta))
+        temp=TZERO*theta
+        pressure=PZERO*delta
+        density=RHOZERO*sigma
+        asound=AZERO*math.sqrt(theta)
+        Itxt.write("%6.1f %8.5E %8.3E %5.1f" %
+            (temp,pressure,density,asound))
+        viscosity=MetricViscosity(theta)
+        kinematicViscosity=viscosity/density
+        Itxt.write("%6.2f %8.2E\n" %
+            (1.0E6*viscosity, kinematicViscosity))
+
+    Itxt.close()
+
+def ShortSITable():
+    Itxt = open('si2py.prt', 'w')
+    Itxt.write(' alt  sigma  delta  theta ')
+    Itxt.write(' temp  press  dens   a    visc   k.visc ratio\n')
+    Itxt.write('  Km                      ')
+    Itxt.write('   K  N/sq.m  kcm  m/sec kg/m-s  sq.m/s  1/m\n')
+
+    for i in range(-1, 41):
+        altKm=0.5*i
+        (sigma, delta, theta) = SimpleAtmosphere(altKm)
+        Itxt.write("%4.1f %6.4f %6.4f %6.4f" %
+            (altKm, sigma, delta, theta))
+        temp=TZERO*theta
+        pressure=PZERO*delta
+        density=RHOZERO*sigma
+        asound=AZERO*math.sqrt(theta)
+        Itxt.write("%6.1f %6.0f %5.3f %5.1f" %
+            (temp, pressure,density,asound))
+        viscosity=MetricViscosity(theta)
+        kinematicViscosity=viscosity/density
+        vratio=asound/kinematicViscosity
+        Itxt.write("%6.2f %8.2E %5.2f\n" %
+            (1.0E6*viscosity, kinematicViscosity, 1.0E-6*vratio))
+
+    Itxt.close()
+
+def Atmosphere(alt):
+    """ Compute temperature, density, and pressure in standard atmosphere.
+    Correct to 86 km.  Only approximate thereafter.
+    Input:
+    alt geometric altitude, km.
+    Return: (sigma, delta, theta)
+    sigma   density/sea-level standard density
+    delta   pressure/sea-level standard pressure
+    theta   temperature/sea-level std. temperature
+    """
+
+    REARTH = 6369.0     # radius of the Earth (km)
+    GMR = 34.163195
+    NTAB = 8            # length of tables
+
+    htab = [ 0.0,  11.0, 20.0, 32.0, 47.0,
+    51.0, 71.0, 84.852 ]
+    ttab = [ 288.15, 216.65, 216.65, 228.65, 270.65,
+    270.65, 214.65, 186.946 ]
+    ptab = [ 1.0, 2.2336110E-1, 5.4032950E-2, 8.5666784E-3, 1.0945601E-3,
+    6.6063531E-4, 3.9046834E-5, 3.68501E-6 ]
+    gtab = [ -6.5, 0.0, 1.0, 2.8, 0, -2.8, -2.0, 0.0 ]
+
+    h = alt*REARTH/(alt+REARTH) # geometric to geopotential altitude
+
+    i=0; j=len(htab)-1
+    while (j > i+1):
+        k = (i+j)//2      # this is floor division in Python 3
+        if h < htab[k]:
+            j = k
+        else:
+            i = k
+    tgrad = gtab[i]     # temp. gradient of local layer
+    tbase = ttab[i]     # base  temp. of local layer
+    deltah=h-htab[i]        # height above local base
+    tlocal=tbase+tgrad*deltah   # local temperature
+    theta = tlocal/ttab[0]  # temperature ratio
+
+    if 0.0 == tgrad:
+        delta=ptab[i]*math.exp(-GMR*deltah/tbase)
+    else:
+        delta=ptab[i]*math.pow(tbase/tlocal, GMR/tgrad)
+    sigma = delta/theta
+    return ( sigma, delta, theta )
+
+def SimpleAtmosphere(alt):
+    """ Compute temperature, density, and pressure in simplified
+    standard atmosphere.
+
+    Correct to 20 km.  Only approximate thereafter.
+
+    Input:
+        alt geometric altitude, km.
+    Return: (sigma, delta, theta)
+        sigma   density/sea-level standard density
+        delta   pressure/sea-level standard pressure
+        theta   temperature/sea-level std. temperature
+    """
+    REARTH = 6369.0     # radius of the Earth (km)
+    GMR = 34.163195     # gas constant
+
+    h = alt*REARTH/(alt+REARTH) # geometric to geopotential altitude
+
+    if h<11.0:      # troposphere
+        theta=(288.15-6.5*h)/288.15
+        delta=math.pow(theta, GMR/6.5)
+    else:       # stratosphere
+        theta=216.65/288.15
+        delta=0.2233611*math.exp(-GMR*(h-11.0)/216.65)
+    sigma = delta/theta
+    return ( sigma, delta, theta )
+
+def MetricViscosity(theta):
+    t=theta*TZERO
+    return BETAVISC*math.sqrt(t*t*t)/(t+SUTH)
+
+def main():
+    print ("Executing ", sys.argv[0])
+    print (greeting)
+    print (author)
+    if modifier != "":
+        print ("Modified by ", modifier)
+        print ("    version ", version)
+    LongUSTable()
+    ShortUSTable()
+    LongSITable()
+    ShortSITable()
+
+    print (farewell)
+    print (finalmess)
+
+main()

scripts/USSA 1976/perturbed_atmoshere.py

@@ -0,0 +1,242 @@
+#!/usr/bin/env Python3
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.lines import Line2D
+from matplotlib import cm
+from matplotlib.colors import LogNorm
+import sys
+from AstroAtmosphere import *
+from scipy import interpolate
+
+'''
+@ Author:		Joost van den Born 
+@ Contact:		born@astron.nl
+@ Description:	This script uses the US Standard Atmosphere and perturbs it 
+				at several points. Then full integration is used to calculate
+				the atmospheric dispersion. This is then compared to when we
+				would have only used the atmospheric conditions directly at 
+				the observatory. This is done to assess whether we should 
+				monitor the weather at every height along the line of sight of
+				the telescope.
+'''
+
+def perturbAtmosphere(h, Temp_perturbed, Press_perturbed):
+	for _i in range(len(h)):
+		if h[_i] >= 1.0 and h[_i] <= 2.0:
+			Tnoise = np.random.normal(0,5,1)
+			Temp_perturbed[_i] = Temp_perturbed[_i] + Tnoise
+			Pnoise = np.random.normal(0,1000,1)
+			Press_perturbed[_i] = Press_perturbed[_i] + Pnoise
+		if h[_i] >= 7.0 and h[_i] <= 10.0:
+			Tnoise = np.random.normal(0,5,1)
+			Temp_perturbed[_i] = Temp_perturbed[_i] + Tnoise
+			Pnoise = np.random.normal(0,500,1)
+			Press_perturbed[_i] = Press_perturbed[_i] + Pnoise
+		if h[_i] >= 25.0 and h[_i] <= 30.0:
+			Tnoise = np.random.normal(0,10,1)
+			Temp_perturbed[_i] = Temp_perturbed[_i] + Tnoise
+			Pnoise = np.random.normal(0,250,1)
+			Press_perturbed[_i] = Press_perturbed[_i] + Pnoise
+		if h[_i] >= 50.0 and h[_i] <= 60.0:
+			Tnoise = np.random.normal(0,15,1)
+			Temp_perturbed[_i] = Temp_perturbed[_i] + Tnoise
+			Pnoise = np.random.normal(0,4,1)
+			Press_perturbed[_i] = Press_perturbed[_i] + Pnoise
+	return Temp_perturbed, Press_perturbed
+
+plt.rc('font', family='serif')
+plt.rc('axes', labelsize='large')
+
+h, _, _, _, Temp, Press, Dens, vSound, _,_ = np.genfromtxt('si1py_modified.txt', unpack=True, skip_header=6)
+
+# Conditions at Cerro Armazones
+T   = 288.15 	# k
+P   = 101324  	# Pa
+H   = 0.0
+xc  = 314   	# ppm
+lat = -24.5983 	# degrees
+h0   = 0  	# m
+l1  = 1.49 		# micron
+l2  = 1.78 		# micron
+zenith = 60 	# degrees
+
+# Calculate the index of refraction
+at  = Observatory()
+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)
+rho = at.rho(p=P, T=T, RH=H, xc=xc)
+
+refr = refraction(lat, h0) 
+refr.setReducedHeight(P, rho)
+disp = dispersion(lat, h0)
+disp.setReducedHeight(P, rho)
+
+
+
+# We must create a cube of weather data:
+# We assume a flat atmosphere for now, so that h does not need to be adjusted
+# x_max = 200 km, gives zenith_max = ~66.7 degrees
+# Every 20 km in horizontal direction a new perturbation and interpolation
+# for every 500 m.
+
+x_max = 200 # km
+
+Temp_cube  = []
+Press_cube = []
+Temp_cube_n = []
+Press_cube_n = []
+Temp_cube.append(Temp)
+Press_cube.append(Press)
+Temp_cube_n.append(Temp)
+Press_cube_n.append(Press)
+
+
+
+for x in range(int(x_max/(2))):
+	# Adding noise
+	# 1-2 km: 2 degree noise, 100 Pa noise
+	# 7-10 km: 3 degree noise, 150 Pa noise
+	# 25-30 km: 5 degree noise, 250 Pa noise
+	Temp_perturbed = Temp.copy() 
+	Press_perturbed = Press.copy()
+
+	Temp_perturbed, Press_perturbed = perturbAtmosphere(h, Temp_perturbed, Press_perturbed)
+	Temp_cube.append(Temp_perturbed)
+	Press_cube.append(Press_perturbed)
+	Temp_cube_n.append(Temp)
+	Press_cube_n.append(Press)
+
+
+f_Temp_n  	= interpolate.interp1d(np.linspace(0,x_max, int(x_max/2)+1), Temp_cube_n, axis=0)
+f_Press_n 	= interpolate.interp1d(np.linspace(0,x_max, int(x_max/2)+1), Press_cube_n, axis=0)
+f_Temp  	= interpolate.interp1d(np.linspace(0,x_max, int(x_max/2)+1), Temp_cube, axis=0)
+f_Press 	= interpolate.interp1d(np.linspace(0,x_max, int(x_max/2)+1), Press_cube, axis=0)
+
+Temp_interpolated_n 	= f_Temp_n(np.linspace(0,x_max, int(x_max/0.5)+1) )
+Press_interpolated_n 	= f_Press_n(np.linspace(0,x_max, int(x_max/0.5)+1) )
+Temp_interpolated 		= f_Temp(np.linspace(0,x_max, int(x_max/0.5)+1) )
+Press_interpolated 		= f_Press(np.linspace(0,x_max, int(x_max/0.5)+1) )
+
+x, y = np.mgrid[0:200:(int(x_max/0.5)+1)*1j, 0:np.max(h):len(h)*1j]
+
+fig = plt.figure()
+ax1 = fig.add_subplot(211)
+im = ax1.pcolor(x,y, Temp_interpolated, cmap=cm.coolwarm)
+fig.colorbar(im, ax=ax1, label='Temperature (K)')
+ax1.set_aspect(1)
+ax1.set_ylabel('Height (km)')
+
+ax2 = fig.add_subplot(212)
+im2 = ax2.pcolor(x,y, Press_interpolated, cmap=cm.coolwarm, norm=LogNorm())
+ax2.set_aspect(1)
+ax2.set_xlabel('Horizontal distance (km)')
+ax2.set_ylabel('Height (km)')
+fig.colorbar(im2, ax=ax2, label='Pressure (Pa)')
+plt.tight_layout()
+plt.savefig('2d_atmosphere.png', bbox_inches='tight')
+plt.savefig('2d_atmosphere.pdf', bbox_inches='tight')
+plt.show(block=False)
+
+
+# Now for a given zenith angle we need to find the temp and press along the line of sight
+x = np.linspace(0,200,(int(x_max/0.5)+1))
+y = h
+
+f_Temp_2d_n  = interpolate.interp2d(y,x,Temp_interpolated_n, kind='linear')
+f_Press_2d_n = interpolate.interp2d(y,x,Press_interpolated_n, kind='linear')
+f_Temp_2d  	 = interpolate.interp2d(y,x,Temp_interpolated, kind='linear')
+f_Press_2d 	 = interpolate.interp2d(y,x,Press_interpolated, kind='linear')
+
+fig2 = plt.figure()
+ax = fig2.add_subplot(121)
+ax.grid(True, alpha=0.3)
+ax.plot(Temp, h, 'k-')
+
+ax2 = fig2.add_subplot(122)
+ax2.grid(True, alpha=0.3)
+ax2.plot(Press, h, 'k-')
+
+
+zenithRange = np.linspace(30,65,100)
+refrDiff, dispDiff = [], []
+
+for zenith in zenithRange:
+	z = np.deg2rad(zenith)
+	y_los = np.linspace(0, 86, len(h))
+	x_los = np.linspace(0, 86*np.tan(z), len(h))
+
+	Temp_los, Press_los = [], []
+	Temp_los_n, Press_los_n = [], []
+	for i in range(len(h)):
+		 Temp_los.append(f_Temp_2d(y_los[i], x_los[i])[0])
+		 Press_los.append(f_Press_2d(y_los[i], x_los[i])[0])
+		 Temp_los_n.append(f_Temp_2d_n(y_los[i], x_los[i])[0])
+		 Press_los_n.append(f_Press_2d_n(y_los[i], x_los[i])[0])
+
+	Temp_los  	= np.asarray(Temp_los)
+	Press_los 	= np.asarray(Press_los)
+	Temp_los_n 	= np.asarray(Temp_los_n)
+	Press_los_n = np.asarray(Press_los_n)
+	rho_los   	= at.rho(p=Press_los, T=Temp_los, RH=H, xc=xc)
+	rho_los_n 	= at.rho(p=Press_los_n, T=Temp_los_n, RH=H, xc=xc)
+
+	ax.plot(Temp_los, h, label='z=%i' %zenith)
+	ax2.plot(Press_los, h, label='z=%i' %zenith)
+
+	refrInt = refr.refractionIntegral(n1, z, useStandardAtmosphere=False, heightData=h, rhoData=rho_los_n)*3600*1000
+	refrInt_pert = refr.refractionIntegral(n1, z, useStandardAtmosphere=False, heightData=h, rhoData=rho_los)*3600*1000
+	refrDiff.append(refrInt_pert - refrInt)
+
+	dispInt = disp.refractionIntegral(n1, n2, z, useStandardAtmosphere=False, heightData=h, rhoData=rho_los_n)*3.6e6
+	dispInt_pert = disp.refractionIntegral(n1, n2, z, useStandardAtmosphere=False, heightData=h, rhoData=rho_los)*3.6e6
+	dispDiff.append(dispInt_pert - dispInt)
+
+ax.legend()
+ax2.legend()
+
+
+
+# fig3 = plt.figure()
+# ax1 = fig3.add_subplot(211)
+# ax2 = fig3.add_subplot(212)
+# ax1.tick_params(axis ='both', direction='in') 
+# ax2.tick_params(axis ='both', direction='in') 
+# ax1.grid(True, alpha=0.3)
+# ax2.grid(True, alpha=0.3)
+# ax1.semilogy(zenithRange, np.abs(refrDiff), 'k-')
+# ax2.semilogy(zenithRange, np.abs(dispDiff), 'k-')
+# ax1.set_xlabel('Zenith (degrees)')
+# ax2.set_xlabel('Zenith (degrees)')
+# ax1.set_ylabel('|$R_{p} - R_{SA}$| (mas)')
+# ax2.set_ylabel('|$\Delta R_{p} - \Delta R_{SA}$| (mas)')
+# plt.tight_layout()
+# plt.savefig('differential_USSA_for_different_z.png', bbox_inches='tight')
+# plt.savefig('differential_USSA_for_different_z.pdf', bbox_inches='tight')
+# plt.show()
+
+fig3 = plt.figure(figsize=(9,7))
+ax1 = fig3.add_subplot(211)
+ax2 = fig3.add_subplot(212)
+ax1.tick_params(axis ='both', direction='in', which='both', labelsize=16) 
+ax2.tick_params(axis ='both', direction='in', which='both', labelsize=16) 
+ax1.grid(True, alpha=0.3)
+ax2.grid(True, alpha=0.3)
+
+
+ax1.plot(zenithRange, np.log10(np.abs(refrDiff)), 'k-', linewidth=2)
+ax2.plot(zenithRange, np.log10(np.abs(dispDiff)), 'k-', linewidth=2)
+ax1.set_xlabel('Zenith (degrees)', fontsize=18)
+ax2.set_xlabel('Zenith (degrees)', fontsize=18)
+ax1.set_ylabel('$\log_{10}($|$R - R_{o}$|$)$ (mas)', fontsize=18)
+ax2.set_ylabel('$\log_{10}($|$\Delta R - \Delta R_{o}$|$)$ (mas)', fontsize=18)
+plt.tight_layout()
+# plt.savefig('differential_USSA_for_different_z.png', bbox_inches='tight')
+# plt.savefig('differential_USSA_for_different_z.pdf', bbox_inches='tight')
+plt.show()
+
+
+
+
+# Comparison of different air columns.
+
+	

scripts/USSA 1976/refraction_height.py

@@ -0,0 +1,91 @@
+#!/usr/bin/env Python3
+import numpy as np
+from pdb import set_trace
+import matplotlib.pyplot as plt
+from matplotlib.lines import Line2D
+import sys
+from AstroAtmosphere import *
+
+'''
+@ Author:		Joost van den Born 
+@ Contact:		born@astron.nl
+@ Description:	We plot the total refraction as a function of height.
+				We find that most refraction happens in the lower atmosphere
+				i.e. <10 km.
+'''
+
+def refractionIntegral2(n, R, z0, n0, R0):
+	z0 = np.deg2rad(z0)
+	return (n0 * R0 * np.sin(z0))/(n*np.sqrt((n*R)**2 - (n0 * R0 * np.sin(z0))**2 ))
+
+def n_for_heights(n0, z0, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None):
+	'''
+		Using data that gives temperature, pressure and density as a function of height, 
+		the gladstone-dale relation is invoked to calculate the path of a light ray 
+		along the atmosphere. This is allows us to also include atmospheric activity in 
+		different atmosphere layers.
+		Recommended use with the US Standard Atmosphere from 1976.
+		USSA1976 data generated from code supplied by http://www.pdas.com/atmos.html
+	'''
+	if useStandardAtmosphere:
+		_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('si1py_modified.txt', unpack=True, skip_header=6)
+	else:
+		_Dens = rhoData
+		_h 	  = heightData
+	z0 		= np.deg2rad(z0)
+	_Dens 	= np.asarray(_Dens)
+	_h 	 	= np.asarray(_h) 
+	_R 		= R0 + _h * 1000
+	
+	# Gladstone-Dale relation between refractive index and density
+	_n_h = 1 + (n0 - 1) / _Dens[0] * _Dens
+	_z_h = np.arcsin( (n0 * R0) / (_n_h * _R) * np.sin(z0) )
+
+	return _h, _n_h, np.degrees(_z_h)
+
+
+h, _, _, _, Temp, Press, Dens, vSound, _,_ = np.genfromtxt('si1py_modified.txt', unpack=True, skip_header=6)
+
+
+# Conditions
+T   = 288.15 	# k
+P   = 101324  	# Pa
+H   = 0.0
+xc  = 314   	# ppm
+lat = -24.5983 	# degrees
+h0   = 0  	# m
+l1  = 1.49 		# micron
+l2  = 1.78 		# micron
+zenith = 60 	# degrees
+
+# Calculate the index of refraction
+at  = Observatory()
+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)
+rho = at.rho(p=P, T=T, RH=H, xc=xc)
+
+refr = refraction(lat, h0) 
+refr.setReducedHeight(P, rho)
+disp = dispersion(lat, h0)
+disp.setReducedHeight(P, rho)
+
+
+fig = plt.figure(figsize=(5,4))
+ax = fig.add_subplot(111)
+ax.grid(True, alpha=0.3)
+ax.set_xlabel('Height (km)')
+ax.set_ylabel('Normalized cumulative refraction')
+ax.set_xlim([0, h[-1]])
+
+for z in [30]:
+	refList = []
+	dispList = []
+	for i in range(len(h)):
+		refList.append(refr.refractionIntegral(n1, z, useStandardAtmosphere=False, heightData=h[:i+1], rhoData=Dens[:i+1])*3600)
+		dispList.append(disp.refractionIntegral(n1, n2, z, useStandardAtmosphere=False, heightData=h[:i+1], rhoData=Dens[:i+1])*3.6e6)
+
+	ax.plot(h, np.asarray(refList)/refList[-1], 'k-', label='$\lambda=1.49$ $\mu$m, z=%i$^{\\circ}$' %z)	
+ax.set_ylim([0,1.1])
+
+plt.tight_layout()
+plt.show()

scripts/USSA 1976/si1py_modified.txt

@@ -0,0 +1,179 @@
+ alt    sigma      delta    theta  temp   press      dens     a   visc   k.visc
+  Km                                 K   N/sq.m    kg/cu.m  m/sec kg/m-s sq.m/s
+-2.0 1.2067E+00 1.2611E+00 1.0451 301.2 1.27783E+05 1.478E+00 347.9 18.51 1.25E-05
+-1.5 1.1522E+00 1.1912E+00 1.0338 297.9 1.20696E+05 1.411E+00 346.0 18.36 1.30E-05
+-1.0 1.0996E+00 1.1244E+00 1.0226 294.7 1.13931E+05 1.347E+00 344.1 18.21 1.35E-05
+-0.5 1.0489E+00 1.0607E+00 1.0113 291.4 1.07478E+05 1.285E+00 342.2 18.05 1.40E-05
+ 0.0 1.0000E+00 1.0000E+00 1.0000 288.1 1.01325E+05 1.225E+00 340.3 17.89 1.46E-05
+ 0.5 9.5288E-01 9.4213E-01 0.9887 284.9 9.54613E+04 1.167E+00 338.4 17.74 1.52E-05
+ 1.0 9.0748E-01 8.8701E-01 0.9774 281.7 8.98763E+04 1.112E+00 336.4 17.58 1.58E-05
+ 1.5 8.6376E-01 8.3454E-01 0.9662 278.4 8.45597E+04 1.058E+00 334.5 17.42 1.65E-05
+ 2.0 8.2168E-01 7.8462E-01 0.9549 275.2 7.95014E+04 1.007E+00 332.5 17.26 1.71E-05
+ 2.5 7.8119E-01 7.3715E-01 0.9436 271.9 7.46917E+04 9.570E-01 330.6 17.10 1.79E-05
+ 3.0 7.4225E-01 6.9204E-01 0.9324 268.7 7.01211E+04 9.093E-01 328.6 16.94 1.86E-05
+ 3.5 7.0482E-01 6.4920E-01 0.9211 265.4 6.57804E+04 8.634E-01 326.6 16.78 1.94E-05
+ 4.0 6.6885E-01 6.0854E-01 0.9098 262.2 6.16604E+04 8.193E-01 324.6 16.61 2.03E-05
+ 4.5 6.3432E-01 5.6997E-01 0.8986 258.9 5.77525E+04 7.770E-01 322.6 16.45 2.12E-05
+ 5.0 6.0117E-01 5.3341E-01 0.8873 255.7 5.40482E+04 7.364E-01 320.5 16.28 2.21E-05
+ 5.5 5.6936E-01 4.9878E-01 0.8760 252.4 5.05393E+04 6.975E-01 318.5 16.12 2.31E-05
+ 6.0 5.3887E-01 4.6600E-01 0.8648 249.2 4.72176E+04 6.601E-01 316.5 15.95 2.42E-05
+ 6.5 5.0964E-01 4.3499E-01 0.8535 245.9 4.40754E+04 6.243E-01 314.4 15.78 2.53E-05
+ 7.0 4.8165E-01 4.0568E-01 0.8423 242.7 4.11052E+04 5.900E-01 312.3 15.61 2.65E-05
+ 7.5 4.5485E-01 3.7799E-01 0.8310 239.5 3.82996E+04 5.572E-01 310.2 15.44 2.77E-05
+ 8.0 4.2921E-01 3.5185E-01 0.8198 236.2 3.56515E+04 5.258E-01 308.1 15.27 2.90E-05
+ 8.5 4.0470E-01 3.2721E-01 0.8085 233.0 3.31541E+04 4.958E-01 306.0 15.10 3.05E-05
+ 9.0 3.8128E-01 3.0398E-01 0.7973 229.7 3.08006E+04 4.671E-01 303.8 14.93 3.20E-05
+ 9.5 3.5891E-01 2.8211E-01 0.7860 226.5 2.85846E+04 4.397E-01 301.7 14.75 3.36E-05
+10.0 3.3756E-01 2.6153E-01 0.7748 223.3 2.64998E+04 4.135E-01 299.5 14.58 3.53E-05
+10.5 3.1720E-01 2.4219E-01 0.7635 220.0 2.45401E+04 3.886E-01 297.4 14.40 3.71E-05
+11.0 2.9780E-01 2.2403E-01 0.7523 216.8 2.26998E+04 3.648E-01 295.2 14.22 3.90E-05
+11.5 2.7545E-01 2.0710E-01 0.7519 216.7 2.09847E+04 3.374E-01 295.1 14.22 4.21E-05
+12.0 2.5464E-01 1.9146E-01 0.7519 216.7 1.93993E+04 3.119E-01 295.1 14.22 4.56E-05
+12.5 2.3541E-01 1.7699E-01 0.7519 216.7 1.79340E+04 2.884E-01 295.1 14.22 4.93E-05
+13.0 2.1763E-01 1.6363E-01 0.7519 216.7 1.65795E+04 2.666E-01 295.1 14.22 5.33E-05
+13.5 2.0119E-01 1.5127E-01 0.7519 216.7 1.53275E+04 2.465E-01 295.1 14.22 5.77E-05
+14.0 1.8600E-01 1.3985E-01 0.7519 216.7 1.41702E+04 2.279E-01 295.1 14.22 6.24E-05
+14.5 1.7196E-01 1.2929E-01 0.7519 216.7 1.31005E+04 2.107E-01 295.1 14.22 6.75E-05
+15.0 1.5898E-01 1.1953E-01 0.7519 216.7 1.21117E+04 1.948E-01 295.1 14.22 7.30E-05
+15.5 1.4698E-01 1.1051E-01 0.7519 216.7 1.11976E+04 1.801E-01 295.1 14.22 7.90E-05
+16.0 1.3589E-01 1.0217E-01 0.7519 216.7 1.03527E+04 1.665E-01 295.1 14.22 8.54E-05
+16.5 1.2564E-01 9.4465E-02 0.7519 216.7 9.57164E+03 1.539E-01 295.1 14.22 9.24E-05
+17.0 1.1616E-01 8.7339E-02 0.7519 216.7 8.84961E+03 1.423E-01 295.1 14.22 9.99E-05
+17.5 1.0740E-01 8.0752E-02 0.7519 216.7 8.18215E+03 1.316E-01 295.1 14.22 1.08E-04
+18.0 9.9302E-02 7.4662E-02 0.7519 216.7 7.56512E+03 1.216E-01 295.1 14.22 1.17E-04
+18.5 9.1815E-02 6.9032E-02 0.7519 216.7 6.99471E+03 1.125E-01 295.1 14.22 1.26E-04
+19.0 8.4893E-02 6.3828E-02 0.7519 216.7 6.46738E+03 1.040E-01 295.1 14.22 1.37E-04
+19.5 7.8494E-02 5.9017E-02 0.7519 216.7 5.97989E+03 9.616E-02 295.1 14.22 1.48E-04
+20.0 7.2578E-02 5.4569E-02 0.7519 216.7 5.52921E+03 8.891E-02 295.1 14.22 1.60E-04
+20.5 6.6979E-02 5.0460E-02 0.7534 217.1 5.11291E+03 8.205E-02 295.4 14.24 1.74E-04
+21.0 6.1807E-02 4.6670E-02 0.7551 217.6 4.72884E+03 7.571E-02 295.7 14.27 1.88E-04
+21.5 5.7045E-02 4.3173E-02 0.7568 218.1 4.37446E+03 6.988E-02 296.0 14.29 2.05E-04
+22.0 5.2660E-02 3.9945E-02 0.7585 218.6 4.04741E+03 6.451E-02 296.4 14.32 2.22E-04
+22.5 4.8622E-02 3.6965E-02 0.7603 219.1 3.74551E+03 5.956E-02 296.7 14.35 2.41E-04
+23.0 4.4902E-02 3.4214E-02 0.7620 219.6 3.46678E+03 5.500E-02 297.0 14.38 2.61E-04
+23.5 4.1474E-02 3.1674E-02 0.7637 220.1 3.20939E+03 5.081E-02 297.4 14.40 2.83E-04
+24.0 3.8316E-02 2.9328E-02 0.7654 220.6 2.97167E+03 4.694E-02 297.7 14.43 3.07E-04
+24.5 3.5404E-02 2.7161E-02 0.7672 221.1 2.75206E+03 4.337E-02 298.1 14.46 3.33E-04
+25.0 3.2721E-02 2.5158E-02 0.7689 221.6 2.54915E+03 4.008E-02 298.4 14.48 3.61E-04
+25.5 3.0246E-02 2.3308E-02 0.7706 222.0 2.36163E+03 3.705E-02 298.7 14.51 3.92E-04
+26.0 2.7964E-02 2.1597E-02 0.7723 222.5 2.18831E+03 3.426E-02 299.1 14.54 4.24E-04
+26.5 2.5858E-02 2.0016E-02 0.7740 223.0 2.02807E+03 3.168E-02 299.4 14.57 4.60E-04
+27.0 2.3916E-02 1.8553E-02 0.7758 223.5 1.87991E+03 2.930E-02 299.7 14.59 4.98E-04
+27.5 2.2124E-02 1.7201E-02 0.7775 224.0 1.74289E+03 2.710E-02 300.1 14.62 5.39E-04
+28.0 2.0470E-02 1.5950E-02 0.7792 224.5 1.61614E+03 2.508E-02 300.4 14.65 5.84E-04
+28.5 1.8943E-02 1.4793E-02 0.7809 225.0 1.49888E+03 2.320E-02 300.7 14.67 6.32E-04
+29.0 1.7533E-02 1.3722E-02 0.7826 225.5 1.39037E+03 2.148E-02 301.0 14.70 6.84E-04
+29.5 1.6231E-02 1.2731E-02 0.7844 226.0 1.28994E+03 1.988E-02 301.4 14.73 7.41E-04
+30.0 1.5028E-02 1.1813E-02 0.7861 226.5 1.19698E+03 1.841E-02 301.7 14.75 8.01E-04
+30.5 1.3917E-02 1.0964E-02 0.7878 227.0 1.11092E+03 1.705E-02 302.0 14.78 8.67E-04
+31.0 1.2891E-02 1.0177E-02 0.7895 227.5 1.03122E+03 1.579E-02 302.4 14.81 9.38E-04
+31.5 1.1942E-02 9.4488E-03 0.7912 228.0 9.57402E+02 1.463E-02 302.7 14.83 1.01E-03
+32.0 1.1065E-02 8.7740E-03 0.7930 228.5 8.89024E+02 1.355E-02 303.0 14.86 1.10E-03
+32.5 1.0228E-02 8.1493E-03 0.7968 229.6 8.25726E+02 1.253E-02 303.8 14.92 1.19E-03
+33.0 9.4469E-03 7.5724E-03 0.8016 231.0 7.67272E+02 1.157E-02 304.7 14.99 1.30E-03
+33.5 8.7297E-03 7.0395E-03 0.8064 232.4 7.13279E+02 1.069E-02 305.6 15.07 1.41E-03
+34.0 8.0709E-03 6.5470E-03 0.8112 233.7 6.63380E+02 9.887E-03 306.5 15.14 1.53E-03
+34.5 7.4654E-03 6.0917E-03 0.8160 235.1 6.17243E+02 9.145E-03 307.4 15.21 1.66E-03
+35.0 6.9085E-03 5.6705E-03 0.8208 236.5 5.74564E+02 8.463E-03 308.3 15.29 1.81E-03
+35.5 6.3961E-03 5.2807E-03 0.8256 237.9 5.35066E+02 7.835E-03 309.2 15.36 1.96E-03
+36.0 5.9245E-03 4.9198E-03 0.8304 239.3 4.98495E+02 7.257E-03 310.1 15.43 2.13E-03
+36.5 5.4901E-03 4.5854E-03 0.8352 240.7 4.64618E+02 6.725E-03 311.0 15.51 2.31E-03
+37.0 5.0899E-03 4.2756E-03 0.8400 242.1 4.33223E+02 6.235E-03 311.9 15.58 2.50E-03
+37.5 4.7209E-03 3.9883E-03 0.8448 243.4 4.04115E+02 5.783E-03 312.8 15.65 2.71E-03
+38.0 4.3806E-03 3.7218E-03 0.8496 244.8 3.77116E+02 5.366E-03 313.7 15.72 2.93E-03
+38.5 4.0666E-03 3.4746E-03 0.8544 246.2 3.52061E+02 4.982E-03 314.6 15.79 3.17E-03
+39.0 3.7767E-03 3.2450E-03 0.8592 247.6 3.28801E+02 4.626E-03 315.4 15.87 3.43E-03
+39.5 3.5089E-03 3.0318E-03 0.8640 249.0 3.07198E+02 4.298E-03 316.3 15.94 3.71E-03
+40.0 3.2616E-03 2.8337E-03 0.8688 250.4 2.87125E+02 3.995E-03 317.2 16.01 4.01E-03
+40.5 3.0329E-03 2.6496E-03 0.8736 251.7 2.68467E+02 3.715E-03 318.1 16.08 4.33E-03
+41.0 2.8214E-03 2.4783E-03 0.8784 253.1 2.51116E+02 3.456E-03 318.9 16.15 4.67E-03
+41.5 2.6257E-03 2.3190E-03 0.8832 254.5 2.34974E+02 3.216E-03 319.8 16.22 5.04E-03
+42.0 2.4445E-03 2.1708E-03 0.8880 255.9 2.19952E+02 2.995E-03 320.7 16.29 5.44E-03
+42.5 2.2768E-03 2.0327E-03 0.8928 257.3 2.05965E+02 2.789E-03 321.5 16.36 5.87E-03
+43.0 2.1214E-03 1.9041E-03 0.8976 258.6 1.92937E+02 2.599E-03 322.4 16.43 6.32E-03
+43.5 1.9774E-03 1.7843E-03 0.9024 260.0 1.80799E+02 2.422E-03 323.3 16.50 6.81E-03
+44.0 1.8438E-03 1.6727E-03 0.9072 261.4 1.69484E+02 2.259E-03 324.1 16.57 7.34E-03
+44.5 1.7199E-03 1.5685E-03 0.9120 262.8 1.58932E+02 2.107E-03 325.0 16.64 7.90E-03
+45.0 1.6050E-03 1.4714E-03 0.9168 264.2 1.49090E+02 1.966E-03 325.8 16.71 8.50E-03
+45.5 1.4983E-03 1.3808E-03 0.9216 265.5 1.39905E+02 1.835E-03 326.7 16.78 9.14E-03
+46.0 1.3992E-03 1.2961E-03 0.9263 266.9 1.31330E+02 1.714E-03 327.5 16.85 9.83E-03
+46.5 1.3071E-03 1.2171E-03 0.9311 268.3 1.23322E+02 1.601E-03 328.4 16.92 1.06E-02
+47.0 1.2215E-03 1.1433E-03 0.9359 269.7 1.15841E+02 1.496E-03 329.2 16.99 1.14E-02
+47.5 1.1437E-03 1.0743E-03 0.9393 270.6 1.08849E+02 1.401E-03 329.8 17.04 1.22E-02
+48.0 1.0748E-03 1.0095E-03 0.9393 270.6 1.02287E+02 1.317E-03 329.8 17.04 1.29E-02
+48.5 1.0100E-03 9.4865E-04 0.9393 270.6 9.61215E+01 1.237E-03 329.8 17.04 1.38E-02
+49.0 9.4912E-04 8.9147E-04 0.9393 270.6 9.03286E+01 1.163E-03 329.8 17.04 1.47E-02
+49.5 8.9193E-04 8.3776E-04 0.9393 270.6 8.48858E+01 1.093E-03 329.8 17.04 1.56E-02
+50.0 8.3819E-04 7.8728E-04 0.9393 270.6 7.97716E+01 1.027E-03 329.8 17.04 1.66E-02
+50.5 7.8770E-04 7.3986E-04 0.9393 270.6 7.49663E+01 9.649E-04 329.8 17.04 1.77E-02
+51.0 7.4026E-04 6.9530E-04 0.9393 270.6 7.04511E+01 9.068E-04 329.8 17.04 1.88E-02
+51.5 6.9630E-04 6.5342E-04 0.9384 270.4 6.62082E+01 8.530E-04 329.7 17.02 2.00E-02
+52.0 6.5758E-04 6.1395E-04 0.9336 269.0 6.22084E+01 8.055E-04 328.8 16.96 2.10E-02
+52.5 6.2085E-04 5.7668E-04 0.9289 267.7 5.84321E+01 7.605E-04 328.0 16.89 2.22E-02
+53.0 5.8599E-04 5.4150E-04 0.9241 266.3 5.48678E+01 7.178E-04 327.1 16.82 2.34E-02
+53.5 5.5293E-04 5.0831E-04 0.9193 264.9 5.15047E+01 6.773E-04 326.3 16.75 2.47E-02
+54.0 5.2158E-04 4.7700E-04 0.9145 263.5 4.83322E+01 6.389E-04 325.4 16.68 2.61E-02
+54.5 4.9187E-04 4.4748E-04 0.9098 262.1 4.53405E+01 6.025E-04 324.6 16.61 2.76E-02
+55.0 4.6370E-04 4.1964E-04 0.9050 260.8 4.25201E+01 5.680E-04 323.7 16.54 2.91E-02
+55.5 4.3702E-04 3.9341E-04 0.9002 259.4 3.98620E+01 5.354E-04 322.9 16.47 3.08E-02
+56.0 4.1175E-04 3.6869E-04 0.8954 258.0 3.73576E+01 5.044E-04 322.0 16.40 3.25E-02
+56.5 3.8782E-04 3.4541E-04 0.8907 256.6 3.49988E+01 4.751E-04 321.2 16.33 3.44E-02
+57.0 3.6516E-04 3.2349E-04 0.8859 255.3 3.27778E+01 4.473E-04 320.3 16.26 3.64E-02
+57.5 3.4373E-04 3.0286E-04 0.8811 253.9 3.06871E+01 4.211E-04 319.4 16.19 3.85E-02
+58.0 3.2344E-04 2.8344E-04 0.8763 252.5 2.87198E+01 3.962E-04 318.6 16.12 4.07E-02
+58.5 3.0426E-04 2.6518E-04 0.8716 251.1 2.68692E+01 3.727E-04 317.7 16.05 4.31E-02
+59.0 2.8612E-04 2.4800E-04 0.8668 249.8 2.51289E+01 3.505E-04 316.8 15.98 4.56E-02
+59.5 2.6897E-04 2.3186E-04 0.8620 248.4 2.34929E+01 3.295E-04 315.9 15.91 4.83E-02
+60.0 2.5276E-04 2.1668E-04 0.8573 247.0 2.19554E+01 3.096E-04 315.1 15.84 5.11E-02
+60.5 2.3746E-04 2.0243E-04 0.8525 245.6 2.05110E+01 2.909E-04 314.2 15.77 5.42E-02
+61.0 2.2300E-04 1.8904E-04 0.8477 244.3 1.91546E+01 2.732E-04 313.3 15.69 5.75E-02
+61.5 2.0935E-04 1.7647E-04 0.8430 242.9 1.78811E+01 2.565E-04 312.4 15.62 6.09E-02
+62.0 1.9647E-04 1.6468E-04 0.8382 241.5 1.66860E+01 2.407E-04 311.5 15.55 6.46E-02
+62.5 1.8431E-04 1.5361E-04 0.8334 240.2 1.55648E+01 2.258E-04 310.7 15.48 6.86E-02
+63.0 1.7285E-04 1.4323E-04 0.8287 238.8 1.45132E+01 2.117E-04 309.8 15.41 7.28E-02
+63.5 1.6204E-04 1.3351E-04 0.8239 237.4 1.35275E+01 1.985E-04 308.9 15.33 7.72E-02
+64.0 1.5185E-04 1.2439E-04 0.8191 236.0 1.26036E+01 1.860E-04 308.0 15.26 8.20E-02
+64.5 1.4225E-04 1.1585E-04 0.8144 234.7 1.17382E+01 1.743E-04 307.1 15.19 8.72E-02
+65.0 1.3321E-04 1.0785E-04 0.8096 233.3 1.09277E+01 1.632E-04 306.2 15.12 9.26E-02
+65.5 1.2470E-04 1.0036E-04 0.8048 231.9 1.01690E+01 1.528E-04 305.3 15.04 9.85E-02
+66.0 1.1668E-04 9.3354E-05 0.8001 230.5 9.45909E+00 1.429E-04 304.4 14.97 1.05E-01
+66.5 1.0914E-04 8.6800E-05 0.7953 229.2 8.79502E+00 1.337E-04 303.5 14.90 1.11E-01
+67.0 1.0204E-04 8.0672E-05 0.7906 227.8 8.17408E+00 1.250E-04 302.6 14.82 1.19E-01
+67.5 9.5371E-05 7.4944E-05 0.7858 226.4 7.59372E+00 1.168E-04 301.7 14.75 1.26E-01
+68.0 8.9101E-05 6.9593E-05 0.7811 225.1 7.05149E+00 1.091E-04 300.7 14.67 1.34E-01
+68.5 8.3209E-05 6.4595E-05 0.7763 223.7 6.54509E+00 1.019E-04 299.8 14.60 1.43E-01
+69.0 7.7675E-05 5.9929E-05 0.7715 222.3 6.07234E+00 9.515E-05 298.9 14.53 1.53E-01
+69.5 7.2478E-05 5.5576E-05 0.7668 221.0 5.63120E+00 8.879E-05 298.0 14.45 1.63E-01
+70.0 6.7601E-05 5.1515E-05 0.7620 219.6 5.21972E+00 8.281E-05 297.1 14.38 1.74E-01
+70.5 6.3026E-05 4.7728E-05 0.7573 218.2 4.83607E+00 7.721E-05 296.1 14.30 1.85E-01
+71.0 5.8734E-05 4.4199E-05 0.7525 216.8 4.47851E+00 7.195E-05 295.2 14.23 1.98E-01
+71.5 5.4712E-05 4.0912E-05 0.7478 215.5 4.14543E+00 6.702E-05 294.3 14.15 2.11E-01
+72.0 5.0905E-05 3.7852E-05 0.7436 214.3 3.83531E+00 6.236E-05 293.4 14.08 2.26E-01
+72.5 4.7296E-05 3.5007E-05 0.7402 213.3 3.54709E+00 5.794E-05 292.8 14.03 2.42E-01
+73.0 4.3927E-05 3.2365E-05 0.7368 212.3 3.27939E+00 5.381E-05 292.1 13.98 2.60E-01
+73.5 4.0786E-05 2.9912E-05 0.7334 211.3 3.03083E+00 4.996E-05 291.4 13.92 2.79E-01
+74.0 3.7856E-05 2.7635E-05 0.7300 210.3 2.80012E+00 4.637E-05 290.7 13.87 2.99E-01
+74.5 3.5125E-05 2.5522E-05 0.7266 209.4 2.58605E+00 4.303E-05 290.1 13.81 3.21E-01
+75.0 3.2580E-05 2.3563E-05 0.7232 208.4 2.38749E+00 3.991E-05 289.4 13.76 3.45E-01
+75.5 3.0209E-05 2.1746E-05 0.7198 207.4 2.20338E+00 3.701E-05 288.7 13.70 3.70E-01
+76.0 2.8001E-05 2.0061E-05 0.7164 206.4 2.03272E+00 3.430E-05 288.0 13.65 3.98E-01
+76.5 2.5946E-05 1.8501E-05 0.7131 205.5 1.87458E+00 3.178E-05 287.4 13.59 4.28E-01
+77.0 2.4033E-05 1.7055E-05 0.7097 204.5 1.72810E+00 2.944E-05 286.7 13.54 4.60E-01
+77.5 2.2253E-05 1.5716E-05 0.7063 203.5 1.59247E+00 2.726E-05 286.0 13.48 4.95E-01
+78.0 2.0597E-05 1.4477E-05 0.7029 202.5 1.46692E+00 2.523E-05 285.3 13.43 5.32E-01
+78.5 1.9058E-05 1.3331E-05 0.6995 201.6 1.35076E+00 2.335E-05 284.6 13.37 5.73E-01
+79.0 1.7627E-05 1.2270E-05 0.6961 200.6 1.24331E+00 2.159E-05 283.9 13.32 6.17E-01
+79.5 1.6298E-05 1.1290E-05 0.6927 199.6 1.14396E+00 1.996E-05 283.2 13.26 6.64E-01
+80.0 1.5063E-05 1.0384E-05 0.6893 198.6 1.05213E+00 1.845E-05 282.5 13.21 7.16E-01
+80.5 1.3917E-05 9.5464E-06 0.6860 197.7 9.67291E-01 1.705E-05 281.8 13.15 7.71E-01
+81.0 1.2853E-05 8.7731E-06 0.6826 196.7 8.88933E-01 1.574E-05 281.1 13.10 8.32E-01
+81.5 1.1866E-05 8.0591E-06 0.6792 195.7 8.16590E-01 1.454E-05 280.4 13.04 8.97E-01
+82.0 1.0950E-05 7.4002E-06 0.6758 194.7 7.49827E-01 1.341E-05 279.7 12.98 9.68E-01
+82.5 1.0101E-05 6.7924E-06 0.6724 193.8 6.88237E-01 1.237E-05 279.0 12.93 1.04E+00
+83.0 9.3145E-06 6.2318E-06 0.6690 192.8 6.31441E-01 1.141E-05 278.3 12.87 1.13E+00
+83.5 8.5856E-06 5.7151E-06 0.6657 191.8 5.79087E-01 1.052E-05 277.6 12.82 1.22E+00
+84.0 7.9106E-06 5.2391E-06 0.6623 190.8 5.30848E-01 9.690E-06 276.9 12.76 1.32E+00
+84.5 7.2857E-06 4.8006E-06 0.6589 189.9 4.86416E-01 8.925E-06 276.2 12.70 1.42E+00
+85.0 6.7074E-06 4.3968E-06 0.6555 188.9 4.45509E-01 8.217E-06 275.5 12.65 1.54E+00
+85.5 6.1724E-06 4.0253E-06 0.6521 187.9 4.07863E-01 7.561E-06 274.8 12.59 1.67E+00
+86.0 5.6777E-06 3.6835E-06 0.6488 186.9 3.73232E-01 6.955E-06 274.1 12.53 1.80E+00

scripts/analytical_adc_fitting.py

@@ -0,0 +1,137 @@
+#!/usr/bin/env Python3
+import numpy as np
+import matplotlib.pyplot as plt
+from AstroAtmosphere import *
+from ADC import *
+import sys
+import scipy.optimize as opt
+
+
+
+'''
+@ Author:		Joost van den Born 
+@ Contact:		born@astron.nl
+@ Description:	This script finds the filter constant using
+				- Direct application of equation 18 and 19
+				- The minimization algorithm implemented in the ADC class
+				and compares these two.
+'''
+
+
+
+
+
+
+
+class filter:
+	def __init__(self, name, l1, l2):
+		self.name 				= name
+		self.lmin 				= np.min([l1,l2])
+		self.lmax 				= np.max([l1,l2])
+		self.bandwidth 			= self.lmax - self.lmin
+		self.optimumAngle 		= None
+		self.residualDispersion = None
+		self.optimumAngleForAllZ = None
+		self.cf 				= None
+
+	def analytic_optimum_position_ADC(self, z, l1, l2, dispersion_object, adc_object):
+		disp = dispersion_object
+		adc = adc_object
+
+		_R = disp.rc / (disp.rc + disp.H)
+		n1  = at.n_tph(l=self.lmin, T=T, p=P, RH=H, xc=xc)
+		n2  = at.n_tph(l=self.lmax, T=T, p=P, RH=H, xc=xc)
+		dn_atm = n2 - n1
+		_dispersion = np.deg2rad(disp.cassini(n1,n2,z))
+
+		_A1 = np.deg2rad(adc.p1.tilt2x_OG - adc.p1.tilt1x_OG)
+		_A2 = np.deg2rad(adc.p2.tilt2x_OG - adc.p2.tilt1x_OG)
+
+		_A3 = -1 * np.deg2rad(adc.p3.tilt2x_OG - adc.p3.tilt1x_OG) # -1 because doublet gets flipped after setting rotation
+		_A4 = -1 * np.deg2rad(adc.p4.tilt2x_OG - adc.p4.tilt1x_OG) # -1 because doublet gets flipped after setting rotation
+
+		_dn_p1 = adc.p1.indexOfRefraction(l1) - adc.p1.indexOfRefraction(l2)
+		_dn_p2 = adc.p2.indexOfRefraction(l1) - adc.p2.indexOfRefraction(l2)
+
+		_pupilRatio = 0.07943/38.542
+		_Psi = (_dn_p1 * (_A1 + _A4) + _dn_p2 * (_A2 + _A3)) * _pupilRatio
+
+		_cf = -1/_pupilRatio * (dn_atm)/((_dn_p1 * (_A1 + _A4) + _dn_p2 * (_A2 + _A3)))
+		_thetaOpt = np.arccos(1/_Psi * _dispersion)
+		return np.degrees(_thetaOpt), _cf
+
+def ADCposAnalytical(z, Cf):
+	return np.arccos(Cf*np.tan(np.deg2rad(z)))*180/np.pi
+
+
+if __name__ == "__main__":
+	Iband = filter('I-band',0.8,0.934)
+	Jband = filter('J-band',1.15,1.345)
+	Hband = filter('H-band',1.49,1.78)
+	Kband = filter('K-band',1.97,2.32)
+
+
+	##########################################################
+	# Atmospheric model
+	# 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()
+
+	# Density of the atmosphere (following CIPM-81/91 equations)
+	rho = at.rho(p=P, T=T, RH=H, xc=xc)
+
+	# Initializing dispersion object
+	disp = dispersion(lat, h)
+	disp.setReducedHeight(P, rho)
+	
+	##########################################################
+
+	adc = ADC()
+
+
+	# Unlock this when the above works
+	# zRange = range(2,62)
+	zRange = np.arange(5,65.1,0.1)
+	optimumAngleFilterList = []
+	fig = plt.figure()
+	ax = fig.add_subplot(111)
+
+
+	for f in [Hband, Iband, Jband, Kband]:
+		# For new filters or if the zenith values are different, the new values are calculated		
+		print ('Starting %s...' %f.name	)
+		# Calculating indices of refraction for l1 and l2
+		n1  = at.n_tph(l=f.lmin, T=T, p=P, RH=H, xc=xc)
+		n2  = at.n_tph(l=f.lmax, T=T, p=P, RH=H, xc=xc)
+		dispersion_obs = disp.cassini(n1, n2, zRange)*3.6e6
+
+		adc.setOrientation(0)
+		theta_opt = [adc.optimize(disp, f.lmin, f.lmax) for disp in dispersion_obs]
+		f.optimumAngleForAllZ = theta_opt
+		# Plotting the data
+		ax.plot(zRange, f.optimumAngleForAllZ, label=f.name)
+
+	for f in [Hband, Iband, Jband, Kband]:
+		popt, pcov = opt.curve_fit(ADCposAnalytical, zRange[:-50], f.optimumAngleForAllZ[:-50], p0=0.5)
+		f.cf = popt[0]
+		ax.plot(zRange, ADCposAnalytical(zRange, f.cf), 'k--', label='best fit %s' %f.name)
+		theta_opt_analytical, cf_analytical = f.analytic_optimum_position_ADC(zRange, f.lmin, f.lmax, disp, adc)
+		ax.plot(zRange, ADCposAnalytical(zRange, cf_analytical), 'b--', label='analytical cf %s' %f.name)
+		print('{} - Best fit: {} - Analytical: {}'.format(f.name, f.cf, cf_analytical))
+
+	# Plot properties
+	ax.set_title('Optimum rotation angle of the first ADC prism', fontsize=9)
+	ax.set_ylabel('Rotation angle (deg)', fontsize=9)
+	ax.set_xlabel('Zenith angle (deg)', fontsize=9)	
+	ax.grid(True, alpha=0.5)
+	ax.legend(fontsize=9)
+
+	plt.show()
+ 	

scripts/compare_models.py

@@ -0,0 +1,299 @@
+#!/usr/bin/env Python
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import host_subplot
+from AstroAtmosphere import *
+import pandas as pd
+
+'''
+@ Author:		Joost van den Born 
+@ Contact:		born@astron.nl
+@ Description:	Modeling the atmospheric refraction using multiple 
+				models of the atmosphere and then comparing them.
+
+				Models included:
+				- Full integration method
+				- Plane-parallel atmosphere
+				- Cassini's model
+				- Tan5 model (A*tan(z) + B*tan^3(z) + C*tan^5(z))
+				- Oriani's theorem (A*tan(z) + B*tan^3(z))
+				- Erf model (from Corbard et al., 2019)
+				- SLALIB
+
+				One of the figures and all three tables are used in the paper.
+'''
+
+
+if __name__ == '__main__':
+	plt.rc('font', family='serif')
+	plt.rc('axes', labelsize='large')
+
+	# Conditions at Cerro Armazones
+	# T   = 279.65 	# k
+	# P   = 71200  	# Pa
+	# H   = 0.22
+	# xc  = 450   	# ppm
+	# Conditions for the standard atmosphere
+	T   = 288.15 	# k
+	P   = 101325  	# Pa
+	H   = 0.0
+	xc  = 314   	# ppm
+	lat = -24.5983 	# degrees
+	h   = 3064  	# m
+	l1  = 1.49 		# micron
+	l2  = 1.78 		# micron
+	zenith = 30 	# degrees
+
+	# Calculate the index of refraction
+	at  = Observatory()
+	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)
+	rho = at.rho(p=P, T=T, RH=H, xc=xc)
+
+	refractionAnalysis  = True
+	dispersionAnalysis  = True
+	tableCompare 		= True
+	errorTable 			= True
+	if dispersionAnalysis:
+		# Initiate the refraction class
+		ref = dispersion(lat, h)
+		ref.setReducedHeight(P, rho)
+
+		# Compare the refraction models to check if the output is okay
+		print ('Plane - parallel atmosphere dispersion: %s mas' %(ref.planeParallel(n1, n2, zenith) * 3.6e6))
+		print ('Cassini dispersion: %s mas' %(ref.cassini(n1, n2, zenith) * 3.6e6))
+		print ('Tan5 dispersion: %s mas' %(ref.tan5(n1, n2, zenith) * 3.6e6))
+		print ('Oriani dispersion: %s mas' %(ref.oriani(n1, n2, zenith) * 3.6e6))
+		print ('Corbard dispersion: %s mas' %(ref.corbard(n1, n2, zenith)*3.6e6))
+		print ('Full integration: %s mas' %(ref.refractionIntegral(n1, n2, zenith)*3.6e6))
+		print('-'*35)
+
+
+		# Now we calculate the dispersion over a wide range of zenith angles
+		# for ideal gas atmosphere
+		zenithRange = np.linspace(5, 85, 1000)
+		ppList 		= ref.planeParallel(n1, n2, zenithRange) * 3.6e6
+		cassiniList = ref.cassini(n1, n2, zenithRange) * 3.6e6
+		tan5List 	= ref.tan5(n1, n2, zenithRange) * 3.6e6
+		orianiList  = ref.oriani(n1, n2, zenithRange) * 3.6e6
+		corbardList = [ref.corbard(n1, n2, z) * 3.6e6 for z in zenithRange]
+		fullIntegrationList = [ref.refractionIntegral(n1, n2, z) * 3.6e6 for z in zenithRange]
+
+		# for isothermal atmosphere
+		cassiniList_isotherm = ref.cassini(n1, n2, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e6
+		tan5List_isotherm 	 = ref.tan5(n1, n2, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e6
+		orianiList_isotherm  = ref.oriani(n1, n2, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e6
+		corbardList_isotherm = [ref.corbard(n1, n2, z, reduced_height=ref.H_isotherm(T)) * 3.6e6 for z in zenithRange]
+
+		# check when the variation between the different models pass different values
+		dispersionWidth = []
+		widthList = []
+		for z in range(len(zenithRange)):
+			compList = [cassiniList[z], tan5List[z], orianiList[z], corbardList[z], cassiniList_isotherm[z], tan5List_isotherm[z], orianiList_isotherm[z], corbardList_isotherm[z], fullIntegrationList[z]]
+			width    = np.max(compList) - np.min(compList)
+			widthList.append(width)
+			dispersionWidth.append((zenithRange[z], width))
+
+		masBorder = np.argmax(np.asarray(widthList)>1)
+		p1masBorder = np.argmax(np.asarray(widthList)>0.1)
+		p01masBorder = np.argmax(np.asarray(widthList)>0.01)
+		p001masBorder = np.argmax(np.asarray(widthList)>0.001)
+
+		plotAbsolute = 1
+		if plotAbsolute:
+			# Plotting it all, with vertical lines where the variation between the models exceed 1, 0.1, 0.01 and 0.001 mas.
+			fig = plt.figure(figsize=(9,5))
+			ax  = fig.add_subplot(111)
+			ax.set_xlim([30,85])
+			ax.set_xlabel('Zenith (degrees)')
+			ax.set_ylabel('Dispersion (mas)')
+			ax.grid(True, alpha=0.3)
+			ax.axvline(dispersionWidth[masBorder][0], linestyle='--', color='k', linewidth=0.5)
+			ax.text(dispersionWidth[masBorder][0]-0.5, 402, '1 mas', ha='center', va='bottom', rotation='vertical')
+			ax.axvline(dispersionWidth[p1masBorder][0], linestyle='--', color='k', linewidth=0.5)
+			ax.text(dispersionWidth[p1masBorder][0]-0.5, 402, '0.1 mas', ha='center', va='bottom', rotation='vertical')
+			ax.axvline(dispersionWidth[p01masBorder][0], linestyle='--', color='k', linewidth=0.5)
+			ax.text(dispersionWidth[p01masBorder][0]-0.5, 402, '10 $\mu$as', ha='center', va='bottom', rotation='vertical')
+			ax.axvline(dispersionWidth[p001masBorder][0], linestyle='--', color='k', linewidth=0.5)
+			ax.text(dispersionWidth[p001masBorder][0]-0.5, 402, '1 $\mu$as', ha='center', va='bottom', rotation='vertical')
+
+			ax.plot(zenithRange, ppList, 'k-', label='Plane-parallel', alpha=0.2)
+			ax.plot(zenithRange, cassiniList, 'g-', label='Cassini')
+			ax.plot(zenithRange, tan5List, 'b-', label='$Atan(z) - Btan^3(z) + Ctan^5(z)$')
+			ax.plot(zenithRange, orianiList, 'y-', label='$Atan(z) - Btan^3(z)$')
+			ax.plot(zenithRange, corbardList, 'r-', label='Error function')
+			ax.plot(zenithRange, fullIntegrationList, 'c-', label='Full integration')
+
+			ax.plot(zenithRange, cassiniList_isotherm, 'g--', label='Cassini (isotherm)')
+			ax.plot(zenithRange, tan5List_isotherm, 'b--', label='$Atan(z) - Btan^3(z) + Ctan^5(z)$ (isotherm)')
+			ax.plot(zenithRange, orianiList_isotherm, 'y--', label='$Atan(z) - Btan^3(z)$ (isotherm)')
+			ax.plot(zenithRange, corbardList_isotherm, 'r--', label='Error function (isotherm)')
+
+
+			ax.tick_params(axis ='both', direction='in') 
+			ax.legend()
+			plt.tight_layout()
+			# plt.show()
+
+		plotRelative = 1
+		if plotRelative:
+			# RELATIVE TO FULL INTEGRATION
+			fig = plt.figure(figsize=(9,6))
+			ax  = fig.add_subplot(111)
+			ax.set_xlim([30,85])
+			
+			ax.set_xlabel('Zenith (degrees)', fontsize=18)
+			ax.set_ylabel(r'$\log_{10}(|(\Delta R - \Delta R_{int})|$ mas$^{-1})$', fontsize=18)
+			
+			ax.axhline(np.log10(0.4), linestyle='-', linewidth=1, color='red')
+			ax.text(30.5, np.log10(0.4)+0.05, 'Residual dispersion requirement', fontsize=12)
+
+			ax.grid(True, which='major', alpha=0.3)
+			ax.grid(True, which='minor', alpha=0.2)
+
+			ax.plot(zenithRange, np.log10(np.abs(ppList-fullIntegrationList)), 					color='k', linestyle='-',  linewidth=2, label='Plane-parallel')
+			ax.plot(zenithRange, np.log10(np.abs(cassiniList-fullIntegrationList)), 			color='g', linestyle=(0, (1, 1)),  linewidth=2, label='Cassini')
+			ax.plot(zenithRange, np.log10(np.abs(orianiList-fullIntegrationList)), 				color='y', linestyle='-.', linewidth=2, label='Oriani')
+			ax.plot(zenithRange, np.log10(np.abs(tan5List-fullIntegrationList)), 				color='b', linestyle=(0, (5, 1)), linewidth=2, label='Oriani incl. fifth order')
+			ax.plot(zenithRange, np.log10(np.abs(np.asarray(corbardList)-fullIntegrationList)), color='r', linestyle=(0, (3, 1, 1, 1, 1, 1)), linewidth=2, label='Error function')
+			ax.plot(zenithRange, np.log10(np.abs(cassiniList-fullIntegrationList)), 			color='g', linestyle=(0, (1, 1)),  linewidth=2)
+
+			ax.set_ylim(top=2.6, bottom=-7.2)
+			ax.tick_params(axis ='both', direction='in', labelsize=16) 
+			ax.legend(fontsize=14, labelspacing=0.3, columnspacing=0.5, handletextpad=0.3,  ncol=1, loc='lower right')
+			plt.tight_layout()
+			# plt.savefig('./figures/dispersionModels_with_SLALIB_minimal2.pdf', transparent=True, bbox_inches='tight')
+			# plt.savefig('./figures/dispersionModels_with_SLALIB_minimal2.png', transparent=True, bbox_inches='tight')
+
+
+		
+
+	if refractionAnalysis:
+		# Initiate the refraction class
+		ref = refraction(lat, h)
+		ref.setReducedHeight(P, rho)
+
+
+		# Compare the refraction models to check if the output is okay
+		print ('Plane - parallel atmosphere refraction: %s arcsec' %(ref.planeParallel(n1, zenith) * 3.6e3))
+		print ('Cassini refraction: %s arcsec' %(ref.cassini(n1, zenith) * 3.6e3))
+		print ('Tan5 refraction: %s arcsec' %(ref.tan5(n1, zenith) * 3.6e3))
+		print ('Oriani refraction: %s arcsec' %(ref.oriani(n1, zenith) * 3.6e3))
+		print ('Corbard refraction: %s arcsec' %(ref.corbard(n1, zenith)*3.6e3))
+		print ('Full integration: %s mas' %(ref.refractionIntegral(n1, zenith)*3.6e3))
+		print('-'*35)
+
+
+		# Now we calculate the dispersion over a wide range of zenith angles
+		# for ideal gas atmosphere
+		zenithRange = np.linspace(0, 88, 1000)
+		ppList 		= ref.planeParallel(n1, zenithRange) * 3.6e3
+		cassiniList = ref.cassini(n1, zenithRange) * 3.6e3
+		tan5List 	= ref.tan5(n1, zenithRange) * 3.6e3
+		orianiList  = ref.oriani(n1, zenithRange) * 3.6e3
+		corbardList = [ref.corbard(n1, z) * 3.6e3 for z in zenithRange]
+		fullIntegrationList = [ref.refractionIntegral(n1, z) * 3.6e3 for z in zenithRange]
+
+		# for isothermal atmosphere
+		cassiniList_isotherm = ref.cassini(n1, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e3
+		tan5List_isotherm 	 = ref.tan5(n1, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e3
+		orianiList_isotherm  = ref.oriani(n1, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e3
+		corbardList_isotherm = [ref.corbard(n1, z, reduced_height=ref.H_isotherm(T)) * 3.6e3 for z in zenithRange]
+
+		# Plotting it all, with vertical lines where the variation between the models exceed 1, 0.1, 0.01 and 0.001 arcsec.
+		fig = plt.figure(figsize=(9,5))
+		ax  = fig.add_subplot(111)
+		ax.set_xlim([80,88])
+		ax.set_xlabel('Zenith (degrees)')
+		ax.set_ylabel('Refraction (arcsec)')
+		ax.grid(True, alpha=0.3)
+
+		ax.plot(zenithRange, ppList, 'k-', label='Plane-parallel', alpha=0.2)
+		ax.plot(zenithRange, cassiniList, 'g-', label='Cassini')
+		ax.plot(zenithRange, tan5List, 'b-', label='$Atan(z) - Btan^3(z) + Ctan^5(z)$')
+		ax.plot(zenithRange, orianiList, 'y-', label='$Atan(z) - Btan^3(z)$')
+		ax.plot(zenithRange, corbardList, 'r-', label='Error function')
+		ax.plot(zenithRange, fullIntegrationList, 'c-', label='Full integration')
+
+
+		ax.tick_params(axis ='both', direction='in') 
+		ax.legend()
+		plt.tight_layout()
+		
+		
+
+	if tableCompare:
+		refTableData = []
+		dispTableData = []
+		ref = refraction(lat, h)
+		ref.setReducedHeight(P, rho)		
+		disp = dispersion(lat, h)
+		disp.setReducedHeight(P, rho)
+		zenithRange = [10, 30, 50, 60, 70, 80, 85]
+		headerList  = np.array(['z', 'Plane-parallel', 'Cassini', 'tan3', 'tan5', 'Error function', 'Full integration'])
+		refTableData.append(headerList)
+		dispTableData.append(headerList)
+		
+		for z in zenithRange:
+			ppValue 		= round(ref.planeParallel(n1, z) * 3.6e3, 3)
+			cassiniValue 	= round(ref.cassini(n1, z) * 3.6e3, 3)
+			tan3Value 		= round(ref.oriani(n1, z)  * 3.6e3, 3)
+			tan5Value 		= round(ref.tan5(n1, z)  * 3.6e3, 3)
+			erfValue 		= round(ref.corbard(n1, z) * 3.6e3, 3)
+			fullIntValue 	= round(ref.refractionIntegral(n1, z) * 3.6e3, 3)
+			refTableData.append([z, ppValue, cassiniValue, tan3Value, tan5Value, erfValue, fullIntValue])
+		
+			ppValue 		= round(disp.planeParallel(n1, n2, z) * 3.6e6, 3)
+			cassiniValue 	= round(disp.cassini(n1, n2, z) * 3.6e6, 3)
+			tan3Value 		= round(disp.oriani(n1, n2, z)  * 3.6e6, 3)
+			tan5Value 		= round(disp.tan5(n1, n2, z)  * 3.6e6, 3)
+			erfValue 		= round(disp.corbard(n1, n2, z) * 3.6e6, 3)
+			fullIntValue 	= round(disp.refractionIntegral(n1, n2, z) * 3.6e6, 3)
+			dispTableData.append([z, ppValue, cassiniValue, tan3Value, tan5Value, erfValue, fullIntValue])
+
+		dfr = pd.DataFrame(np.asarray(refTableData))
+		dfd = pd.DataFrame(np.asarray(dispTableData))
+		print('-'*35,'\nREFRACTION COMPARISON:\n')
+		print(dfr.to_string(index=False))
+		print('-'*35)
+		print('-'*35,'\nDISPERSION COMPARISON:\n')
+		print(dfd.to_string(index=False))
+		print('-'*35)
+			
+		# print(dfr.to_latex(index=False))
+		# print(dfd.to_latex(index=False))
+
+
+	if errorTable:
+		ref = dispersion(lat, h)
+		ref.setReducedHeight(P, rho)
+		dispTableData = []
+		headerList  = np.array(['s_P', 's_T', 's_H', 's_disp', 's_disp(dl=0)'])
+		dispTableData.append(headerList)
+
+		# dPlist = [2, 20, 200, 0,    0,   0, 0,     0,    0,   0,  0, 0,   2,     20,   200]
+		# dTlist = [0, 0,  0,   0.02, 0.2, 2, 0,     0,    0,   0,  0, 0,   0.02,  0.2,  2  ]
+		# dRHlist= [0, 0,  0,   0,    0,   0, 0.002, 0.02, 0.2, 0,  0, 0,   0.002, 0.02, 0.2]
+		# dCO2list=[0, 0,  0,   0,    0,   0, 0,     0,    0,   2, 20, 200, 2,     20,   200]
+		dPlist = [0,    10,    10,   20,  100]
+		dTlist = [0,   0.1,   0.2,  0.2,  1.0]
+		dRHlist= [0,  0.01,  0.01, 0.02, 0.05]
+		dCO2list=[0,     20,     20,    20,    20]
+
+		for i in range(len(dPlist)):
+			# nerr = at.dn_tph(l1, T, P, H, xc, dl=0, dT=dTlist[i], dP=dPlist[i], dRH=dRHlist[i], dCO2=dCO2list[i])
+			disp_error = ref.cassiniError(zenith, l1, l2, T, P, H, xc, dl1=0.001, dl2=0.001, dT=dTlist[i], dP=dPlist[i], dRH=dRHlist[i], 
+				dCO2=dCO2list[i], dz=0) * 3.6e9
+
+			disp_error_dl0 = ref.cassiniError(zenith, l1, l2, T, P, H, xc, dl1=0, dl2=0, dT=dTlist[i], dP=dPlist[i], dRH=dRHlist[i], 
+				dCO2=dCO2list[i], dz=0) * 3.6e9
+
+			dispTableData.append([dPlist[i], dTlist[i], dRHlist[i], disp_error, disp_error_dl0])
+			
+		dfd = pd.DataFrame(np.asarray(dispTableData))
+		
+		print('-'*35,'\nUNCERTAINTY ANALYSIS:\n')
+		print(dfd.to_string(index=False))
+		print('-'*35)
+
+	plt.show()

scripts/compare_models_with_SLALIB.py

@@ -0,0 +1,247 @@
+#!/usr/bin/env Python
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import host_subplot
+from AstroAtmosphere import *
+from pyslalib import slalib
+
+
+'''
+@ Author:		Joost van den Born 
+@ Contact:		born@astron.nl
+@ Description:	This script was used to compare the impact of using the old 
+				B&S equations, and the particular implementation of SLALIB. 
+				The rest of the script is almost a one-to-one copy of other 
+				script that compared the atmospheric models.
+
+				Modeling the atmospheric refraction using multiple models of 
+				the atmosphere and then comparing them.
+
+				Models included:
+				- Full integration method
+				- Plane-parallel atmosphere
+				- Cassini's model
+				- Tan5 model (A*tan(z) + B*tan^3(z) + C*tan^5(z))
+				- Oriani's theorem (A*tan(z) + B*tan^3(z))
+				- Erf model (from Corbard et al., 2019)
+				- SLALIB
+'''
+class atmosphereSLALIB:
+	def __init__(self, T=279.65, P=712, H=0.22):
+		self.name 	= 'SLALIB atmosphere'
+		self.T 		= T
+		self.P 		= P
+		self.H 		= H
+		self.h 		= 3064
+		self.lat 	= -24.5983
+		self.TLR 	= 0.0065
+		self.prec 	= 1e-8
+
+	def atmModel(self, A, B, zenith):
+		zenith = np.radians(zenith)
+		return A * np.tan(zenith) + B * (np.tan(zenith)**3)
+
+	def dispersion(self, l1, l2, z):
+		A1, B1 = slalib.sla_refco(self.h, self.T, self.P, self.H, l1, np.deg2rad(self.lat), self.TLR, self.prec)
+		A2, B2 = slalib.sla_atmdsp(self.T, self.P, self.H, l1, A1, B1, l2)
+		D1 = np.degrees(self.atmModel(A1, B1, z)) * 3600
+		D2 = np.degrees(self.atmModel(A2, B2, z)) * 3600
+		dispersionValue = ((D1 - D2) * 1000)
+		return dispersionValue
+
+	def refractivity(self, l=None):
+		'''
+			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.
+		'''
+		zObs = np.deg2rad(45.06563)
+		h = 0
+		T = self.T
+		P = self.P
+		H = self.H
+		lat = 24
+		tlr = 0.0065
+		eps = 1e-10
+		if not l:
+			l = self.lrange
+		R = []
+		if type(l) == float:
+			l = [l]
+		for wl in l:
+			R.append(slalib.sla_refro(zObs, h, T, P, H, wl, lat, tlr, eps)+1)
+		return np.asarray(R)
+
+
+if __name__ == '__main__':
+	plt.rc('font', family='serif')
+	# plt.rc('xtick', labelsize='x-small')
+	# plt.rc('ytick', labelsize='x-small')
+	plt.rc('axes', labelsize='large')
+
+	# Standard atmospheric conditions
+	# T   = 279.65 	# k
+	# P   = 71200  	# Pa
+	# H   = 0.22
+	# xc  = 450   	# ppm
+	T   = 288.15 	# k
+	P   = 101325  	# Pa
+	H   = 0.0
+	xc  = 314   	# ppm
+	lat = -24.5983 	# degrees
+	h   = 3064  	# m
+	l1  = 1.49 		# micron
+	l2  = 1.78 		# micron
+	zenith = 30 	# degrees
+
+	# Calculate the index of refraction
+	at  = Observatory()
+	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)
+	rho = at.rho(p=P, T=T, RH=H, xc=xc)
+	
+	#######################################################
+	# DIRECT COMPARISON
+	ref = dispersion(lat, h)
+	ref.setReducedHeight(P, rho)
+	print('Ciddor refractivity: ', n1, n2)
+	print ('Cassini dispersion: %s mas' %(ref.cassini(n1, n2, zenith) * 3.6e6))
+	n1c = n1
+	n2c = n2
+
+	refraction = refraction(lat, h)
+	refraction.setReducedHeight(P, rho)
+	refractionValueCiddor = refraction.cassini(n1c, zenith) * 3.6e6 # in mas
+
+	# SLALIB
+	atSLALIB = atmosphereSLALIB(T, P/100, H)
+	n1 = atSLALIB.refractivity(l1)
+	n2 = atSLALIB.refractivity(l2)
+
+	print('B&S refractivity: ', n1, n2)
+	print ('SLALIB model: %s mas' %(atSLALIB.dispersion(l1, l2, zenith) ) )
+
+	# SLALIB refraction:
+	A1, B1 = slalib.sla_refco(atSLALIB.h, atSLALIB.T, atSLALIB.P, atSLALIB.H, l1, np.deg2rad(atSLALIB.lat), atSLALIB.TLR, atSLALIB.prec)
+	refractionValueSLALIB = np.degrees(atSLALIB.atmModel(A1, B1, zenith)) * 3.6e6 # in mas
+
+	print('Differential dispersion = {} mas'.format(atSLALIB.dispersion(l1, l2, zenith) - ref.cassini(n1c, n2c, zenith) * 3.6e6) )
+	print('Differential refraction = {} mas'.format( refractionValueSLALIB - refractionValueCiddor ) )
+
+	#######################################################
+
+	refractionAnalysis  = True
+	dispersionAnalysis  = True
+	tableCompare 		= True
+	if dispersionAnalysis:
+		# Initiate the refraction class
+		ref = dispersion(lat, h)
+		ref.setReducedHeight(P, rho)
+
+		# Compare the refraction models to check if the output is okay
+		print ('Plane - parallel atmosphere dispersion: %s mas' %(ref.planeParallel(n1, n2, zenith) * 3.6e6))
+		print ('Cassini dispersion: %s mas' %(ref.cassini(n1, n2, zenith) * 3.6e6))
+		print ('Tan5 dispersion: %s mas' %(ref.tan5(n1, n2, zenith) * 3.6e6))
+		print ('Oriani dispersion: %s mas' %(ref.oriani(n1, n2, zenith) * 3.6e6))
+		print ('Corbard dispersion: %s mas' %(ref.corbard(n1, n2, zenith)*3.6e6))
+		print ('Full integration: %s mas' %(ref.refractionIntegral(n1, n2, zenith)*3.6e6))
+		print ('SLALIB model: %s mas' %(atSLALIB.dispersion(l1, l2, zenith) ) )
+
+		# Now we calculate the dispersion over a wide range of zenith angles
+		# for ideal gas atmosphere
+		zenithRange = np.linspace(30, 85, 1000)
+		ppList 		= ref.planeParallel(n1, n2, zenithRange) * 3.6e6
+		cassiniList = ref.cassini(n1, n2, zenithRange) * 3.6e6
+		tan5List 	= ref.tan5(n1, n2, zenithRange) * 3.6e6
+		orianiList  = ref.oriani(n1, n2, zenithRange) * 3.6e6
+		corbardList = [ref.corbard(n1, n2, z)[0] * 3.6e6 for z in zenithRange]
+		fullIntegrationList = np.asarray([ref.refractionIntegral(n1, n2, z) * 3.6e6 for z in zenithRange])
+		slalibList	= atSLALIB.dispersion(l1, l2, zenithRange)
+
+		# for isothermal atmosphere
+		cassiniList_isotherm = ref.cassini(n1, n2, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e6
+		tan5List_isotherm 	 = ref.tan5(n1, n2, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e6
+		orianiList_isotherm  = ref.oriani(n1, n2, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e6
+		corbardList_isotherm = [ref.corbard(n1, n2, z, reduced_height=ref.H_isotherm(T))[0] * 3.6e6 for z in zenithRange]
+
+		# check when the variation between the different models pass different values
+		dispersionWidth = []
+		widthList = []
+		for z in range(len(zenithRange)):
+			compList = [cassiniList[z], tan5List[z], orianiList[z], corbardList[z], cassiniList_isotherm[z], tan5List_isotherm[z], orianiList_isotherm[z], corbardList_isotherm[z], slalibList[z]]#, fullIntegrationList[z]]
+			width    = np.max(compList) - np.min(compList)
+			widthList.append(width)
+			dispersionWidth.append((zenithRange[z], width))
+
+		masBorder = np.argmax(np.asarray(widthList)>1)
+		p1masBorder = np.argmax(np.asarray(widthList)>0.1)
+		p01masBorder = np.argmax(np.asarray(widthList)>0.01)
+		p001masBorder = np.argmax(np.asarray(widthList)>0.001)
+
+		plotAbsolute = 1
+		if plotAbsolute:
+			# Plotting it all, with vertical lines where the variation between the models exceed 1, 0.1, 0.01 and 0.001 mas.
+			fig = plt.figure(figsize=(9,5))
+			ax  = fig.add_subplot(111)
+			ax.set_xlim([30,85])
+			ax.set_xlabel('Zenith (degrees)')
+			ax.set_ylabel('Dispersion (mas)')
+			ax.grid(True, alpha=0.3)
+			ax.axvline(dispersionWidth[masBorder][0], linestyle='--', color='k', linewidth=0.5)
+			ax.text(dispersionWidth[masBorder][0]-0.5, 402, '1 mas', ha='center', va='bottom', rotation='vertical')
+			ax.axvline(dispersionWidth[p1masBorder][0], linestyle='--', color='k', linewidth=0.5)
+			ax.text(dispersionWidth[p1masBorder][0]-0.5, 402, '0.1 mas', ha='center', va='bottom', rotation='vertical')
+			ax.axvline(dispersionWidth[p01masBorder][0], linestyle='--', color='k', linewidth=0.5)
+			ax.text(dispersionWidth[p01masBorder][0]-0.5, 402, '10 $\mu$as', ha='center', va='bottom', rotation='vertical')
+			ax.axvline(dispersionWidth[p001masBorder][0], linestyle='--', color='k', linewidth=0.5)
+			ax.text(dispersionWidth[p001masBorder][0]-0.5, 402, '1 $\mu$as', ha='center', va='bottom', rotation='vertical')
+
+			ax.plot(zenithRange, ppList, 'k-', label='Plane-parallel', alpha=0.2)
+			ax.plot(zenithRange, cassiniList, 'g-', label='Cassini')
+			ax.plot(zenithRange, tan5List, 'b-', label='$Atan(z) - Btan^3(z) + Ctan^5(z)$')
+			ax.plot(zenithRange, orianiList, 'y-', label='$Atan(z) - Btan^3(z)$')
+			ax.plot(zenithRange, corbardList, 'r-', label='Error function')
+			ax.plot(zenithRange, fullIntegrationList, 'c-', label='Full integration')
+
+			ax.plot(zenithRange, cassiniList_isotherm, 'g--', label='Cassini (isotherm)')
+			ax.plot(zenithRange, tan5List_isotherm, 'b--', label='$Atan(z) - Btan^3(z) + Ctan^5(z)$ (isotherm)')
+			ax.plot(zenithRange, orianiList_isotherm, 'y--', label='$Atan(z) - Btan^3(z)$ (isotherm)')
+			ax.plot(zenithRange, corbardList_isotherm, 'r--', label='Error function (isotherm)')
+
+
+			ax.tick_params(axis ='both', direction='in') 
+			ax.legend()
+			plt.tight_layout()
+			plt.show()
+
+		plotRelative = 1
+		if plotRelative:
+			# RELATIVE TO FULL INTEGRATION
+			fig = plt.figure(figsize=(9,6))
+			ax  = fig.add_subplot(111)
+			ax.set_xlim([30,85])
+			
+			ax.set_xlabel('Zenith (degrees)', fontsize=18)
+			ax.set_ylabel(r'$\log_{10}(|(\Delta R - \Delta R_{int})|$ mas$^{-1})$', fontsize=18)
+			
+			ax.axhline(np.log10(0.4), linestyle='-', linewidth=1, color='red')
+			ax.text(30.5, np.log10(0.4)+0.05, 'Residual dispersion requirement', fontsize=12)
+
+			ax.grid(True, which='major', alpha=0.3)
+			ax.grid(True, which='minor', alpha=0.2)
+
+			ax.plot(zenithRange, np.log10(np.abs(ppList-fullIntegrationList)), 					color='k', linestyle='-',  linewidth=2, label='Plane-parallel')
+			ax.plot(zenithRange, np.log10(np.abs(cassiniList-fullIntegrationList)), 			color='g', linestyle=(0, (1, 1)),  linewidth=2, label='Cassini')
+			ax.plot(zenithRange, np.log10(np.abs(orianiList-fullIntegrationList)), 				color='y', linestyle='-.', linewidth=2, label='Oriani')
+			ax.plot(zenithRange, np.log10(np.abs(tan5List-fullIntegrationList)), 				color='b', linestyle=(0, (5, 1)), linewidth=2, label='Oriani incl. fifth order')
+			ax.plot(zenithRange, np.log10(np.abs(np.asarray(corbardList)-fullIntegrationList)), color='r', linestyle=(0, (3, 1, 1, 1, 1, 1)), linewidth=2, label='Error function')
+			ax.plot(zenithRange, np.log10(np.abs(np.asarray(slalibList)-fullIntegrationList)), 	color='m', linewidth=2, label='SLALIB')
+			ax.plot(zenithRange, np.log10(np.abs(cassiniList-fullIntegrationList)), 			color='g', linestyle=(0, (1, 1)),  linewidth=2)
+
+			ax.set_ylim(top=2.6, bottom=-7.2)
+			ax.tick_params(axis ='both', direction='in', labelsize=16) 
+			ax.legend(fontsize=14, labelspacing=0.3, columnspacing=0.5, handletextpad=0.3,  ncol=1, loc='lower right')
+			plt.tight_layout()
+			plt.show()
+
+

scripts/compare_refractivity_models.py

@@ -0,0 +1,177 @@
+#!/usr/bin/env Python3
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import host_subplot
+from matplotlib.colors import LogNorm
+from tabulate import tabulate
+import AstroAtmosphere
+
+
+'''
+@ Author:		Joost van den Born 
+@ Contact:		born@astron.nl
+@ 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).
+'''
+
+
+class atmosphericModel:
+
+	def __init__(self, T, P, H, CO2, l1, l2):
+		self.T 		= T
+		self.P 		= P
+		self.H 		= H
+		self.CO2	= CO2
+		self.lrange = np.linspace(l1, l2, 100)
+
+
+	def BarrellAndSears(self, l=None):
+		if not l:
+			l = self.lrange
+		return (AstroAtmosphere.BarrellAndSears(l, T=self.T, p=self.P)-1)
+
+
+	def Ciddor(self, l=None):
+		if not l:
+			l = self.lrange
+		return (AstroAtmosphere.Ciddor(l, CO2=self.CO2)-1)
+
+
+	def Edlen1953(self, l=None):
+		if not l:
+			l = self.lrange
+		return(AstroAtmosphere.Edlen1953(l)-1)
+
+
+	def Edlen1966(self, l=None):
+		if not l:
+			l = self.lrange
+		return(AstroAtmosphere.Edlen1966(l)-1)	
+
+
+	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)
+
+
+	def BonschPotulski(self, l=None):
+		if not l:
+			l = self.lrange
+		return (AstroAtmosphere.BonschPotulski(l, T=self.T, p=self.P, CO2=self.CO2)-1)
+
+	def BirchDowns(self, l=None):
+		if not l:
+			l = self.lrange
+		return (AstroAtmosphere.BirchDowns(l)-1)
+
+	def PeckReeder(self, l=None):
+		if not l:
+			l = self.lrange
+		return (AstroAtmosphere.PeckReeder(l)-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)
+
+
+if __name__ == '__main__':
+
+	# Calculation of the refractivity of the different models
+	# multiplication by 1e8 to make the differentials more
+	# pronounced.
+	at = atmosphericModel(T=288.15, P=101325, H=0, CO2=300, l1=0.3, l2=2.5)
+	BSdata 				= 1e8 * at.BarrellAndSears()
+	CiddorData 			= 1e8 * at.Ciddor()
+	Edlen1953Data 		= 1e8 * at.Edlen1953()
+	Edlen1966Data 		= 1e8 * at.Edlen1966()
+	OwensData 			= 1e8 * at.Owens()
+	BonschPotulskiData 	= 1e8 * at.BonschPotulski()
+	BirchDownsData 		= 1e8 * at.BirchDowns()
+	PeckReederData 		= 1e8 * at.PeckReeder()
+	slalibData 			= 1e8 * at.slalib()
+	
+	# Plotting the results
+	fig = plt.figure(figsize=(9,5))
+	ax2 = fig.add_subplot(111)
+	ax2.grid(True, alpha=0.3)
+	ax2.plot(at.lrange, BSdata - CiddorData, 'b-', label='Barrell & Sears (1939)')
+	ax2.plot(at.lrange, CiddorData - CiddorData, 'k-', label='Ciddor (1996)')
+	ax2.plot(at.lrange, Edlen1953Data - CiddorData, 'r:', label='Edlen (1953)')
+	ax2.plot(at.lrange, Edlen1966Data - CiddorData, 'g--', label='Edlen (1966)')
+	ax2.plot(at.lrange, OwensData - CiddorData, 'm--', label='Owens (1967)')
+	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, 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])
+	plt.tight_layout()
+	plt.show()
+
+
+	
+	# Let's compare the impact on the atmospheric dispersion.
+	# For the dispersion we must compare the delta for two lambda values
+	methodList = 	[	at.BarrellAndSears,
+						at.Ciddor,
+						at.Edlen1953,
+						at.Edlen1966,
+						at.Owens,
+						at.BonschPotulski,
+						at.BirchDowns,
+						at.PeckReeder,
+						at.slalib
+					]
+	
+
+	filters = 	[	[0.800, 0.934], # I-band
+					[1.150, 1.345], # J-band
+					[1.490, 1.780], # H-band
+					[1.970, 2.320]	# K-band
+				]
+	filterDispersion = []
+	
+	# Iterate over the four MICADO broad band filters
+	for f in filters:
+		# Wavelengths of this particular filter
+		l1 = f[0]
+		l2 = f[1]
+
+		dispersionList = []
+		for m in methodList:
+			d1 = m(l=l1)
+			d2 = m(l=l2)
+			deltan = d1 - d2
+			
+			# We assume, for sake of simplicity, a plane parallel atmosphere
+			# and calculate the dispersion in milli arcseconds.
+			dispersion = deltan * np.tan(np.deg2rad(45)) * 206265 * 1000 
+			dispersionList.append(dispersion)
+
+		# We compare the atmospheric dispersion relative to the atmospheric
+		# dispersion as calculated using the Ciddor model.
+		dispersionList = np.asarray(dispersionList) - dispersionList[1]
+		filterDispersion.append(dispersionList)
+
+	# Print the results in a table.
+	methodStringList = ['Barrell & Sears (1939)',
+						'Ciddor (1996)',
+						'Edlen (1953)',
+						'Edlen (1966)',
+						'Bonsch & Potulski (1998)',
+						'Birch & Downs (1993, 1994)',
+						'Peck & Reeder (1972)',
+						'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)']))
+

scripts/differential_dispersion.py

@@ -0,0 +1,70 @@
+#!/usr/bin/env Python3
+import numpy as np
+import matplotlib.pyplot as plt
+import AstroAtmosphere as atm
+from ADC import * 
+
+'''
+	@author: 		Joost van den Born
+	@Contact:		born@astron.nl
+	@Description: 	Caclulates the differential dispersion between the ADC and 
+					the atmosphere following the ciddor and cassini dispersion
+					model.
+					This script generates figure 4 from the paper.
+'''
+
+# Initializing ADC object and defining the filter bands
+adc = ADC()
+bandName = ['I-band', 'J-band', 'H-band', 'K-band']
+bandList = [(0.8, 0.934), (1.15, 1.345), (1.49, 1.78), (1.97, 2.32)]
+
+# Some empty lists for later use
+wavelengthList = []
+atDispersionList = []
+adcDispersionList = []
+
+for i in range(len(bandList)):
+	# We want to find the dispersion of the atmosphere and 
+	# and of the ADC, at z=30 deg.
+	lMin 			= bandList[i][0]
+	lMax 			= bandList[i][1]
+	z 				= 30 # degrees
+	lRange 			= np.linspace(lMin, lMax, 100)
+	wavelengthList.append(lRange)
+
+	# Set ADC position optimized for H-band, z=30 deg.
+	atDispersion 	= atm.quick_dispersion(lMin, lMax, z, conditions='CERRO_ARMAZONES') * 3.6e6
+	adcOptimumAngle = adc.optimize(atDispersion, lMin, lMax) 
+	adc.setOrientation(adcOptimumAngle)
+
+	# Now iterate over the band to see how the differential dispersion changes
+	atDispersionValues  = []
+	adcDispersionValues = []
+	for l in lRange:
+		atDispersionValues.append(atm.quick_dispersion(np.min(lRange), l, z, conditions='CERRO_ARMAZONES')* 3.6e6 )
+		adcDispersionValues.append(adc.calculateDispersion(np.min(lRange),l)[1])
+	atDispersionList.append(atDispersionValues)
+	adcDispersionList.append(adcDispersionValues)
+
+# Calculating the differential dispersion
+atDispersionList 	= np.asarray(atDispersionList)
+adcDispersionList 	= np.asarray(adcDispersionList)	
+deltaDispersion 	= atDispersionList - adcDispersionList
+
+
+
+# Plot the results
+plt.rc('font', family='serif')
+plt.rc('axes', labelsize='large')
+fig = plt.figure(figsize=(9,6))
+
+ax = fig.add_subplot(111)
+linestyleList = [':', '-.', '--', '-']
+for i in range(len(bandList)):
+	ax.plot(wavelengthList[i], deltaDispersion[i], label=bandName[i], linestyle=linestyleList[i], linewidth=2, color='k')
+ax.set_ylabel('Differential $\Delta$R (mas)', fontsize=18)
+ax.set_xlabel('$\lambda$ ($\mu$m)', fontsize=18)
+ax.grid(True, alpha=0.3)
+ax.legend(fontsize=16, loc='lower right')
+ax.tick_params(axis ='both', direction='in', labelsize=16) 
+plt.show()

scripts/spectral_differential_refraction.py

@@ -0,0 +1,324 @@
+#!/usr/bin/env Python
+from __future__ import division
+import numpy as np
+import matplotlib.pyplot as plt
+import AstroAtmosphere as atm
+import sys
+from ADC import *
+
+'''
+@author: 		Joost van den Born
+@Description: 	Calculates differential refraction (dispersion) between stars with a different 
+				effective temperature. Similar to Gubler & Tytler 1998
+'''
+
+def blackbody(l, T):
+	'''
+		Calculates the intensity output of a star with effective temperature T at wavelength l
+	'''
+	h = 6.626070150e-34	# J*s 
+	k = 1.38064852e-23	# J/K
+	c = 299792458		# m/s
+	return 2 * np.pi * h * c**2 / l**5 * 1 / (np.exp( h * c / (l * k * T)) - 1)
+
+def effectiveWavelength(T, lmin, lmax):
+	'''
+		Calculates effective wavelength of a star from the effective temperature.
+	'''
+	lmin = lmin/1e6
+	lmax = lmax/1e6
+	_l = np.linspace(lmin, lmax, 1000)
+	_I = blackbody(_l, T)
+	return np.sum(_I*_l)/np.sum(_I)*1e6
+
+if __name__ == '__main__':
+	# Parameters at Cerro Armazones
+	T   = 279.65 	# k
+	P   = 71200  	# Pa
+	H   = 0.22
+	xc  = 450   	# ppm
+	lat = -24.5983 	# degrees
+	h   = 3064  	# m
+	l1  = 1.49 		# micron
+	l2  = 1.78 		# micron
+	zenith = 60 	# degrees
+	offset = 0#.285/2		# degrees, ADC offset to simulate finite step size
+
+	S2temps = [2500, 4000, 5500, 7000, 8500]
+	S1temps = np.linspace(2500, 1e4, 100)
+
+
+	########################################################################################
+	########################################################################################
+	# 1. We look at the differential refraction through the atmosphere, without an ADC
+	# at = atm.dispersion(lat, h)
+	at = atm.refraction(lat, h)
+	at.setReducedHeight(P, at.rho(P, T, H, xc) )
+
+
+	# Calculating the differential atmospheric refraction for stars of different temperatures
+	data = []
+	for T2 in S2temps:
+		# Properties of star 2
+		l_eff_2 = effectiveWavelength(T2, l1, l2)
+		n2 = at.n_tph(l_eff_2, T, P, H, xc)
+		R2 = at.cassini(n2, zenith)
+
+		# Calculating the differences with star 1
+		dataEntry = []
+		for T1 in S1temps:
+			l_eff_1 = effectiveWavelength(T1, l1, l2)
+			n1 		= at.n_tph(l_eff_1, T, P, H, xc)
+			R1 		= at.cassini(n1, zenith)
+			diffRef = (R1 - R2) * 3.6e9
+			dataEntry.append(diffRef)
+		data.append(dataEntry)
+
+	# Plotting the results
+	plt.rc('font', family='serif')
+	plt.rc('axes', labelsize='large')
+	fig = plt.figure(figsize=(9,6))
+	ax = fig.add_subplot(111)
+	ax2 = ax.twiny()
+	ax.grid(True, alpha=0.3)
+	linestylesList = [':', '-.', '--', '-', (0, (3, 1, 1, 1, 1, 1))]
+	for i in range(len(data)):
+		ax.plot(S1temps, data[i], label=str(S2temps[i]) + ' K', linestyle=linestylesList[i], color='k', linewidth=2)
+	
+	S1temp_ticks = ax.get_xticks()
+	S1temp_tick_loc = ax.get_xticks()
+	
+	l_eff_ticks = []
+	for T1 in S1temp_ticks:
+		l_eff_ticks.append(round( (effectiveWavelength(T1, l1, l2)*1000-1620),1))
+	
+	ax2.set_xticklabels(l_eff_ticks)
+	ax2.set_xlim(np.min(S1temps), np.max(S1temps))
+	ax2.set_xlabel('$\lambda_{eff} - 1620$ (nm) ', fontsize=18)
+	ax2.tick_params(axis ='both', direction='in', labelsize=16) 
+
+	ax.set_xlabel('Temperature of Star 1 (K)', fontsize=18)
+	ax.set_ylabel('$\Delta$R ($\mu$as)', fontsize=18)
+	ax.legend(fontsize=16, labelspacing=0.3)
+	ax.set_xlim(np.min(S1temps), np.max(S1temps))
+	ax.tick_params(axis ='both', direction='in', labelsize=16) 
+	# plt.show()
+
+	########################################################################################
+	########################################################################################
+	
+
+
+	########################################################################################
+	########################################################################################
+	# 2. We look now at the differential refraction through the atmosphere, including an ADC
+	at = atm.refraction(lat, h)
+	at.setReducedHeight(P, at.rho(P, T, H, xc) )
+
+	adc = ADC()
+
+	# Conditions at the band cut-on and cut-off wavelength
+	n_l1 = at.n_tph(l1, T, P, H, xc)
+	n_l2 = at.n_tph(l2, T, P, H, xc)
+
+
+	atd = atm.dispersion(lat, h)
+	atd.setReducedHeight(P, at.rho(P, T, H, xc) )
+	# Optimize ADC for broad band edges at zenith distance.
+	atm_disp_full_band = atd.cassini(n_l1, n_l2, zenith)
+	print(atm_disp_full_band*3.6e6)
+	adcOptimumAngle = adc.optimize(atm_disp_full_band*3.6e6, l1, l2) 
+	adc.setOrientation(adcOptimumAngle+offset)
+
+
+	# Calculating the differential atmospheric refraction for stars of different temperatures
+	data = []
+	for T2 in S2temps:
+		# Properties of star 2
+		l_eff_2 = effectiveWavelength(T2, l1, l2)
+		n2 = at.n_tph(l_eff_2, T, P, H, xc)
+		R2_atm = at.cassini(n2, zenith)
+		R2_adc = adc.calculateDispersion(l_eff_2, l1, refraction=True)
+
+		# Calculating the differences with star 1
+		dataEntry = []
+		for T1 in S1temps:
+			l_eff_1 = effectiveWavelength(T1, l1, l2)
+			n1 		= at.n_tph(l_eff_1, T, P, H, xc)
+			R1_atm	= at.cassini(n1, zenith)
+			R1_adc 	= adc.calculateDispersion(l_eff_1, l1, refraction=True)
+			diffRef = ( (R1_adc[1] - R2_adc[1]) + ((R1_atm - R2_atm) * 3.6e6 ) ) * 1e3
+			dataEntry.append(diffRef)
+		data.append(dataEntry)
+
+	# Plotting the results
+	plt.rc('font', family='serif')
+	plt.rc('axes', labelsize='large')
+	fig = plt.figure(figsize=(9,6))
+	ax = fig.add_subplot(111)
+	ax2 = ax.twiny()
+	ax.grid(True, alpha=0.3)
+	linestylesList = [':', '-.', '--', '-', (0, (3, 1, 1, 1, 1, 1))]
+	for i in range(len(data)):
+		ax.plot(S1temps, data[i], label=str(S2temps[i]) + ' K', linestyle=linestylesList[i], color='k', linewidth=2)
+	
+	S1temp_ticks = ax.get_xticks()
+	S1temp_tick_loc = ax.get_xticks()
+	
+	l_eff_ticks = []
+	for T1 in S1temp_ticks:
+		l_eff_ticks.append(round( (effectiveWavelength(T1, l1, l2)*1000-1620),1))
+	
+	ax2.set_xticklabels(l_eff_ticks)
+	ax2.set_xlim(np.min(S1temps), np.max(S1temps))
+	ax2.set_xlabel('$\lambda_{eff} - 1620$ (nm) ', fontsize=18)
+	ax2.tick_params(axis ='both', direction='in', labelsize=16) 
+
+	ax.set_xlabel('Temperature of Star 1 (K)', fontsize=18)
+	ax.set_ylabel('$\Delta$R ($\mu$as)', fontsize=18)
+	ax.legend(fontsize=16, labelspacing=0.3, loc='lower right')
+	ax.set_xlim(np.min(S1temps), np.max(S1temps))
+	ax.tick_params(axis ='both', direction='in', labelsize=16) 
+	plt.show()
+
+	########################################################################################
+	########################################################################################
+
+
+
+	########################################################################################
+	########################################################################################
+	# 3. Now we look at the differential dispersion, without an ADC
+	at = atm.dispersion(lat, h)
+	at.setReducedHeight(P, at.rho(P, T, H, xc) )
+
+	# Conditions at the band cut-on and cut-off wavelength
+	n_l1 = at.n_tph(l1, T, P, H, xc)
+	n_l2 = at.n_tph(l2, T, P, H, xc)
+
+	# Optimize ADC for broad band edges at zenith distance.
+	atm_disp_full_band = at.cassini(n_l1, n_l2, zenith)
+	print(atm_disp_full_band*3.6e6)
+
+
+	# Check that we can use dispersion? If not, how to calculate refraction...?
+	# We plot below the difference in atmospheric dispersion of the two stars, incl. an ADC.
+	data = []
+	for T2 in S2temps:
+		l_eff_2 = effectiveWavelength(T2, l1, l2)
+		n2 = at.n_tph(l_eff_2, T, P, H, xc)
+		dataEntry = []
+		for T1 in S1temps:
+			l_eff_1 = effectiveWavelength(T1, l1, l2)
+			n1 		= at.n_tph(l_eff_1, T, P, H, xc)
+			atm_disp = at.cassini(n1, n2, zenith) * 3.6e9
+			dataEntry.append(atm_disp)
+		data.append(dataEntry)
+
+	plt.rc('font', family='serif')
+	plt.rc('axes', labelsize='large')
+	fig2 = plt.figure(figsize=(9,6))
+	ax = fig2.add_subplot(111)
+	ax2 = ax.twiny()
+	ax.grid(True, alpha=0.3)
+	linestylesList = [':', '-.', '--', '-', (0, (3, 1, 1, 1, 1, 1))]
+	for i in range(len(data)):
+		ax.plot(S1temps, data[i], label=str(S2temps[i]) + ' K', linestyle=linestylesList[i], color='k', linewidth=2)
+	
+	S1temp_ticks = ax.get_xticks()
+	S1temp_tick_loc = ax.get_xticks()
+	
+	l_eff_ticks = []
+	for T1 in S1temp_ticks:
+		l_eff_ticks.append(round( (effectiveWavelength(T1, l1, l2)*1000-1620),1))
+	
+	# ax2.plot(S1temps, data[0], alpha=0)
+	ax2.set_xticklabels(l_eff_ticks)
+	ax2.set_xlim(np.min(S1temps), np.max(S1temps))
+	ax2.set_xlabel('$\lambda_{eff} - 1620$ (nm) ', fontsize=18)
+	ax2.tick_params(axis ='both', direction='in', labelsize=16) 
+
+	ax.set_xlabel('Temperature of Star 1 (K)', fontsize=18)
+	ax.set_ylabel('$\Delta(\Delta$R)  ($\mu$as), excl. ADC', fontsize=18)
+	ax.legend(fontsize=16, labelspacing=0.3, loc='lower right')
+	ax.set_xlim(np.min(S1temps), np.max(S1temps))
+	ax.tick_params(axis ='both', direction='in', labelsize=16) 
+	# plt.show()
+
+	########################################################################################
+	########################################################################################
+	########################################################################################
+
+
+
+
+	########################################################################################
+	########################################################################################
+	# 4. Now we add the ADC and look again at the differential refraction, but on the detector
+	at = atm.dispersion(lat, h)
+	at.setReducedHeight(P, at.rho(P, T, H, xc) )
+	adc = ADC()
+
+	# Conditions at the band cut-on and cut-off wavelength
+	n_l1 = at.n_tph(l1, T, P, H, xc)
+	n_l2 = at.n_tph(l2, T, P, H, xc)
+
+	# Optimize ADC for broad band edges at zenith distance.
+	atm_disp_full_band = at.cassini(n_l1, n_l2, zenith)
+	print(atm_disp_full_band*3.6e6)
+	adcOptimumAngle = adc.optimize(atm_disp_full_band*3.6e6, l1, l2) 
+	adc.setOrientation(adcOptimumAngle+offset)
+
+
+	# Check that we can use dispersion? If not, how to calculate refraction...?
+	# We plot below the difference in atmospheric dispersion of the two stars, incl. an ADC.
+	data = []
+	for T2 in S2temps:
+		l_eff_2 	= effectiveWavelength(T2, l1, l2)
+		n2 			= at.n_tph(l_eff_2, T, P, H, xc)
+		atm_disp2 	= at.cassini(n_l1, n2, zenith) * 3.6e6
+		adc_disp2 	= adc.calculateDispersion(l1, l_eff_2)[1]
+
+		dataEntry = []
+		for T1 in S1temps:
+			l_eff_1 	= effectiveWavelength(T1, l1, l2)
+			n1 			= at.n_tph(l_eff_1, T, P, H, xc)
+			atm_disp1 	= at.cassini(n_l1, n1, zenith) * 3.6e6
+			adc_disp1 	= adc.calculateDispersion(l1,l_eff_1)[1]
+			diffDisp 	= ((adc_disp1 - adc_disp2) - (atm_disp1 - atm_disp2)) * 1e3
+			dataEntry.append(diffDisp)
+		data.append(dataEntry)
+
+	plt.rc('font', family='serif')
+	plt.rc('axes', labelsize='large')
+	fig2 = plt.figure(figsize=(9,6))
+	ax = fig2.add_subplot(111)
+	ax2 = ax.twiny()
+	ax.grid(True, alpha=0.3)
+	linestylesList = [':', '-.', '--', '-', (0, (3, 1, 1, 1, 1, 1))]
+	for i in range(len(data)):
+		ax.plot(S1temps, data[i], label=str(S2temps[i]) + ' K', linestyle=linestylesList[i], color='k', linewidth=2)
+	
+	S1temp_ticks = ax.get_xticks()
+	S1temp_tick_loc = ax.get_xticks()
+	
+	l_eff_ticks = []
+	for T in S1temp_ticks:
+		l_eff_ticks.append(round( (effectiveWavelength(T, l1, l2)*1000-1620),1))
+	
+	# ax2.plot(S1temps, data[0], alpha=0)
+	ax2.set_xticklabels(l_eff_ticks)
+	ax2.set_xlim(np.min(S1temps), np.max(S1temps))
+	ax2.set_xlabel('$\lambda_{eff} - 1620$ (nm) ', fontsize=18)
+	ax2.tick_params(axis ='both', direction='in', labelsize=16) 
+
+	ax.set_xlabel('Temperature of Star 1 (K)', fontsize=18)
+	ax.set_ylabel('$\Delta(\Delta$R)  ($\mu$as), incl. ADC', fontsize=18)
+	ax.legend(fontsize=16, labelspacing=0.3, loc='lower right')
+	ax.set_xlim(np.min(S1temps), np.max(S1temps))
+	ax.tick_params(axis ='both', direction='in', labelsize=16) 
+	plt.show()
+
+	########################################################################################
+	########################################################################################
+	########################################################################################

scripts/uncertainty_analysis.py

@@ -0,0 +1,353 @@
+#!/usr/bin/env Python
+from __future__ import division
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import host_subplot
+from AstroAtmosphere import *
+from tabulate import tabulate
+
+
+
+'''
+@ Author:		Joost van den Born 
+@ Contact:		born@astron.nl
+@ Description:	An exploration of the parameter space for the Ciddor-Cassini model.
+				The results here are given in table 3 and Appendix A of the paper.
+				
+'''
+
+
+if __name__ == '__main__':
+
+	plt.rc('font', family='serif')
+	plt.rc('axes', labelsize='large')
+
+	at = Observatory()
+	
+	pdrFigs 		= False
+	errorSpace 		= False
+	generate_images = True
+
+
+	T   = 279.65 	# k
+	P   = 71200  	# Pa
+	H   = 0.22
+	xc  = 450   	# ppm
+	lat = -24.5983 	# degrees
+	h   = 3064  	# m
+	l1  = 1.49 		# micron
+	l2  = 1.78 		# micron
+	zenith = 30 	# degrees
+
+	rho = at.rho(p=P, T=T, RH=H, xc=xc)
+	disp = dispersion(lat, h)
+	disp.setReducedHeight(P, rho)
+
+	# Change dl, dT, dP, dRH, dCO2 to the values of interest.
+	n1 = at.n_tph(l=1.49, T=279.65, p=71200, RH=0.22, xc=450)
+	n2 = at.n_tph(l=1.78, T=279.65, p=71200, RH=0.22, xc=450)
+	R1ca = at.dn_tph(l=1.49, T=279.65, p=71200, RH=0.22, xc=450,
+		dl=0.001, dT=0.4, dP=10, dRH=0.01, dCO2=20)
+	R2ca = at.dn_tph(l=1.78, T=279.65, p=71200, RH=0.22, xc=450,
+		dl=0.001, dT=0.4, dP=10, dRH=0.01, dCO2=20)
+	
+
+	# print (n1, at.dn_tph(l=1.49, T=288.0, p=101325, RH=0.0, xc=450, dl=0.0, dT=0.5, dP=0.0, dRH=0.0, dCO2=0) )
+	# print(disp.cassiniError( zenith, l1, l2, T, P, H, xc, dl1=0.001, dl2=0.001, dT=0.4, dP=10, dRH=0.01, dCO2=20 )*3.6e9)
+
+	zenith_true = 30 # degrees
+	dispValue = disp.cassini(n1, n2, zenith_true)
+	defaultDispCA = disp.cassini(n1, n2, zenith_true)*3.6e6
+	errorDispValue = disp.cassiniError(zenith, l1, l2, T, P, H, xc, dl1=0.001, dl2=0.001, dT=0.4, dP=10, dRH=0.01, dCO2=20 )
+	print ('Dispersion between 1.49 and 1.78 micron at z = %i deg: %.3f +/- %.3f mas' %(zenith_true, dispValue*3.6e6, errorDispValue*3.6e6))
+	
+
+	R1s  = at.n_tph(l=1.49, T=288.15, p=101325, RH=0.0, xc=450)
+	R2s  = at.n_tph(l=1.78, T=288.15, p=101325, RH=0.0, xc=450)
+
+	# Varying the parameters in the non-errror formula:
+	if pdrFigs:
+		Trange 			= np.linspace(262, 298, 1000)
+		Prange 			= np.linspace(69000, 73000, 1000)
+		Hrange 			= np.linspace(0.1, 0.3, 1000)
+
+		fig = plt.figure()
+
+		# Variation in T
+		dispListS, dispListCA = [], []
+		for t in Trange:
+			R1ca = at.n_tph(l=1.49, T=t, p=71200, RH=0.22, xc=450)
+			R2ca = at.n_tph(l=1.78, T=t, p=71200, RH=0.22, xc=450)
+			dispValueS = disp.cassini(R1s, R2s, zenith_true) * 3.6e6
+			dispValueCA = disp.cassini(R1ca, R2ca, zenith_true) * 3.6e6
+			dispListS.append(dispValueS)
+			dispListCA.append(dispValueCA)
+		ax1 = fig.add_subplot(131)
+		ax1.grid(True, alpha=0.3)
+		ax1.plot(Trange, np.abs(dispListCA-defaultDispCA), 'r-', label='Cerro Armazones (p=712 mbar, T=6.5C, H=22%)')
+		ax1.set_xlabel('Temperature (K)')
+		ax1.set_ylabel('Dispersion (mas)')
+		ax1.tick_params(axis ='both', direction='in') 
+
+		# Variation in P
+		dispListS, dispListCA = [], []
+		for p in Prange:
+			R1ca = at.n_tph(l=1.49, T=279.65, p=p, RH=0.22, xc=450)
+			R2ca = at.n_tph(l=1.78, T=279.65, p=p, RH=0.22, xc=450)
+			dispValueS = disp.cassini(R1s, R2s, zenith_true) * 3.6e6
+			dispValueCA = disp.cassini(R1ca, R2ca, zenith_true) * 3.6e6
+			dispListS.append(dispValueS)
+			dispListCA.append(dispValueCA)
+		ax2 = fig.add_subplot(132)
+		ax2.grid(True, alpha=0.3)
+		ax2.plot(Prange/100, np.abs(dispListCA-defaultDispCA), 'r-', label='Cerro Armazones (p=712 mbar, T=6.5C, H=22%)')
+		ax2.set_xlabel('Pressure (mbar)')
+		ax2.set_ylabel('Dispersion (mas)')
+		ax2.tick_params(axis ='both', direction='in') 
+
+		# Variation in T
+		dispListS, dispListCA = [], []
+		for h in Hrange:
+			R1ca = at.n_tph(l=1.49, T=279.65, p=71200, RH=h, xc=450)
+			R2ca = at.n_tph(l=1.78, T=279.65, p=71200, RH=h, xc=450)
+			dispValueS = disp.cassini(R1s, R2s, zenith_true) * 3.6e6
+			dispValueCA = disp.cassini(R1ca, R2ca, zenith_true) * 3.6e6
+			dispListS.append(dispValueS)
+			dispListCA.append(dispValueCA)
+		ax3 = fig.add_subplot(133)
+		ax3.grid(True, alpha=0.3)
+		ax3.plot(Hrange*100, np.abs(dispListCA-defaultDispCA), 'r-', label='Cerro Armazones (p=712 mbar, T=6.5C, H=22%)')
+		ax3.set_xlabel('Humidity (%)')
+		ax3.set_ylabel('Dispersion (mas)')
+		ax3.tick_params(axis ='both', direction='in') 
+
+		plt.tight_layout()
+		plt.show()
+
+
+	# Plotting variations in parameters with corresponding errors.	
+	if errorSpace:	
+		Trange 		= np.linspace(-50, 50, 100)
+		Prange 		= np.linspace(-50, 50, 10)
+		Hrange 		= np.linspace(-0.025, 0.025, 10)
+
+		fig = plt.figure()
+		# Variation in T
+		errorRange = [0.1, 0.2, 0.3, 1.0]
+		errorRange = errorRange[::-1]
+		
+		ax1 = fig.add_subplot(111)
+		ax1.grid(True, alpha=0.3)
+
+		for e in range(len(errorRange)):
+			dispList = []
+			errorList = []
+			for t in Trange:
+				R1n = at.n_tph(l=1.49, T=279.65+t, p=71200, RH=0.22, xc=450)
+				R2n = at.n_tph(l=1.78, T=279.65+t, p=71200, RH=0.22, xc=450)
+				dispValue  = disp.cassini(R1n,  R2n,  zenith_true) * 3.6e6	
+				errorValue = disp.cassiniError(zenith_true, l1, l2, T, P, H, xc, dl1=0.0, dl2=0.0, dT=errorRange[e], dP=0, dRH=0, dCO2=0 )  * 3.6e6
+				dispList.append(dispValue)
+				errorList.append(errorValue)
+			dispList = np.asarray(dispList)
+			errorList = np.asarray(errorList)
+			ax1.fill_between(Trange+279.65, dispList-errorList, dispList+errorList, label='$\delta p$ = %s $^{\circ}$C' %(errorRange[e]))
+
+		ax1.plot(Trange+279.65, dispList, 'k-')
+		ax1.set_xlabel('Temperature (K)')
+		ax1.set_ylabel('Dispersion (mas)')
+		ax1.set_xlim([np.min(Trange)+279.65,np.max(Trange)+279.65])
+		ax1.set_ylim([19,32])
+		handles, labels = ax1.get_legend_handles_labels()
+		ax1.legend(handles[::-1], labels[::-1], title='Error in measurements')
+		plt.tight_layout()
+		plt.show()
+
+	
+	# Varying the parameters in the uncertainty formula:
+	if generate_images:
+		################################################
+		# Listing some values and the results in a table
+		pList = [101325, 101325, 71200, 71200, 71200]
+		TList = [288.15, 288.15, 279.65, 279.65, 279.65]
+		HList = [0, 0, 0.22, 0.22, 0.6]
+		CO2List = [450] * len(pList)
+		lList = [1.35] * len(pList)
+		dpList = [0, 20, 20, 10, 20]
+		dTlist = [0, 0.2, 0.2, 0.1, 0.2]
+		dRHlist = [0, 0.02, 0.02, 0.01, 0.02]
+		dCO2list = [0, 20, 20, 10, 20]
+		dlList = [0, 0.001, 0.001, 0.001, 0.001]
+
+		refractionList, sigmaList = [], []
+		for i in range(len(pList)):
+			n = at.n_tph(l=lList[i], T=TList[i], p=pList[i], RH=HList[i], xc=CO2List[i])
+			s = at.dn_tph(l=lList[i], T=TList[i], p=pList[i], RH=HList[i], xc=CO2List[i],
+				dl=dlList[i], dT=dTlist[i], dP=dpList[i], dRH=dRHlist[i], dCO2=dCO2list[i])
+			refractionList.append((n-1)*1e8)
+			sigmaList.append(s*1e8)
+
+		tableData = []
+		for i in range(len(pList)):
+			tableData.append([pList[i], TList[i], HList[i], CO2List[i], lList[i], dpList[i], dTlist[i], dRHlist[i], dCO2list[i], dlList[i], refractionList[i], sigmaList[i]])
+		print('REFRACTION: ')
+		print (tabulate(tableData, headers = ['p', 'T', 'H', 'x_c', 'l', 'dp', 'dT', 'dH', 'dx_c', 'dl', 'n', 's']))
+		
+		#################################################
+		print ('\n\n\n')
+
+		################################################
+		# Listing some values and the results in a table
+		pList = [101325, 101325, 101325, 101325, 101325]
+		TList = [288.15, 288.15, 288.15, 288.15, 288.15]
+		HList = [0.0, 0.0, 0.0, 0.0, 0.0]
+		CO2List = [450] * len(pList)
+		# lList = [1.35] * len(pList)
+		dpList = [0, 10, 10, 20, 100]
+		dTlist = [0, 0.1, 0.2, 0.2, 1.0]
+		dRHlist = [0, 0.01, 0.01, 0.02, 0.05]
+		dCO2list = [0, 20, 20, 20, 20]
+		# dlList = [0, 0.001, 0.001, 0.001, 0.001] # Uncomment this and comment the next 
+												   # line to check the impact with uncertainty in wavelength
+		dlList = [0, 0, 0, 0, 0]
+		dzList = [0] * len(pList)
+
+		dispersionList, errorList = [], []
+		for i in range(len(pList)):
+			zenith_true = 30 # degrees
+			n1  = at.n_tph(l=1.49, T=TList[i], p=pList[i], RH=HList[i], xc=CO2List[i])
+			n2  = at.n_tph(l=1.78, T=TList[i], p=pList[i], RH=HList[i], xc=CO2List[i])
+			dn1 = at.dn_tph(l=1.49, T=TList[i], p=pList[i], RH=HList[i], xc=CO2List[i],
+				dl=dlList[i], dT=dTlist[i], dP=dpList[i], dRH=dRHlist[i], dCO2=dCO2list[i])
+			dn2 = at.dn_tph(l=1.78, T=TList[i], p=pList[i], RH=HList[i], xc=CO2List[i],
+				dl=dlList[i], dT=dTlist[i], dP=dpList[i], dRH=dRHlist[i], dCO2=dCO2list[i])
+			dispValue = disp.cassini(n1, n2, zenith_true) * 3.6e6
+			errorDispValue = disp.cassiniError( zenith_true, 1.49, 1.78, TList[i], pList[i], HList[i], CO2List[i], 
+				dl1=dlList[i], dl2=dlList[i], dT=dTlist[i], dP=dpList[i], dRH=dRHlist[i], dCO2=dCO2list[i], dz=dzList[i] ) * 3.6e9
+			dispersionList.append(dispValue)
+			errorList.append(errorDispValue)
+
+		tableData = []
+		for i in range(len(pList)):
+			tableData.append([pList[i], TList[i], HList[i], CO2List[i], dpList[i], dTlist[i], dRHlist[i], dCO2list[i], dlList[i], dispersionList[i], errorList[i]])
+		print('DISPERSION: ')
+		print (tabulate(tableData, headers = ['p', 'T', 'H', 'x_c', 'dp', 'dT', 'dH', 'dx_c', 'dl', 'dispersion', 'error']))
+		
+		# Variation of the uncertainty in P
+		pRange = np.linspace(0, 100, 1000)
+		errorList = []
+		for p in pRange:
+			errorDispValue = disp.cassiniError( zenith_true, 1.49, 1.78, T=279.65, p=71200, RH=0.22, xc=450, 
+				dl1=0.001, dl2=0.001, dT=0.2, dP=p, dRH=0.01, dCO2=20, dz=0 ) * 3.6e9
+			errorList.append(errorDispValue)
+
+		fig = plt.figure(figsize=(4,3))
+		ax = fig.add_subplot(111)
+		ax.grid(True, alpha=0.3)
+		ax.plot(pRange, errorList, 'k-')
+		ax.set_xlabel('$\sigma_p$ (Pa)')
+		ax.set_ylabel('$\sigma_{\Delta R}$ ($\mu$as)')
+		ax.set_xlim([np.min(pRange), np.max(pRange)])
+		ax.set_ylim([np.min(errorList)-0.05*(np.max(errorList) - np.min(errorList)), np.max(errorList)])
+		ax.tick_params(axis ='both', direction='in') 
+		plt.tight_layout()
+		
+
+		# Variation of the uncertainty in T
+		TRange = np.linspace(0.0, 1, 1000)
+		errorList = []
+		for T in TRange:
+			errorDispValue = disp.cassiniError( zenith_true, 1.49, 1.78, T=279.65, p=71200, RH=0.22, xc=450, 
+				dl1=0.001, dl2=0.001, dT=T, dP=10, dRH=0.01, dCO2=20, dz=0 ) * 3.6e9
+			errorList.append(errorDispValue)
+
+		fig = plt.figure(figsize=(4,3))
+		ax = fig.add_subplot(111)
+		ax.grid(True, alpha=0.3)
+		ax.plot(TRange, errorList, 'k-')
+		ax.set_xlabel('$\sigma_T$ (K)')
+		ax.set_ylabel('$\sigma_{\Delta R}$ ($\mu$as)')
+		ax.set_xlim([np.min(TRange), np.max(TRange)])
+		ax.set_ylim([np.min(errorList)-0.05*(np.max(errorList) - np.min(errorList)), np.max(errorList)])
+		ax.tick_params(axis ='both', direction='in')
+		plt.tight_layout()
+
+		# Variation of the uncertainty in H
+		HRange = np.linspace(0.000, 0.05, 1000)
+		errorList = []
+		for H in HRange:
+			errorDispValue = disp.cassiniError( zenith_true, 1.49, 1.78, T=279.65, p=71200, RH=0.22, xc=450, 
+				dl1=0.001, dl2=0.001, dT=0.2, dP=10, dRH=H, dCO2=20, dz=0 ) * 3.6e9
+			errorList.append(errorDispValue)
+
+		fig = plt.figure(figsize=(4,3))
+		ax = fig.add_subplot(111)
+		ax.grid(True, alpha=0.3)
+		ax.plot(HRange, errorList, 'k-')
+		ax.set_xlabel('$\sigma_H$')
+		ax.set_ylabel('$\sigma_{\Delta R}$ ($\mu$as)')
+		ax.set_xlim([np.min(HRange), np.max(HRange)])
+		ax.set_ylim([np.min(errorList)-0.05*(np.max(errorList) - np.min(errorList)), np.max(errorList)])
+		ax.tick_params(axis ='both', direction='in')
+		plt.tight_layout()
+
+
+		# Variation of the uncertainty in xc
+		CO2Range = np.linspace(0, 200, 1000)
+		errorList = []
+		for CO2 in CO2Range:
+			errorDispValue = disp.cassiniError( zenith_true, 1.49, 1.78, T=279.65, p=71200, RH=0.22, xc=450, 
+				dl1=0.001, dl2=0.001, dT=0.2, dP=10, dRH=0.01, dCO2=CO2, dz=0 ) * 3.6e9
+			errorList.append(errorDispValue)
+
+		fig = plt.figure(figsize=(4,3))
+		ax = fig.add_subplot(111)
+		ax.grid(True, alpha=0.3)
+		ax.plot(CO2Range, errorList, 'k-')
+		ax.set_xlabel('$\sigma_{x_c}$ (ppm)')
+		ax.set_ylabel('$\sigma_{\Delta R}$ ($\mu$as)')
+		ax.set_xlim([np.min(CO2Range), np.max(CO2Range)])
+		ax.set_ylim([np.min(errorList)-0.05*(np.max(errorList) - np.min(errorList)), np.max(errorList)])
+		ax.tick_params(axis ='both', direction='in')
+		plt.tight_layout()
+
+		# Variation of the uncertainty in wavelength
+		lRange = np.linspace(0.0, 0.001, 1000)
+		errorList = []
+		for l in lRange:
+			errorDispValue = disp.cassiniError( zenith_true, 1.49, 1.78, T=279.65, p=71200, RH=0.22, xc=450, 
+				dl1=l, dl2=l, dT=0.2, dP=10, dRH=0.01, dCO2=20, dz=0 ) * 3.6e9
+			errorList.append(errorDispValue)
+
+		fig = plt.figure(figsize=(4,3))
+		ax = fig.add_subplot(111)
+		ax.grid(True, alpha=0.3)
+		ax.plot(lRange*1000, errorList, 'k-')
+		ax.set_xlabel('$\sigma_{\lambda}$ (nm)')
+		ax.set_ylabel('$\sigma_{\Delta R}$ ($\mu$as)')
+		ax.set_xlim([np.min(lRange)*1000, np.max(lRange)*1000])
+		ax.set_ylim([np.min(errorList)-0.05*(np.max(errorList) - np.min(errorList)), np.max(errorList)])
+		ax.tick_params(axis ='both', direction='in')
+		plt.tight_layout()
+		
+		
+		# Variation of the uncertainty in zenith
+		zRange = np.linspace(0, 60/3600, 1000)
+		errorList = []
+		for z in zRange:
+			errorDispValue = disp.cassiniError( zenith_true, 1.49, 1.78, T=279.65, p=71200, RH=0.22, xc=450, 
+				dl1=0.001, dl2=0.001, dT=0.2, dP=10, dRH=0.01, dCO2=20, dz=z ) * 3.6e9
+			errorList.append(errorDispValue)
+
+		fig = plt.figure(figsize=(4,3))
+		ax = fig.add_subplot(111)
+		ax.grid(True, alpha=0.3)
+		ax.plot(zRange*3600, errorList, 'k-')
+		ax.set_xlim([np.min(zRange)*3600, np.max(zRange)*3600])
+		ax.set_ylim([np.min(errorList)-0.05*(np.max(errorList) - np.min(errorList)), np.max(errorList)])
+		ax.set_xlabel('$\sigma_z$ (arcseconds)')
+		ax.set_ylabel('$\sigma_{\Delta R}$ ($\mu$as)')
+		ax.tick_params(axis ='both', direction='in')
+		plt.tight_layout()
+		plt.show()
+		

scripts/uncertainty_analysis_verification.py

@@ -0,0 +1,85 @@
+import numpy as np
+from AstroAtmosphere import *
+
+'''
+	Date: 02 May 2019
+	Author: Joost van den Born
+	Description: This script compares the uncertainty propagation results to a Monte Carlo method.
+'''
+
+
+# Parameters at Cerro Armazones
+T   = 288.15 	# k
+P   = 101325  	# Pa
+H   = 0.0
+xc  = 450   	# ppm
+lat = -24.5983 	# degrees
+h   = 3064  	# m
+l1  = 1.49 		# micron
+l2  = 1.78 		# micron
+zenith = 45 	# degrees
+
+dl1 = 0
+dl2 = 0
+dT  = 0
+dP  = 2
+dRH = 0
+dCO2= 0
+dz = 1/3600 # degrees
+
+# Initializing dispersion model
+at  = Observatory()
+
+dPlist = [2, 20, 200, 0,    0,   0, 0,     0,    0,   0,  0, 0,   2,     20,   200]
+dTlist = [0, 0,  0,   0.02, 0.2, 2, 0,     0,    0,   0,  0, 0,   0.02,  0.2,  2  ]
+dRHlist= [0, 0,  0,   0,    0,   0, 0.002, 0.02, 0.2, 0,  0, 0,   0.002, 0.02, 0.2]
+dCO2list=[0, 0,  0,   0,    0,   0, 0,     0,    0,   2, 20, 200, 2,     20,   200]
+
+for i in range(len(dPlist)):
+	nerr = at.dn_tph(l1, T, P, H, xc, dl=0, dT=dTlist[i], dP=dPlist[i], dRH=dRHlist[i], dCO2=dCO2list[i])
+	print(nerr, dPlist[i], dTlist[i], dRHlist[i], dCO2list[i])
+
+
+
+
+
+
+disp = dispersion(lat, h)
+rho_obs = at.rho(P, T, H, xc)
+disp.setReducedHeight(P, rho_obs)
+zenith_true = 10
+
+# We do a Monte Carlo error analysis on the dispersion.
+print('MONTE-CARLO UNCERTAINTY ANALYSIS')
+print('-'*35)
+
+for j in range(len(dPlist)):
+    sigmaP 	 = dPlist[j] # mbar
+    sigmaT 	 = dTlist[j] # K
+    sigmaH 	 = dRHlist[j]
+    sigmaCO2 = dCO2list[j] # ppm
+
+    N = 10000
+    l1sample = np.random.normal(l1, 0, N)
+    l2sample = np.random.normal(l2, 0, N)
+    Tsample   = np.random.normal(T,   sigmaT,   N)
+    Psample   = np.random.normal(P,   sigmaP,   N)
+    Hsample   = np.random.normal(H,   sigmaH,   N)
+    CO2sample = np.random.normal(xc, sigmaCO2, N)
+    dzsample  = np.random.normal(zenith_true, dz, N)
+    
+    results = []
+    nList = []
+
+    for i in range(N):
+        # Calculate the refraction at l1 and l2
+        R1ca2 = at.n_tph(l=l1sample[i], T=Tsample[i], p=Psample[i], RH=Hsample[i], xc=CO2sample[i])
+        R2ca2 = at.n_tph(l=l2sample[i], T=Tsample[i], p=Psample[i], RH=Hsample[i], xc=CO2sample[i])
+        dispersionValue = disp.cassini(R1ca2, R2ca2, zenith_true)
+        results.append(dispersionValue)
+        nList.append(R1ca2)
+    
+    cassErrMC = np.std(results)*3.6e9
+    nerrMC    = np.std(nList)
+    print(nerrMC, dPlist[j], dTlist[j], dRHlist[j], dCO2list[j])
+

setup.py

@@ -21,14 +21,14 @@ def setup_package():
     # write_version_py()
 
     setup(name='AstroAtmosphere',
-          version='0.1.2',
+          version='1.0',
           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/v_012/atmosphericmodel-v_012.tar.gz',
+          download_url='https://gitlab.astro-wise.org/micado/atmosphericmodel/-/archive/v_10/atmosphericmodel-v_10.tar.gz',
           license='GNU General Public License v3',
           data_files = [('',['LICENSE.txt'])],
           package_dir={'AstroAtmosphere': 'AstroAtmosphere'},