;+
; NAME: ra_mike_pipeline.pro
;
; PURPOSE:
;   Run the Mike pipeline (Rosetta Alice) to process engineering level
;   data into science level data
;
; CATEGORY: data procesing
;
; CALLING SEQUENCE:
;  ra_mike_pipeline, FileList, alldata = alldata, outdir = outdir, Verbose=Verbose,
;    deadflag = deadflag, darkflag = darkflag, flatflag = flatflag, aeffflag = aeffflag,
;    wavefile = wavefile, darkfile = darkfile, flatfile = flatfile, aefffile = aefffile,
;    badflag = badflag, badfile = badfile, rectflag = rectflag, logfile = logfile,
;    kernelfile = kernelfile, loadkernels = loadkernels, movefiles2sci = movefiles2sci
;
; OPTIONAL INPUTS:
;   filelist -- if supplied, a string or string array giving the
;               filenames of the files to be processed.  The string
;               may be a valid UNIX wildcard.
;   darkfile -- filepath to dark count file
;   flatfile -- filepath to flatfield file
;   aefffile -- filepath to effective area calibration file
;   
; KEYWORD PARAMETERS:
;
;   deadflag -- Default behavior is to apply deadtime correction to
;               all data. If set to 0, deadtime correction is not
;               performed.
;
;   darkflag -- Default behavior is to apply dark counts correction to
;               all data. If set to 0, dark count correction is not
;               performed.
;
;   flatflag -- Default behavior is to apply flatfield correction to
;               all data. If set to 0, flatfield correction is not
;               performed.
;
;   aeffflag -- Default behavior is to apply effective area
;               calibration to all data. If set to 0, the effective
;               area calibration  is not  performed.
;
;   badflag  -- Set this keyword to apply the bad pixel
;               mask. CURRENTLY, THIS KEYWORD IS NOT OPERATIONAL
;
;   rectflag -- Set this keyword to apply the rectification correction
;               to the data. CURRENTLY, THIS KEYWORD IS NOT OPERATIONAL
;
;   alldata -- set this keyword to process all eng files in the
;              rosetta data directory.
;
;   ra_data_dir -- set this keyword to the fully-qualified path to the
;                  rosetta data directory. If not provided, the
;                  current directory is used.
;
;   outdir -- set this keyword to the directory in which to output
;             files.
;
;   movefiles2sci --  If this keyword is provided, data are written
;             to the /data/sci/YYYY/MM/ directory, where YYYY and MM
;             are the appropriate Year and Month, respectively, of the
;             data file.
;
; OUTPUTS: The pipeline code produces science level data files
;
; COMMON BLOCKS:
;   mike_data -- contains a structure with various flags for processing
;
; SIDE EFFECTS:
;   none
;
;
; RESTRICTIONS:
;
;
;
; PROCEDURE:
;
;
;
; EXAMPLE:
;
;
;
; MODIFICATION HISTORY:
; The code has been completely rewritten to simplify the Rosetta Alice
; pipeline and significantly improve performance. After this rewrite,
; THE PIPELINE IS NOT COMPATIBLE WITH PREVIOUS VERSIONS.
; Version Date: 
;   v4. 2007 Feb 22  AJS
;   v5. 2007 Jun 13  AJS Added et column to pixel list table
;-

pro ra_mike_pipeline, FileList, alldata = alldata, outdir = outdir, Verbose=Verbose, $
  deadflag = deadflag, darkflag = darkflag, flatflag = flatflag, aeffflag = aeffflag, $
  wavefile = wavefile, darkfile = darkfile, flatfile = flatfile, aefffile = aefffile, $
  badflag = badflag, badfile = badfile, rectflag = rectflag, logfile = logfile, $
  kernelfile = kernelfile, loadkernels = loadkernels, movefiles2sci = movefiles2sci

COMMON MIKE_DATA, mikeflags, logmessages

IF n_params() EQ 0 AND NOT keyword_set(alldata) THEN BEGIN
  print, 'ra_mike_pipeline, FileList, alldata = alldata, outdir = outdir, Verbose=Verbose,'
  print,'  deadflag = deadflag, darkflag = darkflag, flatflag = flatflag, aeffflag = aeffflag,'
  print,'  wavefile = wavefile, darkfile = darkfile, flatfile = flatfile, aefffile = aefffile,'
  print,'  badflag = badflag, badfile = badfile, rectflag = rectflag, logfile = logfile,'
  print,'  kernelfile = kernelfile, loadkernels = loadkernels, movefiles2sci = movefiles2sci'
RETURN
ENDIF

; Default behavior is to apply corrections to the data files
IF n_elements(deadflag) EQ 0 THEN deadflag = 1
IF n_elements(darkflag) EQ 0 THEN darkflag = 1
IF n_elements(aeffflag) EQ 0 THEN aeffflag = 1

IF n_elements(flatflag) EQ 0 THEN flatflag = 0 ; DEFAULT IS NOW NO FLATFIELDING!
IF n_elements(badflag) EQ 0 THEN badflag = 0
IF n_elements(rectflag) EQ 0 THEN rectflag = 0

IF NOT keyword_set(kernelfile) THEN kernelfile = $
  '/private/var/automount/net/web/d0/web/webman/ralice/private/data/' + $
  'kernels/rosetta_kernels_argonath.txt'

; These directories and file references will need to be modified to
; match the directory structure on the local machine. 
calibration_dir = '/private/var/automount/net/web/d0/web/webman/ralice/private/data/cal/'
data_dir = '/private/var/automount/net/web/d0/web/webman/ralice/private/data/'

; Establish the error handler
;catch, error_status

; Wavelength Calibration file
IF NOT keyword_set(wavefile) THEN wavefile = calibration_dir + $
  'wavelength/RA_WAVE_003.FIT'
row_offset = readfits(wavefile, waveheader, /silent)

; Dark Correction File
IF NOT keyword_set(darkfile) THEN darkfile = calibration_dir + $
  'dark/RA_DARK_003.FIT'
IF keyword_set(darkflag) THEN BEGIN
   dark = readfits(darkfile, darkhdr, /silent)
   darkexp = sxpar(darkhdr, 'exptime')
   darkerror = sqrt(dark * darkexp) / darkexp
ENDIF

; Flatfield Correction file
IF NOT keyword_set(flatfile) THEN flatfile = calibration_dir + $
  'flat/RA_FLAT_002.FIT'
IF keyword_set(flatflag) THEN flat = readfits(flatfile, flathdr, /silent)

; Effective area Calibration file
IF NOT keyword_set(aefffile) THEN aefffile = calibration_dir + $
  'aeff/RA_AEFF_004.TAB'
IF keyword_set(aeffflag) THEN readcol, aefffile, aeff_wave, aeff, format = 'F,F', /silent

;;;
;   Set up the various keywords and flags.
;
IF n_elements(mikeflags) EQ 0 THEN $
  MikeFlags = {Version:'4.0 [2007 Feb 22]', $
               Mode:'', $
               Verbose: keyword_set(verbose) ? verbose : 0 , $
               Log:1, $
               Error:0, $
               StartTime:systime(1), $
               kernel:0, $
               nogeometry:0, $
               apdoor:0}

Time0 = systime(1)

; Start the log file
mike_log, mess, /start, /init, /time, verbose = 1

; Get the files to process.
; Use FILE_SEARCH to allow the user to enter wildcards in the filelist string
IF keyword_set(alldata) THEN BEGIN
      filelist = file_search(data_dir + 'eng/', 'ra_*eng.fit', count = nfiles)
   ENDIF ELSE BEGIN
      nfiles = n_elements(filelist)
      IF (nfiles eq 1) THEN filelist = file_search(filelist, count=nfiles)
   ENDELSE

IF (nfiles EQ 0) THEN BEGIN
   mike_log, err=2, 'No files to process.'
   RETURN
ENDIF


; check to see that the CSPICE ICY DLM has been registered with IDL
; either at startup (the RSI preferred way of doing this) or via the
; DLM_REGISTER procedure.
help,/dlm,/brief,out=dlm_help
foo = where(strpos(dlm_help, 'ICY') GE 0, icy_loaded)
IF icy_loaded EQ 0 THEN BEGIN
   mike_log, 'The CSPICE ICY DLM has not been registered with IDL, ' + $
             'so geometrical information cannot be determined.', errnum = 1
   mikeflags.nogeometry = 1
ENDIF

; The kernels required for SPICE to work have not been loaded, so load
; them now. This is accomplished via a SPICE metakernel file
IF NOT keyword_set(mikeflags.kernel) OR keyword_set(loadkernels) THEN BEGIN
   CSPICE_FURNSH, kernelfile
   mikeflags.kernel = 1
ENDIF

;;;
;   Go through each file and process it.
;   If pixel list data, then reform each time group of data into a histogram
; image and process it that way.
;
IF (n_elements(NumIms) NE 0) THEN nfiles = nfiles < NumIms
mike_log, string(nfiles)+' files to process', Verb=1

; Loop over the files
fOR F = 0,(nfiles-1) DO BEGIN

   Filename = Filelist[F]
   Filenum = string(F+1,strtrim(nfiles,2),format='(I5,"/",A)')
   IF keyword_set(verbose) THEN print, '% RA_MIKE_PIPELINE: Now processing file ' + $
     filename + ' ' + filenum

   IF strpos(strupcase(filename), 'SCI') GE 0 THEN BEGIN
      mike_log, err = 2, 'File has already been processed!'
      GOTO, error
   ENDIF

   outfile = repstr(filename, '_eng', '_sci')
   
   ; Test for existence of file
   exists = file_test(filename)
   IF exists EQ 0 THEN BEGIN
      mike_log, 'File ' + filename + ' not found.', err = 2, verbose = verbose
      GOTO, error
   ENDIF

; Make sure we are reading a valid FITS file
   A = bytarr(10)
   openr, Lun, Filename, /get_lun
   readu, Lun, A
   free_lun, Lun
   
   IF (string(A) NE 'SIMPLE  = ') THEN BEGIN
      mike_log, err = 2, 'The file: ' + filename + ' is not a valid FITS file'
      GOTO, error
   ENDIF

; Read in the data
   rdfits_struct, filename, data, /silent
   primary_header = data.hdr0
   sxaddpar, primary_header, 'BITPIX', -32, /savecomment

; Is the aperture door open?If not, it doesn't make much sense to
; apply corrections for dark current, effective area and flatfield.
   IF strupcase(sxpar(primary_header, 'APDOOR')) EQ 'OPEN' THEN mikeflags.apdoor = 1 ELSE $
     mikeflags.apdoor = 0

; Right now skip data that have been binned or windowed. If this
; becomes a popular thing to do, this will have to be revised.
   IF sxpar(primary_header, 'WILOSPEC') NE 0 OR $
     (sxpar(primary_header, 'WIHISPEC') NE 0 AND sxpar(primary_header, 'WIHISPEC') NE 1023) OR $
     sxpar(primary_header, 'WILOSPAT') NE 0 OR $
     (sxpar(primary_header, 'WIHISPAT') NE 0 AND sxpar(primary_header, 'WIHISPAT') NE 31) OR $
     (sxpar(primary_header, 'WICOSPEC') NE 0 AND sxpar(primary_header, 'WICOSPEC') NE 1) OR $
     (sxpar(primary_header, 'WICOSPAT') NE 0 AND sxpar(primary_header, 'WICOSPAT') NE 1) $
     THEN BEGIN
      mike_log, err = 2, 'Windowed and/or binned data currently not allowed'
      GOTO, error
   ENDIF

   mikeflags.mode = strupcase(sxpar(primary_header, 'ACQMODE'))

; Initialize the MIKE header block
   mike_header_init, filename, primary_header
   sxaddpar, primary_header, 'DARKFILE', /savecomment, $
             strmid(darkfile,strpos(darkfile,'/',/reverse_search)+1)
   sxaddpar, primary_header, 'FLATFILE', /savecomment, $
             strmid(flatfile,strpos(flatfile,'/',/reverse_search)+1)
   sxaddpar, primary_header, 'AEFFFILE', /savecomment, $
             strmid(aefffile,strpos(aefffile,'/',/reverse_search)+1)

; Add the geometry keywords to the header
   ra_mike_add_geometry, primary_header

   exptime = float(sxpar(primary_header, 'EXPTIME'))

; Split the pipeline here based on the acquisition mode of the
; observation.
   CASE mikeflags.mode OF 

;=====================================
; HISTOGRAM observations
;=====================================
      'HISTOGRAM': BEGIN
         image = data.im0

; The format of the images as it comes out of LIMA contains the PHD in
; the first 16 elemnts of the image array. This is unneccesary since
; this information is also contained in the PHD extension, and it also
; messes up automatic image scaling routines, since the values of PHD
; bins will be significantly greater than actual image values,
; reducing the available dynamic range. So, set the PHD portion of the
; image to 0. 
         image[0:15, 0] = 0

; If the HISTOGRAM image has saturated, then replace these pixel
; values with NaNs (Not A Numbers).
         saturated = where(image EQ 65535, nsaturated)
         IF nsaturated GT 0 THEN image[saturated] = !values.f_nan

; Create the error image
         error_image = sqrt(data.im0)
         mkhdr, error_header, error_image, /image
         sxaddpar, error_header, 'EXTNAME', 'Errors'

; Divide by the exposure time
         image /= exptime
         error_image /= exptime
         sxaddpar, primary_header, 'BUNIT', 'counts s**-1', /savecomment
         sxaddpar, error_header, 'BUNIT', 'counts s**-1', 'Units of data in the error extension'

; Create the wavelength image
         IF sxpar(primary_header, 'DETSTIM') EQ 'stim' AND $
           sxpar(primary_header, 'EVENTS') GE 100 THEN BEGIN
; Find the average location of the stim pulse.
            stim = float(image[1:40, 22])
            stimindex = findgen(40) + 1
            stimpos = total(stim * stimindex)/total(stim)
            wave_image = ralice_wavecal(stimpos = stimpos, row_offset = row_offset, $
                           hdr = waveheader)
         ENDIF ELSE wave_image = ralice_wavecal(t = sxpar(primary_header, 'T_DELECC'), $
           row_offset = row_offset, hdr = waveheader)

; Dead time correction
         IF keyword_set(deadflag) THEN BEGIN
            deadtime_correction = mike_deadtime(sxpar(primary_header, 'EVENTS')/exptime)
            image *= deadtime_correction
            error_image *= deadtime_correction
            sxaddpar, primary_header, 'DEADFLAG', 'T', ' Dead-time correction applied?'
         ENDIF

; Correction for dark counts. 
         IF keyword_set(darkflag)  AND mikeflags.apdoor THEN BEGIN
            mike_log, '   dark subtraction = YES'
            image = adderr(image, error_image, dark, darkerror, err = error_image, /subtract)
            sxaddpar, primary_header, 'DARKFLAG', 'T', ' Dark count correction applied?'
         ENDIF

; Flatfield Correction
         IF keyword_set(flatflag)  AND mikeflags.apdoor THEN BEGIN
            mike_log, '   flatfield correction = YES'
            image /= flat 
            error_image /= flat 
            sxaddpar, primary_header, 'FLATFLAG', 'T', ' Flatfield correction applied?'
         ENDIF

; Effective Area Calibration         
         IF keyword_set(aeffflag) AND mikeflags.apdoor THEN BEGIN
            aeff_image = interpol(aeff, aeff_wave, wave_image)
            image /= aeff_image
            error_image /= aeff_image 
            sxaddpar, primary_header, 'AEFFFLAG', 'T', ' Effective area correction applied?'
            sxaddpar, primary_header, 'BUNIT', 'photons cm**-2 s**-1'
            sxaddpar, error_header, 'BUNIT', 'photons cm**-2 s**-1'
         ENDIF

; Now write the FITS file and its extensions
         sxaddpar, primary_header, 'BZERO', 0, ' Data are 32 bit floating point numbers'
         writefits, OutFile, Image, primary_Header
         writefits, outfile, error_image, error_header, /app

         mkhdr, Hwave, wave_image, /image
         sxaddpar, Hwave, 'XTENSION', 'IMAGE', 'Image extension: Wavelength Image'
         sxaddpar, Hwave, '',''
         sxaddpar, Hwave, 'COMMENT', 'Alice Wavelength Calibration 2-D Image'
         sxaddpar, Hwave, 'EXTNAME','Wavelengths'
         sxaddpar, Hwave, 'BUNIT', 'Angstroms', 'Units for the wavelength array'
         writefits, OutFile, Wave_Image, Hwave, /app

         HDRtmp = data.hdr1
         sxaddpar, HDRtmp, 'EXTNAME','Pulse Height Distribution'
         writefits, OutFile, data.im1, HDRtmp, /app
         HDRtmp = data.hdr2
         sxaddpar, HDRtmp, 'EXTNAME','Count Rate'
         writefits, OutFile, data.im2, HDRtmp, /app

      END

;=====================================
      'PIXELLIST': BEGIN
;=====================================
         image = data.im0
         
; Create the error image
         error_image = sqrt(image)
         mkhdr, error_header, error_image, /image
         sxaddpar, error_header, 'EXTNAME', 'Errors'

; Create the wavelength image   
         IF sxpar(primary_header, 'DETSTIM') EQ 'stim' AND $
           sxpar(primary_header, 'EVENTS') GE 100 THEN BEGIN
; Find the average location of the stim pulse.
            stim = float(image[1:40, 22])
            stimindex = findgen(40) + 1
            stimpos = total(stim * stimindex)/total(stim)
            wave_image = ralice_wavecal(stimpos = stimpos, row_offset = row_offset, $
                           hdr = waveheader)
         ENDIF ELSE wave_image = ralice_wavecal(t = sxpar(primary_header, 'T_DELECC'), $
           row_offset = row_offset, hdr = waveheader)

; Call the routine deform_pixellist to get information about each event
         pixtable = pixellist_table(data.im1, wave_image = wave_image)

; Add a column containing the ephemeris time of the photon event to
; the pixel list table. 
         time = pixtable[3,*] * sxpar(data.hdr0, 'TIMEHCKC')
         CSPICE_UTC2ET, sxpar(data.hdr0, 'STRTSCET'), starttime
         time += starttime
         pixtable = [pixtable, time]

; Dead time correction
         IF keyword_set(deadflag) THEN BEGIN
            orig_image = image
            image = fltarr(1024, 32)

            cntrate = histogram(pixtable[3,*], binsize = 1, reverse_indices = R)
            hacktime = 2.^(-8 + sxpar(primary_Header,'TIMEHCKR'))
            deadtime_correction = mike_deadtime(cntrate/hacktime)

; Use the arcane and powerful reverse indices to find the number of
; observed photon events between consecutive timehacks. 
            FOR i = 0l, n_elements(cntrate)-1 DO BEGIN
               IF r[i] NE r[i+1] THEN BEGIN
                  inbin = r[r[i]:r[i+1]-1]
; There is probably a better way to do this without using a for
; loop,but I can't think of one right now, so just plow on ahead with
; brutre force...
                  FOR j = 0, n_elements(inbin)-1 DO $
                          image[pixtable[0,inbin[j]], pixtable[1,inbin[j]]] += $
                            deadtime_correction[i]
               ENDIF
            ENDFOR

            notzero = where(orig_image NE 0, count)
            IF count GT 0 THEN BEGIN
               factor = image[notzero]/orig_image[notzero]
               error_image[notzero] *= factor
            ENDIF
            
            sxaddpar, primary_header, 'DEADFLAG', 'T', ' Dead-time correction applied?'
         ENDIF

; Divide by the exposure time
         image /= exptime
         error_image /= exptime
         sxaddpar, primary_header, 'BUNIT', 'counts s**-1'
         sxaddpar, error_header, 'BUNIT', 'counts s**-1', 'Units of data in the error extension'

; Correction for dark counts. 
         IF keyword_set(darkflag)  AND mikeflags.apdoor THEN BEGIN
            mike_log, '   dark subtraction = YES'
            image = adderr(image, error_image, dark, darkerror, err = error_image, /subtract)
            sxaddpar, primary_header, 'DARKFLAG', 'T', ' Dark count correction applied?'

         ENDIF

; Flatfield Correction
         IF keyword_set(flatflag) AND mikeflags.apdoor  THEN BEGIN
            mike_log, '   flatfield correction = YES'
            image /= flat 
            error_image /= flat 
            sxaddpar, primary_header, 'FLATFLAG', 'T', ' Flatfield correction applied?'
         ENDIF

; Effective Area Calibration         
         IF keyword_set(aeffflag) AND mikeflags.apdoor THEN BEGIN
            aeff_image = interpol(aeff, aeff_wave, wave_image)
            image /= aeff_image
            error_image /= aeff_image 
            sxaddpar, primary_header, 'AEFFFLAG', 'T', ' Effective area correction applied?'
            sxaddpar, primary_header, 'BUNIT', 'photons cm**-2 s**-1', /savecomment
            sxaddpar, error_header, 'BUNIT', 'photons cm**-2 s**-1', /savecomment
         ENDIF

; Now write the FITS file and its extensions
         sxaddpar, primary_header, 'BZERO', 0, ' Data are 32 bit floating point numbers'
         writefits, OutFile, Image, primary_Header
         writefits, outfile, error_image, error_header, /app

         mkhdr, Hwave, wave_image, /image
         sxaddpar, Hwave, 'XTENSION', 'IMAGE', 'Image extension: Wavelength Image'
         sxaddpar, Hwave, '',''
         sxaddpar, Hwave, 'COMMENT', 'Alice Wavelength Calibration 2-D Image'
         sxaddpar, Hwave, 'EXTNAME','Wavelengths'
         sxaddpar, Hwave, 'BUNIT', 'Angstroms', 'Units for the wavelength array'
         writefits, OutFile, Wave_Image, Hwave, /app

         cols = 5
         rows = n_elements(pixtable)/cols
         ftcreate, cols, rows, tabhdr, tab
         sxaddpar, tabhdr, 'EXTNAME','Pixel List Event Table'
         ftput, tabhdr, tab, 'X Position', 0, string(reform(pixtable[0,*]), format = '(I4)')
         ftput, tabhdr, tab, 'Y Position', 0, string(reform(pixtable[1,*]), format = '(I3)')
         ftput, tabhdr, tab, 'Wavelength', 0, string(reform(pixtable[2,*]), format = '(F8.2)')
         ftput, tabhdr, tab, 'Timestep',   0, string(reform(pixtable[3,*]), format = '(I6)')
         ftput, tabhdr, tab, 'Time (et)',  0, string(reform(pixtable[4,*]), format = '(F16.3)')
         sxaddpar, tabhdr, 'TFORM1', 'I4      ', /savecomment
         sxaddpar, tabhdr, 'TFORM2', 'I3      ', /savecomment
         sxaddpar, tabhdr, 'TFORM3', 'F8.2    ', /savecomment
         sxaddpar, tabhdr, 'TFORM4', 'I6      ', /savecomment
         sxaddpar, tabhdr, 'TFORM5', 'I      ', /savecomment
         sxaddpar, tabhdr, 'TUNIT1', 'pixels', after = 'TFORM1', /savecomment
         sxaddpar, tabhdr, 'TUNIT2', 'pixels', after = 'TFORM2', /savecomment
         sxaddpar, tabhdr, 'TUNIT3', 'Angstroms', after = 'TFORM3', /savecomment
         sxaddpar, tabhdr, 'TUNIT4', 'time steps', after = 'TFORM4', /savecomment
         sxaddpar, tabhdr, 'TUNIT5', 'seconds past J2000', after = 'TFORM5', /savecomment

; Add a <crlf> pair to the end of each ascii table row. These are
; "stealth" characters in that they are not described by the FITS
; header as being in columns. 
         cr = bytarr(1, rows) + 13b
         lf = bytarr(1, rows) + 10b
         tab = [tab, cr, lf]
         sxaddpar, tabhdr, 'NAXIS1', n_elements(tab[*,0]) , /savecomment
         writefits, outfile, tab, tabhdr, /app

         HDRtmp = data.hdr2
         sxaddpar, HDRtmp, 'EXTNAME','Count Rate'
         writefits, OutFile, data.im2, HDRtmp, /app
      END

;=====================================
      'COUNTRATE':BEGIN
;=====================================
         cnt = data.im0

         cnterr = sqrt(cnt)
         mkhdr, error_header, cnterr, /image
         sxaddpar, error_header, 'EXTNAME','Errors'

; Divide by the exposure time
         cntinterval = sxpar(primary_header, 'crintvc')
         cnt /= cntinterval
         cnterr /= cntinterval
         sxaddpar, primary_header, 'BUNIT', 'counts s**-1', /savecomment
         sxaddpar, error_header, 'BUNIT', 'counts s**-1', 'Units of data in the error extension'

; Dead time correction
         IF keyword_set(deadflag) THEN BEGIN
            deadtime_correction = mike_deadtime(cnt)
            cnt *= deadtime_correction
            cnterr *= deadtime_correction
            sxaddpar, primary_header, 'DEADFLAG', 'T', ' Dead-time correction applied?'
         ENDIF

; Correction for dark counts. 
         IF keyword_set(darkflag)  AND mikeflags.apdoor THEN BEGIN
            totaldark = replicate(total(dark), n_elements(cnt)) 
            totaldark_error = replicate(sqrt(total(darkerror^2)), n_elements(cnt)) 
            mike_log, '   dark subtraction = YES'
            cnt = adderr(cnt, cnterr, totaldark, totaldark_error, err = cnterr, /subtract)
            sxaddpar, primary_header, 'DARKFLAG', 'T', ' Dark count correction applied?'
         ENDIF

         sxaddpar, primary_header, 'EXTEND', 'T', ' file may contain extensions', $
                   before = 'ORIGIN'
         sxaddpar, primary_header, 'BZERO', 0, ' Data are 32 bit floating point numbers'

         writefits, OutFile, cnt, primary_Header
         writefits, outfile, cnterr, error_header, /app
      END
   ENDCASE

   mike_log, 'Processed data written to file: ' + OutFile

   IF keyword_set(movefiles2sci) THEN BEGIN
      year = '20' + strmid(filename, strpos(filename, 'ra_') + 3, 2)
      mon = strmid(filename, strpos(filename, 'ra_') + 5, 2)
      outdir = data_dir + 'sci/' + year + '/' + mon + '/'
   ENDIF

   IF keyword_set(outdir) THEN file_move, outfile, outdir, /overwrite

; The actual error handling routine.If we ever get here, some error
; occurred and the calibrated file could not be created.

; STILL NEEDS WORK!
   catch, error_status
   IF error_status NE 0 THEN BEGIN
      help, calls = calls
      mike_log, err = 2, '%% Fatal Error at ' + calls
      mike_log, !error_state.msg, errnum = 2
      error: mike_log, 'File ' + filename + ' not processed!', errnum = 2
   ENDIF

endfor

dT = systime(1) - Time0
dT = strtrim(string(dT / 60., format='(F20.2)'),2)
mike_log, 'Total processing time = ' + dT + ' minutes', verb=1

IF keyword_set(logfile) THEN BEGIN
   openw, lun, logfile, /get
   printf, lun, logmessages
   free_lun, lun
ENDIF

end