diff options
Diffstat (limited to 'wcs/sph.c')
| -rw-r--r-- | wcs/sph.c | 234 |
1 files changed, 234 insertions, 0 deletions
diff --git a/wcs/sph.c b/wcs/sph.c new file mode 100644 index 0000000..b8ba23d --- /dev/null +++ b/wcs/sph.c @@ -0,0 +1,234 @@ +/*============================================================================ +* +* WCSLIB - an implementation of the FITS WCS proposal. +* Copyright (C) 1995-2002, Mark Calabretta +* +* 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 WCSLIB may be directed to: +* Internet email: mcalabre@atnf.csiro.au +* Postal address: Dr. Mark Calabretta, +* Australia Telescope National Facility, +* P.O. Box 76, +* Epping, NSW, 2121, +* AUSTRALIA +* +*============================================================================= +* +* C routines for the spherical coordinate transformations used by the FITS +* "World Coordinate System" (WCS) convention. +* +* Summary of routines +* ------------------- +* The spherical coordinate transformations are implemented via separate +* functions for the transformation in each direction. +* +* Forward transformation; sphfwd() +* -------------------------------- +* Transform celestial coordinates to the native coordinates of a projection. +* +* Given: +* lng,lat double Celestial longitude and latitude, in degrees. +* eul[5] double Euler angles for the transformation: +* 0: Celestial longitude of the native pole, in +* degrees. +* 1: Celestial colatitude of the native pole, or +* native colatitude of the celestial pole, in +* degrees. +* 2: Native longitude of the celestial pole, in +* degrees. +* 3: cos(eul[1]) +* 4: sin(eul[1]) +* +* Returned: +* phi, double Longitude and latitude in the native coordinate +* theta system of the projection, in degrees. +* +* Function return value: +* int Error status +* 0: Success. +* +* Reverse transformation; sphrev() +* -------------------------------- +* Transform native coordinates of a projection to celestial coordinates. +* +* Given: +* phi, double Longitude and latitude in the native coordinate +* theta system of the projection, in degrees. +* eul[5] double Euler angles for the transformation: +* 0: Celestial longitude of the native pole, in +* degrees. +* 1: Celestial colatitude of the native pole, or +* native colatitude of the celestial pole, in +* degrees. +* 2: Native longitude of the celestial pole, in +* degrees. +* 3: cos(eul[1]) +* 4: sin(eul[1]) +* +* Returned: +* lng,lat double Celestial longitude and latitude, in degrees. +* +* Function return value: +* int Error status +* 0: Success. +* +* Author: Mark Calabretta, Australia Telescope National Facility +* $Id: sph.c,v 2.7 2002/04/03 01:25:29 mcalabre Exp $ +*===========================================================================*/ + +#include <math.h> +#include "wcslib.h" + +#ifndef __STDC__ +#ifndef const +#define const +#endif +#endif + +const double tol = 1.0e-5; + +int sphfwd (lng, lat, eul, phi, theta) + +const double lat, lng, eul[5]; +double *phi, *theta; + +{ + double coslat, coslng, dlng, dphi, sinlat, sinlng, x, y, z; + + coslat = cosdeg (lat); + sinlat = sindeg (lat); + + dlng = lng - eul[0]; + coslng = cosdeg (dlng); + sinlng = sindeg (dlng); + + /* Compute the native longitude. */ + x = sinlat*eul[4] - coslat*eul[3]*coslng; + if (fabs(x) < tol) { + /* Rearrange formula to reduce roundoff errors. */ + x = -cosdeg (lat+eul[1]) + coslat*eul[3]*(1.0 - coslng); + } + y = -coslat*sinlng; + if (x != 0.0 || y != 0.0) { + dphi = atan2deg (y, x); + } else { + /* Change of origin of longitude. */ + dphi = dlng - 180.0; + } + *phi = eul[2] + dphi; + + /* Normalize the native longitude. */ + if (*phi > 180.0) { + *phi -= 360.0; + } else if (*phi < -180.0) { + *phi += 360.0; + } + + /* Compute the native latitude. */ + if (fmod(dlng,180.0) == 0.0) { + *theta = lat + coslng*eul[1]; + if (*theta > 90.0) *theta = 180.0 - *theta; + if (*theta < -90.0) *theta = -180.0 - *theta; + } else { + z = sinlat*eul[3] + coslat*eul[4]*coslng; + /* Use an alternative formula for greater numerical accuracy. */ + if (fabs(z) > 0.99) { + if (z < 0) + *theta = -acosdeg (sqrt(x*x+y*y)); + else + *theta = acosdeg (sqrt(x*x+y*y)); + } else { + *theta = asindeg (z); + } + } + + return 0; +} + +/*-----------------------------------------------------------------------*/ + +int sphrev (phi, theta, eul, lng, lat) + +const double phi, theta, eul[5]; +double *lng, *lat; + +{ + double cosphi, costhe, dlng, dphi, sinphi, sinthe, x, y, z; + + costhe = cosdeg (theta); + sinthe = sindeg (theta); + + dphi = phi - eul[2]; + cosphi = cosdeg (dphi); + sinphi = sindeg (dphi); + + /* Compute the celestial longitude. */ + x = sinthe*eul[4] - costhe*eul[3]*cosphi; + if (fabs(x) < tol) { + /* Rearrange formula to reduce roundoff errors. */ + x = -cosdeg (theta+eul[1]) + costhe*eul[3]*(1.0 - cosphi); + } + y = -costhe*sinphi; + if (x != 0.0 || y != 0.0) { + dlng = atan2deg (y, x); + } else { + /* Change of origin of longitude. */ + dlng = dphi + 180.0; + } + *lng = eul[0] + dlng; + + /* Normalize the celestial longitude. */ + if (eul[0] >= 0.0) { + if (*lng < 0.0) *lng += 360.0; + } else { + if (*lng > 0.0) *lng -= 360.0; + } + + if (*lng > 360.0) { + *lng -= 360.0; + } else if (*lng < -360.0) { + *lng += 360.0; + } + + /* Compute the celestial latitude. */ + if (fmod(dphi,180.0) == 0.0) { + *lat = theta + cosphi*eul[1]; + if (*lat > 90.0) *lat = 180.0 - *lat; + if (*lat < -90.0) *lat = -180.0 - *lat; + } else { + z = sinthe*eul[3] + costhe*eul[4]*cosphi; + + /* Use an alternative formula for greater numerical accuracy. */ + if (fabs(z) > 0.99) { + if (z < 0) + *lat = -acosdeg (sqrt(x*x+y*y)); + else + *lat = acosdeg (sqrt(x*x+y*y)); + } else { + *lat = asindeg (z); + } + } + + return 0; +} +/* Dec 20 1999 Doug Mink - Change cosd() and sind() to cosdeg() and sindeg() + * Dec 20 1999 Doug Mink - Include wcslib.h, which includes wcstrig.h, sph.h + * Dec 20 1999 Doug Mink - Define copysign only if it is not already defined + * + * Jan 5 2000 Doug Mink - Drop copysign + * + * Sep 19 2001 Doug Mink - No change for WCSLIB 2.7 + */ |