idlastro / FITS Astrometry and Calibration: WCSSPH2XY

[Source code]

NAME
WCSSPH2XY
PURPOSE
Convert spherical coordinates to x and y (map) angular coordinates
EXPLANATION
Convert spherical (longitude and latitude -- sky) coordinates to x
and y intermediate world coordinates (still nominally in degrees) in
the projection plane of the map.  This procedure is the inverse of
WCSXY2SPH.    See WCS_DEMO for example of use.
This is a lower level procedure -- given a FITS header, the user will
usually use ADXY which will then call WCSSPH2XY with the appropriate
parameters.
CATEGORY
Mapping and Auxiliary FITS Routine
CALLING SEQUENCE
wcssph2xy, longitude, latitude, x, y, [ map_type , CTYPE = ,
         FACE =, PV1= PV2= , CRVAL = , CRXY = , LONGPOLE = ,
         LATPOLE = , PHI0 = , NORTH_OFFSET =, SOUTH_OFFSET =, BADINDEX =]
INPUT PARAMETERS
longitude - longitude of data, scalar or vector, in degrees
latitude  - latitude of data, same number of elements as longitude,
            in degrees
map_type  - optional positional parameter, numeric scalar (0-26)
          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.
S  Number  Name                       Comments
e   code
-  ------  -----------------------    -----------------------------------
F     0    Default = Plate Carree
P     1    Zenithal perspective       PV2_1 required
N     2    Gnomic                     AZP w/ mu = 0
N     3    Orthographic               PV2_1,PV2_2 optional
G     4    Stereographic              AZP w/ mu = 1
C     5    Zenithal Equidistant
N     6    Zenithal polynomial        PV2_0, PV2_1....PV2_20 possible
A     7    Zenithal equal area
R     8    Airy                       PV2_1 required
P     9    Cylindrical perspective    PV2_1 and PV2_2 required
R    10    Plate Carree
R    11    Mercator
A    12    Cylindrical equal area     PV2_1 required
P    13    Conical perspective        PV2_1 and PV2_2 required
D    14    Conical equidistant        PV2_1 and PV2_2 required
E    15    Conical equal area         PV2_1 and PV2_2 required
O    16    Conical orthomorphic       PV2_1 and PV2_2 required
N    17    Bonne's equal area         PV2_1 required
O    18    Polyconic
L    19    Sanson-Flamsteed  (GLS is allowed as a synonym for SFL)
R    20    Parabolic
T    21    Hammer-Aitoff
L    22    Mollweide
C    23    Cobe Quadrilateralized     convergence of inverse is poor
           Spherical Cube
C    24    Quadrilateralized
           Spherical Cube
C    25    Tangential Spherical Cube
P    26    Slant Zenithal Projection   PV2_1,PV2_2, PV2_3 optional
X    27    HealPix
T    28    HealCart (Cartesian approximation of Healpix)
H    29    HEALPix butterfly projection
OPTIONAL INPUT 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.
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.
CRXY -    2 element vector giving the x and y coordinates of the
          reference point. if this is not set the offset is [0,0].
          Used to implement (x0,y0) in Sect 2.5 of Griesen & Calabretta 2002
          Do not confuse with CRPIX.
arameters simply passed to WCS_ROTATE:
CRVAL - 2 element vector containing standard system coordinates (the
          longitude and latitude) of the reference point
PV1   - Vector of projection parameters associated with longitude
LONGPOLE -  native longitude of standard system's North Pole
LATPOLE  -  "target" native latitude of the standard system's North Pole
arameters intended to enhance invertability:
NORTH_OFFSET - offset (radians) added to input points near north pole.
SOUTH_OFFSET - offset (radians) added to input points near south pole.
OUTPUT PARAMETERS
x - x coordinate of data, same number of elements as longitude, in
        degrees; if CRXY is set, then x will be returned offset by
        crxy(0).  NOTE: x in all map projections increases to the
        left, not the right.
y - y coordinate of data, same number of elements as longitude, in
        degrees; if CRXY is set, y will be returned offset by crxy[1]
OPTIONAL OUTPUT KEYWORD PARAMETERS
FACE - a output variable used for spherical cube projections to
        designate the face of the cube on which the x and y
        coordinates lie.   Will contain the same number of elements as
        X and Y.    Must contain at least 1 arbitrary element on input
        If FACE is NOT defined on input, it is assumed that the
        spherical cube projection is laid out over the whole sky
        in the "sideways T" configuration.
DINDEX - vector, list of transformed points too close to poles.
NOTES
The conventions followed here are described in more detail in
"Representations of Celestial Coordinates in FITS" by Calabretta
and  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 first use WCS_ROTATE to convert
coordinates in one of three standard systems (celestial, galactic,
or ecliptic) into a "native system" of latitude and longitude.  The
latitude and longitude are then converted into x and y coordinates
which depend on the map projection which is performed.   The rotation
from standard to native coordinates can be skipped if one so desires.
This procedure necessitates two basic sections.  The first converts
"standard" coordinates to "native" coordinates while the second converts
"native" coordinates to x and y coordinates.  The first section is
simply a call to WCS_ROTATE, while the second contains the guts of
the code in which all of the map projection is done.  This procedure
can be called in a form similar to AITOFF, EQPOLE, or QDCB by calling
wcssph2xy with a fifth parameter specifying the map projection by
number and by not using any of the keywords related to the map
projection type (e.g. CTYPE).
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 to a map projection
type).  All variables are converted into double precision values and
angular measurements are converted from degrees into radians.
If necessary, longitude values are converted into the range -pi to pi.
Any latitude points close to the  of the poles are mapped to a specific
latitude of  from the pole so that the map transformations become
completely invertible.  The magnitude of this correction is given by
the keywords NORTH_OFFSET and SOUTH_OFFSET and a list of affected
points is optionally returned in the "badindex" output parameter.
The next task of the procedure is to convert the "standard"
coordinates to "native" coordinates by rotating the coordinate system.
This rotation is performed by the procedure WCS_ROTATE and is governed
by the keywords CRVAL and LONGPOLE.   The final task of the WCSSPH2XY
is to take "native" latitude and longitude coordinates and convert
them into x and y coordinates.  Any 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).
Note that a further transformation (using the CD matrix) is required
to convert the (x,y) coordinates to pixel coordinates.
COMMON BLOCKS
none
PROCEDURES CALLED
WCS_ROTATE
ORIGINAL AUTHOR
Rick Balsano   LANL   V1.1     8/31/93
MODIFICATIONS/REVISION LEVEL
2.3     9/15/93  W. Landsman (HSTX) Update quad cube coords, vectorize
                 keywords
2.4     12/29/93 I. Freedman (HSTX) Eliminated LU decomposition
2.5     1/5/93   I. Freedman (HSTX) Offset keywords / bad point index
2.6     Dec 94   Compute pole for transformations where the reference
                pixel is at the native origin    W. Landsman (HSTX)
2.7     May 95  Change internal variable BETA for V4.0 compatibility
2.8     June 95 Change loop indices from integer to long
2.9     3/18/96 Change FACE usage for cube projections to match WCSLIB
                C/FORTRAN software library.
2.10    02/18/99 Fixed implementation of ARC algorithm
2.11    June 2003 Update conic projections, add LATPOLE keyword
 2.12    Aug 2003, N.Rich - Fix pre-V5.5 bug from previous update
2.13    Sep 2003, W. Landsman CTYPE keywords need not be 8 characters
2.14    Jan 2004, W. Landsman don't modify scalars, fix PARabolic code
2.15    Feb 2004, W. Landsman Fix AZP and AIR algorithms
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     Jul 2005 W.Landsman/C. Markwardt Set unprojectable points in
            tangent projection to NaN
3.1.1   Jul 2005 Fixed 3.1 mod to work for scalars
3.2     Dec 2005 Fixed Airy projection for latitude centered at 90 deg
3.3     Aug 2007 R. Munoz, W.Landsman Correct treatment of PV1_2 and
                 PV2_2 parameters
3.4    Oct 2007  Sergey Koposov Support HEALPIX projection
3.4.1  June 2009 Check for range of validity of ZPN polynomial W.L.
3.5    May 2012  Benjamin Alan Weaver, Add nonstandard HEALCART 
                 projection, Allow map_index to be > 25
3.5.1  May 2013  W. Landsman Allow GLS as a synonym for SFL
3.6    Jul 2013  J. P. Leahy added XPH projection, apply polar offsets
                 only for cylindrical & conic projections. 
3.6.1  Dec 2013  W. Landsman Polar offsets done in radians
3.6.2  Jan 2016  W. Landsman Lat and Long can have different size so long 
                 as they have the same number of elements