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}