summaryrefslogtreecommitdiffstats
path: root/funtools/wcs/wcscon.c
diff options
context:
space:
mode:
Diffstat (limited to 'funtools/wcs/wcscon.c')
-rw-r--r--funtools/wcs/wcscon.c2328
1 files changed, 2328 insertions, 0 deletions
diff --git a/funtools/wcs/wcscon.c b/funtools/wcs/wcscon.c
new file mode 100644
index 0000000..6e99bd3
--- /dev/null
+++ b/funtools/wcs/wcscon.c
@@ -0,0 +1,2328 @@
+/*** File wcscon.c
+ *** March 30, 2010
+ *** Doug Mink, Harvard-Smithsonian Center for Astrophysics
+ *** Some subroutines are based on Starlink subroutines by Patrick Wallace
+ *** Copyright (C) 1995-2010
+ *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2 of the License, or (at your option) any later version.
+
+ This library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with this library; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+ Correspondence concerning WCSTools should be addressed as follows:
+ Internet email: jmink@cfa.harvard.edu
+ Postal address: Jessica Mink
+ Smithsonian Astrophysical Observatory
+ 60 Garden St.
+ Cambridge, MA 02138 USA
+
+ * Module: wcscon.c (World Coordinate System conversion)
+ * Purpose: Convert between various sky coordinate systems
+ * Subroutine: wcscon (sys1,sys2,eq1,eq2,theta,phi,epoch)
+ * convert between coordinate systems
+ * Subroutine: wcsconp (sys1,sys2,eq1,eq2,ep1,ep2,dtheta,dphi,ptheta,pphi)
+ * convert coordinates and proper motion between coordinate systems
+ * Subroutine: wcsconv (sys1,sys2,eq1,eq2,ep1,ep2,dtheta,dphi,ptheta,pphi,px,rv)
+ * convert coordinates and proper motion between coordinate systems
+ * Subroutine: wcscsys (cstring) returns code for coordinate system in string
+ * Subroutine: wcsceq (wcstring) returns equinox in years from system string
+ * Subroutine: wcscstr (sys,equinox,epoch) returns system string from equinox
+ * Subroutine: fk524 (ra,dec) Convert J2000(FK5) to B1950(FK4) coordinates
+ * Subroutine: fk524e (ra, dec, epoch) (more accurate for known position epoch)
+ * Subroutine: fk524m (ra,dec,rapm,decpm) exact
+ * Subroutine: fk524pv (ra,dec,rapm,decpm,parallax,rv) more exact
+ * Subroutine: fk425 (ra,dec) Convert B1950(FK4) to J2000(FK5) coordinates
+ * Subroutine: fk425e (ra, dec, epoch) (more accurate for known position epoch)
+ * Subroutine: fk425m (ra, dec, rapm, decpm) exact
+ * Subroutine: fk425pv (ra,dec,rapm,decpm,parallax,rv) more exact
+ * Subroutine: fk42gal (dtheta,dphi) Convert B1950(FK4) to galactic coordinates
+ * Subroutine: fk52gal (dtheta,dphi) Convert J2000(FK5) to galactic coordinates
+ * Subroutine: gal2fk4 (dtheta,dphi) Convert galactic coordinates to B1950(FK4)
+ * Subroutine: gal2fk5 (dtheta,dphi) Convert galactic coordinates to J2000<FK5)
+ * Subroutine: fk42ecl (dtheta,dphi,epoch) Convert B1950(FK4) to ecliptic coordinates
+ * Subroutine: fk52ecl (dtheta,dphi,epoch) Convert J2000(FK5) to ecliptic coordinates
+ * Subroutine: ecl2fk4 (dtheta,dphi,epoch) Convert ecliptic coordinates to B1950(FK4)
+ * Subroutine: ecl2fk5 (dtheta,dphi,epoch) Convert ecliptic coordinates to J2000<FK5)
+ * Subroutine: fk5prec (ep0, ep1, ra, dec) Precession ep0 to ep1, FK5 system
+ * Subroutine: fk4prec (ep0, ep1, ra, dec) Precession ep0 to ep1, FK4 system
+ * Subroutine: d2v3 (rra, rdec, r, pos) RA and Dec in degrees, Distance to Cartesian
+ * Subroutine: v2d3 (pos, rra, rdec, r) Cartesian to RA and Dec in degrees, Distance
+ * Subroutine: s2v3 (rra, rdec, r, pos) RA, Dec, Distance to Cartesian
+ * Subroutine: v2s3 (pos, rra, rdec, r) Cartesian to RA, Dec, Distance
+ * Subroutine: rotmat (axes, rot1, rot2, rot3, matrix) Rotation angles to matrix
+ *
+ * Note: Proper motions are always in RA/Dec degrees/year; no cos(Dec) correction
+ */
+
+#include <math.h>
+#ifndef VMS
+#include <stdlib.h>
+#endif
+#include <stdio.h> /* for fprintf() and sprintf() */
+#include <ctype.h>
+#include <string.h>
+#include "wcs.h"
+
+void fk524(), fk524e(), fk524m(), fk524pv();
+void fk425(), fk425e(), fk425m(), fk425pv();
+void fk42gal(), fk52gal(), gal2fk4(), gal2fk5();
+void fk42ecl(), fk52ecl(), ecl2fk4(), ecl2fk5();
+
+/* Convert from coordinate system sys1 to coordinate system sys2, converting
+ proper motions, too, and adding them if an epoch is specified */
+
+void
+wcsconp (sys1, sys2, eq1, eq2, ep1, ep2, dtheta, dphi, ptheta, pphi)
+
+int sys1; /* Input coordinate system (J2000, B1950, ECLIPTIC, GALACTIC */
+int sys2; /* Output coordinate system (J2000, B1950, ECLIPTIC, GALACTIC */
+double eq1; /* Input equinox (default of sys1 if 0.0) */
+double eq2; /* Output equinox (default of sys2 if 0.0) */
+double ep1; /* Input Besselian epoch in years (for proper motion) */
+double ep2; /* Output Besselian epoch in years (for proper motion) */
+double *dtheta; /* Longitude or right ascension in degrees
+ Input in sys1, returned in sys2 */
+double *dphi; /* Latitude or declination in degrees
+ Input in sys1, returned in sys2 */
+double *ptheta; /* Longitude or right ascension proper motion in RA degrees/year
+ Input in sys1, returned in sys2 */
+double *pphi; /* Latitude or declination proper motion in Dec degrees/year
+ Input in sys1, returned in sys2 */
+
+{
+ void fk5prec(), fk4prec();
+
+ /* Set equinoxes if 0.0 */
+ if (eq1 == 0.0) {
+ if (sys1 == WCS_B1950)
+ eq1 = 1950.0;
+ else
+ eq1 = 2000.0;
+ }
+ if (eq2 == 0.0) {
+ if (sys2 == WCS_B1950)
+ eq2 = 1950.0;
+ else
+ eq2 = 2000.0;
+ }
+
+ /* Set epochs if 0.0 */
+ if (ep1 == 0.0) {
+ if (sys1 == WCS_B1950)
+ ep1 = 1950.0;
+ else
+ ep1 = 2000.0;
+ }
+ if (ep2 == 0.0) {
+ if (sys2 == WCS_B1950)
+ ep2 = 1950.0;
+ else
+ ep2 = 2000.0;
+ }
+
+ if (sys1 == WCS_ICRS && sys2 == WCS_ICRS)
+ eq2 = eq1;
+
+ if (sys1 == WCS_J2000 && sys2 == WCS_ICRS && eq1 == 2000.0) {
+ eq2 = eq1;
+ sys1 = sys2;
+ }
+
+ /* Set systems and equinoxes so that ICRS coordinates are not precessed */
+ if (sys1 == WCS_ICRS && sys2 == WCS_J2000 && eq2 == 2000.0) {
+ eq1 = eq2;
+ sys1 = sys2;
+ }
+
+ /* If systems and equinoxes are the same, add proper motion and return */
+ if (sys2 == sys1 && eq1 == eq2) {
+ if (ep1 != ep2) {
+ if (sys1 == WCS_J2000) {
+ *dtheta = *dtheta + ((ep2 - ep1) * *ptheta);
+ *dphi = *dphi + ((ep2 - ep1) * *pphi);
+ }
+ else if (sys1 == WCS_B1950) {
+ *dtheta = *dtheta + ((ep2 - ep1) * *ptheta);
+ *dphi = *dphi + ((ep2 - ep1) * *pphi);
+ }
+ }
+ if (eq1 != eq2) {
+ if (sys1 == WCS_B1950)
+ fk4prec (eq1, eq2, dtheta, dphi);
+ if (sys1 == WCS_J2000)
+ fk5prec (eq1, 2000.0, dtheta, dphi);
+ }
+ return;
+ }
+
+ /* Precess from input equinox to input system equinox, if necessary */
+ if (sys1 == WCS_B1950 && eq1 != 1950.0)
+ fk4prec (eq1, 1950.0, dtheta, dphi);
+ if (sys1 == WCS_J2000 && eq1 != 2000.0)
+ fk5prec (eq1, 2000.0, dtheta, dphi);
+
+ /* Convert to B1950 FK4 */
+ if (sys2 == WCS_B1950) {
+ if (sys1 == WCS_J2000) {
+ if (*ptheta != 0.0 || *pphi != 0.0) {
+ fk524m (dtheta, dphi, ptheta, pphi);
+ if (ep2 != 1950.0) {
+ *dtheta = *dtheta + ((ep2 - 1950.0) * *ptheta);
+ *dphi = *dphi + ((ep2 - 1950.0) * *pphi);
+ }
+ }
+ else if (ep2 != 1950.0)
+ fk524e (dtheta, dphi, ep2);
+ else
+ fk524 (dtheta, dphi);
+ }
+ else if (sys1 == WCS_GALACTIC)
+ gal2fk4 (dtheta, dphi);
+ else if (sys1 == WCS_ECLIPTIC)
+ ecl2fk4 (dtheta, dphi, ep2);
+ }
+
+ else if (sys2 == WCS_J2000) {
+ if (sys1 == WCS_B1950) {
+ if (*ptheta != 0.0 || *pphi != 0.0) {
+ fk425m (dtheta, dphi, ptheta, pphi);
+ if (ep2 != 2000.0) {
+ *dtheta = *dtheta + ((ep2 - 2000.0) * *ptheta);
+ *dphi = *dphi + ((ep2 - 2000.0) * *pphi);
+ }
+ }
+ else if (ep2 > 0.0)
+ fk425e (dtheta, dphi, ep2);
+ else
+ fk425 (dtheta, dphi);
+ }
+ else if (sys1 == WCS_GALACTIC)
+ gal2fk5 (dtheta, dphi);
+ else if (sys1 == WCS_ECLIPTIC)
+ ecl2fk5 (dtheta, dphi, ep2);
+ }
+
+ else if (sys2 == WCS_GALACTIC) {
+ if (sys1 == WCS_B1950) {
+ if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
+ *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
+ *dphi = *dphi + (*pphi * (ep2 - ep1));
+ }
+ fk42gal (dtheta, dphi);
+ }
+ else if (sys1 == WCS_J2000) {
+ if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
+ *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
+ *dphi = *dphi + (*pphi * (ep2 - ep1));
+ }
+ fk52gal (dtheta, dphi);
+ }
+ else if (sys1 == WCS_ECLIPTIC) {
+ ecl2fk5 (dtheta, dphi, ep2);
+ fk52gal (dtheta, dphi);
+ }
+ }
+
+ else if (sys2 == WCS_ECLIPTIC) {
+ if (sys1 == WCS_B1950) {
+ if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
+ *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
+ *dphi = *dphi + (*pphi * (ep2 - ep1));
+ }
+ if (ep2 > 0.0)
+ fk42ecl (dtheta, dphi, ep2);
+ else
+ fk42ecl (dtheta, dphi, 1950.0);
+ }
+ else if (sys1 == WCS_J2000) {
+ if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
+ *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
+ *dphi = *dphi + (*pphi * (ep2 - ep1));
+ }
+ fk52ecl (dtheta, dphi, ep2);
+ }
+ else if (sys1 == WCS_GALACTIC) {
+ gal2fk5 (dtheta, dphi);
+ fk52ecl (dtheta, dphi, ep2);
+ }
+ }
+
+ /* Precess to desired equinox, if necessary */
+ if (sys2 == WCS_B1950 && eq2 != 1950.0)
+ fk4prec (1950.0, eq2, dtheta, dphi);
+ if (sys2 == WCS_J2000 && eq2 != 2000.0)
+ fk5prec (2000.0, eq2, dtheta, dphi);
+
+ /* Keep latitude/declination between +90 and -90 degrees */
+ if (*dphi > 90.0) {
+ *dphi = 180.0 - *dphi;
+ *dtheta = *dtheta + 180.0;
+ }
+ else if (*dphi < -90.0) {
+ *dphi = -180.0 - *dphi;
+ *dtheta = *dtheta + 180.0;
+ }
+
+ /* Keep longitude/right ascension between 0 and 360 degrees */
+ if (*dtheta > 360.0)
+ *dtheta = *dtheta - 360.0;
+ else if (*dtheta < 0.0)
+ *dtheta = *dtheta + 360.0;
+ return;
+}
+
+
+/* Convert from coordinate system sys1 to coordinate system sys2, converting
+ proper motions, too, and adding them if an epoch is specified */
+
+void
+wcsconv (sys1, sys2, eq1, eq2, ep1, ep2, dtheta, dphi, ptheta, pphi, px, rv)
+
+int sys1; /* Input coordinate system (J2000, B1950, ECLIPTIC, GALACTIC */
+int sys2; /* Output coordinate system (J2000, B1950, ECLIPTIC, GALACTIC */
+double eq1; /* Input equinox (default of sys1 if 0.0) */
+double eq2; /* Output equinox (default of sys2 if 0.0) */
+double ep1; /* Input Besselian epoch in years (for proper motion) */
+double ep2; /* Output Besselian epoch in years (for proper motion) */
+double *dtheta; /* Longitude or right ascension in degrees
+ Input in sys1, returned in sys2 */
+double *dphi; /* Latitude or declination in degrees
+ Input in sys1, returned in sys2 */
+double *ptheta; /* Longitude or right ascension proper motion in degrees/year
+ Input in sys1, returned in sys2 */
+double *pphi; /* Latitude or declination proper motion in degrees/year
+ Input in sys1, returned in sys2 */
+double *px; /* Parallax in arcseconds */
+double *rv; /* Radial velocity in km/sec */
+
+{
+ void fk5prec(), fk4prec();
+
+ /* Set equinoxes if 0.0 */
+ if (eq1 == 0.0) {
+ if (sys1 == WCS_B1950)
+ eq1 = 1950.0;
+ else
+ eq1 = 2000.0;
+ }
+ if (eq2 == 0.0) {
+ if (sys2 == WCS_B1950)
+ eq2 = 1950.0;
+ else
+ eq2 = 2000.0;
+ }
+
+ /* Set epochs if 0.0 */
+ if (ep1 == 0.0) {
+ if (sys1 == WCS_B1950)
+ ep1 = 1950.0;
+ else
+ ep1 = 2000.0;
+ }
+ if (ep2 == 0.0) {
+ if (sys2 == WCS_B1950)
+ ep2 = 1950.0;
+ else
+ ep2 = 2000.0;
+ }
+
+ /* Set systems and equinoxes so that ICRS coordinates are not precessed */
+ if (sys1 == WCS_ICRS && sys2 == WCS_ICRS)
+ eq2 = eq1;
+
+ if (sys1 == WCS_J2000 && sys2 == WCS_ICRS && eq1 == 2000.0) {
+ eq2 = eq1;
+ sys1 = sys2;
+ }
+
+ if (sys1 == WCS_ICRS && sys2 == WCS_J2000 && eq2 == 2000.0) {
+ eq1 = eq2;
+ sys1 = sys2;
+ }
+
+ /* If systems and equinoxes are the same, add proper motion and return */
+ if (sys2 == sys1 && eq1 == eq2) {
+ if (ep1 != ep2) {
+ if (sys1 == WCS_J2000) {
+ *dtheta = *dtheta + ((ep2 - ep1) * *ptheta);
+ *dphi = *dphi + ((ep2 - ep1) * *pphi);
+ }
+ else if (sys1 == WCS_B1950) {
+ *dtheta = *dtheta + ((ep2 - ep1) * *ptheta);
+ *dphi = *dphi + ((ep2 - ep1) * *pphi);
+ }
+ }
+ return;
+ }
+
+ /* Precess from input equinox to input system equinox, if necessary */
+ if (eq1 != eq2) {
+ if (sys1 == WCS_B1950 && eq1 != 1950.0)
+ fk4prec (eq1, 1950.0, dtheta, dphi);
+ if (sys1 == WCS_J2000 && eq1 != 2000.0)
+ fk5prec (eq1, 2000.0, dtheta, dphi);
+ }
+
+ /* Convert to B1950 FK4 */
+ if (sys2 == WCS_B1950) {
+ if (sys1 == WCS_J2000) {
+ if (*ptheta != 0.0 || *pphi != 0.0) {
+ if (*px != 0.0 || *rv != 0.0)
+ fk524pv (dtheta, dphi, ptheta, pphi, px, rv);
+ else
+ fk524m (dtheta, dphi, ptheta, pphi);
+ if (ep1 == 2000.0)
+ ep1 = 1950.0;
+ if (ep2 != 1950.0) {
+ *dtheta = *dtheta + ((ep2 - 1950.0) * *ptheta);
+ *dphi = *dphi + ((ep2 - 1950.0) * *pphi);
+ }
+ }
+ else if (ep2 != 1950.0)
+ fk524e (dtheta, dphi, ep2);
+ else
+ fk524 (dtheta, dphi);
+ }
+ else if (sys1 == WCS_GALACTIC)
+ gal2fk4 (dtheta, dphi);
+ else if (sys1 == WCS_ECLIPTIC)
+ ecl2fk4 (dtheta, dphi, ep2);
+ }
+
+ else if (sys2 == WCS_J2000) {
+ if (sys1 == WCS_B1950) {
+ if (*ptheta != 0.0 || *pphi != 0.0) {
+ if (*px != 0.0 || *rv != 0.0)
+ fk425pv (dtheta, dphi, ptheta, pphi, px, rv);
+ else
+ fk425m (dtheta, dphi, ptheta, pphi);
+ if (ep2 != 2000.0) {
+ *dtheta = *dtheta + ((ep2 - 2000.0) * *ptheta);
+ *dphi = *dphi + ((ep2 - 2000.0) * *pphi);
+ }
+ }
+ else if (ep2 > 0.0)
+ fk425e (dtheta, dphi, ep2);
+ else
+ fk425 (dtheta, dphi);
+ }
+ else if (sys1 == WCS_GALACTIC)
+ gal2fk5 (dtheta, dphi);
+ else if (sys1 == WCS_ECLIPTIC)
+ ecl2fk5 (dtheta, dphi, ep2);
+ }
+
+ else if (sys2 == WCS_GALACTIC) {
+ if (sys1 == WCS_B1950) {
+ if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
+ *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
+ *dphi = *dphi + (*pphi * (ep2 - ep1));
+ }
+ fk42gal (dtheta, dphi);
+ }
+ else if (sys1 == WCS_J2000) {
+ if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
+ *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
+ *dphi = *dphi + (*pphi * (ep2 - ep1));
+ }
+ fk52gal (dtheta, dphi);
+ }
+ else if (sys1 == WCS_ECLIPTIC) {
+ ecl2fk5 (dtheta, dphi, ep2);
+ fk52gal (dtheta, dphi);
+ }
+ }
+
+ else if (sys2 == WCS_ECLIPTIC) {
+ if (sys1 == WCS_B1950) {
+ if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
+ *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
+ *dphi = *dphi + (*pphi * (ep2 - ep1));
+ }
+ if (ep2 > 0.0)
+ fk42ecl (dtheta, dphi, ep2);
+ else
+ fk42ecl (dtheta, dphi, 1950.0);
+ }
+ else if (sys1 == WCS_J2000) {
+ if (ep2 != 0.0 && (*ptheta != 0.0 || *pphi != 0.0)) {
+ *dtheta = *dtheta + (*ptheta * (ep2 - ep1));
+ *dphi = *dphi + (*pphi * (ep2 - ep1));
+ }
+ fk52ecl (dtheta, dphi, ep2);
+ }
+ else if (sys1 == WCS_GALACTIC) {
+ gal2fk5 (dtheta, dphi);
+ fk52ecl (dtheta, dphi, ep2);
+ }
+ }
+
+ /* Precess to desired equinox, if necessary */
+ if (eq1 != eq2) {
+ if (sys2 == WCS_B1950 && eq2 != 1950.0)
+ fk4prec (1950.0, eq2, dtheta, dphi);
+ if (sys2 == WCS_J2000 && eq2 != 2000.0)
+ fk5prec (2000.0, eq2, dtheta, dphi);
+ }
+
+ /* Keep latitude/declination between +90 and -90 degrees */
+ if (*dphi > 90.0) {
+ *dphi = 180.0 - *dphi;
+ *dtheta = *dtheta + 180.0;
+ }
+ else if (*dphi < -90.0) {
+ *dphi = -180.0 - *dphi;
+ *dtheta = *dtheta + 180.0;
+ }
+
+ /* Keep longitude/right ascension between 0 and 360 degrees */
+ if (*dtheta > 360.0)
+ *dtheta = *dtheta - 360.0;
+ else if (*dtheta < 0.0)
+ *dtheta = *dtheta + 360.0;
+ return;
+}
+
+
+/* Convert from coordinate system sys1 to coordinate system sys2 */
+
+void
+wcscon (sys1, sys2, eq1, eq2, dtheta, dphi, epoch)
+
+int sys1; /* Input coordinate system (J2000, B1950, ECLIPTIC, GALACTIC */
+int sys2; /* Output coordinate system (J2000, B1950, ECLIPTIC, GALACTIC */
+double eq1; /* Input equinox (default of sys1 if 0.0) */
+double eq2; /* Output equinox (default of sys2 if 0.0) */
+double *dtheta; /* Longitude or right ascension in degrees
+ Input in sys1, returned in sys2 */
+double *dphi; /* Latitude or declination in degrees
+ Input in sys1, returned in sys2 */
+double epoch; /* Besselian epoch in years */
+
+{
+ void fk5prec(), fk4prec();
+
+ /* Set equinoxes if 0.0 */
+ if (eq1 == 0.0) {
+ if (sys1 == WCS_B1950)
+ eq1 = 1950.0;
+ else
+ eq1 = 2000.0;
+ }
+ if (eq2 == 0.0) {
+ if (sys2 == WCS_B1950)
+ eq2 = 1950.0;
+ else
+ eq2 = 2000.0;
+ }
+
+ /* Set systems and equinoxes so that ICRS coordinates are not precessed */
+ if (sys1 == WCS_ICRS && sys2 == WCS_ICRS)
+ eq2 = eq1;
+
+ if (sys1 == WCS_J2000 && sys2 == WCS_ICRS && eq1 == 2000.0) {
+ eq2 = eq1;
+ sys1 = sys2;
+ }
+
+ if (sys1 == WCS_ICRS && sys2 == WCS_J2000 && eq2 == 2000.0) {
+ eq1 = eq2;
+ sys1 = sys2;
+ }
+
+ /* If systems and equinoxes are the same, return */
+ if (sys2 == sys1 && eq1 == eq2)
+ return;
+
+ /* Precess from input equinox, if necessary */
+ if (eq1 != eq2) {
+ if (sys1 == WCS_B1950 && eq1 != 1950.0)
+ fk4prec (eq1, 1950.0, dtheta, dphi);
+ if (sys1 == WCS_J2000 && eq1 != 2000.0)
+ fk5prec (eq1, 2000.0, dtheta, dphi);
+ }
+
+ /* Convert to B1950 FK4 */
+ if (sys2 == WCS_B1950) {
+ if (sys1 == WCS_J2000) {
+ if (epoch > 0)
+ fk524e (dtheta, dphi, epoch);
+ else
+ fk524 (dtheta, dphi);
+ }
+ else if (sys1 == WCS_GALACTIC)
+ gal2fk4 (dtheta, dphi);
+ else if (sys1 == WCS_ECLIPTIC) {
+ if (epoch > 0)
+ ecl2fk4 (dtheta, dphi, epoch);
+ else
+ ecl2fk4 (dtheta, dphi, 1950.0);
+ }
+ }
+
+ else if (sys2 == WCS_J2000) {
+ if (sys1 == WCS_B1950) {
+ if (epoch > 0)
+ fk425e (dtheta, dphi, epoch);
+ else
+ fk425 (dtheta, dphi);
+ }
+ else if (sys1 == WCS_GALACTIC)
+ gal2fk5 (dtheta, dphi);
+ else if (sys1 == WCS_ECLIPTIC) {
+ if (epoch > 0)
+ ecl2fk5 (dtheta, dphi, epoch);
+ else
+ ecl2fk5 (dtheta, dphi, 2000.0);
+ }
+ }
+
+ else if (sys2 == WCS_GALACTIC) {
+ if (sys1 == WCS_B1950)
+ fk42gal (dtheta, dphi);
+ else if (sys1 == WCS_J2000)
+ fk52gal (dtheta, dphi);
+ else if (sys1 == WCS_ECLIPTIC) {
+ if (epoch > 0)
+ ecl2fk5 (dtheta, dphi, epoch);
+ else
+ ecl2fk5 (dtheta, dphi, 2000.0);
+ fk52gal (dtheta, dphi);
+ }
+ }
+
+ else if (sys2 == WCS_ECLIPTIC) {
+ if (sys1 == WCS_B1950) {
+ if (epoch > 0)
+ fk42ecl (dtheta, dphi, epoch);
+ else
+ fk42ecl (dtheta, dphi, 1950.0);
+ }
+ else if (sys1 == WCS_J2000) {
+ if (epoch > 0)
+ fk52ecl (dtheta, dphi, epoch);
+ else
+ fk52ecl (dtheta, dphi, 2000.0);
+ }
+ else if (sys1 == WCS_GALACTIC) {
+ gal2fk5 (dtheta, dphi);
+ if (epoch > 0)
+ fk52ecl (dtheta, dphi, epoch);
+ else
+ fk52ecl (dtheta, dphi, 2000.0);
+ }
+ }
+
+ /* Precess to desired equinox, if necessary */
+ if (eq1 != eq2) {
+ if (sys2 == WCS_B1950 && eq2 != 1950.0)
+ fk4prec (1950.0, eq2, dtheta, dphi);
+ if (sys2 == WCS_J2000 && eq2 != 2000.0)
+ fk5prec (2000.0, eq2, dtheta, dphi);
+ }
+
+ /* Keep latitude/declination between +90 and -90 degrees */
+ if (*dphi > 90.0) {
+ *dphi = 180.0 - *dphi;
+ *dtheta = *dtheta + 180.0;
+ }
+ else if (*dphi < -90.0) {
+ *dphi = -180.0 - *dphi;
+ *dtheta = *dtheta + 180.0;
+ }
+
+ /* Keep longitude/right ascension between 0 and 360 degrees */
+ if (*dtheta > 360.0)
+ *dtheta = *dtheta - 360.0;
+ else if (*dtheta < 0.0)
+ *dtheta = *dtheta + 360.0;
+
+ return;
+}
+
+
+/* Set coordinate system from string */
+int
+wcscsys (wcstring)
+
+char *wcstring; /* Name of coordinate system */
+{
+ double equinox;
+
+ if (wcstring[0] == 'J' || wcstring[0] == 'j' ||
+ !strcmp (wcstring,"2000") || !strcmp (wcstring, "2000.0") ||
+ !strcmp (wcstring,"ICRS") || !strcmp (wcstring, "icrs") ||
+ !strncmp (wcstring,"FK5",3) || !strncmp (wcstring, "fk5",3))
+ return WCS_J2000;
+
+ if (wcstring[0] == 'B' || wcstring[0] == 'b' ||
+ !strcmp (wcstring,"1950") || !strcmp (wcstring, "1950.0") ||
+ !strncmp (wcstring,"FK4",3) || !strncmp (wcstring, "fk4",3))
+ return WCS_B1950;
+
+ else if (wcstring[0] == 'I' || wcstring[0] == 'i' )
+ return WCS_ICRS;
+
+ else if (wcstring[0] == 'G' || wcstring[0] == 'g' )
+ return WCS_GALACTIC;
+
+ else if (wcstring[0] == 'E' || wcstring[0] == 'e' )
+ return WCS_ECLIPTIC;
+
+ else if (wcstring[0] == 'A' || wcstring[0] == 'a' )
+ return WCS_ALTAZ;
+
+ else if (wcstring[0] == 'N' || wcstring[0] == 'n' )
+ return WCS_NPOLE;
+
+ else if (wcstring[0] == 'L' || wcstring[0] == 'l' )
+ return WCS_LINEAR;
+
+ else if (!strncasecmp (wcstring, "pixel", 5))
+ return WCS_XY;
+
+ else if (wcstring[0] == 'P' || wcstring[0] == 'p' )
+ return WCS_PLANET;
+
+ else if (isnum (wcstring)) {
+ equinox = atof (wcstring);
+ if (equinox > 1980.0)
+ return WCS_J2000;
+ else if (equinox > 1900.0)
+ return WCS_B1950;
+ else
+ return -1;
+ }
+ else
+ return -1;
+}
+
+
+/* Set equinox from string (return 0.0 if not obvious) */
+
+double
+wcsceq (wcstring)
+
+char *wcstring; /* Name of coordinate system */
+{
+ if (wcstring[0] == 'J' || wcstring[0] == 'j' ||
+ wcstring[0] == 'B' || wcstring[0] == 'b')
+ return (atof (wcstring+1));
+ else if (!strncmp (wcstring, "FK4",3) ||
+ !strncmp (wcstring, "fk4",3))
+ return (1950.0);
+ else if (!strncmp (wcstring, "FK5",3) ||
+ !strncmp (wcstring, "fk5",3))
+ return (2000.0);
+ else if (!strncmp (wcstring, "ICRS",4) ||
+ !strncmp (wcstring, "icrs",4))
+ return (2000.0);
+ else if (wcstring[0] == '1' || wcstring[0] == '2')
+ return (atof (wcstring));
+ else
+ return (0.0);
+}
+
+
+/* Set coordinate system type string from system and equinox */
+
+void
+wcscstr (cstr, syswcs, equinox, epoch)
+
+char *cstr; /* Coordinate system string (returned) */
+int syswcs; /* Coordinate system code */
+double equinox; /* Equinox of coordinate system */
+double epoch; /* Epoch of coordinate system */
+{
+
+ char *estr;
+
+ if (syswcs == WCS_XY) {
+ strcpy (cstr, "XY");
+ return;
+ }
+
+ /* Try to figure out coordinate system if it is not set */
+ if (epoch == 0.0)
+ epoch = equinox;
+ if (syswcs < 0) {
+ if (equinox > 0.0) {
+ if (equinox == 2000.0)
+ syswcs = WCS_J2000;
+ else if (equinox == 1950.0)
+ syswcs = WCS_B1950;
+ }
+ else if (epoch > 0.0) {
+ if (epoch > 1980.0) {
+ syswcs = WCS_J2000;
+ equinox = 2000.0;
+ }
+ else {
+ syswcs = WCS_B1950;
+ equinox = 1950.0;
+ }
+ }
+ else
+ syswcs = WCS_J2000;
+ }
+
+ /* Set coordinate system string from system flag and epoch */
+ if (syswcs == WCS_B1950) {
+ if (epoch == 1950.0 || epoch == 0.0)
+ strcpy (cstr, "B1950");
+ else
+ sprintf (cstr, "B%7.2f", equinox);
+ if ((estr = strsrch (cstr,".00")) != NULL) {
+ estr[0] = (char) 0;
+ estr[1] = (char) 0;
+ estr[2] = (char) 0;
+ }
+ }
+ else if (syswcs == WCS_GALACTIC)
+ strcpy (cstr, "galactic");
+ else if (syswcs == WCS_ECLIPTIC)
+ strcpy (cstr, "ecliptic");
+ else if (syswcs == WCS_J2000) {
+ if (epoch == 2000.0 || epoch == 0.0)
+ strcpy (cstr, "J2000");
+ else
+ sprintf (cstr, "J%7.2f", equinox);
+ if ((estr = strsrch (cstr,".00")) != NULL) {
+ estr[0] = (char) 0;
+ estr[1] = (char) 0;
+ estr[2] = (char) 0;
+ }
+ }
+ else if (syswcs == WCS_ICRS) {
+ strcpy (cstr, "ICRS");
+ }
+ else if (syswcs == WCS_PLANET) {
+ strcpy (cstr, "PLANET");
+ }
+ else if (syswcs == WCS_LINEAR || syswcs == WCS_XY) {
+ strcpy (cstr, "LINEAR");
+ }
+ return;
+}
+
+
+/* Constant vector and matrix (by columns)
+ These values were obtained by inverting C.Hohenkerk's forward matrix
+ (private communication), which agrees with the one given in reference
+ 2 but which has one additional decimal place. */
+
+static double a[3] = {-1.62557e-6, -0.31919e-6, -0.13843e-6};
+static double ad[3] = {1.245e-3, -1.580e-3, -0.659e-3};
+static double d2pi = 6.283185307179586476925287; /* two PI */
+static double tiny = 1.e-30; /* small number to avoid arithmetic problems */
+
+/* FK524 convert J2000 FK5 star data to B1950 FK4
+ based on Starlink sla_fk524 by P.T.Wallace 27 October 1987 */
+
+static double emi[6][6] = {
+ { 0.9999256795, /* emi[0][0] */
+ 0.0111814828, /* emi[0][1] */
+ 0.0048590039, /* emi[0][2] */
+ -0.00000242389840, /* emi[0][3] */
+ -0.00000002710544, /* emi[0][4] */
+ -0.00000001177742 }, /* emi[0][5] */
+
+ { -0.0111814828, /* emi[1][0] */
+ 0.9999374849, /* emi[1][1] */
+ -0.0000271771, /* emi[1][2] */
+ 0.00000002710544, /* emi[1][3] */
+ -0.00000242392702, /* emi[1][4] */
+ 0.00000000006585 }, /* emi[1][5] */
+
+ { -0.0048590040, /* emi[2][0] */
+ -0.0000271557, /* emi[2][1] */
+ 0.9999881946, /* emi[2][2] */
+ 0.00000001177742, /* emi[2][3] */
+ 0.00000000006585, /* emi[2][4] */
+ -0.00000242404995 }, /* emi[2][5] */
+
+ { -0.000551, /* emi[3][0] */
+ 0.238509, /* emi[3][1] */
+ -0.435614, /* emi[3][2] */
+ 0.99990432, /* emi[3][3] */
+ 0.01118145, /* emi[3][4] */
+ 0.00485852 }, /* emi[3][5] */
+
+ { -0.238560, /* emi[4][0] */
+ -0.002667, /* emi[4][1] */
+ 0.012254, /* emi[4][2] */
+ -0.01118145, /* emi[4][3] */
+ 0.99991613, /* emi[4][4] */
+ -0.00002717 }, /* emi[4][5] */
+
+ { 0.435730, /* emi[5][0] */
+ -0.008541, /* emi[5][1] */
+ 0.002117, /* emi[5][2] */
+ -0.00485852, /* emi[5][3] */
+ -0.00002716, /* emi[5][4] */
+ 0.99996684 } /* emi[5][5] */
+ };
+
+void
+fk524 (ra,dec)
+
+double *ra; /* Right ascension in degrees (J2000 in, B1950 out) */
+double *dec; /* Declination in degrees (J2000 in, B1950 out) */
+
+{
+ double rapm; /* Proper motion in right ascension */
+ double decpm; /* Proper motion in declination */
+ /* In: deg/jul.yr. Out: deg/trop.yr. */
+
+ rapm = (double) 0.0;
+ decpm = (double) 0.0;
+ fk524m (ra, dec, &rapm, &decpm);
+ return;
+}
+
+void
+fk524e (ra, dec, epoch)
+
+double *ra; /* Right ascension in degrees (J2000 in, B1950 out) */
+double *dec; /* Declination in degrees (J2000 in, B1950 out) */
+double epoch; /* Besselian epoch in years */
+
+{
+ double rapm; /* Proper motion in right ascension */
+ double decpm; /* Proper motion in declination */
+ /* In: deg/jul.yr. Out: deg/trop.yr. */
+
+ rapm = (double) 0.0;
+ decpm = (double) 0.0;
+ fk524m (ra, dec, &rapm, &decpm);
+ *ra = *ra + (rapm * (epoch - 1950.0));
+ *dec = *dec + (decpm * (epoch - 1950.0));
+ return;
+}
+
+void
+fk524m (ra,dec,rapm,decpm)
+
+double *ra; /* Right ascension in degrees (J2000 in, B1950 out) */
+double *dec; /* Declination in degrees (J2000 in, B1950 out) */
+double *rapm; /* Proper motion in right ascension */
+double *decpm; /* Proper motion in declination */
+ /* In: ra/dec deg/jul.yr. Out: ra/dec deg/trop.yr. */
+
+{
+ double parallax = 0.0;
+ double rv = 0.0;
+
+ fk524pv (ra, dec, rapm, decpm, &parallax, &rv);
+ return;
+}
+
+
+void
+fk524pv (ra,dec,rapm,decpm, parallax, rv)
+
+double *ra; /* Right ascension in degrees (J2000 in, B1950 out) */
+double *dec; /* Declination in degrees (J2000 in, B1950 out) */
+double *rapm; /* Proper motion in right ascension */
+double *decpm; /* Proper motion in declination
+ * In: ra/dec degrees/Julian year (not ra*cos(dec))
+ * Out: ra/dec degrees/tropical year */
+double *parallax; /* Parallax (arcsec) */
+double *rv; /* Rradial velocity (km/s, +ve = moving away) */
+
+/* This routine converts stars from the IAU 1976 FK5 Fricke
+ system, to the old Bessel-Newcomb FK4 system, using Yallop's
+ implementation (see ref 2) of a matrix method due to Standish
+ (see ref 3). The numerical values of ref 2 are used canonically.
+
+ * Conversion from other than Julian epoch 2000.0 to other than Besselian
+ epoch 1950.0 will require use of the appropriate precession, proper
+ motion, and e-terms routines before and/or after fk524 is called.
+
+ * In the FK4 catalogue the proper motions of stars within 10 degrees
+ of the poles do not embody the differential e-term effect and should,
+ strictly speaking, be handled in a different manner from stars outside
+ these regions. however, given the general lack of homogeneity of the
+ star data available for routine astrometry, the difficulties of handling
+ positions that may have been determined from astrometric fields spanning
+ the polar and non-polar regions, the likelihood that the differential
+ e-terms effect was not taken into account when allowing for proper motion
+ in past astrometry, and the undesirability of a discontinuity in the
+ algorithm, the decision has been made in this routine to include the
+ effect of differential e-terms on the proper motions for all stars,
+ whether polar or not, at epoch 2000, and measuring on the sky rather
+ than in terms of dra, the errors resulting from this simplification are
+ less than 1 milliarcsecond in position and 1 milliarcsecond per century
+ in proper motion.
+
+ References:
+
+ 1 "Mean and apparent place computations in the new IAU System.
+ I. The transformation of astrometric catalog systems to the
+ equinox J2000.0." Smith, C.A.; Kaplan, G.H.; Hughes, J.A.;
+ Seidelmann, P.K.; Yallop, B.D.; Hohenkerk, C.Y.
+ Astronomical Journal vol. 97, Jan. 1989, p. 265-273.
+
+ 2 "Mean and apparent place computations in the new IAU System.
+ II. Transformation of mean star places from FK4 B1950.0 to
+ FK5 J2000.0 using matrices in 6-space." Yallop, B.D.;
+ Hohenkerk, C.Y.; Smith, C.A.; Kaplan, G.H.; Hughes, J.A.;
+ Seidelmann, P.K.; Astronomical Journal vol. 97, Jan. 1989,
+ p. 274-279.
+
+ 3 Seidelmann, P.K. (ed), 1992. "Explanatory Supplement to
+ the Astronomical Almanac", ISBN 0-935702-68-7.
+
+ 4 "Conversion of positions and proper motions from B1950.0 to the
+ IAU system at J2000.0", Standish, E.M. Astronomy and
+ Astrophysics, vol. 115, no. 1, Nov. 1982, p. 20-22.
+
+ P.T.Wallace Starlink 19 December 1993
+ Doug Mink Smithsonian Astrophysical Observatory 1 November 2000 */
+
+{
+ double r2000,d2000; /* J2000.0 ra,dec (radians) */
+ double r1950,d1950; /* B1950.0 ra,dec (rad) */
+
+ /* Miscellaneous */
+ double ur,ud;
+ double sr, cr, sd, cd, x, y, z, w, wd;
+ double v1[6],v2[6];
+ double xd,yd,zd;
+ double rxyz, rxysq, rxy;
+ double dra,ddec;
+ int i,j;
+ int diag = 0;
+
+ /* Constants */
+ double zero = (double) 0.0;
+ double vf = 21.095; /* Km per sec to AU per tropical century */
+ /* = 86400 * 36524.2198782 / 149597870 */
+
+ /* Convert J2000 RA and Dec from degrees to radians */
+ r2000 = degrad (*ra);
+ d2000 = degrad (*dec);
+
+ /* Convert J2000 RA and Dec proper motion from degrees/year to arcsec/tc */
+ ur = *rapm * 360000.0;
+ ud = *decpm * 360000.0;
+
+ /* Spherical to Cartesian */
+ sr = sin (r2000);
+ cr = cos (r2000);
+ sd = sin (d2000);
+ cd = cos (d2000);
+
+ x = cr * cd;
+ y = sr * cd;
+ z = sd;
+
+ v1[0] = x;
+ v1[1] = y;
+ v1[2] = z;
+
+ if (ur != zero || ud != zero) {
+ v1[3] = -(ur*y) - (cr*sd*ud);
+ v1[4] = (ur*x) - (sr*sd*ud);
+ v1[5] = (cd*ud);
+ }
+ else {
+ v1[3] = zero;
+ v1[4] = zero;
+ v1[5] = zero;
+ }
+
+ /* Convert position + velocity vector to bn system */
+ for (i = 0; i < 6; i++) {
+ w = zero;
+ for (j = 0; j < 6; j++) {
+ w = w + emi[i][j] * v1[j];
+ }
+ v2[i] = w;
+ }
+
+ /* Vector components */
+ x = v2[0];
+ y = v2[1];
+ z = v2[2];
+
+ /* Magnitude of position vector */
+ rxyz = sqrt (x*x + y*y + z*z);
+
+ /* Apply e-terms to position */
+ w = (x * a[0]) + (y * a[1]) + (z * a[2]);
+ x = x + (a[0] * rxyz) - (w * x);
+ y = y + (a[1] * rxyz) - (w * y);
+ z = z + (a[2] * rxyz) - (w * z);
+
+ /* Recompute magnitude of position vector */
+ rxyz = sqrt (x*x + y*y + z*z);
+
+ /* Apply e-terms to position and velocity */
+ x = v2[0];
+ y = v2[1];
+ z = v2[2];
+ w = (x * a[0]) + (y * a[1]) + (z * a[2]);
+ wd = (x * ad[0]) + (y * ad[1]) + (z * ad[2]);
+ x = x + (a[0] * rxyz) - (w * x);
+ y = y + (a[1] * rxyz) - (w * y);
+ z = z + (a[2] * rxyz) - (w * z);
+ xd = v2[3] + (ad[0] * rxyz) - (wd * x);
+ yd = v2[4] + (ad[1] * rxyz) - (wd * y);
+ zd = v2[5] + (ad[2] * rxyz) - (wd * z);
+
+ /* Convert to spherical */
+ rxysq = (x * x) + (y * y);
+ rxy = sqrt (rxysq);
+
+ /* Convert back to spherical coordinates */
+ if (x == zero && y == zero)
+ r1950 = zero;
+ else {
+ r1950 = atan2 (y,x);
+ if (r1950 < zero)
+ r1950 = r1950 + d2pi;
+ }
+ d1950 = atan2 (z,rxy);
+
+ if (rxy > tiny) {
+ ur = (x*yd - y*xd) / rxysq;
+ ud = (zd*rxysq - z * (x*xd + y*yd)) / ((rxysq + z*z) * rxy);
+ }
+
+ if (*parallax > tiny) {
+ *rv = ((x * xd) + (y * yd) + (z * zd)) / (*parallax * vf * rxyz);
+ *parallax = *parallax / rxyz;
+ }
+
+ /* Return results */
+ *ra = raddeg (r1950);
+ *dec = raddeg (d1950);
+ *rapm = ur / 360000.0;
+ *decpm = ud / 360000.0;
+
+ if (diag) {
+ dra = 240.0 * raddeg (r1950 - r2000);
+ ddec = 3600.0 * raddeg (d1950 - d2000);
+ fprintf(stderr,"B1950-J2000: dra= %11.5f sec ddec= %f11.5f arcsec\n",
+ dra, ddec);
+ }
+
+ return;
+}
+
+
+/* Convert B1950.0 FK4 star data to J2000.0 FK5 */
+static double em[6][6] = {
+ { 0.9999256782, /* em[0][0] */
+ -0.0111820611, /* em[0][1] */
+ -0.0048579477, /* em[0][2] */
+ 0.00000242395018, /* em[0][3] */
+ -0.00000002710663, /* em[0][4] */
+ -0.00000001177656 }, /* em[0][5] */
+
+ { 0.0111820610, /* em[1][0] */
+ 0.9999374784, /* em[1][1] */
+ -0.0000271765, /* em[1][2] */
+ 0.00000002710663, /* em[1][3] */
+ 0.00000242397878, /* em[1][4] */
+ -0.00000000006587 }, /* em[1][5] */
+
+ { 0.0048579479, /* em[2][0] */
+ -0.0000271474, /* em[2][1] */
+ 0.9999881997, /* em[2][2] */
+ 0.00000001177656, /* em[2][3] */
+ -0.00000000006582, /* em[2][4] */
+ 0.00000242410173 }, /* em[2][5] */
+
+ { -0.000551, /* em[3][0] */
+ -0.238565, /* em[3][1] */
+ 0.435739, /* em[3][2] */
+ 0.99994704, /* em[3][3] */
+ -0.01118251, /* em[3][4] */
+ -0.00485767 }, /* em[3][5] */
+
+ { 0.238514, /* em[4][0] */
+ -0.002667, /* em[4][1] */
+ -0.008541, /* em[4][2] */
+ 0.01118251, /* em[4][3] */
+ 0.99995883, /* em[4][4] */
+ -0.00002718 }, /* em[4][5] */
+
+ { -0.435623, /* em[5][0] */
+ 0.012254, /* em[5][1] */
+ 0.002117, /* em[5][2] */
+ 0.00485767, /* em[5][3] */
+ -0.00002714, /* em[5][4] */
+ 1.00000956 } /* em[5][5] */
+ };
+
+void
+fk425 (ra, dec)
+
+double *ra; /* Right ascension in degrees (B1950 in, J2000 out) */
+double *dec; /* Declination in degrees (B1950 in, J2000 out) */
+
+{
+double rapm; /* Proper motion in right ascension */
+double decpm; /* Proper motion in declination */
+ /* In: rad/trop.yr. Out: rad/jul.yr. */
+
+ rapm = (double) 0.0;
+ decpm = (double) 0.0;
+ fk425m (ra, dec, &rapm, &decpm);
+ return;
+}
+
+
+void
+fk425e (ra, dec, epoch)
+
+double *ra; /* Right ascension in degrees (B1950 in, J2000 out) */
+double *dec; /* Declination in degrees (B1950 in, J2000 out) */
+double epoch; /* Besselian epoch in years */
+{
+double rapm; /* Proper motion in right ascension */
+double decpm; /* Proper motion in declination */
+ /* In: rad/trop.yr. Out: rad/jul.yr. */
+
+ rapm = (double) 0.0;
+ decpm = (double) 0.0;
+ fk425m (ra, dec, &rapm, &decpm);
+ *ra = *ra + (rapm * (epoch - 2000.0));
+ *dec = *dec + (decpm * (epoch - 2000.0));
+ return;
+}
+
+void
+fk425m (ra, dec, rapm, decpm)
+
+double *ra, *dec; /* Right ascension and declination in degrees
+ input: B1950.0,FK4 returned: J2000.0,FK5 */
+double *rapm, *decpm; /* Proper motion in right ascension and declination
+ input: B1950.0,FK4 returned: J2000.0,FK5
+ ra/dec deg/trop.yr. ra/dec deg/jul.yr. */
+{
+ double parallax = 0.0;
+ double rv = 0.0;
+
+ fk425pv (ra, dec, rapm, decpm, &parallax, &rv);
+ return;
+}
+
+
+void
+fk425pv (ra,dec,rapm,decpm, parallax, rv)
+
+double *ra; /* Right ascension in degrees (J2000 in, B1950 out) */
+double *dec; /* Declination in degrees (J2000 in, B1950 out) */
+double *rapm; /* Proper motion in right ascension */
+double *decpm; /* Proper motion in declination
+ * In: ra/dec degrees/Julian year (not ra*cos(dec))
+ * Out: ra/dec degrees/tropical year */
+double *parallax; /* Parallax (arcsec) */
+double *rv; /* Rradial velocity (km/s, +ve = moving away) */
+
+/* This routine converts stars from the old Bessel-Newcomb FK4 system
+ to the IAU 1976 FK5 Fricke system, using Yallop's implementation
+ (see ref 2) of a matrix method due to Standish (see ref 3). The
+ numerical values of ref 2 are used canonically.
+
+ * Conversion from other than Besselian epoch 1950.0 to other than Julian
+ epoch 2000.0 will require use of the appropriate precession, proper
+ motion, and e-terms routines before and/or after fk425 is called.
+
+ * In the FK4 catalogue the proper motions of stars within 10 degrees
+ of the poles do not embody the differential e-term effect and should,
+ strictly speaking, be handled in a different manner from stars outside
+ these regions. however, given the general lack of homogeneity of the
+ star data available for routine astrometry, the difficulties of handling
+ positions that may have been determined from astrometric fields spanning
+ the polar and non-polar regions, the likelihood that the differential
+ e-terms effect was not taken into account when allowing for proper motion
+ in past astrometry, and the undesirability of a discontinuity in the
+ algorithm, the decision has been made in this routine to include the
+ effect of differential e-terms on the proper motions for all stars,
+ whether polar or not, at epoch 2000, and measuring on the sky rather
+ than in terms of dra, the errors resulting from this simplification are
+ less than 1 milliarcsecond in position and 1 milliarcsecond per century
+ in proper motion.
+
+ References:
+
+ 1 "Mean and apparent place computations in the new IAU System.
+ I. The transformation of astrometric catalog systems to the
+ equinox J2000.0." Smith, C.A.; Kaplan, G.H.; Hughes, J.A.;
+ Seidelmann, P.K.; Yallop, B.D.; Hohenkerk, C.Y.
+ Astronomical Journal vol. 97, Jan. 1989, p. 265-273.
+
+ 2 "Mean and apparent place computations in the new IAU System.
+ II. Transformation of mean star places from FK4 B1950.0 to
+ FK5 J2000.0 using matrices in 6-space." Yallop, B.D.;
+ Hohenkerk, C.Y.; Smith, C.A.; Kaplan, G.H.; Hughes, J.A.;
+ Seidelmann, P.K.; Astronomical Journal vol. 97, Jan. 1989,
+ p. 274-279.
+
+ 3 "Conversion of positions and proper motions from B1950.0 to the
+ IAU system at J2000.0", Standish, E.M. Astronomy and
+ Astrophysics, vol. 115, no. 1, Nov. 1982, p. 20-22.
+
+ P.T.Wallace Starlink 20 December 1993
+ Doug Mink Smithsonian Astrophysical Observatory 7 June 1995 */
+
+{
+ double r1950,d1950; /* B1950.0 ra,dec (rad) */
+ double r2000,d2000; /* J2000.0 ra,dec (rad) */
+
+ /* Miscellaneous */
+ double ur,ud,sr,cr,sd,cd,w,wd;
+ double x,y,z,xd,yd,zd, dra,ddec;
+ double rxyz, rxysq, rxy, rxyzsq, spxy, spxyz;
+ int i,j;
+ int diag = 0;
+
+ double r0[3],rd0[3]; /* star position and velocity vectors */
+ double v1[6],v2[6]; /* combined position and velocity vectors */
+
+ /* Constants */
+ double zero = (double) 0.0;
+ double vf = 21.095; /* Km per sec to AU per tropical century */
+ /* = 86400 * 36524.2198782 / 149597870 */
+
+ /* Convert B1950 RA and Dec from degrees to radians */
+ r1950 = degrad (*ra);
+ d1950 = degrad (*dec);
+
+ /* Convert B1950 RA and Dec proper motion from degrees/year to arcsec/tc */
+ ur = *rapm * 360000.0;
+ ud = *decpm * 360000.0;
+
+ /* Convert direction to Cartesian */
+ sr = sin (r1950);
+ cr = cos (r1950);
+ sd = sin (d1950);
+ cd = cos (d1950);
+ r0[0] = cr * cd;
+ r0[1] = sr * cd;
+ r0[2] = sd;
+
+ /* Convert motion to Cartesian */
+ w = vf * *rv * *parallax;
+ if (ur != zero || ud != zero || (*rv != zero && *parallax != zero)) {
+ rd0[0] = (-sr * cd * ur) - (cr * sd * ud) + (w * r0[0]);
+ rd0[1] = (cr * cd * ur) - (sr * sd * ud) + (w * r0[1]);
+ rd0[2] = (cd * ud) + (w * r0[2]);
+ }
+ else {
+ rd0[0] = zero;
+ rd0[1] = zero;
+ rd0[2] = zero;
+ }
+
+ /* Remove e-terms from position and express as position+velocity 6-vector */
+ w = (r0[0] * a[0]) + (r0[1] * a[1]) + (r0[2] * a[2]);
+ for (i = 0; i < 3; i++)
+ v1[i] = r0[i] - a[i] + (w * r0[i]);
+
+ /* Remove e-terms from proper motion and express as 6-vector */
+ wd = (r0[0] * ad[0]) + (r0[1] * ad[1]) + (r0[2] * ad[2]);
+ for (i = 0; i < 3; i++)
+ v1[i+3] = rd0[i] - ad[i] + (wd * r0[i]);
+
+ /* Alternately: Put proper motion in 6-vector without adding e-terms
+ for (i = 0; i < 3; i++)
+ v1[i+3] = rd0[i]; */
+
+ /* Convert position + velocity vector to FK5 system */
+ for (i = 0; i < 6; i++) {
+ w = zero;
+ for (j = 0; j < 6; j++) {
+ w += em[i][j] * v1[j];
+ }
+ v2[i] = w;
+ }
+
+ /* Vector components */
+ x = v2[0];
+ y = v2[1];
+ z = v2[2];
+ xd = v2[3];
+ yd = v2[4];
+ zd = v2[5];
+
+ /* Magnitude of position vector */
+ rxysq = x*x + y*y;
+ rxy = sqrt (rxysq);
+ rxyzsq = rxysq + z*z;
+ rxyz = sqrt (rxyzsq);
+
+ spxy = (x * xd) + (y * yd);
+ spxyz = spxy + (z * zd);
+
+ /* Convert back to spherical coordinates */
+ if (x == zero && y == zero)
+ r2000 = zero;
+ else {
+ r2000 = atan2 (y,x);
+ if (r2000 < zero)
+ r2000 = r2000 + d2pi;
+ }
+ d2000 = atan2 (z,rxy);
+
+ if (rxy > tiny) {
+ ur = ((x * yd) - (y * xd)) / rxysq;
+ ud = ((zd * rxysq) - (z * spxy)) / (rxyzsq * rxy);
+ }
+
+ if (*parallax > tiny) {
+ *rv = spxyz / (*parallax * rxyz * vf);
+ *parallax = *parallax / rxyz;
+ }
+
+ /* Return results */
+ *ra = raddeg (r2000);
+ *dec = raddeg (d2000);
+ *rapm = ur / 360000.0;
+ *decpm = ud / 360000.0;
+
+ if (diag) {
+ dra = 240.0 * raddeg (r2000 - r1950);
+ ddec = 3600.0 * raddeg (d2000 - d1950);
+ fprintf(stderr,"J2000-B1950: dra= %11.5f sec ddec= %f11.5f arcsec\n",
+ dra, ddec);
+ }
+ return;
+}
+
+int idg=0;
+
+/* l2,b2 system of galactic coordinates
+ * p = 192.25 ra of galactic north pole (mean b1950.0)
+ * q = 62.6 inclination of galactic to mean b1950.0 equator
+ * r = 33 longitude of ascending node
+ * p,q,r are degrees
+
+ * Equatorial to galactic rotation matrix
+ (The Eulerian angles are p, q, 90-r)
+ +cp.cq.sr-sp.cr +sp.cq.sr+cp.cr -sq.sr
+ -cp.cq.cr-sp.sr -sp.cq.cr+cp.sr +sq.cr
+ cp.sq +sp.sq +cq
+ */
+
+static
+double bgal[3][3] =
+ {{-0.066988739415,-0.872755765852,-0.483538914632},
+ {0.492728466075,-0.450346958020, 0.744584633283},
+ {-0.867600811151,-0.188374601723, 0.460199784784}};
+
+/*--- Transform B1950.0 FK4 equatorial coordinates to
+ * IAU 1958 galactic coordinates */
+
+void
+fk42gal (dtheta,dphi)
+
+double *dtheta; /* B1950.0 FK4 right ascension in degrees
+ Galactic longitude (l2) in degrees (returned) */
+double *dphi; /* B1950.0 FK4 declination in degrees
+ Galactic latitude (b2) in degrees (returned) */
+
+/* Input equatorial coordinates are B1950 FK4.
+ Use fk52gal() to convert from j2000.0 coordinates.
+ Reference: Blaauw et al, MNRAS,121,123 (1960) */
+{
+ double pos[3],pos1[3],r,dl,db,rl,rb,rra,rdec,dra,ddec;
+ void v2s3(),s2v3();
+ int i;
+ char *eqcoor, *eqstrn();
+
+ dra = *dtheta;
+ ddec = *dphi;
+ rra = degrad (dra);
+ rdec = degrad (ddec);
+
+ /* remove e-terms */
+ /* call jpabe (rra,rdec,-1,idg) */
+
+ /* Spherical to Cartesian */
+ r = 1.;
+ s2v3 (rra,rdec,r,pos);
+
+ /* rotate to galactic */
+ for (i = 0; i<3; i++) {
+ pos1[i] = pos[0]*bgal[i][0] + pos[1]*bgal[i][1] + pos[2]*bgal[i][2];
+ }
+
+ /* Cartesian to spherical */
+ v2s3 (pos1,&rl,&rb,&r);
+
+ dl = raddeg (rl);
+ db = raddeg (rb);
+ *dtheta = dl;
+ *dphi = db;
+
+ /* Print result if in diagnostic mode */
+ if (idg) {
+ eqcoor = eqstrn (dra,ddec);
+ fprintf (stderr,"FK42GAL: B1950 RA,Dec= %s\n",eqcoor);
+ fprintf (stderr,"FK42GAL: long = %.5f lat = %.5f\n",dl,db);
+ free (eqcoor);
+ }
+
+ return;
+}
+
+
+/*--- Transform IAU 1958 galactic coordinates to B1950.0 'FK4'
+ * equatorial coordinates */
+
+void
+gal2fk4 (dtheta,dphi)
+
+double *dtheta; /* Galactic longitude (l2) in degrees
+ B1950 FK4 RA in degrees (returned) */
+double *dphi; /* Galactic latitude (b2) in degrees
+ B1950 FK4 Dec in degrees (returned) */
+
+/* Output equatorial coordinates are B1950.0 FK4.
+ Use gal2fk5() to convert to J2000 coordinates.
+ Reference: Blaauw et al, MNRAS,121,123 (1960) */
+
+{
+ double pos[3],pos1[3],r,dl,db,rl,rb,rra,rdec,dra,ddec;
+ void v2s3(),s2v3();
+ char *eqcoor, *eqstrn();
+ int i;
+
+ /* spherical to cartesian */
+ dl = *dtheta;
+ db = *dphi;
+ rl = degrad (dl);
+ rb = degrad (db);
+ r = 1.0;
+ s2v3 (rl,rb,r,pos);
+
+ /* rotate to equatorial coordinates */
+ for (i = 0; i < 3; i++) {
+ pos1[i] = pos[0]*bgal[0][i] + pos[1]*bgal[1][i] + pos[2]*bgal[2][i];
+ }
+
+ /* cartesian to spherical */
+ v2s3 (pos1,&rra,&rdec,&r);
+
+/* introduce e-terms */
+/* jpabe (rra,rdec,-1,idg); */
+
+ dra = raddeg (rra);
+ ddec = raddeg (rdec);
+ *dtheta = dra;
+ *dphi = ddec;
+
+ /* print result if in diagnostic mode */
+ if (idg) {
+ fprintf (stderr,"GAL2FK4: long = %.5f lat = %.5f\n",dl,db);
+ eqcoor = eqstrn (dra,ddec);
+ fprintf (stderr,"GAL2FK4: B1950 RA,Dec= %s\n",eqcoor);
+ free (eqcoor);
+ }
+
+ return;
+}
+
+
+/* l2,b2 system of galactic coordinates
+ p = 192.25 ra of galactic north pole (mean b1950.0)
+ q = 62.6 inclination of galactic to mean b1950.0 equator
+ r = 33 longitude of ascending node
+ p,q,r are degrees */
+
+/* Equatorial to galactic rotation matrix
+ The eulerian angles are p, q, 90-r
+ +cp.cq.sr-sp.cr +sp.cq.sr+cp.cr -sq.sr
+ -cp.cq.cr-sp.sr -sp.cq.cr+cp.sr +sq.cr
+ +cp.sq +sp.sq +cq */
+
+static
+double jgal[3][3] =
+ {{-0.054875539726,-0.873437108010,-0.483834985808},
+ {0.494109453312,-0.444829589425, 0.746982251810},
+ {-0.867666135858,-0.198076386122, 0.455983795705}};
+
+/* Transform J2000 equatorial coordinates to IAU 1958 galactic coordinates */
+
+void
+fk52gal (dtheta,dphi)
+
+double *dtheta; /* J2000 right ascension in degrees
+ Galactic longitude (l2) in degrees (returned) */
+double *dphi; /* J2000 declination in degrees
+ Galactic latitude (b2) in degrees (returned) */
+
+/* Rotation matrices by P.T.Wallace, Starlink eqgal and galeq, March 1986 */
+
+/* Input equatorial coordinates are J2000 FK5.
+ Use gal2fk4() if converting from B1950 FK4 coordinates.
+ Reference: Blaauw et al, MNRAS,121,123 (1960) */
+{
+ double pos[3],pos1[3],r,dl,db,rl,rb,rra,rdec,dra,ddec;
+ void v2s3(),s2v3();
+ char *eqcoor, *eqstrn();
+ int i;
+
+ /* Spherical to cartesian */
+ dra = *dtheta;
+ ddec = *dphi;
+ rra = degrad (dra);
+ rdec = degrad (ddec);
+ r = 1.0;
+ (void)s2v3 (rra,rdec,r,pos);
+
+ /* Rotate to galactic */
+ for (i = 0; i < 3; i++) {
+ pos1[i] = pos[0]*jgal[i][0] + pos[1]*jgal[i][1] + pos[2]*jgal[i][2];
+ }
+
+ /* Cartesian to spherical */
+ v2s3 (pos1,&rl,&rb,&r);
+
+ dl = raddeg (rl);
+ db = raddeg (rb);
+ *dtheta = dl;
+ *dphi = db;
+
+ /* Print result if in diagnostic mode */
+ if (idg) {
+ eqcoor = eqstrn (dra,ddec);
+ fprintf (stderr,"FK52GAL: J2000 RA,Dec= %s\n",eqcoor);
+ fprintf (stderr,"FK52GAL: long = %.5f lat = %.5f\n",dl,db);
+ free (eqcoor);
+ }
+
+ return;
+}
+
+
+/*--- Transform IAU 1958 galactic coordinates to J2000 equatorial coordinates */
+
+void
+gal2fk5 (dtheta,dphi)
+
+double *dtheta; /* Galactic longitude (l2) in degrees
+ J2000.0 ra in degrees (returned) */
+double *dphi; /* Galactic latitude (b2) in degrees
+ J2000.0 dec in degrees (returned) */
+
+/* Output equatorial coordinates are J2000.
+ Use gal2fk4() to convert to B1950 coordinates.
+ Reference: Blaauw et al, MNRAS,121,123 (1960) */
+
+{
+ double pos[3],pos1[3],r,dl,db,rl,rb,rra,rdec,dra,ddec;
+ void v2s3(),s2v3();
+ int i;
+ char *eqcoor, *eqstrn();
+
+ /* Spherical to Cartesian */
+ dl = *dtheta;
+ db = *dphi;
+ rl = degrad (dl);
+ rb = degrad (db);
+ r = 1.0;
+ s2v3 (rl,rb,r,pos);
+
+ /* Rotate to equatorial coordinates */
+ for (i = 0; i < 3; i++) {
+ pos1[i] = pos[0]*jgal[0][i] + pos[1]*jgal[1][i] + pos[2]*jgal[2][i];
+ }
+
+ /* Cartesian to Spherical */
+ v2s3 (pos1,&rra,&rdec,&r);
+ dra = raddeg (rra);
+ ddec = raddeg (rdec);
+ *dtheta = dra;
+ *dphi = ddec;
+
+ /* Print result if in diagnostic mode */
+ if (idg) {
+ fprintf (stderr,"GAL2FK5: long = %.5f lat = %.5f\n",dl,db);
+ eqcoor = eqstrn (dra,ddec);
+ fprintf (stderr,"GAL2FK5: J2000 RA,Dec= %s\n",eqcoor);
+ free (eqcoor);
+ }
+
+ return;
+}
+
+
+/* Return string with right ascension in hours and declination in degrees */
+
+char *eqstrn (dra, ddec)
+
+double dra; /* Right ascension in degrees */
+double ddec; /* Declination in degrees */
+
+{
+char *eqcoor; /* ASCII character string of position (returned) */
+char decp;
+int rah,irm,decd,decm;
+double xpos,ypos,xp,yp,ras,decs;
+
+ /* Right ascension to hours, minutes, and seconds */
+ xpos = dra / 15.0;
+ rah = (int) xpos;
+ xp = (double) 60.0 * (xpos - (double) rah);
+ irm = (int) xp;
+ ras = (double) 60.0 * (xp - (double) irm);
+
+ /* Declination to degrees, minutes, seconds */
+ if (ddec < 0) {
+ ypos = -ddec;
+ decp = '-';
+ }
+ else {
+ decp = '+';
+ ypos = ddec;
+ }
+ decd = (int) ypos;
+ yp = (double) 60.0 * (ypos - (double) decd);
+ decm = (int) yp;
+ decs = (double) 60.0 * (yp - (double) decm);
+
+ eqcoor = malloc (32);
+ (void)sprintf (eqcoor,"%02d:%02d:%06.3f %c%02d:%02d:%05.2f",
+ rah,irm,ras,decp,decd,decm,decs);
+ if (eqcoor[6] == ' ')
+ eqcoor[6] = '0';
+ if (eqcoor[20] == ' ')
+ eqcoor[20] = '0';
+
+ return (eqcoor);
+}
+
+
+/* Convert geocentric equatorial rectangular coordinates to
+ right ascension and declination, and distance */
+
+
+/* These routines are based on similar ones in Pat Wallace's slalib package */
+
+/* Convert B1950 right ascension and declination to ecliptic coordinates */
+
+void
+fk42ecl (dtheta, dphi, epoch)
+
+double *dtheta; /* B1950 right ascension in degrees
+ Galactic longitude (l2) in degrees (returned) */
+double *dphi; /* B1950 declination in degrees
+ Galactic latitude (b2) in degrees (returned) */
+double epoch; /* Besselian epoch in years */
+
+{
+ void fk425e(), fk52ecl();
+
+ /* Convert from B1950 to J2000 coordinates */
+ fk425e (dtheta, dphi, epoch);
+
+ /* Convert from J2000 to ecliptic coordinates */
+ fk52ecl (dtheta, dphi, epoch);
+
+ return;
+}
+
+/* Convert J2000 right ascension and declination to ecliptic coordinates */
+
+void
+fk52ecl (dtheta, dphi, epoch)
+
+double *dtheta; /* J2000 right ascension in degrees
+ Galactic longitude (l2) in degrees (returned) */
+double *dphi; /* J2000 declination in degrees
+ Galactic latitude (b2) in degrees (returned) */
+double epoch; /* Besselian epoch in years */
+
+{
+ int i, j;
+ double t, eps0, rphi, rtheta;
+ double v1[3], v2[3], r;
+ double rmat[9], *rmati; /* Rotation matrix */
+
+ void rotmat(), v2s3(), s2v3(), fk5prec();
+
+ /* Precess coordinates from J2000 to epoch */
+ if (epoch != 2000.0)
+ fk5prec (2000.0, epoch, dtheta, dphi);
+
+ /* Convert from degrees to radians */
+ rtheta = degrad (*dtheta);
+ rphi = degrad (*dphi);
+
+ /* Convert RA,Dec to x,y,z */
+ r = 1.0;
+ s2v3 (rtheta, rphi, r, v1);
+
+ /* Interval between basic epoch J2000.0 and current epoch (JC) in centuries*/
+ t = (epoch - 2000.0) * 0.01;
+
+ /* Mean obliquity */
+ eps0 = secrad ((84381.448 + (-46.8150 + (-0.00059 + 0.001813*t) * t) * t));
+
+ /* Form the equatorial to ecliptic rotation matrix (IAU 1980 theory).
+ * References: Murray, C.A., Vectorial Astrometry, section 4.3.
+ * The matrix is in the sense v[ecl] = rmat * v[equ]; the
+ * equator, equinox and ecliptic are mean of date. */
+ rotmat (1, eps0, 0.0, 0.0, rmat);
+
+ /* Multiply position vector by equatoria to eccliptic rotation matrix */
+ rmati = rmat;
+ for (i = 0; i < 3; i++) {
+ v2[i] = 0;
+ for (j = 0; j < 3; j++)
+ v2[i] = v2[i] + (*rmati++ * v1[j]);
+ }
+
+ /* Convert x,y,z to latitude, longitude */
+ v2s3 (v2, &rtheta, &rphi, &r);
+
+ /* Convert from radians to degrees */
+ *dtheta = raddeg (rtheta);
+ *dphi = raddeg (rphi);
+}
+
+
+/* Convert ecliptic coordinates to B1950 right ascension and declination */
+
+void
+ecl2fk4 (dtheta, dphi, epoch)
+
+double *dtheta; /* Galactic longitude (l2) in degrees
+ B1950 right ascension in degrees (returned) */
+double *dphi; /* Galactic latitude (b2) in degrees
+ B1950 declination in degrees (returned) */
+double epoch; /* Besselian epoch in years */
+
+{
+ void ecl2fk5(), fk524e();
+
+ /* Convert from ecliptic to J2000 coordinates */
+ ecl2fk5 (dtheta, dphi, epoch);
+
+ /* Convert from J2000 to B1950 coordinates */
+ fk524e (dtheta, dphi, epoch);
+
+ return;
+}
+
+
+
+/* Convert ecliptic coordinates to J2000 right ascension and declination */
+
+void
+ecl2fk5 (dtheta, dphi, epoch)
+
+double *dtheta; /* Galactic longitude (l2) in degrees
+ J2000 right ascension in degrees (returned) */
+double *dphi; /* Galactic latitude (b2) in degrees
+ J2000 declination in degrees (returned) */
+double epoch; /* Besselian epoch in years */
+
+{
+ int i, j;
+ double rtheta, rphi, v1[3], v2[3];
+ double t, eps0, r;
+ double rmat[9]; /* Rotation matrix */
+ void v2s3(),s2v3(), fk5prec(), rotmat();
+
+ rtheta = degrad (*dtheta);
+ rphi = degrad (*dphi);
+
+ /* Convert RA,Dec to x,y,z */
+ r = 1.0;
+ s2v3 (rtheta, rphi, r, v1);
+
+ /* Interval between basic epoch J2000.0 and current epoch (JC) in centuries*/
+ t = (epoch - 2000.0) * 0.01;
+
+ /* Mean obliquity */
+ eps0 = secrad ((84381.448 + (-46.8150 + (-0.00059 + 0.001813*t) * t) * t));
+
+ /* Form the equatorial to ecliptic rotation matrix (IAU 1980 theory).
+ * References: Murray, C.A., Vectorial Astrometry, section 4.3.
+ * The matrix is in the sense v[ecl] = rmat * v[equ]; the
+ * equator, equinox and ecliptic are mean of date. */
+ rotmat (1, eps0, 0.0, 0.0, rmat);
+
+ /* Multiply position vector by ecliptic to equatorial rotation matrix */
+ for (i = 0; i < 3; i++) {
+ v2[i] = 0;
+ for (j = 0; j < 3; j++)
+ v2[i] = v2[i] + (rmat[3*j + i] * v1[j]);
+ }
+
+ /* Cartesian to spherical */
+ v2s3 (v2, &rtheta, &rphi, &r);
+
+ /* Convert from radians to degrees */
+ *dtheta = raddeg (rtheta);
+ *dphi = raddeg (rphi);
+
+ if (epoch != 2000.0)
+ fk5prec (epoch, 2000.0, dtheta, dphi);
+}
+
+
+/* The following routines are modified from Patrick Wallace's SLALIB */
+
+/* Precess coordinates between epochs in FK4 */
+void
+fk4prec (ep0, ep1, ra, dec)
+
+double ep0; /* Starting Besselian epoch */
+double ep1; /* Ending Besselian epoch */
+double *ra; /* RA in degrees mean equator & equinox of epoch ep0
+ mean equator & equinox of epoch ep1 (returned) */
+double *dec; /* Dec in degrees mean equator & equinox of epoch ep0
+ mean equator & equinox of epoch ep1 (returned) */
+/*
+** Precession - FK4 (Bessel-Newcomb, pre-IAU1976)
+**
+** This routine will not correctly convert between FK4 and FK5
+** For output in FK5, precess to 1950.0 and use fk425() on result.
+**
+** Based on slaPreces(), P.T.Wallace Starlink 22 December 1993
+*/
+{
+ int i, j;
+ double pm[9], *pmi, v1[3], v2[3], rra, rdec, r;
+ void v2s3(),s2v3(), mprecfk4();
+
+ rra = degrad (*ra);
+ rdec = degrad (*dec);
+ r = 1.0;
+
+ /* Generate appropriate precession matrix */
+ mprecfk4 ( ep0, ep1, pm );
+
+ /* Convert RA,Dec to x,y,z */
+ s2v3 (rra, rdec, r, v1);
+
+ /* Multiply position vector by precession matrix */
+ pmi = pm;
+ for (i = 0; i < 3; i++) {
+ v2[i] = 0;
+ for (j = 0; j < 3; j++)
+ v2[i] = v2[i] + (*pmi++ * v1[j]);
+ }
+
+ /* Back to RA,Dec */
+ v2s3 (v2, &rra, &rdec, &r);
+
+ /* Convert from radians to degrees */
+ *ra = raddeg (rra);
+ *dec = raddeg (rdec);
+}
+
+void
+fk5prec (ep0, ep1, ra, dec)
+
+double ep0; /* Starting epoch */
+double ep1; /* Ending epoch */
+double *ra; /* RA in degrees mean equator & equinox of epoch ep0
+ mean equator & equinox of epoch ep1 (returned) */
+double *dec; /* Dec in degrees mean equator & equinox of epoch ep0
+ mean equator & equinox of epoch ep1 (returned) */
+/*
+** Precession - FK5 (Fricke, post-IAU1976)
+**
+** This routine will not correctly convert between FK5 and FK4.
+** For output in FK4, precess to 2000.0 and use fk524() on result.
+**
+** Based on slaPreces(), P.T.Wallace Starlink 22 December 1993
+*/
+{
+ int i, j;
+ double pm[9], *pmi, v1[3], v2[3], rra, rdec, r;
+ void v2s3(),s2v3(), mprecfk5();
+
+ rra = degrad (*ra);
+ rdec = degrad (*dec);
+ r = 1.0;
+
+ /* Generate appropriate precession matrix */
+ mprecfk5 (ep0, ep1, pm);
+
+ /* Convert RA,Dec to x,y,z */
+ s2v3 (rra, rdec, r, v1);
+
+ /* Multiply position vector by precession matrix */
+ pmi = pm;
+ for (i = 0; i < 3; i++) {
+ v2[i] = 0;
+ for (j = 0; j < 3; j++)
+ v2[i] = v2[i] + ( v1[j] * *pmi++ );
+ }
+
+ /* Back to RA,Dec */
+ v2s3 (v2, &rra, &rdec, &r);
+
+ /* Convert from radians to degrees */
+ *ra = raddeg (rra);
+ *dec = raddeg (rdec);
+ return;
+}
+
+
+void
+mprecfk4 (bep0, bep1, rmatp)
+
+double bep0; /* Beginning Besselian epoch */
+double bep1; /* Ending Besselian epoch */
+double rmatp[9]; /* 3x3 Precession matrix (returned) */
+
+/*
+** Generate the matrix of precession between two epochs,
+** using the old, pre-IAU1976, Bessel-Newcomb model, using
+** Kinoshita's formulation (double precision)
+**
+** The matrix is in the sense v(bep1) = rmatp * v(bep0)
+**
+** Reference:
+** Kinoshita, H. (1975) 'Formulas for precession', SAO Special
+** Report No. 364, Smithsonian Institution Astrophysical
+** Observatory, Cambridge, Massachusetts.
+**
+** Based on slaPrebn() by P.T.Wallace Starlink 30 October 1993
+*/
+{
+ double bigt, t, tas2r, w, zeta, z, theta;
+ void rotmat();
+
+ /* Interval between basic epoch B1850.0 and beginning epoch in TC */
+ bigt = ( bep0 - 1850.0 ) / 100.0;
+
+ /* Interval over which precession required, in tropical centuries */
+ t = ( bep1 - bep0 ) / 100.0;
+
+ /* Euler angles */
+ tas2r = secrad (t);
+ w = 2303.5548 + ( 1.39720 + 0.000059 * bigt ) * bigt;
+ zeta = (w + ( 0.30242 - 0.000269 * bigt + 0.017996 * t ) * t ) * tas2r;
+ z = (w + ( 1.09478 + 0.000387 * bigt + 0.018324 * t ) * t ) * tas2r;
+ theta = ( 2005.1125 + ( - 0.85294 - 0.000365* bigt ) * bigt +
+ ( - 0.42647 - 0.000365 * bigt - 0.041802 * t ) * t ) * tas2r;
+
+ /* Rotation matrix */
+ rotmat (323, -zeta, theta, -z, rmatp);
+ return;
+}
+
+
+void
+mprecfk5 (ep0, ep1, rmatp)
+
+double ep0; /* Beginning epoch */
+double ep1; /* Ending epoch */
+double rmatp[9]; /* 3x3 Precession matrix (returned) */
+
+/*
+** Form the matrix of precession between two epochs (IAU 1976, FK5).
+** Notes:
+** 1) The epochs are TDB (loosely ET) Julian epochs.
+** 2) The matrix is in the sense v(ep1) = rmatp * v(ep0) .
+**
+** References:
+** Lieske,J.H., 1979. Astron. Astrophys.,73,282.
+** equations (6) & (7), p283.
+** Kaplan,G.H., 1981. USNO circular no. 163, pa2.
+**
+** Based on slaPrec(), P.T.Wallace Starlink 31 October 1993
+*/
+{
+ double t0, t, tas2r, w, zeta, z, theta;
+ void rotmat();
+
+ /* Interval between basic epoch J2000.0 and beginning epoch (JC) */
+ t0 = ( ep0 - 2000.0 ) / 100.0;
+
+ /* Interval over which precession required (JC) */
+ t = ( ep1 - ep0 ) / 100.0;
+
+ /* Euler angles */
+ tas2r = secrad (t);
+ w = 2306.2181 + ( ( 1.39656 - ( 0.000139 * t0 ) ) * t0 );
+ zeta = (w + ( ( 0.30188 - 0.000344 * t0 ) + 0.017998 * t ) * t ) * tas2r;
+ z = (w + ( ( 1.09468 + 0.000066 * t0 ) + 0.018203 * t ) * t ) * tas2r;
+ theta = ( ( 2004.3109 + ( - 0.85330 - 0.000217 * t0 ) * t0 )
+ + ( ( -0.42665 - 0.000217 * t0 ) - 0.041833 * t ) * t ) * tas2r;
+
+ /* Rotation matrix */
+ rotmat (323, -zeta, theta, -z, rmatp);
+ return;
+}
+
+
+/* Make 3-D rotation matrix from up to three rotations */
+
+void
+rotmat (axes, rot1, rot2, rot3, matrix)
+
+int axes; /* Axes about which coordinates are rotated (1=x, 2=y, 3=z) */
+double rot1; /* First rotation in degrees */
+double rot2; /* Second rotation in degrees */
+double rot3; /* Third rotation in degrees */
+double *matrix; /* 3x3 rotation matrix (returned) */
+
+{
+ int i, j, k, naxis, iaxes, iaxis;
+ double rot[3], srot, crot, *mati, w, wm[9], *wmi, matn[9];
+ int axis[3];
+
+ /* Initial final rotation matrix */
+ mati = matrix;
+ for (i = 0; i < 3; i++) {
+ for (j=0; j < 3; j++) {
+ if (i == j)
+ *mati++ = 1.0;
+ else
+ *mati++ = 0.0;
+ }
+ }
+
+ /* Separate digits of rotation axis string and count rotations */
+ naxis = 0;
+ iaxes = axes;
+ axis[0] = iaxes / 100;
+ if (axis[0] > 0) {
+ naxis++;
+ iaxes = iaxes - (100 * axis[0]);
+ }
+ axis[naxis] = iaxes / 10;
+ if (axis[naxis] > 0) {
+ iaxes = iaxes - (10 * axis[naxis]);
+ naxis++;
+ }
+ axis[naxis] = iaxes;
+ if (axis[naxis] > 0)
+ naxis++;
+
+ /* Set up rotation angles */
+ rot[0] = rot1;
+ rot[1] = rot2;
+ rot[2] = rot3;
+
+ /* For each digit of axis string, set up matrix */
+ for (iaxis = 0; iaxis < naxis; iaxis++) {
+
+ /* Initialize current rotation matrix */
+ mati = matn;
+ for (i = 0; i < 3; i++) {
+ for (j=0; j < 3; j++) {
+ if (i == j)
+ *mati++ = 1.0;
+ else
+ *mati++ = 0.0;
+ }
+ }
+
+ srot = sin (rot[iaxis]);
+ crot = cos (rot[iaxis]);
+
+ /* Matrix for rotation in X */
+ if (axis[iaxis] == 1) {
+ matn[4] = crot;
+ matn[5] = srot;
+ matn[7] = -srot;
+ matn[8] = crot;
+ }
+
+ /* Matrix for rotation in Y */
+ else if (axis[iaxis] == 2) {
+ matn[0] = crot;
+ matn[2] = -srot;
+ matn[6] = srot;
+ matn[8] = crot;
+ }
+
+ /* Matrix for rotation in Z */
+ else {
+ matn[0] = crot;
+ matn[1] = srot;
+ matn[3] = -srot;
+ matn[4] = crot;
+ }
+
+ /* Multiply existing rotation matrix by new rotation matrix */
+ for (i = 0; i < 3; i++) {
+ for (j = 0; j < 3; j++) {
+ w = 0.0;
+ for (k = 0; k < 3; k++)
+ w+= matn[3*i + k] * matrix[3*k + j];
+ wm[3*i + j] = w;
+ }
+ }
+
+ /* Update output matrix */
+ mati = matrix;
+ wmi = wm;
+ for (i = 0; i < 9; i++) {
+ *mati++ = *wmi++;
+ }
+ }
+ return;
+}
+
+
+/* The following routines are from Doug Mink's Fortran ephemeris library */
+
+/* Convert right ascensiona and declination in degrees and distance to
+ geocentric equatorial rectangular coordinates */
+
+void
+d2v3 (rra,rdec,r,pos)
+
+double rra; /* Right ascension in degrees */
+double rdec; /* Declination in degrees */
+double r; /* Distance to object in same units as pos */
+double pos[3]; /* x,y,z geocentric equatorial position of object (returned) */
+{
+ s2v3 (degrad (rra), degrad (rdec), r, pos);
+
+ return;
+}
+
+
+/* Convert right ascension, declination, and distance to
+ geocentric equatorial rectangular coordinates */
+
+void
+s2v3 (rra,rdec,r,pos)
+
+double rra; /* Right ascension in radians */
+double rdec; /* Declination in radians */
+double r; /* Distance to object in same units as pos */
+double pos[3]; /* x,y,z geocentric equatorial position of object (returned) */
+{
+ pos[0] = r * cos (rra) * cos (rdec);
+ pos[1] = r * sin (rra) * cos (rdec);
+ pos[2] = r * sin (rdec);
+
+ return;
+}
+
+
+/* Convert geocentric equatorial rectangular coordinates to
+ right ascension and declination in degrees and distance */
+
+void
+v2d3 (pos,rra,rdec,r)
+
+double pos[3]; /* x,y,z geocentric equatorial position of object */
+double *rra; /* Right ascension in degrees (returned) */
+double *rdec; /* Declination in degrees (returned) */
+double *r; /* Distance to object in same units as pos (returned) */
+{
+ v2s3 (pos, rra, rdec, r);
+ *rra = raddeg (*rra);
+ *rdec = raddeg (*rdec);
+ return;
+}
+
+/* Convert geocentric equatorial rectangular coordinates to
+ right ascension, declination, and distance */
+
+void
+v2s3 (pos,rra,rdec,r)
+
+double pos[3]; /* x,y,z geocentric equatorial position of object */
+double *rra; /* Right ascension in radians (returned) */
+double *rdec; /* Declination in radians (returned) */
+double *r; /* Distance to object in same units as pos (returned) */
+{
+ double x,y,z,rxy,rxy2,z2;
+
+ x = pos[0];
+ y = pos[1];
+ z = pos[2];
+
+ *rra = atan2 (y, x);
+
+ /* Keep RA within 0 to 2pi range */
+ if (*rra < 0.0)
+ *rra = *rra + (2.0 * PI);
+ if (*rra > 2.0 * PI)
+ *rra = *rra - (2.0 * PI);
+
+ rxy2 = x*x + y*y;
+ rxy = sqrt (rxy2);
+ *rdec = atan2 (z, rxy);
+
+ z2 = z * z;
+ *r = sqrt (rxy2 + z2);
+
+ return;
+}
+
+/*
+ * Nov 6 1995 Include stdlib.h instead of malloc.h
+ * Apr 1 1996 Add arbitrary epoch precession
+ * Apr 26 1996 Add FK4 <-> FK5 subroutines for use when epoch is known
+ * Aug 6 1996 Clean up after lint
+ * Nov 4 1996 Break SLA subroutines into separate file slasubs.c
+ * Dec 9 1996 Change arguments to degrees in FK4 and FK5 precession programs
+ * Dec 10 1996 All subroutine arguments are degrees except vector conversions
+ *
+ * Mar 20 1997 Drop unused variables after lint
+ *
+ * Apr 14 1998 Add ecliptic coordinate conversions and general conversion routines
+ * Apr 23 1998 Add LINEAR coordinate system
+ * Apr 28 1998 Change coordinate system flags to WCS_*
+ * Apr 28 1998 Return -1 from wcscsys if not a legal coordinate system
+ * May 7 1998 Keep theta within 0 to 2pi in ecl2fk5()
+ * May 13 1998 Add wcsceq()
+ * May 13 1998 Add equinox arguments to wcscon()
+ * Jun 24 1998 Set J2000 from ICRS in wcscsys()
+ * Jul 9 1998 Include stdio.h for fprintf() and sprintf() declarations
+ * Sep 17 1998 Add wcscstr() to get coordinate string
+ * Sep 21 1998 Fix bug in wcscstr() which returned B2000 instead of J2000
+ * Sep 21 1998 Add subroutine to convert proper motions, too.
+ * Oct 21 1998 In wcscstr(), drop .00 from returned string
+ * Nov 18 1998 Rename jpcop() v2s3() and jpcon() s2v3() (spherical to vector)
+ * Dec 2 1998 Add PLANET coordinate system to wcscsys() and wcscstr()
+ *
+ * Mar 10 2000 Precess coordinates correctly from other than 1950.0 and 2000.0
+ * Mar 10 2000 Set coordinate system to J2000 or B1950 if string is numeric
+ * Mar 14 2000 Clean up code in fk524m() and fk425m()
+ * May 31 2000 Add proper motion correctly if proper motion precessed
+ * Jun 26 2000 Add some support for WCS_XY image coordinates
+ * Sep 14 2000 Return -1 from wcscsys if equinox is less than 1900.0
+ * Oct 31 2000 Add proper motion after fk425 or fk524 from system epoch
+ * Oct 31 2000 Fix proper motion units in fk524p() and fk425p()
+ * Nov 6 2000 Update fk425 and fk524 algorithms to include parallax and rv
+ *
+ * Jan 11 2001 Print all messages to stderr
+ * Mar 21 2001 Move braces around bgal[] and jgal[] matrix initialization
+ *
+ * Feb 13 2002 Fix precession units problem in ecl2fk5() and fk52ecl()
+ *
+ * Apr 13 2005 Replace all sla_lib calls with local code
+ * Nov 1 2005 Add WCS_ICRS, and unprecessable system
+ *
+ * Jan 5 2006 Fix bugs in precession subroutines mprecxxx()
+ * May 3 2006 Drop declarations of unused variables suggested by Robert Lupton
+ * Oct 6 2006 If pixel coordinates, set system to WCS_XY in wcscsys()
+ * Oct 30 2006 Add LINEAR and ICRS to wcscstr() returns
+ *
+ * Aug 15 2007 Clean up code in rotmat()
+ * Nov 8 2007 In wcsconp, make it clear that proper motion is in spherical coordinates
+ *
+ * Mar 29 2010 Fix bug in computing the magnitude of the e-terms in fk524()
+ * Mar 30 2010 Drop ep1 assignment after line 178 in wcsconp()
+ */