mirror of
https://gitlab.astro-wise.org/micado/atmosphericmodel.git
synced 2026-10-29 20:51:05 +00:00
Add error propagation to the models
This commit is contained in:
Joost van den Born committed
1 parent
10bf0fa17d
commit
a81c17a4b4
5 files changed
+646
-60
No files matched your search
+497
-34
@@ -17,14 +17,13 @@ from scipy import integrate
|
||||
'''
|
||||
|
||||
|
||||
class Ciddor1996:
|
||||
def __init__(self, p=101325, T=288.15, H=0, xc=450):
|
||||
self.p = p # Pressure in Pa
|
||||
self.T = T # Temperature in Kelvin
|
||||
self.H = H # Relative humidity on a scale from 0 to 1
|
||||
self.xc = xc # ppm of CO2
|
||||
pass
|
||||
|
||||
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):
|
||||
'''
|
||||
@@ -191,15 +190,216 @@ class Ciddor1996:
|
||||
return np.sqrt( (1 + 2 * _Ltotal) / (1 - _Ltotal) )
|
||||
|
||||
|
||||
class error(Ciddor1996):
|
||||
def __init__(self, p=101325, T=288.15, H=0, xc=450, dT=0.2, dP=20, dH=0.02, dCO2=20, dl=0.001):
|
||||
Ciddor1996.__init__(self, p, T, H)
|
||||
self.dT = 0#dT # in K
|
||||
self.dP = 0#dP # in Pa
|
||||
self.dH = 0#dH # relative humidity between 0 and 1
|
||||
self.dCO2 = 0#dCO2 # in ppm
|
||||
self.dl = 0#dl # in micron
|
||||
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):
|
||||
'''
|
||||
@@ -384,6 +584,7 @@ class error(Ciddor1996):
|
||||
|
||||
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:
|
||||
@@ -443,24 +644,27 @@ class error(Ciddor1996):
|
||||
_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:
|
||||
def __init__(self, lat, h, P, rho):
|
||||
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.rc() * 1000 # Radius of curvature of the earth at the observer
|
||||
self.H = self.H(P, rho) # Reduced height of the atmosphere assuming ideal gas law
|
||||
# self.H2 = self.H_isotherm(T) # Reduced height of isotherm atmosphere
|
||||
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.
|
||||
'''
|
||||
# _theta1 = np.degrees( np.arcsin( np.sin( np.deg2rad(zenith) / n1) ) )
|
||||
# _theta2 = np.degrees( np.arcsin( np.sin( np.deg2rad(zenith) / n2) ) )
|
||||
_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)
|
||||
@@ -555,7 +759,7 @@ class dispersion:
|
||||
if not R0:
|
||||
R0 = self.rc
|
||||
if useStandardAtmosphere:
|
||||
_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('1976USSA.txt', unpack=True, skip_header=6)
|
||||
_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('./US1976 Standard Atmosphere/si1py_modified.txt', unpack=True, skip_header=6)
|
||||
else:
|
||||
_Dens = rhoData
|
||||
_h = heightData
|
||||
@@ -589,7 +793,7 @@ class dispersion:
|
||||
return np.exp(x**2) * _int[0]
|
||||
|
||||
|
||||
def rc(self):
|
||||
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
|
||||
@@ -603,7 +807,7 @@ class dispersion:
|
||||
return _rc0 + self.h/1000
|
||||
|
||||
|
||||
def H(self, P, rho):
|
||||
def set_H(self, P, rho):
|
||||
'''
|
||||
Calculates the reduced height of the atmosphere, assuming ideal gas law.
|
||||
'''
|
||||
@@ -634,21 +838,262 @@ class dispersion:
|
||||
|
||||
return (_kb*T) / (_m * _g0_local)
|
||||
|
||||
class refraction:
|
||||
def __init__(self, lat, h, P, rho):
|
||||
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.rc() * 1000 # Radius of curvature of the earth at the observer
|
||||
self.H = self.H(P, rho) # Reduced height of the atmosphere assuming ideal gas law
|
||||
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.
|
||||
'''
|
||||
# _theta1 = np.degrees( np.arcsin( np.sin( np.deg2rad(zenith) / n) ) )
|
||||
_R = np.arcsin(n * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
|
||||
return np.degrees(_R)
|
||||
|
||||
@@ -732,7 +1177,7 @@ class refraction:
|
||||
if not R0:
|
||||
R0 = self.rc
|
||||
if useStandardAtmosphere:
|
||||
_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('1976USSA.txt', unpack=True, skip_header=6)
|
||||
_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('./US1976 Standard Atmosphere/si1py_modified.txt', unpack=True, skip_header=6)
|
||||
else:
|
||||
_Dens = rhoData
|
||||
_h = heightData
|
||||
@@ -760,7 +1205,7 @@ class refraction:
|
||||
|
||||
|
||||
|
||||
def rc(self):
|
||||
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
|
||||
@@ -775,7 +1220,7 @@ class refraction:
|
||||
|
||||
|
||||
|
||||
def H(self, P, rho):
|
||||
def set_H(self, P, rho):
|
||||
'''
|
||||
Calculates the reduced height of the atmosphere, assuming ideal gas law.
|
||||
'''
|
||||
@@ -806,3 +1251,21 @@ class refraction:
|
||||
|
||||
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)
|
||||
Reference in new issue
Block a user