rfa 0.5.9

use crate::{rfam::*, utils::*};
use crate::prec_nut::pwf06::*;
use crate::vector_matrix::{build_rotations::{rx::*, rz::*}, init::ir::*, matrix_vec_prod::*, spherical_cartesian::s2pv::*};
///  Approximate geocentric position and velocity of the Moon.
///
///  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     double[2][3]   Moon p,v, GCRS (AU, AU/d, Note 5)
///
///  Notes:
///
///  1) The TT date date1+date2 is a Julian Date, apportioned in any
///     convenient way between the two arguments.  For example,
///     JD(TT)=2450123.7 could be expressed in any of these ways, among
///     others:
///
///            date1          date2
///
///         2450123.7           0.0       (JD method)
///         2451545.0       -1421.3       (J2000 method)
///         2400000.5       50123.2       (MJD method)
///         2450123.5           0.2       (date & time method)
///
///     The JD method is the most natural and convenient to use in cases
///     where the loss of several decimal digits of resolution is
///     acceptable.  The J2000 method is best matched to the way the
///     argument is handled internally and will deliver the optimum
///     resolution.  The MJD method and the date & time methods are both
///     good compromises between resolution and convenience.  The limited
///     accuracy of the present algorithm is such that any of the methods
///     is satisfactory.
///
///  2) This function is a full implEMentation of the algorithm
///     published by Meeus (see reference) except that the light-time
///     correction to the Moon's mean longitude has been omitted.
///
///  3) Comparisons with ELP/MPP02 over the interval 1950-2100 gave RMS
///     errors of 2.9 arcsec in geocentric direction, 6.1 km in position
///     and 36 mm/s in velocity.  The worst case errors were 18.3 arcsec
///     in geocentric direction, 31.7 km in position and 172 mm/s in
///     velocity.
///
///  4) The original algorithm is expressed in terms of "dynamical time",
///     which can either be TDB or TT without any significant change in
///     accuracy.  UT cannot be used without incurring significant errors
///     (30 arcsec in the present era) due to the Moon's 0.5 arcsec/sec
///     movEMent.
///
///  5) The result is with respect to the GCRS (the same as J2000.0 mean
///     equator and equinox to within 23 mas).
///
///  6) Velocity is obtained by a complete analytical differentiation
///     of the Meeus model.
///
///  7) The Meeus algorithm generates position and velocity in mean
///     ecliptic coordinates of date, which the present function then
///     rotates into GCRS.  Because the ecliptic systEM is precessing,
///     there is a coupling between this spin (about 1.4 degrees per
///     century) and the Moon position that produces a small velocity
///     contribution.  In the present function this effect is neglected
///     as it corresponds to a maximum difference of less than 3 mm/s and
///     increases the RMS error by only 0.4%.
///
///  References:
///
///     Meeus, J., Astronomical Algorithms, 2nd edition, Willmann-Bell,
///     1998, p337.
///
///     Simon, J.L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
///     Francou, G. & Laskar, J., Astron.Astrophys., 1994, 282, 663
///
///  Defined in erfam.h:
///     URSA_DAU           astronomical unit (m)
///     URSA_DJC           days per Julian century
///     URSA_DJ00          reference epoch (J2000.0), Julian Date
///     URSA_DD2R          degrees to radians
///
///  Called:
///     eraS2pv      spherical coordinates to pv-vector
///     eraPfw06     bias-precession F-W angles, IAU 2006
///     eraIr        initialize r-matrix to identity
///     eraRz        rotate around Z-axis
///     eraRx        rotate around X-axis
///     eraRxpv      product of r-matrix and pv-vector
///
///  This revision:  2021 May 11

pub fn moon98(date1: f64, date2: f64, pv: &mut [[f64; 3]; 2])
{
    /*
    ///  Coefficients for fundamental arguments:
    ///
    ///  . Powers of time in Julian centuries
    ///  . Units are degrees.
    */
    
    /* 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;
    
    /*
    /// Other arguments
    */
    
    /* 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;
    
    /*
    /// Coefficients for Moon longitude and distance series
    */
       struct Termlr {
          nd:  i32,           /* multiple of D  in argument           */
          nem: i32,          /*     "    "  M   "    "               */
          nemp: i32,         /*     "    "  M'  "    "               */
          nf: i32,           /*     "    "  F   "    "               */
          coefl: f64,     /* coefficient of L sine argument (deg) */
          coefr: f64,     /* coefficient of R cosine argument (m) */
       }

       impl Termlr {
        pub const fn new( nd:  i32, nem: i32, nemp: i32, nf: i32, coefl: f64,  coefr: f64) ->Self{
            Termlr {nd:nd, nem, nemp, nf: nf, coefl: coefl, coefr: coefr}
        }
      }
    
const TLR:[Termlr;60]=[Termlr::new(0,  0,  1,  0,  6.288774, -20905355.0),
                       Termlr::new(2,  0, -1,  0,  1.274027,  -3699111.0),
                       Termlr::new(2,  0,  0,  0,  0.658314,  -2955968.0),
                       Termlr::new(0,  0,  2,  0,  0.213618,   -569925.0),
                       Termlr::new(0,  1,  0,  0, -0.185116,     48888.0),
                       Termlr::new(0,  0,  0,  2, -0.114332,     -3149.0),
                       Termlr::new(2,  0, -2,  0,  0.058793,    246158.0),
                       Termlr::new(2, -1, -1,  0,  0.057066,   -152138.0),
                       Termlr::new(2,  0,  1,  0,  0.053322,   -170733.0),
                       Termlr::new(2, -1,  0,  0,  0.045758,   -204586.0),
                       Termlr::new(0,  1, -1,  0, -0.040923,   -129620.0),
                       Termlr::new(1,  0,  0,  0, -0.034720,    108743.0),
                       Termlr::new(0,  1,  1,  0, -0.030383,    104755.0),
                       Termlr::new(2,  0,  0, -2,  0.015327,     10321.0),
                       Termlr::new(0,  0,  1,  2, -0.012528,         0.0),
                       Termlr::new(0,  0,  1, -2,  0.010980,     79661.0),
                       Termlr::new(4,  0, -1,  0,  0.010675,    -34782.0),
                       Termlr::new(0,  0,  3,  0,  0.010034,    -23210.0),
                       Termlr::new(4,  0, -2,  0,  0.008548,    -21636.0),
                       Termlr::new(2,  1, -1,  0, -0.007888,     24208.0),
                       Termlr::new(2,  1,  0,  0, -0.006766,     30824.0),
                       Termlr::new(1,  0, -1,  0, -0.005163,     -8379.0),
                       Termlr::new(1,  1,  0,  0,  0.004987,    -16675.0),
                       Termlr::new(2, -1,  1,  0,  0.004036,    -12831.0),
                       Termlr::new(2,  0,  2,  0,  0.003994,    -10445.0),
                       Termlr::new(4,  0,  0,  0,  0.003861,    -11650.0),
                       Termlr::new(2,  0, -3,  0,  0.003665,     14403.0),
                       Termlr::new(0,  1, -2,  0, -0.002689,     -7003.0),
                       Termlr::new(2,  0, -1,  2, -0.002602,         0.0),
                       Termlr::new(2, -1, -2,  0,  0.002390,     10056.0),
                       Termlr::new(1,  0,  1,  0, -0.002348,      6322.0),
                       Termlr::new(2, -2,  0,  0,  0.002236,     -9884.0),
                       Termlr::new(0,  1,  2,  0, -0.002120,      5751.0),
                       Termlr::new(0,  2,  0,  0, -0.002069,         0.0),
                       Termlr::new(2, -2, -1,  0,  0.002048,     -4950.0),
                       Termlr::new(2,  0,  1, -2, -0.001773,      4130.0),
                       Termlr::new(2,  0,  0,  2, -0.001595,         0.0),
                       Termlr::new(4, -1, -1,  0,  0.001215,     -3958.0),
                       Termlr::new(0,  0,  2,  2, -0.001110,         0.0),
                       Termlr::new(3,  0, -1,  0, -0.000892,      3258.0),
                       Termlr::new(2,  1,  1,  0, -0.000810,      2616.0),
                       Termlr::new(4, -1, -2,  0,  0.000759,     -1897.0),
                       Termlr::new(0,  2, -1,  0, -0.000713,     -2117.0),
                       Termlr::new(2,  2, -1,  0, -0.000700,      2354.0),
                       Termlr::new(2,  1, -2,  0,  0.000691,         0.0),
                       Termlr::new(2, -1,  0, -2,  0.000596,         0.0),
                       Termlr::new(4,  0,  1,  0,  0.000549,     -1423.0),
                       Termlr::new(0,  0,  4,  0,  0.000537,     -1117.0),
                       Termlr::new(4, -1,  0,  0,  0.000520,     -1571.0),
                       Termlr::new(1,  0, -2,  0, -0.000487,     -1739.0),
                       Termlr::new(2,  1,  0, -2, -0.000399,         0.0),
                       Termlr::new(0,  0,  2, -2, -0.000381,     -4421.0),
                       Termlr::new(1,  1,  1,  0,  0.000351,         0.0),
                       Termlr::new(3,  0, -2,  0, -0.000340,         0.0),
                       Termlr::new(4,  0, -3,  0,  0.000330,         0.0),
                       Termlr::new(2, -1,  2,  0,  0.000327,         0.0),
                       Termlr::new(0,  2,  1,  0, -0.000323,      1165.0),
                       Termlr::new(1,  1, -1,  0,  0.000299,         0.0),
                       Termlr::new(2,  0,  3,  0,  0.000294,         0.0),
                       Termlr::new(2,  0, -1, -2,  0.000000,      8752.0)];
    
    const NLR: usize = TLR.len();
    
    /*
    /// Coefficients for Moon latitude series
    */
       struct Termb {
          nd:  i32,           /* multiple of D  in argument           */
          nem: i32,          /*     "    "  M   "    "               */
          nemp: i32,         /*     "    "  M'  "    "               */
          nf: i32,           /*     "    "  F   "    "               */
          coefb: f64,        /* coefficient of B sine argument (deg) */
       }

       impl Termb {
           pub const fn new( nd:  i32, nem: i32, nemp: i32, nf: i32, coefb: f64) ->Self{
                Termb{nd:nd, nem, nemp, nf: nf, coefb: coefb}
           }
       }
    
const TB:[Termb;60]=[ Termb::new(0,  0,  0,  1,  5.128122),
                      Termb::new(0,  0,  1,  1,  0.280602),
                      Termb::new(0,  0,  1, -1,  0.277693),
                      Termb::new(2,  0,  0, -1,  0.173237),
                      Termb::new(2,  0, -1,  1,  0.055413),
                      Termb::new(2,  0, -1, -1,  0.046271),
                      Termb::new(2,  0,  0,  1,  0.032573),
                      Termb::new(0,  0,  2,  1,  0.017198),
                      Termb::new(2,  0,  1, -1,  0.009266),
                      Termb::new(0,  0,  2, -1,  0.008822),
                      Termb::new(2, -1,  0, -1,  0.008216),
                      Termb::new(2,  0, -2, -1,  0.004324),
                      Termb::new(2,  0,  1,  1,  0.004200),
                      Termb::new(2,  1,  0, -1, -0.003359),
                      Termb::new(2, -1, -1,  1,  0.002463),
                      Termb::new(2, -1,  0,  1,  0.002211),
                      Termb::new(2, -1, -1, -1,  0.002065),
                      Termb::new(0,  1, -1, -1, -0.001870),
                      Termb::new(4,  0, -1, -1,  0.001828),
                      Termb::new(0,  1,  0,  1, -0.001794),
                      Termb::new(0,  0,  0,  3, -0.001749),
                      Termb::new(0,  1, -1,  1, -0.001565),
                      Termb::new(1,  0,  0,  1, -0.001491),
                      Termb::new(0,  1,  1,  1, -0.001475),
                      Termb::new(0,  1,  1, -1, -0.001410),
                      Termb::new(0,  1,  0, -1, -0.001344),
                      Termb::new(1,  0,  0, -1, -0.001335),
                      Termb::new(0,  0,  3,  1,  0.001107),
                      Termb::new(4,  0,  0, -1,  0.001021),
                      Termb::new(4,  0, -1,  1,  0.000833),
                      Termb::new(0,  0,  1, -3,  0.000777),
                      Termb::new(4,  0, -2,  1,  0.000671),
                      Termb::new(2,  0,  0, -3,  0.000607),
                      Termb::new(2,  0,  2, -1,  0.000596),
                      Termb::new(2, -1,  1, -1,  0.000491),
                      Termb::new(2,  0, -2,  1, -0.000451),
                      Termb::new(0,  0,  3, -1,  0.000439),
                      Termb::new(2,  0,  2,  1,  0.000422),
                      Termb::new(2,  0, -3, -1,  0.000421),
                      Termb::new(2,  1, -1,  1, -0.000366),
                      Termb::new(2,  1,  0,  1, -0.000351),
                      Termb::new(4,  0,  0,  1,  0.000331),
                      Termb::new(2, -1,  1,  1,  0.000315),
                      Termb::new(2, -2,  0, -1,  0.000302),
                      Termb::new(0,  0,  1,  3, -0.000283),
                      Termb::new(2,  1,  1, -1, -0.000229),
                      Termb::new(1,  1,  0, -1,  0.000223),
                      Termb::new(1,  1,  0,  1,  0.000223),
                      Termb::new(0,  1, -2, -1, -0.000220),
                      Termb::new(2,  1, -1, -1, -0.000220),
                      Termb::new(1,  0,  1,  1, -0.000185),
                      Termb::new(2, -1, -2, -1,  0.000181),
                      Termb::new(0,  1,  2,  1, -0.000177),
                      Termb::new(4,  0, -2, -1,  0.000176),
                      Termb::new(4, -1, -1, -1,  0.000166),
                      Termb::new(1,  0,  1, -1, -0.000164),
                      Termb::new(4,  0,  1, -1,  0.000132),
                      Termb::new(1,  0, -1, -1, -0.000119),
                      Termb::new(4, -1,  0, -1,  0.000115),
                      Termb::new(2, -2,  0,  1,  0.000107)];
    
       const NB: usize = TB.len();
    
/* ------------------------------------------------------------------ */
    
    /* Centuries since J2000.0 */
      let t = ((date1 - URSA_DJ00) + date2) / URSA_DJC;
    
    /* --------------------- */
    /* Fundamental arguments */
    /* --------------------- */
    
    /* Arguments (radians) and derivatives (radians per Julian century)
       for the current date. */
    
    /* Moon's mean longitude. */
       let elp = URSA_DD2R * fmod ( ELP0
                       + ( ELP1
                       + ( ELP2
                       + ( ELP3
                       +   ELP4 * t ) * t ) * t ) * t, 360.0 );
       let delp = URSA_DD2R * (     ELP1
                       + ( ELP2 * 2.0
                       + ( ELP3 * 3.0
                       +   ELP4 * 4.0 * t ) * t ) * t );
    
    /* Moon's mean elongation. */
       let d = URSA_DD2R * fmod ( D0
                     + ( D1
                     + ( D2
                     + ( D3
                     +   D4 * t ) * t ) * t ) * t, 360.0 );
       let dd = URSA_DD2R * (     D1
                     + ( D2 * 2.0
                     + ( D3 * 3.0
                     +   D4 * 4.0 * t ) * t ) * t );
    
    /* Sun's mean anomaly. */
       let em = URSA_DD2R * fmod ( EM0
                      + ( EM1
                      + ( EM2
                      + ( EM3
                      +   EM4 * t ) * t ) * t ) * t, 360.0 );
       let dem = URSA_DD2R * (     EM1
                      + ( EM2 * 2.0
                      + ( EM3 * 3.0
                      +   EM4 * 4.0 * t ) * t ) * t );
    
    /* Moon's mean anomaly. */
       let emp = URSA_DD2R * fmod ( EMP0
                       + ( EMP1
                       + ( EMP2
                       + ( EMP3
                       +   EMP4 * t ) * t ) * t ) * t, 360.0 );
       let demp = URSA_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 = URSA_DD2R * fmod ( F0
                     + ( F1
                     + ( F2
                     + ( F3
                     +   F4 * t ) * t ) * t ) * t, 360.0 );
       let df = URSA_DD2R * (     F1
                     + ( F2 * 2.0
                     + ( F3 * 3.0
                     +   F4 * 4.0 * t ) * t ) * t );
    
    /* Meeus further arguments. */
       let a1 = URSA_DD2R * ( A10 + A11*t );
       let da1 = URSA_DD2R * AL1;
       let a2 = URSA_DD2R * ( A20 + A21*t );
       let da2 = URSA_DD2R * A21;
       let a3 = URSA_DD2R * ( A30 + A31*t );
       let da3 = URSA_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 * sin(a1)
           + AL2 * sin(elpmf)
           + AL3 * sin(a2);
       let mut vdel = AL1 * cos(a1) * da1
            + AL2 * cos(elpmf) * delpmf
            + AL3 * cos(a2) * 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 * sin(elp)
          + AB2 * sin(a3)
          + AB3 * sin(a1mf)
          + AB4 * sin(a1pf)
          + AB5 * sin(dlpmp)
          + AB6 * sin(slpmp);
       let mut vdb = AB1 * cos(elp) * delp
           + AB2 * cos(a3) * da3
           + AB3 * cos(a1mf) * da1mf
           + AB4 * cos(a1pf) * da1pf
           + AB5 * cos(dlpmp) * (delp-demp)
           + AB6 * cos(slpmp) * (delp+demp);
    
    /* ----------------- */
    /* Series expansions */
    /* ----------------- */
    
    /* Longitude and distance plus derivatives. */
       for n in 0..NLR {
          let dn = TLR[n].nd as f64;
          let i = TLR[n].nem;
          let emn = TLR[n].nem as f64;
          let empn = TLR[n].nemp as f64;
          let tfn = TLR[n].nf as f64;
          let en: f64; let den: f64;
          match i.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 + tfn*f;
          let darg = dn*dd + emn*dem + empn*demp + tfn*df;
          let farg = sin(arg);
          let v = farg * en;
          let dv = cos(arg)*darg*en + farg*den;
          let coeff = TLR[n].coefl;
          vel += coeff * v;
          vdel += coeff * dv;
          let farg = cos(arg);
          let v = farg * en;
          let dv = -sin(arg)*darg*en + farg*den;
          let coeff = TLR[n].coefr;
          vr += coeff * v;
          vdr += coeff * dv;
       }
       let el = elp + URSA_DD2R*vel;
       let del = ( delp + URSA_DD2R*vdel ) / URSA_DJC;
       let r = ( vr + R0 ) / URSA_DAU;
       let dr = vdr / URSA_DAU / URSA_DJC;
    
    /* Latitude plus derivative. */
       for n in 0..NB {
          let dn = TB[n].nd as f64;
          let i =  TB[n].nem;
          let emn = TB[n].nem as f64;
          let empn = TB[n].nemp as f64;
          let tfn = TB[n].nf as f64;
          let en: f64; let den: f64;
          match i.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 + tfn*f;
          let darg = dn*dd + emn*dem + empn*demp + tfn*df;
          let farg = sin(arg);
          let v = farg * en;
          let dv = cos(arg)*darg*en + farg*den;
          let coeff = TB[n].coefb;
          vb += coeff * v;
          vdb += coeff * dv;
       }
       let b = vb * URSA_DD2R;
       let db = vdb * URSA_DD2R / URSA_DJC;
    
    /* ------------------------------ */
    /* Transformation into final form */
    /* ------------------------------ */
    let mut local_pv = [[0.0; 3]; 2];
    /* Longitude, latitude to x, y, z (AU). */
       s2pv ( el, b, r, del, db, dr, &mut local_pv );
    let mut gamb = 0.0; let mut phib =0.0; let mut psib = 0.0; let mut epsa =0.0;
    /* IAU 2006 Fukushima-Williams bias+precession angles. */
       pfw06 ( date1, date2, &mut gamb, &mut phib, &mut psib, &mut epsa );
    let mut rm = [[0.0; 3];3];
    /* 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). */
       rxpv ( &rm, &local_pv, pv );
    
    /* Finished. */
   
}