; REGRD
;---------------------------------------------------------------
;! Regrids an image from one co-ordinate frame to another
;# TASK IMAGING COORDINATES
;-----------------------------------------------------------------------
;;  Copyright (C) 1995, 2002
;;  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
;-----------------------------------------------------------------------
;---------------------------------------------------------------
REGRD     LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
REGRD: Transforms image coordinate system and geometry
USERID     -32000.0     32000.0    User ID.  0 => current user,
                                     32000 => any user.
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
APARM                              Input map parameters, (if not
                                   specified in the map header)
                                   1) epoch of mean coordinates
                                   2) epoch prefix
                                     1: Julian    (eg J2000.0)
                                     2: Besselian (eg B1950.0)
                                     3: Besselian without
                                        E-terms   (eg b1950.0)
BPARM                              Output map parameters
                                   1) coordinate system
                                     1: equatorial
                                     2: galactic
                                     3: ecliptic
                                   2) epoch of mean coordinates
                                   3) epoch prefix
                                     1: Julian    (eg J2000.0)
                                     2: Besselian (eg B1950.0)
                                     3: Besselian without
                                        E-terms   (eg b1950.0)
                                   4) geometry
                                     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
                                   1-5): first axis (see HELP)
                                     1: hour   (or degree)
                                     2: minute (or arcmin)
                                     3: second (or arcsec)
                                     4: reference pixel
                                     5: coord increment (arcsec)
                                   6-10): second axis similarly

          See HELP for important information concerning
          transforming coordinates, and usage of CPARM.
----------------------------------------------------------------
REGRD
Type: Task
Use:  REGRD will regrid an image from one coordinate frame
      and geometry to another.  This includes precession as
      described below.

                        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 REGRD implies that the
      coordinates are referenced to the new IAU1976/FK5 system.

      Specifying a Besselian epoch to REGRD 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 cooordinates 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.

Adverbs:
  USERID......User ID of owner of image.  0 => current user,
                32000 => any user.
  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.
  APARM.......Input map parameters to be used if not defined in
              the image header.
              1) epoch of mean equatorial or ecliptic
                 coordinates, e.g. 1950, 2000.
              2) epoch prefix,
                  1 = Julian    (as J2000.0)
                  2 = Besselian (as B1950.0)
                  3 = Besselian without E-terms (eg b1950.0)
              If not defined in the header or APARMs, the epoch
              defaults to 2000.0.  The prefix defaults to B if
              the epoch is 1950.0, otherwise to J,
  BPARM.......Coordinate frame and geometry of the output map
              1) coordinate frame
                 1: equatorial (mean of epoch)
                 2: galactic
                 3: ecliptic (mean of epoch)
                 anything else defaults to the input map value.
              2) epoch of mean equatorial or ecliptic
                 coordinates, e.g. 1950, 2000.
                 If negative or zero, defaults to the input
                 value if defined, otherwise 2000.0
              3) epoch prefix,
                  1 = Julian    (as J2000.0)
                  2 = Besselian (as B1950.0)
                  3 = Besselian without E-terms (eg b1950.0)
                 Anything else defaults to the input value if
                 defined, otherwise if the epoch is 1950.0 it
                 defaults to B, otherwise to J.
              4) geometry
                 1: SIN, sine (orthographic)
                 2: TAN, tangent (gnomonic)
                 3: ARC, arc (zenithal equidistant)
                 4: NCP, north celestial pole tangent
                 5: STG, stereographic
                 6: AIT, Aitov
                 7: GLS, global sinusoid (Sanson-Flamsteed)
                 8: MER, Mercator
              5) output blanking control
                 0: "magic" blanking
                 1: zero
  CPARM.......Output axis specification: all zero causes the input
              reference pixel to be used (as transformed).
              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
                    (TRANSFORMED) COORDINATES OF THE CENTRE OF
                    THE INPUT MAP WILL BE USED.
                 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 map is assumed.
                 5: coordinate increment (arcsec per pixel,
                    should be negative).  If zero, the input
                    coordinate increment will be scaled
                    according to the input and output image
                    sizes.
              6-10) apply to the second axis as for the first
                    except the range is -90 to +90 degrees
----------------------------------------------------------------
REGRD:   Transforms image coordinate system and geometry
Author:  Mark Calabretta
Related tasks: 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 (equatorial FK4 coordinates only)
   b) perform the spherical coordinate rotation specified by
      three Euler angles
   c) apply E-terms (equatorial FK4 coordinates only)
3) compute the pixel coordinates on the input map
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 (steps 2a and 2c,
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.
----------------------------------------------------------------
