"""IAU 1980 Nutation model

Functions
---------
.. autosummary::
    :toctree:
    :template: clean_function.rst

    nutation

Data
----
.. autosummary::
   :toctree:

   coeff80
"""

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"""

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
