FACT++  1.0
palFk524.c
Go to the documentation of this file.
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 
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.0E-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 }
int i
Definition: db_dim_client.c:21
static const double PAL__VF
Definition: palmac.h:106
#define PAL__PMF
Definition: palmac.h:110
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)
Definition: palFk524.c:132
static const double PAL__D2PI
Definition: palmac.h:66