/*
Copyright (c) 2022 Peter Hristozov / Петър Христозов

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/

//! ICRF/ITRF reference system as in IAU 2006 resolutions.

pub mod iau06_r {

    use crate::datetime::mjd;
    use crate::math::angle;
    use crate::math::matrix::mat_r;
    use crate::math::scalar::scl;
    use std::f64;

    /**
    # Arguments
    # Returns
    # References*/
    pub fn gam0() -> f64 {
        -0.052928 * angle::ASECRAD
    }

    /**
    # Arguments
    # Returns
    # References*/
    pub fn phi0() -> f64 {
        84381.412819 * angle::ASECRAD
    }

    /**
    # Arguments
    # Returns
    # References*/
    pub fn psi0() -> f64 {
        -0.041775 * angle::ASECRAD
    }

    /**
    # Arguments
    # Returns
    # References*/
    pub fn eps0() -> f64 {
        84381.406 * angle::ASECRAD
    }

    /// Fukushima-Williams angle gamma_GCRS.
    /// # Arguments
    /// * `dtc` TT, Julian centuries since J2000.0
    /// # Returns
    /// GCRS RA difference of ecliptic poles in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn gamb(dtc: f64) -> f64 {
        let coef = [
            -0.052928,
            10.556378,
            0.4932044,
            -0.00031238,
            -0.000002788,
            0.0000000260,
        ];
        scl::polynom(dtc, &coef) * angle::ASECRAD
    }

    /// Fukushima-Williams angle phi_GCRS.
    /// # Arguments
    /// * `dtc` TT, Julian centuries since J2000.0
    /// # Returns
    /// GCRS codeclination of ecliptic pole in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn phib(dtc: f64) -> f64 {
        let coef = [
            84381.412819,
            -46.811016,
            0.0511268,
            0.00053289,
            -0.000000440,
            -0.0000000176,
        ];
        scl::polynom(dtc, &coef) * angle::ASECRAD
    }

    /// Fukushima-Williams angle psi_GCRS.
    /// # Arguments
    /// * `dtc` TT, Julian centuries since J2000.0
    /// # Returns
    /// longitude difference of equator poles, GCRS in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn psib(dtc: f64) -> f64 {
        let coef = [
            -0.041775,
            5038.481484,
            1.5584175,
            -0.00018522,
            -0.000026452,
            -0.0000000148,
        ];
        scl::polynom(dtc, &coef) * angle::ASECRAD
    }

    /// Fukushima-Williams angle gamma_J2000.
    /// # Arguments
    /// * `dtc` TT, Julian centuries since J2000.0
    /// # Returns
    /// J2000.0 RA difference of ecliptic poles in radians
    /// # References
    /// Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351
    pub fn gam(dtc: f64) -> f64 {
        let coef = [
            0.0,
            10.556403,
            0.4932044,
            -0.00031238,
            -0.000002788,
            0.0000000260,
        ];
        scl::polynom(dtc, &coef) * angle::ASECRAD
    }

    /// Fukushima-Williams angle phi_J2000.
    /// # Arguments
    /// * `dtc` TT, Julian centuries since J2000.0
    /// # Returns
    /// J2000.0 codeclination of ecliptic pole in radians
    /// # References
    /// Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351
    pub fn phi(dtc: f64) -> f64 {
        let coef = [
            84381.406,
            -46.811015,
            0.0511269,
            0.00053289,
            -0.000000440,
            -0.0000000176,
        ];
        scl::polynom(dtc, &coef) * angle::ASECRAD
    }

    /// Fukushima-Williams angle psi_J2000.
    /// # Arguments
    /// * `dtc` TT, Julian centuries since J2000.0
    /// # Returns
    /// longitude difference of equator poles, J2000.0 in radians
    /// # References
    /// Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351
    pub fn psi(dtc: f64) -> f64 {
        let coef = [
            0.0,
            5038.481507,
            1.5584176,
            -0.00018522,
            -0.000026452,
            -0.0000000148,
        ];
        scl::polynom(dtc, &coef) * angle::ASECRAD
    }

    /// Mean obliquity of the ecliptic, IAU 2006 precession model.
    /// # Arguments
    /// * `dtc` TT, Julian centuries since J2000.0
    /// # Returns
    /// obliquity of the ecliptic in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn epsa(dtc: f64) -> f64 {
        let coef = [
            84381.406,
            -46.836769,
            -0.0001831,
            0.00200340,
            -0.000000576,
            -0.0000000434,
        ];
        scl::polynom(dtc, &coef) * angle::ASECRAD
    }

    /// Form rotation matrix given the Fukushima-Williams angles.
    /// # Arguments
    /// * `gam` F-W angle gamma (radians)
    /// * `phi` F-W angle phi (radians)
    /// * `psi` F-W angle psi (radians)
    /// * `eps` F-W angle epsilon (radians)
    /// # Returns
    /// The matrix representing the combined effects of
    /// frame bias, precession and nutation is:
    /// # References
    /// Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351
    pub fn fwm(gam: f64, phi: f64, psi: f64, eps: f64) -> mat_r::MatR {
        mat_r::mul(
            &mat_r::rot_x(-eps),
            &mat_r::mul(
                &mat_r::rot_z(-psi),
                &mat_r::mul(&mat_r::rot_x(phi), &mat_r::rot_z(gam)),
            ),
        )
    }

    /// Form frame bias matrix.
    /// # Arguments
    /// # Returns
    /// frame bias matrix that transforms vectors from GCRS
    /// to mean J2000.0 reference system
    /// # References
    /// Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855
    pub fn gcrs_j2000() -> mat_r::MatR {
        fwm(gam0(), phi0(), psi0(), eps0())
    }

    /// Form frame bias and precession matrix.
    /// # Arguments
    /// * `tt-mjd` Terrestrial time as Modified Julian Date
    /// # Returns
    /// bias-precession matrix that transforms vectors from GCRS
    /// to mean of date reference system
    /// # References
    /// Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855
    pub fn gcrs_mean(tt_mjd: f64) -> mat_r::MatR {
        let dtc = (tt_mjd - mjd::J2000) / mjd::DJCENT;
        let gam = gamb(dtc);
        let phi = phib(dtc);
        let psi = psib(dtc);
        let eps = epsa(dtc);
        fwm(gam, phi, psi, eps)
    }

    /// Form precession matrix.
    /// # Arguments
    /// * `tt-mjd` Terrestrial time as Modified Julian Date
    /// # Returns
    /// precession matrix that transforms vectors from J2000
    /// to mean of date reference system
    /// # References
    /// Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855
    pub fn j2000_mean(tt_mjd: f64) -> mat_r::MatR {
        let dtc = (tt_mjd - mjd::J2000) / mjd::DJCENT;
        let gam = gam(dtc);
        let phi = phi(dtc);
        let psi = psi(dtc);
        let eps = epsa(dtc);
        fwm(gam, phi, psi, eps)
    }

    /// Form frame bias and precession matrix.
    /// # Arguments
    /// * `tt-mjd` Terrestrial time as Modified Julian Date
    /// # Returns
    /// bias-precession matrix that transforms vectors from GCRS
    /// to mean of date ecliptic system
    /// # References
    /// Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855
    pub fn gcrs_mean_ecl(tt_mjd: f64) -> mat_r::MatR {
        let dtc = (tt_mjd - mjd::J2000) / mjd::DJCENT;
        let gam = gamb(dtc);
        let phi = phib(dtc);
        let psi = psib(dtc);
        mat_r::mul(
            &mat_r::rot_z(-psi),
            &mat_r::mul(&mat_r::rot_x(phi), &mat_r::rot_z(gam)),
        )
    }

    /// Form frame bias, precession and nutation matrix.
    /// # Arguments
    /// * `tt-mjd` Terrestrial time as Modified Julian Date
    /// * `dpsi` nutation in longitude in radians
    /// * `deps` nutation in obliquity in radians
    /// # Returns
    /// bias-precession-nutation matrix that transforms vectors from GCRS
    /// to true of date reference system
    /// # References
    /// Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855
    pub fn gcrs_true(tt_mjd: f64, dpsi: f64, deps: f64) -> mat_r::MatR {
        let dtc = (tt_mjd - mjd::J2000) / mjd::DJCENT;
        let gam = gamb(dtc);
        let phi = phib(dtc);
        let psi = psib(dtc) + dpsi;
        let eps = epsa(dtc) + deps;
        fwm(gam, phi, psi, eps)
    }

    /// Form precession and nutation matrix.
    /// # Arguments
    /// * `tt-mjd` Terrestrial time as Modified Julian Date
    /// * `dpsi` nutation in longitude in radians
    /// * `deps` nutation in obliquity in radians
    /// # Returns
    /// precession-nutation matrix that transforms vectors from J2000
    /// to true of date reference system
    /// # References
    /// Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855
    pub fn j2000_true(tt_mjd: f64, dpsi: f64, deps: f64) -> mat_r::MatR {
        let dtc = (tt_mjd - mjd::J2000) / mjd::DJCENT;
        let gam = gam(dtc);
        let phi = phi(dtc);
        let psi = psi(dtc) + dpsi;
        let eps = epsa(dtc) + deps;
        fwm(gam, phi, psi, eps)
    }

    /// Form frame bias, precession and nutation matrix.
    /// # Arguments
    /// * `tt-mjd` Terrestrial time as Modified Julian Date
    /// * `dpsi` nutation in longitude in radians
    /// # Returns
    /// bias-precession-nutation matrix that transforms vectors from GCRS
    /// to true of date ecliptic system
    /// # References
    /// Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855
    pub fn gcrs_true_ecl(tt_mjd: f64, dpsi: f64) -> mat_r::MatR {
        let dtc = (tt_mjd - mjd::J2000) / mjd::DJCENT;
        let gam = gamb(dtc);
        let phi = phib(dtc);
        let psi = psib(dtc) + dpsi;
        mat_r::mul(
            &mat_r::rot_z(-psi),
            &mat_r::mul(&mat_r::rot_x(phi), &mat_r::rot_z(gam)),
        )
    }

    /// Extract from the bias-precession-nutation matrix the X,Y coordinates
    /// of the Celestial Intermediate Pole.
    /// # Arguments
    /// * `bpn` celestial-to-true matrix
    /// # Returns
    /// Celestial Intermediate Pole
    /// * `sxy2` s+x*y/2 series where s: CIO locator; x, y: CIP coordinates;
    /// # References
    /// Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351
    pub fn bpn_xy(bpn: &mat_r::MatR) -> (f64, f64) {
        let x = bpn.r20;
        let y = bpn.r21;
        (x, y)
    }

    /// The CIO locator s, positioning the Celestial Intermediate Origin on
    /// the equator of the Celestial Intermediate Pole, given the CIP's X,Y
    /// coordinates.  Compatible with IAU 2006/2000A precession-nutation.
    /// # Arguments
    /// * `tt_mjd` Terrestrial time as Modified Julian Date
    /// * `x`, `y` CIP coordinates
    /// # Returtns
    /// * the CIO locator s in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn s_lo(tt_mjd: f64, x: f64, y: f64) -> f64 {
        let dtc = (tt_mjd - mjd::J2000) / mjd::DJCENT;
        let om = faom(dtc);
        let som = om.sin();
        let s2om = (om + om).sin();
        (94.00e-6 - 2640.73e-6 * som - 63.53e-6 * s2om
            + (3808.65e-6
                + (-122.68e-6
                    + 743.52e-6 * som
                    + (-72574.11e-6 + (27.98e-6 + 15.62e-6 * dtc) * dtc) * dtc)
                    * dtc)
                * dtc)
            * angle::ASECRAD
            - x * y / 2.0
    }

    /// Form the celestial to intermediate-frame-of-date matrix given the CIP
    /// X,Y and the CIO locator s.
    /// # Arguments
    /// * `x` Celestial Intermediate Pole x component
    /// * `y` Celestial Intermediate Pole y component
    /// * `s` the CIO locator in radians
    /// # Returns
    /// celestial-to-intermediate matrix
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn gcrs_cirs(x: f64, y: f64, s: f64) -> mat_r::MatR {
        let r2 = x * x + y * y;
        let e = y.atan2(x);
        let d = (r2 / (1.0 - r2)).sqrt().atan();
        mat_r::mul(
            &mat_r::rot_z(-e - s),
            &mat_r::mul(&mat_r::rot_y(d), &mat_r::rot_z(e)),
        )
    }

    /// Assemble the celestial to terrestrial matrix from CIO-based
    /// components (the celestial-to-intermediate matrix, the Earth Rotation
    /// Angle and the polar motion matrix).
    /// # Arguments
    /// * `gcrs_cirs` celestial-to-intermediate matrix
    /// * `era` Earth rotation angle in radians
    /// * `pom` polar-motion matrix
    /// # Returns
    /// celestial-to-terrestrial matrix
    /// \[TRS\] = pom * Rz(era) * gcrs_cirs * \[CRS\]
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn gcrs_itrs_cio(gcrs_cirs: &mat_r::MatR, era: f64, pom: &mat_r::MatR) -> mat_r::MatR {
        mat_r::mul(&pom, &mat_r::mul(&mat_r::rot_z(era), &gcrs_cirs))
    }

    /// Assemble the celestial to terrestrial matrix from equinox-based
    /// components (the celestial-to-true matrix, the Greenwich Apparent
    /// Sidereal Time and the polar motion matrix).
    /// # Arguments
    /// * `bpn` celestial-to-true matrix
    /// * `gst` Greenwich (apparent) Sidereal Time in radians
    /// * `pom` polar-motion matrix
    /// # Returns
    /// celestial-to-terrestrial matrix
    /// \[TRS\] = pom * Rz(gst) * bpn * \[CRS\]
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn gcrs_itrs_eqx(bpn: &mat_r::MatR, gst: f64, pom: &mat_r::MatR) -> mat_r::MatR {
        mat_r::mul(&pom, &mat_r::mul(&mat_r::rot_z(gst), &bpn))
    }

    /// Earth rotation angle (IAU 2000 model).
    /// # Arguments
    /// * `ut1_mjd` UT1 as a Modified Julian Date
    /// # Returns
    /// Earth rotation angle in radians, range 0-2pi
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn era(ut_mjd: f64) -> f64 {
        let dt = ut_mjd - mjd::J2000;
        let f = 0.5 + ut_mjd % 1.0;
        angle::norm_dblpi(angle::DBLPI * (f + 0.7790572732640 + 0.00273781191135448 * dt))
    }

    /// Equation of the origins, given the classical NPB matrix and the
    /// quantity s.
    /// # Arguments
    /// * `bpn` classical nutation x precession x bias matrix
    /// * `s` the quantity s (the CIO locator) in radians
    /// # Returns
    /// the equation of the origins in radians
    /// # References
    /// Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981
    pub fn eo(bpn: &mat_r::MatR, s: f64) -> f64 {
        let (x, y) = bpn_xy(&bpn);
        let ax = x / (1.0 + bpn.r22);
        let xs = 1.0 - ax * x;
        let ys = -ax * y;
        let zs = -x;
        let p = bpn.r00 * xs + bpn.r01 * ys + bpn.r02 * zs;
        let q = bpn.r10 * xs + bpn.r11 * ys + bpn.r12 * zs;
        if (p != 0.0) || (q != 0.0) {
            s - q.atan2(p)
        } else {
            s
        }
    }

    /// Greenwich mean sidereal time (consistent with IAU 2006 precession).
    /// # Arguments
    /// * `ut1_mjd` UT1 as a Modified Julian Date
    /// * `tt_mjd` TT as a Modified Julian Date
    /// # Returns
    /// Greenwich mean sidereal time in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn gmst(ut1_mjd: f64, tt_mjd: f64) -> f64 {
        let dtc = (tt_mjd - mjd::J2000) / mjd::DJCENT;
        let coef = [
            0.014506,
            4612.156534,
            1.3915817,
            -0.00000044,
            -0.000029956,
            -0.0000000368,
        ];
        angle::norm_dblpi(era(ut1_mjd) + scl::polynom(dtc, &coef) * angle::ASECRAD)
    }

    /// Greenwich apparent sidereal time, IAU 2006, given the NPB matrix.
    /// # Arguments
    /// * `ut1_mjd` UT1 as Modified Julian Date
    /// * `bpn` nutation x precession x bias matrix
    /// * `s` CIO locator in radians
    /// # Returns
    /// Greenwich apparent sidereal time in radians
    /// # References
    /// Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981
    pub fn gst(ut1_mjd: f64, bpn: &mat_r::MatR, s: f64) -> f64 {
        let eo = eo(&bpn, s);
        let era = era(ut1_mjd);
        angle::norm_dblpi(era - eo)
    }

    /// Form the matrix of polar motion for a given date, IAU 2000.
    /// # Arguments
    /// * `tt_mjd`
    /// * `xp`, `yp` coordinates of the pole in radians
    /// # Returns
    /// polar-motion matrix in transformation \[TRS\] = pom * \[CIP\]
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn pom(tt_mjd: f64, xp: f64, yp: f64) -> mat_r::MatR {
        let dtc = (tt_mjd - mjd::J2000) / mjd::DJCENT;
        let sp = -47e-6 * dtc * angle::ASECRAD;
        mat_r::mul(
            &mat_r::rot_x(-yp),
            &mat_r::mul(&mat_r::rot_y(-xp), &&mat_r::rot_z(sp)),
        )
    }

    /// Fundamental argument, mean anomaly of the Moon.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// L mean anomaly of the Moon in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fal(dtc: f64) -> f64 {
        let coef = [
            485868.249036,
            1717915923.2178,
            31.8792,
            0.051635,
            -0.00024470,
        ];
        scl::polynom(dtc, &coef) % angle::TURNAS * angle::ASECRAD
    }

    /// Fundamental argument, mean anomaly of the Sun.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Lp mean anomaly of the Sun in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn falp(dtc: f64) -> f64 {
        let coef = [
            1287104.793048,
            129596581.0481,
            -0.5532,
            0.000136,
            -0.00001149,
        ];
        scl::polynom(dtc, &coef) % angle::TURNAS * angle::ASECRAD
    }

    /// Fundamental argument, mean longitude of the Moon minus mean longitude of the ascending node.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// F mean longitude of the Moon minus mean longitude of the ascending node in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn faf(dtc: f64) -> f64 {
        let coef = [
            335779.526232,
            1739527262.8478,
            -12.7512,
            -0.001037,
            0.00000417,
        ];
        scl::polynom(dtc, &coef) % angle::TURNAS * angle::ASECRAD
    }

    /// Fundamental argument, mean elongation of the Moon from the Sun.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// D mean elongation of the Moon from the Sun in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fad(dtc: f64) -> f64 {
        let coef = [
            1072260.703692,
            1602961601.2090,
            -6.3706,
            0.006593,
            -0.00003169,
        ];
        scl::polynom(dtc, &coef) % angle::TURNAS * angle::ASECRAD
    }

    /// Fundamental argument, mean longitude of the Moon's ascending node.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Omega mean longitude of the Moon's ascending node in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn faom(dtc: f64) -> f64 {
        let coef = [450160.398036, -6962890.5431, 7.4722, 0.007702, -0.00005939];
        scl::polynom(dtc, &coef) % angle::TURNAS * angle::ASECRAD
    }

    /// Fundamental argument, mean longitude of Mercury.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Mercury in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fame(dtc: f64) -> f64 {
        (4.402608842 + 2608.7903141574 * dtc) % angle::DBLPI
    }

    /// Fundamental argument, mean longitude of Venus.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Venus in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fave(dtc: f64) -> f64 {
        (3.176146697 + 1021.3285546211 * dtc) % angle::DBLPI
    }

    /// Fundamental argument, mean longitude of Earth.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Earth in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fae(dtc: f64) -> f64 {
        (1.753470314 + 628.3075849991 * dtc) % angle::DBLPI
    }

    /// Fundamental argument, mean longitude of Mars.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Mars in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fama(dtc: f64) -> f64 {
        (6.203480913 + 334.0612426700 * dtc) % angle::DBLPI
    }

    /// Fundamental argument, mean longitude of Jupiter.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Jupiter in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn faju(dtc: f64) -> f64 {
        (0.599546497 + 52.9690962641 * dtc) % angle::DBLPI
    }

    /// Fundamental argument, mean longitude of Saturn.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Saturn in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fasa(dtc: f64) -> f64 {
        (0.874016757 + 21.3299104960 * dtc) % angle::DBLPI
    }

    /// Fundamental argument, mean longitude of Uranus.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Uranus in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn faur(dtc: f64) -> f64 {
        (5.481293872 + 7.4781598567 * dtc) % angle::DBLPI
    }

    /// Fundamental argument, mean longitude of Neptune.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Neptune in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fane(dtc: f64) -> f64 {
        (5.311886287 + 3.8133035638 * dtc) % angle::DBLPI
    }

    /// Fundamental argument, general accumulated precession in longitude.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// General accumulated precession in longitude in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fapa(dtc: f64) -> f64 {
        (0.024381750 + 0.00000538691 * dtc) * dtc
    }

    /// Nutation, IAU 2000B model.
    /// # Arguments
    /// * `tt_mjd` terrestrial time as Modified Julian Date
    /// # Returns
    /// The nutation components in longitude and obliquity in radians
    /// and with respect to the equinox and ecliptic of date.
    /// # References
    /// McCarthy, D.D. & Luzum, B.J., "An abridged model of the
    /// precession-nutation of the celestial pole", Cel.Mech.Dyn.Astron.
    /// 85, 37-49 (2003)
    pub fn nut00b(tt_mjd: f64) -> (f64, f64) {
        struct Terms(i8, i8, i8, i8, i8, f64, f64, f64, f64, f64, f64);
        let terms_for_nutation = [
            Terms(
                0,
                0,
                0,
                0,
                1,
                -172064161.0,
                -174666.0,
                33386.0,
                92052331.0,
                9086.0,
                15377.0,
            ),
            Terms(
                0,
                0,
                2,
                -2,
                2,
                -13170906.0,
                -1675.0,
                -13696.0,
                5730336.0,
                -3015.0,
                -4587.0,
            ),
            Terms(
                0, 0, 2, 0, 2, -2276413.0, -234.0, 2796.0, 978459.0, -485.0, 1374.0,
            ),
            Terms(
                0, 0, 0, 0, 2, 2074554.0, 207.0, -698.0, -897492.0, 470.0, -291.0,
            ),
            Terms(
                0, 1, 0, 0, 0, 1475877.0, -3633.0, 11817.0, 73871.0, -184.0, -1924.0,
            ),
            Terms(
                0, 1, 2, -2, 2, -516821.0, 1226.0, -524.0, 224386.0, -677.0, -174.0,
            ),
            Terms(1, 0, 0, 0, 0, 711159.0, 73.0, -872.0, -6750.0, 0.0, 358.0),
            Terms(
                0, 0, 2, 0, 1, -387298.0, -367.0, 380.0, 200728.0, 18.0, 318.0,
            ),
            Terms(
                1, 0, 2, 0, 2, -301461.0, -36.0, 816.0, 129025.0, -63.0, 367.0,
            ),
            Terms(
                0, -1, 2, -2, 2, 215829.0, -494.0, 111.0, -95929.0, 299.0, 132.0,
            ),
            Terms(0, 0, 2, -2, 1, 128227.0, 137.0, 181.0, -68982.0, -9.0, 39.0),
            Terms(-1, 0, 2, 0, 2, 123457.0, 11.0, 19.0, -53311.0, 32.0, -4.0),
            Terms(-1, 0, 0, 2, 0, 156994.0, 10.0, -168.0, -1235.0, 0.0, 82.0),
            Terms(1, 0, 0, 0, 1, 63110.0, 63.0, 27.0, -33228.0, 0.0, -9.0),
            Terms(-1, 0, 0, 0, 1, -57976.0, -63.0, -189.0, 31429.0, 0.0, -75.0),
            Terms(-1, 0, 2, 2, 2, -59641.0, -11.0, 149.0, 25543.0, -11.0, 66.0),
            Terms(1, 0, 2, 0, 1, -51613.0, -42.0, 129.0, 26366.0, 0.0, 78.0),
            Terms(-2, 0, 2, 0, 1, 45893.0, 50.0, 31.0, -24236.0, -10.0, 20.0),
            Terms(0, 0, 0, 2, 0, 63384.0, 11.0, -150.0, -1220.0, 0.0, 29.0),
            Terms(0, 0, 2, 2, 2, -38571.0, -1.0, 158.0, 16452.0, -11.0, 68.0),
            Terms(0, -2, 2, -2, 2, 32481.0, 0.0, 0.0, -13870.0, 0.0, 0.0),
            Terms(-2, 0, 0, 2, 0, -47722.0, 0.0, -18.0, 477.0, 0.0, -25.0),
            Terms(2, 0, 2, 0, 2, -31046.0, -1.0, 131.0, 13238.0, -11.0, 59.0),
            Terms(1, 0, 2, -2, 2, 28593.0, 0.0, -1.0, -12338.0, 10.0, -3.0),
            Terms(-1, 0, 2, 0, 1, 20441.0, 21.0, 10.0, -10758.0, 0.0, -3.0),
            Terms(2, 0, 0, 0, 0, 29243.0, 0.0, -74.0, -609.0, 0.0, 13.0),
            Terms(0, 0, 2, 0, 0, 25887.0, 0.0, -66.0, -550.0, 0.0, 11.0),
            Terms(0, 1, 0, 0, 1, -14053.0, -25.0, 79.0, 8551.0, -2.0, -45.0),
            Terms(-1, 0, 0, 2, 1, 15164.0, 10.0, 11.0, -8001.0, 0.0, -1.0),
            Terms(0, 2, 2, -2, 2, -15794.0, 72.0, -16.0, 6850.0, -42.0, -5.0),
            Terms(0, 0, -2, 2, 0, 21783.0, 0.0, 13.0, -167.0, 0.0, 13.0),
            Terms(1, 0, 0, -2, 1, -12873.0, -10.0, -37.0, 6953.0, 0.0, -14.0),
            Terms(0, -1, 0, 0, 1, -12654.0, 11.0, 63.0, 6415.0, 0.0, 26.0),
            Terms(-1, 0, 2, 2, 1, -10204.0, 0.0, 25.0, 5222.0, 0.0, 15.0),
            Terms(0, 2, 0, 0, 0, 16707.0, -85.0, -10.0, 168.0, -1.0, 10.0),
            Terms(1, 0, 2, 2, 2, -7691.0, 0.0, 44.0, 3268.0, 0.0, 19.0),
            Terms(-2, 0, 2, 0, 0, -11024.0, 0.0, -14.0, 104.0, 0.0, 2.0),
            Terms(0, 1, 2, 0, 2, 7566.0, -21.0, -11.0, -3250.0, 0.0, -5.0),
            Terms(0, 0, 2, 2, 1, -6637.0, -11.0, 25.0, 3353.0, 0.0, 14.0),
            Terms(0, -1, 2, 0, 2, -7141.0, 21.0, 8.0, 3070.0, 0.0, 4.0),
            Terms(0, 0, 0, 2, 1, -6302.0, -11.0, 2.0, 3272.0, 0.0, 4.0),
            Terms(1, 0, 2, -2, 1, 5800.0, 10.0, 2.0, -3045.0, 0.0, -1.0),
            Terms(2, 0, 2, -2, 2, 6443.0, 0.0, -7.0, -2768.0, 0.0, -4.0),
            Terms(-2, 0, 0, 2, 1, -5774.0, -11.0, -15.0, 3041.0, 0.0, -5.0),
            Terms(2, 0, 2, 0, 1, -5350.0, 0.0, 21.0, 2695.0, 0.0, 12.0),
            Terms(0, -1, 2, -2, 1, -4752.0, -11.0, -3.0, 2719.0, 0.0, -3.0),
            Terms(0, 0, 0, -2, 1, -4940.0, -11.0, -21.0, 2720.0, 0.0, -9.0),
            Terms(-1, -1, 0, 2, 0, 7350.0, 0.0, -8.0, -51.0, 0.0, 4.0),
            Terms(2, 0, 0, -2, 1, 4065.0, 0.0, 6.0, -2206.0, 0.0, 1.0),
            Terms(1, 0, 0, 2, 0, 6579.0, 0.0, -24.0, -199.0, 0.0, 2.0),
            Terms(0, 1, 2, -2, 1, 3579.0, 0.0, 5.0, -1900.0, 0.0, 1.0),
            Terms(1, -1, 0, 0, 0, 4725.0, 0.0, -6.0, -41.0, 0.0, 3.0),
            Terms(-2, 0, 2, 0, 2, -3075.0, 0.0, -2.0, 1313.0, 0.0, -1.0),
            Terms(3, 0, 2, 0, 2, -2904.0, 0.0, 15.0, 1233.0, 0.0, 7.0),
            Terms(0, -1, 0, 2, 0, 4348.0, 0.0, -10.0, -81.0, 0.0, 2.0),
            Terms(1, -1, 2, 0, 2, -2878.0, 0.0, 8.0, 1232.0, 0.0, 4.0),
            Terms(0, 0, 0, 1, 0, -4230.0, 0.0, 5.0, -20.0, 0.0, -2.0),
            Terms(-1, -1, 2, 2, 2, -2819.0, 0.0, 7.0, 1207.0, 0.0, 3.0),
            Terms(-1, 0, 2, 0, 0, -4056.0, 0.0, 5.0, 40.0, 0.0, -2.0),
            Terms(0, -1, 2, 2, 2, -2647.0, 0.0, 11.0, 1129.0, 0.0, 5.0),
            Terms(-2, 0, 0, 0, 1, -2294.0, 0.0, -10.0, 1266.0, 0.0, -4.0),
            Terms(1, 1, 2, 0, 2, 2481.0, 0.0, -7.0, -1062.0, 0.0, -3.0),
            Terms(2, 0, 0, 0, 1, 2179.0, 0.0, -2.0, -1129.0, 0.0, -2.0),
            Terms(-1, 1, 0, 1, 0, 3276.0, 0.0, 1.0, -9.0, 0.0, 0.0),
            Terms(1, 1, 0, 0, 0, -3389.0, 0.0, 5.0, 35.0, 0.0, -2.0),
            Terms(1, 0, 2, 0, 0, 3339.0, 0.0, -13.0, -107.0, 0.0, 1.0),
            Terms(-1, 0, 2, -2, 1, -1987.0, 0.0, -6.0, 1073.0, 0.0, -2.0),
            Terms(1, 0, 0, 0, 2, -1981.0, 0.0, 0.0, 854.0, 0.0, 0.0),
            Terms(-1, 0, 0, 1, 0, 4026.0, 0.0, -353.0, -553.0, 0.0, -139.0),
            Terms(0, 0, 2, 1, 2, 1660.0, 0.0, -5.0, -710.0, 0.0, -2.0),
            Terms(-1, 0, 2, 4, 2, -1521.0, 0.0, 9.0, 647.0, 0.0, 4.0),
            Terms(-1, 1, 0, 1, 1, 1314.0, 0.0, 0.0, -700.0, 0.0, 0.0),
            Terms(0, -2, 2, -2, 1, -1283.0, 0.0, 0.0, 672.0, 0.0, 0.0),
            Terms(1, 0, 2, 2, 1, -1331.0, 0.0, 8.0, 663.0, 0.0, 4.0),
            Terms(-2, 0, 2, 2, 2, 1383.0, 0.0, -2.0, -594.0, 0.0, -2.0),
            Terms(-1, 0, 0, 0, 2, 1405.0, 0.0, 4.0, -610.0, 0.0, 2.0),
            Terms(1, 1, 2, -2, 2, 1290.0, 0.0, 0.0, -556.0, 0.0, 0.0),
        ];

        let dtc = (tt_mjd - mjd::J2000) / mjd::DJCENT;
        // Full polynomials for fundamental arguments.
        let el = (485868.249036 + 1717915923.2178 * dtc) % angle::TURNAS * angle::ASECRAD;
        let elp = (1287104.79305 + 129596581.0481 * dtc) % angle::TURNAS * angle::ASECRAD;
        let f = (335779.526232 + 1739527262.8478 * dtc) % angle::TURNAS * angle::ASECRAD;
        let d = (1072260.70369 + 1602961601.2090 * dtc) % angle::TURNAS * angle::ASECRAD;
        let om = (450160.398036 - 6962890.5431 * dtc) % angle::TURNAS * angle::ASECRAD;

        let mut dp = 0.0;
        let mut de = 0.0;

        for x in terms_for_nutation.iter() {
            let arg = (x.0 as f64) * el
                + (x.1 as f64) * elp
                + (x.2 as f64) * f
                + (x.3 as f64) * d
                + (x.4 as f64) * om;
            let (s, c) = arg.sin_cos();
            dp += (x.5 + x.6 * dtc) * s + x.7 * c;
            de += (x.8 + x.9 * dtc) * c + x.10 * s;
        }
        let dpsils = dp * 1e-7 * angle::ASECRAD;
        let depsls = de * 1e-7 * angle::ASECRAD;
        let dpsip = -0.135e-3 * angle::ASECRAD;
        let depsp = 0.388e-3 * angle::ASECRAD;

        (dpsils + dpsip, depsls + depsp)
    }
}

pub mod iau06_rv {

    use crate::datetime::mjd;
    use crate::math::angle;
    use crate::math::matrix::mat_rv;
    use crate::math::scalar::scl_sv::{self, SclSv};
    use std::f64;

    /// Fukushima-Williams angle gamma_GCRS.
    /// # Arguments
    /// * `dtc` TT, Julian centuries since J2000.0
    /// # Returns
    /// GCRS RA difference of ecliptic poles in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn gamb(dtc: f64) -> SclSv {
        let coef = [
            -0.052928,
            10.556378,
            0.4932044,
            -0.00031238,
            -0.000002788,
            0.0000000260,
        ];
        scl_sv::scale2(
            angle::ASECRAD,
            angle::ASECRAD / mjd::DJCENT,
            &scl_sv::polynom(dtc, &coef),
        )
    }

    /// Fukushima-Williams angle phi_GCRS.
    /// # Arguments
    /// * `dtc` TT, Julian centuries since J2000.0
    /// # Returns
    /// GCRS codeclination of ecliptic pole in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn phib(dtc: f64) -> SclSv {
        let coef = [
            84381.412819,
            -46.811016,
            0.0511268,
            0.00053289,
            -0.000000440,
            -0.0000000176,
        ];
        scl_sv::scale2(
            angle::ASECRAD,
            angle::ASECRAD / mjd::DJCENT,
            &scl_sv::polynom(dtc, &coef),
        )
    }

    /// Fukushima-Williams angle psi_GCRS.
    /// # Arguments
    /// * `dtc` TT, Julian centuries since J2000.0
    /// # Returns
    /// longitude difference of equator poles, GCRS in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn psib(dtc: f64) -> SclSv {
        let coef = [
            -0.041775,
            5038.481484,
            1.5584175,
            -0.00018522,
            -0.000026452,
            -0.0000000148,
        ];
        scl_sv::scale2(
            angle::ASECRAD,
            angle::ASECRAD / mjd::DJCENT,
            &scl_sv::polynom(dtc, &coef),
        )
    }

    /// Fukushima-Williams angle gamma_J2000.
    /// # Arguments
    /// * `dtc` TT, Julian centuries since J2000.0
    /// # Returns
    /// J2000.0 RA difference of ecliptic poles in radians
    /// # References
    /// Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351
    pub fn gam(dtc: f64) -> SclSv {
        let coef = [
            0.0,
            10.556403,
            0.4932044,
            -0.00031238,
            -0.000002788,
            0.0000000260,
        ];
        scl_sv::scale2(
            angle::ASECRAD,
            angle::ASECRAD / mjd::DJCENT,
            &scl_sv::polynom(dtc, &coef),
        )
    }

    /// Fukushima-Williams angle phi_J2000.
    /// # Arguments
    /// * `dtc` TT, Julian centuries since J2000.0
    /// # Returns
    /// J2000.0 codeclination of ecliptic pole in radians
    /// # References
    /// Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351
    pub fn phi(dtc: f64) -> SclSv {
        let coef = [
            84381.406,
            -46.811015,
            0.0511269,
            0.00053289,
            -0.000000440,
            -0.0000000176,
        ];
        scl_sv::scale2(
            angle::ASECRAD,
            angle::ASECRAD / mjd::DJCENT,
            &scl_sv::polynom(dtc, &coef),
        )
    }

    /// Fukushima-Williams angle psi_J2000.
    /// # Arguments
    /// * `dtc` TT, Julian centuries since J2000.0
    /// # Returns
    /// longitude difference of equator poles, J2000.0 in radians
    /// # References
    /// Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351
    pub fn psi(dtc: f64) -> SclSv {
        let coef = [
            0.0,
            5038.481507,
            1.5584176,
            -0.00018522,
            -0.000026452,
            -0.0000000148,
        ];
        scl_sv::scale2(
            angle::ASECRAD,
            angle::ASECRAD / mjd::DJCENT,
            &scl_sv::polynom(dtc, &coef),
        )
    }

    /// Mean obliquity of the ecliptic, IAU 2006 precession model.
    /// # Arguments
    /// * `dtc` TT, Julian centuries since J2000.0
    /// # Returns
    /// obliquity of the ecliptic in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn epsa(dtc: f64) -> SclSv {
        let coef = [
            84381.406,
            -46.836769,
            -0.0001831,
            0.00200340,
            -0.000000576,
            -0.0000000434,
        ];
        scl_sv::scale2(
            angle::ASECRAD,
            angle::ASECRAD / mjd::DJCENT,
            &scl_sv::polynom(dtc, &coef),
        )
    }

    /// Fundamental argument, mean anomaly of the Moon.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// L mean anomaly of the Moon in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fal(dtc: f64) -> SclSv {
        let coef = [
            485868.249036,
            1717915923.2178,
            31.8792,
            0.051635,
            -0.00024470,
        ];
        scl_sv::scale2(
            angle::ASECRAD,
            angle::ASECRAD / mjd::DJCENT,
            &scl_sv::poly_norm(dtc, &coef, angle::TURNAS),
        )
    }

    /// Fundamental argument, mean anomaly of the Sun.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Lp mean anomaly of the Sun in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn falp(dtc: f64) -> SclSv {
        let coef = [
            1287104.793048,
            129596581.0481,
            -0.5532,
            0.000136,
            -0.00001149,
        ];
        scl_sv::scale2(
            angle::ASECRAD,
            angle::ASECRAD / mjd::DJCENT,
            &scl_sv::poly_norm(dtc, &coef, angle::TURNAS),
        )
    }

    /// Fundamental argument, mean longitude of the Moon minus mean longitude of the ascending node.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// F mean longitude of the Moon minus mean longitude of the ascending node in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn faf(dtc: f64) -> SclSv {
        let coef = [
            335779.526232,
            1739527262.8478,
            -12.7512,
            -0.001037,
            0.00000417,
        ];
        scl_sv::scale2(
            angle::ASECRAD,
            angle::ASECRAD / mjd::DJCENT,
            &scl_sv::poly_norm(dtc, &coef, angle::TURNAS),
        )
    }

    /// Fundamental argument, mean elongation of the Moon from the Sun.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// D mean elongation of the Moon from the Sun in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fad(dtc: f64) -> SclSv {
        let coef = [
            1072260.703692,
            1602961601.2090,
            -6.3706,
            0.006593,
            -0.00003169,
        ];
        scl_sv::scale2(
            angle::ASECRAD,
            angle::ASECRAD / mjd::DJCENT,
            &scl_sv::poly_norm(dtc, &coef, angle::TURNAS),
        )
    }

    /// Fundamental argument, mean longitude of the Moon's ascending node.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Omega mean longitude of the Moon's ascending node in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn faom(dtc: f64) -> SclSv {
        let coef = [450160.398036, -6962890.5431, 7.4722, 0.007702, -0.00005939];
        scl_sv::scale2(
            angle::ASECRAD,
            angle::ASECRAD / mjd::DJCENT,
            &scl_sv::poly_norm(dtc, &coef, angle::TURNAS),
        )
    }

    /// Fundamental argument, mean longitude of Mercury.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Mercury in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fame(dtc: f64) -> SclSv {
        SclSv {
            s: (4.402608842 + 2608.7903141574 * dtc) % angle::DBLPI,
            v: 2608.7903141574 / mjd::DJCENT,
        }
    }

    /// Fundamental argument, mean longitude of Venus.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Venus in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fave(dtc: f64) -> SclSv {
        SclSv {
            s: (3.176146697 + 1021.3285546211 * dtc) % angle::DBLPI,
            v: 1021.3285546211 / mjd::DJCENT,
        }
    }

    /// Fundamental argument, mean longitude of Earth.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Earth in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fae(dtc: f64) -> SclSv {
        SclSv {
            s: (1.753470314 + 628.3075849991 * dtc) % angle::DBLPI,
            v: 628.3075849991 / mjd::DJCENT,
        }
    }

    /// Fundamental argument, mean longitude of Mars.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Mars in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fama(dtc: f64) -> SclSv {
        SclSv {
            s: (6.203480913 + 334.0612426700 * dtc) % angle::DBLPI,
            v: 334.0612426700 / mjd::DJCENT,
        }
    }

    /// Fundamental argument, mean longitude of Jupiter.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Jupiter in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn faju(dtc: f64) -> SclSv {
        SclSv {
            s: (0.599546497 + 52.9690962641 * dtc) % angle::DBLPI,
            v: 52.9690962641 / mjd::DJCENT,
        }
    }

    /// Fundamental argument, mean longitude of Saturn.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Saturn in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fasa(dtc: f64) -> SclSv {
        SclSv {
            s: (0.874016757 + 21.3299104960 * dtc) % angle::DBLPI,
            v: 21.3299104960 / mjd::DJCENT,
        }
    }

    /// Fundamental argument, mean longitude of Uranus.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Uranus in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn faur(dtc: f64) -> SclSv {
        SclSv {
            s: (5.481293872 + 7.4781598567 * dtc) % angle::DBLPI,
            v: 7.4781598567 / mjd::DJCENT,
        }
    }

    /// Fundamental argument, mean longitude of Neptune.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// Mean longitude of Neptune in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fane(dtc: f64) -> SclSv {
        SclSv {
            s: (5.311886287 + 3.8133035638 * dtc) % angle::DBLPI,
            v: 3.8133035638 / mjd::DJCENT,
        }
    }

    /// Fundamental argument, general accumulated precession in longitude.
    /// # Arguments
    /// * `dtc` TDB, Julian centuries since J2000.0
    /// # Returns
    /// General accumulated precession in longitude in radians
    /// # References
    /// Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
    /// IERS Technical Note No. 36, BKG (2010)
    pub fn fapa(dtc: f64) -> SclSv {
        SclSv {
            s: (0.024381750 + 0.00000538691 * dtc) * dtc,
            v: (0.024381750 + 2.0 * 0.00000538691 * dtc) / mjd::DJCENT,
        }
    }

    /// Form rotation matrix given the Fukushima-Williams angles.
    /// # Arguments
    /// * `gam` F-W angle gamma (radians)
    /// * `phi` F-W angle phi (radians)
    /// * `psi` F-W angle psi (radians)
    /// * `eps` F-W angle epsilon (radians)
    /// # Returns
    /// The matrix representing the combined effects of
    /// frame bias, precession and nutation is:
    /// # References
    /// Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351
    pub fn fwm(gam: SclSv, phi: SclSv, psi: SclSv, eps: SclSv) -> mat_rv::MatRv {
        mat_rv::mul(
            &mat_rv::rot_x(&scl_sv::scale(-1.0, &eps)),
            &mat_rv::mul(
                &mat_rv::rot_z(&scl_sv::scale(-1.0, &psi)),
                &mat_rv::mul(&mat_rv::rot_x(&phi), &mat_rv::rot_z(&gam)),
            ),
        )
    }

    /// Form frame bias and precession matrix.
    /// # Arguments
    /// * `tt-mjd` Terrestrial time as Modified Julian Date
    /// # Returns
    /// bias-precession matrix that transforms vectors from GCRS
    /// to mean of date reference system
    /// # References
    /// Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855
    pub fn gcrs_mean(tt_mjd: f64) -> mat_rv::MatRv {
        let dtc = (tt_mjd - mjd::J2000) / mjd::DJCENT;
        let gam = gamb(dtc);
        let phi = phib(dtc);
        let psi = psib(dtc);
        let eps = epsa(dtc);
        fwm(gam, phi, psi, eps)
    }

    /// Form precession matrix.
    /// # Arguments
    /// * `tt-mjd` Terrestrial time as Modified Julian Date
    /// # Returns
    /// precession matrix that transforms vectors from J2000
    /// to mean of date reference system
    /// # References
    /// Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855
    pub fn j2000_mean(tt_mjd: f64) -> mat_rv::MatRv {
        let dtc = (tt_mjd - mjd::J2000) / mjd::DJCENT;
        let gam = gam(dtc);
        let phi = phi(dtc);
        let psi = psi(dtc);
        let eps = epsa(dtc);
        fwm(gam, phi, psi, eps)
    }
}
