sofars/astro/pvstar.rs
1use crate::{
2 consts::{DAU, DAYSEC, DC, DJY, DR2AS},
3 vm::{anp, pdp, pm, pmp, pn, ppp, pv2s, sxp},
4};
5
6/// Space motion pv−vector to star catalog data
7///
8/// Convert star position+velocity vector to catalog coordinates.
9///
10/// This function is part of the International Astronomical Union's
11/// SOFA (Standards of Fundamental Astronomy) software collection.
12///
13/// Status: support function.
14///
15/// Given (Note 1):
16/// ```
17/// pv double[2][3] pv-vector (au, au/day)
18/// ```
19/// Returned (Note 2):
20/// ```
21/// ra double right ascension (radians)
22/// dec double declination (radians)
23/// pmr double RA proper motion (radians/year)
24/// pmd double Dec proper motion (radians/year)
25/// px double parallax (arcsec)
26/// rv double radial velocity (km/s, positive = receding)
27/// ```
28/// Returned (function value):
29/// ```
30/// int status:
31/// 0 = OK
32/// -1 = superluminal speed (Note 5)
33/// -2 = null position vector
34/// ```
35/// Notes:
36///
37/// 1) The specified pv-vector is the coordinate direction (and its rate
38/// of change) for the date at which the light leaving the star
39/// reached the solar-system barycenter.
40///
41/// 2) The star data returned by this function are "observables" for an
42/// imaginary observer at the solar-system barycenter. Proper motion
43/// and radial velocity are, strictly, in terms of barycentric
44/// coordinate time, TCB. For most practical applications, it is
45/// permissible to neglect the distinction between TCB and ordinary
46/// "proper" time on Earth (TT/TAI). The result will, as a rule, be
47/// limited by the intrinsic accuracy of the proper-motion and
48/// radial-velocity data; moreover, the supplied pv-vector is likely
49/// to be merely an intermediate result (for example generated by the
50/// function iauStarpv), so that a change of time unit will cancel
51/// out overall.
52///
53/// In accordance with normal star-catalog conventions, the object's
54/// right ascension and declination are freed from the effects of
55/// secular aberration. The frame, which is aligned to the catalog
56/// equator and equinox, is Lorentzian and centered on the SSB.
57///
58/// Summarizing, the specified pv-vector is for most stars almost
59/// identical to the result of applying the standard geometrical
60/// "space motion" transformation to the catalog data. The
61/// differences, which are the subject of the Stumpff paper cited
62/// below, are:
63///
64/// (i) In stars with significant radial velocity and proper motion,
65/// the constantly changing light-time distorts the apparent proper
66/// motion. Note that this is a classical, not a relativistic,
67/// effect.
68///
69/// (ii) The transformation complies with special relativity.
70///
71/// 3) Care is needed with units. The star coordinates are in radians
72/// and the proper motions in radians per Julian year, but the
73/// parallax is in arcseconds; the radial velocity is in km/s, but
74/// the pv-vector result is in au and au/day.
75///
76/// 4) The proper motions are the rate of change of the right ascension
77/// and declination at the catalog epoch and are in radians per Julian
78/// year. The RA proper motion is in terms of coordinate angle, not
79/// true angle, and will thus be numerically larger at high
80/// declinations.
81///
82/// 5) Straight-line motion at constant speed in the inertial frame is
83/// assumed. If the speed is greater than or equal to the speed of
84/// light, the function aborts with an error status.
85///
86/// 6) The inverse transformation is performed by the function iauStarpv.
87///
88/// Called:
89/// ```
90/// iauPn decompose p-vector into modulus and direction
91/// iauPdp scalar product of two p-vectors
92/// iauSxp multiply p-vector by scalar
93/// iauPmp p-vector minus p-vector
94/// iauPm modulus of p-vector
95/// iauPpp p-vector plus p-vector
96/// iauPv2s pv-vector to spherical
97/// iauAnp normalize angle into range 0 to 2pi
98/// ```
99/// Reference:
100///
101/// Stumpff, P., 1985, Astron.Astrophys. 144, 232-240.
102pub fn pvstar(pv: &[[f64; 3]; 2]) -> Result<[f64; 6], i32> {
103 let mut pv = *pv;
104 let d: f64;
105 let w: f64;
106 let del: f64;
107
108 // Isolate the radial component of the velocity (au/day, inertial).
109 let (_, pu) = pn(&pv[0]);
110 let vr = pdp(&pu, &pv[1]);
111 let ur = sxp(vr, &pu);
112
113 // Isolate the transverse component of the velocity (au/day, inertial).
114 let ut = pmp(&pv[1], &ur);
115 let vt = pm(ut);
116
117 // Special-relativity dimensionless parameters.
118 let bett = vt / DC;
119 let betr = vr / DC;
120
121 // The observed-to-inertial correction terms.
122 d = 1.0 + betr;
123 w = betr * betr + bett * bett;
124 if d == 0.0 || w > 1.0 {
125 return Err(-1);
126 }
127 del = -w / ((1.0 - w).sqrt() + 1.0);
128
129 // Scale inertial tangential velocity vector into observed (au/d).
130 let ust = sxp(1.0 / d, &ut);
131
132 // Compute observed radial velocity vector (au/d).
133 let usr = sxp(DC * (betr - del) / d, &pu);
134
135 // Combine the two to obtain the observed velocity vector.
136 pv[1] = ppp(&usr, &ust);
137
138 // Cartesian to spherical.
139 let (mut ra, dec, r, rad, decd, rd) = pv2s(&pv);
140 if r == 0.0 {
141 return Err(-2);
142 }
143
144 // Return RA in range 0 to 2pi.
145 ra = anp(ra);
146
147 // Return proper motions in radians per year.
148 let pmr = rad * DJY;
149 let pmd = decd * DJY;
150
151 // Return parallax in arcsec.
152 let px = DR2AS / r;
153
154 // Return radial velocity in km/s.
155 let rv = 1e-3 * rd * DAU / DAYSEC;
156
157 // Success.
158 Ok([ra, dec, pmr, pmd, px, rv])
159}