"""IAU 1980 Nutation model
"""
from math import fmod, fsum
import numpy as np
coeff80 = np.array([[ 0, 0, 0, 0, 1, -171996.0, -174.2, 92025.0, 8.9],
[ 0, 0, 0, 0, 2, 2062.0, 0.2, -895.0, 0.5],
[-2, 0, 2, 0, 1, 46.0, 0.0, -24.0, 0.0],
[ 2, 0, -2, 0, 0, 11.0, 0.0, 0.0, 0.0],
[-2, 0, 2, 0, 2, -3.0, 0.0, 1.0, 0.0],
[ 1, -1, 0, -1, 0, -3.0, 0.0, 0.0, 0.0],
[ 0, -2, 2, -2, 1, -2.0, 0.0, 1.0, 0.0],
[ 2, 0, -2, 0, 1, 1.0, 0.0, 0.0, 0.0],
[ 0, 0, 2, -2, 2, -13187.0, -1.6, 5736.0, -3.1],
[ 0, 1, 0, 0, 0, 1426.0, -3.4, 54.0, -0.1],
[ 0, 1, 2, -2, 2, -517.0, 1.2, 224.0, -0.6],
[ 0, -1, 2, -2, 2, 217.0, -0.5, -95.0, 0.3],
[ 0, 0, 2, -2, 1, 129.0, 0.1, -70.0, 0.0],
[ 2, 0, 0, -2, 0, 48.0, 0.0, 1.0, 0.0],
[ 0, 0, 2, -2, 0, -22.0, 0.0, 0.0, 0.0],
[ 0, 2, 0, 0, 0, 17.0, -0.1, 0.0, 0.0],
[ 0, 1, 0, 0, 1, -15.0, 0.0, 9.0, 0.0],
[ 0, 2, 2, -2, 2, -16.0, 0.1, 7.0, 0.0],
[ 0, -1, 0, 0, 1, -12.0, 0.0, 6.0, 0.0],
[-2, 0, 0, 2, 1, -6.0, 0.0, 3.0, 0.0],
[ 0, -1, 2, -2, 1, -5.0, 0.0, 3.0, 0.0],
[ 2, 0, 0, -2, 1, 4.0, 0.0, -2.0, 0.0],
[ 0, 1, 2, -2, 1, 4.0, 0.0, -2.0, 0.0],
[ 1, 0, 0, -1, 0, -4.0, 0.0, 0.0, 0.0],
[ 2, 1, 0, -2, 0, 1.0, 0.0, 0.0, 0.0],
[ 0, 0, -2, 2, 1, 1.0, 0.0, 0.0, 0.0],
[ 0, 1, -2, 2, 0, -1.0, 0.0, 0.0, 0.0],
[ 0, 1, 0, 0, 2, 1.0, 0.0, 0.0, 0.0],
[-1, 0, 0, 1, 1, 1.0, 0.0, 0.0, 0.0],
[ 0, 1, 2, -2, 0, -1.0, 0.0, 0.0, 0.0],
[ 0, 0, 2, 0, 2, -2274.0, -0.2, 977.0, -0.5],
[ 1, 0, 0, 0, 0, 712.0, 0.1, -7.0, 0.0],
[ 0, 0, 2, 0, 1, -386.0, -0.4, 200.0, 0.0],
[ 1, 0, 2, 0, 2, -301.0, 0.0, 129.0, -0.1],
[ 1, 0, 0, -2, 0, -158.0, 0.0, -1.0, 0.0],
[-1, 0, 2, 0, 2, 123.0, 0.0, -53.0, 0.0],
[ 0, 0, 0, 2, 0, 63.0, 0.0, -2.0, 0.0],
[ 1, 0, 0, 0, 1, 63.0, 0.1, -33.0, 0.0],
[-1, 0, 0, 0, 1, -58.0, -0.1, 32.0, 0.0],
[-1, 0, 2, 2, 2, -59.0, 0.0, 26.0, 0.0],
[ 1, 0, 2, 0, 1, -51.0, 0.0, 27.0, 0.0],
[ 0, 0, 2, 2, 2, -38.0, 0.0, 16.0, 0.0],
[ 2, 0, 0, 0, 0, 29.0, 0.0, -1.0, 0.0],
[ 1, 0, 2, -2, 2, 29.0, 0.0, -12.0, 0.0],
[ 2, 0, 2, 0, 2, -31.0, 0.0, 13.0, 0.0],
[ 0, 0, 2, 0, 0, 26.0, 0.0, -1.0, 0.0],
[-1, 0, 2, 0, 1, 21.0, 0.0, -10.0, 0.0],
[-1, 0, 0, 2, 1, 16.0, 0.0, -8.0, 0.0],
[ 1, 0, 0, -2, 1, -13.0, 0.0, 7.0, 0.0],
[-1, 0, 2, 2, 1, -10.0, 0.0, 5.0, 0.0],
[ 1, 1, 0, -2, 0, -7.0, 0.0, 0.0, 0.0],
[ 0, 1, 2, 0, 2, 7.0, 0.0, -3.0, 0.0],
[ 0, -1, 2, 0, 2, -7.0, 0.0, 3.0, 0.0],
[ 1, 0, 2, 2, 2, -8.0, 0.0, 3.0, 0.0],
[ 1, 0, 0, 2, 0, 6.0, 0.0, 0.0, 0.0],
[ 2, 0, 2, -2, 2, 6.0, 0.0, -3.0, 0.0],
[ 0, 0, 0, 2, 1, -6.0, 0.0, 3.0, 0.0],
[ 0, 0, 2, 2, 1, -7.0, 0.0, 3.0, 0.0],
[ 1, 0, 2, -2, 1, 6.0, 0.0, -3.0, 0.0],
[ 0, 0, 0, -2, 1, -5.0, 0.0, 3.0, 0.0],
[ 1, -1, 0, 0, 0, 5.0, 0.0, 0.0, 0.0],
[ 2, 0, 2, 0, 1, -5.0, 0.0, 3.0, 0.0],
[ 0, 1, 0, -2, 0, -4.0, 0.0, 0.0, 0.0],
[ 1, 0, -2, 0, 0, 4.0, 0.0, 0.0, 0.0],
[ 0, 0, 0, 1, 0, -4.0, 0.0, 0.0, 0.0],
[ 1, 1, 0, 0, 0, -3.0, 0.0, 0.0, 0.0],
[ 1, 0, 2, 0, 0, 3.0, 0.0, 0.0, 0.0],
[ 1, -1, 2, 0, 2, -3.0, 0.0, 1.0, 0.0],
[-1, -1, 2, 2, 2, -3.0, 0.0, 1.0, 0.0],
[-2, 0, 0, 0, 1, -2.0, 0.0, 1.0, 0.0],
[ 3, 0, 2, 0, 2, -3.0, 0.0, 1.0, 0.0],
[ 0, -1, 2, 2, 2, -3.0, 0.0, 1.0, 0.0],
[ 1, 1, 2, 0, 2, 2.0, 0.0, -1.0, 0.0],
[-1, 0, 2, -2, 1, -2.0, 0.0, 1.0, 0.0],
[ 2, 0, 0, 0, 1, 2.0, 0.0, -1.0, 0.0],
[ 1, 0, 0, 0, 2, -2.0, 0.0, 1.0, 0.0],
[ 3, 0, 0, 0, 0, 2.0, 0.0, 0.0, 0.0],
[ 0, 0, 2, 1, 2, 2.0, 0.0, -1.0, 0.0],
[-1, 0, 0, 0, 2, 1.0, 0.0, -1.0, 0.0],
[ 1, 0, 0, -4, 0, -1.0, 0.0, 0.0, 0.0],
[-2, 0, 2, 2, 2, 1.0, 0.0, -1.0, 0.0],
[-1, 0, 2, 4, 2, -2.0, 0.0, 1.0, 0.0],
[ 2, 0, 0, -4, 0, -1.0, 0.0, 0.0, 0.0],
[ 1, 1, 2, -2, 2, 1.0, 0.0, -1.0, 0.0],
[ 1, 0, 2, 2, 1, -1.0, 0.0, 1.0, 0.0],
[-2, 0, 2, 4, 2, -1.0, 0.0, 1.0, 0.0],
[-1, 0, 4, 0, 2, 1.0, 0.0, 0.0, 0.0],
[ 1, -1, 0, -2, 0, 1.0, 0.0, 0.0, 0.0],
[ 2, 0, 2, -2, 1, 1.0, 0.0, -1.0, 0.0],
[ 2, 0, 2, 2, 2, -1.0, 0.0, 0.0, 0.0],
[ 1, 0, 0, 2, 1, -1.0, 0.0, 0.0, 0.0],
[ 0, 0, 4, -2, 2, 1.0, 0.0, 0.0, 0.0],
[ 3, 0, 2, -2, 2, 1.0, 0.0, 0.0, 0.0],
[ 1, 0, 2, -2, 0, -1.0, 0.0, 0.0, 0.0],
[ 0, 1, 2, 0, 1, 1.0, 0.0, 0.0, 0.0],
[-1, -1, 0, 2, 1, 1.0, 0.0, 0.0, 0.0],
[ 0, 0, -2, 0, 1, -1.0, 0.0, 0.0, 0.0],
[ 0, 0, 2, -1, 2, -1.0, 0.0, 0.0, 0.0],
[ 0, 1, 0, 2, 0, -1.0, 0.0, 0.0, 0.0],
[ 1, 0, -2, -2, 0, -1.0, 0.0, 0.0, 0.0],
[ 0, -1, 2, 0, 1, -1.0, 0.0, 0.0, 0.0],
[ 1, 1, 0, -2, 1, -1.0, 0.0, 0.0, 0.0],
[ 1, 0, -2, 2, 0, -1.0, 0.0, 0.0, 0.0],
[ 2, 0, 0, 2, 0, 1.0, 0.0, 0.0, 0.0],
[ 0, 0, 2, 4, 2, -1.0, 0.0, 0.0, 0.0],
[ 0, 1, 0, 1, 0, 1.0, 0.0, 0.0, 0.0]
])
"""IAU 1980 nutation coefficients"""
[docs]
def nutation(TT_JC, const, nTerms=106):
'''Calculate dPsi and dEps for IAU80 nutation model'''
nl = coeff80[:, 0].astype(int)
nlp = coeff80[:, 1].astype(int)
nf = coeff80[:, 2].astype(int)
nd = coeff80[:, 3].astype(int)
nom = coeff80[:, 4].astype(int)
sp = coeff80[:, 5].astype(float)
spt = coeff80[:, 6].astype(float)
ce = coeff80[:, 7].astype(float)
cet = coeff80[:, 8].astype(float)
TT = TT_JC
nTerms = nTerms if nTerms <= 106 else 106
# Mean longitude of Moon (from mean long of Moon's perigee)
Mmoon = fmod((485866.733 + (715922.633 + (31.310 + 0.064 * TT) * TT) * TT)
* const.arcsec + fmod(1325.0 * TT, 1.0) * const.twopi,
const.twopi)
# Mean longitude of Sun (from mean long of Sun's perigee)
Msun = fmod((1287099.804 + (1292581.224 + (-0.577 - 0.012 * TT) * TT) * TT)
* const.arcsec + fmod(99.0 * TT, 1.0) * const.twopi,
const.twopi)
# Offset in mean longitude of Moon from mean lunar node
uMmoon = fmod((335778.877 + (295263.137 + (-13.257 + 0.011 * TT) * TT) * TT)
* const.arcsec + fmod(1342.0 * TT, 1.0) * const.twopi,
const.twopi)
# Mean elongation of Moon
Dsun = fmod((1072261.307 + (1105601.328 + (-6.891 + 0.019 * TT) * TT) * TT)
* const.arcsec + fmod(1236.0 * TT, 1.0) * const.twopi,
const.twopi)
# Longitude of the mean ascending node of Moon's orbit (from the mean equinox of date)
Nmoon = fmod((450160.280 + (-482890.539 + (7.455 + 0.008 * TT) * TT) * TT)
* const.arcsec + fmod(-5.0 * TT, 1.0) * const.twopi,
const.twopi)
psiparts = []
epsparts = []
# Sum nutation terms, in reverse order for best numerical accuracy
for idx in range(nTerms)[::-1]:
arg = (nl[idx] * Mmoon + nlp[idx] * Msun + nf[idx] * uMmoon
+ nd[idx] * Dsun + nom[idx] * Nmoon)
sterm = sp[idx] + spt[idx] * TT
cterm = ce[idx] + cet[idx] * TT
psiparts.append(sterm * np.sin(arg))
epsparts.append(cterm * np.cos(arg))
psiSum = fsum([pp for pp in psiparts if pp != 0])
epsSum = fsum([ep for ep in epsparts if ep != 0])
# Convert to radians (from 1/10 milli-arcseconds)
dPsi = psiSum / 10000
dEps = epsSum / 10000
return dPsi, dEps
