;+
; NAME:
;  comb_spe_flat
;
; PURPOSE: (one line)
;  combine spe_flat files
;
; DESCRIPTION:
;  combine spe_flat files
;
; CATEGORY:
;  SPE
;  
; CALLING SEQUENCE:
;  comb_spe_flat, ffilenames, flat, ferr, fbad, gain, $
;                    silent=silent, $
;                    aexptime=aexptime, ainfo=ainfo
; INPUTS:
;  ffilenames : filenames of SPE flatfield average images
;
; OPTIONAL INPUT PARAMETERS:
;  -------- behavior of the routine ------------
;  silent : 1 if the routine is to perform silently. Default = 0,
;           print out statements
; 
;  -------- description of the file the zero file will be applied to ------
;  aexptime : exposure time in seconds
;  ainfo: string (filename), array of strings (FITS header) or 
;    structure (SPE header)
;  
; OUTPUTS:
;  flat - the average flat field image
;  ferr - the error in the average flat image 
;  fbad - bad pixel image for avg (0 = good)
;  gain - electrons per ADU
;
; SIDE EFFECTS:
;
; EXAMPLE:
;
; MODIFICATION HISTORY:
;  Written  14 Jun 2007 Leslie Ann Young SwRI
;  Modified 13 Sep 2007 LAY
;-

pro  comb_spe_flat,ffilenames, flat, ferr, fbad, gain, $
                   silent=silent, ainfo=ainfo, $
                   method=method, clip=clip, $
                   caldir=caldir, outfilename = outfilename, pre=pre, $
                   overridefilename = overridefilename, $
                   nowrite = nowrite
                   

    n = n_elements(ffilenames)
    if n eq 1 then ffilenames = reform(ffilenames, 1)

    ; read the first one and make the larger arrays
    i = 0
    fn = ffilenames[i]
    a = get_spe_flat(fn, e, b, fitshead=h, $
      silent=silent, ainfo=ainfo)
    nxnyn, a, nx0, ny0

    aarr = fltarr(nx0,ny0,n)
    earr = fltarr(nx0,ny0,n)
    barr = fltarr(nx0,ny0,n)
    harr = strarr( n_elements(h), n)

    if n eq 1 then begin
        aarr = reform(aarr, nx0,ny0,n)
        earr = reform(earr, nx0,ny0,n)
        barr = reform(barr, nx0,ny0,n)
        harr = reform(harr, n_elements(h), n)
    endif


    for i = 0, n-1 do begin
        fn = ffilenames[i]
        if not isfile(fn) then begin
            print, 'comb_spe_flat : ', fn, ' not found'
            stop
        endif
        a = get_spe_flat(fn, e, b, fitsh=h, silent=silent)
        nxnyn, a, nx, ny
        if nx eq nx0 and ny eq ny0 then begin
            aarr[*,*,i] = a
            barr[*,*,i] = b
            earr[*,*,i] = e
            harr[*,i] = h
        endif
    end
    garr = fltarr(n) ; gain
    for i=0,n-1 do garr[i] = sxpar(harr[*,i],'GAIN')

    if keyword_set(clip) then begin
        barr = barr or (aarr lt clip[0]) or (aarr gt clip[1])
    end

    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_spe_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

        'weight_by_nimage': begin
            nimages = lonarr(n)
            for i = 0, n-1 do nimages[i] = sxpar(harr[*,i],'NIMAGES')
            wght = fltarr(nx,ny,n)
            for i = 0, n-1 do wght[*,*,i] = nimages[i] / total(nimages)
            wght = wght * (barr eq 0)
            wghttot = total(wght, 3)
            gd = where(wghttot gt 0, ngd)
            if ngd eq 0 then begin
                print, 'comb_spe_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

        'weight_by_error': begin
            wght = 1./earr^2
            wght = wght * (barr eq 0)
            wghttot = total(wght, 3)
            gd = where(wghttot gt 0, ngd)
            if ngd eq 0 then begin
                print, 'comb_spe_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

        'median' : begin
            wght = float((barr eq 0))
            wghttot = total(wght, 3)
            gd = where(wghttot gt 0, ngd)
            if ngd eq 0 then begin
                print, 'comb_spe_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) 
            flat = fltarr(nx,ny)
            for i = 0, nx-1 do begin
                for j = 0, ny-1 do begin
                    gd = where (barr[i,j,*] eq 0, ngd)
                    if ngd gt 0 then begin
                        flat[i,j] = median(aarr[i,j,gd], /even)
                    endif
                endfor
            endfor
            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

        'medianw' : begin
            wght = float((barr eq 0))
            wghttot = total(wght, 3)
            gd = where(wghttot gt 0, ngd)
            if ngd eq 0 then begin
                print, 'comb_spe_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 = fltarr(nx,ny)
            for i = 0, nx-1 do begin
                for j = 0, ny-1 do begin
                    gd = where (barr[i,j,*] eq 0, ngd)
                    if ngd gt 0 then begin
                        flat[i,j] = medianw(aarr[i,j,gd], wght[i,j,gd])
                    endif
                endfor
            endfor
            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

    fbad = ( total( reform(barr eq 0,nx,ny,n), 3) eq 0 ) ; NONE are good 

    if not keyword_set(nowrite) then begin

    ;--------------------------------
    ; Get or create the output filename
    ;--------------------------------

        if not keyword_set(caldir) then caldir = 'cal/' 
        if not keyword_set(pre) then pre = 'super' 
        if keyword_set(overridefilename) then begin
            flatfn = outfilename
        endif else begin
            flatfn = caldir + pre + '_flat.fits'
            outfilename = flatfn
        endelse
        if not silent then print, flatfn
    
    ;---------------------------------------------
    ; write it out
    ;---------------------------------------------

        mkhdr, hdr, flat, /EXTEND
        iend = where( strmid(hdr,0,8) eq 'END     ')
        header = spe_head(h, /justspe) ; the lines between FIRSTSPE and LASTSPE
        header = [ hdr[0:iend-1], header, hdr[iend] ]
        sxaddpar, header, 'OBJECT', $
          'median image of flat-field files', before='FIRSTSPE'
        sxaddpar, header, $
          'NFILES', n, 'Number of files combined', $
          before='FIRSTSPE'
        for i = 0, n-1 do begin
            keyword = 'INFILE' + strtrim(string(i),2)
            sxaddpar, header, keyword, ffilenames[i], before='FIRSTSPE'
        endfor

        sxaddpar, header, $
          'GAIN', gain, $
          'Derived gain, (e/DN)', $
          before='FIRSTSPE'
        sxaddpar, header, $
          'COMBMETH', method, 'Method for combining the flat-field files', $
          before='FIRSTSPE'
        WRITEFITS, flatfn, flat, header

        mkhdr, header, ferr, /EXTEND 
        sxaddpar, header, 'OBJECT', 'error in median image'
        WRITEFITS, flatfn, ferr, /append
        
        mkhdr, header, fbad, /EXTEND
        sxaddpar, header, 'OBJECT', 'bad pixel list (0=good)'
        WRITEFITS, flatfn, fbad, /append

    endif


end