;+
; NAME:
;  comb_dark
;
; PURPOSE: (one line)
;  combine zero-flux (dark+bias) averaged image from image cube
;
; DESCRIPTION:
;  
;
; CATEGORY:
;  IMAGES
;  
; CALLING SEQUENCE:
;
; INPUTS:
;  cube : 
;
; OPTIONAL INPUT PARAMETERS:
;
; OUTPUTS:
;
; OPTIONAL OUTPUTS:
;
; EXAMPLE:
;
; MODIFICATION HISTORY:
;  HISTORY OF comb_spe_flat.pro
;  Written  14 Jun 2007 Leslie Ann Young SwRI
;  Modified 13 Sep 2007 LAY
;  HISTORY OF comb_dark
;  2010-09-22 --- convert from comb_spe_flat
;-

pro  comb_flat, aarr, flat, ferr, fbad, err=earr, qual=barr, method=method
    
    nxnyn, aarr, nx,ny,n

    ; make the larger arrays
    if not keyword_set(earr) then begin
       earr = fltarr(nx,ny,n) + 1.
    endif
    if not keyword_set(barr) then begin
       barr = bytarr(nx,ny,n)
    endif

    ; reform to cubes if needed
    aar2d = (size(aarr, /n_dim) eq 2)
    if aar2d then aarr = reform(aarr,nx,ny,1)
    ear2d = (size(earr, /n_dim) eq 2)
    if ear2d then earr = reform(earr,nx,ny,1)
    bar2d = (size(barr, /n_dim) eq 2)
    if bar2d then barr = reform(barr,nx,ny,1)
    
    if not keyword_set(method) then method = 'average'

    case method of
       'average': begin
          wght = reform(float((barr eq 0)), nx, ny, n)
          wghttot = total(wght, 3)
          gd = where(wghttot gt 0, ngd)
          if ngd eq 0 then begin
             print, 'comb_flat : all points are bad'
             flat = fltarr(nx,ny)
             ferr = fltarr(nx,ny)
             fbad = bytarr(nx,ny) + 1
             return
          endif
          wghttot = wghttot > min(wghttot[gd])
          for i = 0, n-1 do wght[*,*,i] = wght[*,*,i] / wghttot
          flat = total( reform(aarr * wght,nx,ny,n) ,3) 
          ferr = sqrt(total( reform(wght^2 * earr^2,nx,ny,n) , 3) )
          fbad = ( total(barr,3) gt 0 )
;          gain = total( total(total(wght,1),1) * garr ) / total(wght)
       end
       
    endcase

    ; reform from cubes if needed
    if aar2d then aar = reform(aarr,nx,ny)
    if ear2d then aar = reform(earr,nx,ny)
    if bar2d then aar = reform(barr,nx,ny)
    
 end