Merge commit '339420dd5dd874c41f6bab5808291fb4036dd022' as 'funtools'

1 files changed, 474 insertions, 0 deletions

diff --git a/funtools/wcs/cel.c b/funtools/wcs/cel.c
new file mode 100644
index 0000000..9109284
--- /dev/null
+++ b/funtools/wcs/cel.c
@@ -0,0 +1,474 @@
+/*=============================================================================
+*
+*   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
+ */