// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.
//! Precession and nutation code.
pub use nut00a;
use crate*;
/// Mean obliquity of the ecliptic, IAU 2006 precession model. (`eraObl06`)
///
/// Given:
/// * `date1`,`date2`: TT as a 2-part Julian Date (Note 1)
///
/// Returned:
/// * obliquity of the ecliptic (radians, Note 2)
///
/// # 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.
///
/// 2) The result is the angle between the ecliptic and mean equator of date
/// `date1+date2`.
///
/// # References:
///
/// * Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351
///
/// Precession angles, IAU 2006, equinox based. (`eraP06e`)
///
/// Given:
/// * `date1`,`date2`: TT as a 2-part Julian Date (Note 1)
///
/// Returned (see Note 2):
/// * `eps0`: epsilon_0
/// * `psia`: psi_A
/// * `oma`: omega_A
/// * `bpa`: P_A
/// * `bqa`: Q_A
/// * `pia`: pi_A
/// * `bpia`: Pi_A
/// * `epsa`: obliquity epsilon_A
/// * `chia`: chi_A
/// * `za`: z_A
/// * `zetaa`: zeta_A
/// * `thetaa`: theta_A
/// * `pa`: p_A
/// * `gam`: F-W angle gamma_J2000
/// * `phi`: F-W angle phi_J2000
/// * `psi`: F-W angle psi_J2000
///
/// 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.
///
/// 2) This function returns the set of equinox based angles for the Capitaine
/// et al. "P03" precession theory, adopted by the IAU in
/// 2006. The angles are set out in Table 1 of Hilton et al. (2006):
///
/// * `eps0` `epsilon_0` obliquity at J2000.0
/// * `psia` `psi_A` luni-solar precession
/// * `oma` `omega_A` inclination of equator wrt J2000.0 ecliptic
/// * `bpa` `P_A` ecliptic pole x, J2000.0 ecliptic triad
/// * `bqa` `Q_A` ecliptic pole -y, J2000.0 ecliptic triad
/// * `pia` `pi_A` angle between moving and J2000.0 ecliptics
/// * `bpia` `Pi_A` longitude of ascending node of the ecliptic
/// * `epsa` `epsilon_A` obliquity of the ecliptic
/// * `chia` `chi_A` planetary precession
/// * `za` `z_A` equatorial precession: -3rd 323 Euler angle
/// * `zetaa` `zeta_A` equatorial precession: -1st 323 Euler angle
/// * `thetaa` `theta_A` equatorial precession: 2nd 323 Euler angle
/// * `pa` `p_A` general precession (n.b. see below)
/// * `gam` `gamma_J2000` J2000.0 RA difference of ecliptic poles
/// * `phi` `phi_J2000` J2000.0 codeclination of ecliptic pole
/// * `psi` `psi_J2000` longitude difference of equator poles, J2000.0
///
/// The returned values are all radians.
///
/// Note that the `t^5` coefficient in the series for `p_A` from Capitaine et
/// al. (2003) is incorrectly signed in Hilton et al. (2006).
///
/// 3) Hilton et al. (2006) Table 1 also contains angles that depend on models
/// distinct from the P03 precession theory itself, namely the IAU 2000A
/// frame bias and nutation. The quoted polynomials are used in other ERFA
/// functions:
///
/// * `eraXy06` contains the polynomial parts of the `X` and `Y` series.
///
/// * `eraS06` contains the polynomial part of the `s+XY/2` series.
///
/// * `eraPfw06`: implements the series for the Fukushima-Williams angles
/// that are with respect to the GCRS pole (i.e. the variants that include
/// frame bias).
///
/// 4) The IAU resolution stipulated that the choice of parameterization was
/// left to the user, and so an IAU compliant precession implementation can
/// be constructed using various combinations of the angles returned by the
/// present function.
///
/// 5) The parameterization used by ERFA is the version of the Fukushima-
/// Williams angles that refers directly to the GCRS pole. These angles may
/// be calculated by calling the function
/// [`eraPfw06`](precession_angles_fw06). ERFA also supports the direct
/// computation of the CIP GCRS `X`,`Y` by series, available by calling
/// `eraXy06`.
///
/// 6) The agreement between the different parameterizations is at the 1
/// microarcsecond level in the present era.
///
/// 7) When constructing a precession formulation that refers to the GCRS pole
/// rather than the dynamical pole, it may (depending on the choice of
/// angles) be necessary to introduce the frame bias explicitly.
///
/// 8) It is permissible to re-use the same variable in the returned arguments.
/// The quantities are stored in the stated order.
///
/// References:
///
/// * Capitaine, N., Wallace, P.T. & Chapront, J., 2003, Astron.Astrophys., 412,
/// 567
///
/// * Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351
///
/// Precession angles, IAU 2006 (Fukushima-Williams 4-angle formulation).
/// (`eraPfw06`)
///
/// Given:
/// * `date1`,`date2`: TT as a 2-part Julian Date (Note 1)
///
/// Returned:
/// * `gamb`: F-W angle gamma_bar (radians)
/// * `phib`: F-W angle phi_bar (radians)
/// * `psib`: F-W angle psi_bar (radians)
/// * `epsa`: F-W angle epsilon_A (radians)
///
/// # 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.
///
/// 2) Naming the following points:
///
/// * `e` = J2000.0 ecliptic pole,
/// * `p` = GCRS pole,
/// * `E` = mean ecliptic pole of date,
/// * `P` = mean pole of date,
///
/// the four Fukushima-Williams angles are as follows:
///
/// * `gamb` = gamma_bar = epE
/// * `phib` = phi_bar = pE
/// * `psib` = psi_bar = pEP
/// * `epsa` = epsilon_A = EP
///
/// 3) The matrix representing the combined effects of frame bias and precession
/// is:
///
/// `PxB = R_1(-epsa).R_3(-psib).R_1(phib).R_3(gamb)``
///
/// 4) The matrix representing the combined effects of frame bias, precession
/// and nutation is simply:
///
/// `NxPxB = R_1(-epsa-dE).R_3(-psib-dP).R_1(phib).R_3(gamb)``
///
/// where `dP` and `dE` are the nutation components with respect to the
/// ecliptic of date.
///
/// References:
///
/// * Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351
///
/// IAU 2000A nutation with adjustments to match the IAU 2006 precession.
/// (`eraNut06a`)
///
/// Given:
/// * `date1`,`date2`: TT as a 2-part Julian Date (Note 1)
///
/// Returned:
/// * `dpsi`,`deps`: nutation, luni-solar + planetary (Note 2)
///
/// # 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.
///
/// 2) The nutation components in longitude and obliquity are in radians and
/// with respect to the mean equinox and ecliptic of date, IAU 2006
/// precession model (Hilton et al. 2006, Capitaine et al. 2005).
///
/// 3) The function first computes the IAU 2000A nutation, then applies
/// adjustments for (i) the consequences of the change in obliquity from the
/// IAU 1980 ecliptic to the IAU 2006 ecliptic and (ii) the secular variation
/// in the Earth's dynamical form factor J2.
///
/// 4) The present function provides classical nutation, complementing the IAU
/// 2000 frame bias and IAU 2006 precession. It delivers a pole which is at
/// current epochs accurate to a few tens of microarcseconds, apart from the
/// free core nutation.
///
/// # References:
///
/// * Chapront, J., Chapront-Touze, M. & Francou, G. 2002, Astron.Astrophys.
/// 387, 700
/// * Lieske, J.H., Lederle, T., Fricke, W. & Morando, B. 1977,
/// Astron.Astrophys. 58, 1-16
/// * Mathews, P.M., Herring, T.A., Buffet, B.A. 2002, J.Geophys.Res. 107, B4.
/// The MHB_2000 code itself was obtained on 9th September 2002 from
/// <ftp://maia.usno.navy.mil/conv2000/chapter5/IAU2000A>.
/// * Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M., Francou,
/// G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683
/// * Souchay, J., Loysel, B., Kinoshita, H., Folgueira, M. 1999,
/// Astron.Astrophys.Supp.Ser. 135, 111
/// * Wallace, P.T., "Software for Implementing the IAU 2000 Resolutions", in
/// IERS Workshop 5.1 (2002)
///
/// Form rotation matrix given the Fukushima-Williams angles. (`eraFw2m`)
///
/// Given:
/// * `gamb`: F-W angle gamma_bar (radians)
/// * `phib`: F-W angle phi_bar (radians)
/// * `psi`: F-W angle psi (radians)
/// * `eps`: F-W angle epsilon (radians)
///
/// Returned:
/// * `r`: rotation matrix
///
/// # Notes:
///
/// 1) Naming the following points:
///
/// * `e` = J2000.0 ecliptic pole,
/// * `p` = GCRS pole,
/// * `E` = ecliptic pole of date,
/// * `P` = CIP,
///
/// the four Fukushima-Williams angles are as follows:
///
/// * `gamb = gamma = epE`,
/// * `phib = phi = pE`,
/// * `psi = psi = pEP`,
/// * `eps = epsilon = EP`
///
/// 2) The matrix representing the combined effects of frame bias, precession
/// and nutation is:
///
/// `NxPxB = R_1(-eps).R_3(-psi).R_1(phib).R_3(gamb)`
///
/// 3) The present function can construct three different matrices, depending on
/// which angles are supplied as the arguments `gamb`, `phib`, `psi` and
/// `eps`:
///
/// - To obtain the nutation x precession x frame bias matrix, first
/// generate the four precession angles known conventionally as
/// `gamma_bar`, `phi_bar`, `psi_bar` and `epsilon_A`, then generate the
/// nutation components `Dpsi` and `Depsilon` and add them to `psi_bar`
/// and `epsilon_A`, and finally call the present function using those
/// four angles as arguments.
///
/// - To obtain the precession x frame bias matrix, generate the four
/// precession angles and call the present function.
///
/// - To obtain the frame bias matrix, generate the four precession angles
/// for date J2000.0 and call the present function.
///
/// The nutation-only and precession-only matrices can if necessary be
/// obtained by combining these three appropriately.
///
/// # References:
///
/// * Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855
/// * Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351
///
/// Form the matrix of precession-nutation for a given date (including frame
/// bias), equinox based, IAU 2006 precession and IAU 2000A nutation models.
/// (`eraPnm06a`)
///
/// Given:
/// * `date1`,`date2`: TT as a 2-part Julian Date (Note 1)
///
/// Returned:
/// * `rbpn`: bias-precession-nutation matrix (Note 2)
///
/// # 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.
///
/// 2) The matrix operates in the sense `V(date) = rbpn * V(GCRS)`, where the
/// p-vector `V(date)` is with respect to the true equatorial triad of date
/// `date1+date2` and the p-vector `V(GCRS)` is with respect to the
/// Geocentric Celestial Reference System (IAU, 2000).
///
/// # References:
///
/// * Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855.
///
/// Extract from the bias-precession-nutation matrix the X,Y coordinates of the
/// Celestial Intermediate Pole. (`eraBpn2xy`)
///
/// Given:
/// * `rbpn`: celestial-to-true matrix (Note 1)
///
/// Returned:
/// * `x`,`y`: Celestial Intermediate Pole (Note 2)
///
/// # Notes:
///
/// 1) The matrix `rbpn` transforms vectors from GCRS to true equator (and CIO
/// or equinox) of date, and therefore the Celestial Intermediate Pole unit
/// vector is the bottom row of the matrix.
///
/// 2) `x` and `y` are components of the Celestial Intermediate Pole unit vector
/// in the Geocentric Celestial Reference System.
///
/// # Reference:
///
/// * "Expressions for the Celestial Intermediate Pole and Celestial Ephemeris
/// Origin consistent with the IAU 2000A precession- nutation model",
/// Astron.Astrophys. 400, 1145-1154 (2003)
///
/// n.b. The celestial ephemeris origin (CEO) was renamed "celestial
/// intermediate origin" (CIO) by IAU 2006 Resolution 2.
///
/// Equation of the origins, given the classical NPB matrix and the quantity
/// `s`. (`eraEors`)
///
/// Given:
/// * `rnpb` classical nutation x precession x bias matrix
/// * `s`: the quantity s (the CIO locator) in radians
///
/// Returned:
/// * the equation of the origins in radians
///
/// # Notes:
///
/// 1) The equation of the origins is the distance between the true equinox and
/// the celestial intermediate origin and, equivalently, the difference
/// between Earth rotation angle and Greenwich apparent sidereal time
/// (ERA-GST). It comprises the precession (since J2000.0) in right
/// ascension plus the equation of the equinoxes (including the small
/// correction terms).
///
/// 2) The algorithm is from Wallace & Capitaine (2006).
///
/// # References:
///
/// * Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855
/// * Wallace, P. & Capitaine, N., 2006, Astron.Astrophys. 459, 981
///
/// Precession matrix (including frame bias) from GCRS to a specified date, IAU
/// 2006 model. (`eraPmat06`)
///
/// Given:
/// * `date1`,`date2`: TT as a 2-part Julian Date (Note 1)
///
/// Returned:
/// * `rbp`: bias-precession matrix (Note 2)
///
/// # 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.
///
/// 2) The matrix operates in the sense `V(date) = rbp * V(GCRS)`, where the
/// p-vector `V(GCRS)` is with respect to the Geocentric Celestial Reference
/// System (IAU, 2000) and the p-vector `V(date)` is with respect to the mean
/// equatorial triad of the given date.
///
/// # References:
///
/// * Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855
///
/// * IAU: Trans. International Astronomical Union, Vol. XXIVB; Proc. 24th
/// General Assembly, Manchester, UK. Resolutions B1.3, B1.6. (2000)
///
/// * Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981
///