sofars/pnp/
bp00.rs

1use crate::consts::{DAS2R, DJ00, DJC};
2use crate::pnp::{bi00, pr00};
3use crate::vm::{ir, rx, rxr, ry, rz};
4
5///  Frame bias and precession matrices, IAU 2000.
6///
7///  Frame bias and precession, IAU 2000.
8///
9///  This function is part of the International Astronomical Union's
10///  SOFA (Standards of Fundamental Astronomy) software collection.
11///
12///  Status:  canonical model.
13///
14///  Given:
15///     date1,date2  double         TT as a 2-part Julian Date (Note 1)
16///
17///  Returned (function value):
18///     (rb, rp, rbp) ([[f64; 3]; 3], [[f64; 3]; 3], [[f64; 3]; 3])
19///
20///     rb           frame bias matrix (Note 2)
21///     rp           precession matrix (Note 3)
22///     rbp          bias-precession matrix (Note 4)
23///
24///  Notes:
25///
26///  1) The TT date date1+date2 is a Julian Date, apportioned in any
27///     convenient way between the two arguments.  For example,
28///     JD(TT)=2450123.7 could be expressed in any of these ways,
29///     among others:
30///
31///             date1         date2
32///
33///         2450123.7           0.0       (JD method)
34///         2451545.0       -1421.3       (J2000 method)
35///         2400000.5       50123.2       (MJD method)
36///         2450123.5           0.2       (date & time method)
37///
38///     The JD method is the most natural and convenient to use in
39///     cases where the loss of several decimal digits of resolution
40///     is acceptable.  The J2000 method is best matched to the way
41///     the argument is handled internally and will deliver the
42///     optimum resolution.  The MJD method and the date & time methods
43///     are both good compromises between resolution and convenience.
44///
45///  2) The matrix rb transforms vectors from GCRS to mean J2000.0 by
46///     applying frame bias.
47///
48///  3) The matrix rp transforms vectors from J2000.0 mean equator and
49///     equinox to mean equator and equinox of date by applying
50///     precession.
51///
52///  4) The matrix rbp transforms vectors from GCRS to mean equator and
53///     equinox of date by applying frame bias then precession.  It is
54///     the product rp x rb.
55///
56///  Called:
57///     iauBi00      frame bias components, IAU 2000
58///     iauPr00      IAU 2000 precession adjustments
59///     iauIr        initialize r-matrix to identity
60///     iauRx        rotate around X-axis
61///     iauRy        rotate around Y-axis
62///     iauRz        rotate around Z-axis
63///     iauCr        copy r-matrix
64///     iauRxr       product of two r-matrices
65///
66///  Reference:
67///     "Expressions for the Celestial Intermediate Pole and Celestial
68///     Ephemeris Origin consistent with the IAU 2000A precession-
69///     nutation model", Astron.Astrophys. 400, 1145-1154 (2003)
70///
71///     n.b. The celestial ephemeris origin (CEO) was renamed "celestial
72///          intermediate origin" (CIO) by IAU 2006 Resolution 2.
73///
74///  This revision:  2021 May 11
75///
76///  SOFA release 2023-10-11
77///
78///  Copyright (C) 2023 IAU SOFA Board.  See notes at end.
79///
80pub fn bp00(date1: f64, date2: f64) -> ([[f64; 3]; 3], [[f64; 3]; 3], [[f64; 3]; 3]) {
81    /* J2000.0 obliquity (Lieske et al. 1977) */
82    const EPS0: f64 = 84381.448 * DAS2R;
83
84    /* Interval between fundamental epoch J2000.0 and current date (JC). */
85    let t = ((date1 - DJ00) + date2) / DJC;
86
87    /* Frame bias. */
88    let (dpsibi, depsbi, dra0) = bi00();
89
90    /* Precession angles (Lieske et al. 1977) */
91    let psia77 = (5038.7784 + (-1.07259 + (-0.001147) * t) * t) * t * DAS2R;
92    let oma77 = EPS0 + ((0.05127 + (-0.007726) * t) * t) * t * DAS2R;
93    let chia = (10.5526 + (-2.38064 + (-0.001125) * t) * t) * t * DAS2R;
94
95    /* Apply IAU 2000 precession corrections. */
96    let (dpsipr, depspr) = pr00(date1, date2);
97    let psia = psia77 + dpsipr;
98    let oma = oma77 + depspr;
99
100    /* Frame bias matrix: GCRS to J2000.0. */
101    let mut rb = [[0.0; 3]; 3];
102    ir(&mut rb);
103    rz(dra0, &mut rb);
104    ry(dpsibi * EPS0.sin(), &mut rb);
105    rx(-depsbi, &mut rb);
106
107    /* Precession matrix: J2000.0 to mean of date. */
108    let mut rp = [[0.0; 3]; 3];
109    ir(&mut rp);
110    rx(EPS0, &mut rp);
111    rz(-psia, &mut rp);
112    rx(-oma, &mut rp);
113    rz(chia, &mut rp);
114
115    /* Bias-precession matrix: GCRS to mean of date. */
116    let mut rbp = [[0.0; 3]; 3];
117    rxr(&rp, &rb, &mut rbp);
118
119    (rb, rp, rbp)
120}
sofars/pnp/bp00.rs

sofars/pnp/
bp00.rs
