/astrom                                                   April 2008

This directory contain IDL procedures for manipulating FITS images which include
header keywords specifying the coordinate system.   The FITS header must follow
the world coordinate (WCS) specification in the paper "Representations of
Celestial Coordinates in FITS"  by Eric Greisen and Mark Calabretta (2002, A&A,
395, 1061).   This paper and other WCS documentation may be obtained from    
http://fits.gsfc.nasa.gov/fits_wcs.html.  All projections in that paper are 
implemented as of May 2004      The IDL code has been  validated using the test
images on http://www.atnf.csiro.au/people/mcalabre/WCS/sample_data/index.html 

Bill Thompson has written a different IDL implementation of the WCS
specification which also includes support for solar and spectroscopic
coordinates.    His routines are available at the SolarSoftWare (SSW) library
at     http://sohowww.nascom.nasa.gov/solarsoft/gen/idl/wcs/

The FITS convention defines the center of the first pixel in the
array to have coordinates (1,1).    This differs from the IDL convention where
the first pixel has coordinates (0,0).    Whenever an (X,Y) position is
extracted or inserted into a FITS header then this difference must be accounted
for.

Several procedures use an IDL structure that contains astrometry information in
the following tags:

  .NAXIS - 2 element long vector giving dimensions of the images
  .CD	-  2 x 2 array containing the astrometry parameters CD1_1 CD1_2
	   in DEGREES/PIXEL				      CD2_1 CD2_2
  .CDELT - 2 element vector giving physical increment at the reference pixel
  .CRPIX - 2 element vector giving X and Y coordinates of reference pixel
	   (def = NAXIS/2) in FITS convention (first pixel is 1,1)
  .CRVAL - 2 element double precision vector giving R.A. and DEC of 
	 reference pixel in DEGREES
  .CTYPE - 2 element string vector giving projection types, default
	 ['RA---TAN','DEC--TAN']
  .LONGPOLE - scalar longitude of north pole (default = 180)
  .LATPOLE - scalar giving native latitude of the celestial pole default=0)
  .PV2 - Vector of parameters (PV2_1, PV2_2...) needed in some projections 
  .DISTORT - Optional substructure giving distortion parameters.   Currently
         only implemented for the Spitzer simple imaging polynomial (SIP) see
         http://ssc.spitzer.caltech.edu/postbcd/doc/shupeADASS.pdf

The procedures EXTAST and PUTAST can be used, respectively, to extract
astrometry information from a FITS header and place it into an astrometry
structure, and to put the information from an astrometry structure into a FITS
header.  

      IDL> extast, hdr, astr     ;Extract astrometry from a FITS header
      IDL> putast, hdr, astr     ;Put astrometry info into a FITS header 

(1) Following the discussion in Section 5 of the Greisen & Calabretta paper,
this software supports three WCS formats: (1) the PC matrix with the separate
CDELT parameters (preferred), (2) the CD matrix without CDELT parameters (IRAF
standard),  and (3) the old-style AIPS (CDELT + CROTA) notation. The CDELT tag
is not used in format (2), and is normally set to [1.0, 1.0].    However, for
an AIPS-style header, the CDELT values are stored in the CDELT tag, and the
CROTA value is stored in the CD matrix     

cd = [ [ cos(crota), -sin(crota) ] , [ sin(crota), cos(crota)] ] 

If either element of the CDELT tag differs from 1.0 then this signals that the
astrometry came from an AIPS type header.    Programs using the astrometry
structure should always use the following code 

 if N_elements(CDELT) GE 2 then if (cdelt[0] NE 1.0) then begin
        cd[0,0] = cd[0,0]*cdelt[0] & cd[0,1] =  cd[0,1]*cdelt[0]
        cd[1,1] = cd[1,1]*cdelt[1] & cd[1,0] =  cd[1,0]*cdelt[1]
 endif


(2) Earlier draft versions of the World Coordinate System paper represented the
CD matrix in a FITS header with names such as 'CD001001' rather than 'CD1_1', and
'PC001001' instead of 'PC1_1'.   This notation is no longer recognized, but the
procedure FITS_CD_FIX can be used to convert the FITS header to modern notation.

(3) The CRPIX value is stored in the astrometry structure the same way as it is
in the FITS header, i.e. in FORTRAN (first pixel is 1) convention.    

(4) The procedures HREBIN and HCONGRID will update the astrometry, including the
SIP distortion coefficients,  when expanding or contracting an image.   However,
note that HROT or HROTATE currently will not update any SIP coefficients if
present. 

Special astrometry procedures are required for the Digitized Sky Survey Images
which were released on CD-Rom in early 1994.    Information about the Guidestar
images can be obtained from http://archive.stsci.edu/dss/
and the Guidestar images can also be obtained from the SKYVIEW facility
at http://skyview.gsfc.nasa.gov/skyview.html.  The Schmidt plates used in
this survey have a highly nonlinear plate solution that is not covered
by the Greisen & Calabretta conventions.     The procedures GSSSEXTAST, 
GSSSXYAD and GSSSADXY are the Guidestar survey analogues of the procedures
EXTAST, XY2AD, and ADXY for standard astrometry.     All the astrometry
procedures in the library will test for a Guidestar image header (by looking
for the 'PPO1' keyword) and call the appropriate procedures.    The procedure
GSSS_STDAST will convert the astrometry in a guidestar image header to 
an standard tangent projection with very slightly degraded accuracy.   

A couple of procedures in other directories also use the FITS world coordinate
system including
  IMCONTOUR - (in /astro) Contour plots with astronomical labeling (either 
               RA,Dec or arc distance from the image center
  IMDBASE - (in /database) Find all catalog sources within the field of an 
            image (given a FITS header with astrometry)

ADD_DISTORT - Add a distortion structure into a FITS header
AD2XY - Use astrometry structure to convert celestial to pixel coordinates 
ADXY  - Use FITS header to convert celestial (RA,Dec) to pixel coordinates
CONS_DEC() - Obtain the X and Y coordinates of a line of constant 
        declination
CONS_RA() - Obtain the X and Y coordinates of a line of constant right
        ascension
EXTAST- EXTract ASTrometry parameters from a FITS header into an IDL structure
FITS_CD_FIX - Update obsolete representations of the CD matrix in a FITS header 
GET_EQUINOX() - Return a numeric equinox value from a FITS header
GETROT - GET ROTation and plate scale from a FITS header
GSSS_STDAST - Insert the closest tangent projection astrometry into an STScI
         Guidestar Survey Image
GSSSADXY - Convert RA, Dec to pixel coordinates for an STScI survey image
GSSSEXTAST - Extract astrometry parameters from an STScI Survey Image
GSSSXYAD - Convert pixel coordinates to RA, Dec for an STScI survey image
HASTROM - Rotate, Congrid, and/or shift an image until astrometry matches
        that in a reference FITS header.  Used to align images.
HBOXAVE - Boxaverage an image and update astrometry in a FITS header
HCONGRID - CONGRID an image and update astrometry in a FITS header
HEULER - Convert between Galactic, celestial and ecliptic coordinates in FITS a
header
HEXTRACT - Extract a subimage and update astrometry in a FITS header
HPRECESS - Precess the astrometry in a FITS header to a new equinox.
HREBIN - REBIN an image and update the astrometry in a FITS header
HREVERSE - Reverse an image about either dimension and update astrometry
        in a FITS header                                    
HROT  - Rotate an image and update astrometry in a FITS header.
HROTATE - Apply IDL ROTATE function and update astrometry in a FITS header
MAKE_ASTR - Build an astrometry structure from input parameter values
PRECESS_CD - Precess coordinate description (CD) matrix in a FITS header 
          to a new equinox.   Called by HPRECESS
PUTAST - Put astrometry parameters (e.g. rotation, plate scale) into a FITS header. 
STARAST - Obtain an exact astrometry solution given the coordinates and
        plate position of 2 or 3 stars.
UPDATE_DISTORT - Update SIP distortion coefficients for a linear transformation	
WCS_GETPOLE -  Compute the coordinates of the native pole for non-polar projection
WCSSPH2XY - Convert between longitude,latitude to X,Y angular coordinates for
        25 different map projection types
WCSXY2SPH - Inverse of WCSSPH2XY
WCS_DEMO - Demo program for WCSSPH2XY and WCSXY2SPH
WCS_ROTATE - Rotate between standard (e.g. celestial) and native coordinates
WFPC2_METRIC -  Compute the distortion in a WFPC2 image and return coordinates
XYAD  - Use FITS header to convert pixel (X,Y) to celestial(RA, Dec) coordinates
XY2AD - Use astrometry structure to convert pixel to celestial coordinates
XYXY -  Convert X,Y values on one image to X,Y values in another image
         using astrometry in the FITS headers
