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