1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
use crate::rfam::*;
use super::ltp::*;
/// Long-term precession matrix, including ICRS frame bias.
///
/// Given:
/// epj double Julian epoch (TT)
///
/// Returned:
/// rpb double[3][3] precession-bias matrix, J2000.0 to date
///
/// Notes:
///
/// 1) The matrix is in the sense
///
/// P_date = rpb x P_ICRS,
///
/// where P_ICRS is a vector in the Geocentric Celestial Reference
/// System, and P_date is the vector with respect to the Celestial
/// Intermediate Reference System at that date but with nutation
/// neglected.
///
/// 2) A first order frame bias formulation is used, of sub-
/// microarcsecond accuracy compared with a full 3D rotation.
///
/// 3) The Vondrak et al. (2011, 2012) 400 millennia precession model
/// agrees with the IAU 2006 precession at J2000.0 and stays within
/// 100 microarcseconds during the 20th and 21st centuries. It is
/// accurate to a few arcseconds throughout the historical period,
/// worsening to a few tenths of a degree at the end of the
/// +/- 200,000 year time span.
///
/// References:
///
/// Vondrak, J., Capitaine, N. and Wallace, P., 2011, New precession
/// expressions, valid for long time intervals, Astron.Astrophys. 534,
/// A22
///
/// Vondrak, J., Capitaine, N. and Wallace, P., 2012, New precession
/// expressions, valid for long time intervals (Corrigendum),
/// Astron.Astrophys. 541, C1
///
/// This revision: 2021 May 11
pub fn ltpb(epj: f64, rpb: &mut[[f64;3];3])
{
/* Frame bias (IERS Conventions 2010, Eqs. 5.21 and 5.33) */
const DX: f64 = -0.016617 * URSA_DAS2R;
const DE: f64 = -0.0068192 * URSA_DAS2R;
const DR: f64 = -0.0146 * URSA_DAS2R;
let mut rp = [[0.0; 3];3];
/* Precession matrix. */
ltp(epj, &mut rp);
/* Apply the bias. */
for i in 0..3 {
rpb[i][0] = rp[i][0] - rp[i][1]*DR + rp[i][2]*DX;
rpb[i][1] = rp[i][0]*DR + rp[i][1] + rp[i][2]*DE;
rpb[i][2] = -rp[i][0]*DX - rp[i][1]*DE + rp[i][2];
}
/* Finished. */
}