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

#![allow(non_snake_case)]
use crate::consts::{D2PI, DAYSEC};
use crate::coords::gd2gc;
use crate::pnp::pom00;
use crate::vm::trxp;

///  Observatory position and velocity
///
///  Position and velocity of a terrestrial observing station.
///
///  This function is part of the International Astronomical Union's
///  SOFA (Standards of Fundamental Astronomy) software collection.
///
///  Status:  support function.
///
///  Given:
///  ```
///     elong   double       longitude (radians, east +ve, Note 1)
///     phi     double       latitude (geodetic, radians, Note 1)
///     hm      double       height above ref. ellipsoid (geodetic, m)
///     xp,yp   double       coordinates of the pole (radians, Note 2)
///     sp      double       the TIO locator s' (radians, Note 2)
///     theta   double       Earth rotation angle (radians, Note 3)
///  ```
///  Returned:
///  ```
///     pv      double[2][3] position/velocity vector (m, m/s, CIRS)
///  ```
///  Notes:
///
///  1) The terrestrial coordinates are with respect to the WGS84
///     reference ellipsoid.
///
///  2) xp and yp are the coordinates (in radians) of the Celestial
///     Intermediate Pole with respect to the International Terrestrial
///     Reference System (see IERS Conventions), measured along the
///     meridians 0 and 90 deg west respectively.  sp is the TIO locator
///     s', in radians, which positions the Terrestrial Intermediate
///     Origin on the equator.  For many applications, xp, yp and
///     (especially) sp can be set to zero.
///
///  3) If theta is Greenwich apparent sidereal time instead of Earth
///     rotation angle, the result is with respect to the true equator
///     and equinox of date, i.e. with the x-axis at the equinox rather
///     than the celestial intermediate origin.
///
///  4) The velocity units are meters per UT1 second, not per SI second.
///     This is unlikely to have any practical consequences in the modern
///     era.
///
///  5) No validation is performed on the arguments.  Error cases that
///     could lead to arithmetic exceptions are trapped by the iauGd2gc
///     function, and the result set to zeros.
///
///  References:
///
///     McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
///     IERS Technical Note No. 32, BKG (2004)
///
///     Urban, S. & Seidelmann, P. K. (eds), Explanatory Supplement to
///     the Astronomical Almanac, 3rd ed., University Science Books
///     (2013), Section 7.4.3.3.
///
///  Called:
///  ```
///     iauGd2gc     geodetic to geocentric transformation
///     iauPom00     polar motion matrix
///     iauTrxp      product of transpose of r-matrix and p-vector
///  ```
pub fn pvtob(
    elong: f64,
    phi: f64,
    hm: f64,
    xp: f64,
    yp: f64,
    sp: f64,
    theta: f64,
    pv: &mut [[f64; 3]; 2],
) {
    /* Earth rotation rate in radians per UT1 second */
    let OM = 1.00273781191135448 * D2PI / DAYSEC;

    /* Geodetic to geocentric transformation (WGS84). */
    let xyzm = gd2gc(1, elong, phi, hm).unwrap_or([0.0; 3]);

    /* Polar motion and TIO position. */
    let mut xyz = [0.0; 3];
    let rpm = pom00(xp, yp, sp);
    trxp(&rpm, &xyzm, &mut xyz);
    let x = xyz[0];
    let y = xyz[1];
    let z = xyz[2];

    /* Functions of ERA. */
    let s = theta.sin();
    let c = theta.cos();

    /* Position. */
    pv[0][0] = c * x - s * y;
    pv[0][1] = s * x + c * y;
    pv[0][2] = z;

    /* Velocity. */
    pv[1][0] = OM * (-s * x - c * y);
    pv[1][1] = OM * (c * x - s * y);
    pv[1][2] = 0.0;
}