Given two lists of source on field, find the dx,dy offset between lists.

  x1 - X coordinate from list 1, in pixels.
  y1 - Y coordinate from list 1, in pixels.
  x2 - X coordinate from list 2, in pixels.
  y2 - Y coordinate from list 2, in pixels.

  NX - maximum extent in X to consider (default is max([x1,x2]))
  NY - maximum extent in Y to consider (default is max([y1,y2]))
  MAXERR - maximum error allowed in initial spread test of position.
  FNDRAD - Size of the aperture used for the final offset measurement.
              DEFAULT VALUE = 12 pixels
           The default is provided based on its historical value.  In most
              cases this appears to work pretty well and should generally
              be left alone.  However, some data have been seen to get
              confused with a value that is this big.  Changing this value
              will require knowing a better value for a specific dataset.
  SCALEFAC - Scaling factor to apply on the initial crude offset.  The
               default is 1.0.  This control is used for images where the
               pixel scale is very oversampled and very small on an absolute
               astrometric basis.  One particular case where this was
               needed is in Magellan IMACS f/4 data where the image scale
               is 0.111 arcsec/pixel.  With seeing of 1 arcsec the offset
               calculation does not get a good correlation peak.  Binning
               the result makes the peak sharper and easier to find.  For
               this case a scalefactor of 0.5 or 0.3 worked quite well.

  xoff - X offset (2-1) between positions in each list.
  yoff - Y offset (2-1) between positions in each list.
  error - Code, set if something went wrong in correlating the lists.
            0 - everything appears to be good.
            1 - failure during input validation
            2 - spread in the initial x offset is too big (>maxerr)
            3 - spread in the initial y offset is too big (>maxerr)
            4 - correlation spot has negative "flux" or fwhm
            5 - Final pass on x offset excluded all points in robomean
            6 - Final pass on y offset excluded all points in robomean
            7 - All of the final pass x offsets were bigger than 1.5*xsize
    FOM - Figure of merit, a number than can be used (differentially) to
           measure how good the spatial correlation is.  This number is
           approximately the fraction of objects in the shortest list that
           ended up spatially correlated.  A number close to 1 should be
   INDEX- index into list 2 for points in list1, ie, list2[index[i]] is the
         closest, or one of a group of closest points, in list 2 to the
         ith element of list 1, given the xoff, yoff determined.
         If SPATIAL is specified, elements of the index may be invalidated
         by setting to -1- these represent invalid matches from list 1.
         On return, if error is set, the index output should be ignored.
 SPATIAL- Filtering parameters to frmdxdy - a vector or scalar of either 1 
           or 2 elements- the first is the max distance in pixels from the
           mean correlation dx,dy for a match to be valid, and the 2nd is
           the maximum threshold in sigma from the mean correlation dx,dy
           for a match to be valid. If SPATIAL is specified, invalid matches
           will be excluded (via an initial distance test and robomean) from 
           both the IDX 
           outputs and the final dx,dy result- if not specified, the distance
           and sigma criteria will be defaulted by frmdxdy.  In this case, 
           invalid matches will be excluded from the final dx,dy result but
           INCLUDED in the IDX output. Note that if only the first element of
           SPATIAL is specified, the second is defaulted to 3.0. The default
           assumed for spatial[0] is 3.0.



 It is conventional (and faster) for list 1 to be the shorter
 of the two lists.   Success is independent of the order in which lists 
 are presented, although if list 1 is longer than list 2, the index generated
 will not be unique (many->1).


  99/03/22, Written by Marc W. Buie, Lowell Observatory
  2005/06/21, MWB, changed called to robomean to trap errors.
  2007/11/21, MWB, merged with alternate code buried in
  2009/07/23, MWB, modified so that x,y input arrays do not have to be
                    positivie definite.
  2009/07/24, MWB, added XOUT,YOUT optional output.
  2010/02/14, MWB, merged with alternate version from Peter Collins, this
                     brings in the INDEX and SPATIAL keywords.
  2010/07/19, MWB, minor tweak to ensure that the error flag is set for
                     all cases of premature return.  Added FNDRAD keyword.
  2012/12/03, MWB, error code 3 never returned, fixed.
                     Added SCALEFAC keyword