rfa 0.5.9

use crate::rfam::*;
use super::{pwf06::*, fw2m::*, };
use crate::vector_matrix::{copy_ext::*, matrix_ops::{rxr::*, tr::*,}};
///  Precession-nutation, IAU 2006 model:  a multi-purpose function,
///  supporting classical (equinox-based) use directly and CIO-based use
///  indirectly.
///
///  Given:
///     date1,date2  double          TT as a 2-part Julian Date (Note 1)
///     dpsi,deps    double          nutation (Note 2)
///
///  Returned:
///     epsa         double          mean obliquity (Note 3)
///     rb           double[3][3]    frame bias matrix (Note 4)
///     rp           double[3][3]    precession matrix (Note 5)
///     rbp          double[3][3]    bias-precession matrix (Note 6)
///     rn           double[3][3]    nutation matrix (Note 7)
///     rbpn         double[3][3]    GCRS-to-true matrix (Notes 8,9)
///
///  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 caller is responsible for providing the nutation components;
///      they are in longitude and obliquity, in radians and are with
///      respect to the equinox and ecliptic of date.  For high-accuracy
///      applications, free core nutation should be included as well as
///      any other relevant corrections to the position of the CIP.
///
///  3)  The returned mean obliquity is consistent with the IAU 2006
///      precession.
///
///  4)  The matrix rb transforms vectors from GCRS to J2000.0 mean
///      equator and equinox by applying frame bias.
///
///  5)  The matrix rp transforms vectors from J2000.0 mean equator and
///      equinox to mean equator and equinox of date by applying
///      precession.
///
///  6)  The matrix rbp transforms vectors from GCRS to mean equator and
///      equinox of date by applying frame bias then precession.  It is
///      the product rp x rb.
///
///  7)  The matrix rn transforms vectors from mean equator and equinox
///      of date to true equator and equinox of date by applying the
///      nutation (luni-solar + planetary).
///
///  8)  The matrix rbpn transforms vectors from GCRS to true equator and
///      equinox of date.  It is the product rn x rbp, applying frame
///      bias, precession and nutation in that order.
///
///  9)  The X,Y,Z coordinates of the Celestial Intermediate Pole are
///      elements (3,1-3) of the GCRS-to-true matrix, i.e. rbpn[2][0-2].
///
///  10) It is permissible to re-use the same array in the returned
///      arguments.  The arrays are filled in the stated order.
///
///  Called:
///     eraPfw06     bias-precession F-W angles, IAU 2006
///     eraFw2m      F-W angles to r-matrix
///     eraCr        copy r-matrix
///     eraTr        transpose r-matrix
///     eraRxr       product of two r-matrices
///
///  References:
///
///     Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855
///
///     Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981
///
///  This revision:  2021 May 11
pub fn pn06(date1: f64, date2: f64, dpsi: f64, deps: f64,
    epsa: &mut f64,
    rb: &mut [[f64; 3];3], rp: &mut [[f64; 3];3], rbp: &mut [[f64; 3];3],
    rn: &mut [[f64; 3];3], rbpn: &mut [[f64; 3];3])
{
    let mut gamb =0.0; let mut phib =0.0; let mut psib =0.0; let mut eps =0.0;
    let mut r1 = [[0.0; 3]; 3]; let mut r2 = [[0.0; 3]; 3]; let mut rt = [[0.0; 3]; 3];
 
 
 /* Bias-precession Fukushima-Williams angles of J2000.0 = frame bias. */
    pfw06(URSA_DJM0, URSA_DJM00, &mut gamb, &mut phib, &mut psib, &mut eps);
 
 /* B matrix. */
    fw2m(gamb, phib, psib, eps, &mut r1);
    cr(&r1, rb);
 
 /* Bias-precession Fukushima-Williams angles of date. */
    pfw06(date1, date2, &mut gamb, &mut phib, &mut psib, &mut eps);
 
 /* Bias-precession matrix. */
    fw2m(gamb, phib, psib, eps, &mut r2);
    cr(&r2, rbp);
 
 /* Solve for precession matrix. */
    tr(&r1, &mut rt);
    rxr(&r2, &rt, rp);
 
 /* Equinox-based bias-precession-nutation matrix. */
    fw2m(gamb, phib, psib + dpsi, eps + deps, &mut r1);
    cr(&r1, rbpn);
 
 /* Solve for nutation matrix. */
    tr(&r2, &mut rt);
    rxr(&r1, &rt, rn);
 
 /* Obliquity, mean of date. */
    *epsa = eps;
 
 /* Finished. */
 
 }