1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
|
/*
*+
* Name:
* palFk524
* Purpose:
* Convert J2000.0 FK5 star data to B1950.0 FK4.
* Language:
* Starlink ANSI C
* Type of Module:
* Library routine
* Invocation:
* palFk524( double r2000, double d2000, double dr2000, double dd2000,
* double p2000, double v2000, double *r1950, double *d1950,
* double *dr1950, double *dd1950, double *p1950, double *v1950 )
* Arguments:
* r2000 = double (Given)
* J2000.0 FK5 RA (radians).
* d2000 = double (Given)
* J2000.0 FK5 Dec (radians).
* dr2000 = double (Given)
* J2000.0 FK5 RA proper motion (rad/Jul.yr)
* dd2000 = double (Given)
* J2000.0 FK5 Dec proper motion (rad/Jul.yr)
* p2000 = double (Given)
* J2000.0 FK5 parallax (arcsec)
* v2000 = double (Given)
* J2000.0 FK5 radial velocity (km/s, +ve = moving away)
* r1950 = double * (Returned)
* B1950.0 FK4 RA (radians).
* d1950 = double * (Returned)
* B1950.0 FK4 Dec (radians).
* dr1950 = double * (Returned)
* B1950.0 FK4 RA proper motion (rad/Jul.yr)
* dd1950 = double * (Returned)
* B1950.0 FK4 Dec proper motion (rad/Jul.yr)
* p1950 = double * (Returned)
* B1950.0 FK4 parallax (arcsec)
* v1950 = double * (Returned)
* B1950.0 FK4 radial velocity (km/s, +ve = moving away)
* Description:
* This function converts stars from the IAU 1976, FK5, Fricke
* system, to the Bessel-Newcomb, FK4 system. The precepts
* of Smith et al (Ref 1) are followed, using the implementation
* by Yallop et al (Ref 2) of a matrix method due to Standish.
* Kinoshita's development of Andoyer's post-Newcomb precession is
* used. The numerical constants from Seidelmann et al (Ref 3) are
* used canonically.
* Notes:
* - The proper motions in RA are dRA/dt rather than
* cos(Dec)*dRA/dt, and are per year rather than per century.
* - Note that conversion from Julian epoch 2000.0 to Besselian
* epoch 1950.0 only is provided for. Conversions involving
* other epochs will require use of the appropriate precession,
* proper motion, and E-terms routines before and/or after
* FK524 is called.
* - In the FK4 catalogue the proper motions of stars within
* 10 degrees of the poles do not embody the differential
* E-term effect and should, strictly speaking, be handled
* in a different manner from stars outside these regions.
* However, given the general lack of homogeneity of the star
* data available for routine astrometry, the difficulties of
* handling positions that may have been determined from
* astrometric fields spanning the polar and non-polar regions,
* the likelihood that the differential E-terms effect was not
* taken into account when allowing for proper motion in past
* astrometry, and the undesirability of a discontinuity in
* the algorithm, the decision has been made in this routine to
* include the effect of differential E-terms on the proper
* motions for all stars, whether polar or not. At epoch 2000,
* and measuring on the sky rather than in terms of dRA, the
* errors resulting from this simplification are less than
* 1 milliarcsecond in position and 1 milliarcsecond per
* century in proper motion.
*
* References:
* - Smith, C.A. et al, 1989. "The transformation of astrometric
* catalog systems to the equinox J2000.0". Astron.J. 97, 265.
* - Yallop, B.D. et al, 1989. "Transformation of mean star places
* from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space".
* Astron.J. 97, 274.
* - Seidelmann, P.K. (ed), 1992. "Explanatory Supplement to
* the Astronomical Almanac", ISBN 0-935702-68-7.
* Authors:
* PTW: Pat Wallace (STFC)
* DSB: David Berry (JAC, Hawaii)
* {enter_new_authors_here}
* History:
* 2012-02-13 (DSB):
* Initial version with documentation taken from Fortran SLA
* Adapted with permission from the Fortran SLALIB library.
* {enter_further_changes_here}
* Copyright:
* Copyright (C) 1995 Rutherford Appleton Laboratory
* Copyright (C) 2012 Science and Technology Facilities Council.
* All Rights Reserved.
* Licence:
* This program is free software: you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation, either
* version 3 of the License, or (at your option) any later
* version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General
* License along with this program. If not, see
* <http://www.gnu.org/licenses/>.
* Bugs:
* {note_any_bugs_here}
*-
*/
#include "pal.h"
#include "palmac.h"
#include "math.h"
void palFk524( double r2000, double d2000, double dr2000, double dd2000,
double p2000, double v2000, double *r1950, double *d1950,
double *dr1950, double *dd1950, double *p1950, double *v1950 ){
/* Local Variables; */
double r, d, ur, ud, px, rv;
double sr, cr, sd, cd, x, y, z, w;
double v1[ 6 ], v2[ 6 ];
double xd, yd, zd;
double rxyz, wd, rxysq, rxy;
int i, j;
/* Small number to avoid arithmetic problems. */
static const double tiny = 1.0E-30;
/* Canonical constants (see references). Constant vector and matrix. */
double a[ 6 ] = { -1.62557E-6, -0.31919E-6, -0.13843E-6,
+1.245E-3, -1.580E-3, -0.659E-3 };
double emi[ 6 ][ 6 ] = {
{ 0.9999256795, 0.0111814828, 0.0048590039,
-0.00000242389840, -0.00000002710544, -0.00000001177742},
{-0.0111814828, 0.9999374849, -0.0000271771,
0.00000002710544, -0.00000242392702, 0.00000000006585 },
{-0.0048590040, -0.0000271557, 0.9999881946,
0.00000001177742, 0.00000000006585, -0.00000242404995 },
{-0.000551, 0.238509, -0.435614,
0.99990432, 0.01118145, 0.00485852 },
{-0.238560, -0.002667, 0.012254,
-0.01118145, 0.99991613, -0.00002717},
{ 0.435730, -0.008541, 0.002117,
-0.00485852, -0.00002716, 0.99996684 } };
/* Pick up J2000 data (units radians and arcsec/JC). */
r = r2000;
d = d2000;
ur = dr2000*PAL__PMF;
ud = dd2000*PAL__PMF;
px = p2000;
rv = v2000;
/* Spherical to Cartesian. */
sr = sin( r );
cr = cos( r );
sd = sin( d );
cd = cos( d );
x = cr*cd;
y = sr*cd;
z = sd;
w = PAL__VF*rv*px;
v1[ 0 ] = x;
v1[ 1 ] = y;
v1[ 2 ] = z;
v1[ 3 ] = -ur*y - cr*sd*ud + w*x;
v1[ 4 ] = ur*x - sr*sd*ud + w*y;
v1[ 5 ] = cd*ud + w*z;
/* Convert position+velocity vector to BN system. */
for( i = 0; i < 6; i++ ) {
w = 0.0;
for( j = 0; j < 6; j++ ) {
w += emi[ i ][ j ]*v1[ j ];
}
v2[ i ] = w;
}
/* Position vector components and magnitude. */
x = v2[ 0 ];
y = v2[ 1 ];
z = v2[ 2 ];
rxyz = sqrt( x*x + y*y + z*z );
/* Apply E-terms to position. */
w = x*a[ 0 ] + y*a[ 1 ] + z*a[ 2 ];
x += a[ 0 ]*rxyz - w*x;
y += a[ 1 ]*rxyz - w*y;
z += a[ 2 ]*rxyz - w*z;
/* Recompute magnitude. */
rxyz = sqrt( x*x + y*y + z*z );
/* Apply E-terms to both position and velocity. */
x = v2[ 0 ];
y = v2[ 1 ];
z = v2[ 2 ];
w = x*a[ 0 ] + y*a[ 1 ] + z*a[ 2 ];
wd = x*a[ 3 ] + y*a[ 4 ] + z*a[ 5 ];
x += a[ 0 ]*rxyz - w*x;
y += a[ 1 ]*rxyz - w*y;
z += a[ 2 ]*rxyz - w*z;
xd = v2[ 3 ] + a[ 3 ]*rxyz - wd*x;
yd = v2[ 4 ] + a[ 4 ]*rxyz - wd*y;
zd = v2[ 5 ] + a[ 5 ]*rxyz - wd*z;
/* Convert to spherical. */
rxysq = x*x + y*y;
rxy = sqrt( rxysq );
if( x == 0.0 && y == 0.0 ) {
r = 0.0;
} else {
r = atan2( y, x );
if( r < 0.0 ) r += PAL__D2PI;
}
d = atan2( z, rxy );
if( rxy > tiny ) {
ur = ( x*yd - y*xd )/rxysq;
ud = ( zd*rxysq - z*( x*xd + y*yd ) )/( ( rxysq + z*z )*rxy );
}
/* Radial velocity and parallax. */
if( px > tiny ) {
rv = ( x*xd + y*yd + z*zd )/( px*PAL__VF*rxyz );
px /= rxyz;
}
/* Return results. */
*r1950 = r;
*d1950 = d;
*dr1950 = ur/PAL__PMF;
*dd1950 = ud/PAL__PMF;
*p1950 = px;
*v1950 = rv;
}
|
