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
67
68
69
70
71
72
73
74
use super::{ltpequ::*, ltpecl::*};
use crate::vector_matrix::{vector_ops::{pn::*, pxp::*, }};
/// Long-term precession matrix.
///
/// Given:
/// epj double Julian epoch (TT)
///
/// Returned:
/// rp double[3][3] precession matrix, J2000.0 to date
///
/// Notes:
///
/// 1) The matrix is in the sense
///
/// P_date = rp x P_J2000,
///
/// where P_J2000 is a vector with respect to the J2000.0 mean
/// equator and equinox and P_date is the same vector with respect to
/// the equator and equinox of epoch epj.
///
/// 2) 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.
///
/// Called:
/// eraLtpequ equator pole, long term
/// eraLtpecl ecliptic pole, long term
/// eraPxp vector product
/// eraPn normalize vector
///
/// 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 ltp(epj: f64, rp: &mut[[f64;3];3])
{
let mut peqr = [0.0; 3]; let mut pecl = [0.0; 3]; let mut v = [0.0; 3]; let mut w =0.0;
let mut eqx = [0.0; 3];
/* Equator pole (bottom row of matrix). */
ltpequ(epj, &mut peqr);
/* Ecliptic pole. */
ltpecl(epj, &mut pecl);
/* Equinox (top row of matrix). */
pxp(&peqr, &pecl, &mut v);
pn(&v, &mut w, &mut eqx);
/* Middle row of matrix. */
pxp(&peqr, &eqx, &mut v);
/* Assemble the matrix. */
for i in 0..3 {
rp[0][i] = eqx[i];
rp[1][i] = v[i];
rp[2][i] = peqr[i];
}
/* Finished. */
}