#!/usr/bin/env Python3
import numpy as np
from scipy import integrate
from pdb import set_trace

'''
@ Author:		Joost van den Born 
@ Contact:		born@astron.nl
@ Description:	This script offers a complete model of the atmospheric refraction
				and dispersion. It uses the Ciddor1996 paper to calculate the 
				index of refraction for a wide range of conditions, with corresponding
				error propagation if needed.
				The actual refraction is calculated by any of the following models:
				- Plane parallel atmosphere
				- Cassini spherical homogeneous atmosphere
				- Other methods described in Corbard et al., 2019.
'''


class Observatory:
	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
		self.dH 	= dH	# relative humidity between 0 and 1
		self.dCO2 	= dCO2	# in ppm
		self.dl 	= dl 	# in micron

	def n_as(self, l):
		'''
			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!
		'''
		_s = 1 / l
		_k0 = 238.0185 	# um^-2
		_k1 = 5792105  	# um^-2
		_k2 = 57.362		# um^-2
		_k3 = 167917		# um^-2
		return (_k1 / (_k0 - _s**2) + _k3 / (_k2 - _s**2)) / 1e8 + 1


	def n_axs(self, l, xc):
		'''
			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)
		'''
		_n_as = self.n_as(l)
		return ( (_n_as - 1) * (1 + 0.534e-6 * (xc - 450) ) ) + 1


	def n_ws(self, l):
		'''
			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
		'''
		_s  = 1 / l
		_cf = 1.022
		_w0 = 295.235 	# um^-2
		_w1 = 2.6422	# um^-2
		_w2 = -0.032380	# um^-4
		_w3 = 0.004028	# um^-6
		return (_cf * (_w0 + _w1 * _s**2 + _w2 * _s**4 + _w3 * _s**6) ) / 1e8 + 1


	def f(self, p, T):
		'''
			Enhancement factor of water vapour in air
			t in degrees celsius
			p in Pascal (Pa)
		'''
		_t = T - 273.15
		_a = 1.00062
		_b = 3.14e-8 	# Pa^-1
		_c = 5.6e-7		# deg C^-2
		return _a + _b * p + _c * _t**2


	def svp(self, T):
		'''
			Saturation vapour pressure in Pa.
		'''
		_A =  1.2378847e-5 # K^-2
		_B = -1.9121316e-2 # K^-2
		_C = 33.93711047
		_D = -6.3431645e3  # K^-2
		return np.exp(_A*T**2 + _B*T + _C + _D/T)


	def svp2(self, T):
		'''
			Calculation of svp from IAPWS formula, based on documentation by NIST toolbox.
		'''
		K1 =  1.16705214528e+03
		K2 = -7.24213167032e+05
		K3 = -1.70738469401e+01
		K4 =  1.20208247025e+04
		K5 = -3.23255503223e+06
		K6 =  1.49151086135e+01
		K7 = -4.82326573616e+03
		K8 =  4.05113405421e+05
		K9 = -2.38555575678e-01
		K10 = 6.50175348448e+02
		_W = T + K9 / (T - K10)
		_A = _W**2 + K1*_W + K2
		_B = K3*_W**2 + K4*_W + K5
		_C = K6*_W**2 + K7*_W + K8
		_X = -_B + np.sqrt(_B**2 - 4*_A*_C)
		return (2*_C / _X)**4 * 1e6


	def xw(self, p, T, RH):
		'''
			Molar fraction of water vapour in moist air.
		'''
		_f 		= self.f(p, T)
		_svp 	= self.svp(T)
		return _f * RH * _svp / p


	def Z(self, p, T, RH, xw=''):
		'''
			Compressibility of moist air
		'''
		_a0 =  1.58123e-6 # K Pa^-1
		_a1 = -2.9331e-8  # Pa^-1
		_a2 =  1.1043e-10 # K^-1 Pa^-1
		_b0 =  5.707e-6   # K Pa^-1
		_b1 = -2.051e-8   # Pa^-1
		_c0 =  1.9898e-4  # K Pa^-1
		_c1 = -2.376e-6   # Pa^-1
		_d  =  1.83e-11   # K^2 Pa^-2
		_e  = -0.765e-8   # K^2 Pa^-2
		_t  = T - 273.15
		if xw == '':
			_xw = self.xw(p, T, RH)
		else:
			_xw = xw
		Z = 1 - (p/T) * (_a0 + _a1*_t + _a2*_t**2 + (_b0 + _b1*_t)*_xw + (_c0 + _c1*_t)*_xw**2) + (p/T)**2 * (_d + _e*_xw**2)
		return Z


	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
		'''
		_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
		_R  = 8.314510										# J mol^-1 K^-1, gas constant
		if xw == '':
			_xw = self.xw(p, T, RH)
		else: 
			_xw = xw
		_Z  = self.Z(p, T, RH, xw=_xw)
		return ( ( p * _Ma ) / (_Z * _R * T) ) * ( 1 - _xw * (1 - _Mw/_Ma)) 


	def n_tph(self, l, T, p, RH, xc):
		'''
			Calculates the refraction of the atmosphere for any conditions.
		'''
		_n_a 	 = self.n_axs(l, xc)
		_n_w 	 = self.n_ws(l)
		_La 	 = (_n_a**2 - 1) / (_n_a**2 + 2)
		_Lw 	 = (_n_w**2 - 1) / (_n_w**2 + 2)
		
		# Calculate the density of the dry air components, using the partial pressure of the dry air, i.e. p_dry = p_atm - p_water
		_rho_a 	 = self.rho(p=p-RH*self.svp(T), T=T, RH=0, xc=xc)				# Dry air component for the actual conditions.
		_rho_axs = self.rho(p=101325, T=288.15, RH=0, xc=xc) 					# Dry air component at standard conditions.
		
		if RH == 0:
			# nan is returned for _rho_w if we put in a (partial) pressure of 0
			_rho_w = 0
		else:
			# Partial pressure of water vapour is p_water = h*svp
			_rho_w 	 = self.rho(p=RH*self.svp(T), T=T, RH=RH, xc=xc)			# Moist air component for the actual conditions.
		# _rho_ws	 = self.rho(p=1333, T=293.15, RH=1/1.7564043320974863, xc=xc) 	# Moist air component at standard conditions.
		_rho_ws	 = self.rho(p=1333, T=293.15, RH=1333/self.svp(293.15), xc=xc) 	# Moist air component at standard conditions.
		
		_Ltotal  = ( _rho_a / _rho_axs ) * _La + ( _rho_w / _rho_ws ) * _Lw

		# print 'Dry air component = ', (_rho_a/_rho_axs)*(_n_a-1)*1e8
		# print 'Wet component = ', ( _rho_w / _rho_ws ) * (_n_w - 1)*1e8
		# print 'Sum of components = ', (_rho_a/_rho_axs)*(_n_a-1)*1e8+( _rho_w / _rho_ws ) * (_n_w - 1)*1e8
		return np.sqrt( (1 + 2 * _Ltotal) / (1 - _Ltotal) )


	def dn_tph(self, l, T, p, RH, xc, dl=None, dT=None, dP=None, dRH=None, dCO2=None):
		'''
			Calculates the uncertainty in the refractive index, following the Ciddor 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.
		
			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
		'''

		if dP == None:
			dP = self.dP
		if dT == None:
			dT = self.dT
		if dRH == None:
			dRH = self.dH
		if dCO2 == None:
			dCO2 = self.dCO2

		# constants
		_cf = 1.022
		_w0 = 295.235 	# um^-2
		_w1 = 2.6422	# um^-2
		_w2 = -0.032380	# um^-4
		_w3 = 0.004028	# um^-6

		_k0 = 238.0185 		# um^-2
		_k1 = 5792105  		# um^-2
		_k2 = 57.362		# um^-2
		_k3 = 167917		# um^-2	

		_A =  1.2378847e-5 # K^-2
		_B = -1.9121316e-2 # K^-2
		_C = 33.93711047
		_D = -6.3431645e3  # K^-2

		_Mw  = 0.018015   # kg/mol
		_R   = 8.314510	  # J mol^-1 K^-1, gas constant

		_a0 =  1.58123e-6 # K Pa^-1
		_a1 = -2.9331e-8  # Pa^-1
		_a2 =  1.1043e-10 # K^-1 Pa^-1
		_b0 =  5.707e-6   # K Pa^-1
		_b1 = -2.051e-8   # Pa^-1
		_c0 =  1.9898e-4  # K Pa^-1
		_c1 = -2.376e-6   # Pa^-1
		_d  =  1.83e-11   # K^2 Pa^-2
		_e  = -0.765e-8   # K^2 Pa^-2
		_t  = T - 273.15

		_aa = 1.00062
		_bb = 3.14e-8 	# Pa^-1
		_cc = 5.6e-7	# deg C^-2

		_RH_a = 0
		_RH_w = RH

		# Calculated values
		# _rho_obs = self.rho(p, T, RH, xc)
		# _ref	 = refraction(lat, h, p, _rho_obs)
		# _alpha   = _ref.rc * np.sin(np.deg2rad(z)) / (_ref.rc + _ref.H)
		_p_a 	 = p - RH*self.svp(T) 
		_p_w     = RH*self.svp(T)
		_x_w_a	 = self.xw(_p_a, T, _RH_a)
		_x_w_w 	 = self.f(_p_w, T)
		_Z_a 	 = self.Z(p=_p_a, T=T, RH=_RH_a)
		_Z_w 	 = self.Z(p=_p_w, T=T, RH=_RH_w, xw=_x_w_w)
		_Ma  	 = 1e-3 * (28.9635 + 12.011e-6 * (xc - 400) )	# Molar mass in kg/mol of dry air containing xc ppm of CO2
		_p_sv 	 = self.svp(T)
		_rho_a 	 = self.rho(p=p-RH*self.svp(T), T=T, RH=0, xc=xc)
		_rho_axs = self.rho(p=101325, T=288.15, RH=0, xc=xc)
		if RH == 0.0:
			# nan is returned for _rho_w if we put in a (partial) pressure of 0
			_rho_w = 0
		else:
			_rho_w = self.rho(p=RH*self.svp(T), T=T, RH=RH, xc=xc) 
		_rho_ws  = self.rho(p=1333, T=293.15, RH=1333/self.svp(293.15), xc=xc)
		_n_a1 = self.n_axs(l, xc)
		_n_w1 = self.n_ws(l)

		_L1_a = (_n_a1**2 - 1) / (_n_a1**2 + 2)  # Dry air component of L
		_L1_w = (_n_w1**2 - 1) / (_n_w1**2 + 2)
		_L1   = ( _rho_a / _rho_axs ) * _L1_a + ( _rho_w / _rho_ws ) * _L1_w # Overall L
		_n1   = np.sqrt( (1 + 2 * _L1) / (1 - _L1) )

		# Partial derivatives for derivation to l1
		_dn_wdl  = _cf * 1e-8 * (-1 * _w1 * 2/l**3 - _w2 * 4/l**5 - _w3 * 6/l**7)
		_dLdn_w1 = _rho_w / _rho_ws * 6 * _n_w1 / (_n_w1**2 + 2)**2
		_dn_adl  = 1e-8 * (1 + 0.534e-6 * (xc - 450)) * ( - (2 * _k1 * l ) / ((_k0 * l**2 - 1)**2) - (2 * _k3 * l) / ((_k2 * l**2 - 1)**2) )
		_dLdn_a1 = _rho_a / _rho_axs * 6 * _n_a1 / (_n_a1**2 + 2)**2
		_dndL1   = ( ( (1 + 2*_L1) / (1 - _L1)**2 ) + (2 / (1-_L1)) ) / (2 * np.sqrt( (1 + 2*_L1) / (1 - _L1) ) ) 
		_dndl   = _dndL1 * ( _dLdn_a1 * _dn_adl + _dLdn_w1 * _dn_wdl)

		# # Partial derivatives for derivation to l2


		# Partial derivatives for derivation to T
		# Dry air parts:
		_drho_adp_a  = _Ma / (_Z_a * _R * T) * (1 - _x_w_a * (1 - _Mw / _Ma))
		_dp_adp_sv 	 = - RH
		_dp_svdT 	 = (2*_A*T + _B - _D/(T**2)) * np.exp(_A * T**2 + _B * T + _C + _D/T )
		_drho_adTZ	 = - _p_a * _Ma / (_R * ( _Z_a * T)**2) * (1 - _x_w_a * (1 - _Mw / _Ma))
		_dTZdT_a	 = _p_a / T * (_a2 * _t**2 + _a1 * _t + _a0 + (_b1*_t + _b0) * _x_w_a + (_c0 + _c1*_t) * _x_w_a**2 ) - _p_a * (_a1 + 2*_a2*_t + _b1*_x_w_a + _c1*_x_w_a**2) - 2*(_p_a/T)**2 * (_d + _e * _x_w_a**2) + _Z_a
		_dTZdxw_a	 = 2*_e*_p_a**2 / T * _x_w_a - _p_a * (_b1*_t + _b0 + 2 * (_c0 + _c1*_t) * _x_w_a)
		_dxw_adf	 = _RH_a * _p_sv / _p_a
		_dfdT 		 = 2*_cc*_t
		_dxw_adp_sv	 = self.f(_p_a, T) * _RH_a / _p_a
		_drho_adxw_a = 0 # as rho_a has not water content, i.e. rho = p*M/ZRT


		# Moist air parts:
		_drho_wdp_w  = _Ma / (_Z_w * _R * T) * (1 - _x_w_w * (1 - _Mw / _Ma))
		_dp_wdp_sv 	 = RH
		_dp_svdT 	 = (2*_A*T + _B - _D/(T**2)) * np.exp(_A * T**2 + _B * T + _C + _D/T )
		_drho_wdTZ	 = - _p_w * _Ma / (_R * (_Z_w * T)**2) * (1 - _x_w_w * (1 - _Mw / _Ma))
		_dTZdT_w	 = _p_w / T * (_a2 * _t**2 + _a1 * _t + _a0 + (_b1*_t + _b0) * _x_w_w + (_c0 + _c1*_t) * _x_w_w**2 ) - _p_w * (_a1 + 2*_a2*_t + _b1*_x_w_w + _c1*_x_w_w**2) - 2*(_p_w/T)**2 * (_d + _e * _x_w_w**2) + _Z_w
		_dTZdxw_w	 = 2*_e*_p_w**2 / T * _x_w_w - _p_w * (_b1*_t + _b0 + 2 * (_c0 + _c1*_t) * _x_w_w)
		_dxw_wdf	 = 1
		_dfdT 		 = 2*_cc*_t
		_dxw_wdp_sv	 = self.f(_p_w, T) / _p_sv
		_drho_wdxw_w = _p_w * _Ma / (_Z_w * _R * T) * (_Mw / _Ma - 1)

		_tA =  _drho_adp_a * _dp_adp_sv * _dp_svdT +  _drho_adTZ * (_dTZdT_a + _dTZdxw_a * (_dxw_adf * _dfdT + _dxw_adp_sv * _dp_svdT) ) + _drho_adxw_a * ( _dxw_adf * _dfdT + _dxw_adp_sv * _dp_svdT )
		_tB =  _drho_wdp_w * _dp_wdp_sv * _dp_svdT +  _drho_wdTZ * (_dTZdT_w + _dTZdxw_w * (_dxw_wdf * _dfdT + _dxw_wdp_sv * _dp_svdT) ) + _drho_wdxw_w * ( _dxw_wdf * _dfdT + _dxw_wdp_sv * _dp_svdT )
		_dL1drho_a 	 = _L1_a / _rho_axs
		_dL1drho_w 	 = _L1_w / _rho_ws
		
		_dndT = _tA * (_dndL1 * _dL1drho_a) + _tB * (_dndL1 * _dL1drho_w)

		# Partial derivatives for derivation to p
		# Dry air parts:
		_drho_adZ 	= - _p_a * _Ma / ( _Z_a**2 * _R * T) * (1 - _x_w_a * (1 - _Mw / _Ma))
		_dZdp_a 	= - 1 / T * (_a0 + _a1 * _t + _a2 * _t**2 + (_b0 + _b1*_t)*_x_w_a + (_c0 + _c1*_t)*_x_w_a**2 ) + 2*_p_a/(T**2) * (_d + _e * _x_w_a**2)
		_dfdp 		= _bb
		_dxw_adp_a 	= - self.f(_p_a, T) * _RH_a * _p_sv / (_p_a)**2 

		# Moist air parts
		_dZdp_w 	= - 1 / T * (_a0 + _a1 * _t + _a2 * _t**2 + (_b0 + _b1*_t)*_x_w_w + (_c0 + _c1*_t)*_x_w_w**2 ) + 2*_p_w/(T**2) * (_d + _e * _x_w_w**2)
		# Not necessary, since dpwdp=0

		_pA = _drho_adp_a + _drho_adZ * _dZdp_a + _drho_adxw_a * (_dxw_adf * _dfdp + _dxw_adp_a)
		_dndP = _pA * (_dndL1 * _dL1drho_a)

		
		# Partial derivatives for derivation to H
		# Dry air parts:
		_dp_adH 	= - _p_sv
		_dZdxw_a 	= (_p_a / T)**2 * 2 * _e * _x_w_a - (_p_a / T) * ( (_b0 + _b1 * _t) + 2 * (_c0 + _c1 * _t) * _x_w_a )
		_dxw_adH 	= self.f(_p_a, T) * _p_sv / _p_a

		# Moist air parts:
		_dp_wdH 	= _p_sv
		_dZdxw_w 	= (_p_w / T)**2 * 2 * _e * _x_w_w - (_p_w / T) * ( (_b0 + _b1 * _t) + 2 * (_c0 + _c1 * _t) * _x_w_w )
		_dxw_wdH 	= 0
		_dxw_wdp_w 	= 0		
		_drho_wdZ 	= - _p_w * _Ma / ( _Z_w**2 * _R * T) * (1 - _x_w_w * (1 - _Mw / _Ma))

		_hA = _drho_adp_a * _dp_adH + _drho_adZ * ( _dZdp_a * _dp_adH + _dZdxw_a * (_dxw_adH + _dxw_adp_a * _dp_adH + _dxw_adf * _dfdp * _dp_adH) ) + _drho_adxw_a * (_dxw_adH + _dxw_adp_a * _dp_adH + _dxw_adf * _dfdp * _dp_adH)
		_hB = _drho_wdp_w * _dp_wdH + _drho_wdZ * ( _dZdp_w * _dp_wdH + _dZdxw_w * ( _bb * _p_sv ) ) + _drho_wdxw_w * ( _bb * _p_sv )
		_dndRH = _hA * (_dndL1 * _dL1drho_a) + _hB * (_dndL1 * _dL1drho_w)


		# Partial derivatives for derivation to xc
		# Dry air parts:
		_drho_adMa   = _p_a * (1 - _x_w_a) / (_Z_a * _R * T)
		_dL1drho_axs = - (_rho_a * _L1_a)/(_rho_axs**2)
		_drho_axsdMa = _drho_adMa
		_dL1dL_a	 = _rho_a / _rho_axs
		_dL_adn_a1	 = 6 * _n_a1 / (_n_a1**2 + 2)**2

		# Moist air parts:
		_drho_wdMa   = _p_w * (1 - _x_w_w) / (_Z_w * _R * T)
		_dL1drho_ws = - (_rho_w * _L1_w)/(_rho_ws**2)
		_drho_wsdMa = _drho_wdMa

		# Other parts:
		_dMadx_c   = 12.011e-9
		_dn_adx_c1 = (self.n_as(l) - 1) * 0.534e-6

		_dL1dx_c = _dL1drho_a * _drho_adMa * _dMadx_c + _dL1drho_w * _drho_wdMa * _dMadx_c + _dL1drho_axs * _drho_axsdMa * _dMadx_c + _dL1drho_ws * _drho_wsdMa * _dMadx_c + _dL1dL_a * _dL_adn_a1 * _dn_adx_c1
		_dndCO2 = _dndL1 * _dL1dx_c

		# Final error:
		_nError = np.sqrt( (_dndl * dl)**2 + (_dndT * dT)**2  + (_dndP * dP)**2 + (_dndRH * dRH)**2 + (_dndCO2 * dCO2)**2 )

		return _nError

###
### THIS ERROR CLASS IS OUTDATED
### Use Observatory.dn_tph(), dispersion.cassiniError() or refraction.cassiniError() instead.
###
class error(Observatory):
	def __init__(self, dT=0.2, dP=20, dH=0.02, dCO2=20, dl=0.001):
		Observatory.__init__(self)
		self.dT 	= dT 	# in K
		self.dP 	= dP 	# in Pa
		self.dH 	= dH	# relative humidity between 0 and 1
		self.dCO2 	= dCO2	# in ppm
		self.dl 	= dl 	# in micron
		print('THIS ERROR CLASS IS OUTDATED\nUse Observatory.dn_tph(), dispersion.cassiniError() or refraction.cassiniError() instead.')

	def dn_as(self, l, dl=None):
		'''
			Formula that calculate the error on the refractive index for standard 
			atmospheric conditions, i.e. 15 C, 101,325 Pa, 0% RH, 450 ppm CO2.
			It returns dn_as, not d(n_as - 1) * 1e8!
		'''
		if not dl:
			dl = self.dl

		_k0 = 238.0185 		# um^-2
		_k1 = 5792105  		# um^-2
		_k2 = 57.362		# um^-2
		_k3 = 167917		# um^-2	
		_dn_as = np.sqrt((( - (2 * _k1 * l ) / ((_k0 * l**2 - 1)**2) - (2 * _k3 * l) / ((_k2 * l**2 - 1)**2) )*dl)**2)
		return _dn_as / 1e8


	def dn_axs(self, l, xc, dCO2=None, dl=None):
		'''
			Formula that calculates the error on the refractive index for standard 
			atmospheric	conditions, i.e. 15 C, 101,325 Pa, 0% RH, but with xc ppm of CO2.
			It returns dn_axs, not d(n_axs - 1)
		'''
		if not dCO2:
			dCO2 = self.dCO2
		if not dl:
			dl = self.dl
		_dn_as = self.dn_as(l, dl)
		_n_as  = self.n_as(l)
		_dn_axs = np.sqrt( ( (_n_as - 1) * 0.534e-6 * dCO2 )**2 + ( (0.534e-6 * (xc - 450) + 1) * _dn_as )**2 )
		return _dn_axs


	def dn_ws(self, l, dl=None):
		'''
			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
		'''
		if not dl:
			dl = self.dl
		_cf = 1.022
		_w0 = 295.235 	# um^-2
		_w1 = 2.6422	# um^-2
		_w2 = -0.032380	# um^-4
		_w3 = 0.004028	# um^-6
		
		_dndl = _cf * (-1 * _w1 * 2/l**3 - _w2 * 4/l**5 - _w3 * 6/l**7)
		return np.sqrt( (_dndl * dl)**2 ) / 1e8


	def df(self, p, T, dP=None, dT=None):
		'''
			Uncertainty in enhancement factor of
			water vapour in air.
			t in degrees celsius
		'''
		if not dP:
			dP = self.dP
		if not dT:
			dT = self.dT

		_t = T - 273.15
		_b = 3.14e-8 	# Pa^-1
		_c = 5.6e-7		# deg C^-2
		return np.sqrt( (_b*dP)**2 + (2*_c*_t*dT)**2 )


	def dsvp(self, T, dT=None):
		'''
			Uncertainty in saturation vapour pressure in Pa.
		'''
		if not dT:
			dT = self.dT
		_A =  1.2378847e-5 # K^-2
		_B = -1.9121316e-2 # K^-2
		_C = 33.93711047
		_D = -6.3431645e3  # K^-2
		_p1 = 2*_A*T + _B - _D/(T**2)
		_p2 = np.exp(_A*T**2 + _B*T + _C + _D/T)
		return np.sqrt( (_p1 * _p2 * dT)**2 )


	def dxw(self, p, T, RH, dP=None, dT=None, dRH=None):
		'''
			Uncertainty in molar fraction of water vapour in moist air.
		'''
		if not dP:
			dP = self.dP
		if not dT:
			dT = self.dT
		if not dRH:
			dRH = self.dH

		_f 		= self.f(p, T)
		_svp 	= self.svp(T)

		_df 	= self.df(p=p, T=T, dP=dP, dT=dT)
		_dsvp	= self.dsvp(T=T, dT=dT)
		_p1		= (RH * _svp / p) * _df
		_p2		= (_f * _svp / p) * dRH
		_p3		= (_f * RH  / p) * _dsvp
		_p4		= (-1 * _f * RH * _svp / p**2) * dP
		return np.sqrt( _p1**2 + _p2**2 + _p3**2 + _p4**2)


	def dZ(self, p, T, RH, xw='', dP=None, dT=None, dRH=None):
		'''
			Uncertainty in compressibility of moist air
		'''
		if not dP:
			dP = self.dP
		if not dT:
			dT = self.dT
		if not dRH:
			dRH = self.dH
		_dt  = dT
		_dxw = self.dxw(p=p, T=T, RH=RH, dP=dP, dT=dT, dRH=dRH)

		_a0 =  1.58123e-6 # K Pa^-1
		_a1 = -2.9331e-8  # Pa^-1
		_a2 =  1.1043e-10 # K^-1 Pa^-1
		_b0 =  5.707e-6   # K Pa^-1
		_b1 = -2.051e-8   # Pa^-1
		_c0 =  1.9898e-4  # K Pa^-1
		_c1 = -2.376e-6   # Pa^-1
		_d  =  1.83e-11   # K^2 Pa^-2
		_e  = -0.765e-8   # K^2 Pa^-2
		_t  = T - 273.15
		if xw == '':
			_xw = self.xw(p, T, RH)
		else:
			_xw = xw

		_dZdp  = -1/T * ( _a0 + _a1*_t + _a2*_t**2 + (_b0 + _b1*_t)*_xw + (_c0 + _c1*_t)*_xw**2) + (2*p/T**2) * (_d + _e*_xw**2)
		_dZdT  = p/T**2 * ( _a0 + _a1*_t + _a2*_t**2 + (_b0 + _b1*_t)*_xw + (_c0 + _c1*_t)*_xw**2) - 2*p**2/T**3 * (_d + _e*_xw**2)
		_dZdt  = -1*p/T * (_a1 + 2*_a2*_t + _b1*_xw + _c1*_xw**2)
		_dZdxw = -1*p/T * ( (_b0 + _b1*_t) + 2*(_c0 + _c1*_t)*_xw ) + 2 * (p/T)**2 * _e * _xw
		_dZ    = np.sqrt( (_dZdp * dP)**2 + (_dZdT * dT)**2 + (_dZdt * _dt)**2 + (_dZdxw * _dxw)**2 )
		return _dZ


	def drho(self, p, T, RH, xc, xw='', dP=None, dT=None, dRH=None, dCO2=None):
		'''
			Calculates the uncertainty in the density of air for the given parameters.
			T in Kelvin
			p in Pascal
			RH between 0 and 1
			xc in ppm
		'''
		if not dP:
			dP = self.dP
		if not dT:
			dT = self.dT
		if not dRH:
			dRH = self.dH
		if not dCO2:
			dCO2 = self.dCO2

		_Mw  = 0.018015 										# kg/mol
		_R   = 8.314510										# J mol^-1 K^-1, gas constant
		_Ma  = 1e-3 * (28.9635 + 12.011e-6 * (xc - 400) )	# Molar mass in kg/mol of dry air containing xc ppm of CO2
		if xw == '':
			_xw = self.xw(p, T, RH)
		else: 
			_xw = xw
		_Z   = self.Z(p=p, T=T, RH=RH, xw=_xw)

		_dMa = np.sqrt( (12.011e-9*dCO2)**2 )
		_dZ  = self.dZ(p=p, T=T, RH=RH, xw=_xw, dP=dP, dT=dT, dRH=dRH)
		_dxw = self.dxw(p=p, T=T, RH=RH, dP=dP, dT=dT, dRH=dRH)

		_drhodp  = _Ma / (_Z * _R * T) * (1 - _xw * (1 - _Mw / _Ma))
		_drhodMa =   p / (_Z * _R * T) * (1 - _xw)
		_drhodZ  = -1 * (p * _Ma) / (_Z**2 * _R * T) * (1 - _xw * (1 - _Mw / _Ma))
		_drhodT  = -1 * (p * _Ma) / (_Z * _R * T**2) * (1 - _xw * (1 - _Mw / _Ma))
		_drhodxw = p / (_Z * _R * T) * (_Mw - _Ma)

		return np.sqrt( (_drhodp * dP)**2 + (_drhodMa * _dMa)**2 + (_drhodZ * _dZ)**2 + (_drhodT * dT)**2 + (_drhodxw * _dxw)**2 )


	def dn_tph(self, l, T, p, RH, xc, dl=None, dT=None, dP=None, dRH=None, dCO2=None):
		'''
			[DEPRECATED]
			Calculates the uncertainty in the refraction of the atmosphere for any conditions.
		'''
		if not dP:
			dP = self.dP
		if not dT:
			dT = self.dT
		if not dRH:
			dRH = self.dH
		if not dCO2:
			dCO2 = self.dCO2

		_n_a 	 = self.n_axs(l, xc)
		_n_w 	 = self.n_ws(l)
		_dn_a    = self.dn_axs(l=l, xc=xc, dl=dl, dCO2=dCO2)
		_dn_w    = self.dn_ws(l=l, dl=dl)

		_La 	 = (_n_a**2 - 1) / (_n_a**2 + 2)
		_Lw 	 = (_n_w**2 - 1) / (_n_w**2 + 2)
		_dLa	 = np.sqrt( ( (6*_n_a) / (_n_a**2+2)**2 * _dn_a)**2 )
		_dLw	 = np.sqrt( ( (6*_n_w) / (_n_w**2+2)**2 * _dn_w)**2 )
		
		# Calculate the density of the dry air components, using the partial pressure of the dry air, i.e. p_dry = p_atm - p_water
		_pw 	 = RH*self.svp(T)
		_pa 	 = p - _pw
		_rho_a 	 = self.rho(p=_pa, T=T, RH=0, xc=xc)				# Dry air component for the actual conditions.
		_rho_axs = self.rho(p=101325, T=288.15, RH=0, xc=xc) 		# Dry air component at standard conditions.

		_dpa 	 = np.sqrt( (dP)**2 + (-1*self.svp(T) * dRH)**2 + (-1 * RH * self.dsvp(T=T, dT=dT) )**2 )
		_dpw 	 = np.sqrt( (RH * self.dsvp(T=T, dT=dT))**2 + (self.svp(T) * dRH)**2 )
		_drho_a  = self.drho(p=_pa, T=T, RH=RH, xc=xc, dP=_dpa, dT=dT, dRH=dRH, dCO2=dCO2)
		_drho_axs= 0
		
		if RH == 0:
			# nan is returned for _rho_w if we put in a (partial) pressure of 0
			_rho_w = 0
			_drho_w = 0
		else:
			# Partial pressure of water vapour is p_water = h*svp
			_rho_w 	 = self.rho(p=_pw, T=T, RH=RH, xc=xc)			# Moist air component for the actual conditions.
			_drho_w  = self.drho(p=_pw, T=T, RH=RH, xc=xc, dP=_dpw, dT=dT, dRH=dRH, dCO2=dCO2)
		_rho_ws	 = self.rho(p=1333, T=293.15, RH=1333/self.svp(293.15), xc=xc) 	# Moist air component at standard conditions.
		_drho_ws = 0

		
		_dLdL_a 	= _rho_a / _rho_axs 
		_dLdL_w 	= _rho_w / _rho_ws
		_dLdrho_a 	= _La / _rho_axs
		_dLdrho_axs = - (_rho_a * _La) / (_rho_axs**2)
		_dLdrho_w 	= _Lw / _rho_ws
		_dLdrho_ws 	= - (_rho_w * _Lw) / _rho_ws**2

		_Ltotal  = ( _rho_a / _rho_axs ) * _La + ( _rho_w / _rho_ws ) * _Lw
		_dLtotal = np.sqrt( (_dLdL_a * _dLa)**2 + (_dLdL_w * _dLw)**2 + (_dLdrho_a * _drho_a)**2 
			+ (_dLdrho_axs * _drho_axs)**2 + (_dLdrho_w * _drho_w)**2 + (_dLdrho_ws * _drho_ws)**2 )

		_nLL = np.sqrt( (1 + 2 * _Ltotal) / (1 - _Ltotal) )
		_dnLL = np.sqrt( ( (3*np.sqrt(1 - _Ltotal )) / ( 2 * (1 - _Ltotal)**2 * np.sqrt(1 + 2*_Ltotal) ) * _dLtotal )**2 )
		
		return (_nLL, _dnLL)
###
### END OF OUTDATED ERROR CLASS
###
	


class dispersion(Observatory):
	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 planeParallel(self, n1, n2, zenith):
		'''
			Refraction based on a flat atmosphere, easily derived from Snell's law.
		'''
		_R1 = np.arcsin(n1 * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
		_R2 = np.arcsin(n2 * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
		return np.degrees(_R1) - np.degrees(_R2)



	def cassini(self, n1, n2, zenith, reduced_height=None):
		'''
			Refraction of spherical atmosphere, derived from geometry using Sine law
		'''
		if reduced_height:
			H = reduced_height
		else:
			H = self.H
		
		_R1 = np.arcsin(n1 * self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H)) - np.arcsin(self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H) )
		_R2 = np.arcsin(n2 * self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H)) - np.arcsin(self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H) )	
		return np.degrees(_R1) - np.degrees(_R2)



	def tan5(self, n1, n2, zenith, reduced_height=None):
		'''
			Same as oriani's theorem, but including a higher order term.
			Found in Corbard et al., 2019.
		'''
		if reduced_height:
			H = reduced_height
		else:
			H = self.H
		
		_rList = []
		for _n in [n1, n2]:
			_a = _n - 1
			_b = H / self.rc
			_R = _a * (1 - _b) * np.tan( np.deg2rad(zenith) )  - _a * (_b - _a/2) * np.tan( np.deg2rad(zenith) )**3 + 3 * _a * (_b - _a/2)**2 * np.tan( np.deg2rad(zenith) )**5
			_rList.append(np.degrees(_R))
		return _rList[0] - _rList[1]



	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.
		'''
		if reduced_height:
			H = reduced_height
		else:
			H = self.H

		_rList = []
		for _n in [n1, n2]:
			_a = _n - 1
			_b = H / self.rc
			_R = _a * (1 - _b) * np.tan( np.deg2rad(zenith) ) - _a * (_b - _a/2) * np.tan( np.deg2rad(zenith) )**3
			_rList.append(np.degrees(_R))
		return _rList[0] - _rList[1]



	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'.
		'''
		if reduced_height:
			H = reduced_height
		else:
			H = self.H

		_rList = []
		for _n in [n1, n2]:
			_a = _n - 1
			_b = H / self.rc
			_R = _a * ( (2 - _a) / (np.sqrt( 2*_b - _a )) ) * np.sin( np.deg2rad(zenith) ) * self.psi( (np.cos(np.deg2rad(zenith))) / np.sqrt(2*_b - _a) )
			_rList.append(np.degrees(_R))
		return _rList[0] - _rList[1]



	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 
			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.
		'''
		if not R0:
			R0 = self.rc
		if useStandardAtmosphere:
			_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('./US1976 Standard Atmosphere/si1py_modified.txt', unpack=True, skip_header=6)
		else:
			_Dens = rhoData
			_h 	  = heightData 

		
		_rList 					= []
		z0 						= np.deg2rad(z0)
		
		# Integration of the refraction integral over the atmosphere layers.
		for _n in [n1, n2]:
			_Dens 	= np.asarray(_Dens)
			_h 	 	= np.asarray(_h) 
			_R 		= R0 + _h * 1000
			
			# Gladstone-Dale relation between refractive index and density
			_n_h = 1 + (_n - 1) / _Dens[0] * _Dens 

			# Calculation of the zenith angle for given height and refractive index, from the refractive invariant
			_z 		= np.arcsin( (_n * R0) / (_n_h * _R) * np.sin(z0) )
			
			# Calculation of log space derivative of n and r
			_lnn = np.log(_n_h)
			_lnr = np.log(_R)
			_dlnn = np.asarray([_lnn[i+1] - _lnn[i] for i in range(len(_lnn)-1)])
			_dlnr = np.asarray([_lnr[i+1] - _lnr[i] for i in range(len(_lnr)-1)])
			_dz   = np.asarray([(_z[i] - _z[i+1]) for i in range(len(_z)-1)])

			# Calculation of the integrand in eq. 3 in Auer & Standish
			_Integrand = -1*(_dlnn/_dlnr) / (1 + _dlnn/_dlnr) * _dz
			_rList.append(np.degrees(np.sum(_Integrand)))
		return _rList[0] - _rList[1]



	def psi(self, x):
		_f   = lambda a: np.exp(-a**2)
		_int = integrate.quad(_f, x, np.inf)
		return np.exp(x**2) * _int[0]


	def set_rc(self):
		'''
			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.
		'''
		_A   = 6378.137 # km
		_B   = 6356.752 # km
		_phi = np.deg2rad(self.lat)
		
		_rc0 = (_A*_B)**2 / (_A**2 * np.cos(_phi)**2 + _B**2 * np.sin(_phi)**2)**(3/2)
		return _rc0 + self.h/1000


	def set_H(self, P, rho):
		'''
			Calculates the reduced height of the atmosphere, assuming ideal gas law.
		'''
		_phi = np.deg2rad(self.lat)
		_c1 = 5.2790414e-3
		_c2 = 2.32718e-5
		_c3 = 1.262e-7
		_c4 = 7e-10
		_g0 = 9.780327 # m/s^2
		_g0_local = _g0 * (1 + _c1 * np.sin(_phi)**2 + _c2 * np.sin(_phi)**4 + _c3 * np.sin(_phi)**6 + _c4 * np.sin(_phi)**8)
		_g = _g0_local - (3.0877e-6 - 4.3e-9 * np.sin(_phi)**2) * self.h + 7.2e-13 * self.h**2
		return P / (rho * _g)

	def H_isotherm(self, T):
		'''
			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(self.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 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.

			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.

			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
		'''
		if lat == None:
			lat = self.lat
		else:
			self.lat = lat
		if h == None:
			h 	= self.h
		else:
			self.h = h

		# constants
		_cf = 1.022
		_w0 = 295.235 	# um^-2
		_w1 = 2.6422	# um^-2
		_w2 = -0.032380	# um^-4
		_w3 = 0.004028	# um^-6

		_k0 = 238.0185 		# um^-2
		_k1 = 5792105  		# um^-2
		_k2 = 57.362		# um^-2
		_k3 = 167917		# um^-2	

		_A =  1.2378847e-5 # K^-2
		_B = -1.9121316e-2 # K^-2
		_C = 33.93711047
		_D = -6.3431645e3  # K^-2

		_Mw  = 0.018015   # kg/mol
		_R   = 8.314510	  # J mol^-1 K^-1, gas constant

		_a0 =  1.58123e-6 # K Pa^-1
		_a1 = -2.9331e-8  # Pa^-1
		_a2 =  1.1043e-10 # K^-1 Pa^-1
		_b0 =  5.707e-6   # K Pa^-1
		_b1 = -2.051e-8   # Pa^-1
		_c0 =  1.9898e-4  # K Pa^-1
		_c1 = -2.376e-6   # Pa^-1
		_d  =  1.83e-11   # K^2 Pa^-2
		_e  = -0.765e-8   # K^2 Pa^-2
		_t  = T - 273.15

		_aa = 1.00062
		_bb = 3.14e-8 	# Pa^-1
		_cc = 5.6e-7	# deg C^-2

		_RH_a = 0
		_RH_w = RH

		# Calculated values
		_rho_obs = self.rho(p, T, RH, xc)
		self.rc  = self.set_rc() * 1000
		self.H   = self.set_H(p, _rho_obs) 
		_alpha   = self.rc * np.sin(np.deg2rad(z)) / (self.rc + self.H)
		_p_a 	 = p - RH*self.svp(T) 
		_p_w     = RH*self.svp(T)
		_x_w_a	 = self.xw(_p_a, T, _RH_a)
		_x_w_w 	 = self.f(_p_w, T)
		_Z_a 	 = self.Z(p=_p_a, T=T, RH=_RH_a)
		_Z_w 	 = self.Z(p=_p_w, T=T, RH=_RH_w, xw=_x_w_w)
		_Ma  	 = 1e-3 * (28.9635 + 12.011e-6 * (xc - 400) )	# Molar mass in kg/mol of dry air containing xc ppm of CO2
		_p_sv 	 = self.svp(T)
		_rho_a 	 = self.rho(p=p-RH*self.svp(T), T=T, RH=0, xc=xc)
		_rho_axs = self.rho(p=101325, T=288.15, RH=0, xc=xc)
		if RH == 0.0:
			# nan is returned for _rho_w if we put in a (partial) pressure of 0
			_rho_w = 0
		else:
			_rho_w = self.rho(p=RH*self.svp(T), T=T, RH=RH, xc=xc) 
		_rho_ws  = self.rho(p=1333, T=293.15, RH=1333/self.svp(293.15), xc=xc)
		_n_a1 	 = self.n_axs(l1, xc)
		_n_w1 	 = self.n_ws(l1)
		_n_a2 	 = self.n_axs(l2, xc)
		_n_w2 	 = self.n_ws(l2)

		_L1_a 	 = (_n_a1**2 - 1) / (_n_a1**2 + 2)  # Dry air component of L
		_L1_w 	 = (_n_w1**2 - 1) / (_n_w1**2 + 2)
		_L1   	 = ( _rho_a / _rho_axs ) * _L1_a + ( _rho_w / _rho_ws ) * _L1_w # Overall L
		_L2_a 	 = (_n_a2**2 - 1) / (_n_a2**2 + 2) # Dry air component of L
		_L2_w 	 = (_n_w2**2 - 1) / (_n_w2**2 + 2)
		_L2    	 = ( _rho_a / _rho_axs ) * _L2_a + ( _rho_w / _rho_ws ) * _L2_w # Overall L

		_n1   	 = np.sqrt( (1 + 2 * _L1) / (1 - _L1) )
		_n2   	 = np.sqrt( (1 + 2 * _L2) / (1 - _L2) )

		# Partial derivatives for derivation to l1
		_dn_wdl1 = _cf * 1e-8 * (-1 * _w1 * 2/l1**3 - _w2 * 4/l1**5 - _w3 * 6/l1**7)
		_dLdn_w1 = _rho_w / _rho_ws * 6 * _n_w1 / (_n_w1**2 + 2)**2
		_dn_adl1 = 1e-8 * (1 + 0.534e-6 * (xc - 450)) * ( - (2 * _k1 * l1 ) / ((_k0 * l1**2 - 1)**2) - (2 * _k3 * l1) / ((_k2 * l1**2 - 1)**2) )
		_dLdn_a1 = _rho_a / _rho_axs * 6 * _n_a1 / (_n_a1**2 + 2)**2
		_dndL1   = ( ( (1 + 2*_L1) / (1 - _L1)**2 ) + (2 / (1-_L1)) ) / (2 * np.sqrt( (1 + 2*_L1) / (1 - _L1) ) ) 
		_dddn1   = _alpha / np.sqrt( 1 - _alpha**2 * _n1**2 )
		_dddl1   = _dddn1 * _dndL1 * ( _dLdn_a1 * _dn_adl1 + _dLdn_w1 * _dn_wdl1)

		# Partial derivatives for derivation to l2
		_dn_wdl2 = _cf * 1e-8 * (-1 * _w1 * 2/l2**3 - _w2 * 4/l2**5 - _w3 * 6/l2**7)
		_dLdn_w2 = _rho_w / _rho_ws * 6 * _n_w2 / (_n_w2**2 + 2)**2
		_dn_adl2 = 1e-8 * (1 + 0.534e-6 * (xc - 450)) * ( - (2 * _k1 * l2 ) / ((_k0 * l2**2 - 1)**2) - (2 * _k3 * l2) / ((_k2 * l2**2 - 1)**2) )
		_dLdn_a2 = _rho_a / _rho_axs * 6 * _n_a2 / (_n_a2**2 + 2)**2
		_dndL2   = ( ( (1 + 2*_L2) / (1 - _L2)**2 ) + (2 / (1-_L2)) ) / (2 * np.sqrt( (1 + 2*_L2) / (1 - _L2) ) ) 
		_dddn2   = _alpha / np.sqrt( 1 - _alpha**2 * _n2**2 )
		_dddl2   = _dddn2 * _dndL2 * ( _dLdn_a2 * _dn_adl2 + _dLdn_w2 * _dn_wdl2)


		# Partial derivatives for derivation to T
		# Dry air parts:
		_drho_adp_a  = _Ma / (_Z_a * _R * T) * (1 - _x_w_a * (1 - _Mw / _Ma))
		_dp_adp_sv 	 = - RH
		_dp_svdT 	 = (2*_A*T + _B - _D/(T**2)) * np.exp(_A * T**2 + _B * T + _C + _D/T )
		_drho_adTZ	 = - _p_a * _Ma / (_R * ( _Z_a * T)**2) * (1 - _x_w_a * (1 - _Mw / _Ma))
		_dTZdT_a	 = _p_a / T * (_a2 * _t**2 + _a1 * _t + _a0 + (_b1*_t + _b0) * _x_w_a + (_c0 + _c1*_t) * _x_w_a**2 ) - _p_a * (_a1 + 2*_a2*_t + _b1*_x_w_a + _c1*_x_w_a**2) - 2*(_p_a/T)**2 * (_d + _e * _x_w_a**2) + _Z_a
		_dTZdxw_a	 = 2*_e*_p_a**2 / T * _x_w_a - _p_a * (_b1*_t + _b0 + 2 * (_c0 + _c1*_t) * _x_w_a)
		_dxw_adf	 = _RH_a * _p_sv / _p_a
		_dfdT 		 = 2*_cc*_t
		_dxw_adp_sv	 = self.f(_p_a, T) * _RH_a / _p_a
		_drho_adxw_a = 0 # as rho_a has not water content, i.e. rho = p*M/ZRT


		# Moist air parts:
		_drho_wdp_w  = _Ma / (_Z_w * _R * T) * (1 - _x_w_w * (1 - _Mw / _Ma))
		_dp_wdp_sv 	 = RH
		_dp_svdT 	 = (2*_A*T + _B - _D/(T**2)) * np.exp(_A * T**2 + _B * T + _C + _D/T )
		_drho_wdTZ	 = - _p_w * _Ma / (_R * (_Z_w * T)**2) * (1 - _x_w_w * (1 - _Mw / _Ma))
		_dTZdT_w	 = _p_w / T * (_a2 * _t**2 + _a1 * _t + _a0 + (_b1*_t + _b0) * _x_w_w + (_c0 + _c1*_t) * _x_w_w**2 ) - _p_w * (_a1 + 2*_a2*_t + _b1*_x_w_w + _c1*_x_w_w**2) - 2*(_p_w/T)**2 * (_d + _e * _x_w_w**2) + _Z_w
		_dTZdxw_w	 = 2*_e*_p_w**2 / T * _x_w_w - _p_w * (_b1*_t + _b0 + 2 * (_c0 + _c1*_t) * _x_w_w)
		_dxw_wdf	 = 1 # _RH_w * _p_sv / _p_w
		_dfdT 		 = 2*_cc*_t
		_dxw_wdp_sv	 = self.f(_p_w, T) / _p_sv # * _RH_w / _p_w
		_drho_wdxw_w = _p_w * _Ma / (_Z_w * _R * T) * (_Mw / _Ma - 1)

		_tA =  _drho_adp_a * _dp_adp_sv * _dp_svdT +  _drho_adTZ * (_dTZdT_a + _dTZdxw_a * (_dxw_adf * _dfdT + _dxw_adp_sv * _dp_svdT) ) + _drho_adxw_a * ( _dxw_adf * _dfdT + _dxw_adp_sv * _dp_svdT )
		_tB =  _drho_wdp_w * _dp_wdp_sv * _dp_svdT +  _drho_wdTZ * (_dTZdT_w + _dTZdxw_w * (_dxw_wdf * _dfdT + _dxw_wdp_sv * _dp_svdT) ) + _drho_wdxw_w * ( _dxw_wdf * _dfdT + _dxw_wdp_sv * _dp_svdT )
		_dL1drho_a 	 = _L1_a / _rho_axs
		_dL2drho_a 	 = _L2_a / _rho_axs
		_dL1drho_w 	 = _L1_w / _rho_ws
		_dL2drho_w 	 = _L2_w / _rho_ws
		
		_dddT = _tA * (_dddn1 * _dndL1 * _dL1drho_a - _dddn2 * _dndL2 * _dL2drho_a) + _tB * (_dddn1 * _dndL1 * _dL1drho_w - _dddn2 * _dndL2 * _dL2drho_w)

		# Partial derivatives for derivation to p
		# Dry air parts:
		_drho_adZ 	= - _p_a * _Ma / ( _Z_a**2 * _R * T) * (1 - _x_w_a * (1 - _Mw / _Ma))
		_dZdp_a 	= - 1 / T * (_a0 + _a1 * _t + _a2 * _t**2 + (_b0 + _b1*_t)*_x_w_a + (_c0 + _c1*_t)*_x_w_a**2 ) + 2*_p_a/(T**2) * (_d + _e * _x_w_a**2)
		_dfdp 		= _bb
		_dxw_adp_a 	= - self.f(_p_a, T) * _RH_a * _p_sv / (_p_a)**2 

		# Moist air parts
		_dZdp_w 	= - 1 / T * (_a0 + _a1 * _t + _a2 * _t**2 + (_b0 + _b1*_t)*_x_w_w + (_c0 + _c1*_t)*_x_w_w**2 ) + 2*_p_w/(T**2) * (_d + _e * _x_w_w**2)
		# Not necessary, since dpwdp=0

		_pA = _drho_adp_a + _drho_adZ * _dZdp_a + _drho_adxw_a * (_dxw_adf * _dfdp + _dxw_adp_a)
		_dddP = _pA * (_dddn1 * _dndL1 * _dL1drho_a - _dddn2 * _dndL2 * _dL2drho_a)

		
		# Partial derivatives for derivation to H
		# Dry air parts:
		_dp_adH 	= - _p_sv
		_dZdxw_a 	= (_p_a / T)**2 * 2 * _e * _x_w_a - (_p_a / T) * ( (_b0 + _b1 * _t) + 2 * (_c0 + _c1 * _t) * _x_w_a )
		_dxw_adH 	= self.f(_p_a, T) * _p_sv / _p_a

		# Moist air parts:
		_dp_wdH 	= _p_sv
		_dZdxw_w 	= (_p_w / T)**2 * 2 * _e * _x_w_w - (_p_w / T) * ( (_b0 + _b1 * _t) + 2 * (_c0 + _c1 * _t) * _x_w_w )
		_dxw_wdH 	= 0
		_dxw_wdp_w 	= 0		
		_drho_wdZ 	= - _p_w * _Ma / ( _Z_w**2 * _R * T) * (1 - _x_w_w * (1 - _Mw / _Ma))

		_hA = _drho_adp_a * _dp_adH + _drho_adZ * ( _dZdp_a * _dp_adH + _dZdxw_a * (_dxw_adH + _dxw_adp_a * _dp_adH + _dxw_adf * _dfdp * _dp_adH) ) + _drho_adxw_a * (_dxw_adH + _dxw_adp_a * _dp_adH + _dxw_adf * _dfdp * _dp_adH)
		_hB = _drho_wdp_w * _dp_wdH + _drho_wdZ * ( _dZdp_w * _dp_wdH + _dZdxw_w * ( _bb * _p_sv ) ) + _drho_wdxw_w * ( _bb * _p_sv )
		_dddRH = _hA * (_dddn1 * _dndL1 * _dL1drho_a - _dddn2 * _dndL2 * _dL2drho_a) + _hB * (_dddn1 * _dndL1 * _dL1drho_w - _dddn2 * _dndL2 * _dL2drho_w)


		# Partial derivatives for derivation to xc
		# Dry air parts:
		_drho_adMa   = _p_a * (1 - _x_w_a) / (_Z_a * _R * T)
		_dL1drho_axs = - (_rho_a * _L1_a)/(_rho_axs**2)
		_dL2drho_axs = - (_rho_a * _L2_a)/(_rho_axs**2)
		_drho_axsdMa = _drho_adMa
		_dL1dL_a	 = _rho_a / _rho_axs
		_dL_adn_a1	 = 6 * _n_a1 / (_n_a1**2 + 2)**2
		_dL2dL_a	 = _dL1dL_a
		_dL_adn_a2	 = 6 * _n_a2 / (_n_a2**2 + 2)**2

		# Moist air parts:
		_drho_wdMa   = _p_w * (1 - _x_w_w) / (_Z_w * _R * T)
		_dL1drho_ws = - (_rho_w * _L1_w)/(_rho_ws**2)
		_dL2drho_ws = - (_rho_w * _L2_w)/(_rho_ws**2)
		_drho_wsdMa = _drho_wdMa

		# Other parts:
		_dMadx_c   = 12.011e-9
		_dn_adx_c1 = (self.n_as(l1) - 1) * 0.534e-6
		_dn_adx_c2 = (self.n_as(l2) - 1) * 0.534e-6

		_dL1dx_c = _dL1drho_a * _drho_adMa * _dMadx_c + _dL1drho_w * _drho_wdMa * _dMadx_c + _dL1drho_axs * _drho_axsdMa * _dMadx_c + _dL1drho_ws * _drho_wsdMa * _dMadx_c + _dL1dL_a * _dL_adn_a1 * _dn_adx_c1
		_dL2dx_c = _dL2drho_a * _drho_adMa * _dMadx_c + _dL2drho_w * _drho_wdMa * _dMadx_c + _dL2drho_axs * _drho_axsdMa * _dMadx_c + _dL2drho_ws * _drho_wsdMa * _dMadx_c + _dL2dL_a * _dL_adn_a2 * _dn_adx_c2
		_dddCO2 = _dddn1 * _dndL1 * _dL1dx_c - _dddn2 * _dndL2 * _dL2dx_c


		# Partial derivatives for derivation to z
		_zenith = np.deg2rad(z)
		_dz   	= np.deg2rad(dz)
		_R 		= self.rc
		_H 		= self.H
		_dddz 	= (_R * np.cos(_zenith)) / (_R + _H) * ( _n1 / np.sqrt(1-(_n1**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2) - _n2 / np.sqrt(1-(_n2**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2) )		

		# Final error:
		_dispersionError = np.sqrt( (_dddl1 * dl1)**2 + (_dddl2 * dl2)**2 + (_dddT * dT)**2  + (_dddP * dP)**2 + (_dddRH * dRH)**2 + (_dddCO2 * dCO2)**2  + (_dddz * _dz)**2)

		return np.degrees(_dispersionError)





class refraction(Observatory):
	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
		# 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 planeParallel(self, n, zenith):
		'''
			Refraction based on a flat atmosphere, easily derived from Snell's law.
		'''
		_R = np.arcsin(n * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
		return np.degrees(_R)



	def cassini(self, n1, zenith, reduced_height=None):
		'''
			Refraction of spherical atmosphere, derived from geometry using Sine law
		'''
		if reduced_height:
			H = reduced_height
		else:
			H = self.H
		
		_R1 = np.arcsin(n1 * self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H)) - np.arcsin(self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H) )
		return np.degrees(_R1)



	def tan5(self, n1, zenith, reduced_height=None):
		'''
			Same as oriani's theorem, but including a higher order term.
			Found in Corbard et al., 2019.
		'''
		if reduced_height:
			H = reduced_height
		else:
			H = self.H
		
		_a = n1 - 1
		_b = H / self.rc
		_R = _a * (1 - _b) * np.tan( np.deg2rad(zenith) )  - _a * (_b - _a/2) * np.tan( np.deg2rad(zenith) )**3 + 3 * _a * (_b - _a/2)**2 * np.tan( np.deg2rad(zenith) )**5
		return np.degrees(_R)



	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.
		'''
		if reduced_height:
			H = reduced_height
		else:
			H = self.H

		_a = n1 - 1
		_b = H / self.rc
		_R = _a * (1 - _b) * np.tan( np.deg2rad(zenith) ) - _a * (_b - _a/2) * np.tan( np.deg2rad(zenith) )**3
		return np.degrees(_R)



	def corbard(self, n1, 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'.
		'''
		if reduced_height:
			H = reduced_height
		else:
			H = self.H

		_a = n1 - 1
		_b = H / self.rc
		_R = _a * ( (2 - _a) / (np.sqrt( 2*_b - _a )) ) * np.sin( np.deg2rad(zenith) ) * self.psi( (np.cos(np.deg2rad(zenith))) / np.sqrt(2*_b - _a) )
		return np.degrees(_R)



	def refractionIntegral(self, 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 not R0:
			R0 = self.rc
		if useStandardAtmosphere:
			_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('./US1976 Standard Atmosphere/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 

		# Calculation of the zenith angle for given height and refractive index, from the refractive invariant
		_z 		= np.arcsin( (n0 * R0) / (_n_h * _R) * np.sin(z0) )
		
		# Calculation of log space derivative of n and r
		_lnn = np.log(_n_h)
		_lnr = np.log(_R)
		_dlnn = np.asarray([_lnn[i+1] - _lnn[i] for i in range(len(_lnn)-1)])
		_dlnr = np.asarray([_lnr[i+1] - _lnr[i] for i in range(len(_lnr)-1)])
		_dz   = np.asarray([(_z[i] - _z[i+1]) for i in range(len(_z)-1)])

		# Calculation of the integrand in eq. 3 in Auer & Standish
		_Integrand = -1*(_dlnn/_dlnr) / (1 + _dlnn/_dlnr) * _dz


		return np.degrees(np.sum(_Integrand))



	def psi(self, x):
		_f   = lambda a: np.exp(-a**2)
		_int = integrate.quad(_f, x, np.inf)
		return np.exp(x**2) * _int[0]



	def set_rc(self):
		'''
			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.
		'''
		_A   = 6378.137 # km
		_B   = 6356.752 # km
		_phi = np.deg2rad(self.lat)
		
		_rc0 = (_A*_B)**2 / (_A**2 * np.cos(_phi)**2 + _B**2 * np.sin(_phi)**2)**(3/2)
		return _rc0 + self.h/1000



	def set_H(self, P, rho):
		'''
			Calculates the reduced height of the atmosphere, assuming ideal gas law.
		'''
		_phi = np.deg2rad(self.lat)
		_c1 = 5.2790414e-3
		_c2 = 2.32718e-5
		_c3 = 1.262e-7
		_c4 = 7e-10
		_g0 = 9.780327 # m/s^2
		_g0_local = _g0 * (1 + _c1 * np.sin(_phi)**2 + _c2 * np.sin(_phi)**4 + _c3 * np.sin(_phi)**6 + _c4 * np.sin(_phi)**8)
		_g = _g0_local - (3.0877e-6 - 4.3e-9 * np.sin(_phi)**2) * self.h + 7.2e-13 * self.h**2
		return P / (rho * _g)

	def H_isotherm(self, T):
		'''
			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(self.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 cassiniError(self, n, dn, zenith, dz=0, reduced_height=None):
		'''
		Calculates the uncertainty in the refraction angle in degrees.
		'''
		if reduced_height:
			_H = reduced_height
		else:
			_H = self.H
		_zenith = np.deg2rad(zenith)
		_dz = np.deg2rad(dz)
		_R = self.rc
		_dRdn = (_R * np.sin(_zenith)) / ((_R + _H) * np.sqrt(1 - ((n**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2)))

		_dRdz = (_R * np.cos(_zenith)) / (_R + _H) * ( n / np.sqrt(1-(n**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2) - 1 / np.sqrt(1-(_R**2 * np.sin(_zenith)**2)/(_R + _H)**2) )

		_dR = np.sqrt( (_dRdn * dn)**2 + (_dRdz * _dz)**2)

		return np.degrees(_dR)