| 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 | }
|
|---|
