NAME:
    findmax
 PURPOSE: (one line)
    Find the interpolated local maximum in a 2-D array.
 DESCRIPTION:

 CATEGORY:
    Numerical
 CALLING SEQUENCE:
    findmax, x, y, f, xmax, ymax, fmax, $
             DELTA=in_delta, EPS=in_eps, GOLD=in_gold
 INPUTS:
    x, y   : Position of the initial guess.
    f      : The 2-D function array.
 OPTIONAL INPUT PARAMETERS:

 KEYWORD PARAMETERS:
    DELTA   = Half-width of the box containing the desired maximum.
              Default is 1.0 pixel.
    EPS     = Stop criterion. Default=1.0E-5.
    GOLD    = Pad value on DELTA.  Default is 1.0E-4.
 OUTPUTS:
    xmax, ymax  : Position of computed maximum.
    fmax        : Computed function maximum.
 COMMON BLOCKS:

 SIDE EFFECTS:

 RESTRICTIONS:

 PROCEDURE:
    Binary search (2D).
    A guess for the location of the maximum is chosen.  The external function
 sint2d is called to obtain interpolated function values at two symmetric
 points along each axis (x and y).  For each axis, the two points are used
 to determine which way to shift the location of the maximum.  IF the
 function values at the two points are not equal, the location of the
 maximum is shifted by half the previous amount in the indicated direction,
 and a new set of four points are computed at half the offset used previously.
 This process continues until the offset falls below some small threshold
 value (epsilon).
 MODIFICATION HISTORY:
    Written by Doug Loucks, Lowell Observatory, September, 1993.
    2007/06/21, MWB, promote x,y inputs to float (internally only).