Merge commit '23c7930d27fe11c4655e1291a07a821dbbaba78d' as 'funtools'

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()
+ */