/*
* sph.c
*
* Spherical coordinate transformations.
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* This file part of: AstrOmatic WCS library
*
* Copyright: (C) 2000-2010 IAP/CNRS/UPMC
* (C) 1995-1999 Mark Calabretta
*
* Authors: Emmanuel Bertin (this version)
* Mark Calabretta (original version)
*
* Licenses: GNU General Public License
*
* AstrOmatic software is free software: you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3 of the
* License, or (at your option) any later version.
* AstrOmatic software 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 General Public License for more details.
* You should have received a copy of the GNU General Public License
* along with AstrOmatic software.
* If not, see .
*
* Last modified: 10/10/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/*============================================================================
*
* WCSLIB - an implementation of the FITS WCS proposal.
* Copyright (C) 1995-1999, Mark Calabretta
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Library 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 Library
* General Public License for more details.
*
* You should have received a copy of the GNU Library General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 675 Mass Ave, Cambridge, MA 02139, 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 1.1.1.1 2002/03/15 16:33:26 bertin Exp $
*===========================================================================*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_MATHIMF_H
#include
#else
#include
#endif
#include "wcstrig.h"
#include "sph.h"
#ifndef __STDC__
#ifndef const
#define const
#endif
#endif
#define wcs_copysign(X, Y) ((Y) < 0.0 ? -fabs(X) : fabs(X))
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 = wcs_cosd(lat);
sinlat = wcs_sind(lat);
dlng = lng - eul[0];
coslng = wcs_cosd(dlng);
sinlng = wcs_sind(dlng);
/* Compute the native longitude. */
x = sinlat*eul[4] - coslat*eul[3]*coslng;
if (fabs(x) < tol) {
/* Rearrange formula to reduce roundoff errors. */
x = -wcs_cosd(lat+eul[1]) + coslat*eul[3]*(1.0 - coslng);
}
y = -coslat*sinlng;
if (x != 0.0 || y != 0.0) {
dphi = wcs_atan2d(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;
if (fabs(z) > 0.99) {
/* Use an alternative formula for greater numerical accuracy. */
*theta = wcs_copysign(wcs_acosd(sqrt(x*x+y*y)), z);
} else {
*theta = wcs_asind(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 = wcs_cosd(theta);
sinthe = wcs_sind(theta);
dphi = phi - eul[2];
cosphi = wcs_cosd(dphi);
sinphi = wcs_sind(dphi);
/* Compute the celestial longitude. */
x = sinthe*eul[4] - costhe*eul[3]*cosphi;
if (fabs(x) < tol) {
/* Rearrange formula to reduce roundoff errors. */
x = -wcs_cosd(theta+eul[1]) + costhe*eul[3]*(1.0 - cosphi);
}
y = -costhe*sinphi;
if (x != 0.0 || y != 0.0) {
dlng = wcs_atan2d(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;
if (fabs(z) > 0.99) {
/* Use an alternative formula for greater numerical accuracy. */
*lat = wcs_copysign(wcs_acosd(sqrt(x*x+y*y)), z);
} else {
*lat = wcs_asind(z);
}
}
return 0;
}