Skip to main content

sofars/astro/
ab.rs

1use crate::consts::SRS;
2use crate::vm::pdp;
3
4///  Apply stellar aberration
5///
6///  Apply aberration to transform natural direction into proper
7///  direction.
8///
9///  This function is part of the International Astronomical Union's
10///  SOFA (Standards of Fundamental Astronomy) software collection.
11///
12///  Status:  support function.
13///
14///  Given:
15///    pnat    double[3]   natural direction to the source (unit vector)
16///    v       double[3]   observer barycentric velocity in units of c
17///    s       double      distance between the Sun and the observer (au)
18///    bm1     double      sqrt(1-|v|^2): reciprocal of Lorenz factor
19///
20///  Returned:
21///    ppr     double[3]   proper direction to source (unit vector)
22///
23///  Notes:
24///
25///  1) The algorithm is based on Expr. (7.40) in the Explanatory
26///     Supplement (Urban & Seidelmann 2013), but with the following
27///     changes:
28///
29///     o  Rigorous rather than approximate normalization is applied.
30///
31///     o  The gravitational potential term from Expr. (7) in
32///        Klioner (2003) is added, taking into account only the Sun's
33///        contribution.  This has a maximum effect of about
34///        0.4 microarcsecond.
35///
36///  2) In almost all cases, the maximum accuracy will be limited by the
37///     supplied velocity.  For example, if the SOFA iauEpv00 function is
38///     used, errors of up to 5 microarcseconds could occur.
39///
40///  References:
41///
42///     Urban, S. & Seidelmann, P. K. (eds), Explanatory Supplement to
43///     the Astronomical Almanac, 3rd ed., University Science Books
44///     (2013).
45///
46///     Klioner, Sergei A., "A practical relativistic model for micro-
47///     arcsecond astrometry in space", Astr. J. 125, 1580-1597 (2003).
48///
49///  Called:
50///     iauPdp       scalar product of two p-vectors
51pub fn ab(pnat: &[f64; 3], v: &[f64; 3], s: f64, bm1: f64) -> [f64; 3] {
52    let mut p = [0.0; 3];
53    let mut r2 = 0.0;
54
55    let pdv = pdp(pnat, v);
56    let w1 = 1.0 + pdv / (1.0 + bm1);
57    let w2 = SRS / s; // SRS constant value
58
59    for i in 0..3 {
60        let w = pnat[i] * bm1 + w1 * v[i] + w2 * (v[i] - pdv * pnat[i]);
61        p[i] = w;
62        r2 += w * w;
63    }
64
65    let r = r2.sqrt();
66    let mut ppr = [0.0; 3];
67    for i in 0..3 {
68        ppr[i] = p[i] / r;
69    }
70
71    ppr
72}