;+
; NAME:
;  make_spe_cal1
;
; PURPOSE: (one line)
;  make a calibrated from spe files
;
; DESCRIPTION:
;  Given a filenames of SPE zero-flux files, make averaged flat-field file
;
; CATEGORY:
;  IMAGES
;  
; CALLING SEQUENCE:
;  make_spe_cal1,filename, cal, err, bad, $
;                    silent=silent, $
;                    gain=gain, znoise=znoise, $
;                    zfilename = zfilename, ffilename = ffilename, $
;                    caldir=caldir, outfilename = outfilename, $
;                    overrideoutfilename = overrideoutfilename
; INPUTS:
;  filename : filename of SPE image cube with zero-flux (dark+bias) 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 input file ------
;  
;  zfilename = zero flux file
;  znoise : noise per pixel in a zero-flux image (eg readnoise for
;    negligible dark), in DN -- overrides znoise in zfilename
;  
;  ffilename = flat field file
;  gain : electron/dn, used for calculating expected noise.  Default =
;     estimate from scatter.
;
;  -------- output file ------------
;  caldir : data directory for cal files (with terminal '/').  Default = 'cal/'
;  overrideoutfilename : if set (1), use outfilename as output file name.
;     Default = 0, construct the output filename from the input file name
;  outputfilename : if overrideoutfilename is set, this is the output
;      filename used
;  nowrite : if set, don't write the output file
;  
; OUTPUTS:
;  cal - the average zero-flux image
;  err - the error in the average zero-flux image (NOT error per
;              image per pixel)
;  bad - bad pixel image for avg (0 = good)
;
; OPTIONAL OUTPUTS:
;
; SIDE EFFECTS:
;  For each file, writes 
;  cal/<basefilename>_cal.fits
;
;  Primary Data Header PDU = D 
;  calibrated data = (raw - zero)/flat
;
;  TIME_REF
;  DATE-OBS 
;  TIME-OBS
;  STAR_REF
;  RA_STAR
;  DEC_STAR
;  OBSERVER
;  SITENAME= 'LOWELL-HALL'
;  SITEABBR= 'LWL-1.1'
;  LONGITUD= '111:32:10.6W'
;  LATITUDE= '35:05:48.7N'
;  ALTITUDE=        2.20900000000 
;  SITE_REF= 'p445.3_sites_v4.tk ' 
;
;ORIENT   =                 91.0 / orientation in degrees clockwise from E-W      
;
;  Extension 1 - gain array G
;  gain (e/dn) * flat
;
;  Extension 2 - zero noise array ZERR
;  sqrt( znoise^2 + zero_err^2 ) / flat
;  where znoise is the readnoise (or readoise + darknoise)
;  in a single pixel of a single frame, and zero_err
;  is the error in zero, from the file zfilename
; 
;  Extension 3 - flat error FERR
;  flat_err/flat
;  where flat_err is the error in the flatfield, from
;  the file ffilename
;
;  With D, G, ZERR, and FERR, we can calculate the error for
;  any frame
;  
;  err^2 = D/G + ZERR^2 + (D * FERR)^2
;
;  Extension 4 - bad pixel array
;  NX, by NY, 0 = good, 1 = bad
;
;  Extension 5 - bad image list
;  N, 0 = good, 1 = bad
;
; EXAMPLE:
;
; MODIFICATION HISTORY:
;  Written  Oct 2 2007 Leslie Ann Young SwRI
;-

pro make_spe_cal1,filename, cal, cerr, cbad, $
                  silent=silent, $
                  time_ref = time_ref, $
                  utc0 = utc0, ipp=ipp, $
                  date_obs = date_obs, time_obs = time_obs, $
                  crval3 = crval3, $
                  target = target, $
                  star = star, $
                  observer = observer, $
                  filter=filter, $
                  camera=camera, $
                  comments = comments, $
                  lc_abbr=lc_abbr, $
                  gain=gain, znoise=znoise, $
                  zfilename = zfilename, ffilename = ffilename, $
                  caldir=caldir, outfilename = outfilename, $
                  overrideoutfilename = overrideoutfilename, $
                  nowrite = nowrite

    if not keyword_set(silent) then silent = 0b
    if not keyword_set(caldir) then caldir = 'cal/'
    if not keyword_set(observer) then observer = ''
    
    ;--------------------------------
    ; Read the file using filename
    ;--------------------------------
    fullfn = filename[0]
    if not silent then begin
        print, 'make_spe_cal1: calibrating ', fullfn
    endif
    fn = file_basename(fullfn)
    len = strlen(fn)
    if strmatch(fn,'*.gz', /fold) then begin
        basefn = strmid(fn,0,len-7)
    endif else begin
        basefn = strmid(fn,0,len-4)
    endelse
    if not isfile(fullfn) then begin
        print, 'make_spe_cal1 : ', fullfn, ' does not exist'
        return
    endif
    read_princeton_gz,fullfn,cal,head=spehead
    nxnyn, cal, nx, ny, nimages

    ;--------------------------------
    ; ALLOCATE THE ERROR AND BADPIX ARRAYS HERE
    ;--------------------------------
    cerr = fltarr(nx,ny,nimages)
    cbad = bytarr(nx,ny,nimages)

    ;--------------------------------
    ; CAMERA INFORMATION
    ;--------------------------------
    if not keyword_set(filter) then filter = ''
    if not keyword_set (camera) then begin
        camera = ''
        case spehead.dettype of
            94 : camtype='micromax'
            0 : camtype='photonmax'
            else: camtype='unknown'
        endcase
    endif else begin
        case strupcase(camera) of
            'DOROTHEA': camtype='micromax'
            'ANSEL': camtype='micromax'
            'WALKER': camtype='micromax'
            'HENRI': camtype='micromax'
            'DOC': camtype='photonmax'
            'GJON': camtype='photonmax'
            'EADWEARD': camtype='photonmax'
            else: camtype = 'unknown'
        endcase
    endelse
    if camtype eq 'unknown' then begin
        print, 'make_spe_cal1 : WARNING unknown camera type'
    endif
    if not silent then print, 'CAMERA = ', camera, ', CAMTYPE = ', camtype

    ;--------------------------------
    ; TIMING INFORMATION
    ;--------------------------------

    ; get the start time
    ; if it's not passed, then get the time from the header
    ;   and set time_ref to 'SPE header'
    if not keyword_set(utc0) then begin
        if keyword_set(ipp) then begin
            print, 'make_spe_cal1 : WARNING interpulse period set, but'+$
              ' utc0 was not passed'
        endif
        date_spe = spehead.date
        time_spe = spehead.EXPERIMENTTIMEUTC
        utc0 = strmid(date_spe,5,4) + ' ' + $
          strmid(date_spe,2,3) + ' ' + $
          strmid(date_spe,0,2) + ' ' + $
          strmid(time_spe,0,2) + ':' + $
          strmid(time_spe,2,2) + ':' + $
          strmid(time_spe,4,2) 
        time_ref = 'SPE header, set from data acquisition computer'
    endif else begin
        if not keyword_set(time_ref) then time_ref = ''
    endelse
    et0 = utc2et(utc0)
    date_obs = strmid(et2utc(et0, fo='ISOC'),0,10)
    time_obs = strmid(et2utc(et0, fo='ISOC'),11,12)

    if keyword_set(ipp) then begin
        ; ipp set - triggered exposure
        case camtype of 
            'micromax' : begin
                exptime = spehead.exp_sec
                crpix3 = 1
                crval3 = spehead.exp_sec / 2. 
                cdelt3 = ipp
            end

            'photonmax' : begin
                exptime = ipp
                crpix3 = 1
                crval3 = spehead.exp_sec - ipp/ 2.
                cdelt3 = ipp
            end

            'unknown' : begin
                exptime = ipp
                crpix3 = 1
                crval3 = spehead.exp_sec / 2.
                cdelt3 = ipp
            end
            
        endcase
    endif else begin
        exptime = spehead.exp_sec
        crpix3 = 1
        crval3 = exptime/2.
        if spehead.readouttime/1e3 lt exptime then begin
            ; frame transfer
            cdelt3 = exptime
        endif else begin
            cdelt3 = exptime + spehead.readouttime/1e3
        endelse
    endelse
    if not silent then begin
        print, 'EXPTIME = ', float(exptime)
    end

    ;--------------------------------
    ; STAR INFORMATION
    ;--------------------------------
    if keyword_set(star) then begin
        star_ref = star.poscat
        ra_star = rastrf(star.ra)
        dec_star = decstrf(star.dec)
    endif else begin
        star_ref = ''
        ra_star = ''
        dec_star = ''
    endelse

    ;--------------------------------
    ; SITE INFORMATION
    ;--------------------------------
    sitename = ''
    longitud = ''
    latitude = ''
    altitude = ''
    site_ref = ''
    if keyword_set(lc_abbr) then begin
        kernels = naif_loadedkernels()
        itk = where(strmatch(kernels, '*.tk'), ntk)
        if ntk eq 0 then begin
            print, 'make_spe_cal1 ; WARNING ', $
              'no text kernel loaded'
        endif else begin
            for i = 0, ntk - 1 do begin
                site_ref = site_ref + kernels[itk[i]] + ' '
            endfor
            ; get the naifid and site name
            cspice_bodn2c, lc_abbr, obsid, found
            if found ne 1 then begin
                print, 'make_spe_cal1 : WARNING ', $
                  lc_abbr, ' not found'
            endif else begin
                cspice_bodc2n, obsid, site, found
                if found ne 1 then begin
                    print, 'make_spe_cal1 : WARNING ', $
                      obsid, ' not found'
                endif else begin
                    sitename = site
                    deg = !dpi/180.d
                    cspice_gdpool, site+'_LATLON', 0, 3, cvals, found
                    latitude = latstrf(cvals[0] * deg)
                    longitud = lonstrf(cvals[1] * deg)
                    altitude = cvals[2]
                endelse
            endelse
        endelse
    endif else begin
        lc_abbr = ''
    endelse

    ;--------------------------------
    ; Target
    ;--------------------------------
    if not keyword_set(target) then target = ''

    ;--------------------------------
    ; Subtract the zero-flux
    ;--------------------------------
    if keyword_set(zfilename) then begin
        if zfilename eq 'auto' then begin
            zfn = caldir + strtrim(string(round(exptime*1000)),2)+'_zero.fits'
        endif else begin
            zfn = zfilename
        endelse
        zero = get_spe_zero(zfn, zerr, zbad, znoise_zfn, $
          silent=silent, ainfo=h)
    end else begin
        zfn = ''
        zero = fltarr(nx, ny)
        zerr = fltarr(nx, ny)
        zbad = bytarr(nx,ny)
        znoise_zfn = 0.
    endelse
    cal = float(cal)
    for i = 0, nimages-1 do begin
        cal[*,*,i] = cal[*,*,i] - zero
    endfor
    if n_elements(znoise) eq 0 then begin
        znoise = znoise_zfn
    endif
    
    ;--------------------------------
    ; Divide by the flat field
    ;--------------------------------
    ffn = ''
    flat = fltarr(nx, ny) + 1.
    ferr = fltarr(nx, ny)
    fbad = bytarr(nx,ny)
    gain_ffn = 1.
    if keyword_set(ffilename) then begin
        if  ffilename ne '' then begin
            ffn = ffilename
            flat = get_spe_flat(ffn, ferr, fbad, gain_ffn, $
                                silent=silent, ainfo=spehead)
        endif
    endif
    if n_elements(gain) eq 0 then begin
        gain = gain_ffn
    endif
    for i = 0, nimages-1 do begin
        cal[*,*,i] = cal[*,*,i]/flat
    endfor
    
    ;--------------------------------
    ; Make the error arrays
    ;--------------------------------
    gain_arr = gain * flat
    znoise_arr = sqrt( znoise^2 + zerr^2 ) / flat
    ferr_arr = ferr / flat

    ;cerr = fltarr(nx,ny,nimages)
    for i = 0, nimages-1 do begin
        cerr[*,*,i] = $
          (cal[*,*,i]>0)/gain_arr + $
          znoise_arr^2 + $
          (cal[*,*,i] * ferr_arr)^2
    endfor
    if min(cerr) lt 0. then stop
    cerr[*,*,*] = sqrt(cerr[*,*,*])

    ;--------------------------------
    ; Make the bad pixel arrays
    ;--------------------------------
    bad_arr = (zbad ne 0) or (fbad ne 0)
    bad_list = bytarr(nimages)
    if n_elements(imbad) gt 0 then begin
        bad_list[imbad] = 1
    endif

    ; this is maybe too esoteric?
    ;cbad = $
    ;  reform(reform(bad_arr eq 0,nx*long(ny)) # (bad_list eq 0), nx, ny, n) $
    ;  eq 0
    ;cbad = bytarr(nx,ny,nimages)
    for i = 0, nimages-1 do begin
        cbad[*,*,i] = (bad_arr ne 0) or (bad_list[i] ne 0)
    endfor

    if not keyword_set(nowrite) then begin

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

        if keyword_set(overridefilename) then begin
            flatfn = outfilename
        endif else begin
            flatfn = caldir + basefn+'_cal.fits'
            outfilename = flatfn
        endelse
        if not silent then print, flatfn

    ;---------------------------------------------
    ; write it out
    ;---------------------------------------------
        mkhdr, header, cal, /EXTEND
        header = spe_head2fits(spehead, infits = header)
        sxaddpar, header, 'FIRSTOBS', '', $
          'Placeholder for adding comments in order', $
          before='FIRSTSPE'
        sxaddpar, header, 'COMMENT', $
          'Created by make_spe_cal1', $
          before = 'FIRSTOBS'
        sxaddpar, header, 'COMMENT', $
          'Primary Data Unit is the calibrated file', $
          before = 'FIRSTOBS'
        sxaddpar, header, 'COMMENT', $
          'Extensions are used by get_spe_cal to make ',$
          before = 'FIRSTOBS'
        sxaddpar, header, 'COMMENT', $
          ' error and bad-pixel mask arrays ', $
          before = 'FIRSTOBS'
        for i = 0, n_elements(comments)-1 do begin
            sxaddpar, header, 'COMMENT', comments[i], $
              before='FIRSTSPE'
        endfor
        sxaddpar, header, 'TIME_REF', time_ref, $
          before='FIRSTSPE'
        sxaddpar, header, 'DATE-OBS', date_obs, $
          before='FIRSTSPE'
        sxaddpar, header, 'TIME-OBS', time_obs, $
          before='FIRSTSPE'
        sxaddpar, header, 'CDELT3', cdelt3, $
          'Interval between exposures', $
          before='FIRSTSPE'          
        sxaddpar, header, 'CRPIX3', crpix3, $
          before='FIRSTSPE'          
        sxaddpar, header, 'CRVAL3', crval3, $
          'Offset of midtime of 1st exposure from TIME-OBS', $
          before='FIRSTSPE'          
        sxaddpar, header, 'EXPTIME', $
          exptime, 'Exposure time in seconds', $
          before='FIRSTSPE'
        sxaddpar, header, 'FILTER', filter, $
          before='FIRSTSPE'
        sxaddpar, header, 'CAMERA', camera, $
          before='FIRSTSPE'
        sxaddpar, header, 'CAMTYPE', camtype, $
          before='FIRSTSPE'
        sxaddpar, header, 'OBSERVER', observer, $
          before='FIRSTSPE'
        sxaddpar, header, 'RA_STAR', ra_star, $
          before='FIRSTSPE'
        sxaddpar, header, 'DEC_STAR', dec_star, $
          before='FIRSTSPE'
        sxaddpar, header, 'STAR_REF', star_ref, $
          before='FIRSTSPE'
        sxaddpar, header, 'TARGET', target, $
          before='FIRSTSPE'
        sxaddpar, header, 'SITENAME', sitename, $
          before='FIRSTSPE'
        sxaddpar, header, 'LC_ABBR', lc_abbr, $
          before='FIRSTSPE'
        sxaddpar, header, 'LONGITUD', longitud, 'East longitude', $
          before='FIRSTSPE'
        sxaddpar, header, 'LATITUDE', latitude, 'Latitude', $
          before='FIRSTSPE'
        sxaddpar, header, 'ALTITUDE', altitude, 'altitude in km', $
          before='FIRSTSPE'
        sxaddpar, header, 'SITE_REF', site_ref, $
          before='FIRSTSPE'
        sxaddpar, header, 'INFILE', fullfn, before='FIRSTSPE'
        sxaddpar, header, 'ZEROFILE', strend(zfn,68), before='FIRSTSPE'
        sxaddpar, header, 'FLATFILE', strend(ffn,68), before='FIRSTSPE'
        sxaddpar, header, 'OBJECT', $
          'calibrated file', before='FIRSTSPE'
        sxaddpar, header, $
          'READN_DN', znoise, $
          'noise per pixel, read noise + dark noise (DN)', $
          before='FIRSTSPE'
        sxaddpar, header, $
          'GAIN', gain, $
          'gain, electrons per DN', $
          before='FIRSTSPE'
        sxaddpar, header, 'LASTOBS', '', $
          'Placeholder for adding comments in order', $
          before='FIRSTSPE'


        WRITEFITS, flatfn, cal, header

        mkhdr, header, gain_arr, /EXTEND 
        sxaddpar, header, 'OBJECT', 'gain', $
          'electrons per calibrated DN'
        WRITEFITS, flatfn, gain_arr, /append
        
        mkhdr, header, znoise_arr, /EXTEND 
        sxaddpar, header, 'OBJECT', 'read noise', $
          'calibrated DN per pixel'
        WRITEFITS, flatfn, znoise_arr, /append
        
        mkhdr, header, ferr_arr, /EXTEND 
        sxaddpar, header, 'OBJECT', 'fractional error in flat'
        WRITEFITS, flatfn, ferr_arr, /append
        
        mkhdr, header, bad_arr, /EXTEND 
        sxaddpar, header, 'OBJECT', 'array of bad pixels'
        WRITEFITS, flatfn, bad_arr, /append
        
        mkhdr, header, bad_list, /EXTEND 
        sxaddpar, header, 'OBJECT', 'list of bad frames'
        WRITEFITS, flatfn, bad_list, /append
        
    endif

end