sofars/pnp/
pn00.rs

1#![allow(unused_variables)]
2use crate::vm::{cr, rxr};
3
4use super::{bp00, numat, obl80, pr00};
5
6///  Bias/precession/nutation results, IAU 2000
7///
8///  Precession-nutation, IAU 2000 model:  a multi-purpose function,
9///  supporting classical (equinox-based) use directly and CIO-based
10///  use indirectly.
11///
12///  This function is part of the International Astronomical Union's
13///  SOFA (Standards of Fundamental Astronomy) software collection.
14///
15///  Status:  support function.
16///
17///  Given:
18///     date1,date2  double          TT as a 2-part Julian Date (Note 1)
19///     dpsi,deps    double          nutation (Note 2)
20///
21///  Returned:
22///     epsa         double          mean obliquity (Note 3)
23///     rb           double[3][3]    frame bias matrix (Note 4)
24///     rp           double[3][3]    precession matrix (Note 5)
25///     rbp          double[3][3]    bias-precession matrix (Note 6)
26///     rn           double[3][3]    nutation matrix (Note 7)
27///     rbpn         double[3][3]    GCRS-to-true matrix (Note 8)
28///
29///  Notes:
30///
31///  1) The TT date date1+date2 is a Julian Date, apportioned in any
32///     convenient way between the two arguments.  For example,
33///     JD(TT)=2450123.7 could be expressed in any of these ways,
34///     among others:
35///
36///            date1          date2
37///
38///         2450123.7           0.0       (JD method)
39///         2451545.0       -1421.3       (J2000 method)
40///         2400000.5       50123.2       (MJD method)
41///         2450123.5           0.2       (date & time method)
42///
43///     The JD method is the most natural and convenient to use in
44///     cases where the loss of several decimal digits of resolution
45///     is acceptable.  The J2000 method is best matched to the way
46///     the argument is handled internally and will deliver the
47///     optimum resolution.  The MJD method and the date & time methods
48///     are both good compromises between resolution and convenience.
49///
50///  2) The caller is responsible for providing the nutation components;
51///     they are in longitude and obliquity, in radians and are with
52///     respect to the equinox and ecliptic of date.  For high-accuracy
53///     applications, free core nutation should be included as well as
54///     any other relevant corrections to the position of the CIP.
55///
56///  3) The returned mean obliquity is consistent with the IAU 2000
57///     precession-nutation models.
58///
59///  4) The matrix rb transforms vectors from GCRS to J2000.0 mean
60///     equator and equinox by applying frame bias.
61///
62///  5) The matrix rp transforms vectors from J2000.0 mean equator and
63///     equinox to mean equator and equinox of date by applying
64///     precession.
65///
66///  6) The matrix rbp transforms vectors from GCRS to mean equator and
67///     equinox of date by applying frame bias then precession.  It is
68///     the product rp x rb.
69///
70///  7) The matrix rn transforms vectors from mean equator and equinox of
71///     date to true equator and equinox of date by applying the nutation
72///     (luni-solar + planetary).
73///
74///  8) The matrix rbpn transforms vectors from GCRS to true equator and
75///     equinox of date.  It is the product rn x rbp, applying frame
76///     bias, precession and nutation in that order.
77///
78///  9) It is permissible to re-use the same array in the returned
79///     arguments.  The arrays are filled in the order given.
80///
81///  Called:
82///     iauPr00      IAU 2000 precession adjustments
83///     iauObl80     mean obliquity, IAU 1980
84///     iauBp00      frame bias and precession matrices, IAU 2000
85///     iauCr        copy r-matrix
86///     iauNumat     form nutation matrix
87///     iauRxr       product of two r-matrices
88///
89///  Reference:
90///
91///     Capitaine, N., Chapront, J., Lambert, S. and Wallace, P.,
92///     "Expressions for the Celestial Intermediate Pole and Celestial
93///     Ephemeris Origin consistent with the IAU 2000A precession-
94///     nutation model", Astron.Astrophys. 400, 1145-1154 (2003)
95///
96///     n.b. The celestial ephemeris origin (CEO) was renamed "celestial
97///          intermediate origin" (CIO) by IAU 2006 Resolution 2.
98pub fn pn00(
99    date1: f64,
100    date2: f64,
101    dpsi: f64,
102    deps: f64,
103    epsa: &mut f64,
104    rb: &mut [[f64; 3]; 3],
105    rp: &mut [[f64; 3]; 3],
106    rbp: &mut [[f64; 3]; 3],
107    rn: &mut [[f64; 3]; 3],
108    rbpn: &mut [[f64; 3]; 3],
109) {
110    let rbpw = &mut [[0.0; 3]; 3];
111    let rnw = &mut [[0.0; 3]; 3];
112
113    /* IAU 2000 precession-rate adjustments. */
114    let (dpsipr, depspr) = &mut pr00(date1, date2);
115
116    /* Mean obliquity, consistent with IAU 2000 precession-nutation. */
117    *epsa = obl80(date1, date2) + *depspr;
118
119    /* Frame bias and precession matrices and their product. */
120    bp00(date1, date2, rb, rp, rbpw);
121
122    cr(rbpw, rbp);
123
124    /* Nutation matrix. */
125    numat(*epsa, dpsi, deps, rnw);
126    cr(rnw, rn);
127
128    /* Bias-precession-nutation matrix (classical). */
129    rxr(rnw, rbpw, rbpn);
130}