PRO annstats, image, x0, y0, radius, tot, npix, mean, sigma, stderr, $
    mask=mask, verbose=verbose      
;+
; NAME:
;    ANNSTATS
; PURPOSE:
;    Integrate flux in annular apertures and get std. dev. in each
;    annulus.  User responsible for subtracting sky first.  See 
;    notes (1) and (2) below.
; CALLING SEQUENCE:
;    annstats, image, x0, y0, radius, tot, npix, sigma, stderr, $
;           mask=mask, /verbose
; INPUT PARAMETERS:
;    image      2-D image array
;    x0, y0     center of circular annuli (need not be integral)
;    radius     vector of annulus radii
; OUTPUT PARAMETERS:
;    For each annulus,
;       tot     =   total of pixel values
;       npix    =   number of pixels 
;       mean    =   mean value per pixel
;       sigma   =   standard deviation (scatter) among pixels
;       stderr  =   standard error of the mean
; OPTIONAL INPUT KEYWORDS:
;    mask       Must be same size as image.  If a pixel in
;               the mask is 1, the corresponding pixel in the image is
;               counted.  If 0, the corresponding pixel in the image is
;               ignored (in all results).  See note (3).
;    verbose    Use of /verbose will make the routine tell you what stage
;               it's at.  See note (1).
; COMMON BLOCKS:
;    none
; NOTES:
;   (1) If you aren't sure you've set up right, use the /verbose keyword,
;       because the routine is fairly slow.
;   (2) Mask is intended to mask out stars or garbage.  Depending on your
;       application, you might be better off modifying the input image.
; SIDE EFFECTS:  none
; PROCEDURE:   Similar to IDL/UIT/DAOPHOT aperture routine.
; MODIFICATION HISTORY:
;  Modified from CIRCINT.  R.S.Hill, Hughes STX Corp., 5 August 1993.
;  Message fixed.  One aperture possible.  RSH, HSTX, 12-Oct-1993.
;-
on_error,1
IF n_params(0) LT 1 THEN BEGIN
   message,'Calling sequence:  annstats,image,x0,y0,' $
          +'radius,tot,npix,mean,sigma,stderr,mask=mask,/verbose'
ENDIF
IF NOT keyword_set(verbose) THEN verbose=0
sz       = size(image)
bell     = string(07b)
IF keyword_set(mask) THEN BEGIN
   szm = size(mask)
   IF (szm(1) NE sz(1)) OR (szm(2) NE sz(2)) THEN BEGIN
      message,'Mask image must be same size as flux image.'
   ENDIF
ENDIF
nrad       = n_elements(radius)                                             
IF nrad GT 1 THEN BEGIN
   order_test = radius(1:nrad-1) - radius(0:nrad-2)
   ww         = where((order_test LE 0),count)
   IF count GT 0 THEN $                             
      message,'Radii must be specified in ascending order.' 
ENDIF
ww  = where((radius LT 0.5), count)
IF count GT 0 THEN $
   message,'All radii must be 0.5 pixels or greater' $
      +' (first element of returned arrays is for central disc).'
IF verbose NE 0 THEN $ 
   message,/inf,'Beginning routine to integrate in circular annuli ...'
maxrad   = float(max(radius))
maxrad2  = maxrad^2
distance = replicate(2.0*maxrad2, sz(1), sz(2)) ; Like FLTARR with init
xdist2   = (findgen(sz(1)) - x0)^2
ydist2   = (findgen(sz(2)) - y0)^2
xmin = min(where(xdist2 LE maxrad2)) - 1
xmax = max(where(xdist2 LE maxrad2)) + 1
ymin = min(where(ydist2 LE maxrad2)) - 1
ymax = max(where(ydist2 LE maxrad2)) + 1
FOR ix=xmin,xmax DO FOR iy=ymin,ymax DO $ 
   distance(ix,iy) = xdist2(ix) + ydist2(iy)
distance = sqrt(distance)-0.5
tot      = fltarr(nrad)
npix     = tot
mean     = tot
sigma    = tot
stderr   = tot
;
IF verbose NE 0 THEN message,/inf,'Initialization complete.  Now integrating.'
FOR ir=0, nrad-1 DO BEGIN
   outer = radius(ir)
   IF ir GT 0 THEN inner = radius(ir-1) ELSE inner = -0.5
   IF verbose NE 0 THEN BEGIN
      message,/inf,'Now doing annulus of radii ' $
         +strtrim(inner,2)+' and '$ 
         +strtrim(outer,2)+' pixels'
      IF inner LT 0 THEN message,/inf,'     (central disc)'
   ENDIF
   IF keyword_set(mask) THEN BEGIN
      within = (distance LE outer) $
      AND (distance GE (inner)-1) AND mask
   ENDIF ELSE within = distance LE outer $
      AND (distance GE inner-1)
   thisapr = distance(where(within))
   fractno = ((outer-thisapr) < 1.0 >0.0 ) ;Fraction of pixels to count
   fractni = ((thisapr+1-inner) < 1.0 >0.0 ) ;Forced to 1 for central disc
   fractn  = fractno * fractni
   subimage= image(where(within))
   tot(ir)    = total(subimage*fractn)
   npix(ir)   = total(fractn)
   mean(ir)   = tot(ir)/npix(ir)
   deviation  = subimage - mean(ir)
   variance   = deviation^2
   weighted_variance = fractn*variance
   average_variance  = total(weighted_variance)/(npix(ir)-1.0)
   sigma(ir) = sqrt(average_variance)
   stderr(ir) = sigma(ir)/sqrt(npix(ir))
   subimage=0 & thisapr=0 & fractno=0 & fractni=0 & fractn=0 & within=0
   deviation=0 & variance=0 & weighted_variance=0
ENDFOR
IF verbose NE 0 THEN message,/inf,bell+'Done.'
RETURN
END