/*** File wcscon.c
 *** June 9, 2016
 *** Doug Mink, Harvard-Smithsonian Center for Astrophysics
 *** Some subroutines are based on Starlink subroutines by Patrick Wallace
 *** Copyright (C) 1995-2016
 *** 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) == 1 || isnum (wcstring) == 2) {
	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()
 *
 * Jun  9 2016	Fix isnum() tests for added coloned times and dashed dates
 */