diff options
Diffstat (limited to 'funtools/wcs/wcscon.c')
| -rw-r--r-- | funtools/wcs/wcscon.c | 2328 |
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, ¶llax, &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, ¶llax, &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() + */ |