pro sxt_sumxy,index,data,sumX,sumY,unc_data,satpix,silent=silent
;+
;  NAME:
;    sxt_sumxy
;  PURPOSE:
;    Sum SXT X-ray images by sumX by sumY pixels.
;    Preserve decompression uncertainty and saturated pixel information.
;
;    Warning:  This routine clobbers the input variables.
;
;  CALLING SEQUENCE:
;    sxt_sumxy,index,data,sumX
;    sxt_sumxy,index,data,sumX,sumY
;    sxt_sumxy,index,data,sumX,sumY,unc_data
;    sxt_sumxy,index,data,sumX,sumY,unc_data,satpix
;  INPUTS/OUTPUTS:
;    index	= SXT index structure
;    data	= Image cube			(overwrites the input variables)
;  INPUTS:
;    sumX	= Number of pixels to sum together
;  OPTIONAL INPUTS:
;    sumY	= If missing, sumY will be set equal to sumX
;  OPTIONAL INPUT/OUTPUTS
;    unc_data	= Decompression uncertainties
;    satpix	= Saturated pixels (can be a sparse structure)
;		  (Max value of satpix is 255.)
;  RESTRICTIONS:
;    data must be decompressed and 2-d or 3-d
;    Currently, the index record is not updated.
;
;  MODIFICATION HISTORY:
;    9-Mar-93, J. R. Lemen, LPARL
;   16-mar-93, JRL, Fixed up round-off scheme for unc_data
;-

; -------------------------------------------------------------
;  Check to see if data is byte-type.
; -------------------------------------------------------------
sz = size(data) & sz_type = sz(sz(0)+1)
if sz_type eq 1 then begin
  if not keyword_set(silent) then begin
      print,'****  Error in sxt_sumxy:  data must be decompressed (not byte type)'
      tbeep & help,data & print,'      No summation was performed'
  endif
  return
endif

; -------------------------------------------------------------
;  Check that data is 2-d or 3-d
; -------------------------------------------------------------

if (sz(0) ne 2) and (sz(0) ne 3) then begin
    print,'**** Error in sxt_sumxy:  data must be 2-d or 3-d'
    tbeep & help,data
    return
endif

; -------------------------------------------------------------
;  Check that sumX and sumY are integral multiples
; -------------------------------------------------------------

n_X = sz(1) & n_Y = sz(2) & nimg = n_elements(data(0,0,*))
if n_params() lt 4 then sumY = sumX


if ((N_x mod sumX) ne 0) or ((N_y mod sumY) ne 0) then begin
   print,'***  Warning in sxt_sumxy:  ' & tbeep
   print,'     Summation must be an integral factor of input array size'
   print,'     No summation performed'
   return
endif

if (sumX gt 1) or (sumY gt 1) then begin
      if not keyword_set(silent) then 	$
	print,'***  sxt_sumxy: Sum mode (X,Y) = ',strtrim(sumX,2),', ',strtrim(sumY,2)
      NsumX = N_x / sumX & NsumY = N_y / sumY
;  Make the output array integer*4 if summing 4x4 pixels or more pixels:
      sz(1) = NsumX & sz(2) = NsumY
      if (sumX gt 2) and (sumY gt 2) and (sz(sz(0)+1) eq 2) then 	$
			sz(sz(0)+1) = 3			; Force to be long
      new_data = make_array(size=sz)
;  Set up new saturated and uncertainty arrays:
      if n_params() ge 5 then begin
	sz1 = size(unc_data) & sz1(1) = NsumX & sz1(2) = NsumY
	unc_data2 = make_array(size=sz1)
; Round off if the unc_data is not floating:
        if sz1(sz1(0)+1) le 3 then round_off1 = 0.5 else round_off1 = 0.
      endif
      if n_params() ge 6 then begin
	sz2 = size(satpix)   & sz2(1) = NsumX & sz2(2) = NsumY
	satpix2 =   make_array(size=sz2)  
      endif
      npix = float(sumX*sumY)			; Number of pixels in macro pix
;  Begin the summation loop here (how do you do this without loops?)
      for i=0,nimg-1 do begin			; Loop over images
        for j=0,NsumX-1 do begin		; Loop over columns
           j1 = j*sumX & j2 = j1 + sumX-1
           for k=0,NsumY-1 do begin		; Loop over rows
              k1 = k*sumY & k2 = k1 + sumY-1	; Extract the data to sum
              new_data(j,k,i) = total(data(j1:j2,k1:k2,i))
              if n_params() ge 5 then unc_data2(j,k,i) =        $
                 sqrt(total(float(unc_data(j1:j2,k1:k2,i))^2))/npix + round_off1
              if n_params() ge 6 then satpix2(j,k,i) = 		$
						total(satpix(j1:j2,k1:k2,i))<255
           endfor
        endfor
      endfor					; Loop over final image
      data = 0					; Clear memory
      data = new_data				; Return in data
      new_data = 0				; Clear memory
      if n_params() ge 5 then unc_data = unc_data2	; Uncertainties
      if n_params() ge 6 then satpix = satpix2		; saturated pixels
endif						; (sumX gt 1) or (sumY gt 1)

end