sofars 0.6.2

use crate::consts::{DAU, DD2R, DJ00, DJC};
use crate::pnp::pfw06;
use crate::vm::{ir, rx, rxpv, rz, s2pv};

///  Approximate geocentric position and velocity of the Moon.
///
///  This function is part of the International Astronomical Union's
///  SOFA (Standards of Fundamental Astronomy) software collection.
///
///  Status:  support function.
///
///  n.b. Not IAU-endorsed and without canonical status.
///
///  Given:
///     date1  double         TT date part A (Notes 1,4)
///     date2  double         TT date part B (Notes 1,4)
///
///  Returned:
///     pv     [[f64; 3]; 2]  Moon p,v, GCRS (au, au/d, Note 5)
///
pub fn moon98(date1: f64, date2: f64) -> [[f64; 3]; 2] {
    /* Moon's mean longitude (wrt mean equinox and ecliptic of date) */
    const ELP0: f64 = 218.31665436; /* Simon et al. (1994). */
    const ELP1: f64 = 481267.88123421;
    const ELP2: f64 = -0.0015786;
    const ELP3: f64 = 1.0 / 538841.0;
    const ELP4: f64 = -1.0 / 65194000.0;

    /* Moon's mean elongation */
    const D0: f64 = 297.8501921;
    const D1: f64 = 445267.1114034;
    const D2: f64 = -0.0018819;
    const D3: f64 = 1.0 / 545868.0;
    const D4: f64 = 1.0 / 113065000.0;

    /* Sun's mean anomaly */
    const EM0: f64 = 357.5291092;
    const EM1: f64 = 35999.0502909;
    const EM2: f64 = -0.0001536;
    const EM3: f64 = 1.0 / 24490000.0;
    const EM4: f64 = 0.0;

    /* Moon's mean anomaly */
    const EMP0: f64 = 134.9633964;
    const EMP1: f64 = 477198.8675055;
    const EMP2: f64 = 0.0087414;
    const EMP3: f64 = 1.0 / 69699.0;
    const EMP4: f64 = -1.0 / 14712000.0;

    /* Mean distance of the Moon from its ascending node */
    const F0: f64 = 93.2720950;
    const F1: f64 = 483202.0175233;
    const F2: f64 = -0.0036539;
    const F3: f64 = 1.0 / 3526000.0;
    const F4: f64 = 1.0 / 863310000.0;

    /* Meeus A_1, due to Venus (deg) */
    const A10: f64 = 119.75;
    const A11: f64 = 131.849;

    /* Meeus A_2, due to Jupiter (deg) */
    const A20: f64 = 53.09;
    const A21: f64 = 479264.290;

    /* Meeus A_3, due to sidereal motion of the Moon in longitude (deg) */
    const A30: f64 = 313.45;
    const A31: f64 = 481266.484;

    /* Coefficients for Meeus additive terms (deg) */
    const AL1: f64 = 0.003958;
    const AL2: f64 = 0.001962;
    const AL3: f64 = 0.000318;
    const AB1: f64 = -0.002235;
    const AB2: f64 = 0.000382;
    const AB3: f64 = 0.000175;
    const AB4: f64 = 0.000175;
    const AB5: f64 = 0.000127;
    const AB6: f64 = -0.000115;

    /* Fixed term in distance (m) */
    const R0: f64 = 385000560.0;

    /* Coefficients for (dimensionless) E factor */
    const E1: f64 = -0.002516;
    const E2: f64 = -0.0000074;

    struct TermLR {
        nd: i32,
        nem: i32,
        nemp: i32,
        nf: i32,
        coefl: f64,
        coefr: f64,
    }

    const TLR: [TermLR; 60] = [
        TermLR {
            nd: 0,
            nem: 0,
            nemp: 1,
            nf: 0,
            coefl: 6.288774,
            coefr: -20905355.0,
        },
        TermLR {
            nd: 2,
            nem: 0,
            nemp: -1,
            nf: 0,
            coefl: 1.274027,
            coefr: -3699111.0,
        },
        TermLR {
            nd: 2,
            nem: 0,
            nemp: 0,
            nf: 0,
            coefl: 0.658314,
            coefr: -2955968.0,
        },
        TermLR {
            nd: 0,
            nem: 0,
            nemp: 2,
            nf: 0,
            coefl: 0.213618,
            coefr: -569925.0,
        },
        TermLR {
            nd: 0,
            nem: 1,
            nemp: 0,
            nf: 0,
            coefl: -0.185116,
            coefr: 48888.0,
        },
        TermLR {
            nd: 0,
            nem: 0,
            nemp: 0,
            nf: 2,
            coefl: -0.114332,
            coefr: -3149.0,
        },
        TermLR {
            nd: 2,
            nem: 0,
            nemp: -2,
            nf: 0,
            coefl: 0.058793,
            coefr: 246158.0,
        },
        TermLR {
            nd: 2,
            nem: -1,
            nemp: -1,
            nf: 0,
            coefl: 0.057066,
            coefr: -152138.0,
        },
        TermLR {
            nd: 2,
            nem: 0,
            nemp: 1,
            nf: 0,
            coefl: 0.053322,
            coefr: -170733.0,
        },
        TermLR {
            nd: 2,
            nem: -1,
            nemp: 0,
            nf: 0,
            coefl: 0.045758,
            coefr: -204586.0,
        },
        TermLR {
            nd: 0,
            nem: 1,
            nemp: -1,
            nf: 0,
            coefl: -0.040923,
            coefr: -129620.0,
        },
        TermLR {
            nd: 1,
            nem: 0,
            nemp: 0,
            nf: 0,
            coefl: -0.034720,
            coefr: 108743.0,
        },
        TermLR {
            nd: 0,
            nem: 1,
            nemp: 1,
            nf: 0,
            coefl: -0.030383,
            coefr: 104755.0,
        },
        TermLR {
            nd: 2,
            nem: 0,
            nemp: 0,
            nf: -2,
            coefl: 0.015327,
            coefr: 10321.0,
        },
        TermLR {
            nd: 0,
            nem: 0,
            nemp: 1,
            nf: 2,
            coefl: -0.012528,
            coefr: 0.0,
        },
        TermLR {
            nd: 0,
            nem: 0,
            nemp: 1,
            nf: -2,
            coefl: 0.010980,
            coefr: 79661.0,
        },
        TermLR {
            nd: 4,
            nem: 0,
            nemp: -1,
            nf: 0,
            coefl: 0.010675,
            coefr: -34782.0,
        },
        TermLR {
            nd: 0,
            nem: 0,
            nemp: 3,
            nf: 0,
            coefl: 0.010034,
            coefr: -23210.0,
        },
        TermLR {
            nd: 4,
            nem: 0,
            nemp: -2,
            nf: 0,
            coefl: 0.008548,
            coefr: -21636.0,
        },
        TermLR {
            nd: 2,
            nem: 1,
            nemp: -1,
            nf: 0,
            coefl: -0.007888,
            coefr: 24208.0,
        },
        TermLR {
            nd: 2,
            nem: 1,
            nemp: 0,
            nf: 0,
            coefl: -0.006766,
            coefr: 30824.0,
        },
        TermLR {
            nd: 1,
            nem: 0,
            nemp: -1,
            nf: 0,
            coefl: -0.005163,
            coefr: -8379.0,
        },
        TermLR {
            nd: 1,
            nem: 1,
            nemp: 0,
            nf: 0,
            coefl: 0.004987,
            coefr: -16675.0,
        },
        TermLR {
            nd: 2,
            nem: -1,
            nemp: 1,
            nf: 0,
            coefl: 0.004036,
            coefr: -12831.0,
        },
        TermLR {
            nd: 2,
            nem: 0,
            nemp: 2,
            nf: 0,
            coefl: 0.003994,
            coefr: -10445.0,
        },
        TermLR {
            nd: 4,
            nem: 0,
            nemp: 0,
            nf: 0,
            coefl: 0.003861,
            coefr: -11650.0,
        },
        TermLR {
            nd: 2,
            nem: 0,
            nemp: -3,
            nf: 0,
            coefl: 0.003665,
            coefr: 14403.0,
        },
        TermLR {
            nd: 0,
            nem: 1,
            nemp: -2,
            nf: 0,
            coefl: -0.002689,
            coefr: -7003.0,
        },
        TermLR {
            nd: 2,
            nem: 0,
            nemp: -1,
            nf: 2,
            coefl: -0.002602,
            coefr: 0.0,
        },
        TermLR {
            nd: 2,
            nem: -1,
            nemp: -2,
            nf: 0,
            coefl: 0.002390,
            coefr: 10056.0,
        },
        TermLR {
            nd: 1,
            nem: 0,
            nemp: 1,
            nf: 0,
            coefl: -0.002348,
            coefr: 6322.0,
        },
        TermLR {
            nd: 2,
            nem: -2,
            nemp: 0,
            nf: 0,
            coefl: 0.002236,
            coefr: -9884.0,
        },
        TermLR {
            nd: 0,
            nem: 1,
            nemp: 2,
            nf: 0,
            coefl: -0.002120,
            coefr: 5751.0,
        },
        TermLR {
            nd: 0,
            nem: 2,
            nemp: 0,
            nf: 0,
            coefl: -0.002069,
            coefr: 0.0,
        },
        TermLR {
            nd: 2,
            nem: -2,
            nemp: -1,
            nf: 0,
            coefl: 0.002048,
            coefr: -4950.0,
        },
        TermLR {
            nd: 2,
            nem: 0,
            nemp: 1,
            nf: -2,
            coefl: -0.001773,
            coefr: 4130.0,
        },
        TermLR {
            nd: 2,
            nem: 0,
            nemp: 0,
            nf: 2,
            coefl: -0.001595,
            coefr: 0.0,
        },
        TermLR {
            nd: 4,
            nem: -1,
            nemp: -1,
            nf: 0,
            coefl: 0.001215,
            coefr: -3958.0,
        },
        TermLR {
            nd: 0,
            nem: 0,
            nemp: 2,
            nf: 2,
            coefl: -0.001110,
            coefr: 0.0,
        },
        TermLR {
            nd: 3,
            nem: 0,
            nemp: -1,
            nf: 0,
            coefl: -0.000892,
            coefr: 3258.0,
        },
        TermLR {
            nd: 2,
            nem: 1,
            nemp: 1,
            nf: 0,
            coefl: -0.000810,
            coefr: 2616.0,
        },
        TermLR {
            nd: 4,
            nem: -1,
            nemp: -2,
            nf: 0,
            coefl: 0.000759,
            coefr: -1897.0,
        },
        TermLR {
            nd: 0,
            nem: 2,
            nemp: -1,
            nf: 0,
            coefl: -0.000713,
            coefr: -2117.0,
        },
        TermLR {
            nd: 2,
            nem: 2,
            nemp: -1,
            nf: 0,
            coefl: -0.000700,
            coefr: 2354.0,
        },
        TermLR {
            nd: 2,
            nem: 1,
            nemp: -2,
            nf: 0,
            coefl: 0.000691,
            coefr: 0.0,
        },
        TermLR {
            nd: 2,
            nem: -1,
            nemp: 0,
            nf: -2,
            coefl: 0.000596,
            coefr: 0.0,
        },
        TermLR {
            nd: 4,
            nem: 0,
            nemp: 1,
            nf: 0,
            coefl: 0.000549,
            coefr: -1423.0,
        },
        TermLR {
            nd: 0,
            nem: 0,
            nemp: 4,
            nf: 0,
            coefl: 0.000537,
            coefr: -1117.0,
        },
        TermLR {
            nd: 4,
            nem: -1,
            nemp: 0,
            nf: 0,
            coefl: 0.000520,
            coefr: -1571.0,
        },
        TermLR {
            nd: 1,
            nem: 0,
            nemp: -2,
            nf: 0,
            coefl: -0.000487,
            coefr: -1739.0,
        },
        TermLR {
            nd: 2,
            nem: 1,
            nemp: 0,
            nf: -2,
            coefl: -0.000399,
            coefr: 0.0,
        },
        TermLR {
            nd: 0,
            nem: 0,
            nemp: 2,
            nf: -2,
            coefl: -0.000381,
            coefr: -4421.0,
        },
        TermLR {
            nd: 1,
            nem: 1,
            nemp: 1,
            nf: 0,
            coefl: 0.000351,
            coefr: 0.0,
        },
        TermLR {
            nd: 3,
            nem: 0,
            nemp: -2,
            nf: 0,
            coefl: -0.000340,
            coefr: 0.0,
        },
        TermLR {
            nd: 4,
            nem: 0,
            nemp: -3,
            nf: 0,
            coefl: 0.000330,
            coefr: 0.0,
        },
        TermLR {
            nd: 2,
            nem: -1,
            nemp: 2,
            nf: 0,
            coefl: 0.000327,
            coefr: 0.0,
        },
        TermLR {
            nd: 0,
            nem: 2,
            nemp: 1,
            nf: 0,
            coefl: -0.000323,
            coefr: 1165.0,
        },
        TermLR {
            nd: 1,
            nem: 1,
            nemp: -1,
            nf: 0,
            coefl: 0.000299,
            coefr: 0.0,
        },
        TermLR {
            nd: 2,
            nem: 0,
            nemp: 3,
            nf: 0,
            coefl: 0.000294,
            coefr: 0.0,
        },
        TermLR {
            nd: 2,
            nem: 0,
            nemp: -1,
            nf: -2,
            coefl: 0.000000,
            coefr: 8752.0,
        },
    ];

    struct TermB {
        nd: i32,
        nem: i32,
        nemp: i32,
        nf: i32,
        coefb: f64,
    }

    const TB: [TermB; 60] = [
        TermB {
            nd: 0,
            nem: 0,
            nemp: 0,
            nf: 1,
            coefb: 5.128122,
        },
        TermB {
            nd: 0,
            nem: 0,
            nemp: 1,
            nf: 1,
            coefb: 0.280602,
        },
        TermB {
            nd: 0,
            nem: 0,
            nemp: 1,
            nf: -1,
            coefb: 0.277693,
        },
        TermB {
            nd: 2,
            nem: 0,
            nemp: 0,
            nf: -1,
            coefb: 0.173237,
        },
        TermB {
            nd: 2,
            nem: 0,
            nemp: -1,
            nf: 1,
            coefb: 0.055413,
        },
        TermB {
            nd: 2,
            nem: 0,
            nemp: -1,
            nf: -1,
            coefb: 0.046271,
        },
        TermB {
            nd: 2,
            nem: 0,
            nemp: 0,
            nf: 1,
            coefb: 0.032573,
        },
        TermB {
            nd: 0,
            nem: 0,
            nemp: 2,
            nf: 1,
            coefb: 0.017198,
        },
        TermB {
            nd: 2,
            nem: 0,
            nemp: 1,
            nf: -1,
            coefb: 0.009266,
        },
        TermB {
            nd: 0,
            nem: 0,
            nemp: 2,
            nf: -1,
            coefb: 0.008822,
        },
        TermB {
            nd: 2,
            nem: -1,
            nemp: 0,
            nf: -1,
            coefb: 0.008216,
        },
        TermB {
            nd: 2,
            nem: 0,
            nemp: -2,
            nf: -1,
            coefb: 0.004324,
        },
        TermB {
            nd: 2,
            nem: 0,
            nemp: 1,
            nf: 1,
            coefb: 0.004200,
        },
        TermB {
            nd: 2,
            nem: 1,
            nemp: 0,
            nf: -1,
            coefb: -0.003359,
        },
        TermB {
            nd: 2,
            nem: -1,
            nemp: -1,
            nf: 1,
            coefb: 0.002463,
        },
        TermB {
            nd: 2,
            nem: -1,
            nemp: 0,
            nf: 1,
            coefb: 0.002211,
        },
        TermB {
            nd: 2,
            nem: -1,
            nemp: -1,
            nf: -1,
            coefb: 0.002065,
        },
        TermB {
            nd: 0,
            nem: 1,
            nemp: -1,
            nf: -1,
            coefb: -0.001870,
        },
        TermB {
            nd: 4,
            nem: 0,
            nemp: -1,
            nf: -1,
            coefb: 0.001828,
        },
        TermB {
            nd: 0,
            nem: 1,
            nemp: 0,
            nf: 1,
            coefb: -0.001794,
        },
        TermB {
            nd: 0,
            nem: 0,
            nemp: 0,
            nf: 3,
            coefb: -0.001749,
        },
        TermB {
            nd: 0,
            nem: 1,
            nemp: -1,
            nf: 1,
            coefb: -0.001565,
        },
        TermB {
            nd: 1,
            nem: 0,
            nemp: 0,
            nf: 1,
            coefb: -0.001491,
        },
        TermB {
            nd: 0,
            nem: 1,
            nemp: 1,
            nf: 1,
            coefb: -0.001475,
        },
        TermB {
            nd: 0,
            nem: 1,
            nemp: 1,
            nf: -1,
            coefb: -0.001410,
        },
        TermB {
            nd: 0,
            nem: 1,
            nemp: 0,
            nf: -1,
            coefb: -0.001344,
        },
        TermB {
            nd: 1,
            nem: 0,
            nemp: 0,
            nf: -1,
            coefb: -0.001335,
        },
        TermB {
            nd: 0,
            nem: 0,
            nemp: 3,
            nf: 1,
            coefb: 0.001107,
        },
        TermB {
            nd: 4,
            nem: 0,
            nemp: 0,
            nf: -1,
            coefb: 0.001021,
        },
        TermB {
            nd: 4,
            nem: 0,
            nemp: -1,
            nf: 1,
            coefb: 0.000833,
        },
        TermB {
            nd: 0,
            nem: 0,
            nemp: 1,
            nf: -3,
            coefb: 0.000777,
        },
        TermB {
            nd: 4,
            nem: 0,
            nemp: -2,
            nf: 1,
            coefb: 0.000671,
        },
        TermB {
            nd: 2,
            nem: 0,
            nemp: 0,
            nf: -3,
            coefb: 0.000607,
        },
        TermB {
            nd: 2,
            nem: 0,
            nemp: 2,
            nf: -1,
            coefb: 0.000596,
        },
        TermB {
            nd: 2,
            nem: -1,
            nemp: 1,
            nf: -1,
            coefb: 0.000491,
        },
        TermB {
            nd: 2,
            nem: 0,
            nemp: -2,
            nf: 1,
            coefb: -0.000451,
        },
        TermB {
            nd: 0,
            nem: 0,
            nemp: 3,
            nf: -1,
            coefb: 0.000439,
        },
        TermB {
            nd: 2,
            nem: 0,
            nemp: 2,
            nf: 1,
            coefb: 0.000422,
        },
        TermB {
            nd: 2,
            nem: 0,
            nemp: -3,
            nf: -1,
            coefb: 0.000421,
        },
        TermB {
            nd: 2,
            nem: 1,
            nemp: -1,
            nf: 1,
            coefb: -0.000366,
        },
        TermB {
            nd: 2,
            nem: 1,
            nemp: 0,
            nf: 1,
            coefb: -0.000351,
        },
        TermB {
            nd: 4,
            nem: 0,
            nemp: 0,
            nf: 1,
            coefb: 0.000331,
        },
        TermB {
            nd: 2,
            nem: -1,
            nemp: 1,
            nf: 1,
            coefb: 0.000315,
        },
        TermB {
            nd: 2,
            nem: -2,
            nemp: 0,
            nf: -1,
            coefb: 0.000302,
        },
        TermB {
            nd: 0,
            nem: 0,
            nemp: 1,
            nf: 3,
            coefb: -0.000283,
        },
        TermB {
            nd: 2,
            nem: 1,
            nemp: 1,
            nf: -1,
            coefb: -0.000229,
        },
        TermB {
            nd: 1,
            nem: 1,
            nemp: 0,
            nf: -1,
            coefb: 0.000223,
        },
        TermB {
            nd: 1,
            nem: 1,
            nemp: 0,
            nf: 1,
            coefb: 0.000223,
        },
        TermB {
            nd: 0,
            nem: 1,
            nemp: -2,
            nf: -1,
            coefb: -0.000220,
        },
        TermB {
            nd: 2,
            nem: 1,
            nemp: -1,
            nf: -1,
            coefb: -0.000220,
        },
        TermB {
            nd: 1,
            nem: 0,
            nemp: 1,
            nf: 1,
            coefb: -0.000185,
        },
        TermB {
            nd: 2,
            nem: -1,
            nemp: -2,
            nf: -1,
            coefb: 0.000181,
        },
        TermB {
            nd: 0,
            nem: 1,
            nemp: 2,
            nf: 1,
            coefb: -0.000177,
        },
        TermB {
            nd: 4,
            nem: 0,
            nemp: -2,
            nf: -1,
            coefb: 0.000176,
        },
        TermB {
            nd: 4,
            nem: -1,
            nemp: -1,
            nf: -1,
            coefb: 0.000166,
        },
        TermB {
            nd: 1,
            nem: 0,
            nemp: 1,
            nf: -1,
            coefb: -0.000164,
        },
        TermB {
            nd: 4,
            nem: 0,
            nemp: 1,
            nf: -1,
            coefb: 0.000132,
        },
        TermB {
            nd: 1,
            nem: 0,
            nemp: -1,
            nf: -1,
            coefb: -0.000119,
        },
        TermB {
            nd: 4,
            nem: -1,
            nemp: 0,
            nf: -1,
            coefb: 0.000115,
        },
        TermB {
            nd: 2,
            nem: -2,
            nemp: 0,
            nf: 1,
            coefb: 0.000107,
        },
    ];

    let mut rm = [[0.0; 3]; 3];

    /* Time since J2000.0, Julian centuries. */
    let t = ((date1 - DJ00) + date2) / DJC;

    /* Fundamental arguments (Simon et al. 1994). */

    /* Moon's mean longitude. */
    let elp = DD2R * (ELP0 + (ELP1 + (ELP2 + (ELP3 + ELP4 * t) * t) * t) * t).rem_euclid(360.0);
    let delp = DD2R * (ELP1 + (ELP2 * 2.0 + (ELP3 * 3.0 + ELP4 * 4.0 * t) * t) * t);

    /* Moon's mean elongation. */
    let d = DD2R * (D0 + (D1 + (D2 + (D3 + D4 * t) * t) * t) * t).rem_euclid(360.0);
    let dd = DD2R * (D1 + (D2 * 2.0 + (D3 * 3.0 + D4 * 4.0 * t) * t) * t);

    /* Sun's mean anomaly. */
    let em = DD2R * (EM0 + (EM1 + (EM2 + (EM3 + EM4 * t) * t) * t) * t).rem_euclid(360.0);
    let dem = DD2R * (EM1 + (EM2 * 2.0 + (EM3 * 3.0 + EM4 * 4.0 * t) * t) * t);

    /* Moon's mean anomaly. */
    let emp = DD2R * (EMP0 + (EMP1 + (EMP2 + (EMP3 + EMP4 * t) * t) * t) * t).rem_euclid(360.0);
    let demp = DD2R * (EMP1 + (EMP2 * 2.0 + (EMP3 * 3.0 + EMP4 * 4.0 * t) * t) * t);

    /* Mean distance of the Moon from its ascending node. */
    let f = DD2R * (F0 + (F1 + (F2 + (F3 + F4 * t) * t) * t) * t).rem_euclid(360.0);
    let df = DD2R * (F1 + (F2 * 2.0 + (F3 * 3.0 + F4 * 4.0 * t) * t) * t);

    /* Meeus further arguments. */
    let a1 = DD2R * (A10 + A11 * t);
    let da1 = DD2R * A11;
    let a2 = DD2R * (A20 + A21 * t);
    let da2 = DD2R * A21;
    let a3 = DD2R * (A30 + A31 * t);
    let da3 = DD2R * A31;

    /* E-factor, and square. */
    let e = 1.0 + (E1 + E2 * t) * t;
    let de = E1 + 2.0 * E2 * t;
    let esq = e * e;
    let desq = 2.0 * e * de;

    /* Use the Meeus additive terms (deg) to start off the summations. */
    let elpmf = elp - f;
    let delpmf = delp - df;
    let mut vel = AL1 * a1.sin() + AL2 * elpmf.sin() + AL3 * a2.sin();
    let mut vdel = AL1 * a1.cos() * da1 + AL2 * elpmf.cos() * delpmf + AL3 * a2.cos() * da2;

    let mut vr = 0.0;
    let mut vdr = 0.0;

    let a1mf = a1 - f;
    let da1mf = da1 - df;
    let a1pf = a1 + f;
    let da1pf = da1 + df;
    let dlpmp = elp - emp;
    let slpmp = elp + emp;
    let mut vb = AB1 * elp.sin()
        + AB2 * a3.sin()
        + AB3 * a1mf.sin()
        + AB4 * a1pf.sin()
        + AB5 * dlpmp.sin()
        + AB6 * slpmp.sin();
    let mut vdb = AB1 * elp.cos() * delp
        + AB2 * a3.cos() * da3
        + AB3 * a1mf.cos() * da1mf
        + AB4 * a1pf.cos() * da1pf
        + AB5 * dlpmp.cos() * (delp - demp)
        + AB6 * slpmp.cos() * (delp + demp);

    /* ----------------- */
    /* Series expansions */
    /* ----------------- */

    /* Longitude and distance plus derivatives. */
    for n in (0..TLR.len()).rev() {
        let dn = TLR[n].nd as f64;
        let i_val = TLR[n].nem;
        let emn = i_val as f64;
        let empn = TLR[n].nemp as f64;
        let fn_val = TLR[n].nf as f64;
        let en;
        let den;
        match i_val.abs() {
            1 => {
                en = e;
                den = de;
            }
            2 => {
                en = esq;
                den = desq;
            }
            _ => {
                en = 1.0;
                den = 0.0;
            }
        }
        let arg = dn * d + emn * em + empn * emp + fn_val * f;
        let darg = dn * dd + emn * dem + empn * demp + fn_val * df;
        let mut farg = arg.sin();
        let mut v = farg * en;
        let mut dv = arg.cos() * darg * en + farg * den;
        let mut coeff = TLR[n].coefl;
        vel += coeff * v;
        vdel += coeff * dv;
        farg = arg.cos();
        v = farg * en;
        dv = -arg.sin() * darg * en + farg * den;
        coeff = TLR[n].coefr;
        vr += coeff * v;
        vdr += coeff * dv;
    }
    let el = elp + DD2R * vel;
    let del = (delp + DD2R * vdel) / DJC;
    let r = (vr + R0) / DAU;
    let dr = vdr / DAU / DJC;

    /* Latitude plus derivative. */
    for n in (0..TB.len()).rev() {
        let dn = TB[n].nd as f64;
        let i_val = TB[n].nem;
        let emn = i_val as f64;
        let empn = TB[n].nemp as f64;
        let fn_val = TB[n].nf as f64;
        let en;
        let den;
        match i_val.abs() {
            1 => {
                en = e;
                den = de;
            }
            2 => {
                en = esq;
                den = desq;
            }
            _ => {
                en = 1.0;
                den = 0.0;
            }
        }
        let arg = dn * d + emn * em + empn * emp + fn_val * f;
        let darg = dn * dd + emn * dem + empn * demp + fn_val * df;
        let farg = arg.sin();
        let v = farg * en;
        let dv = arg.cos() * darg * en + farg * den;
        let coeff = TB[n].coefb;
        vb += coeff * v;
        vdb += coeff * dv;
    }
    let b = vb * DD2R;
    let db = vdb * DD2R / DJC;

    /* ------------------------------ */
    /* Transformation into final form */
    /* ------------------------------ */

    /* Longitude, latitude to x, y, z (au). */
    let pv = s2pv(el, b, r, del, db, dr);

    /* IAU 2006 Fukushima-Williams bias+precession angles. */
    let (gamb, phib, psib, _) = pfw06(date1, date2);

    /* Mean ecliptic coordinates to GCRS rotation matrix. */
    ir(&mut rm);
    rz(psib, &mut rm);
    rx(-phib, &mut rm);
    rz(-gamb, &mut rm);

    /* Rotate the Moon position and velocity into GCRS (Note 6). */
    let mut res = [[0.0; 3]; 2];
    rxpv(&rm, &pv, &mut res);

    res
}