1 /*** File wcscon.c
2 *** June 9, 2016
3 *** Doug Mink, Harvard-Smithsonian Center for Astrophysics
4 *** Some subroutines are based on Starlink subroutines by Patrick Wallace
5 *** Copyright (C) 1995-2016
6 *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA
7
8 This library is free software; you can redistribute it and/or
9 modify it under the terms of the GNU Lesser General Public
10 License as published by the Free Software Foundation; either
11 version 2 of the License, or (at your option) any later version.
12
13 This library is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 Lesser General Public License for more details.
17
18 You should have received a copy of the GNU Lesser General Public
19 License along with this library; if not, write to the Free Software
20 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21
22 Correspondence concerning WCSTools should be addressed as follows:
23 Internet email: jmink@cfa.harvard.edu
24 Postal address: Jessica Mink
25 Smithsonian Astrophysical Observatory
26 60 Garden St.
27 Cambridge, MA 02138 USA
28
29 * Module: wcscon.c (World Coordinate System conversion)
30 * Purpose: Convert between various sky coordinate systems
31 * Subroutine: wcscon (sys1,sys2,eq1,eq2,theta,phi,epoch)
32 * convert between coordinate systems
33 * Subroutine: wcsconp (sys1,sys2,eq1,eq2,ep1,ep2,dtheta,dphi,ptheta,pphi)
34 * convert coordinates and proper motion between coordinate systems
35 * Subroutine: wcsconv (sys1,sys2,eq1,eq2,ep1,ep2,dtheta,dphi,ptheta,pphi,px,rv)
36 * convert coordinates and proper motion between coordinate systems
37 * Subroutine: wcscsys (cstring) returns code for coordinate system in string
38 * Subroutine: wcsceq (wcstring) returns equinox in years from system string
39 * Subroutine: wcscstr (sys,equinox,epoch) returns system string from equinox
40 * Subroutine: fk524 (ra,dec) Convert J2000(FK5) to B1950(FK4) coordinates
41 * Subroutine: fk524e (ra, dec, epoch) (more accurate for known position epoch)
42 * Subroutine: fk524m (ra,dec,rapm,decpm) exact
43 * Subroutine: fk524pv (ra,dec,rapm,decpm,parallax,rv) more exact
44 * Subroutine: fk425 (ra,dec) Convert B1950(FK4) to J2000(FK5) coordinates
45 * Subroutine: fk425e (ra, dec, epoch) (more accurate for known position epoch)
46 * Subroutine: fk425m (ra, dec, rapm, decpm) exact
47 * Subroutine: fk425pv (ra,dec,rapm,decpm,parallax,rv) more exact
48 * Subroutine: fk42gal (dtheta,dphi) Convert B1950(FK4) to galactic coordinates
49 * Subroutine: fk52gal (dtheta,dphi) Convert J2000(FK5) to galactic coordinates
50 * Subroutine: gal2fk4 (dtheta,dphi) Convert galactic coordinates to B1950(FK4)
51 * Subroutine: gal2fk5 (dtheta,dphi) Convert galactic coordinates to J2000<FK5)
52 * Subroutine: fk42ecl (dtheta,dphi,epoch) Convert B1950(FK4) to ecliptic coordinates
53 * Subroutine: fk52ecl (dtheta,dphi,epoch) Convert J2000(FK5) to ecliptic coordinates
54 * Subroutine: ecl2fk4 (dtheta,dphi,epoch) Convert ecliptic coordinates to B1950(FK4)
55 * Subroutine: ecl2fk5 (dtheta,dphi,epoch) Convert ecliptic coordinates to J2000<FK5)
56 * Subroutine: fk5prec (ep0, ep1, ra, dec) Precession ep0 to ep1, FK5 system
57 * Subroutine: fk4prec (ep0, ep1, ra, dec) Precession ep0 to ep1, FK4 system
58 * Subroutine: d2v3 (rra, rdec, r, pos) RA and Dec in degrees, Distance to Cartesian
59 * Subroutine: v2d3 (pos, rra, rdec, r) Cartesian to RA and Dec in degrees, Distance
60 * Subroutine: s2v3 (rra, rdec, r, pos) RA, Dec, Distance to Cartesian
61 * Subroutine: v2s3 (pos, rra, rdec, r) Cartesian to RA, Dec, Distance
62 * Subroutine: rotmat (axes, rot1, rot2, rot3, matrix) Rotation angles to matrix
63 *
64 * Note: Proper motions are always in RA/Dec degrees/year; no cos(Dec) correction
65 */
66
67 #include <math.h>
68 #ifndef VMS
69 #include <stdlib.h>
70 #endif
71 #include <stdio.h> /* for fprintf() and sprintf() */
72 #include <ctype.h>
73 #include <string.h>
74 #include "wcs.h"
75
76 void fk524(), fk524e(), fk524m(), fk524pv();
77 void fk425(), fk425e(), fk425m(), fk425pv();
78 void fk42gal(), fk52gal(), gal2fk4(), gal2fk5();
79 void fk42ecl(), fk52ecl(), ecl2fk4(), ecl2fk5();
80
81 /* Convert from coordinate system sys1 to coordinate system sys2, converting
82 proper motions, too, and adding them if an epoch is specified */
83
84 void
wcsconp(sys1,sys2,eq1,eq2,ep1,ep2,dtheta,dphi,ptheta,pphi)85 wcsconp (sys1, sys2, eq1, eq2, ep1, ep2, dtheta, dphi, ptheta, pphi)
86
87 int sys1; /* Input coordinate system (J2000, B1950, ECLIPTIC, GALACTIC */
88 int sys2; /* Output coordinate system (J2000, B1950, ECLIPTIC, GALACTIC */
89 double eq1; /* Input equinox (default of sys1 if 0.0) */
90 double eq2; /* Output equinox (default of sys2 if 0.0) */
91 double ep1; /* Input Besselian epoch in years (for proper motion) */
92 double ep2; /* Output Besselian epoch in years (for proper motion) */
93 double *dtheta; /* Longitude or right ascension in degrees
94 Input in sys1, returned in sys2 */
95 double *dphi; /* Latitude or declination in degrees
96 Input in sys1, returned in sys2 */
97 double *ptheta; /* Longitude or right ascension proper motion in RA degrees/year
98 Input in sys1, returned in sys2 */
99 double *pphi; /* Latitude or declination proper motion in Dec degrees/year
100 Input in sys1, returned in sys2 */
101
102 {
103 void fk5prec(), fk4prec();
104
105 /* Set equinoxes if 0.0 */
106 if (eq1 == 0.0) {
107 if (sys1 == WCS_B1950)
108 eq1 = 1950.0;
109 else
110 eq1 = 2000.0;
111 }
112 if (eq2 == 0.0) {
113 if (sys2 == WCS_B1950)
114 eq2 = 1950.0;
115 else
116 eq2 = 2000.0;
117 }
118
119 /* Set epochs if 0.0 */
120 if (ep1 == 0.0) {
121 if (sys1 == WCS_B1950)
122 ep1 = 1950.0;
123 else
124 ep1 = 2000.0;
125 }
126 if (ep2 == 0.0) {
127 if (sys2 == WCS_B1950)
128 ep2 = 1950.0;
129 else
130 ep2 = 2000.0;
131 }
132
133 if (sys1 == WCS_ICRS && sys2 == WCS_ICRS)
134 eq2 = eq1;
135
136 if (sys1 == WCS_J2000 && sys2 == WCS_ICRS && eq1 == 2000.0) {
137 eq2 = eq1;
138 sys1 = sys2;
139 }
140
141 /* Set systems and equinoxes so that ICRS coordinates are not precessed */
142 if (sys1 == WCS_ICRS && sys2 == WCS_J2000 && eq2 == 2000.0) {
143 eq1 = eq2;
144 sys1 = sys2;
145 }
146
147 /* If systems and equinoxes are the same, add proper motion and return */
148 if (sys2 == sys1 && eq1 == eq2) {
149 if (ep1 != ep2) {
150 if (sys1 == WCS_J2000) {
151 *dtheta = *dtheta + ((ep2 - ep1) * *ptheta);
152 *dphi = *dphi + ((ep2 - ep1) * *pphi);
153 }
154 else if (sys1 == WCS_B1950) {
155 *dtheta = *dtheta + ((ep2 - ep1) * *ptheta);
156 *dphi = *dphi + ((ep2 - ep1) * *pphi);
157 }
158 }
159 if (eq1 != eq2) {
160 if (sys1 == WCS_B1950)
161 fk4prec (eq1, eq2, dtheta, dphi);
162 if (sys1 == WCS_J2000)
163 fk5prec (eq1, 2000.0, dtheta, dphi);
164 }
165 return;
166 }
167
168 /* Precess from input equinox to input system equinox, if necessary */
169 if (sys1 == WCS_B1950 && eq1 != 1950.0)
170 fk4prec (eq1, 1950.0, dtheta, dphi);
171 if (sys1 == WCS_J2000 && eq1 != 2000.0)
172 fk5prec (eq1, 2000.0, dtheta, dphi);
173
174 /* Convert to B1950 FK4 */
175 if (sys2 == WCS_B1950) {
176 if (sys1 == WCS_J2000) {
177 if (*ptheta != 0.0 || *pphi != 0.0) {
178 fk524m (dtheta, dphi, ptheta, pphi);
179 if (ep2 != 1950.0) {
180 *dtheta = *dtheta + ((ep2 - 1950.0) * *ptheta);
181 *dphi = *dphi + ((ep2 - 1950.0) * *pphi);
182 }
183 }
184 else if (ep2 != 1950.0)
185 fk524e (dtheta, dphi, ep2);
186 else
187 fk524 (dtheta, dphi);
188 }
189 else if (sys1 == WCS_GALACTIC)
190 gal2fk4 (dtheta, dphi);
191 else if (sys1 == WCS_ECLIPTIC)
192 ecl2fk4 (dtheta, dphi, ep2);
193 }
194
195 else if (sys2 == WCS_J2000) {
196 if (sys1 == WCS_B1950) {
197 if (*ptheta != 0.0 || *pphi != 0.0) {
198 fk425m (dtheta, dphi, ptheta, pphi);
199 if (ep2 != 2000.0) {
200 *dtheta = *dtheta + ((ep2 - 2000.0) * *ptheta);
201 *dphi = *dphi + ((ep2 - 2000.0) * *pphi);
202 }
203 }
204 else if (ep2 > 0.0)
205 fk425e (dtheta, dphi, ep2);
206 else
207 fk425 (dtheta, dphi);
208 }
209 else if (sys1 == WCS_GALACTIC)
210 gal2fk5 (dtheta, dphi);
211 else if (sys1 == WCS_ECLIPTIC)
212 ecl2fk5 (dtheta, dphi, ep2);
213 }
214
215 else if (sys2 == WCS_GALACTIC) {
216 if (sys1 == WCS_B1950) {
217 if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
218 *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
219 *dphi = *dphi + (*pphi * (ep2 - ep1));
220 }
221 fk42gal (dtheta, dphi);
222 }
223 else if (sys1 == WCS_J2000) {
224 if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
225 *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
226 *dphi = *dphi + (*pphi * (ep2 - ep1));
227 }
228 fk52gal (dtheta, dphi);
229 }
230 else if (sys1 == WCS_ECLIPTIC) {
231 ecl2fk5 (dtheta, dphi, ep2);
232 fk52gal (dtheta, dphi);
233 }
234 }
235
236 else if (sys2 == WCS_ECLIPTIC) {
237 if (sys1 == WCS_B1950) {
238 if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
239 *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
240 *dphi = *dphi + (*pphi * (ep2 - ep1));
241 }
242 if (ep2 > 0.0)
243 fk42ecl (dtheta, dphi, ep2);
244 else
245 fk42ecl (dtheta, dphi, 1950.0);
246 }
247 else if (sys1 == WCS_J2000) {
248 if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
249 *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
250 *dphi = *dphi + (*pphi * (ep2 - ep1));
251 }
252 fk52ecl (dtheta, dphi, ep2);
253 }
254 else if (sys1 == WCS_GALACTIC) {
255 gal2fk5 (dtheta, dphi);
256 fk52ecl (dtheta, dphi, ep2);
257 }
258 }
259
260 /* Precess to desired equinox, if necessary */
261 if (sys2 == WCS_B1950 && eq2 != 1950.0)
262 fk4prec (1950.0, eq2, dtheta, dphi);
263 if (sys2 == WCS_J2000 && eq2 != 2000.0)
264 fk5prec (2000.0, eq2, dtheta, dphi);
265
266 /* Keep latitude/declination between +90 and -90 degrees */
267 if (*dphi > 90.0) {
268 *dphi = 180.0 - *dphi;
269 *dtheta = *dtheta + 180.0;
270 }
271 else if (*dphi < -90.0) {
272 *dphi = -180.0 - *dphi;
273 *dtheta = *dtheta + 180.0;
274 }
275
276 /* Keep longitude/right ascension between 0 and 360 degrees */
277 if (*dtheta > 360.0)
278 *dtheta = *dtheta - 360.0;
279 else if (*dtheta < 0.0)
280 *dtheta = *dtheta + 360.0;
281 return;
282 }
283
284
285 /* Convert from coordinate system sys1 to coordinate system sys2, converting
286 proper motions, too, and adding them if an epoch is specified */
287
288 void
wcsconv(sys1,sys2,eq1,eq2,ep1,ep2,dtheta,dphi,ptheta,pphi,px,rv)289 wcsconv (sys1, sys2, eq1, eq2, ep1, ep2, dtheta, dphi, ptheta, pphi, px, rv)
290
291 int sys1; /* Input coordinate system (J2000, B1950, ECLIPTIC, GALACTIC */
292 int sys2; /* Output coordinate system (J2000, B1950, ECLIPTIC, GALACTIC */
293 double eq1; /* Input equinox (default of sys1 if 0.0) */
294 double eq2; /* Output equinox (default of sys2 if 0.0) */
295 double ep1; /* Input Besselian epoch in years (for proper motion) */
296 double ep2; /* Output Besselian epoch in years (for proper motion) */
297 double *dtheta; /* Longitude or right ascension in degrees
298 Input in sys1, returned in sys2 */
299 double *dphi; /* Latitude or declination in degrees
300 Input in sys1, returned in sys2 */
301 double *ptheta; /* Longitude or right ascension proper motion in degrees/year
302 Input in sys1, returned in sys2 */
303 double *pphi; /* Latitude or declination proper motion in degrees/year
304 Input in sys1, returned in sys2 */
305 double *px; /* Parallax in arcseconds */
306 double *rv; /* Radial velocity in km/sec */
307
308 {
309 void fk5prec(), fk4prec();
310
311 /* Set equinoxes if 0.0 */
312 if (eq1 == 0.0) {
313 if (sys1 == WCS_B1950)
314 eq1 = 1950.0;
315 else
316 eq1 = 2000.0;
317 }
318 if (eq2 == 0.0) {
319 if (sys2 == WCS_B1950)
320 eq2 = 1950.0;
321 else
322 eq2 = 2000.0;
323 }
324
325 /* Set epochs if 0.0 */
326 if (ep1 == 0.0) {
327 if (sys1 == WCS_B1950)
328 ep1 = 1950.0;
329 else
330 ep1 = 2000.0;
331 }
332 if (ep2 == 0.0) {
333 if (sys2 == WCS_B1950)
334 ep2 = 1950.0;
335 else
336 ep2 = 2000.0;
337 }
338
339 /* Set systems and equinoxes so that ICRS coordinates are not precessed */
340 if (sys1 == WCS_ICRS && sys2 == WCS_ICRS)
341 eq2 = eq1;
342
343 if (sys1 == WCS_J2000 && sys2 == WCS_ICRS && eq1 == 2000.0) {
344 eq2 = eq1;
345 sys1 = sys2;
346 }
347
348 if (sys1 == WCS_ICRS && sys2 == WCS_J2000 && eq2 == 2000.0) {
349 eq1 = eq2;
350 sys1 = sys2;
351 }
352
353 /* If systems and equinoxes are the same, add proper motion and return */
354 if (sys2 == sys1 && eq1 == eq2) {
355 if (ep1 != ep2) {
356 if (sys1 == WCS_J2000) {
357 *dtheta = *dtheta + ((ep2 - ep1) * *ptheta);
358 *dphi = *dphi + ((ep2 - ep1) * *pphi);
359 }
360 else if (sys1 == WCS_B1950) {
361 *dtheta = *dtheta + ((ep2 - ep1) * *ptheta);
362 *dphi = *dphi + ((ep2 - ep1) * *pphi);
363 }
364 }
365 return;
366 }
367
368 /* Precess from input equinox to input system equinox, if necessary */
369 if (eq1 != eq2) {
370 if (sys1 == WCS_B1950 && eq1 != 1950.0)
371 fk4prec (eq1, 1950.0, dtheta, dphi);
372 if (sys1 == WCS_J2000 && eq1 != 2000.0)
373 fk5prec (eq1, 2000.0, dtheta, dphi);
374 }
375
376 /* Convert to B1950 FK4 */
377 if (sys2 == WCS_B1950) {
378 if (sys1 == WCS_J2000) {
379 if (*ptheta != 0.0 || *pphi != 0.0) {
380 if (*px != 0.0 || *rv != 0.0)
381 fk524pv (dtheta, dphi, ptheta, pphi, px, rv);
382 else
383 fk524m (dtheta, dphi, ptheta, pphi);
384 if (ep1 == 2000.0)
385 ep1 = 1950.0;
386 if (ep2 != 1950.0) {
387 *dtheta = *dtheta + ((ep2 - 1950.0) * *ptheta);
388 *dphi = *dphi + ((ep2 - 1950.0) * *pphi);
389 }
390 }
391 else if (ep2 != 1950.0)
392 fk524e (dtheta, dphi, ep2);
393 else
394 fk524 (dtheta, dphi);
395 }
396 else if (sys1 == WCS_GALACTIC)
397 gal2fk4 (dtheta, dphi);
398 else if (sys1 == WCS_ECLIPTIC)
399 ecl2fk4 (dtheta, dphi, ep2);
400 }
401
402 else if (sys2 == WCS_J2000) {
403 if (sys1 == WCS_B1950) {
404 if (*ptheta != 0.0 || *pphi != 0.0) {
405 if (*px != 0.0 || *rv != 0.0)
406 fk425pv (dtheta, dphi, ptheta, pphi, px, rv);
407 else
408 fk425m (dtheta, dphi, ptheta, pphi);
409 if (ep2 != 2000.0) {
410 *dtheta = *dtheta + ((ep2 - 2000.0) * *ptheta);
411 *dphi = *dphi + ((ep2 - 2000.0) * *pphi);
412 }
413 }
414 else if (ep2 > 0.0)
415 fk425e (dtheta, dphi, ep2);
416 else
417 fk425 (dtheta, dphi);
418 }
419 else if (sys1 == WCS_GALACTIC)
420 gal2fk5 (dtheta, dphi);
421 else if (sys1 == WCS_ECLIPTIC)
422 ecl2fk5 (dtheta, dphi, ep2);
423 }
424
425 else if (sys2 == WCS_GALACTIC) {
426 if (sys1 == WCS_B1950) {
427 if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
428 *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
429 *dphi = *dphi + (*pphi * (ep2 - ep1));
430 }
431 fk42gal (dtheta, dphi);
432 }
433 else if (sys1 == WCS_J2000) {
434 if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
435 *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
436 *dphi = *dphi + (*pphi * (ep2 - ep1));
437 }
438 fk52gal (dtheta, dphi);
439 }
440 else if (sys1 == WCS_ECLIPTIC) {
441 ecl2fk5 (dtheta, dphi, ep2);
442 fk52gal (dtheta, dphi);
443 }
444 }
445
446 else if (sys2 == WCS_ECLIPTIC) {
447 if (sys1 == WCS_B1950) {
448 if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
449 *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
450 *dphi = *dphi + (*pphi * (ep2 - ep1));
451 }
452 if (ep2 > 0.0)
453 fk42ecl (dtheta, dphi, ep2);
454 else
455 fk42ecl (dtheta, dphi, 1950.0);
456 }
457 else if (sys1 == WCS_J2000) {
458 if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
459 *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
460 *dphi = *dphi + (*pphi * (ep2 - ep1));
461 }
462 fk52ecl (dtheta, dphi, ep2);
463 }
464 else if (sys1 == WCS_GALACTIC) {
465 gal2fk5 (dtheta, dphi);
466 fk52ecl (dtheta, dphi, ep2);
467 }
468 }
469
470 /* Precess to desired equinox, if necessary */
471 if (eq1 != eq2) {
472 if (sys2 == WCS_B1950 && eq2 != 1950.0)
473 fk4prec (1950.0, eq2, dtheta, dphi);
474 if (sys2 == WCS_J2000 && eq2 != 2000.0)
475 fk5prec (2000.0, eq2, dtheta, dphi);
476 }
477
478 /* Keep latitude/declination between +90 and -90 degrees */
479 if (*dphi > 90.0) {
480 *dphi = 180.0 - *dphi;
481 *dtheta = *dtheta + 180.0;
482 }
483 else if (*dphi < -90.0) {
484 *dphi = -180.0 - *dphi;
485 *dtheta = *dtheta + 180.0;
486 }
487
488 /* Keep longitude/right ascension between 0 and 360 degrees */
489 if (*dtheta > 360.0)
490 *dtheta = *dtheta - 360.0;
491 else if (*dtheta < 0.0)
492 *dtheta = *dtheta + 360.0;
493 return;
494 }
495
496
497 /* Convert from coordinate system sys1 to coordinate system sys2 */
498
499 void
wcscon(sys1,sys2,eq1,eq2,dtheta,dphi,epoch)500 wcscon (sys1, sys2, eq1, eq2, dtheta, dphi, epoch)
501
502 int sys1; /* Input coordinate system (J2000, B1950, ECLIPTIC, GALACTIC */
503 int sys2; /* Output coordinate system (J2000, B1950, ECLIPTIC, GALACTIC */
504 double eq1; /* Input equinox (default of sys1 if 0.0) */
505 double eq2; /* Output equinox (default of sys2 if 0.0) */
506 double *dtheta; /* Longitude or right ascension in degrees
507 Input in sys1, returned in sys2 */
508 double *dphi; /* Latitude or declination in degrees
509 Input in sys1, returned in sys2 */
510 double epoch; /* Besselian epoch in years */
511
512 {
513 void fk5prec(), fk4prec();
514
515 /* Set equinoxes if 0.0 */
516 if (eq1 == 0.0) {
517 if (sys1 == WCS_B1950)
518 eq1 = 1950.0;
519 else
520 eq1 = 2000.0;
521 }
522 if (eq2 == 0.0) {
523 if (sys2 == WCS_B1950)
524 eq2 = 1950.0;
525 else
526 eq2 = 2000.0;
527 }
528
529 /* Set systems and equinoxes so that ICRS coordinates are not precessed */
530 if (sys1 == WCS_ICRS && sys2 == WCS_ICRS)
531 eq2 = eq1;
532
533 if (sys1 == WCS_J2000 && sys2 == WCS_ICRS && eq1 == 2000.0) {
534 eq2 = eq1;
535 sys1 = sys2;
536 }
537
538 if (sys1 == WCS_ICRS && sys2 == WCS_J2000 && eq2 == 2000.0) {
539 eq1 = eq2;
540 sys1 = sys2;
541 }
542
543 /* If systems and equinoxes are the same, return */
544 if (sys2 == sys1 && eq1 == eq2)
545 return;
546
547 /* Precess from input equinox, if necessary */
548 if (eq1 != eq2) {
549 if (sys1 == WCS_B1950 && eq1 != 1950.0)
550 fk4prec (eq1, 1950.0, dtheta, dphi);
551 if (sys1 == WCS_J2000 && eq1 != 2000.0)
552 fk5prec (eq1, 2000.0, dtheta, dphi);
553 }
554
555 /* Convert to B1950 FK4 */
556 if (sys2 == WCS_B1950) {
557 if (sys1 == WCS_J2000) {
558 if (epoch > 0)
559 fk524e (dtheta, dphi, epoch);
560 else
561 fk524 (dtheta, dphi);
562 }
563 else if (sys1 == WCS_GALACTIC)
564 gal2fk4 (dtheta, dphi);
565 else if (sys1 == WCS_ECLIPTIC) {
566 if (epoch > 0)
567 ecl2fk4 (dtheta, dphi, epoch);
568 else
569 ecl2fk4 (dtheta, dphi, 1950.0);
570 }
571 }
572
573 else if (sys2 == WCS_J2000) {
574 if (sys1 == WCS_B1950) {
575 if (epoch > 0)
576 fk425e (dtheta, dphi, epoch);
577 else
578 fk425 (dtheta, dphi);
579 }
580 else if (sys1 == WCS_GALACTIC)
581 gal2fk5 (dtheta, dphi);
582 else if (sys1 == WCS_ECLIPTIC) {
583 if (epoch > 0)
584 ecl2fk5 (dtheta, dphi, epoch);
585 else
586 ecl2fk5 (dtheta, dphi, 2000.0);
587 }
588 }
589
590 else if (sys2 == WCS_GALACTIC) {
591 if (sys1 == WCS_B1950)
592 fk42gal (dtheta, dphi);
593 else if (sys1 == WCS_J2000)
594 fk52gal (dtheta, dphi);
595 else if (sys1 == WCS_ECLIPTIC) {
596 if (epoch > 0)
597 ecl2fk5 (dtheta, dphi, epoch);
598 else
599 ecl2fk5 (dtheta, dphi, 2000.0);
600 fk52gal (dtheta, dphi);
601 }
602 }
603
604 else if (sys2 == WCS_ECLIPTIC) {
605 if (sys1 == WCS_B1950) {
606 if (epoch > 0)
607 fk42ecl (dtheta, dphi, epoch);
608 else
609 fk42ecl (dtheta, dphi, 1950.0);
610 }
611 else if (sys1 == WCS_J2000) {
612 if (epoch > 0)
613 fk52ecl (dtheta, dphi, epoch);
614 else
615 fk52ecl (dtheta, dphi, 2000.0);
616 }
617 else if (sys1 == WCS_GALACTIC) {
618 gal2fk5 (dtheta, dphi);
619 if (epoch > 0)
620 fk52ecl (dtheta, dphi, epoch);
621 else
622 fk52ecl (dtheta, dphi, 2000.0);
623 }
624 }
625
626 /* Precess to desired equinox, if necessary */
627 if (eq1 != eq2) {
628 if (sys2 == WCS_B1950 && eq2 != 1950.0)
629 fk4prec (1950.0, eq2, dtheta, dphi);
630 if (sys2 == WCS_J2000 && eq2 != 2000.0)
631 fk5prec (2000.0, eq2, dtheta, dphi);
632 }
633
634 /* Keep latitude/declination between +90 and -90 degrees */
635 if (*dphi > 90.0) {
636 *dphi = 180.0 - *dphi;
637 *dtheta = *dtheta + 180.0;
638 }
639 else if (*dphi < -90.0) {
640 *dphi = -180.0 - *dphi;
641 *dtheta = *dtheta + 180.0;
642 }
643
644 /* Keep longitude/right ascension between 0 and 360 degrees */
645 if (*dtheta > 360.0)
646 *dtheta = *dtheta - 360.0;
647 else if (*dtheta < 0.0)
648 *dtheta = *dtheta + 360.0;
649
650 return;
651 }
652
653
654 /* Set coordinate system from string */
655 int
wcscsys(wcstring)656 wcscsys (wcstring)
657
658 char *wcstring; /* Name of coordinate system */
659 {
660 double equinox;
661
662 if (wcstring[0] == 'J' || wcstring[0] == 'j' ||
663 !strcmp (wcstring,"2000") || !strcmp (wcstring, "2000.0") ||
664 !strcmp (wcstring,"ICRS") || !strcmp (wcstring, "icrs") ||
665 !strncmp (wcstring,"FK5",3) || !strncmp (wcstring, "fk5",3))
666 return WCS_J2000;
667
668 if (wcstring[0] == 'B' || wcstring[0] == 'b' ||
669 !strcmp (wcstring,"1950") || !strcmp (wcstring, "1950.0") ||
670 !strncmp (wcstring,"FK4",3) || !strncmp (wcstring, "fk4",3))
671 return WCS_B1950;
672
673 else if (wcstring[0] == 'I' || wcstring[0] == 'i' )
674 return WCS_ICRS;
675
676 else if (wcstring[0] == 'G' || wcstring[0] == 'g' )
677 return WCS_GALACTIC;
678
679 else if (wcstring[0] == 'E' || wcstring[0] == 'e' )
680 return WCS_ECLIPTIC;
681
682 else if (wcstring[0] == 'A' || wcstring[0] == 'a' )
683 return WCS_ALTAZ;
684
685 else if (wcstring[0] == 'N' || wcstring[0] == 'n' )
686 return WCS_NPOLE;
687
688 else if (wcstring[0] == 'L' || wcstring[0] == 'l' )
689 return WCS_LINEAR;
690
691 else if (!strncasecmp (wcstring, "pixel", 5))
692 return WCS_XY;
693
694 else if (wcstring[0] == 'P' || wcstring[0] == 'p' )
695 return WCS_PLANET;
696
697 else if (isnum (wcstring) == 1 || isnum (wcstring) == 2) {
698 equinox = atof (wcstring);
699 if (equinox > 1980.0)
700 return WCS_J2000;
701 else if (equinox > 1900.0)
702 return WCS_B1950;
703 else
704 return -1;
705 }
706 else
707 return -1;
708 }
709
710
711 /* Set equinox from string (return 0.0 if not obvious) */
712
713 double
wcsceq(wcstring)714 wcsceq (wcstring)
715
716 char *wcstring; /* Name of coordinate system */
717 {
718 if (wcstring[0] == 'J' || wcstring[0] == 'j' ||
719 wcstring[0] == 'B' || wcstring[0] == 'b')
720 return (atof (wcstring+1));
721 else if (!strncmp (wcstring, "FK4",3) ||
722 !strncmp (wcstring, "fk4",3))
723 return (1950.0);
724 else if (!strncmp (wcstring, "FK5",3) ||
725 !strncmp (wcstring, "fk5",3))
726 return (2000.0);
727 else if (!strncmp (wcstring, "ICRS",4) ||
728 !strncmp (wcstring, "icrs",4))
729 return (2000.0);
730 else if (wcstring[0] == '1' || wcstring[0] == '2')
731 return (atof (wcstring));
732 else
733 return (0.0);
734 }
735
736
737 /* Set coordinate system type string from system and equinox */
738
739 void
wcscstr(cstr,syswcs,equinox,epoch)740 wcscstr (cstr, syswcs, equinox, epoch)
741
742 char *cstr; /* Coordinate system string (returned) */
743 int syswcs; /* Coordinate system code */
744 double equinox; /* Equinox of coordinate system */
745 double epoch; /* Epoch of coordinate system */
746 {
747
748 char *estr;
749
750 if (syswcs == WCS_XY) {
751 strcpy (cstr, "XY");
752 return;
753 }
754
755 /* Try to figure out coordinate system if it is not set */
756 if (epoch == 0.0)
757 epoch = equinox;
758 if (syswcs < 0) {
759 if (equinox > 0.0) {
760 if (equinox == 2000.0)
761 syswcs = WCS_J2000;
762 else if (equinox == 1950.0)
763 syswcs = WCS_B1950;
764 }
765 else if (epoch > 0.0) {
766 if (epoch > 1980.0) {
767 syswcs = WCS_J2000;
768 equinox = 2000.0;
769 }
770 else {
771 syswcs = WCS_B1950;
772 equinox = 1950.0;
773 }
774 }
775 else
776 syswcs = WCS_J2000;
777 }
778
779 /* Set coordinate system string from system flag and epoch */
780 if (syswcs == WCS_B1950) {
781 if (epoch == 1950.0 || epoch == 0.0)
782 strcpy (cstr, "B1950");
783 else
784 sprintf (cstr, "B%7.2f", equinox);
785 if ((estr = strsrch (cstr,".00")) != NULL) {
786 estr[0] = (char) 0;
787 estr[1] = (char) 0;
788 estr[2] = (char) 0;
789 }
790 }
791 else if (syswcs == WCS_GALACTIC)
792 strcpy (cstr, "galactic");
793 else if (syswcs == WCS_ECLIPTIC)
794 strcpy (cstr, "ecliptic");
795 else if (syswcs == WCS_J2000) {
796 if (epoch == 2000.0 || epoch == 0.0)
797 strcpy (cstr, "J2000");
798 else
799 sprintf (cstr, "J%7.2f", equinox);
800 if ((estr = strsrch (cstr,".00")) != NULL) {
801 estr[0] = (char) 0;
802 estr[1] = (char) 0;
803 estr[2] = (char) 0;
804 }
805 }
806 else if (syswcs == WCS_ICRS) {
807 strcpy (cstr, "ICRS");
808 }
809 else if (syswcs == WCS_PLANET) {
810 strcpy (cstr, "PLANET");
811 }
812 else if (syswcs == WCS_LINEAR || syswcs == WCS_XY) {
813 strcpy (cstr, "LINEAR");
814 }
815 return;
816 }
817
818
819 /* Constant vector and matrix (by columns)
820 These values were obtained by inverting C.Hohenkerk's forward matrix
821 (private communication), which agrees with the one given in reference
822 2 but which has one additional decimal place. */
823
824 static double a[3] = {-1.62557e-6, -0.31919e-6, -0.13843e-6};
825 static double ad[3] = {1.245e-3, -1.580e-3, -0.659e-3};
826 static double d2pi = 6.283185307179586476925287; /* two PI */
827 static double tiny = 1.e-30; /* small number to avoid arithmetic problems */
828
829 /* FK524 convert J2000 FK5 star data to B1950 FK4
830 based on Starlink sla_fk524 by P.T.Wallace 27 October 1987 */
831
832 static double emi[6][6] = {
833 { 0.9999256795, /* emi[0][0] */
834 0.0111814828, /* emi[0][1] */
835 0.0048590039, /* emi[0][2] */
836 -0.00000242389840, /* emi[0][3] */
837 -0.00000002710544, /* emi[0][4] */
838 -0.00000001177742 }, /* emi[0][5] */
839
840 { -0.0111814828, /* emi[1][0] */
841 0.9999374849, /* emi[1][1] */
842 -0.0000271771, /* emi[1][2] */
843 0.00000002710544, /* emi[1][3] */
844 -0.00000242392702, /* emi[1][4] */
845 0.00000000006585 }, /* emi[1][5] */
846
847 { -0.0048590040, /* emi[2][0] */
848 -0.0000271557, /* emi[2][1] */
849 0.9999881946, /* emi[2][2] */
850 0.00000001177742, /* emi[2][3] */
851 0.00000000006585, /* emi[2][4] */
852 -0.00000242404995 }, /* emi[2][5] */
853
854 { -0.000551, /* emi[3][0] */
855 0.238509, /* emi[3][1] */
856 -0.435614, /* emi[3][2] */
857 0.99990432, /* emi[3][3] */
858 0.01118145, /* emi[3][4] */
859 0.00485852 }, /* emi[3][5] */
860
861 { -0.238560, /* emi[4][0] */
862 -0.002667, /* emi[4][1] */
863 0.012254, /* emi[4][2] */
864 -0.01118145, /* emi[4][3] */
865 0.99991613, /* emi[4][4] */
866 -0.00002717 }, /* emi[4][5] */
867
868 { 0.435730, /* emi[5][0] */
869 -0.008541, /* emi[5][1] */
870 0.002117, /* emi[5][2] */
871 -0.00485852, /* emi[5][3] */
872 -0.00002716, /* emi[5][4] */
873 0.99996684 } /* emi[5][5] */
874 };
875
876 void
fk524(ra,dec)877 fk524 (ra,dec)
878
879 double *ra; /* Right ascension in degrees (J2000 in, B1950 out) */
880 double *dec; /* Declination in degrees (J2000 in, B1950 out) */
881
882 {
883 double rapm; /* Proper motion in right ascension */
884 double decpm; /* Proper motion in declination */
885 /* In: deg/jul.yr. Out: deg/trop.yr. */
886
887 rapm = (double) 0.0;
888 decpm = (double) 0.0;
889 fk524m (ra, dec, &rapm, &decpm);
890 return;
891 }
892
893 void
fk524e(ra,dec,epoch)894 fk524e (ra, dec, epoch)
895
896 double *ra; /* Right ascension in degrees (J2000 in, B1950 out) */
897 double *dec; /* Declination in degrees (J2000 in, B1950 out) */
898 double epoch; /* Besselian epoch in years */
899
900 {
901 double rapm; /* Proper motion in right ascension */
902 double decpm; /* Proper motion in declination */
903 /* In: deg/jul.yr. Out: deg/trop.yr. */
904
905 rapm = (double) 0.0;
906 decpm = (double) 0.0;
907 fk524m (ra, dec, &rapm, &decpm);
908 *ra = *ra + (rapm * (epoch - 1950.0));
909 *dec = *dec + (decpm * (epoch - 1950.0));
910 return;
911 }
912
913 void
fk524m(ra,dec,rapm,decpm)914 fk524m (ra,dec,rapm,decpm)
915
916 double *ra; /* Right ascension in degrees (J2000 in, B1950 out) */
917 double *dec; /* Declination in degrees (J2000 in, B1950 out) */
918 double *rapm; /* Proper motion in right ascension */
919 double *decpm; /* Proper motion in declination */
920 /* In: ra/dec deg/jul.yr. Out: ra/dec deg/trop.yr. */
921
922 {
923 double parallax = 0.0;
924 double rv = 0.0;
925
926 fk524pv (ra, dec, rapm, decpm, ¶llax, &rv);
927 return;
928 }
929
930
931 void
fk524pv(ra,dec,rapm,decpm,parallax,rv)932 fk524pv (ra,dec,rapm,decpm, parallax, rv)
933
934 double *ra; /* Right ascension in degrees (J2000 in, B1950 out) */
935 double *dec; /* Declination in degrees (J2000 in, B1950 out) */
936 double *rapm; /* Proper motion in right ascension */
937 double *decpm; /* Proper motion in declination
938 * In: ra/dec degrees/Julian year (not ra*cos(dec))
939 * Out: ra/dec degrees/tropical year */
940 double *parallax; /* Parallax (arcsec) */
941 double *rv; /* Rradial velocity (km/s, +ve = moving away) */
942
943 /* This routine converts stars from the IAU 1976 FK5 Fricke
944 system, to the old Bessel-Newcomb FK4 system, using Yallop's
945 implementation (see ref 2) of a matrix method due to Standish
946 (see ref 3). The numerical values of ref 2 are used canonically.
947
948 * Conversion from other than Julian epoch 2000.0 to other than Besselian
949 epoch 1950.0 will require use of the appropriate precession, proper
950 motion, and e-terms routines before and/or after fk524 is called.
951
952 * In the FK4 catalogue the proper motions of stars within 10 degrees
953 of the poles do not embody the differential e-term effect and should,
954 strictly speaking, be handled in a different manner from stars outside
955 these regions. however, given the general lack of homogeneity of the
956 star data available for routine astrometry, the difficulties of handling
957 positions that may have been determined from astrometric fields spanning
958 the polar and non-polar regions, the likelihood that the differential
959 e-terms effect was not taken into account when allowing for proper motion
960 in past astrometry, and the undesirability of a discontinuity in the
961 algorithm, the decision has been made in this routine to include the
962 effect of differential e-terms on the proper motions for all stars,
963 whether polar or not, at epoch 2000, and measuring on the sky rather
964 than in terms of dra, the errors resulting from this simplification are
965 less than 1 milliarcsecond in position and 1 milliarcsecond per century
966 in proper motion.
967
968 References:
969
970 1 "Mean and apparent place computations in the new IAU System.
971 I. The transformation of astrometric catalog systems to the
972 equinox J2000.0." Smith, C.A.; Kaplan, G.H.; Hughes, J.A.;
973 Seidelmann, P.K.; Yallop, B.D.; Hohenkerk, C.Y.
974 Astronomical Journal vol. 97, Jan. 1989, p. 265-273.
975
976 2 "Mean and apparent place computations in the new IAU System.
977 II. Transformation of mean star places from FK4 B1950.0 to
978 FK5 J2000.0 using matrices in 6-space." Yallop, B.D.;
979 Hohenkerk, C.Y.; Smith, C.A.; Kaplan, G.H.; Hughes, J.A.;
980 Seidelmann, P.K.; Astronomical Journal vol. 97, Jan. 1989,
981 p. 274-279.
982
983 3 Seidelmann, P.K. (ed), 1992. "Explanatory Supplement to
984 the Astronomical Almanac", ISBN 0-935702-68-7.
985
986 4 "Conversion of positions and proper motions from B1950.0 to the
987 IAU system at J2000.0", Standish, E.M. Astronomy and
988 Astrophysics, vol. 115, no. 1, Nov. 1982, p. 20-22.
989
990 P.T.Wallace Starlink 19 December 1993
991 Doug Mink Smithsonian Astrophysical Observatory 1 November 2000 */
992
993 {
994 double r2000,d2000; /* J2000.0 ra,dec (radians) */
995 double r1950,d1950; /* B1950.0 ra,dec (rad) */
996
997 /* Miscellaneous */
998 double ur,ud;
999 double sr, cr, sd, cd, x, y, z, w, wd;
1000 double v1[6],v2[6];
1001 double xd,yd,zd;
1002 double rxyz, rxysq, rxy;
1003 double dra,ddec;
1004 int i,j;
1005 int diag = 0;
1006
1007 /* Constants */
1008 double zero = (double) 0.0;
1009 double vf = 21.095; /* Km per sec to AU per tropical century */
1010 /* = 86400 * 36524.2198782 / 149597870 */
1011
1012 /* Convert J2000 RA and Dec from degrees to radians */
1013 r2000 = degrad (*ra);
1014 d2000 = degrad (*dec);
1015
1016 /* Convert J2000 RA and Dec proper motion from degrees/year to arcsec/tc */
1017 ur = *rapm * 360000.0;
1018 ud = *decpm * 360000.0;
1019
1020 /* Spherical to Cartesian */
1021 sr = sin (r2000);
1022 cr = cos (r2000);
1023 sd = sin (d2000);
1024 cd = cos (d2000);
1025
1026 x = cr * cd;
1027 y = sr * cd;
1028 z = sd;
1029
1030 v1[0] = x;
1031 v1[1] = y;
1032 v1[2] = z;
1033
1034 if (ur != zero || ud != zero) {
1035 v1[3] = -(ur*y) - (cr*sd*ud);
1036 v1[4] = (ur*x) - (sr*sd*ud);
1037 v1[5] = (cd*ud);
1038 }
1039 else {
1040 v1[3] = zero;
1041 v1[4] = zero;
1042 v1[5] = zero;
1043 }
1044
1045 /* Convert position + velocity vector to bn system */
1046 for (i = 0; i < 6; i++) {
1047 w = zero;
1048 for (j = 0; j < 6; j++) {
1049 w = w + emi[i][j] * v1[j];
1050 }
1051 v2[i] = w;
1052 }
1053
1054 /* Vector components */
1055 x = v2[0];
1056 y = v2[1];
1057 z = v2[2];
1058
1059 /* Magnitude of position vector */
1060 rxyz = sqrt (x*x + y*y + z*z);
1061
1062 /* Apply e-terms to position */
1063 w = (x * a[0]) + (y * a[1]) + (z * a[2]);
1064 x = x + (a[0] * rxyz) - (w * x);
1065 y = y + (a[1] * rxyz) - (w * y);
1066 z = z + (a[2] * rxyz) - (w * z);
1067
1068 /* Recompute magnitude of position vector */
1069 rxyz = sqrt (x*x + y*y + z*z);
1070
1071 /* Apply e-terms to position and velocity */
1072 x = v2[0];
1073 y = v2[1];
1074 z = v2[2];
1075 w = (x * a[0]) + (y * a[1]) + (z * a[2]);
1076 wd = (x * ad[0]) + (y * ad[1]) + (z * ad[2]);
1077 x = x + (a[0] * rxyz) - (w * x);
1078 y = y + (a[1] * rxyz) - (w * y);
1079 z = z + (a[2] * rxyz) - (w * z);
1080 xd = v2[3] + (ad[0] * rxyz) - (wd * x);
1081 yd = v2[4] + (ad[1] * rxyz) - (wd * y);
1082 zd = v2[5] + (ad[2] * rxyz) - (wd * z);
1083
1084 /* Convert to spherical */
1085 rxysq = (x * x) + (y * y);
1086 rxy = sqrt (rxysq);
1087
1088 /* Convert back to spherical coordinates */
1089 if (x == zero && y == zero)
1090 r1950 = zero;
1091 else {
1092 r1950 = atan2 (y,x);
1093 if (r1950 < zero)
1094 r1950 = r1950 + d2pi;
1095 }
1096 d1950 = atan2 (z,rxy);
1097
1098 if (rxy > tiny) {
1099 ur = (x*yd - y*xd) / rxysq;
1100 ud = (zd*rxysq - z * (x*xd + y*yd)) / ((rxysq + z*z) * rxy);
1101 }
1102
1103 if (*parallax > tiny) {
1104 *rv = ((x * xd) + (y * yd) + (z * zd)) / (*parallax * vf * rxyz);
1105 *parallax = *parallax / rxyz;
1106 }
1107
1108 /* Return results */
1109 *ra = raddeg (r1950);
1110 *dec = raddeg (d1950);
1111 *rapm = ur / 360000.0;
1112 *decpm = ud / 360000.0;
1113
1114 if (diag) {
1115 dra = 240.0 * raddeg (r1950 - r2000);
1116 ddec = 3600.0 * raddeg (d1950 - d2000);
1117 fprintf(stderr,"B1950-J2000: dra= %11.5f sec ddec= %f11.5f arcsec\n",
1118 dra, ddec);
1119 }
1120
1121 return;
1122 }
1123
1124
1125 /* Convert B1950.0 FK4 star data to J2000.0 FK5 */
1126 static double em[6][6] = {
1127 { 0.9999256782, /* em[0][0] */
1128 -0.0111820611, /* em[0][1] */
1129 -0.0048579477, /* em[0][2] */
1130 0.00000242395018, /* em[0][3] */
1131 -0.00000002710663, /* em[0][4] */
1132 -0.00000001177656 }, /* em[0][5] */
1133
1134 { 0.0111820610, /* em[1][0] */
1135 0.9999374784, /* em[1][1] */
1136 -0.0000271765, /* em[1][2] */
1137 0.00000002710663, /* em[1][3] */
1138 0.00000242397878, /* em[1][4] */
1139 -0.00000000006587 }, /* em[1][5] */
1140
1141 { 0.0048579479, /* em[2][0] */
1142 -0.0000271474, /* em[2][1] */
1143 0.9999881997, /* em[2][2] */
1144 0.00000001177656, /* em[2][3] */
1145 -0.00000000006582, /* em[2][4] */
1146 0.00000242410173 }, /* em[2][5] */
1147
1148 { -0.000551, /* em[3][0] */
1149 -0.238565, /* em[3][1] */
1150 0.435739, /* em[3][2] */
1151 0.99994704, /* em[3][3] */
1152 -0.01118251, /* em[3][4] */
1153 -0.00485767 }, /* em[3][5] */
1154
1155 { 0.238514, /* em[4][0] */
1156 -0.002667, /* em[4][1] */
1157 -0.008541, /* em[4][2] */
1158 0.01118251, /* em[4][3] */
1159 0.99995883, /* em[4][4] */
1160 -0.00002718 }, /* em[4][5] */
1161
1162 { -0.435623, /* em[5][0] */
1163 0.012254, /* em[5][1] */
1164 0.002117, /* em[5][2] */
1165 0.00485767, /* em[5][3] */
1166 -0.00002714, /* em[5][4] */
1167 1.00000956 } /* em[5][5] */
1168 };
1169
1170 void
fk425(ra,dec)1171 fk425 (ra, dec)
1172
1173 double *ra; /* Right ascension in degrees (B1950 in, J2000 out) */
1174 double *dec; /* Declination in degrees (B1950 in, J2000 out) */
1175
1176 {
1177 double rapm; /* Proper motion in right ascension */
1178 double decpm; /* Proper motion in declination */
1179 /* In: rad/trop.yr. Out: rad/jul.yr. */
1180
1181 rapm = (double) 0.0;
1182 decpm = (double) 0.0;
1183 fk425m (ra, dec, &rapm, &decpm);
1184 return;
1185 }
1186
1187
1188 void
fk425e(ra,dec,epoch)1189 fk425e (ra, dec, epoch)
1190
1191 double *ra; /* Right ascension in degrees (B1950 in, J2000 out) */
1192 double *dec; /* Declination in degrees (B1950 in, J2000 out) */
1193 double epoch; /* Besselian epoch in years */
1194 {
1195 double rapm; /* Proper motion in right ascension */
1196 double decpm; /* Proper motion in declination */
1197 /* In: rad/trop.yr. Out: rad/jul.yr. */
1198
1199 rapm = (double) 0.0;
1200 decpm = (double) 0.0;
1201 fk425m (ra, dec, &rapm, &decpm);
1202 *ra = *ra + (rapm * (epoch - 2000.0));
1203 *dec = *dec + (decpm * (epoch - 2000.0));
1204 return;
1205 }
1206
1207 void
fk425m(ra,dec,rapm,decpm)1208 fk425m (ra, dec, rapm, decpm)
1209
1210 double *ra, *dec; /* Right ascension and declination in degrees
1211 input: B1950.0,FK4 returned: J2000.0,FK5 */
1212 double *rapm, *decpm; /* Proper motion in right ascension and declination
1213 input: B1950.0,FK4 returned: J2000.0,FK5
1214 ra/dec deg/trop.yr. ra/dec deg/jul.yr. */
1215 {
1216 double parallax = 0.0;
1217 double rv = 0.0;
1218
1219 fk425pv (ra, dec, rapm, decpm, ¶llax, &rv);
1220 return;
1221 }
1222
1223
1224 void
fk425pv(ra,dec,rapm,decpm,parallax,rv)1225 fk425pv (ra,dec,rapm,decpm, parallax, rv)
1226
1227 double *ra; /* Right ascension in degrees (J2000 in, B1950 out) */
1228 double *dec; /* Declination in degrees (J2000 in, B1950 out) */
1229 double *rapm; /* Proper motion in right ascension */
1230 double *decpm; /* Proper motion in declination
1231 * In: ra/dec degrees/Julian year (not ra*cos(dec))
1232 * Out: ra/dec degrees/tropical year */
1233 double *parallax; /* Parallax (arcsec) */
1234 double *rv; /* Rradial velocity (km/s, +ve = moving away) */
1235
1236 /* This routine converts stars from the old Bessel-Newcomb FK4 system
1237 to the IAU 1976 FK5 Fricke system, using Yallop's implementation
1238 (see ref 2) of a matrix method due to Standish (see ref 3). The
1239 numerical values of ref 2 are used canonically.
1240
1241 * Conversion from other than Besselian epoch 1950.0 to other than Julian
1242 epoch 2000.0 will require use of the appropriate precession, proper
1243 motion, and e-terms routines before and/or after fk425 is called.
1244
1245 * In the FK4 catalogue the proper motions of stars within 10 degrees
1246 of the poles do not embody the differential e-term effect and should,
1247 strictly speaking, be handled in a different manner from stars outside
1248 these regions. however, given the general lack of homogeneity of the
1249 star data available for routine astrometry, the difficulties of handling
1250 positions that may have been determined from astrometric fields spanning
1251 the polar and non-polar regions, the likelihood that the differential
1252 e-terms effect was not taken into account when allowing for proper motion
1253 in past astrometry, and the undesirability of a discontinuity in the
1254 algorithm, the decision has been made in this routine to include the
1255 effect of differential e-terms on the proper motions for all stars,
1256 whether polar or not, at epoch 2000, and measuring on the sky rather
1257 than in terms of dra, the errors resulting from this simplification are
1258 less than 1 milliarcsecond in position and 1 milliarcsecond per century
1259 in proper motion.
1260
1261 References:
1262
1263 1 "Mean and apparent place computations in the new IAU System.
1264 I. The transformation of astrometric catalog systems to the
1265 equinox J2000.0." Smith, C.A.; Kaplan, G.H.; Hughes, J.A.;
1266 Seidelmann, P.K.; Yallop, B.D.; Hohenkerk, C.Y.
1267 Astronomical Journal vol. 97, Jan. 1989, p. 265-273.
1268
1269 2 "Mean and apparent place computations in the new IAU System.
1270 II. Transformation of mean star places from FK4 B1950.0 to
1271 FK5 J2000.0 using matrices in 6-space." Yallop, B.D.;
1272 Hohenkerk, C.Y.; Smith, C.A.; Kaplan, G.H.; Hughes, J.A.;
1273 Seidelmann, P.K.; Astronomical Journal vol. 97, Jan. 1989,
1274 p. 274-279.
1275
1276 3 "Conversion of positions and proper motions from B1950.0 to the
1277 IAU system at J2000.0", Standish, E.M. Astronomy and
1278 Astrophysics, vol. 115, no. 1, Nov. 1982, p. 20-22.
1279
1280 P.T.Wallace Starlink 20 December 1993
1281 Doug Mink Smithsonian Astrophysical Observatory 7 June 1995 */
1282
1283 {
1284 double r1950,d1950; /* B1950.0 ra,dec (rad) */
1285 double r2000,d2000; /* J2000.0 ra,dec (rad) */
1286
1287 /* Miscellaneous */
1288 double ur,ud,sr,cr,sd,cd,w,wd;
1289 double x,y,z,xd,yd,zd, dra,ddec;
1290 double rxyz, rxysq, rxy, rxyzsq, spxy, spxyz;
1291 int i,j;
1292 int diag = 0;
1293
1294 double r0[3],rd0[3]; /* star position and velocity vectors */
1295 double v1[6],v2[6]; /* combined position and velocity vectors */
1296
1297 /* Constants */
1298 double zero = (double) 0.0;
1299 double vf = 21.095; /* Km per sec to AU per tropical century */
1300 /* = 86400 * 36524.2198782 / 149597870 */
1301
1302 /* Convert B1950 RA and Dec from degrees to radians */
1303 r1950 = degrad (*ra);
1304 d1950 = degrad (*dec);
1305
1306 /* Convert B1950 RA and Dec proper motion from degrees/year to arcsec/tc */
1307 ur = *rapm * 360000.0;
1308 ud = *decpm * 360000.0;
1309
1310 /* Convert direction to Cartesian */
1311 sr = sin (r1950);
1312 cr = cos (r1950);
1313 sd = sin (d1950);
1314 cd = cos (d1950);
1315 r0[0] = cr * cd;
1316 r0[1] = sr * cd;
1317 r0[2] = sd;
1318
1319 /* Convert motion to Cartesian */
1320 w = vf * *rv * *parallax;
1321 if (ur != zero || ud != zero || (*rv != zero && *parallax != zero)) {
1322 rd0[0] = (-sr * cd * ur) - (cr * sd * ud) + (w * r0[0]);
1323 rd0[1] = (cr * cd * ur) - (sr * sd * ud) + (w * r0[1]);
1324 rd0[2] = (cd * ud) + (w * r0[2]);
1325 }
1326 else {
1327 rd0[0] = zero;
1328 rd0[1] = zero;
1329 rd0[2] = zero;
1330 }
1331
1332 /* Remove e-terms from position and express as position+velocity 6-vector */
1333 w = (r0[0] * a[0]) + (r0[1] * a[1]) + (r0[2] * a[2]);
1334 for (i = 0; i < 3; i++)
1335 v1[i] = r0[i] - a[i] + (w * r0[i]);
1336
1337 /* Remove e-terms from proper motion and express as 6-vector */
1338 wd = (r0[0] * ad[0]) + (r0[1] * ad[1]) + (r0[2] * ad[2]);
1339 for (i = 0; i < 3; i++)
1340 v1[i+3] = rd0[i] - ad[i] + (wd * r0[i]);
1341
1342 /* Alternately: Put proper motion in 6-vector without adding e-terms
1343 for (i = 0; i < 3; i++)
1344 v1[i+3] = rd0[i]; */
1345
1346 /* Convert position + velocity vector to FK5 system */
1347 for (i = 0; i < 6; i++) {
1348 w = zero;
1349 for (j = 0; j < 6; j++) {
1350 w += em[i][j] * v1[j];
1351 }
1352 v2[i] = w;
1353 }
1354
1355 /* Vector components */
1356 x = v2[0];
1357 y = v2[1];
1358 z = v2[2];
1359 xd = v2[3];
1360 yd = v2[4];
1361 zd = v2[5];
1362
1363 /* Magnitude of position vector */
1364 rxysq = x*x + y*y;
1365 rxy = sqrt (rxysq);
1366 rxyzsq = rxysq + z*z;
1367 rxyz = sqrt (rxyzsq);
1368
1369 spxy = (x * xd) + (y * yd);
1370 spxyz = spxy + (z * zd);
1371
1372 /* Convert back to spherical coordinates */
1373 if (x == zero && y == zero)
1374 r2000 = zero;
1375 else {
1376 r2000 = atan2 (y,x);
1377 if (r2000 < zero)
1378 r2000 = r2000 + d2pi;
1379 }
1380 d2000 = atan2 (z,rxy);
1381
1382 if (rxy > tiny) {
1383 ur = ((x * yd) - (y * xd)) / rxysq;
1384 ud = ((zd * rxysq) - (z * spxy)) / (rxyzsq * rxy);
1385 }
1386
1387 if (*parallax > tiny) {
1388 *rv = spxyz / (*parallax * rxyz * vf);
1389 *parallax = *parallax / rxyz;
1390 }
1391
1392 /* Return results */
1393 *ra = raddeg (r2000);
1394 *dec = raddeg (d2000);
1395 *rapm = ur / 360000.0;
1396 *decpm = ud / 360000.0;
1397
1398 if (diag) {
1399 dra = 240.0 * raddeg (r2000 - r1950);
1400 ddec = 3600.0 * raddeg (d2000 - d1950);
1401 fprintf(stderr,"J2000-B1950: dra= %11.5f sec ddec= %f11.5f arcsec\n",
1402 dra, ddec);
1403 }
1404 return;
1405 }
1406
1407 int idg=0;
1408
1409 /* l2,b2 system of galactic coordinates
1410 * p = 192.25 ra of galactic north pole (mean b1950.0)
1411 * q = 62.6 inclination of galactic to mean b1950.0 equator
1412 * r = 33 longitude of ascending node
1413 * p,q,r are degrees
1414
1415 * Equatorial to galactic rotation matrix
1416 (The Eulerian angles are p, q, 90-r)
1417 +cp.cq.sr-sp.cr +sp.cq.sr+cp.cr -sq.sr
1418 -cp.cq.cr-sp.sr -sp.cq.cr+cp.sr +sq.cr
1419 cp.sq +sp.sq +cq
1420 */
1421
1422 static
1423 double bgal[3][3] =
1424 {{-0.066988739415,-0.872755765852,-0.483538914632},
1425 {0.492728466075,-0.450346958020, 0.744584633283},
1426 {-0.867600811151,-0.188374601723, 0.460199784784}};
1427
1428 /*--- Transform B1950.0 FK4 equatorial coordinates to
1429 * IAU 1958 galactic coordinates */
1430
1431 void
fk42gal(dtheta,dphi)1432 fk42gal (dtheta,dphi)
1433
1434 double *dtheta; /* B1950.0 FK4 right ascension in degrees
1435 Galactic longitude (l2) in degrees (returned) */
1436 double *dphi; /* B1950.0 FK4 declination in degrees
1437 Galactic latitude (b2) in degrees (returned) */
1438
1439 /* Input equatorial coordinates are B1950 FK4.
1440 Use fk52gal() to convert from j2000.0 coordinates.
1441 Reference: Blaauw et al, MNRAS,121,123 (1960) */
1442 {
1443 double pos[3],pos1[3],r,dl,db,rl,rb,rra,rdec,dra,ddec;
1444 void v2s3(),s2v3();
1445 int i;
1446 char *eqcoor, *eqstrn();
1447
1448 dra = *dtheta;
1449 ddec = *dphi;
1450 rra = degrad (dra);
1451 rdec = degrad (ddec);
1452
1453 /* remove e-terms */
1454 /* call jpabe (rra,rdec,-1,idg) */
1455
1456 /* Spherical to Cartesian */
1457 r = 1.;
1458 s2v3 (rra,rdec,r,pos);
1459
1460 /* rotate to galactic */
1461 for (i = 0; i<3; i++) {
1462 pos1[i] = pos[0]*bgal[i][0] + pos[1]*bgal[i][1] + pos[2]*bgal[i][2];
1463 }
1464
1465 /* Cartesian to spherical */
1466 v2s3 (pos1,&rl,&rb,&r);
1467
1468 dl = raddeg (rl);
1469 db = raddeg (rb);
1470 *dtheta = dl;
1471 *dphi = db;
1472
1473 /* Print result if in diagnostic mode */
1474 if (idg) {
1475 eqcoor = eqstrn (dra,ddec);
1476 fprintf (stderr,"FK42GAL: B1950 RA,Dec= %s\n",eqcoor);
1477 fprintf (stderr,"FK42GAL: long = %.5f lat = %.5f\n",dl,db);
1478 free (eqcoor);
1479 }
1480
1481 return;
1482 }
1483
1484
1485 /*--- Transform IAU 1958 galactic coordinates to B1950.0 'FK4'
1486 * equatorial coordinates */
1487
1488 void
gal2fk4(dtheta,dphi)1489 gal2fk4 (dtheta,dphi)
1490
1491 double *dtheta; /* Galactic longitude (l2) in degrees
1492 B1950 FK4 RA in degrees (returned) */
1493 double *dphi; /* Galactic latitude (b2) in degrees
1494 B1950 FK4 Dec in degrees (returned) */
1495
1496 /* Output equatorial coordinates are B1950.0 FK4.
1497 Use gal2fk5() to convert to J2000 coordinates.
1498 Reference: Blaauw et al, MNRAS,121,123 (1960) */
1499
1500 {
1501 double pos[3],pos1[3],r,dl,db,rl,rb,rra,rdec,dra,ddec;
1502 void v2s3(),s2v3();
1503 char *eqcoor, *eqstrn();
1504 int i;
1505
1506 /* spherical to cartesian */
1507 dl = *dtheta;
1508 db = *dphi;
1509 rl = degrad (dl);
1510 rb = degrad (db);
1511 r = 1.0;
1512 s2v3 (rl,rb,r,pos);
1513
1514 /* rotate to equatorial coordinates */
1515 for (i = 0; i < 3; i++) {
1516 pos1[i] = pos[0]*bgal[0][i] + pos[1]*bgal[1][i] + pos[2]*bgal[2][i];
1517 }
1518
1519 /* cartesian to spherical */
1520 v2s3 (pos1,&rra,&rdec,&r);
1521
1522 /* introduce e-terms */
1523 /* jpabe (rra,rdec,-1,idg); */
1524
1525 dra = raddeg (rra);
1526 ddec = raddeg (rdec);
1527 *dtheta = dra;
1528 *dphi = ddec;
1529
1530 /* print result if in diagnostic mode */
1531 if (idg) {
1532 fprintf (stderr,"GAL2FK4: long = %.5f lat = %.5f\n",dl,db);
1533 eqcoor = eqstrn (dra,ddec);
1534 fprintf (stderr,"GAL2FK4: B1950 RA,Dec= %s\n",eqcoor);
1535 free (eqcoor);
1536 }
1537
1538 return;
1539 }
1540
1541
1542 /* l2,b2 system of galactic coordinates
1543 p = 192.25 ra of galactic north pole (mean b1950.0)
1544 q = 62.6 inclination of galactic to mean b1950.0 equator
1545 r = 33 longitude of ascending node
1546 p,q,r are degrees */
1547
1548 /* Equatorial to galactic rotation matrix
1549 The eulerian angles are p, q, 90-r
1550 +cp.cq.sr-sp.cr +sp.cq.sr+cp.cr -sq.sr
1551 -cp.cq.cr-sp.sr -sp.cq.cr+cp.sr +sq.cr
1552 +cp.sq +sp.sq +cq */
1553
1554 static
1555 double jgal[3][3] =
1556 {{-0.054875539726,-0.873437108010,-0.483834985808},
1557 {0.494109453312,-0.444829589425, 0.746982251810},
1558 {-0.867666135858,-0.198076386122, 0.455983795705}};
1559
1560 /* Transform J2000 equatorial coordinates to IAU 1958 galactic coordinates */
1561
1562 void
fk52gal(dtheta,dphi)1563 fk52gal (dtheta,dphi)
1564
1565 double *dtheta; /* J2000 right ascension in degrees
1566 Galactic longitude (l2) in degrees (returned) */
1567 double *dphi; /* J2000 declination in degrees
1568 Galactic latitude (b2) in degrees (returned) */
1569
1570 /* Rotation matrices by P.T.Wallace, Starlink eqgal and galeq, March 1986 */
1571
1572 /* Input equatorial coordinates are J2000 FK5.
1573 Use gal2fk4() if converting from B1950 FK4 coordinates.
1574 Reference: Blaauw et al, MNRAS,121,123 (1960) */
1575 {
1576 double pos[3],pos1[3],r,dl,db,rl,rb,rra,rdec,dra,ddec;
1577 void v2s3(),s2v3();
1578 char *eqcoor, *eqstrn();
1579 int i;
1580
1581 /* Spherical to cartesian */
1582 dra = *dtheta;
1583 ddec = *dphi;
1584 rra = degrad (dra);
1585 rdec = degrad (ddec);
1586 r = 1.0;
1587 (void)s2v3 (rra,rdec,r,pos);
1588
1589 /* Rotate to galactic */
1590 for (i = 0; i < 3; i++) {
1591 pos1[i] = pos[0]*jgal[i][0] + pos[1]*jgal[i][1] + pos[2]*jgal[i][2];
1592 }
1593
1594 /* Cartesian to spherical */
1595 v2s3 (pos1,&rl,&rb,&r);
1596
1597 dl = raddeg (rl);
1598 db = raddeg (rb);
1599 *dtheta = dl;
1600 *dphi = db;
1601
1602 /* Print result if in diagnostic mode */
1603 if (idg) {
1604 eqcoor = eqstrn (dra,ddec);
1605 fprintf (stderr,"FK52GAL: J2000 RA,Dec= %s\n",eqcoor);
1606 fprintf (stderr,"FK52GAL: long = %.5f lat = %.5f\n",dl,db);
1607 free (eqcoor);
1608 }
1609
1610 return;
1611 }
1612
1613
1614 /*--- Transform IAU 1958 galactic coordinates to J2000 equatorial coordinates */
1615
1616 void
gal2fk5(dtheta,dphi)1617 gal2fk5 (dtheta,dphi)
1618
1619 double *dtheta; /* Galactic longitude (l2) in degrees
1620 J2000.0 ra in degrees (returned) */
1621 double *dphi; /* Galactic latitude (b2) in degrees
1622 J2000.0 dec in degrees (returned) */
1623
1624 /* Output equatorial coordinates are J2000.
1625 Use gal2fk4() to convert to B1950 coordinates.
1626 Reference: Blaauw et al, MNRAS,121,123 (1960) */
1627
1628 {
1629 double pos[3],pos1[3],r,dl,db,rl,rb,rra,rdec,dra,ddec;
1630 void v2s3(),s2v3();
1631 int i;
1632 char *eqcoor, *eqstrn();
1633
1634 /* Spherical to Cartesian */
1635 dl = *dtheta;
1636 db = *dphi;
1637 rl = degrad (dl);
1638 rb = degrad (db);
1639 r = 1.0;
1640 s2v3 (rl,rb,r,pos);
1641
1642 /* Rotate to equatorial coordinates */
1643 for (i = 0; i < 3; i++) {
1644 pos1[i] = pos[0]*jgal[0][i] + pos[1]*jgal[1][i] + pos[2]*jgal[2][i];
1645 }
1646
1647 /* Cartesian to Spherical */
1648 v2s3 (pos1,&rra,&rdec,&r);
1649 dra = raddeg (rra);
1650 ddec = raddeg (rdec);
1651 *dtheta = dra;
1652 *dphi = ddec;
1653
1654 /* Print result if in diagnostic mode */
1655 if (idg) {
1656 fprintf (stderr,"GAL2FK5: long = %.5f lat = %.5f\n",dl,db);
1657 eqcoor = eqstrn (dra,ddec);
1658 fprintf (stderr,"GAL2FK5: J2000 RA,Dec= %s\n",eqcoor);
1659 free (eqcoor);
1660 }
1661
1662 return;
1663 }
1664
1665
1666 /* Return string with right ascension in hours and declination in degrees */
1667
eqstrn(dra,ddec)1668 char *eqstrn (dra, ddec)
1669
1670 double dra; /* Right ascension in degrees */
1671 double ddec; /* Declination in degrees */
1672
1673 {
1674 char *eqcoor; /* ASCII character string of position (returned) */
1675 char decp;
1676 int rah,irm,decd,decm;
1677 double xpos,ypos,xp,yp,ras,decs;
1678
1679 /* Right ascension to hours, minutes, and seconds */
1680 xpos = dra / 15.0;
1681 rah = (int) xpos;
1682 xp = (double) 60.0 * (xpos - (double) rah);
1683 irm = (int) xp;
1684 ras = (double) 60.0 * (xp - (double) irm);
1685
1686 /* Declination to degrees, minutes, seconds */
1687 if (ddec < 0) {
1688 ypos = -ddec;
1689 decp = '-';
1690 }
1691 else {
1692 decp = '+';
1693 ypos = ddec;
1694 }
1695 decd = (int) ypos;
1696 yp = (double) 60.0 * (ypos - (double) decd);
1697 decm = (int) yp;
1698 decs = (double) 60.0 * (yp - (double) decm);
1699
1700 eqcoor = malloc (32);
1701 (void)sprintf (eqcoor,"%02d:%02d:%06.3f %c%02d:%02d:%05.2f",
1702 rah,irm,ras,decp,decd,decm,decs);
1703 if (eqcoor[6] == ' ')
1704 eqcoor[6] = '0';
1705 if (eqcoor[20] == ' ')
1706 eqcoor[20] = '0';
1707
1708 return (eqcoor);
1709 }
1710
1711
1712 /* Convert geocentric equatorial rectangular coordinates to
1713 right ascension and declination, and distance */
1714
1715
1716 /* These routines are based on similar ones in Pat Wallace's slalib package */
1717
1718 /* Convert B1950 right ascension and declination to ecliptic coordinates */
1719
1720 void
fk42ecl(dtheta,dphi,epoch)1721 fk42ecl (dtheta, dphi, epoch)
1722
1723 double *dtheta; /* B1950 right ascension in degrees
1724 Galactic longitude (l2) in degrees (returned) */
1725 double *dphi; /* B1950 declination in degrees
1726 Galactic latitude (b2) in degrees (returned) */
1727 double epoch; /* Besselian epoch in years */
1728
1729 {
1730 void fk425e(), fk52ecl();
1731
1732 /* Convert from B1950 to J2000 coordinates */
1733 fk425e (dtheta, dphi, epoch);
1734
1735 /* Convert from J2000 to ecliptic coordinates */
1736 fk52ecl (dtheta, dphi, epoch);
1737
1738 return;
1739 }
1740
1741 /* Convert J2000 right ascension and declination to ecliptic coordinates */
1742
1743 void
fk52ecl(dtheta,dphi,epoch)1744 fk52ecl (dtheta, dphi, epoch)
1745
1746 double *dtheta; /* J2000 right ascension in degrees
1747 Galactic longitude (l2) in degrees (returned) */
1748 double *dphi; /* J2000 declination in degrees
1749 Galactic latitude (b2) in degrees (returned) */
1750 double epoch; /* Besselian epoch in years */
1751
1752 {
1753 int i, j;
1754 double t, eps0, rphi, rtheta;
1755 double v1[3], v2[3], r;
1756 double rmat[9], *rmati; /* Rotation matrix */
1757
1758 void rotmat(), v2s3(), s2v3(), fk5prec();
1759
1760 /* Precess coordinates from J2000 to epoch */
1761 if (epoch != 2000.0)
1762 fk5prec (2000.0, epoch, dtheta, dphi);
1763
1764 /* Convert from degrees to radians */
1765 rtheta = degrad (*dtheta);
1766 rphi = degrad (*dphi);
1767
1768 /* Convert RA,Dec to x,y,z */
1769 r = 1.0;
1770 s2v3 (rtheta, rphi, r, v1);
1771
1772 /* Interval between basic epoch J2000.0 and current epoch (JC) in centuries*/
1773 t = (epoch - 2000.0) * 0.01;
1774
1775 /* Mean obliquity */
1776 eps0 = secrad ((84381.448 + (-46.8150 + (-0.00059 + 0.001813*t) * t) * t));
1777
1778 /* Form the equatorial to ecliptic rotation matrix (IAU 1980 theory).
1779 * References: Murray, C.A., Vectorial Astrometry, section 4.3.
1780 * The matrix is in the sense v[ecl] = rmat * v[equ]; the
1781 * equator, equinox and ecliptic are mean of date. */
1782 rotmat (1, eps0, 0.0, 0.0, rmat);
1783
1784 /* Multiply position vector by equatoria to eccliptic rotation matrix */
1785 rmati = rmat;
1786 for (i = 0; i < 3; i++) {
1787 v2[i] = 0;
1788 for (j = 0; j < 3; j++)
1789 v2[i] = v2[i] + (*rmati++ * v1[j]);
1790 }
1791
1792 /* Convert x,y,z to latitude, longitude */
1793 v2s3 (v2, &rtheta, &rphi, &r);
1794
1795 /* Convert from radians to degrees */
1796 *dtheta = raddeg (rtheta);
1797 *dphi = raddeg (rphi);
1798 }
1799
1800
1801 /* Convert ecliptic coordinates to B1950 right ascension and declination */
1802
1803 void
ecl2fk4(dtheta,dphi,epoch)1804 ecl2fk4 (dtheta, dphi, epoch)
1805
1806 double *dtheta; /* Galactic longitude (l2) in degrees
1807 B1950 right ascension in degrees (returned) */
1808 double *dphi; /* Galactic latitude (b2) in degrees
1809 B1950 declination in degrees (returned) */
1810 double epoch; /* Besselian epoch in years */
1811
1812 {
1813 void ecl2fk5(), fk524e();
1814
1815 /* Convert from ecliptic to J2000 coordinates */
1816 ecl2fk5 (dtheta, dphi, epoch);
1817
1818 /* Convert from J2000 to B1950 coordinates */
1819 fk524e (dtheta, dphi, epoch);
1820
1821 return;
1822 }
1823
1824
1825
1826 /* Convert ecliptic coordinates to J2000 right ascension and declination */
1827
1828 void
ecl2fk5(dtheta,dphi,epoch)1829 ecl2fk5 (dtheta, dphi, epoch)
1830
1831 double *dtheta; /* Galactic longitude (l2) in degrees
1832 J2000 right ascension in degrees (returned) */
1833 double *dphi; /* Galactic latitude (b2) in degrees
1834 J2000 declination in degrees (returned) */
1835 double epoch; /* Besselian epoch in years */
1836
1837 {
1838 int i, j;
1839 double rtheta, rphi, v1[3], v2[3];
1840 double t, eps0, r;
1841 double rmat[9]; /* Rotation matrix */
1842 void v2s3(),s2v3(), fk5prec(), rotmat();
1843
1844 rtheta = degrad (*dtheta);
1845 rphi = degrad (*dphi);
1846
1847 /* Convert RA,Dec to x,y,z */
1848 r = 1.0;
1849 s2v3 (rtheta, rphi, r, v1);
1850
1851 /* Interval between basic epoch J2000.0 and current epoch (JC) in centuries*/
1852 t = (epoch - 2000.0) * 0.01;
1853
1854 /* Mean obliquity */
1855 eps0 = secrad ((84381.448 + (-46.8150 + (-0.00059 + 0.001813*t) * t) * t));
1856
1857 /* Form the equatorial to ecliptic rotation matrix (IAU 1980 theory).
1858 * References: Murray, C.A., Vectorial Astrometry, section 4.3.
1859 * The matrix is in the sense v[ecl] = rmat * v[equ]; the
1860 * equator, equinox and ecliptic are mean of date. */
1861 rotmat (1, eps0, 0.0, 0.0, rmat);
1862
1863 /* Multiply position vector by ecliptic to equatorial rotation matrix */
1864 for (i = 0; i < 3; i++) {
1865 v2[i] = 0;
1866 for (j = 0; j < 3; j++)
1867 v2[i] = v2[i] + (rmat[3*j + i] * v1[j]);
1868 }
1869
1870 /* Cartesian to spherical */
1871 v2s3 (v2, &rtheta, &rphi, &r);
1872
1873 /* Convert from radians to degrees */
1874 *dtheta = raddeg (rtheta);
1875 *dphi = raddeg (rphi);
1876
1877 if (epoch != 2000.0)
1878 fk5prec (epoch, 2000.0, dtheta, dphi);
1879 }
1880
1881
1882 /* The following routines are modified from Patrick Wallace's SLALIB */
1883
1884 /* Precess coordinates between epochs in FK4 */
1885 void
fk4prec(ep0,ep1,ra,dec)1886 fk4prec (ep0, ep1, ra, dec)
1887
1888 double ep0; /* Starting Besselian epoch */
1889 double ep1; /* Ending Besselian epoch */
1890 double *ra; /* RA in degrees mean equator & equinox of epoch ep0
1891 mean equator & equinox of epoch ep1 (returned) */
1892 double *dec; /* Dec in degrees mean equator & equinox of epoch ep0
1893 mean equator & equinox of epoch ep1 (returned) */
1894 /*
1895 ** Precession - FK4 (Bessel-Newcomb, pre-IAU1976)
1896 **
1897 ** This routine will not correctly convert between FK4 and FK5
1898 ** For output in FK5, precess to 1950.0 and use fk425() on result.
1899 **
1900 ** Based on slaPreces(), P.T.Wallace Starlink 22 December 1993
1901 */
1902 {
1903 int i, j;
1904 double pm[9], *pmi, v1[3], v2[3], rra, rdec, r;
1905 void v2s3(),s2v3(), mprecfk4();
1906
1907 rra = degrad (*ra);
1908 rdec = degrad (*dec);
1909 r = 1.0;
1910
1911 /* Generate appropriate precession matrix */
1912 mprecfk4 ( ep0, ep1, pm );
1913
1914 /* Convert RA,Dec to x,y,z */
1915 s2v3 (rra, rdec, r, v1);
1916
1917 /* Multiply position vector by precession matrix */
1918 pmi = pm;
1919 for (i = 0; i < 3; i++) {
1920 v2[i] = 0;
1921 for (j = 0; j < 3; j++)
1922 v2[i] = v2[i] + (*pmi++ * v1[j]);
1923 }
1924
1925 /* Back to RA,Dec */
1926 v2s3 (v2, &rra, &rdec, &r);
1927
1928 /* Convert from radians to degrees */
1929 *ra = raddeg (rra);
1930 *dec = raddeg (rdec);
1931 }
1932
1933 void
fk5prec(ep0,ep1,ra,dec)1934 fk5prec (ep0, ep1, ra, dec)
1935
1936 double ep0; /* Starting epoch */
1937 double ep1; /* Ending epoch */
1938 double *ra; /* RA in degrees mean equator & equinox of epoch ep0
1939 mean equator & equinox of epoch ep1 (returned) */
1940 double *dec; /* Dec in degrees mean equator & equinox of epoch ep0
1941 mean equator & equinox of epoch ep1 (returned) */
1942 /*
1943 ** Precession - FK5 (Fricke, post-IAU1976)
1944 **
1945 ** This routine will not correctly convert between FK5 and FK4.
1946 ** For output in FK4, precess to 2000.0 and use fk524() on result.
1947 **
1948 ** Based on slaPreces(), P.T.Wallace Starlink 22 December 1993
1949 */
1950 {
1951 int i, j;
1952 double pm[9], *pmi, v1[3], v2[3], rra, rdec, r;
1953 void v2s3(),s2v3(), mprecfk5();
1954
1955 rra = degrad (*ra);
1956 rdec = degrad (*dec);
1957 r = 1.0;
1958
1959 /* Generate appropriate precession matrix */
1960 mprecfk5 (ep0, ep1, pm);
1961
1962 /* Convert RA,Dec to x,y,z */
1963 s2v3 (rra, rdec, r, v1);
1964
1965 /* Multiply position vector by precession matrix */
1966 pmi = pm;
1967 for (i = 0; i < 3; i++) {
1968 v2[i] = 0;
1969 for (j = 0; j < 3; j++)
1970 v2[i] = v2[i] + ( v1[j] * *pmi++ );
1971 }
1972
1973 /* Back to RA,Dec */
1974 v2s3 (v2, &rra, &rdec, &r);
1975
1976 /* Convert from radians to degrees */
1977 *ra = raddeg (rra);
1978 *dec = raddeg (rdec);
1979 return;
1980 }
1981
1982
1983 void
mprecfk4(bep0,bep1,rmatp)1984 mprecfk4 (bep0, bep1, rmatp)
1985
1986 double bep0; /* Beginning Besselian epoch */
1987 double bep1; /* Ending Besselian epoch */
1988 double rmatp[9]; /* 3x3 Precession matrix (returned) */
1989
1990 /*
1991 ** Generate the matrix of precession between two epochs,
1992 ** using the old, pre-IAU1976, Bessel-Newcomb model, using
1993 ** Kinoshita's formulation (double precision)
1994 **
1995 ** The matrix is in the sense v(bep1) = rmatp * v(bep0)
1996 **
1997 ** Reference:
1998 ** Kinoshita, H. (1975) 'Formulas for precession', SAO Special
1999 ** Report No. 364, Smithsonian Institution Astrophysical
2000 ** Observatory, Cambridge, Massachusetts.
2001 **
2002 ** Based on slaPrebn() by P.T.Wallace Starlink 30 October 1993
2003 */
2004 {
2005 double bigt, t, tas2r, w, zeta, z, theta;
2006 void rotmat();
2007
2008 /* Interval between basic epoch B1850.0 and beginning epoch in TC */
2009 bigt = ( bep0 - 1850.0 ) / 100.0;
2010
2011 /* Interval over which precession required, in tropical centuries */
2012 t = ( bep1 - bep0 ) / 100.0;
2013
2014 /* Euler angles */
2015 tas2r = secrad (t);
2016 w = 2303.5548 + ( 1.39720 + 0.000059 * bigt ) * bigt;
2017 zeta = (w + ( 0.30242 - 0.000269 * bigt + 0.017996 * t ) * t ) * tas2r;
2018 z = (w + ( 1.09478 + 0.000387 * bigt + 0.018324 * t ) * t ) * tas2r;
2019 theta = ( 2005.1125 + ( - 0.85294 - 0.000365* bigt ) * bigt +
2020 ( - 0.42647 - 0.000365 * bigt - 0.041802 * t ) * t ) * tas2r;
2021
2022 /* Rotation matrix */
2023 rotmat (323, -zeta, theta, -z, rmatp);
2024 return;
2025 }
2026
2027
2028 void
mprecfk5(ep0,ep1,rmatp)2029 mprecfk5 (ep0, ep1, rmatp)
2030
2031 double ep0; /* Beginning epoch */
2032 double ep1; /* Ending epoch */
2033 double rmatp[9]; /* 3x3 Precession matrix (returned) */
2034
2035 /*
2036 ** Form the matrix of precession between two epochs (IAU 1976, FK5).
2037 ** Notes:
2038 ** 1) The epochs are TDB (loosely ET) Julian epochs.
2039 ** 2) The matrix is in the sense v(ep1) = rmatp * v(ep0) .
2040 **
2041 ** References:
2042 ** Lieske,J.H., 1979. Astron. Astrophys.,73,282.
2043 ** equations (6) & (7), p283.
2044 ** Kaplan,G.H., 1981. USNO circular no. 163, pa2.
2045 **
2046 ** Based on slaPrec(), P.T.Wallace Starlink 31 October 1993
2047 */
2048 {
2049 double t0, t, tas2r, w, zeta, z, theta;
2050 void rotmat();
2051
2052 /* Interval between basic epoch J2000.0 and beginning epoch (JC) */
2053 t0 = ( ep0 - 2000.0 ) / 100.0;
2054
2055 /* Interval over which precession required (JC) */
2056 t = ( ep1 - ep0 ) / 100.0;
2057
2058 /* Euler angles */
2059 tas2r = secrad (t);
2060 w = 2306.2181 + ( ( 1.39656 - ( 0.000139 * t0 ) ) * t0 );
2061 zeta = (w + ( ( 0.30188 - 0.000344 * t0 ) + 0.017998 * t ) * t ) * tas2r;
2062 z = (w + ( ( 1.09468 + 0.000066 * t0 ) + 0.018203 * t ) * t ) * tas2r;
2063 theta = ( ( 2004.3109 + ( - 0.85330 - 0.000217 * t0 ) * t0 )
2064 + ( ( -0.42665 - 0.000217 * t0 ) - 0.041833 * t ) * t ) * tas2r;
2065
2066 /* Rotation matrix */
2067 rotmat (323, -zeta, theta, -z, rmatp);
2068 return;
2069 }
2070
2071
2072 /* Make 3-D rotation matrix from up to three rotations */
2073
2074 void
rotmat(axes,rot1,rot2,rot3,matrix)2075 rotmat (axes, rot1, rot2, rot3, matrix)
2076
2077 int axes; /* Axes about which coordinates are rotated (1=x, 2=y, 3=z) */
2078 double rot1; /* First rotation in degrees */
2079 double rot2; /* Second rotation in degrees */
2080 double rot3; /* Third rotation in degrees */
2081 double *matrix; /* 3x3 rotation matrix (returned) */
2082
2083 {
2084 int i, j, k, naxis, iaxes, iaxis;
2085 double rot[3], srot, crot, *mati, w, wm[9], *wmi, matn[9];
2086 int axis[3];
2087
2088 /* Initial final rotation matrix */
2089 mati = matrix;
2090 for (i = 0; i < 3; i++) {
2091 for (j=0; j < 3; j++) {
2092 if (i == j)
2093 *mati++ = 1.0;
2094 else
2095 *mati++ = 0.0;
2096 }
2097 }
2098
2099 /* Separate digits of rotation axis string and count rotations */
2100 naxis = 0;
2101 iaxes = axes;
2102 axis[0] = iaxes / 100;
2103 if (axis[0] > 0) {
2104 naxis++;
2105 iaxes = iaxes - (100 * axis[0]);
2106 }
2107 axis[naxis] = iaxes / 10;
2108 if (axis[naxis] > 0) {
2109 iaxes = iaxes - (10 * axis[naxis]);
2110 naxis++;
2111 }
2112 axis[naxis] = iaxes;
2113 if (axis[naxis] > 0)
2114 naxis++;
2115
2116 /* Set up rotation angles */
2117 rot[0] = rot1;
2118 rot[1] = rot2;
2119 rot[2] = rot3;
2120
2121 /* For each digit of axis string, set up matrix */
2122 for (iaxis = 0; iaxis < naxis; iaxis++) {
2123
2124 /* Initialize current rotation matrix */
2125 mati = matn;
2126 for (i = 0; i < 3; i++) {
2127 for (j=0; j < 3; j++) {
2128 if (i == j)
2129 *mati++ = 1.0;
2130 else
2131 *mati++ = 0.0;
2132 }
2133 }
2134
2135 srot = sin (rot[iaxis]);
2136 crot = cos (rot[iaxis]);
2137
2138 /* Matrix for rotation in X */
2139 if (axis[iaxis] == 1) {
2140 matn[4] = crot;
2141 matn[5] = srot;
2142 matn[7] = -srot;
2143 matn[8] = crot;
2144 }
2145
2146 /* Matrix for rotation in Y */
2147 else if (axis[iaxis] == 2) {
2148 matn[0] = crot;
2149 matn[2] = -srot;
2150 matn[6] = srot;
2151 matn[8] = crot;
2152 }
2153
2154 /* Matrix for rotation in Z */
2155 else {
2156 matn[0] = crot;
2157 matn[1] = srot;
2158 matn[3] = -srot;
2159 matn[4] = crot;
2160 }
2161
2162 /* Multiply existing rotation matrix by new rotation matrix */
2163 for (i = 0; i < 3; i++) {
2164 for (j = 0; j < 3; j++) {
2165 w = 0.0;
2166 for (k = 0; k < 3; k++)
2167 w+= matn[3*i + k] * matrix[3*k + j];
2168 wm[3*i + j] = w;
2169 }
2170 }
2171
2172 /* Update output matrix */
2173 mati = matrix;
2174 wmi = wm;
2175 for (i = 0; i < 9; i++) {
2176 *mati++ = *wmi++;
2177 }
2178 }
2179 return;
2180 }
2181
2182
2183 /* The following routines are from Doug Mink's Fortran ephemeris library */
2184
2185 /* Convert right ascensiona and declination in degrees and distance to
2186 geocentric equatorial rectangular coordinates */
2187
2188 void
d2v3(rra,rdec,r,pos)2189 d2v3 (rra,rdec,r,pos)
2190
2191 double rra; /* Right ascension in degrees */
2192 double rdec; /* Declination in degrees */
2193 double r; /* Distance to object in same units as pos */
2194 double pos[3]; /* x,y,z geocentric equatorial position of object (returned) */
2195 {
2196 s2v3 (degrad (rra), degrad (rdec), r, pos);
2197
2198 return;
2199 }
2200
2201
2202 /* Convert right ascension, declination, and distance to
2203 geocentric equatorial rectangular coordinates */
2204
2205 void
s2v3(rra,rdec,r,pos)2206 s2v3 (rra,rdec,r,pos)
2207
2208 double rra; /* Right ascension in radians */
2209 double rdec; /* Declination in radians */
2210 double r; /* Distance to object in same units as pos */
2211 double pos[3]; /* x,y,z geocentric equatorial position of object (returned) */
2212 {
2213 pos[0] = r * cos (rra) * cos (rdec);
2214 pos[1] = r * sin (rra) * cos (rdec);
2215 pos[2] = r * sin (rdec);
2216
2217 return;
2218 }
2219
2220
2221 /* Convert geocentric equatorial rectangular coordinates to
2222 right ascension and declination in degrees and distance */
2223
2224 void
v2d3(pos,rra,rdec,r)2225 v2d3 (pos,rra,rdec,r)
2226
2227 double pos[3]; /* x,y,z geocentric equatorial position of object */
2228 double *rra; /* Right ascension in degrees (returned) */
2229 double *rdec; /* Declination in degrees (returned) */
2230 double *r; /* Distance to object in same units as pos (returned) */
2231 {
2232 v2s3 (pos, rra, rdec, r);
2233 *rra = raddeg (*rra);
2234 *rdec = raddeg (*rdec);
2235 return;
2236 }
2237
2238 /* Convert geocentric equatorial rectangular coordinates to
2239 right ascension, declination, and distance */
2240
2241 void
v2s3(pos,rra,rdec,r)2242 v2s3 (pos,rra,rdec,r)
2243
2244 double pos[3]; /* x,y,z geocentric equatorial position of object */
2245 double *rra; /* Right ascension in radians (returned) */
2246 double *rdec; /* Declination in radians (returned) */
2247 double *r; /* Distance to object in same units as pos (returned) */
2248 {
2249 double x,y,z,rxy,rxy2,z2;
2250
2251 x = pos[0];
2252 y = pos[1];
2253 z = pos[2];
2254
2255 *rra = atan2 (y, x);
2256
2257 /* Keep RA within 0 to 2pi range */
2258 if (*rra < 0.0)
2259 *rra = *rra + (2.0 * PI);
2260 if (*rra > 2.0 * PI)
2261 *rra = *rra - (2.0 * PI);
2262
2263 rxy2 = x*x + y*y;
2264 rxy = sqrt (rxy2);
2265 *rdec = atan2 (z, rxy);
2266
2267 z2 = z * z;
2268 *r = sqrt (rxy2 + z2);
2269
2270 return;
2271 }
2272
2273 /*
2274 * Nov 6 1995 Include stdlib.h instead of malloc.h
2275 * Apr 1 1996 Add arbitrary epoch precession
2276 * Apr 26 1996 Add FK4 <-> FK5 subroutines for use when epoch is known
2277 * Aug 6 1996 Clean up after lint
2278 * Nov 4 1996 Break SLA subroutines into separate file slasubs.c
2279 * Dec 9 1996 Change arguments to degrees in FK4 and FK5 precession programs
2280 * Dec 10 1996 All subroutine arguments are degrees except vector conversions
2281 *
2282 * Mar 20 1997 Drop unused variables after lint
2283 *
2284 * Apr 14 1998 Add ecliptic coordinate conversions and general conversion routines
2285 * Apr 23 1998 Add LINEAR coordinate system
2286 * Apr 28 1998 Change coordinate system flags to WCS_*
2287 * Apr 28 1998 Return -1 from wcscsys if not a legal coordinate system
2288 * May 7 1998 Keep theta within 0 to 2pi in ecl2fk5()
2289 * May 13 1998 Add wcsceq()
2290 * May 13 1998 Add equinox arguments to wcscon()
2291 * Jun 24 1998 Set J2000 from ICRS in wcscsys()
2292 * Jul 9 1998 Include stdio.h for fprintf() and sprintf() declarations
2293 * Sep 17 1998 Add wcscstr() to get coordinate string
2294 * Sep 21 1998 Fix bug in wcscstr() which returned B2000 instead of J2000
2295 * Sep 21 1998 Add subroutine to convert proper motions, too.
2296 * Oct 21 1998 In wcscstr(), drop .00 from returned string
2297 * Nov 18 1998 Rename jpcop() v2s3() and jpcon() s2v3() (spherical to vector)
2298 * Dec 2 1998 Add PLANET coordinate system to wcscsys() and wcscstr()
2299 *
2300 * Mar 10 2000 Precess coordinates correctly from other than 1950.0 and 2000.0
2301 * Mar 10 2000 Set coordinate system to J2000 or B1950 if string is numeric
2302 * Mar 14 2000 Clean up code in fk524m() and fk425m()
2303 * May 31 2000 Add proper motion correctly if proper motion precessed
2304 * Jun 26 2000 Add some support for WCS_XY image coordinates
2305 * Sep 14 2000 Return -1 from wcscsys if equinox is less than 1900.0
2306 * Oct 31 2000 Add proper motion after fk425 or fk524 from system epoch
2307 * Oct 31 2000 Fix proper motion units in fk524p() and fk425p()
2308 * Nov 6 2000 Update fk425 and fk524 algorithms to include parallax and rv
2309 *
2310 * Jan 11 2001 Print all messages to stderr
2311 * Mar 21 2001 Move braces around bgal[] and jgal[] matrix initialization
2312 *
2313 * Feb 13 2002 Fix precession units problem in ecl2fk5() and fk52ecl()
2314 *
2315 * Apr 13 2005 Replace all sla_lib calls with local code
2316 * Nov 1 2005 Add WCS_ICRS, and unprecessable system
2317 *
2318 * Jan 5 2006 Fix bugs in precession subroutines mprecxxx()
2319 * May 3 2006 Drop declarations of unused variables suggested by Robert Lupton
2320 * Oct 6 2006 If pixel coordinates, set system to WCS_XY in wcscsys()
2321 * Oct 30 2006 Add LINEAR and ICRS to wcscstr() returns
2322 *
2323 * Aug 15 2007 Clean up code in rotmat()
2324 * Nov 8 2007 In wcsconp, make it clear that proper motion is in spherical coordinates
2325 *
2326 * Mar 29 2010 Fix bug in computing the magnitude of the e-terms in fk524()
2327 * Mar 30 2010 Drop ep1 assignment after line 178 in wcsconp()
2328 *
2329 * Jun 9 2016 Fix isnum() tests for added coloned times and dashed dates
2330 */
2331