;+
; 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_zero.pro
;  Written  14 Jun 2007 Leslie Ann Young SwRI
;  Modified 13 Sep 2007 LAY
;  HISTORY OF comb_dark
;  2010-08-30 --- convert from comb_spe_zero
;-

pro  comb_dark, aarr, zero, zerr, zbad, err=earr, qual=barr, exparr=exparr, exptime=exptime, $
                bias=bias, berr=biaserr, dark=dark, derr=darkerr, 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(exparr) then exparr = replicate(1.,n)

    if not keyword_set(exptime) then exptime = mean(exparr)

    if not keyword_set(method) then method = 'average'

    case method of
       'average': begin
          zero = total(aarr,3)/n
          zerr = sqrt(total( reform(earr^2,nx,ny,n), 3)/n)
          zbad = ( total(barr,3) gt 0 )
;          znoise = mean(narr)
       end
       
       'bias+dark': begin
            
          i = 0
          w = 1/earr[*,*,i]^2.
          s = w
          sx = exparr[i] * w
          sy = aarr[*,*,i] * w
          sxx = exparr[i]^2. * w
          sxy = exparr[i] * aarr[*,*,i] * w

          for i = 1, n-1 do begin
             w = 1/earr[*,*,i]^2.
             s = s + w
             sx = sx + exparr[i] * w
             sy = sy + aarr[*,*,i] * w
             sxx = sxx + exparr[i]^2. * w
             sxy = sxy + exparr[i] * aarr[*,*,i] * w
          endfor
          
          del = s * sxx - sx^2
          
          if min(del) gt 0. then begin
          ; -----------------------
          ; normal equations for weighted fit to a line
          ; -----------------------
             bias = (sxx * sy - sx * sxy)/del
             dark = (s * sxy - sx * sy)/del
             biaserr = sqrt( sxx/del )
             darkerr = sqrt( s/del )
             cov = -sx/del
          endif else begin
          ; -----------------------
          ; del = 0 if one file or all times are equal.
          ; assume all of the signal is due to the bias
          ; -----------------------
             bias = sy/s
             dark = fltarr(nx,ny)
             biaserr = sqrt(1/s)
             darkerr = fltarr(nx,ny)
             cov = fltarr(nx,ny)
          endelse
          
          zero = bias + dark * exptime
          zerr = sqrt( biaserr^2 + exptime^2 * dark^2 + exptime * cov)
          zbad = ( total(barr,3) gt 0 )
          help, bias, dark
       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