idlastro / FITS Astrometry and Calibration: WCSXY2SPH

[Source code]

NAME
WCSXY2SPH
PURPOSE
Convert x and y (map) coordinates to spherical coordinates
EXPLANATION
To convert x and y (map) coordinates to spherical (longitude and
latitude or sky) coordinates.    This procedure is the inverse of
WCSSPH2XY.
his is a lower level procedure -- given a FITS header, the user will
sually use XYAD which will then call WCSXY2SPH with the appropriate
arameters.
CATEGORY
Mapping and Auxilliary FITS Routine
CALLING SEQUENCE
wcsxy2sph, x, y, longitude, latitude, [map_type], [ CTYPE = ,$
       FACE = , PV1 =, PV2 = ,CRVAL =, CRXY =, LONGPOLE=, LATPOLE=]
INPUT PARAMETERS
x - x coordinate of data, scalar or vector, in degrees, NOTE: x
        increases to the left, not the right
y - y coordinate of data, same number of elements as x, in degrees
map_type - optional positional parameter, scalar corresponding to a
        particular map projection.  This is not a FITS standard, it is
        simply put in to allow function similar to that of less general
        map projection procedures (eg AITOFF).  The following list gives
        the map projection types and their respective numbers.
 Number  Name                       Comments
  code
 ------  -----------------------    -----------------------------------
    0    Default = Plate Carree
    1    Zenithal perspective       pv2_1 required
    2    Gnomic                     AZP w/ pv2_1 = 0
    3    Orthographic               pv2_1, pv2_2 optional
    4    Stereographic              AZP w/ pv2_1 = 1
    5    Zenithal Equidistant
    6    Zenithal polynomial        PV2_0, PV2_1....PV2_20 possible
    7    Zenithal equal area
    8    Airy                       pv2_1 required
    9    Cylindrical perspective    pv2_1 and pv2_2 required
   10    Plate Carree
   11    Mercator
   12    Cylindrical equal area     pv2_1 required
   13    Conical perspective        pv2_1 and pv2_2 required
   14    Conical equidistant        pv2_1 and pv2_2 required
   15    Conical equal area         pv2_1 and pv2_2 required
   16    Conical orthomorphic       pv2_1 and pv2_2 required
   17    Bonne's equal area         pv2_1 required
   18    Polyconic
   19    Sanson-Flamsteed (GLS is allowed as a synonym for SFL)
   20    Parabolic
   21    Hammer-Aitoff
   22    Mollweide
   23    Cobe Quadrilateralized     inverse converges poorly
         Spherical Cube
   24    Quadrilateralized
         Spherical Cube
   25    Tangential Spherical Cube
   26    Slant Zenithal perspective  PV2_1,PV2_2, PV2_3 optional
   27    HEALPix projection          pv2_1 and pv2_2 optional
   28    HealCart (Cartesian approximation of Healpix)
   29    HEALPix butterfly projection (centred on a pole)
OPTIONAL KEYWORD PARAMETERS
CTYPE - One, two, or three element vector containing 8 character
        strings corresponding to the CTYPE1, CTYPE2, and CTYPE3
        FITS keywords:
        CTYPE[0] - first four characters specify standard system
        ('RA--','GLON' or 'ELON' for right ascension, galactic
        longitude or ecliptic longitude respectively), second four
        letters specify the type of map projection (eg '-AIT' for
        Aitoff projection)
        CTYPE[1] - first four characters specify standard system
        ('DEC-','GLAT' or 'ELAT' for declination, galactic latitude
        or ecliptic latitude respectively; these must match
        the appropriate system of ctype1), second four letters of
        ctype2 must match second four letters of ctype1.
        CTYPE[2] - if present must be the 8 character string,'CUBEFACE',
        only used for spherical cube projections to identify an axis
        as containing the face on which each x and y pair of
        coordinates lie.
FACE - a input variable used for spherical cube projections to
        designate the face of the cube on which the x and y
        coordinates lie.   Must contain the same number of elements
        as X and Y.
CRVAL - 2 element vector containing standard system coordinates (the
        longitude and latitude) of the reference point
CRXY - 2 element vector giving the x and y coordinates of the
        reference point, if this is not set the offset of the x
        coordinate is assumed to be 0. 
        See Calabretta & Griesen Sec 2.5.
PV2  - Vector of projection parameter associated with latitude axis
      PV2 will have up to 21 elements for the ZPN projection, up to 3
      for the SIN projection and no more than 2 for any other
      projection.   The first element corresponds to PV2_1, the
      second to PV2_2, etc.
ameters simply passed to WCS_ROTATE:
VAL - 2 element vector containing standard system coordinates (the
        longitude and latitude) of the reference point
1   - Vector of projection parameters associated with longitude
NGPOLE -  native longitude of standard system's North Pole
TPOLE  -  "target" native latitude of the standard system's North Pole
OUTPUT PARAMETERS
longitude - longitude of data, same number of elements as x, in degrees
latitude - latitude of data, same number of elements as x, in degrees
Longitude and latitude will be set to NaN, wherever elements of X,Y
have no corresponding longitude, latitude values.
NOTES
The conventions followed here are described in more detail in the paper
Representations of Celestial Coordinates in FITS" by Calabretta &
Greisen (2002, A&A, 395, 1077, also see
http://fits.gsfc.nasa.gov/fits_wcs.html).   The general scheme
outlined in that article is to convert x and y coordinates into a
"native" longitude and latitude and then rotate the system into one of
three generally recognized systems (celestial, galactic or ecliptic).
This procedure necessitates two basic sections.  The first converts
x and y coordinates to "native" coordinates while the second converts
"native" to "standard" coordinates.  The first section contains the
guts of the code in which all of the map projection is done.  The
second step is performed by WCS_ROTATE and only involves rotation of
coordinate systems.  WCSXY2SPH can be called in a form similar to
AITOFF, EQPOLE, or QDCB by calling wcsxy2sph with a fifth parameter
specifying the map projection by number and by not using any of the
keywords related to the map projection type (eg ctype1 and ctyp2).
PROCEDURE
The first task of the procedure is to do general error-checking to
make sure the procedure was called correctly and none of the
parameters or keywords conflict.  This is particularly important
because the procedure can be called in two ways (either using
FITS-type keywords or using a number corresponding a map projection
type).  All variables are converted into double precision values.
The second task of the procedure is to take x and y coordinates and
convert them into "native" latitude and longitude coordinates.
Map-specific error-checking is done at this time.  All of the
equations were obtained from "Representations of Celestial
Coordinates in FITS" and cases needing special attention are handled
appropriately (see the comments with individual map projections for
more information on special cases).     WCS_ROTATE is then called to
convert the "native" coordinates to "standard" coordinates by rotating
the coordinate system.  This rotation is governed by the keywords
CRVAL, and LONGPOLE.  The transformation is a straightforward
application of euler angles.  Finally, longitude values are converted
into the range from 0 to 360 degrees.
COMMON BLOCKS
none
PROCEDURES CALLED
WCS_ROTATE
ORIGINAL AUTHOR
Rick Balsano  LANL    8/31/93   V1.1
MODIFICATIONS/REVISION LEVEL
1.2 9/12/93 W. Landsman Vectorized CRXY, CRVAL, CTYPE
1.3 29/12/93 I. Freedman Eliminated LU decomposition
1.4 22/09/94 W. Landsman If scalar input, then scalar output
1.5 02/03/05 W. Landsman Change variable name BETA for V4.0 compatibility
1.6 06/07/05 W. Landsman Change loop index from integer to long
1.7 02/18/99 W. Landsman Fixed implementation of ARC algorithm
1.8 June 2003 W. Landsman Update conic projections, add LATPOLE keyword
1.81 Sep 2003 W. Landsman Avoid divide by zero
1.82 Sep 2003 W. Landsman CTYPE keywords need not be 8 characters
1.83 Sep 2003 W. Landsman Preserve input array sizes
1.9 Jan 2004 W. Landsman don't modify scalars, fix PARabolic code
2.0 Feb 2004 W. Landsman Fix AIR and AZP projections
2.1 Feb 2004 W. Landsman Fix tangent projection for matrix input
3.0 May 2004 W. Landsman Support extended SIN (=NCP), slant zenithal
(SZP), and zenithal polynomial (ZPN) projections, use
 PV2 keyword vector instead of PROJP1, PROJP2
3.1 May 2004 W. Landsman/J. Ballet Handle NaN values, flag invalid output
for AITOFF projection
3.1.1 Dec 2005 W. Landsman/W. Thompson Fixed problem with Airy projection
centered on 90 degree latitude
3.1.2 May 2006 W. Landsman/Y.Sato Fix problem selecting the correct root
for the ZPN projection
3.2 Aug 2007 W. Landsman Correct treatment of PVi_j parameters
3.3 Oct 2007 Sergey Koposov Support HEALPIX projection
3.4 May 2012 Benjamin Alan Weaver, Add nonstandard HEALCART
projection, Allow map_index to be > 25
3.4.1 May 2013 W. Landsman Allow GLS as a synonym for SFL
3.5 Jul 2013 J. P. Leahy Add nonstandard XPH projection and
improved HEALPix support; changed sign of CRXY
for consistency with WCSSPH2XY; introduced PV1
3.5.1 Dec 2013 W. Landsman Return scalar for scalar input for ZPN proj.