mirror of
https://gitlab.astro-wise.org/micado/atmosphericmodel.git
synced 2026-10-29 20:51:05 +00:00
Added Mathars exponential refraction model
This commit is contained in:
yosata committed
1 parent
e5e20e3584
commit
b024ae5d3a
4 files changed
+768
-517
No files matched your search
+456
-329
@@ -14,6 +14,7 @@ from .ciddorModel import Observatory
|
|||||||
The actual refraction is calculated by any of the following models:
|
The actual refraction is calculated by any of the following models:
|
||||||
- Plane parallel atmosphere
|
- Plane parallel atmosphere
|
||||||
- Cassini spherical homogeneous atmosphere
|
- Cassini spherical homogeneous atmosphere
|
||||||
|
- Mathar's Barometric exponential model
|
||||||
- Other methods described in Corbard et al., 2019.
|
- Other methods described in Corbard et al., 2019.
|
||||||
'''
|
'''
|
||||||
|
|
||||||
@@ -79,44 +80,6 @@ class dispersion(Observatory):
|
|||||||
self.lat = lat # Latitude of the observer in degrees
|
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.rc = self.set_rc() * 1000 # Radius of curvature of the earth at the observer
|
||||||
|
|
||||||
def setReducedHeight(self, p, rho):
|
|
||||||
'''
|
|
||||||
Sets the reduced atmospheric height.
|
|
||||||
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
p : float
|
|
||||||
Atmospheric pressure in Pascal
|
|
||||||
rho : float
|
|
||||||
Density of the atmospheric air at the observer
|
|
||||||
'''
|
|
||||||
self.H = self.set_H(p, rho) # Reduced height of the atmosphere assuming ideal gas law
|
|
||||||
|
|
||||||
|
|
||||||
def planeParallel(self, n1, n2, zenith):
|
|
||||||
'''
|
|
||||||
Refraction based on a flat atmosphere, easily derived from Snell's law.
|
|
||||||
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
n1 : float
|
|
||||||
Refractive index at the shortest wavelength of interest
|
|
||||||
n2 : float
|
|
||||||
Refractive index at the longest wavelength of interest
|
|
||||||
zenith : float
|
|
||||||
Angle of observation in degrees.
|
|
||||||
zenith = 0 means the observation is directly overhead.
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
dR : float
|
|
||||||
Atmospheric dispersion in degrees
|
|
||||||
'''
|
|
||||||
_R1 = np.arcsin(n1 * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
|
|
||||||
_R2 = np.arcsin(n2 * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
|
|
||||||
return np.degrees(_R1) - np.degrees(_R2)
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
def cassini(self, n1, n2, zenith, reduced_height=None):
|
def cassini(self, n1, n2, zenith, reduced_height=None):
|
||||||
'''
|
'''
|
||||||
@@ -150,296 +113,6 @@ class dispersion(Observatory):
|
|||||||
return np.degrees(_R1) - np.degrees(_R2)
|
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.
|
|
||||||
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
n1 : float
|
|
||||||
Refractive index at the shortest wavelength of interest
|
|
||||||
n2 : float
|
|
||||||
Refractive index at the longest wavelength of interest
|
|
||||||
zenith : float
|
|
||||||
Angle of observation in degrees.
|
|
||||||
zenith = 0 means the observation is directly overhead.
|
|
||||||
reduced_height : float (optional)
|
|
||||||
Reduced height of the atmosphere. Default value is calculated
|
|
||||||
from the object attributes.
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
dR : float
|
|
||||||
Atmospheric dispersion in degrees
|
|
||||||
'''
|
|
||||||
if reduced_height:
|
|
||||||
H = reduced_height
|
|
||||||
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.
|
|
||||||
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
n1 : float
|
|
||||||
Refractive index at the shortest wavelength of interest
|
|
||||||
n2 : float
|
|
||||||
Refractive index at the longest wavelength of interest
|
|
||||||
zenith : float
|
|
||||||
Angle of observation in degrees.
|
|
||||||
zenith = 0 means the observation is directly overhead.
|
|
||||||
reduced_height : float (optional)
|
|
||||||
Reduced height of the atmosphere. Default value is calculated
|
|
||||||
from the object attributes.
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
dR : float
|
|
||||||
Atmospheric dispersion in degrees
|
|
||||||
'''
|
|
||||||
if reduced_height:
|
|
||||||
H = reduced_height
|
|
||||||
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'.
|
|
||||||
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
n1 : float
|
|
||||||
Refractive index at the shortest wavelength of interest
|
|
||||||
n2 : float
|
|
||||||
Refractive index at the longest wavelength of interest
|
|
||||||
zenith : float
|
|
||||||
Angle of observation in degrees.
|
|
||||||
zenith = 0 means the observation is directly overhead.
|
|
||||||
reduced_height : float (optional)
|
|
||||||
Reduced height of the atmosphere. Default value is calculated
|
|
||||||
from the object attributes.
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
dR : float
|
|
||||||
Atmospheric dispersion in degrees
|
|
||||||
'''
|
|
||||||
if reduced_height:
|
|
||||||
H = reduced_height
|
|
||||||
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.
|
|
||||||
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
n1 : float
|
|
||||||
Refractive index at the shortest wavelength of interest
|
|
||||||
n2 : float
|
|
||||||
Refractive index at the longest wavelength of interest
|
|
||||||
z0 : float
|
|
||||||
Angle of observation in degrees.
|
|
||||||
zenith = 0 means the observation is directly overhead.
|
|
||||||
R0 : float (optional)
|
|
||||||
Local radius of curvature of the Earth in km.
|
|
||||||
useStandardAtmosphere : bool
|
|
||||||
Flag to use the included USSA1976 data
|
|
||||||
heightData : 1D Array
|
|
||||||
If "useStandardAtmosphere" is False, the function will ask
|
|
||||||
for alternative atmospheric data to integrate over. This
|
|
||||||
gives the altitudinal points
|
|
||||||
rhoData : 1D Array
|
|
||||||
If useStandardAtmosphere is False, the function will ask
|
|
||||||
for alternative atmospheric data to integrate over. This
|
|
||||||
gives the atmospheric density at the points defined by
|
|
||||||
"heightData".
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
dR : float
|
|
||||||
Atmospheric dispersion in degrees
|
|
||||||
'''
|
|
||||||
if not R0:
|
|
||||||
R0 = self.rc
|
|
||||||
if useStandardAtmosphere:
|
|
||||||
_datafileLocation = os.path.join(os.path.dirname(__file__), 'data/1976USSA.txt')
|
|
||||||
_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt(_datafileLocation, 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):
|
|
||||||
'''
|
|
||||||
Submethod for corbard(). Described explicitly in Corbard et al.,2019.
|
|
||||||
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
x : float
|
|
||||||
Input of the function
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
y : float
|
|
||||||
Output of the function
|
|
||||||
'''
|
|
||||||
_f = lambda a: np.exp(-a**2)
|
|
||||||
_int = integrate.quad(_f, x, np.inf)
|
|
||||||
return np.exp(x**2) * _int[0]
|
|
||||||
|
|
||||||
|
|
||||||
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.
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
rc : float
|
|
||||||
Local radius of curvature of the Earth in kilometers.
|
|
||||||
'''
|
|
||||||
_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.
|
|
||||||
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
p : float
|
|
||||||
Atmospheric pressure in Pascal
|
|
||||||
rho : float
|
|
||||||
Density of the atmospheric air at the observer
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
H : float
|
|
||||||
Reduced height of the atmosphere in meters
|
|
||||||
'''
|
|
||||||
_phi = np.deg2rad(self.lat)
|
|
||||||
_c1 = 5.2790414e-3
|
|
||||||
_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.
|
|
||||||
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
T : float
|
|
||||||
Atmospheric temperature in Kelvin
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
H : float
|
|
||||||
Reduced height of the atmosphere in meters
|
|
||||||
'''
|
|
||||||
_kb = 1.38064852e-23 # kg m2 s-2 K-1
|
|
||||||
_m = 4.809651698e-26 # kg (weight per molecule of air, assuming dry air)
|
|
||||||
|
|
||||||
_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,
|
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):
|
dRH=0.02, dCO2=20, dz=0, lat=None, h=None):
|
||||||
'''
|
'''
|
||||||
@@ -705,4 +378,458 @@ class dispersion(Observatory):
|
|||||||
# Final error:
|
# 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)
|
_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)
|
return np.degrees(_dispersionError)
|
||||||
|
|
||||||
|
|
||||||
|
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'.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
n1 : float
|
||||||
|
Refractive index at the shortest wavelength of interest
|
||||||
|
n2 : float
|
||||||
|
Refractive index at the longest wavelength of interest
|
||||||
|
zenith : float
|
||||||
|
Angle of observation in degrees.
|
||||||
|
zenith = 0 means the observation is directly overhead.
|
||||||
|
reduced_height : float (optional)
|
||||||
|
Reduced height of the atmosphere. Default value is calculated
|
||||||
|
from the object attributes.
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
dR : float
|
||||||
|
Atmospheric dispersion in degrees
|
||||||
|
'''
|
||||||
|
if reduced_height:
|
||||||
|
H = reduced_height
|
||||||
|
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 gammaCoefficients(self, n0, reduced_height=None):
|
||||||
|
'''
|
||||||
|
Submethod for matharExponential(). Calculates the gamma
|
||||||
|
coefficients.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
n0 : float
|
||||||
|
Refractive index at the wavelength of interest at the observer
|
||||||
|
reduced_height : float (optional)
|
||||||
|
Reduced height of the atmosphere. Default value is calculated
|
||||||
|
from the object attributes.
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
g : array
|
||||||
|
List of values, corresponding to the gamma coefficients.
|
||||||
|
g[0] corresponds to g_1, g[1] to g_3, g[2] to g_5, and so on.
|
||||||
|
'''
|
||||||
|
if reduced_height:
|
||||||
|
H = reduced_height
|
||||||
|
else:
|
||||||
|
H = self.H
|
||||||
|
|
||||||
|
_chi0 = n0**2 - 1
|
||||||
|
_Khat = H / self.rc
|
||||||
|
|
||||||
|
_g1 = (1 - 5/24 * _Khat * _chi0**2 - 3/4 * _chi0 - 3/8 * _Khat**2 * _chi0 \
|
||||||
|
- _Khat - 6 * _Khat**3 + 5/8 * _chi0**2 + 63/128 * _chi0**4 + 3/8 \
|
||||||
|
* _Khat * _chi0 + 35/256 * _Khat * _chi0**3 + 2 * _Khat**2 + 5/36 \
|
||||||
|
* _Khat**2 * _chi0**2 - 35/64 * _chi0**3 + 9/16 * _Khat**3 * _chi0 \
|
||||||
|
- 120 * _Khat**5 + 24 * _Khat**4 - 231 / 512 * _chi0**5 - 63 / 640 \
|
||||||
|
* _Khat * _chi0**4 - 35 / 512 * _Khat**2 * _chi0**3 - 5 / 36 * _Khat**3\
|
||||||
|
* _chi0**2 - 9 / 8 * _Khat**4 * _chi0 )
|
||||||
|
|
||||||
|
_g3 = (35/96 * _chi0**3 - 5/12 * _chi0**2 + 333 * _Khat**4 * _chi0 - 2400 \
|
||||||
|
* _Khat**5 + 3605/4608 * _Khat**2 * _chi0**3 + 63/256 * _Khat * _chi0**4\
|
||||||
|
- 819/16 * _Khat**3 * _chi0 - 35/18 * _Khat**2 * _chi0**2 - 35/96 \
|
||||||
|
* _Khat * _chi0**3 + 69/8 * _Khat**2 * _chi0 + 5/8 * _Khat * _chi0**2 \
|
||||||
|
- 3/2 * _Khat * _chi0 + 1355/216 * _Khat**3 * _chi0**2 + 10 * _Khat**2 \
|
||||||
|
+ 336 * _Khat**4 - 2 * _Khat + _chi0/2 - 21/64 * _chi0**4 - 54 * _Khat**3 \
|
||||||
|
+ 77/256 * _chi0**5 ) * 1/2
|
||||||
|
|
||||||
|
_g5 = (-96 * _Khat**3 + 8 * _Khat**2 - 7/24 * _chi0**3 + 65/4 * _Khat**2 \
|
||||||
|
* _chi0**2 + 35/48 * _Khat * _chi0**3 + 30 * _Khat**2 * _chi0 - 5/3 \
|
||||||
|
* _Khat * _chi0**2 - 3 * _Khat * _chi0 + 20091/8 * _Khat**4 * _chi0 \
|
||||||
|
- 32605/216 * _Khat**3 * _chi0**2 - 35/9 * _Khat**2 * _chi0**3 + 984 \
|
||||||
|
* _Khat**4 + _chi0**2/3 - 10200 * _Khat**5 - 77/320 * _chi0**5 + 21/80 \
|
||||||
|
* _chi0**4 - 1089/4 * _Khat**3 * _chi0 - 7/16 * _Khat * _chi0**4 ) * 3/8
|
||||||
|
|
||||||
|
_g7 = (-17520 * _Khat**5 - 9/40 * _chi0**4 - 8755/12 * _Khat**3 * _chi0**2 \
|
||||||
|
+ 1056 * _Khat**4 + _chi0**3/4 + 23517/4 * _Khat**4 * _chi0 - 11/3 \
|
||||||
|
* _Khat * _chi0**2 + 21 * _Khat**2 * _chi0 - 7/4 * _Khat * _chi0**3 \
|
||||||
|
+ 505/9 * _Khat**2 * _chi0**2 - 801/2 * _Khat**3 * _chi0 + 33/160 \
|
||||||
|
* _chi0**5 - 48 * _Khat**3 + 63/80 * _Khat * _chi0**4 + 2345/96 \
|
||||||
|
* _Khat**2 * _chi0**3) * 5/16
|
||||||
|
|
||||||
|
_g9 = (_chi0**4 / 5 + 384 * _Khat**4 - 13440 * _Khat**5 + 5562 * _Khat**4 \
|
||||||
|
* _chi0 + 3125/36 * _Khat**2 * _chi0**3 - 11/60 * _chi0**5 - 25/6 \
|
||||||
|
* _Khat * _chi0**3 - 26950/27 * _Khat**3 * _chi0**2 + 340/9 * _Khat**2 \
|
||||||
|
* _chi0**2 - 180 * _Khat**3 * _chi0 - 9/5 * _Khat * _chi0**4) * 35/128
|
||||||
|
|
||||||
|
g = np.array([_g1, _g3, _g5, _g7, _g9])
|
||||||
|
g*= n0*_chi0/2
|
||||||
|
|
||||||
|
return g
|
||||||
|
|
||||||
|
|
||||||
|
def H_isotherm(self, T):
|
||||||
|
'''
|
||||||
|
Calculates the reduced height of the atmosphere,
|
||||||
|
assuming an isotherm distribution.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
T : float
|
||||||
|
Atmospheric temperature in Kelvin
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
H : float
|
||||||
|
Reduced height of the atmosphere in meters
|
||||||
|
'''
|
||||||
|
_kb = 1.38064852e-23 # kg m2 s-2 K-1
|
||||||
|
_m = 4.809651698e-26 # kg (weight per molecule of air, assuming dry air)
|
||||||
|
|
||||||
|
_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 matharExponential(self, n1, n2, zenith, reduced_height=None):
|
||||||
|
'''
|
||||||
|
Model on atmospheric refraction described by Richard Mathar in
|
||||||
|
"A Barometric Exponential Model of the Atmosphere’s Refractive Index: ZenithAngles and Second Order Aberration in the Entrance Pupil" (2020)
|
||||||
|
ArXiv id: 2004.11808 (https://arxiv.org/abs/2004.11808)
|
||||||
|
|
||||||
|
The code in this function is adapted from eqs 3, 16, 17 and 18 in this
|
||||||
|
work. It is an adaptation of the RefrExpo.cxx C++ file written by Mathar.
|
||||||
|
|
||||||
|
This model assumes the refractive index, n, is related to the susceptibility,
|
||||||
|
chi, following: n(r) = sqrt(1+chi(r)) (eq. 9).
|
||||||
|
|
||||||
|
R. Mathar then assumes that the susceptibility is proportional to the
|
||||||
|
air density and reduced height, following eq. 10:
|
||||||
|
chi(r) = chi0 * exp(-(r-rho)/K).
|
||||||
|
|
||||||
|
The refraction is described using a Taylor expansion; R in terms of
|
||||||
|
the odd powers of the tangent of the observed zenith angle (eq. 3), not
|
||||||
|
unlike Oriani's theorem. This requires coefficients g_i.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
n1 : float
|
||||||
|
Refractive index at the wavelength of interest at the observer
|
||||||
|
n2 : float
|
||||||
|
Refractive index at the longest wavelength of interest at the
|
||||||
|
observer
|
||||||
|
zenith : float
|
||||||
|
Angle of observation in degrees.
|
||||||
|
zenith = 0 means the observation is directly overhead.
|
||||||
|
reduced_height : float (optional)
|
||||||
|
Reduced height of the atmosphere. Default value is calculated
|
||||||
|
from the object attributes.
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
dR : float
|
||||||
|
Atmospheric refraction in degrees
|
||||||
|
'''
|
||||||
|
if reduced_height:
|
||||||
|
_H = reduced_height
|
||||||
|
else:
|
||||||
|
_H = self.H
|
||||||
|
|
||||||
|
_zenith = np.deg2rad(zenith)
|
||||||
|
_rc = self.rc # rho
|
||||||
|
|
||||||
|
# Atmospheric refraction for the first wavelength
|
||||||
|
_gamma1 = self.gammaCoefficients(n1, reduced_height=_H)
|
||||||
|
_R1 = sum([_gamma1[i]*np.tan(_zenith)**(2*i+1) for i in range(len(_gamma1))])
|
||||||
|
|
||||||
|
# Atmospheric refraction for the second wavelength
|
||||||
|
_gamma2 = self.gammaCoefficients(n2, reduced_height=_H)
|
||||||
|
_R2 = sum([_gamma2[i]*np.tan(_zenith)**(2*i+1) for i in range(len(_gamma2))])
|
||||||
|
|
||||||
|
# The atmospheric dispersion in radians
|
||||||
|
_dR = _R1 - _R2
|
||||||
|
|
||||||
|
return np.degrees(_dR)
|
||||||
|
|
||||||
|
|
||||||
|
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.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
n1 : float
|
||||||
|
Refractive index at the shortest wavelength of interest
|
||||||
|
n2 : float
|
||||||
|
Refractive index at the longest wavelength of interest
|
||||||
|
zenith : float
|
||||||
|
Angle of observation in degrees.
|
||||||
|
zenith = 0 means the observation is directly overhead.
|
||||||
|
reduced_height : float (optional)
|
||||||
|
Reduced height of the atmosphere. Default value is calculated
|
||||||
|
from the object attributes.
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
dR : float
|
||||||
|
Atmospheric dispersion in degrees
|
||||||
|
'''
|
||||||
|
if reduced_height:
|
||||||
|
H = reduced_height
|
||||||
|
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 planeParallel(self, n1, n2, zenith):
|
||||||
|
'''
|
||||||
|
Refraction based on a flat atmosphere, easily derived from Snell's law.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
n1 : float
|
||||||
|
Refractive index at the shortest wavelength of interest
|
||||||
|
n2 : float
|
||||||
|
Refractive index at the longest wavelength of interest
|
||||||
|
zenith : float
|
||||||
|
Angle of observation in degrees.
|
||||||
|
zenith = 0 means the observation is directly overhead.
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
dR : float
|
||||||
|
Atmospheric dispersion in degrees
|
||||||
|
'''
|
||||||
|
_R1 = np.arcsin(n1 * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
|
||||||
|
_R2 = np.arcsin(n2 * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
|
||||||
|
return np.degrees(_R1) - np.degrees(_R2)
|
||||||
|
|
||||||
|
|
||||||
|
def psi(self, x):
|
||||||
|
'''
|
||||||
|
Submethod for corbard(). Described explicitly in Corbard et al.,2019.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
x : float
|
||||||
|
Input of the function
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
y : float
|
||||||
|
Output of the function
|
||||||
|
'''
|
||||||
|
_f = lambda a: np.exp(-a**2)
|
||||||
|
_int = integrate.quad(_f, x, np.inf)
|
||||||
|
return np.exp(x**2) * _int[0]
|
||||||
|
|
||||||
|
|
||||||
|
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.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
n1 : float
|
||||||
|
Refractive index at the shortest wavelength of interest
|
||||||
|
n2 : float
|
||||||
|
Refractive index at the longest wavelength of interest
|
||||||
|
z0 : float
|
||||||
|
Angle of observation in degrees.
|
||||||
|
zenith = 0 means the observation is directly overhead.
|
||||||
|
R0 : float (optional)
|
||||||
|
Local radius of curvature of the Earth in km.
|
||||||
|
useStandardAtmosphere : bool
|
||||||
|
Flag to use the included USSA1976 data
|
||||||
|
heightData : 1D Array
|
||||||
|
If "useStandardAtmosphere" is False, the function will ask
|
||||||
|
for alternative atmospheric data to integrate over. This
|
||||||
|
gives the altitudinal points
|
||||||
|
rhoData : 1D Array
|
||||||
|
If useStandardAtmosphere is False, the function will ask
|
||||||
|
for alternative atmospheric data to integrate over. This
|
||||||
|
gives the atmospheric density at the points defined by
|
||||||
|
"heightData".
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
dR : float
|
||||||
|
Atmospheric dispersion in degrees
|
||||||
|
'''
|
||||||
|
if not R0:
|
||||||
|
R0 = self.rc
|
||||||
|
if useStandardAtmosphere:
|
||||||
|
_datafileLocation = os.path.join(os.path.dirname(__file__), 'data/1976USSA.txt')
|
||||||
|
_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt(_datafileLocation, 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 setReducedHeight(self, p, rho):
|
||||||
|
'''
|
||||||
|
Sets the reduced atmospheric height.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
p : float
|
||||||
|
Atmospheric pressure in Pascal
|
||||||
|
rho : float
|
||||||
|
Density of the atmospheric air at the observer
|
||||||
|
'''
|
||||||
|
self.H = self.set_H(p, rho) # Reduced height of the atmosphere assuming ideal gas law
|
||||||
|
|
||||||
|
|
||||||
|
def set_H(self, p, rho):
|
||||||
|
'''
|
||||||
|
Calculates the reduced height of the atmosphere, assuming ideal gas law.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
p : float
|
||||||
|
Atmospheric pressure in Pascal
|
||||||
|
rho : float
|
||||||
|
Density of the atmospheric air at the observer
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
H : float
|
||||||
|
Reduced height of the atmosphere in meters
|
||||||
|
'''
|
||||||
|
_phi = np.deg2rad(self.lat)
|
||||||
|
_c1 = 5.2790414e-3
|
||||||
|
_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 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.
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
rc : float
|
||||||
|
Local radius of curvature of the Earth in kilometers.
|
||||||
|
'''
|
||||||
|
_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 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.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
n1 : float
|
||||||
|
Refractive index at the shortest wavelength of interest
|
||||||
|
n2 : float
|
||||||
|
Refractive index at the longest wavelength of interest
|
||||||
|
zenith : float
|
||||||
|
Angle of observation in degrees.
|
||||||
|
zenith = 0 means the observation is directly overhead.
|
||||||
|
reduced_height : float (optional)
|
||||||
|
Reduced height of the atmosphere. Default value is calculated
|
||||||
|
from the object attributes.
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
dR : float
|
||||||
|
Atmospheric dispersion in degrees
|
||||||
|
'''
|
||||||
|
if reduced_height:
|
||||||
|
H = reduced_height
|
||||||
|
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]
|
||||||
+302
-184
@@ -9,11 +9,11 @@ from .ciddorModel import Observatory
|
|||||||
@ Contact: born@astron.nl
|
@ Contact: born@astron.nl
|
||||||
@ Description: This script offers a complete model of the atmospheric refraction
|
@ Description: This script offers a complete model of the atmospheric refraction
|
||||||
and dispersion. It uses the Ciddor1996 paper to calculate the
|
and dispersion. It uses the Ciddor1996 paper to calculate the
|
||||||
index of refraction for a wide range of conditions, with corresponding
|
index of refraction for a wide range of conditions, with corresponding error propagation if needed.
|
||||||
error propagation if needed.
|
|
||||||
The actual refraction is calculated by any of the following models:
|
The actual refraction is calculated by any of the following models:
|
||||||
- Plane parallel atmosphere
|
- Plane parallel atmosphere
|
||||||
- Cassini spherical homogeneous atmosphere
|
- Cassini spherical homogeneous atmosphere
|
||||||
|
- Mathar's Barometric exponential model
|
||||||
- Other methods described in Corbard et al., 2019.
|
- Other methods described in Corbard et al., 2019.
|
||||||
'''
|
'''
|
||||||
|
|
||||||
@@ -80,41 +80,6 @@ class refraction(Observatory):
|
|||||||
self.rc = self.set_rc() * 1000 # Radius of curvature of the earth at the observer
|
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
|
# self.H2 = self.H_isotherm(T) # Reduced height of isotherm atmosphere
|
||||||
|
|
||||||
def setReducedHeight(self, p, rho):
|
|
||||||
'''
|
|
||||||
Sets the reduced atmospheric height.
|
|
||||||
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
p : float
|
|
||||||
Atmospheric pressure in Pascal
|
|
||||||
rho : float
|
|
||||||
Density of the atmospheric air at the observer
|
|
||||||
'''
|
|
||||||
self.H = self.set_H(p, rho) # Reduced height of the atmosphere assuming ideal gas law
|
|
||||||
|
|
||||||
|
|
||||||
def planeParallel(self, n1, zenith):
|
|
||||||
'''
|
|
||||||
Refraction based on a flat atmosphere, easily derived from Snell's law.
|
|
||||||
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
n1 : float
|
|
||||||
Refractive index at the wavelength of interest
|
|
||||||
zenith : float
|
|
||||||
Angle of observation in degrees.
|
|
||||||
zenith = 0 means the observation is directly overhead.
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
R : float
|
|
||||||
Atmospheric refraction in degrees
|
|
||||||
'''
|
|
||||||
_R = np.arcsin(n * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
|
|
||||||
return np.degrees(_R)
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
def cassini(self, n1, zenith, reduced_height=None):
|
def cassini(self, n1, zenith, reduced_height=None):
|
||||||
'''
|
'''
|
||||||
@@ -145,72 +110,46 @@ class refraction(Observatory):
|
|||||||
return np.degrees(_R1)
|
return np.degrees(_R1)
|
||||||
|
|
||||||
|
|
||||||
|
def cassiniError(self, n, dn, zenith, dz=0, reduced_height=None):
|
||||||
def tan5(self, n1, zenith, reduced_height=None):
|
|
||||||
'''
|
'''
|
||||||
Same as oriani's theorem, but including a higher order term.
|
Calculates the uncertainty in the refraction angle in degrees.
|
||||||
Found in Corbard et al., 2019.
|
See Observatory().n_tph() and Observatory().dn_tph() on how to
|
||||||
|
calculate the refractive index and it's uncertainty.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
n1 : float
|
n : float
|
||||||
Refractive index at the wavelength of interest
|
Refractive index of air
|
||||||
zenith : float
|
dn : float
|
||||||
Angle of observation in degrees.
|
Uncertainty in the refractive index of air
|
||||||
zenith = 0 means the observation is directly overhead.
|
zenith : float
|
||||||
reduced_height : float (optional)
|
Angle of observation in degrees.
|
||||||
Reduced height of the atmosphere. Default value is calculated
|
zenith = 0 means the observation is directly overhead.
|
||||||
from the object attributes.
|
dz : float (optional)
|
||||||
|
Uncertainty in the zenith angle, in degrees.
|
||||||
|
reduced_height : float (optional)
|
||||||
|
Reduced height of the atmosphere. Default value is calculated
|
||||||
|
from the object attributes.
|
||||||
|
|
||||||
Returns
|
Returns
|
||||||
-------
|
-------
|
||||||
R : float
|
dR : float
|
||||||
Atmospheric refraction in degrees
|
Uncertainty in the atmospheric refraction in degrees.
|
||||||
'''
|
'''
|
||||||
if reduced_height:
|
if reduced_height:
|
||||||
H = reduced_height
|
_H = reduced_height
|
||||||
else:
|
else:
|
||||||
H = self.H
|
_H = self.H
|
||||||
|
_zenith = np.deg2rad(zenith)
|
||||||
_a = n1 - 1
|
_dz = np.deg2rad(dz)
|
||||||
_b = H / self.rc
|
_R = 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
|
_dRdn = (_R * np.sin(_zenith)) / ((_R + _H) * np.sqrt(1 - ((n**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2)))
|
||||||
return np.degrees(_R)
|
|
||||||
|
|
||||||
|
_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)
|
||||||
|
|
||||||
def oriani(self, n1, zenith, reduced_height=None):
|
return np.degrees(_dR)
|
||||||
'''
|
|
||||||
Classic refraction formula with two tan(z) terms. This does not
|
|
||||||
assume any information about the structure of the atmosphere.
|
|
||||||
Found with Oriani's theorem.
|
|
||||||
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
n1 : float
|
|
||||||
Refractive index at the wavelength of interest
|
|
||||||
zenith : float
|
|
||||||
Angle of observation in degrees.
|
|
||||||
zenith = 0 means the observation is directly overhead.
|
|
||||||
reduced_height : float (optional)
|
|
||||||
Reduced height of the atmosphere. Default value is calculated
|
|
||||||
from the object attributes.
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
R : float
|
|
||||||
Atmospheric refraction in degrees
|
|
||||||
'''
|
|
||||||
if reduced_height:
|
|
||||||
H = reduced_height
|
|
||||||
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):
|
def corbard(self, n1, zenith, reduced_height=None):
|
||||||
@@ -246,8 +185,231 @@ class refraction(Observatory):
|
|||||||
return np.degrees(_R)
|
return np.degrees(_R)
|
||||||
|
|
||||||
|
|
||||||
|
def gammaCoefficients(self, n0, reduced_height=None):
|
||||||
|
'''
|
||||||
|
Submethod for matharExponential(). Calculates the gamma
|
||||||
|
coefficients.
|
||||||
|
|
||||||
def refractionIntegral(self, n1, z0, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None):
|
Parameters
|
||||||
|
----------
|
||||||
|
n0 : float
|
||||||
|
Refractive index at the wavelength of interest at the observer
|
||||||
|
reduced_height : float (optional)
|
||||||
|
Reduced height of the atmosphere. Default value is calculated
|
||||||
|
from the object attributes.
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
g : array
|
||||||
|
List of values, corresponding to the gamma coefficients.
|
||||||
|
g[0] corresponds to g_1, g[1] to g_3, g[2] to g_5, and so on.
|
||||||
|
'''
|
||||||
|
|
||||||
|
if reduced_height:
|
||||||
|
H = reduced_height
|
||||||
|
else:
|
||||||
|
H = self.H
|
||||||
|
|
||||||
|
_chi0 = n0**2 - 1
|
||||||
|
_Khat = H / self.rc
|
||||||
|
|
||||||
|
_g1 = (1 - 5/24 * _Khat * _chi0**2 - 3/4 * _chi0 - 3/8 * _Khat**2 * _chi0 \
|
||||||
|
- _Khat - 6 * _Khat**3 + 5/8 * _chi0**2 + 63/128 * _chi0**4 + 3/8 \
|
||||||
|
* _Khat * _chi0 + 35/256 * _Khat * _chi0**3 + 2 * _Khat**2 + 5/36 \
|
||||||
|
* _Khat**2 * _chi0**2 - 35/64 * _chi0**3 + 9/16 * _Khat**3 * _chi0 \
|
||||||
|
- 120 * _Khat**5 + 24 * _Khat**4 - 231 / 512 * _chi0**5 - 63 / 640 \
|
||||||
|
* _Khat * _chi0**4 - 35 / 512 * _Khat**2 * _chi0**3 - 5 / 36 * _Khat**3\
|
||||||
|
* _chi0**2 - 9 / 8 * _Khat**4 * _chi0 )
|
||||||
|
|
||||||
|
_g3 = (35/96 * _chi0**3 - 5/12 * _chi0**2 + 333 * _Khat**4 * _chi0 - 2400 \
|
||||||
|
* _Khat**5 + 3605/4608 * _Khat**2 * _chi0**3 + 63/256 * _Khat * _chi0**4\
|
||||||
|
- 819/16 * _Khat**3 * _chi0 - 35/18 * _Khat**2 * _chi0**2 - 35/96 \
|
||||||
|
* _Khat * _chi0**3 + 69/8 * _Khat**2 * _chi0 + 5/8 * _Khat * _chi0**2 \
|
||||||
|
- 3/2 * _Khat * _chi0 + 1355/216 * _Khat**3 * _chi0**2 + 10 * _Khat**2 \
|
||||||
|
+ 336 * _Khat**4 - 2 * _Khat + _chi0/2 - 21/64 * _chi0**4 - 54 * _Khat**3 \
|
||||||
|
+ 77/256 * _chi0**5 ) * 1/2
|
||||||
|
|
||||||
|
_g5 = (-96 * _Khat**3 + 8 * _Khat**2 - 7/24 * _chi0**3 + 65/4 * _Khat**2 \
|
||||||
|
* _chi0**2 + 35/48 * _Khat * _chi0**3 + 30 * _Khat**2 * _chi0 - 5/3 \
|
||||||
|
* _Khat * _chi0**2 - 3 * _Khat * _chi0 + 20091/8 * _Khat**4 * _chi0 \
|
||||||
|
- 32605/216 * _Khat**3 * _chi0**2 - 35/9 * _Khat**2 * _chi0**3 + 984 \
|
||||||
|
* _Khat**4 + _chi0**2/3 - 10200 * _Khat**5 - 77/320 * _chi0**5 + 21/80 \
|
||||||
|
* _chi0**4 - 1089/4 * _Khat**3 * _chi0 - 7/16 * _Khat * _chi0**4 ) * 3/8
|
||||||
|
|
||||||
|
_g7 = (-17520 * _Khat**5 - 9/40 * _chi0**4 - 8755/12 * _Khat**3 * _chi0**2 \
|
||||||
|
+ 1056 * _Khat**4 + _chi0**3/4 + 23517/4 * _Khat**4 * _chi0 - 11/3 \
|
||||||
|
* _Khat * _chi0**2 + 21 * _Khat**2 * _chi0 - 7/4 * _Khat * _chi0**3 \
|
||||||
|
+ 505/9 * _Khat**2 * _chi0**2 - 801/2 * _Khat**3 * _chi0 + 33/160 \
|
||||||
|
* _chi0**5 - 48 * _Khat**3 + 63/80 * _Khat * _chi0**4 + 2345/96 \
|
||||||
|
* _Khat**2 * _chi0**3) * 5/16
|
||||||
|
|
||||||
|
_g9 = (_chi0**4 / 5 + 384 * _Khat**4 - 13440 * _Khat**5 + 5562 * _Khat**4 \
|
||||||
|
* _chi0 + 3125/36 * _Khat**2 * _chi0**3 - 11/60 * _chi0**5 - 25/6 \
|
||||||
|
* _Khat * _chi0**3 - 26950/27 * _Khat**3 * _chi0**2 + 340/9 * _Khat**2 \
|
||||||
|
* _chi0**2 - 180 * _Khat**3 * _chi0 - 9/5 * _Khat * _chi0**4) * 35/128
|
||||||
|
|
||||||
|
g = np.array([_g1, _g3, _g5, _g7, _g9])
|
||||||
|
g*= n0*_chi0/2
|
||||||
|
|
||||||
|
return g
|
||||||
|
|
||||||
|
|
||||||
|
def H_isotherm(self, T):
|
||||||
|
'''
|
||||||
|
Calculates the reduced height of the atmosphere,
|
||||||
|
assuming an isotherm distribution.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
T : float
|
||||||
|
Atmospheric temperature in Kelvin
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
H : float
|
||||||
|
Reduced height of the atmosphere in meters
|
||||||
|
'''
|
||||||
|
_kb = 1.38064852e-23 # kg m2 s-2 K-1
|
||||||
|
_m = 4.809651698e-26 # kg (weight per molecule of air, assuming dry air)
|
||||||
|
|
||||||
|
_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 matharExponential(self, n1, zenith, reduced_height=None):
|
||||||
|
'''
|
||||||
|
Model on atmospheric refraction described by Richard Mathar in
|
||||||
|
"A Barometric Exponential Model of the Atmosphere’s Refractive Index: ZenithAngles and Second Order Aberration in the Entrance Pupil" (2020)
|
||||||
|
ArXiv id: 2004.11808 (https://arxiv.org/abs/2004.11808)
|
||||||
|
|
||||||
|
The code in this function is adapted from eqs 3, 16, 17 and 18 in this
|
||||||
|
work. It is an adaptation of the RefrExpo.cxx C++ file written by Mathar.
|
||||||
|
|
||||||
|
This model assumes the refractive index, n, is related to the susceptibility,
|
||||||
|
chi, following: n(r) = sqrt(1+chi(r)) (eq. 9).
|
||||||
|
|
||||||
|
R. Mathar then assumes that the susceptibility is proportional to the
|
||||||
|
air density and reduced height, following eq. 10:
|
||||||
|
chi(r) = chi0 * exp(-(r-rho)/K).
|
||||||
|
|
||||||
|
The refraction is described using a Taylor expansion; R in terms of
|
||||||
|
the odd powers of the tangent of the observed zenith angle (eq. 3), not
|
||||||
|
unlike Oriani's theorem. This requires coefficients g_i.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
n1 : float
|
||||||
|
Refractive index at the wavelength of interest at the observer
|
||||||
|
zenith : float
|
||||||
|
Angle of observation in degrees.
|
||||||
|
zenith = 0 means the observation is directly overhead.
|
||||||
|
reduced_height : float (optional)
|
||||||
|
Reduced height of the atmosphere. Default value is calculated
|
||||||
|
from the object attributes.
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
R : float
|
||||||
|
Atmospheric refraction in degrees
|
||||||
|
'''
|
||||||
|
|
||||||
|
if reduced_height:
|
||||||
|
_H = reduced_height
|
||||||
|
else:
|
||||||
|
_H = self.H
|
||||||
|
|
||||||
|
_zenith = np.deg2rad(zenith)
|
||||||
|
_rc = self.rc # rho
|
||||||
|
|
||||||
|
# Atmospheric refraction in radians
|
||||||
|
_gamma = self.gammaCoefficients(n1, reduced_height=_H)
|
||||||
|
_R = sum([_gamma[i]*np.tan(_zenith)**(2*i+1) for i in range(len(_gamma))])
|
||||||
|
|
||||||
|
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.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
n1 : float
|
||||||
|
Refractive index at the wavelength of interest
|
||||||
|
zenith : float
|
||||||
|
Angle of observation in degrees.
|
||||||
|
zenith = 0 means the observation is directly overhead.
|
||||||
|
reduced_height : float (optional)
|
||||||
|
Reduced height of the atmosphere. Default value is calculated
|
||||||
|
from the object attributes.
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
R : float
|
||||||
|
Atmospheric refraction in degrees
|
||||||
|
'''
|
||||||
|
if reduced_height:
|
||||||
|
H = reduced_height
|
||||||
|
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 planeParallel(self, n1, zenith):
|
||||||
|
'''
|
||||||
|
Refraction based on a flat atmosphere, easily derived from Snell's law.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
n1 : float
|
||||||
|
Refractive index at the wavelength of interest
|
||||||
|
zenith : float
|
||||||
|
Angle of observation in degrees.
|
||||||
|
zenith = 0 means the observation is directly overhead.
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
R : float
|
||||||
|
Atmospheric refraction in degrees
|
||||||
|
'''
|
||||||
|
_R = np.arcsin(n * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
|
||||||
|
return np.degrees(_R)
|
||||||
|
|
||||||
|
|
||||||
|
def psi(self, x):
|
||||||
|
'''
|
||||||
|
Submethod for corbard(). Described explicitly in Corbard et al.,2019.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
x : float
|
||||||
|
Input of the function
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
y : float
|
||||||
|
Output of the function
|
||||||
|
'''
|
||||||
|
_f = lambda a: np.exp(-a**2)
|
||||||
|
_int = integrate.quad(_f, x, np.inf)
|
||||||
|
return np.exp(x**2) * _int[0]
|
||||||
|
|
||||||
|
|
||||||
|
def refractionIntegral(self, n1, z0, R0=None, useStandardAtmosphere=True,
|
||||||
|
heightData=None, rhoData=None):
|
||||||
'''
|
'''
|
||||||
Using data that gives temperature, pressure and density as a
|
Using data that gives temperature, pressure and density as a
|
||||||
function of height, the gladstone-dale relation is invoked to
|
function of height, the gladstone-dale relation is invoked to
|
||||||
@@ -295,7 +457,7 @@ class refraction(Observatory):
|
|||||||
_Dens = rhoData
|
_Dens = rhoData
|
||||||
_h = heightData
|
_h = heightData
|
||||||
|
|
||||||
z0 = np.deg2rad(z0)
|
z0 = np.deg2rad(z0)
|
||||||
_Dens = np.asarray(_Dens)
|
_Dens = np.asarray(_Dens)
|
||||||
_h = np.asarray(_h)
|
_h = np.asarray(_h)
|
||||||
_R = R0 + _h * 1000
|
_R = R0 + _h * 1000
|
||||||
@@ -320,46 +482,19 @@ class refraction(Observatory):
|
|||||||
return np.degrees(np.sum(_Integrand))
|
return np.degrees(np.sum(_Integrand))
|
||||||
|
|
||||||
|
|
||||||
|
def setReducedHeight(self, p, rho):
|
||||||
def psi(self, x):
|
|
||||||
'''
|
'''
|
||||||
Submethod for corbard(). Described explicitly in Corbard et al.,2019.
|
Sets the reduced atmospheric height.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
x : float
|
p : float
|
||||||
Input of the function
|
Atmospheric pressure in Pascal
|
||||||
|
rho : float
|
||||||
Returns
|
Density of the atmospheric air at the observer
|
||||||
-------
|
|
||||||
y : float
|
|
||||||
Output of the function
|
|
||||||
'''
|
'''
|
||||||
_f = lambda a: np.exp(-a**2)
|
self.H = self.set_H(p, rho) # Reduced height of the atmosphere assuming ideal gas law
|
||||||
_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.
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
rc : float
|
|
||||||
Local radius of curvature of the Earth in kilometers.
|
|
||||||
'''
|
|
||||||
_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):
|
def set_H(self, p, rho):
|
||||||
'''
|
'''
|
||||||
@@ -387,71 +522,54 @@ class refraction(Observatory):
|
|||||||
_g = _g0_local - (3.0877e-6 - 4.3e-9 * np.sin(_phi)**2) * self.h + 7.2e-13 * self.h**2
|
_g = _g0_local - (3.0877e-6 - 4.3e-9 * np.sin(_phi)**2) * self.h + 7.2e-13 * self.h**2
|
||||||
return p / (rho * _g)
|
return p / (rho * _g)
|
||||||
|
|
||||||
def H_isotherm(self, T):
|
|
||||||
'''
|
|
||||||
Calculates the reduced height of the atmosphere,
|
|
||||||
assuming an isotherm distribution.
|
|
||||||
|
|
||||||
Parameters
|
def set_rc(self):
|
||||||
----------
|
'''
|
||||||
T : float
|
Returns the radius of the curvature at the observer in km,
|
||||||
Atmospheric temperature in Kelvin
|
taking into account the latitude of the observer and the corresponding
|
||||||
|
ellipsoidality of the earth.
|
||||||
|
|
||||||
Returns
|
Returns
|
||||||
-------
|
-------
|
||||||
H : float
|
rc : float
|
||||||
Reduced height of the atmosphere in meters
|
Local radius of curvature of the Earth in kilometers.
|
||||||
'''
|
'''
|
||||||
_kb = 1.38064852e-23 # kg m2 s-2 K-1
|
_A = 6378.137 # km
|
||||||
_m = 4.809651698e-26 # kg (weight per molecule of air, assuming dry air)
|
_B = 6356.752 # km
|
||||||
|
|
||||||
_phi = np.deg2rad(self.lat)
|
_phi = np.deg2rad(self.lat)
|
||||||
_c1 = 5.2790414e-3
|
|
||||||
_c2 = 2.32718e-5
|
_rc0 = (_A*_B)**2 / (_A**2 * np.cos(_phi)**2 + _B**2 * np.sin(_phi)**2)**(3/2)
|
||||||
_c3 = 1.262e-7
|
return _rc0 + self.h/1000
|
||||||
_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):
|
def tan5(self, n1, zenith, reduced_height=None):
|
||||||
'''
|
'''
|
||||||
Calculates the uncertainty in the refraction angle in degrees.
|
Same as oriani's theorem, but including a higher order term.
|
||||||
See Observatory().n_tph() and Observatory().dn_tph() on how to
|
Found in Corbard et al., 2019.
|
||||||
calculate the refractive index and it's uncertainty.
|
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
n : float
|
n1 : float
|
||||||
Refractive index of air
|
Refractive index at the wavelength of interest
|
||||||
dn : float
|
zenith : float
|
||||||
Uncertainty in the refractive index of air
|
Angle of observation in degrees.
|
||||||
zenith : float
|
zenith = 0 means the observation is directly overhead.
|
||||||
Angle of observation in degrees.
|
reduced_height : float (optional)
|
||||||
zenith = 0 means the observation is directly overhead.
|
Reduced height of the atmosphere. Default value is calculated
|
||||||
dz : float (optional)
|
from the object attributes.
|
||||||
Uncertainty in the zenith angle, in degrees.
|
|
||||||
reduced_height : float (optional)
|
|
||||||
Reduced height of the atmosphere. Default value is calculated
|
|
||||||
from the object attributes.
|
|
||||||
|
|
||||||
Returns
|
Returns
|
||||||
-------
|
-------
|
||||||
dR : float
|
R : float
|
||||||
Uncertainty in the atmospheric refraction in degrees.
|
Atmospheric refraction in degrees
|
||||||
'''
|
'''
|
||||||
if reduced_height:
|
if reduced_height:
|
||||||
_H = reduced_height
|
H = reduced_height
|
||||||
else:
|
else:
|
||||||
_H = self.H
|
H = self.H
|
||||||
_zenith = np.deg2rad(zenith)
|
|
||||||
_dz = np.deg2rad(dz)
|
_a = n1 - 1
|
||||||
_R = self.rc
|
_b = H / self.rc
|
||||||
_dRdn = (_R * np.sin(_zenith)) / ((_R + _H) * np.sqrt(1 - ((n**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2)))
|
_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)
|
||||||
|
|
||||||
_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)
|
|
||||||
@@ -2,7 +2,7 @@
|
|||||||
|
|
||||||
**Note:**
|
**Note:**
|
||||||
It would be appreciated if a reference to the following work, for which this package was originally build, is included whenever this code is used for a publication:
|
It would be appreciated if a reference to the following work, for which this package was originally build, is included whenever this code is used for a publication:
|
||||||
*Quantification of the expected residual dispersion of the MICADO Near-IR imaging instrument*, van den Born & Jellema, 2020, MNRAS. Accepted 2020 June 23. Received 2020 June 23; in original form 2020 January 24.
|
*Quantification of the expected residual dispersion of the MICADO Near-IR imaging instrument*, van den Born & Jellema, 2020, MNRAS, DOI: [10.1093/mnras/staa1870](https://doi.org/10.1093/mnras/staa1870).
|
||||||
|
|
||||||
|
|
||||||
**Author info**: *Joost van den Born*, *born@astron.nl*
|
**Author info**: *Joost van den Born*, *born@astron.nl*
|
||||||
@@ -41,7 +41,7 @@ disp = dispersion(lat, h)
|
|||||||
disp.setReducedHeight(P, rho)
|
disp.setReducedHeight(P, rho)
|
||||||
|
|
||||||
# Calculation of the atmopheric dipsersion
|
# Calculation of the atmopheric dipsersion
|
||||||
atm_dispersion = ref.cassini(n1, z0) * 3.6e6
|
atm_dispersion = disp.cassini(n1, z0) * 3.6e6
|
||||||
print ('The dispersion is %.03f milli arc seconds' %(atm_dispersion))
|
print ('The dispersion is %.03f milli arc seconds' %(atm_dispersion))
|
||||||
```
|
```
|
||||||
|
|
||||||
@@ -95,6 +95,7 @@ The available models in this package:
|
|||||||
- `oriani()` - Oriani's theorem ( Atan(z) + Btan^3(z) )
|
- `oriani()` - Oriani's theorem ( Atan(z) + Btan^3(z) )
|
||||||
- `tan5()` - Oriani's theorem expanded ( Atan(z) + Btan^3(z) + Ctan^5(z) )
|
- `tan5()` - Oriani's theorem expanded ( Atan(z) + Btan^3(z) + Ctan^5(z) )
|
||||||
- `corbard()` - Error function refraction model
|
- `corbard()` - Error function refraction model
|
||||||
|
- `matharExponential()` - Barometric exponential model by R. Mathar.
|
||||||
|
|
||||||
|
|
||||||
Input parameters for the spherical atmosphere models (all but the first):
|
Input parameters for the spherical atmosphere models (all but the first):
|
||||||
|
|||||||
+7
-2
@@ -9,9 +9,11 @@ class ADC:
|
|||||||
self.p2 = prism('S-FPL51', 15, wedge[1], wedge[2])
|
self.p2 = prism('S-FPL51', 15, wedge[1], wedge[2])
|
||||||
self.p3 = prism('S-FPL51', 15, wedge[2], wedge[1])
|
self.p3 = prism('S-FPL51', 15, wedge[2], wedge[1])
|
||||||
self.p4 = prism('S-LAH71', 15, wedge[1], wedge[0])
|
self.p4 = prism('S-LAH71', 15, wedge[1], wedge[0])
|
||||||
self.theta = 90
|
self.theta = 90
|
||||||
|
self.FR = 0
|
||||||
|
|
||||||
def setOrientation(self, theta):
|
def setOrientation(self, theta):
|
||||||
|
theta = float(theta)
|
||||||
self.p1.setPrismOrientation(theta)
|
self.p1.setPrismOrientation(theta)
|
||||||
self.p2.setPrismOrientation(theta)
|
self.p2.setPrismOrientation(theta)
|
||||||
self.p3.setPrismOrientation(180-1*theta)
|
self.p3.setPrismOrientation(180-1*theta)
|
||||||
@@ -19,10 +21,12 @@ class ADC:
|
|||||||
self.theta = theta
|
self.theta = theta
|
||||||
|
|
||||||
def setFieldRotation(self, FR):
|
def setFieldRotation(self, FR):
|
||||||
|
FR = float(FR)
|
||||||
self.p1.setPrismOrientation(self.theta-FR)
|
self.p1.setPrismOrientation(self.theta-FR)
|
||||||
self.p2.setPrismOrientation(self.theta-FR)
|
self.p2.setPrismOrientation(self.theta-FR)
|
||||||
self.p3.setPrismOrientation(180-1*self.theta-FR)
|
self.p3.setPrismOrientation(180-1*self.theta-FR)
|
||||||
self.p4.setPrismOrientation(180-1*self.theta-FR)
|
self.p4.setPrismOrientation(180-1*self.theta-FR)
|
||||||
|
self.FR = FR
|
||||||
|
|
||||||
def calculateDispersion(self, l1, l2, hx=0, hy=0, refraction=False):
|
def calculateDispersion(self, l1, l2, hx=0, hy=0, refraction=False):
|
||||||
'''
|
'''
|
||||||
@@ -96,6 +100,7 @@ class prism:
|
|||||||
self.tilt1y = 0
|
self.tilt1y = 0
|
||||||
self.tilt2x = tilt2
|
self.tilt2x = tilt2
|
||||||
self.tilt2y = 0
|
self.tilt2y = 0
|
||||||
|
self.theta = 90
|
||||||
|
|
||||||
# Remember original tilts for reset:
|
# Remember original tilts for reset:
|
||||||
self.tilt1x_OG = tilt1
|
self.tilt1x_OG = tilt1
|
||||||
@@ -116,7 +121,7 @@ class prism:
|
|||||||
def setPrismOrientation(self, theta):
|
def setPrismOrientation(self, theta):
|
||||||
self.tilt1x, self.tilt1y = tiltRotation(self.tilt1x_OG, self.tilt1y_OG, theta)
|
self.tilt1x, self.tilt1y = tiltRotation(self.tilt1x_OG, self.tilt1y_OG, theta)
|
||||||
self.tilt2x, self.tilt2y = tiltRotation(self.tilt2x_OG, self.tilt2y_OG, theta)
|
self.tilt2x, self.tilt2y = tiltRotation(self.tilt2x_OG, self.tilt2y_OG, theta)
|
||||||
|
self.theta = theta
|
||||||
|
|
||||||
|
|
||||||
def tiltRotation(tx,ty, theta):
|
def tiltRotation(tx,ty, theta):
|
||||||
|
|||||||
Reference in new issue
Block a user