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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114

use crate::consts::{DJ00, DJY};
use crate::star::fk5hip;
use crate::vm::{anp, pv2s, rv2m, rxp, rxr, s2c, sxp, trxp, pxp};

/// Transform a Hipparcos star position into FK5 J2000.0, assuming
/// zero Hipparcos proper motion.
///
/// Status:  support function.
///
/// Given:
///    rh            f64    Hipparcos RA (radians)
///    dh            f64    Hipparcos Dec (radians)
///    date1,date2   f64    TDB date (Note 1)
///
/// Returned (all FK5, equinox J2000.0, date date1+date2):
///    r5            f64    RA (radians)
///    d5            f64    Dec (radians)
///    dr5           f64    RA proper motion (rad/year, Note 4)
///    dd5           f64    Dec proper motion (rad/year, Note 4)
///
/// 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:
///
///           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 proper motion in RA is dRA/dt rather than cos(Dec)*dRA/dt.
///
/// 3) The FK5 to Hipparcos transformation is modeled as a pure rotation
///    and spin;  zonal errors in the FK5 catalog are not taken into
///    account.
///
/// 4) It was the intention that Hipparcos should be a close
///    approximation to an inertial frame, so that distant objects have
///    zero proper motion;  such objects have (in general) non-zero
///    proper motion in FK5, and this function returns those fictitious
///    proper motions.
///
/// 5) The position returned by this function is in the FK5 J2000.0
///    reference system but at date date1+date2.
///
/// 6) See also fk52h, h2fk5, fk5hz.
///
/// Called:
///    s2c       spherical coordinates to unit vector
///    fk5hip    FK5 to Hipparcos rotation and spin
///    rxp       product of r-matrix and p-vector
///    sxp       multiply p-vector by scalar
///    rv2m      r-vector to r-matrix
///    rxr       product of two r-matrices
///    trxp      product of transpose of r-matrix and p-vector
///    pxp       vector product of two p-vectors
///    pv2s      pv-vector to spherical
///    anp       normalize angle into range 0 to 2pi
///
/// Reference:
///    F.Mignard & M.Froeschle, 2000, Astron.Astrophys. 354, 732-739.
pub fn hfk5z(rh: f64, dh: f64, date1: f64, date2: f64) -> (f64, f64, f64, f64) {
    let mut sh = [0.0; 3];
    let mut rst = [[0.0; 3]; 3];
    let mut r5ht = [[0.0; 3]; 3];
    let mut pv5e = [[0.0; 3]; 2];

    /* Time interval from fundamental epoch J2000.0 to given date (JY). */
    let t = ((date1 - DJ00) + date2) / DJY;

    /* Hipparcos barycentric position vector (normalized). */
    let ph = s2c(rh, dh);

    /* FK5 to Hipparcos orientation matrix and spin vector. */
    let (r5h, s5h) = fk5hip();

    /* Rotate the spin into the Hipparcos system. */
    rxp(&r5h, &s5h, &mut sh);

    /* Accumulated Hipparcos wrt FK5 spin over that interval. */
    let vst = sxp(t, &s5h);

    /* Express the accumulated spin as a rotation matrix. */
    rv2m(&vst, &mut rst);

    /* Rotation matrix:  accumulated spin, then FK5 to Hipparcos. */
    rxr(&r5h, &rst, &mut r5ht);

    /* De-orient & de-spin the Hipparcos position into FK5 J2000.0. */
    trxp(&r5ht, &ph, &mut pv5e[0]);

    /* Apply spin to the position giving a space motion. */
    let vv = pxp(&sh, &ph);

    /* De-orient & de-spin the Hipparcos space motion into FK5 J2000.0. */
    trxp(&r5ht, &vv, &mut pv5e[1]);

    /* FK5 position/velocity pv-vector to spherical. */
    let (w, d5, _, dr5, dd5, _) = pv2s(&pv5e);
    let r5 = anp(w);

    (r5, d5, dr5, dd5)
}