path: root/funtools/wcs/cel.c

/*=============================================================================
*
*   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 which implement the FITS World Coordinate System (WCS)
*   convention.
*
*   Summary of routines
*   -------------------
*   These routines are provided as drivers for the lower level spherical
*   coordinate transformation and projection routines.  There are separate
*   driver routines for the forward, celfwd(), and reverse, celrev(),
*   transformations.
*
*   An initialization routine, celset(), computes intermediate values from
*   the transformation parameters but need not be called explicitly - see the
*   explanation of cel.flag below.
*
*
*   Initialization routine; celset()
*   --------------------------------
*   Initializes members of a celprm data structure which hold intermediate
*   values.  Note that this routine need not be called directly; it will be
*   invoked by celfwd() and celrev() if the "flag" structure member is
*   anything other than a predefined magic value.
*
*   Given:
*      pcode[4] const char
*                        WCS projection code (see below).
*
*   Given and returned:
*      cel      celprm*  Spherical coordinate transformation parameters
*                        (see below).
*      prj      prjprm*  Projection parameters (usage is described in the
*                        prologue to "proj.c").
*
*   Function return value:
*               int      Error status
*                           0: Success.
*                           1: Invalid coordinate transformation parameters.
*                           2: Ill-conditioned coordinate transformation
*                              parameters.
*
*   Forward transformation; celfwd()
*   --------------------------------
*   Compute (x,y) coordinates in the plane of projection from celestial
*   coordinates (lng,lat).
*
*   Given:
*      pcode[4] const char
*                        WCS projection code (see below).
*      lng,lat  const double
*                        Celestial longitude and latitude of the projected
*                        point, in degrees.
*
*   Given and returned:
*      cel      celprm*  Spherical coordinate transformation parameters
*                        (see below).
*
*   Returned:
*      phi,     double*  Longitude and latitude in the native coordinate
*      theta             system of the projection, in degrees.
*
*   Given and returned:
*      prj      prjprm*  Projection parameters (usage is described in the
*                        prologue to "proj.c").
*
*   Returned:
*      x,y      double*  Projected coordinates, "degrees".
*
*   Function return value:
*               int      Error status
*                           0: Success.
*                           1: Invalid coordinate transformation parameters.
*                           2: Invalid projection parameters.
*                           3: Invalid value of (lng,lat).
*
*   Reverse transformation; celrev()
*   --------------------------------
*   Compute the celestial coordinates (lng,lat) of the point with projected
*   coordinates (x,y).
*
*   Given:
*      pcode[4] const char
*                        WCS projection code (see below).
*      x,y      const double
*                        Projected coordinates, "degrees".
*
*   Given and returned:
*      prj      prjprm*  Projection parameters (usage is described in the
*                        prologue to "proj.c").
*
*   Returned:
*      phi,     double*  Longitude and latitude in the native coordinate
*      theta             system of the projection, in degrees.
*
*   Given and returned:
*      cel      celprm*  Spherical coordinate transformation parameters
*                        (see below).
*
*   Returned:
*      lng,lat  double*  Celestial longitude and latitude of the projected
*                        point, in degrees.
*
*   Function return value:
*               int      Error status
*                           0: Success.
*                           1: Invalid coordinate transformation parameters.
*                           2: Invalid projection parameters.
*                           3: Invalid value of (x,y).
*
*   Coordinate transformation parameters
*   ------------------------------------
*   The celprm struct consists of the following:
*
*      int flag
*         The celprm struct contains pointers to the forward and reverse
*         projection routines as well as intermediaries computed from the
*         reference coordinates (see below).  Whenever the projection code
*         (pcode) or any of ref[4] are set or changed then this flag must be
*         set to zero to signal the initialization routine, celset(), to
*         redetermine the function pointers and recompute intermediaries.
*         Once this has been done pcode itself is ignored.
*
*      double ref[4]
*         The first pair of values should be set to the celestial longitude
*         and latitude (usually right ascension and declination) of the
*         reference point of the projection.  These are given by the CRVALn
*         keywords in FITS.
*
*         The second pair of values are the native longitude of the celestial
*         pole and the celestial latitude of the native pole and correspond to
*         FITS keywords LONPOLE and LATPOLE.
*
*         LONPOLE defaults to 0 degrees if the celestial latitude of the
*         reference point of the projection is greater than the native
*         latitude, otherwise 180 degrees.  (This is the condition for the
*         celestial latitude to increase in the same direction as the native
*         latitude at the reference point.)  ref[2] may be set to 999.0 to
*         indicate that the correct default should be substituted.
*
*         In some circumstances the celestial latitude of the native pole may
*         be determined by the first three values only to within a sign and
*         LATPOLE is used to choose between the two solutions.  LATPOLE is
*         set in ref[3] and the solution closest to this value is used to
*         reset ref[3].  It is therefore legitimate, for example, to set
*         ref[3] to 999.0 to choose the more northerly solution - the default
*         if the LATPOLE card is omitted from the FITS header.  For the
*         special case where the reference point of the projection is at
*         native latitude zero, its celestial latitude is zero, and
*         LONPOLE = +/- 90 then the celestial latitude of the pole is not
*         determined by the first three reference values and LATPOLE
*         specifies it completely.
*
*   The remaining members of the celprm struct are maintained by the
*   initialization routines and should not be modified.  This is done for the
*   sake of efficiency and to allow an arbitrary number of contexts to be
*   maintained simultaneously.
*
*      double euler[5]
*         Euler angles and associated intermediaries derived from the
*         coordinate reference values.
*
*
*   WCS projection codes
*   --------------------
*   Zenithals/azimuthals:
*      AZP: zenithal/azimuthal perspective
*      TAN: gnomonic
*      STG: stereographic
*      SIN: synthesis (generalized orthographic)
*      ARC: zenithal/azimuthal equidistant
*      ZPN: zenithal/azimuthal polynomial
*      ZEA: zenithal/azimuthal equal area
*      AIR: Airy
*
*   Cylindricals:
*      CYP: cylindrical perspective
*      CEA: cylindrical equal area
*      CAR: Cartesian
*      MER: Mercator
*
*   Pseudo-cylindricals:
*      SFL: Sanson-Flamsteed
*      PAR: parabolic
*      MOL: Mollweide
*
*   Conventional:
*      AIT: Hammer-Aitoff
*
*   Conics:
*      COP: conic perspective
*      COD: conic equidistant
*      COE: conic equal area
*      COO: conic orthomorphic
*
*   Polyconics:
*      BON: Bonne
*      PCO: polyconic
*
*   Quad-cubes:
*      TSC: tangential spherical cube
*      CSC: COBE quadrilateralized spherical cube
*      QSC: quadrilateralized spherical cube
*
*   Author: Mark Calabretta, Australia Telescope National Facility
*   $Id: cel.c,v 2.14 2002/04/03 01:25:29 mcalabre Exp $
*===========================================================================*/

#include <math.h>
#include <string.h>
#include "wcslib.h"

/* Map error number to error message for each function. */
const char *celset_errmsg[] = {
   0,
   "Invalid coordinate transformation parameters",
   "Ill-conditioned coordinate transformation parameters"};

const char *celfwd_errmsg[] = {
   0,
   "Invalid coordinate transformation parameters",
   "Invalid projection parameters",
   "Invalid value of (lng,lat)"};

const char *celrev_errmsg[] = {
   0,
   "Invalid coordinate transformation parameters",
   "Invalid projection parameters",
   "Invalid value of (x,y)"};
 

int
celset(pcode, cel, prj)

const char pcode[4];
struct celprm *cel;
struct prjprm *prj;

{
   int dophip;
   const double tol = 1.0e-10;
   double clat0, cphip, cthe0, slat0, sphip, sthe0;
   double latp, latp1, latp2;
   double u, v, x, y, z;

   /* Initialize the projection driver routines. */
   if (prjset(pcode, prj)) {
      return 1;
   }

   /* Set default for native longitude of the celestial pole? */
   dophip = (cel->ref[2] == 999.0);

   /* Compute celestial coordinates of the native pole. */
   if (prj->theta0 == 90.0) {
      /* Reference point is at the native pole. */

      if (dophip) {
         /* Set default for longitude of the celestial pole. */
         cel->ref[2] = 180.0;
      }

      latp = cel->ref[1];
      cel->ref[3] = latp;

      cel->euler[0] = cel->ref[0];
      cel->euler[1] = 90.0 - latp;
   } else {
      /* Reference point away from the native pole. */

      /* Set default for longitude of the celestial pole. */
      if (dophip) {
         cel->ref[2] = (cel->ref[1] < prj->theta0) ? 180.0 : 0.0;
      }

      clat0 = cosdeg (cel->ref[1]);
      slat0 = sindeg (cel->ref[1]);
      cphip = cosdeg (cel->ref[2]);
      sphip = sindeg (cel->ref[2]);
      cthe0 = cosdeg (prj->theta0);
      sthe0 = sindeg (prj->theta0);

      x = cthe0*cphip;
      y = sthe0;
      z = sqrt(x*x + y*y);
      if (z == 0.0) {
         if (slat0 != 0.0) {
            return 1;
         }

         /* latp determined by LATPOLE in this case. */
         latp = cel->ref[3];
      } else {
         if (fabs(slat0/z) > 1.0) {
            return 1;
         }

         u = atan2deg (y,x);
         v = acosdeg (slat0/z);

         latp1 = u + v;
         if (latp1 > 180.0) {
            latp1 -= 360.0;
         } else if (latp1 < -180.0) {
            latp1 += 360.0;
         }

         latp2 = u - v;
         if (latp2 > 180.0) {
            latp2 -= 360.0;
         } else if (latp2 < -180.0) {
            latp2 += 360.0;
         }

         if (fabs(cel->ref[3]-latp1) < fabs(cel->ref[3]-latp2)) {
            if (fabs(latp1) < 90.0+tol) {
               latp = latp1;
            } else {
               latp = latp2;
            }
         } else {
            if (fabs(latp2) < 90.0+tol) {
               latp = latp2;
            } else {
               latp = latp1;
            }
         }

         cel->ref[3] = latp;
      }

      cel->euler[1] = 90.0 - latp;

      z = cosdeg (latp)*clat0;
      if (fabs(z) < tol) {
         if (fabs(clat0) < tol) {
            /* Celestial pole at the reference point. */
            cel->euler[0] = cel->ref[0];
            cel->euler[1] = 90.0 - prj->theta0;
         } else if (latp > 0.0) {
            /* Celestial pole at the native north pole.*/
            cel->euler[0] = cel->ref[0] + cel->ref[2] - 180.0;
            cel->euler[1] = 0.0;
         } else if (latp < 0.0) {
            /* Celestial pole at the native south pole. */
            cel->euler[0] = cel->ref[0] - cel->ref[2];
            cel->euler[1] = 180.0;
         }
      } else {
         x = (sthe0 - sindeg (latp)*slat0)/z;
         y =  sphip*cthe0/clat0;
         if (x == 0.0 && y == 0.0) {
            return 1;
         }
         cel->euler[0] = cel->ref[0] - atan2deg (y,x);
      }

      /* Make euler[0] the same sign as ref[0]. */
      if (cel->ref[0] >= 0.0) {
         if (cel->euler[0] < 0.0) cel->euler[0] += 360.0;
      } else {
         if (cel->euler[0] > 0.0) cel->euler[0] -= 360.0;
      }
   }

   cel->euler[2] = cel->ref[2];
   cel->euler[3] = cosdeg (cel->euler[1]);
   cel->euler[4] = sindeg (cel->euler[1]);
   cel->flag = CELSET;

   /* Check for ill-conditioned parameters. */
   if (fabs(latp) > 90.0+tol) {
      return 2;
   }

   return 0;
}

/*--------------------------------------------------------------------------*/

int
celfwd(pcode, lng, lat, cel, phi, theta, prj, x, y)

const char pcode[4];
const double lng, lat;
struct celprm *cel;
double *phi, *theta;
struct prjprm *prj;
double *x, *y;

{
   int    err;

   if (cel->flag != CELSET) {
      if (celset(pcode, cel, prj)) return 1;
   }

   /* Compute native coordinates. */
   sphfwd(lng, lat, cel->euler, phi, theta);

   /* Apply forward projection. */
   if ((err = prj->prjfwd(*phi, *theta, prj, x, y))) {
      return err == 1 ? 2 : 3;
   }

   return 0;
}

/*--------------------------------------------------------------------------*/

int
celrev(pcode, x, y, prj, phi, theta, cel, lng, lat)

const char pcode[4];
const double x, y;
struct prjprm *prj;
double *phi, *theta;
struct celprm *cel;
double *lng, *lat;

{
   int    err;

   if (cel->flag != CELSET) {
      if(celset(pcode, cel, prj)) return 1;
   }

   /* Apply reverse projection. */
   if ((err = prj->prjrev(x, y, prj, phi, theta))) {
      return err == 1 ? 2 : 3;
   }

   /* Compute native coordinates. */
   sphrev(*phi, *theta, cel->euler, lng, lat);

   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 wcsmath.h and cel.h
 *
 * Dec 18 2000	Doug Mink - Include string.h for strcmp()
 *
 * Mar 20 2001	Doug Mink - Add () around err assignments in if statements
 * Sep 19 2001	Doug Mink - Add above changes to WCSLIB-2.7 cel.c
 *
 * Mar 12 2002	Doug Mink - Add changes to WCSLIB-2.8.2 cel.c
 */