1 /*
2 * - - - - - - - - - - -
3 * g a l _ p v s t a r
4 * - - - - - - - - - - -
5 *
6 * This routine is part of the General Astrodynamics Library
7 *
8 * Description:
9 *
10 * Convert star position+velocity vector to catalog coordinates.
11 *
12 * This routine is an independent translation of a FORTRAN routine
13 * that is part of IAU's SOFA software collection.
14 *
15 * Status:
16 *
17 * support routine.
18 *
19 * Given (Note 1):
20 * pv d[2][3] pv-vector (AU, AU/day)
21 *
22 * Returned (Note 2):
23 * *ra d right ascension (radians)
24 * *dec d declination (radians)
25 * *pmr d RA proper motion (radians/year)
26 * *pmd d Dec proper motion (radians/year)
27 * *px d parallax (arcsec)
28 * *rv d radial velocity (km/s, positive = receding)
29 * gal_pvstar i status:
30 * 0 = OK
31 * -1 = superluminal speed (Note 5)
32 * -2 = null position vector
33 *
34 * Notes:
35 *
36 * 1) The specified pv-vector is the coordinate direction (and its rate
37 * of change) for the epoch at which the light leaving the star
38 * reached the solar-system barycenter.
39 *
40 * 2) The star data returned by this routine are "observables" for an
41 * imaginary observer at the solar-system barycenter. Proper motion
42 * and radial velocity are, strictly, in terms of barycentric
43 * coordinate time, TCB. For most practical applications, it is
44 * permissible to neglect the distinction between TCB and ordinary
45 * "proper" time on Earth (TT/TAI). The result will, as a rule, be
46 * limited by the intrinsic accuracy of the proper-motion and radial-
47 * velocity data; moreover, the supplied pv-vector is likely to be
48 * merely an intermediate result (for example generated by the
49 * routine gal_starpv), so that a change of time unit will cancel
50 * out overall.
51 *
52 * In accordance with normal star-catalog conventions, the object's
53 * right ascension and declination are freed from the effects of
54 * secular aberration. The frame, which is aligned to the catalog
55 * equator and equinox, is Lorentzian and centered on the SSB.
56 *
57 * Summarizing, the specified pv-vector is for most stars almost
58 * identical to the result of applying the standard geometrical
59 * "space motion" transformation to the catalog data. The
60 * differences, which are the subject of the Stumpff paper cited
61 * below, are:
62 *
63 * (i) In stars with significant radial velocity and proper motion,
64 * the constantly changing light-time distorts the apparent proper
65 * motion. Note that this is a classical, not a relativistic,
66 * effect.
67 *
68 * (ii) The transformation complies with special relativity.
69 *
70 * 3) Care is needed with units. The star coordinates are in radians
71 * and the proper motions in radians per Julian year, but the
72 * parallax is in arcseconds; the radial velocity is in km/s, but
73 * the pv-vector result is in AU and AU/day.
74 *
75 * 4) The proper motions are the rate of change of the right ascension
76 * and declination at the catalog epoch and are in radians per Julian
77 * year. The RA proper motion is in terms of coordinate angle, not
78 * true angle, and will thus be numerically larger at high
79 * declinations.
80 *
81 * 5) Straight-line motion at constant speed in the inertial frame is
82 * assumed. If the speed is greater than or equal to the speed of
83 * light, the routine aborts with an error status.
84 *
85 * 6) The inverse transformation is performed by the routine gal_starpv.
86 *
87 * Called:
88 *
89 * gal_pn decompose p-vector into modulus and direction
90 * gal_pdp scalar product of two p-vectors
91 * gal_sxp multiply p-vector by scalar
92 * gal_pmp p-vector minus p-vector
93 * gal_pm modulus of p-vector
94 * gal_ppp p-vector plus p-vector
95 * gal_pv2s pv-vector to spherical
96 * gal_anp normalize angle into range 0 to 2pi
97 *
98 * References:
99 *
100 * Stumpff, P., Astron.Astrophys. 144, 232-240 (1985).
101 *
102 * This revision:
103 *
104 * 2006 November 13 ( c version 2008 June 8 )
105 *
106 *
107 * Copyright (C) 2008 Paul C. L. Willmott. See notes at end.
108 *
109 *-----------------------------------------------------------------------
110 */
111
112 #include <math.h>
113 #include "gal_const.h"
114 #include "gal_pvstar.h"
115 #include "gal_pn.h"
116 #include "gal_pdp.h"
117 #include "gal_sxp.h"
118 #include "gal_pmp.h"
119 #include "gal_pm.h"
120 #include "gal_ppp.h"
121 #include "gal_pv2s.h"
122 #include "gal_anp.h"
123
124 int
gal_pvstar(double pv[2][3],double * ra,double * dec,double * pmr,double * pmd,double * px,double * rv)125 gal_pvstar
126 (
127 double pv[2][3],
128 double *ra,
129 double *dec,
130 double *pmr,
131 double *pmd,
132 double *px,
133 double *rv
134 )
135 {
136
137 const double AUM = 149597870e3 ; /* AU (m) */
138 const double C = GAL_D2S / 499.004782e0 ; /* Speed of light (AU per day) */
139
140 double r, vr, vt, bett, betr, d, w, del, a, rad, decd, rd,
141 x[3], ur[3], ut[3], usr[3], ust[3] ;
142 int j ;
143
144 /*
145 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146 */
147
148 /*
149 * Isolate the radial component of the velocity (AU/day, inertial).
150 */
151
152 gal_pn ( &pv[0][0], &r, x ) ;
153 vr = gal_pdp ( x, &pv[1][0] ) ;
154 gal_sxp ( vr, x, ur ) ;
155
156 /*
157 * Isolate the transverse component of the velocity (AU/day, inertial).
158 */
159
160 gal_pmp ( &pv[1][0], ur, ut ) ;
161 vt = gal_pm ( ut ) ;
162
163 /*
164 * Special-relativity dimensionless parameters.
165 */
166
167 bett = vt / C ;
168 betr = vr / C ;
169
170 /*
171 * The inertial-to-observed correction terms.
172 */
173
174 d = 1.0 + betr ;
175 w = 1.0 - betr * betr - bett * bett ;
176 if ( d == 0.0 || w < 0.0 ) {
177 j = -1 ;
178 }
179 else {
180
181 del = sqrt ( w ) - 1.0 ;
182
183 /*
184 * Apply relativistic correction factor to radial velocity component.
185 */
186
187 if ( betr != 0.0 ) {
188 w = ( betr - del ) / ( betr * d ) ;
189 }
190 else {
191 w = 1.0 ;
192 }
193 gal_sxp ( w, ur, usr ) ;
194
195 /*
196 * Apply relativistic correction factor to tangential velocity component.
197 */
198
199 gal_sxp ( 1.0 / d, ut, ust ) ;
200
201 /*
202 * Combine the two to obtain the observed velocity vector (AU/day).
203 */
204
205 gal_ppp ( usr, ust, &pv[1][0] ) ;
206
207 /*
208 * Cartesian to spherical.
209 */
210
211 gal_pv2s ( pv, &a, dec, &r, &rad, &decd, &rd ) ;
212 if ( r == 0.0 ) {
213 j = -2 ;
214 }
215 else {
216
217 /*
218 * Return RA in range 0 to 2pi.
219 */
220
221 *ra = gal_anp ( a ) ;
222
223 /*
224 * Return proper motions in radians per year.
225 */
226
227 *pmr = rad * GAL_DJY ;
228 *pmd = decd * GAL_DJY ;
229
230 /*
231 * Return parallax in arcsec.
232 */
233
234 *px = GAL_R2AS / r ;
235
236 /*
237 * Return radial velocity in km/s.
238 */
239
240 *rv = 1e-3 * rd * AUM / GAL_D2S ;
241
242 /*
243 * OK status.
244 */
245
246 j = 0 ;
247
248 }
249
250 }
251
252 /*
253 * Exit.
254 */
255
256 return j ;
257
258 /*
259 * Finished.
260 */
261
262 }
263
264 /*
265 * gal - General Astrodynamics Library
266 * Copyright (C) 2008 Paul C. L. Willmott
267 *
268 * This program is free software; you can redistribute it and/or modify
269 * it under the terms of the GNU General Public License as published by
270 * the Free Software Foundation; either version 2 of the License, or
271 * (at your option) any later version.
272 *
273 * This program is distributed in the hope that it will be useful,
274 * but WITHOUT ANY WARRANTY; without even the implied warranty of
275 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
276 * GNU General Public License for more details.
277 *
278 * You should have received a copy of the GNU General Public License along
279 * with this program; if not, write to the Free Software Foundation, Inc.,
280 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
281 *
282 * Contact:
283 *
284 * Paul Willmott
285 * vp9mu@amsat.org
286 */
287