PURPOSE:   (one line only)
  Fit a numerical PSF to one or more sources in an image.
  The PSF must be provided as an input to this routine.  Suitable
    PSFs can be generated with PSFSTACK.  For this routine to work, the
    resulting PSF must be amenable to shifting by interpolation.  Thus
    an under-sampled PSF will not work with this tool.
  Generally speaking, if you ask to fit more than one source at a time,
    the sources should have overlapping PSFs.  Otherwise, you are better
    off with separate calls.  This is NOT enforced.
  CCD data processing
  image  - Input image to be fitted.
  psf    - Numerical PSF to fit to image.  This must be normalized to have
              unit volume.
  region - Fitting region for chi-sq calculation.  Provide
            [x1,x2,y1,y2] in image index coordinates.
  xin    - Starting x location for source(s)
  yin    - Starting y location for source(s)
  SILENT - Flag, if set will suppress all diagnostic output
  PSFMAX - Two element vector containing the x,y position of the peak of
             the input PSF.  The position is in the native array coordinates
             of the psf array.  If not provided, this will be computed.
             It will save time for many fits with the same psf to provide
             this value.  If computed internally, the values will be returned
             to this keyword.
  MEANSKY - Background sky value (of entire image).  If this is provided
             and calcsky is not set then this value is used for sky.  If
             calcsky is set then sky is computed and returned in this variable.
  CALCSKY - Flag, if set forces this routine to compute sky background.
             By using this flag you get a sky value that is global for the
             entire image based on up to 60000 randomly selected pixels.
             This flag is ignored if you use SKY1 and SKY2.
  GAIN   - Gain of the image, e-/ADU.  This is used to generate an array
             of uncertainties for each pixel based on photon statistics.
  RDNOISE - Readout noise of CCD (default=10)  [e-]
  NMAX    - Maximum number of iterations for the amoeba call.  Default=5000
             If this limit is reached, the source will be tagged as
             not fittable, ie., chisq = -1.
  EXPTIME - Exposure time, in seconds, default=1
  SKY1  - Inner radius of the sky annulus (unit=pixels), default=0
  SKY2  - Outer radius of the sky annulus, default=0

  PSFFACTOR- Area of constraint factor relative to the size of the PSF.
                Pixels within the distance FWHMPSF*PSFFACTOR of the starting
                location of the source(s) are used to guide the fit.
                The default value is 1.0 and seems to work pretty well.
                This value cannot be arbitrarly large or small.  In general
                it needs to be between 1.0 and 2.0.  If it is too small it
                may give bad answers.  If it is too big the program will
                likely crash.

  CSTART - optional vector, length must match xin,yin.  This gives starting
              values for the source flux that override the automatic values
              computed in this routine.  Any value that is greater than 0 is
              taken to be the starting value to use.  For those that are
              not specified (ie., <=0), the ones that are known are computed
              and subtracted prior to computing the automatic starting

  xout   - Fitted x location for source(s)
  yout   - Fitted y location for source(s)
  counts - Object flux (photons per second)
  CHISQ  - final goodness of fit for the returned values
              If this value is negative, the return values are meaningless.
  FLERR  - Uncertainty (1 sigma) of object flux.
  FLUX   - Object flux (photons per second)
  MAG    - Optional return of the instrumental magnitude
  ERR    - Optional return of the magnitude error
  FWHM   - FWHM from basphote prior to fitting
  PSF_FWHM - FWHM of PSF from basphote, same parameters as FWHM for calculation
  PSFFIT_COM, for internal use only
  Written by Marc W. Buie, Southwest Research Institute, 2012/09/28
  2012/10/04, MWB, added weighted fitting.
  2012/12/06, MWB, fixed a problem with the MEANSKY and CALCSKY keywords
  2013/03/05, MWB, added NMAX keyword and proper trap of bad fit from amoeba
  2013/03/29, BLE, Added new input keywords EXPTIME, SKY1, and SKY2.
                   Calculate and output new keywords FLUX, FLERR, MAG, ERR.
                   The PSFTHRES input keyword was also added.
  2013/04/14, MWB, tweak sigma array to avoid negative flux values
  2013/06/19, MWB, removed PSFTHRES keyword.  Added PSFFACTOR plus numerous
                     improvements including better starting values.
                     CSTART keyword was added.
  2013/07/12, MWB, added FWHM and PSF_FWHM output keywords
  2014/01/16, MWB, if number of inputs is 1 then outputs are scalars.