; SKYVE
;---------------------------------------------------------------
;! Regrids a DSS image from one co-ordinate frame to another
;# TASK IMAGING COORDINATES
;-----------------------------------------------------------------------
;;  Copyright (C) 1995, 2009
;;  Associated Universities, Inc. Washington DC, USA.
;;
;;  This program is free software; you can redistribute it and/or
;;  modify it under the terms of the GNU General Public License as
;;  published by the Free Software Foundation; either version 2 of
;;  the License, or (at your option) any later version.
;;
;;  This program 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 General Public License for more details.
;;
;;  You should have received a copy of the GNU General Public
;;  License along with this program; if not, write to the Free
;;  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge,
;;  MA 02139, USA.
;;
;;  Correspondence concerning AIPS should be addressed as follows:
;;         Internet email: aipsmail@nrao.edu.
;;         Postal address: AIPS Project Office
;;                         National Radio Astronomy Observatory
;;                         520 Edgemont Road
;;                         Charlottesville, VA 22903-2475 USA
;-----------------------------------------------------------------------
SKYVE     LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SKYVE: Transforms DSS image coordinate system and projection
INNAME                             Input name
INCLASS                            Input class
INSEQ           0.0      9999.0    Input sequence, 0 -> high
INDISK          0.0         9.0    Input disk, 0 -> any
OUTNAME                            Output name
OUTCLASS                           Output class
OUTSEQ         -1.0      9999.0    Output sequence
                                     0 -> highest unique
OUTDISK         0.0         9.0    Output disk
                                     0 -> highest with room
IMSIZE          0.0     16384.0    Output image size
BPARM                              Output map parameters
                                   1) coordinate system
                                     1: equatorial (default)
                                     2: galactic
                                     3: ecliptic
                                   2) epoch of mean coordinates
                                      (default: 2000.0)
                                   3) epoch prefix (default: J)
                                     1: Julian    (eg J2000.0)
                                     2: Besselian (eg B1950.0)
                                     3: Besselian without
                                        E-terms   (eg b1950.0)
                                   4) projection (default: NCP)
                                     1: SIN     5: STG
                                     2: TAN     6: AIT
                                     3: ARC     7: GLS
                                     4: NCP     8: MER
                                   5) blanking control
                                     0: "magic blanking"
                                     1: zeros
CPARM                              Output axis specification

                                   See HELP for important
                                   information concerning the
                                   usage of CPARM.

                                   1-5): first axis
                                     1: hour   (or degree)
                                     2: minute (or arcmin)
                                     3: second (or arcsec)
                                     4: reference pixel
                                     5: coord increment (arcsec)
                                   6-10): second axis similarly
----------------------------------------------------------------
SKYVE
Type: Task
Use:  SKYVE will regrid a Digitized Sky Survey (DSS) image to a
      coordinate frame and projection recognized by AIPS.

      The DSS is based on photographic material obtained using
      the UK Schmidt Telescope operated by the Royal Observatory
      Edinburgh (RGO), with funding from the UK Science and
      Engineering Research Council (SERC) until 1988 June, and
      thereafter by the Anglo-Australian Observatory (AAO).  The
      DSS was produced by the Space Telescope Science Institute
      (STScI) under U.S. Government grant NAG W-2166.

                            SKYVE USAGE

      Digitized Sky Survey images may be extracted from the CD
      set as FITS files by a program called 'getimage'. This is
      supplied with the CD set (on CD #61).  The FITS file
      created by 'getimage' may be read into AIPS using IMLOD.

      The DSS image coordinate system is encoded as a set of
      plate solution coefficients in FITS header cards which
      are not generally recognized by AIPS.  IMLOD stores these
      unrecognized header cards in the history file associated
      with the AIPS image.  These header cards are listed in
      the EXPLAIN file.

      SKYVE retrieves the plate solution parameters from the
      history file and regrids the image into a coordinate
      system recognized by AIPS.

      Why would you want to do this?

        a) So you can use AIPS to measure positions from the
           optical image.  Verification tests done on a
           selection of about 70 quasars with VLBI positions
           by Martin Anderson (ATNF/UWS) using MAXFIT show that
           an rms accuracy of better than 0.5 arcsec is
           achievable.

        b) So that you can overlay radio images on top of the
           optical image.

      If you want to overlay images then be sure to extract an
      optical image slightly larger than the radio image to
      avoid edge effects when the optical image is regridded.

      You must also exercise caution in changing the pixel
      spacing in the optical image to be much greater than that
      of the original DSS image (approximately 1.7 arcsec).  If
      you make the spacing too great then stars and other small
      scale objects may be skipped over.  If you already have
      the radio image then:

        a) If the pixel spacing in the radio image is much
           less than 1.7 arcsec then you should regrid the
           optical image to the same spacing as the radio image.
           You should be able to do this directly with SKYVE
           using the IMSIZE, BPARM, and CPARM adverbs to match
           the coordinate system and projection of the radio
           image.

        b) If the pixel spacing in the radio image is much
           greater than 1.7 arcsec then you should regrid the
           radio image to the same spacing as the optical image.
           This may be done with REGRD or HGEOM.

      If you don't already have the radio image then you should
      synthesize it with a cell spacing to match the optical
      image.

                        COORDINATE SYSTEMS

      Coordinate transformations between the IAU1976 and
      Bessel-Newcomb systems are done with full precision
      assuming zero proper motion, parallax, and recessional
      velocity at J2000.0

      Specifying a Julian epoch 'J' to SKYVE implies that the
      output coordinates are referenced to the new IAU1976/FK5
      system.

      Specifying a Besselian epoch to SKYVE implies that the
      coordinates are referenced to the old Bessel-Newcomb/FK4
      system.  An epoch prefix of 'B' indicates the convention
      that the coordinates include the effect of the E-terms,
      whereas 'b' indicates that they have already been removed.

      FK4 catalogue coordinates were not corrected for the
      elliptic terms of aberration (E-terms) except for
      positions within 10 degrees of the pole.  Most earlier
      catalogues did not correct for them.

      The default behaviour here is to assume that the E-terms
      are included in all Besselian coordinates (including near
      the pole).  This can be defeated if it is known that the
      input coordinates have already been corrected, or if it is
      required that the output coordinates not contain them.

      See the EXPLAIN section for a brief description of the
      algorithm.

                           STATISTICS

      SKYVE reports statistics of the pixel displacements
      (output map pixel coordinate minus input map pixel
      coordinate) for the regridding operation.

      The mean and rms for pixel displacements in X and Y are
      reported, and also the correlation coefficient.  If all
      SKYVE defaults are adopted - same image size, J2000.0
      equatorial coordinates, same centre coordinates and pixel
      increment - then the mean shift should be approximately
      zero, the rms should be a few pixels, and the correlation
      coefficient much less than unity.

      However, these statistics do not directly account for a
      net rotation of the image, and this is usually the main
      systematic difference between the output and input maps.

Adverbs:
  INNAME......Input image name,  standard defaults.
  INCLASS.....Input image class, standard defaults.
  INSEQ.......Input image sequence number, 0 -> highest.
  INDISK......Input disk drive number, 0 -> any.
  OUTNAME.....Output image name,  standard defaults.
  OUTCLASS....Output image class, standard defaults.
  OUTSEQ......Output image sequence number, 0 -> highest unique
  OUTDISK.....Output disk drive number, 0 -> highest with
              space.
  IMSIZE......Output image size (pixels), maximum 16384.
              Defaults to the input image size if negative or
              zero.
  BPARM.......Coordinate frame and projection of the output map

              1) Coordinate frame
                 0: -> 1
                 1: equatorial (mean of epoch)
                 2: galactic
                 3: ecliptic (mean of epoch)
                 Anything else produces an error.

              2) Epoch of mean equatorial or ecliptic
                 coordinates, e.g. 1950, 2000.
                 Defaults to 2000.0 if negative or zero.

              3) Epoch prefix
                  1: "J" - Julian    (as J2000.0)
                  2: "B" - Besselian (as B1950.0)
                  3: "b" - Besselian without E-terms (eg b1950.0)
                 Anything else defaults to "J".

              4) Spherical projection (geometry)
                 0: -> 4
                 1: SIN, sine (orthographic)
                 2: TAN, tangent (gnomonic)
                 3: ARC, arc (zenithal equidistant)
                 4: NCP, north celestial pole tangent
                 5: STG, stereographic
                 6: AIT, Hammer-Aitoff
                 7: GLS, global sinusoid (Sanson-Flamsteed)
                 8: MER, Mercator
                 Anything else produces an error.

              5) output blanking control
                 0: "magic" blanking
                 1: zero
                 (less than or equal to 0.5 -> 0;
                  greather than 0.5 -> 1)

  CPARM.......Output axis specification
              1-5)  Apply to the first axis

              1-3)  Specify the coordinate reference pixel.

                    IF THE VALUE SPECIFIED IS OUTSIDE THE RANGE
                    -24HR TO +24HR, OR -360 TO +360 DEGREES, THE
                    COORDINATES OF THE CENTRE OF THE INPUT MAP
                    WILL BE USED (transformed to the coordinate
                    system of the output map if necessary).

                 1: hour   for equatorial, degree for the others
                 2: minute for equatorial, arcmin for the others
                 3: second for equatorial, arcsec for the others

                 4: Coordinate reference pixel.  If zero, the
                    centre of the output map is assumed.

                 5: Coordinate increment (arcsec per pixel,
                    should be negative).  If zero, the DSS plate
                    scale is assumed.

              6-10) Apply to the second axis as for the first
                    except the range is -90 to +90 degrees, and
                    the coordinate increment (if non-zero)
                    should normally be positive.

                    THE ABSOLUTE VALUES OF CPARM(6:8) ARE USED
                    TO COMPUTE THE DECLINATION (LATITUDE); IF
                    ANY OR ALL OF CPARM(6:8) ARE NEGATIVE THE
                    DECLINATION (LATITUDE) IS NEGATED.
----------------------------------------------------------------
SKYVE:   Transforms DSS image coordinate system and projection
Author and documenter:  Mark Calabretta, ATNF
Related tasks: REGRD, GEOM, HGEOM, LGEOM, PGEOM, COMB

Algorithm
~~~~~~~~~
   For each pixel in the output image:
1) Compute its sky coordinates.
2) Transform to sky coordinates on the input map -
   a) remove E-terms (only when the output map is in equatorial
      FK4 coordinates).
   b) perform the spherical coordinate rotation specified by
      three Euler angles.
3) Compute the pixel coordinates on the input map.  This
   requires iterative inversion of the DSS plate solution
   equations.
4) Interpolate the pixel value from the nearest nine pixels -
   a) quadratic interpolation in X for each of the three rows.
   b) quadratic interpolation in Y of the result.

   Parameters used for the transformation (step 2) are recorded
in the history file.

E-terms
~~~~~~~
   The E-terms are recomputed for every pixel (step 2a,
computation of the E-terms is not done iteratively).
1) In RA:  (E1*COS(RA) + E2*SIN(RA))/COS(DEC)
2) In DEC: (E2*COS(RA) - E1*SIN(RA))*SIN(DEC) + E3*COS(DEC)

Euler angles
~~~~~~~~~~~~
   For a spherical coordinate rotation from SYSTEM1 to SYSTEM2:
1) PHI0: Longitude of the ascending node in SYSTEM1.
   Of the two points of intersection of the equators of SYSTEM1
   and SYSTEM2, the ascending node is the one where the equator
   of SYSTEM2 crosses from south to north as viewed in SYSTEM1.
2) THETA: The angle between the poles of the two systems.
   Positive for a positive rotation about the ascending node.
3) PHI: Longitude of the ascending node in SYSTEM2.

Blanking
~~~~~~~~
   Blank pixels are fully accounted for in the sense that one
blank pixel in the input map produces only one blank pixel in
the output map.  The basic criteria is that the output pixel
will be blank if and only if the pixel (P0) on the input map
nearest the position computed at step 3 above is blank.

   If P0 is not blank and any of the eight pixels surrounding
it are, then the quadratic interpolation reduces to a linear or,
if necessary, a constant interpolation or extrapolation.  In the
worst case where all of the neighbouring pixels are blank, the
"interpolated" value would be the value at P0.

DSS FITS header cards
~~~~~~~~~~~~~~~~~~~~~
     CNPIX1,           The DSS pixel coordinates of the bottom
 and CNPIX2            left-hand corner of the bottom left-hand
                       pixel of the image extracted by
                       'getimage'.
     PLTSCALE          Approximate plate scale, arcsec/mm.
     XPIXELSZ,         Plate pixel size in X and Y, micron.
 and YPIXELSZ
     PLTRAH,           J2000.0 right ascension of the plate
 and PLTRAM,           centre (hours, minutes, and seconds).
 and PLTRAS
     PLTDECSN,         J2000.0 declination of the plate centre
 and PLTDECD,          (sign, degrees, arcmin, and arcsec).
 and PLTDECM,
 and PLTDECS
     PPO3,             Plate centre offsets, micron.
 and PP06
     AMDX1, ...        Plate solution coefficients for the xi
  to AMDX13            standard plate coordinate.
     AMDY1, ...        Plate solution coefficients for the eta
  to AMDY13            standard plate coordinate.
----------------------------------------------------------------
