yosata/AtmosphericModel
1
0 Fork 0
mirror of https://gitlab.astro-wise.org/micado/atmosphericmodel.git synced 2026-10-29 20:51:05 +00:00
Code Issues Projects Releases Wiki Activity
Files
eee6181d26d210bc0a5ed80a71b4e569ea303311
AtmosphericModel/atmosphericModel.py
T
Joost van den Born 10bf0fa17d First commit
2019-04-04 11:47:43 +02:00

809 lines
24 KiB
Python
Raw Blame History

#!/usr/bin/env Python3
import numpy as np
from scipy import integrate
'''
@ 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 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
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) )
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_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):
'''
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)
class dispersion:
def __init__(self, lat, h, P, rho):
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
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)
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
'''
if not R0:
R0 = self.rc
if useStandardAtmosphere:
_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('1976USSA.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]:
_refractionIntegralList = []
_nList = []
for _i in range(len(_h)):
_R = R0 + _h[_i]*1000
_n_h = 1 + (_n - 1) / _Dens[0] * _Dens[_i] # Gladstone-Dale relation between refractive index and density
_nList.append(_n_h)
_refractionIntegralList.append( (_n * R0 * np.sin(z0))/(_n_h*np.sqrt((_n_h*_R)**2 - (_n * R0 * np.sin(z0))**2 )) )
_dn = (_nList[0] - _nList[-1]) / ((len(_h)-1))
_rList.append(np.sum(np.degrees(np.asarray(_refractionIntegralList) * _dn)))
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 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 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)
class refraction:
def __init__(self, lat, h, P, rho):
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
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)
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('1976USSA.txt', unpack=True, skip_header=6)
else:
_Dens = rhoData
_h = heightData
_refractionIntegralList = []
_nList = []
z0 = np.deg2rad(z0)
# Integration of the refraction integral over the atmosphere layers.
for _i in range(len(_h)):
_R = R0 + _h[_i]*1000
_n_h = 1 + (n0 - 1) / _Dens[0] * _Dens[_i] # Gladstone-Dale relation between refractive index and density
_nList.append(_n_h)
_refractionIntegralList.append( (n0 * R0 * np.sin(z0))/(_n_h*np.sqrt((_n_h*_R)**2 - (n0 * R0 * np.sin(z0))**2 )) )
_dn = (_nList[0] - _nList[-1]) / ((len(_h)-1))
return np.sum(np.degrees(np.asarray(_refractionIntegralList) * _dn))
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 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 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)
Reference in new issue View Git Blame Copy Permalink