sofars 0.6.0

Pure Rust implementation of the IAU SOFA library
Documentation 
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
75
76
77
78
79

use super::prec76;
use crate::consts::DJ00;
use crate::vm::{cr, ir, ry, rz};

/// Precession matrix from J2000.0 to a specified date, IAU 1976 model.
///
/// This function is part of the International Astronomical Union's
/// SOFA (Standards of Fundamental Astronomy) software collection.
///
/// Status:  support function.
///
/// Given:
///    date1,date2 f64          ending date, TT (Note 1)
///
/// Returned:
///    rmatp       f64[3][3]    precession matrix, J2000.0 -> date1+date2
///
/// 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:
///
/// ```text
///           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 matrix operates in the sense V(date) = RMATP * V(J2000),
///    where the p-vector V(J2000) is with respect to the mean
///    equatorial triad of epoch J2000.0 and the p-vector V(date)
///    is with respect to the mean equatorial triad of the given
///    date.
///
/// 3) Though the matrix method itself is rigorous, the precession
///    angles are expressed through canonical polynomials which are
///    valid only for a limited time span.  In addition, the IAU 1976
///    precession rate is known to be imperfect.  The absolute accuracy
///    of the present formulation is better than 0.1 arcsec from
///    1960AD to 2040AD, better than 1 arcsec from 1640AD to 2360AD,
///    and remains below 3 arcsec for the whole of the period
///    500BC to 3000AD.  The errors exceed 10 arcsec outside the
///    range 1200BC to 3900AD, exceed 100 arcsec outside 4200BC to
///    5600AD and exceed 1000 arcsec outside 6800BC to 8200AD.
///
/// References:
///
///    Lieske, J.H., 1979, Astron.Astrophys. 73, 282.
///     equations (6) & (7), p283.
///
///    Kaplan,G.H., 1981. USNO circular no. 163, pA2.
pub fn pmat76(date1: f64, date2: f64) -> [[f64; 3]; 3] {
    let mut wmat = [[0.0; 3]; 3];
    let mut rmatp = [[0.0; 3]; 3];

    /* Precession Euler angles, J2000.0 to specified date. */
    let (zeta, z, theta) = prec76(DJ00, 0.0, date1, date2);

    /* Form the rotation matrix. */
    ir(&mut wmat);
    rz(-zeta, &mut wmat);
    ry(theta, &mut wmat);
    rz(-z, &mut wmat);
    cr(&wmat, &mut rmatp);

    rmatp
}