source: trunk/MagicSoft/slalib/fk524.c@ 8123

Last change on this file since 8123 was 731, checked in by tbretz, 25 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 8.1 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaFk524 ( double r2000, double d2000, double dr2000,
4 double dd2000, double p2000, double v2000,
5 double *r1950, double *d1950, double *dr1950,
6 double *dd1950, double *p1950, double *v1950 )
7/*
8** - - - - - - - - -
9** s l a F k 5 2 4
10** - - - - - - - - -
11**
12** Convert J2000.0 FK5 star data to B1950.0 FK4.
13**
14** (double precision)
15**
16** This routine converts stars from the new, IAU 1976, FK5, Fricke
17** system, to the old, Bessel-Newcomb, FK4 system. The precepts
18** of Smith et al (Ref 1) are followed, using the implementation
19** by Yallop et al (Ref 2) of a matrix method due to Standish.
20** Kinoshita's development of Andoyer's post-Newcomb precession is
21** used. The numerical constants from Seidelmann et al (Ref 3) are
22** used canonically.
23**
24** Given: (all J2000.0,FK5)
25** r2000,d2000 double J2000.0 RA,Dec (rad)
26** dr2000,dd2000 double J2000.0 proper motions (rad/Jul.yr)
27** p2000 double parallax (arcsec)
28** v2000 double radial velocity (km/s, +ve = moving away)
29**
30** Returned: (all B1950.0,FK4)
31** *r1950,*d1950 double B1950.0 RA,Dec (rad)
32** *dr1950,*dd1950 double B1950.0 proper motions (rad/trop.yr)
33** *p1950 double parallax (arcsec)
34** *v1950 double radial velocity (km/s, +ve = moving away)
35**
36** Notes:
37**
38** 1) The proper motions in RA are dRA/dt rather than
39** cos(Dec)*dRA/dt, and are per year rather than per century.
40**
41** 2) Note that conversion from Julian epoch 2000.0 to Besselian
42** epoch 1950.0 only is provided for. Conversions involving
43** other epochs will require use of the appropriate precession,
44** proper motion, and E-terms routines before and/or after
45** FK524 is called.
46**
47** 3) In the FK4 catalogue the proper motions of stars within
48** 10 degrees of the poles do not embody the differential
49** E-term effect and should, strictly speaking, be handled
50** in a different manner from stars outside these regions.
51** However, given the general lack of homogeneity of the star
52** data available for routine astrometry, the difficulties of
53** handling positions that may have been determined from
54** astrometric fields spanning the polar and non-polar regions,
55** the likelihood that the differential E-terms effect was not
56** taken into account when allowing for proper motion in past
57** astrometry, and the undesirability of a discontinuity in
58** the algorithm, the decision has been made in this routine to
59** include the effect of differential E-terms on the proper
60** motions for all stars, whether polar or not. At epoch 2000,
61** and measuring on the sky rather than in terms of dRA, the
62** errors resulting from this simplification are less than
63** 1 milliarcsecond in position and 1 milliarcsecond per
64** century in proper motion.
65**
66** References:
67**
68** 1 Smith, C.A. et al, 1989. "The transformation of astrometric
69** catalog systems to the equinox J2000.0". Astron.J. 97, 265.
70**
71** 2 Yallop, B.D. et al, 1989. "Transformation of mean star places
72** from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space".
73** Astron.J. 97, 274.
74**
75** 3 Seidelmann, P.K. (ed), 1992. "Explanatory Supplement to
76** the Astronomical Almanac", ISBN 0-935702-68-7.
77**
78** Defined in slamac.h: D2PI
79**
80** Last revision: 20 December 1993
81**
82** Copyright P.T.Wallace. All rights reserved.
83*/
84{
85
86/* Miscellaneous */
87 double r, d, ur, ud, px, rv;
88 double sr, cr, sd, cd, x, y, z, w;
89 double v1[6], v2[6];
90 double xd, yd, zd;
91 double rxyz, wd, rxysq, rxy;
92 int i,j;
93
94/* Radians per year to arcsec per century */
95 static double pmf = 100.0 * 60.0 * 60.0 * 360.0 / D2PI;
96
97/* Small number to avoid arithmetic problems */
98 static double tiny = 1.0e-30;
99
100/*
101** Canonical constants (see references)
102*/
103
104/*
105** km per sec to AU per tropical century
106** = 86400 * 36524.2198782 / 1.49597870e8
107*/
108 static double vf = 21.095;
109
110/* Constant vector and matrix (by rows) */
111 static double a[6] = { -1.62557e-6, -0.31919e-6, -0.13843e-6,
112 1.245e-3, -1.580e-3, -0.659e-3 };
113
114 static double emi[6][6] =
115 {
116 { 0.9999256795, /* emi[0][0] */
117 0.0111814828, /* emi[0][1] */
118 0.0048590039, /* emi[0][2] */
119 -0.00000242389840, /* emi[0][3] */
120 -0.00000002710544, /* emi[0][4] */
121 -0.00000001177742 }, /* emi[0][5] */
122
123 { -0.0111814828, /* emi[1][0] */
124 0.9999374849, /* emi[1][1] */
125 -0.0000271771, /* emi[1][2] */
126 0.00000002710544, /* emi[1][3] */
127 -0.00000242392702, /* emi[1][4] */
128 0.00000000006585 }, /* emi[1][5] */
129
130 { -0.0048590040, /* emi[2][0] */
131 -0.0000271557, /* emi[2][1] */
132 0.9999881946, /* emi[2][2] */
133 0.00000001177742, /* emi[2][3] */
134 0.00000000006585, /* emi[2][4] */
135 -0.00000242404995 }, /* emi[2][5] */
136
137 { -0.000551, /* emi[3][0] */
138 0.238509, /* emi[3][1] */
139 -0.435614, /* emi[3][2] */
140 0.99990432, /* emi[3][3] */
141 0.01118145, /* emi[3][4] */
142 0.00485852 }, /* emi[3][5] */
143
144 { -0.238560, /* emi[4][0] */
145 -0.002667, /* emi[4][1] */
146 0.012254, /* emi[4][2] */
147 -0.01118145, /* emi[4][3] */
148 0.99991613, /* emi[4][4] */
149 -0.00002717 }, /* emi[4][5] */
150
151 { 0.435730, /* emi[5][0] */
152 -0.008541, /* emi[5][1] */
153 0.002117, /* emi[5][2] */
154 -0.00485852, /* emi[5][3] */
155 -0.00002716, /* emi[5][4] */
156 0.99996684 } /* emi[5][5] */
157 };
158
159/* Pick up J2000 data (units radians and arcsec/JC) */
160 r = r2000;
161 d = d2000;
162 ur = dr2000 * pmf;
163 ud = dd2000 * pmf;
164 px = p2000;
165 rv = v2000;
166
167/* Spherical to Cartesian */
168 sr = sin ( r );
169 cr = cos ( r );
170 sd = sin ( d );
171 cd = cos ( d );
172
173 x = cr * cd;
174 y = sr * cd;
175 z = sd;
176
177 w = vf * rv * px;
178
179 v1[0] = x;
180 v1[1] = y;
181 v1[2] = z;
182
183 v1[3] = - ur * y - cr * sd * ud + w * x;
184 v1[4] = ur * x - sr * sd * ud + w * y;
185 v1[5] = cd * ud + w * z;
186
187/* Convert position+velocity vector to BN system */
188 for ( i = 0; i < 6; i++ ) {
189 w = 0.0;
190 for ( j = 0; j < 6; j++ ) {
191 w += emi[i][j] * v1[j];
192 }
193 v2[i] = w;
194 }
195
196/* Position vector components and magnitude */
197 x = v2[0];
198 y = v2[1];
199 z = v2[2];
200 rxyz = sqrt ( x * x + y * y + z * z );
201
202/* Include e-terms */
203 w = x * a[0] + y * a[1] + z * a[2];
204 x += a[0] * rxyz - w * x;
205 y += a[1] * rxyz - w * y;
206 z += a[2] * rxyz - w * z;
207
208/* Recompute magnitude */
209 rxyz = sqrt ( x * x + y * y + z * z );
210
211/* Apply E-terms to both position and velocity */
212 x = v2[0];
213 y = v2[1];
214 z = v2[2];
215 w = x * a[0] + y * a[1] + z * a[2];
216 wd = x * a[3] + y * a[4] + z * a[5];
217 x += a[0] * rxyz - w * x;
218 y += a[1] * rxyz - w * y;
219 z += a[2] * rxyz - w * z;
220 xd = v2[3] + a[3] * rxyz - wd * x;
221 yd = v2[4] + a[4] * rxyz - wd * y;
222 zd = v2[5] + a[5] * rxyz - wd * z;
223
224/* Convert to spherical */
225 rxysq = x * x + y * y;
226 rxy = sqrt ( rxysq );
227
228 if ( ( x == 0.0 ) && ( y == 0.0 ) ) {
229 r = 0.0;
230 } else {
231 r = atan2 ( y, x );
232 if ( r < 0.0 ) r += D2PI;
233 }
234 d = atan2 ( z, rxy );
235
236 if (rxy > tiny) {
237 ur = ( x * yd - y * xd ) / rxysq;
238 ud = ( zd * rxysq - z * ( x * xd + y * yd ) ) /
239 ( ( rxysq + z * z ) * rxy );
240 }
241
242/* Radial velocity and parallax */
243 if ( px > tiny )
244 {
245 rv = ( x * xd + y * yd + z * zd ) / ( px * vf * rxyz );
246 px = px / rxyz;
247 }
248
249/* Return results */
250 *r1950 = r;
251 *d1950 = d;
252 *dr1950 = ur / pmf;
253 *dd1950 = ud / pmf;
254 *v1950 = rv;
255 *p1950 = px;
256}
Note: See TracBrowser for help on using the repository browser.