;+
; NAME:
;    sumannerr
; PURPOSE: (one line)
;    Integrate over an annulus.
; DESCRIPTION:
;
;    This procedure computes the first two image moments of the pixels
;    contained within the input annulus AND THEIR ERRORS.  
;    The position and inner and outer
;    radii define the annulus.  Each pixel in the input image is then
;    assigned a weighting value which is its areal overlap between the pixel
;    and the annulus.  This weight varies from 0 to 1 and is the precise
;    analytic area of overlap (computed by pixwt.pro).  Of course, this
;    program doesn't do the computation for all pixels, only those near
;    the edge.
;
;    UNLIKE IN SUMANN, the moments are NOT into two components
;
; CATEGORY:
;    CCD data processing
; CALLING SEQUENCE:
;    Sumannerr, image, imageerr,xcen, ycen, $
;            inradius, outradius, back, backerr, totweight, $
;            sum, xmom, ymom, $
;            sumerr, xmomerr, ymomerr
;
; INPUTS:
;    image       : CCD image array.
;    imageerr    : std deviation of image
;    xcen,ycen   : Center of annulus.
;    inradius    : Radius of inner circle.
;    outradius   : Radius of outer circle.
;    back        : Background to subtract from each pixel.
;    backerr     : error in background
;
; OPTIONAL INPUT PARAMETERS:
;
; KEYWORD PARAMETERS:
;
; OUTPUTS:
;    totweight          : Area of annulus.
;    sum     : Sums of positive and negative pixels.
;    xmomm   : x moments relative to xcen.
;    ymom    : y moments relative to ycen.
;    sumerr  : formal error in sum
;    xmomerr : formal error in xmom
;    ymomerr : formal error in ymom
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
;    Ported by Doug Loucks, Lowell Observatory, 1992 Oct, from the
;    C-language version written by Marc Buie.
;    April, 1993. DWL. Replaced the inner FOR loop with vectors.
;    98/09/21, MWB, optimizations
;    98/09/29, MWB, fixed nasty bug introduced during optimization.
;    06/01/31, LAY, removed pos and neg sums, added errors, based on subann
;-
PRO Sumannerr, image, imageerr, xcen, ycen, inradius, outradius, $
               back, backerr, $
               totweight, sum, xmom, ymom, sumerr, xmomerr, ymomerr

   totweight = 0.0

   sum  = 0.0
   xmom = 0.0
   ymom = 0.0
   
   sumerr  = 0.0
   xmomerr = 0.0
   ymomerr = 0.0
   
   image_size = SIZE( image )

   r2 = outradius + 2
   npts = LONG( !PI * r2*r2 )
   x = INTARR( npts, /NOZERO )
   y = INTARR( npts, /NOZERO )
   flags = INTARR( npts )

   ;  As the pixel location vectors (x, y) are assembled, the vector flags
   ;  is constructed with the following values:
   ;
   ;     0  : Trivial.  pixwt not needed.
   ;     1  : Trivial.  pixwt not needed.
   ;    -1  : On or near outer circle only.
   ;    -2  : On or near inner circle only.
   ;    -4  : On or near both circles.
   ;
   ;  Then the flags are used to calculate the final weights for each pixel
   ;  location.
   ;
   IF ( 0 LE inradius ) AND ( inradius LT outradius ) THEN BEGIN
      ; Compute the y-range limits.
      r2 = outradius
      r1 = inradius
      r2sq = r2 * r2
      r1sq = r1 * r1
      yp = ycen - r2
      outy0 = FIX( yp )
      IF outy0+0.5 LT yp THEN outy0=outy0+1
      yp = ycen + r2
      outy3 = FIX( yp + 1 )
      IF outy3-0.5 GT yp THEN outy3=outy3-1

      base = 0
      FOR yp = outy0, outy3 DO BEGIN
         cgetrng, ycen, xcen, r2, yp, outx0, outx1, outx2, outx3
         cgetrng, ycen, xcen, r1, yp, inx0, inx1, inx2, inx3
         n = outx3 - outx0 + 1

         ;  Assume on outer circle.
         xx = outx0 + INDGEN( n )
         yy = REPLICATE( yp, n )
         f  = REPLICATE( -1, n )  ; -1 = on or near outer circle

         ;  Entirely inside outer circle.   (1)
         t = WHERE( outx1 LE xx AND xx LT outx2, count )
         IF count GT 0 THEN f[t] = 1

         ;  On inner circle.
         t = WHERE( (inx0 LE xx AND xx LT inx1) OR $
                    (inx2 LE xx AND xx LT inx3), count )
         IF count GT 0 THEN f[t]=f[t]-3

         ;  Entirely inside inner circle.  (0)
         t = WHERE( inx1 LE xx AND xx LT inx2, count )
         IF count GT 0 THEN f[t]=0
         
         x[ base : base+n-1 ] = xx
         y[ base : base+n-1 ] = yy
         flags[ base : base+n-1 ] = f
         base = base + n
      ENDFOR

      ; This weeds out all non-contributing pixels as well as the extra
      ;   storage (npts-base).
      z = WHERE( flags NE 0, count )
      IF count GT 0 THEN BEGIN
         x = x[z]
         y = y[z]
         flags = flags[z]
      ENDIF

      ; Keep only those pixels which are actually in the image array.
      xmin = MIN( x, MAX=xmax )
      ymin = MIN( y, MAX=ymax )
      IF xmin LT 0 OR xmax GE image_size[1] OR $
         ymin LT 0 OR ymax GE image_size[2] THEN BEGIN
         z = WHERE( x GE 0 AND x LE image_size[ 1 ] AND $
                    y GE 0 AND y LE image_size[ 2 ], count )
         IF count EQ 0 THEN RETURN
         x = x[z]
         y = y[z]
         flags = flags[z]
      ENDIF

      wts = replicate(1.0,n_elements(flags))

      z = WHERE( flags EQ -1, count )
      IF count GT 0 THEN wts[z] = pixwt(xcen,ycen,outradius,x[z],y[z]) > 0.

      z = WHERE( flags EQ -2, count )
      IF count GT 0 THEN wts[z] = (1.0 - pixwt(xcen,ycen,inradius,x[z],y[z])) > 0.

      z = WHERE( flags EQ -4, count )
      IF count GT 0 THEN BEGIN
         wts[z] = (pixwt( xcen, ycen, outradius, x[z], y[z] ) - $
                  pixwt( xcen, ycen, inradius,  x[z], y[z] ) ) > 0.
      ENDIF

      totweight = TOTAL( wts, /DOUBLE )

      val = image[ x, y ]
      vxw = ( val - back ) * wts
      dx = x - xcen
      dy = y - ycen
      dx2 = dx^2.d
      dy2 = dy^2.d

      sum  = TOTAL( vxw, /DOUBLE )
      xmom = TOTAL( vxw * dx, /DOUBLE )
      ymom = TOTAL( vxw * dy, /DOUBLE )

      var = (imageerr[ x, y ])^2.
      wts2 = wts^2.

      sumerr  = sqrt( total( var*wts2, /DOUBLE) + backerr^2 * totweight^2 )
      xmomerr = sqrt( total( var*wts2*dx2, /DOUBLE) + $
                      backerr^2 * TOTAL(wts*dx,/double)^2 )
      ymomerr = sqrt( total( var*wts2*dy2, /DOUBLE) + $
                      backerr^2 * TOTAL(wts*dy,/double)^2 )
      
   ENDIF

END