1 /*
2 * - - - - - - - - - - -
3 * g a l _ s t a r p v
4 * - - - - - - - - - - -
5 *
6 * This routine is part of the General Astrodynamics Library
7 *
8 * Description:
9 *
10 * Convert star catalog coordinates to position+velocity vector.
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 * ra d right ascension (radians)
21 * dec d declination (radians)
22 * pmr d RA proper motion (radians/year)
23 * pmd d Dec proper motion (radians/year)
24 * px d parallax (arcseconds)
25 * rv d radial velocity (km/s, positive = receding)
26 *
27 * Returned (Note 2):
28 * pv d[2][3] pv-vector (AU, AU/day)
29 * gal_starpv i status:
30 * 0 = no warnings
31 * 1 = distance overridden (Note 6)
32 * 2 = excessive velocity (Note 7)
33 * 4 = solution didn't converge (Note 8)
34 * else = binary logical OR of the above
35 *
36 * Notes:
37 *
38 * 1) The star data accepted by this routine are "observables" for an
39 * imaginary observer at the solar-system barycenter. Proper motion
40 * and radial velocity are, strictly, in terms of barycentric
41 * coordinate time, TCB. For most practical applications, it is
42 * permissible to neglect the distinction between TCB and ordinary
43 * "proper" time on Earth (TT/TAI). The result will, as a rule, be
44 * limited by the intrinsic accuracy of the proper-motion and radial-
45 * velocity data; moreover, the pv-vector is likely to be merely an
46 * intermediate result, so that a change of time unit would cancel
47 * out overall.
48 *
49 * In accordance with normal star-catalog conventions, the object's
50 * right ascension and declination are freed from the effects of
51 * secular aberration. The frame, which is aligned to the catalog
52 * equator and equinox, is Lorentzian and centered on the SSB.
53 *
54 * 2) The resulting position and velocity pv-vector is with respect to
55 * the same frame and, like the catalog coordinates, is freed from
56 * the effects of secular aberration. Should the "coordinate
57 * direction", where the object was located at the catalog epoch, be
58 * required, it may be obtained by calculating the magnitude of the
59 * position vector pv[0][0-2] dividing by the speed of light in AU/day
60 * to give the light-time, and then multiplying the space velocity
61 * pv[1][0-2] by this light-time and adding the result to pv[0][0-2].
62 *
63 * Summarizing, the pv-vector returned is for most stars almost
64 * identical to the result of applying the standard geometrical
65 * "space motion" transformation. The differences, which are the
66 * subject of the Stumpff paper referenced below, are:
67 *
68 * (i) In stars with significant radial velocity and proper motion,
69 * the constantly changing light-time distorts the apparent proper
70 * motion. Note that this is a classical, not a relativistic,
71 * effect.
72 *
73 * (ii) The transformation complies with special relativity.
74 *
75 * 3) Care is needed with units. The star coordinates are in radians
76 * and the proper motions in radians per Julian year, but the
77 * parallax is in arcseconds; the radial velocity is in km/s, but
78 * the pv-vector result is in AU and AU/day.
79 *
80 * 4) The ra proper motion is in terms of coordinate angle, not true
81 * angle. If the catalog uses arcseconds for both ra and dec proper
82 * motions, the ra proper motion will need to be divided by cos(dec)
83 * before use.
84 *
85 * 5) Straight-line motion at constant speed, in the inertial frame,
86 * is assumed.
87 *
88 * 6) An extremely small (or zero or negative) parallax is interpreted
89 * to mean that the object is on the "celestial sphere", the radius
90 * of which is an arbitrary (large) value (see the constant PXMIN).
91 * When the distance is overridden in this way, the status, initially
92 * zero, has 1 added to it.
93 *
94 * 7) If the space velocity is a significant fraction of c (see the
95 * constant VMAX), it is arbitrarily set to zero. When this action
96 * occurs, 2 is added to the status.
97 *
98 * 8) The relativistic adjustment involves an iterative calculation.
99 * If the process fails to converge within a set number (IMAX) of
100 * iterations, 4 is added to the status.
101 *
102 * 9) The inverse transformation is performed by the routine gal_PVSTAR.
103 *
104 * Called:
105 *
106 * gal_s2pv spherical coordinates to pv-vector
107 * gal_pm modulus of p-vector
108 * gal_zp zero p-vector
109 * gal_pn decompose p-vector into modulus and direction
110 * gal_pdp scalar product of two p-vectors
111 * gal_sxp multiply p-vector by scalar
112 * gal_pmp p-vector minus p-vector
113 * gal_ppp p-vector plus p-vector
114 *
115 * References:
116 *
117 * Stumpff, P., Astron.Astrophys. 144, 232-240 (1985).
118 *
119 * This revision:
120 *
121 * 2006 November 13 ( c version 2008 June 8 )
122 *
123 *
124 * Copyright (C) 2008 Paul C. L. Willmott. See notes at end.
125 *
126 *-----------------------------------------------------------------------
127 */
128
129 #include <math.h>
130 #include "gal_const.h"
131 #include "gal_starpv.h"
132 #include "gal_s2pv.h"
133 #include "gal_pm.h"
134 #include "gal_zp.h"
135 #include "gal_pn.h"
136 #include "gal_pdp.h"
137 #include "gal_sxp.h"
138 #include "gal_pmp.h"
139 #include "gal_ppp.h"
140
141 int
gal_starpv(double ra,double dec,double pmr,double pmd,double px,double rv,double pv[2][3])142 gal_starpv
143 (
144 double ra,
145 double dec,
146 double pmr,
147 double pmd,
148 double px,
149 double rv,
150 double pv[2][3]
151 )
152 {
153
154 /*
155 * Smallest allowed parallax
156 */
157
158 const double PXMIN = 1e-7 ;
159
160 /*
161 * Largest allowed speed (fraction of c)
162 */
163
164 const double VMAX = 0.5 ;
165
166 /*
167 * AU (meters)
168 */
169
170 const double AUM = 149597870e3 ;
171
172 /*
173 * Speed of light (AU per day)
174 */
175
176 const double C = GAL_D2S / 499.004782 ;
177
178 /*
179 * Maximum number of iterations for relativistic solution
180 */
181
182 int i, iwarn ;
183
184 const int IMAX = 100 ;
185
186 double w, r, rd, rad, decd, v, vsr, vst, betst, betsr, bett,
187 betr, od, odel, dd, ddel, odd, oddel, d, del,
188 x[3], usr[3], ust[3], ur[3], ut[3] ;
189
190 /*
191 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
192 */
193
194 /*
195 * Distance (AU).
196 */
197
198 if ( px >= PXMIN ) {
199 w = px ;
200 iwarn = 0 ;
201 }
202 else {
203 w = PXMIN ;
204 iwarn = 1 ;
205 }
206 r = GAL_R2AS / w ;
207
208 /*
209 * Radial velocity (AU/day).
210 */
211
212 rd = GAL_D2S * rv * 1e3 / AUM ;
213
214 /*
215 * Proper motion (radian/day).
216 */
217
218 rad = pmr / GAL_DJY ;
219 decd = pmd / GAL_DJY ;
220
221 /*
222 * To pv-vector (AU,AU/day).
223 */
224
225 gal_s2pv ( ra, dec, r, rad, decd, rd, pv ) ;
226
227 /*
228 * If excessive velocity, arbitrarily set it to zero.
229 */
230
231 v = gal_pm ( &pv[1][0] ) ;
232 if ( v / C > VMAX ) {
233 gal_zp ( &pv[1][0] ) ;
234 iwarn += 2 ;
235 }
236
237 /*
238 * Isolate the radial component of the velocity (AU/day).
239 */
240
241 gal_pn ( &pv[0][0], &w, x ) ;
242 vsr = gal_pdp ( x, &pv[1][0] ) ;
243 gal_sxp ( vsr, x, usr ) ;
244
245 /*
246 * Isolate the transverse component of the velocity (AU/day).
247 */
248
249 gal_pmp ( &pv[1][0], usr, ust ) ;
250 vst = gal_pm ( ust ) ;
251
252 /*
253 * Special-relativity dimensionless parameters.
254 */
255
256 betsr = vsr / C ;
257 betst = vst / C ;
258
259 /*
260 * Determine the inertial-to-observed relativistic correction terms.
261 */
262
263 bett = betst ;
264 betr = betsr ;
265 for ( i = 0; i < IMAX; i++ ) {
266 d = 1.0 + betr ;
267 del = sqrt ( 1.0 - betr * betr - bett * bett ) - 1.0 ;
268 betr = d * betsr + del ;
269 bett = d * betst ;
270 if ( i > 0 ) {
271 dd = fabs ( d - od ) ;
272 ddel = fabs ( del - odel ) ;
273 if ( i > 1 && dd == odd && ddel == oddel ) {
274 break ;
275 }
276 if ( i >= IMAX - 1 ) {
277 iwarn += 4 ;
278 }
279 odd = dd ;
280 oddel = ddel ;
281 }
282 od = d ;
283 odel = del ;
284 }
285
286 /*
287 * Replace observed radial velocity with inertial value.
288 */
289
290 if ( betsr != 0.0 ) {
291 w = d + del / betsr ;
292 }
293 else {
294 w = 1.0 ;
295 }
296 gal_sxp ( w, usr, ur ) ;
297
298 /*
299 * Replace observed tangential velocity with inertial value.
300 */
301
302 gal_sxp ( d, ust, ut ) ;
303
304 /*
305 * Combine the two to obtain the inertial space velocity.
306 */
307
308 gal_ppp ( ur, ut, &pv[1][0] ) ;
309
310 /*
311 * Return the status.
312 */
313
314 return iwarn ;
315
316 /*
317 * Finished.
318 */
319
320 }
321
322 /*
323 * gal - General Astrodynamics Library
324 * Copyright (C) 2008 Paul C. L. Willmott
325 *
326 * This program is free software; you can redistribute it and/or modify
327 * it under the terms of the GNU General Public License as published by
328 * the Free Software Foundation; either version 2 of the License, or
329 * (at your option) any later version.
330 *
331 * This program is distributed in the hope that it will be useful,
332 * but WITHOUT ANY WARRANTY; without even the implied warranty of
333 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
334 * GNU General Public License for more details.
335 *
336 * You should have received a copy of the GNU General Public License along
337 * with this program; if not, write to the Free Software Foundation, Inc.,
338 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
339 *
340 * Contact:
341 *
342 * Paul Willmott
343 * vp9mu@amsat.org
344 */
345