NAME:
  psfstack
 PURPOSE:
  Generate an average numerical psf by stacking observed images.
 DESCRIPTION:

 CATEGORY:
  CCD data processing

 CALLING SEQUENCE:
  psfstack,image,x,y,psf

 INPUTS:
  image - Input image from which psf is collected.
  x     - Scalar or vector list of positions of stars (x coordinate).
  y     - Scalar or vector list of positions of stars (y coordinate).

 OPTIONAL INPUT PARAMETERS:

 KEYWORD INPUT PARAMETERS:
  APODIZE - Flag, if set turns on an apodization constraint that forces the
               PSF decrease to zero by the edge.  This is implemented by 
               looking at the mean radial profile of the final PSF.  In
               this profile, the routine looks for a local minimum.  From
               that point to the edge (corners, really) the PSF is reduced
               linearly from that radial distance to the corners so that
               the edge goes to zero.  It is best to not use this but if
               you have problems with the PSF turning up at the edge, this can
               help make it better.
  DW      - half-width of box centered around PSF, final size is 2*DW+1,
             default size is DW=9 pixels.  It is really important that
             DW be larger than FWHM.  Usually, DW ~ 2*FWHM
  SILENT  - Flag, if set suppresses all printed output.
  SKY1    - Inner radius of the sky annulus (in pixels), default=10.
  SKY2    - Outer radius of the sky annulus (in pixels), default=40.  The sky
               value is computed for every source to get a local sky value.
  SNRTHRESH - Signal-to-noise ratio threshold.  Any source pointed to by x,y
                 that has a lower SNR than this will be ignored.  The default
                 is 25.
  TOPCLIP - Fraction of the brightest pixels that are excluded before doing
               the initial median on the psf stack.  This should be a number
               between 0 and 1, though if you set it to 1 everything would be
               excluded.  The default is 0.

 OUTPUTS:
  psf   - Final average psf, normalized to unit _volume_

 KEYWORD OUTPUT PARAMETERS:
  FLUX  - Fraction of PSF volume interior to FWHM from center of PSF.
  FWHM  - FWHM of final stacked PSF
  XCEN  - Centroided X location of center of PSF within the psf array
  YCEN  - Centroided Y location of center of PSF within the psf array

 COMMON BLOCKS:

 SIDE EFFECTS:

 RESTRICTIONS:

 PROCEDURE:

 MODifICATION HISTORY:
  97/10/22, Written by Marc W. Buie, Lowell Observatory
  2013/03/30, MWB, added APODIZE keyword
  2013/03/31, MWB, added SKY1, SKY2, TOPCLIP, FWHM, XCEN, YCEN, FLUX
                      also some tweaks to the algorithm that leads to
                      a significantly improved PSF.
  2014/01/17, MWB, changed the procedure that leads to the FWHM calculation.
                      Was using an aperture of radius=DW, now uses DW/2.
  2019/12/20, MWB, incorporated new routine, mkxyarr