;+
; 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, flatfile = flatfile, aefffile = aefffile,
; badflag = badflag, badfile = badfile, linflag = linflag, logfile = logfile,
; kernelfile = kernelfile, loadkernels = loadkernels, movefiles = movefiles
;
; 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.
; 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
;
; linflag -- Set this keyword to apply the linearization correction
; to the data.
;
; 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.
;
; movefiles -- 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
;
; LIST OF POSSIBLE DATA QUALITY FLAGS:
; 0: Data are nominal
; 1: Data are missing packets
; 2: Aperture door motion event during exposure
; 4: Stim position/Lyman alpha wavelength correction not available,
; wavelength errors of 5-10 A are possible.
;
; PROCEDURE:
;
;
;
; EXAMPLE:
;
;
;
; MODIFICATION HISTORY:
; The code has been completely rewritten to simplify the Rosetta Alice
; pipeline and significantly improve performance. After this rewrite, (v4)
; 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
; v6. 2008 Jan 10 AJS Use HK Event file to check aperture door state
; v7. 2008 Jan 23 AJS Use non-paralyzable deadtime correction
; v8. 2008 Feb 11 AJS Use proper HV-level for dark subtraction
; v9. 2008 Sep 02 AJS - Default to using post PC8 effective
; area. Correction for 3.8 kV level in Commissioning
; and Mars data. New dark count subtraction
; algorithm. Stim pixel wavelength corrections
; v10. 2008 Nov 25 AJS - Change error handling to use Poisson error bars when
; raw image counts are lower than 2500. Use Gaussian
; (sqrt) error bars when counts are greater than
; 2500. Also, run deadtime correction before calculating
; error bars.
; v11. 2009 Feb 08 AJS - Use high-order polynomial wavelength solution,
; RA_AEFF_007.TAB, hybrid model deadtime. Create
; linearized data products. Handle deadtime
; associated errors like the ideal non-paralyzable
; case. Correct Countrate files for dark counts
; v12. 2009 Feb 22 AJS - Update version flag in data filenames and fixed bug
; that prevented Aeff correction on some files.
; v13. 2009 Apr 14 AJS - Fixed bug that prevented the LIN file header from
; being updated.
; v14. 2009 Nov 10 AJS - Update to generalize how wavelength solution is
; calculated and make code more machine independent
; v15. 2009 Dec 14 AJS - Use PC10 Aeff and newly derived corrections for
; other epochs, correct count rate extension
; interval, add wavefile flag
; v16. 2009 Dec 18 AJS - Fixed bug with LIN files, wherein the derived flux
; levels were off by the ratio of the pixel size in
; the "natural" frame to the pixel size in the linear frame.
; v17. 2010 Feb 10 AJS - Update to produce "Version 3" data for all data
; through ESB3
; v18. 2010 Feb 15 AJS - Fix bug that prevented the headers for the primary
; data unit and error image extension from being
; written
; v19. 2010 Mar 26 AJS - Update ra_mike_add_geometry.pro to use in-flight
; pixel look vectors
; v20. 2010 Mar 27 AJS - Added comment to wavelength extension
;
; v21. 2010 Jul 08 AJS - Changed how effective area correction is applied.
;
; v22. 2010 Jul 13 AJS - Fixed bug that resulted in certain files being
; reported erroneously as darks. Removed correction
; for CRINTVC error that has been fixed in v17 of
; LIMA.
; v23. 2010 Jul 15 AJS - Added optional sky subtraction
pro ra_mike_pipeline, infiles, alldata = alldata, outdir = outdir, Verbose=Verbose, $
deadflag = deadflag, darkflag = darkflag, flatflag = flatflag, aeffflag = aeffflag, $
wavefile = wavefile, flatfile = flatfile, skyflag = skyflag, skyfile = skyfile, skyindex = skyindex, $
badflag = badflag, badfile = badfile, linflag = linflag, logfile = logfile, $
kernelfile = kernelfile, loadkernels = loadkernels, movefiles = movefiles, $
hvflag = hvflag, hvfile38 = hvfile38, hvfile37 = hvfile37, no_linearize = no_linearize, $
version = version
COMMON MIKE_DATA, mikeflags, logmessages
IF n_params() EQ 0 AND NOT keyword_set(alldata) THEN BEGIN
print, 'ra_mike_pipeline, infiles, alldata = alldata, outdir = outdir, Verbose=Verbose, $'
print, ' deadflag = deadflag, darkflag = darkflag, flatflag = flatflag, aeffflag = aeffflag, $'
print, ' wavefile = wavefile, flatfile = flatfile, skyflag = skyflag, skyfile = skyfile, skyindex = skyindex, $'
print, ' badflag = badflag, badfile = badfile, linflag = linflag, logfile = logfile, $'
print, ' kernelfile = kernelfile, loadkernels = loadkernels, movefiles = movefiles, $'
print, ' hvflag = hvflag, hvfile38 = hvfile38, hvfile37 = hvfile37, no_linearize = no_linearize, $'
print, ' version = version'
RETURN
ENDIF
; This is the calibration version number, not the pipeline version number
IF NOT keyword_set(version) THEN version = '3' ELSE version = strtrim(version, 2)
; 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(hvflag) EQ 0 THEN hvflag = 0
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 n_elements(skyflag) EQ 0 THEN skyflag = 0
; What machine are we on?
spawn, 'hostname', Host
Host = strsplit(host, '.', /extract)
; These directories and file references will need to be modified to
; match the directory structure on the local machine.
CASE strlowcase(Host[0]) OF
'argonath': BEGIN
calibration_dir = '/Volumes/d0/joel/projects/rosetta/ADA/data/cal/'
data_dir = '/Volumes/d0/joel/projects/rosetta/ADA/data/'
IF NOT keyword_set(kernelfile) THEN $
kernelfile = '/Volumes/d0/joel/projects/rosetta/ADA/data/spice/rosetta_kernels_argonath.txt'
END
'ariel': BEGIN
calibration_dir = '/Users/joel/projects/rosetta/ADA/data/cal/'
data_dir = '/User/joel/projects/rosetta/ADA/data/'
IF NOT keyword_set(kernelfile) THEN $
kernelfile = '/Volumes/d0/joel/projects/rosetta/ADA/data/spice/rosetta_kernels_argonath.txt'
END
'rioja': BEGIN
calibration_dir = '/Users/steffl/data/ralice/data/cal/'
data_dir = '/Users/steffl/data/ralice/data/'
IF NOT keyword_set(kernelfile) THEN $
kernelfile = '/Users/steffl/kernels/rosetta/rosetta_kernels.txt'
END
'rosetta': BEGIN
calibration_dir = '/home/soc/rdat/data/cal/'
data_dir = '/home/soc/rdat/data/'
IF NOT keyword_set(kernelfile) THEN $
kernelfile = '/home/soc/rdat/data/spice/rosetta_kernels_rosetta.txt'
END
ELSE: BEGIN
print, '% RA_MIKE_PIPELINE: Unknown host name. Unable to determine required filepaths.'
print, '% RA_MIKE_PIPELINE: Returning.'
return
END
ENDCASE
; Wavelength Calibration file
IF NOT keyword_set(wavefile) THEN wavefile = calibration_dir + $
'wavelength/RA_WAVE_006.FIT'
; 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)
; Sky Background file
IF NOT keyword_set(skyfile) THEN skyfile = calibration_dir + $
'sky/RA_SKY_3.9KV_001.FIT'
IF keyword_set(skyflag) THEN BEGIN
rdfits_struct, skyfile, sky, /silent
skyerr = float(sky.im1)
sky = float(sky.im0)
IF NOT keyword_set(skyindex) THEN skyindex = [6, 7, 8, 20, 21, 22]
ENDIF
;;;
; Set up the various keywords and flags.
;
IF n_elements(mikeflags) EQ 0 THEN $
MikeFlags = {Version:'23 [2010 Jul 15]', $
Mode:'', $
Verbose: keyword_set(verbose) ? verbose : 0 , $
Log:1, $
Error:0, $
StartTime:systime(1), $
kernel:0, $
nogeometry:0, $
apdoor:0}
Time0 = systime(1)
; 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
; 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')
filelist = filelist[where(strpos(filelist, 'ra_temp') EQ -1)]
nfiles = n_elements(filelist)
ENDIF ELSE BEGIN
nfilelist = n_elements(infiles)
nfiles = 0
filelist = ''
FOR i = 0, nfilelist-1 DO BEGIN
IF file_test(infiles[i]) THEN BEGIN
filelist_i = infiles[i]
nfiles_i = 1
ENDIF ELSE filelist_i = file_search(data_dir + 'eng/', infiles[i], count=nfiles_i)
filelist = [filelist, filelist_i]
nfiles += nfiles_i
ENDFOR
IF nfiles GT 0 THEN filelist = filelist[1:*]
ENDELSE
IF (nfiles EQ 0) THEN BEGIN
mike_log, err=2, 'No files to process.'
RETURN
ENDIF
; Find the (chronologically) first and last of the files that will be
; processed, in order to find the range of data directories spanned by
; the data. Since it is theoretically possible that a data file could
; span multiple months, add one month before and one month after this
; range.
file_basenames = file_basename(filelist)
file_basenames = file_basenames[sort(file_basenames)]
first = file_basenames[0]
last = file_basenames[nfiles-1]
first_year = fix('20' + strmid(first, 3, 2))
first_month = fix(strmid(first, 5, 2)) - 1
IF first_month LT 1 THEN BEGIN
first_year -= 1
first_month = 12
ENDIF
last_year = fix('20' + strmid(last, 3, 2))
last_month = fix(strmid(last, 5, 2)) + 1
IF last_month GT 12 THEN BEGIN
last_year += 1
last_month = 1
ENDIF
; Find all of the available HK Event files in the time span of the data
mike_log, 'Searching for valid HK Event files'
hkeventfiles = ''
nhkeventfiles = 0l
FOR year = first_year, last_year DO BEGIN
yearstr = strtrim(year, 2)
FOR month = 1, 12 DO BEGIN
IF NOT ((year EQ first_year AND month LT first_month) OR $
(year EQ last_year AND month GT last_month)) THEN BEGIN
monthstr = month LT 10 ? '0'+strtrim(month,2) : strtrim(month,2)
hkevfiles = file_search(data_dir + 'eng/' + yearstr + '/' + monthstr + '/', $
'*evnt_eng.txt', count = nhkevfiles)
IF nhkevfiles GT 0 THEN BEGIN
hkeventfiles = [hkeventfiles, hkevfiles]
nhkeventfiles += nhkevfiles
ENDIF
ENDIF
ENDFOR
ENDFOR
IF n_elements(hkeventfiles) GT 1 THEN hkeventfiles = hkeventfiles[1:*]
;Initialize variables
hkformat = "(F15.3,A24,2F15.3,I2,I6,A7,A60)"
line_scetr = 0d
line_scetc = ''
line_insttime = 0d
line_toffset = 0d
line_instsync = 0b
line_eidsubt = 0l
line_ident = ''
line_descript = ''
; Apdoor contains flags for the aperture door motion event at time,
; scetc
scetc = ''
apdoor = 0
valid_hkevtimes = [0d, 0d]
; Loop over the HK event files
FOR i = 0, nhkeventfiles-1 DO BEGIN
mike_log, 'Now reading ' + hkeventfiles[i]
openr, lun, hkeventfiles[i], /get
first_time = ''
last_time = ''
dataflag = 0
WHILE NOT eof(lun) DO BEGIN
; Read in the header info
line = ''
WHILE dataflag EQ 0 DO BEGIN
readf, lun, line
IF line EQ 'START DATA' THEN dataflag = 1
ENDWHILE
readf, lun, line_Scetr, line_ScetC, line_InstTime, line_TOffset, line_InstSync, $
line_EidSubt, line_Ident, line_Descript, format = hkformat
IF first_time EQ '' THEN first_time = line_scetc
last_time = line_scetc
; Search for Aperture Door Motion events. If found, determine what
; kind of motion event it is.
IF line_ident EQ ' ADMOT:' THEN BEGIN
scetc = [scetc, line_scetc]
CASE 1 OF
strpos(line_descript, 'closed') GE 0 : line_apdoor = 0
strpos(line_descript, 'open') GE 0 : line_apdoor = 1
strpos(line_descript, 'between') GE 0 : line_apdoor = 2
strpos(line_descript, 'error') GE 0 : line_apdoor = 3
ENDCASE
apdoor = [apdoor, line_apdoor]
ENDIF
ENDWHILE
; Get the first and last times covered by this HK event file
CSPICE_UTC2ET, first_time, first_et
CSPICE_UTC2ET, last_time, last_et
valid_hkevtimes = [[valid_hkevtimes], [first_et, last_et]]
free_lun, lun
ENDFOR
IF nhkeventfiles GT 0 THEN BEGIN
scetc = scetc[1:*]
apdoor = apdoor[1:*]
valid_hkevtimes = valid_hkevtimes[*, 1:*]
ENDIF
et_hkev = dblarr(n_elements(scetc))
FOR i = 0, n_elements(scetc)-1 do BEGIN
CSPICE_UTC2ET, scetc[i], et
et_hkev[i] = et
ENDFOR
; Create the 2D array of row offsets in pixel units and initialize the
; wavelength solution coefficients.
rdfits_struct, wavefile, wave, /silent
wavesoln = wave.im0
poly2 = wave.im1
poly1 = wave.im2
row_offset = wave.im3
row_offset = (fltarr(1024) + 1.) # row_offset
wave_comments = [$
'Note that the actual wavelength range covered by the Alice UV ', $
'spectrometer is 700-2050 Angstroms. When a charge pulse strikes the ', $
'readout anodes, the detector electronics assigns it a 10-bit X value', $
'and a 5-bit Y-value. The mapping of physical space on the detector to', $
'data space does not use the full 10-bit X range. Instead, events on the', $
'detector are mapped into columns 130-900 (approximately) in data space.', $
'The regions in data space on either side of this range do not have any', $
'physical meaning, and are conceptually similar to overscan areas found', $
'in CCD data.']
; Define the nominal linear wavelength scale
xindex = findgen(1024)
lin_wavelen = float(poly(xindex, poly1))
lin_wave_image = lin_wavelen # (fltarr(32)+1.)
yindex = (fltarr(1024) + 1.) # findgen(32)
; Call poisson_errorbars.pro to get vector of Poisson errors
poisson_errors = dblarr(2501)
FOR i = 0, 2500 DO poisson_errors[i] = poisson_errorbar(i, /double)
poisson_errors = FLOAT(poisson_errors)
;;;
; Loop over the files
IF (n_elements(NumIms) NE 0) THEN nfiles = nfiles < NumIms
mike_log, string(nfiles)+' files to process', Verb=1
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 ' + $
file_basename(filename) + ' ' + filenum
IF strpos(strupcase(filename), 'SCI') GE 0 THEN BEGIN
mike_log, err = 2, 'File has already been processed!'
GOTO, error
ENDIF
oldversion = strmid(filename, 8, 1, /reverse_offset)
outfile = repstr(file_basename(filename), oldversion + '_eng', version + '_sci')
linfile = repstr(file_basename(filename), oldversion + '_eng', version + '_lin')
; 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 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
exptime = float(sxpar(primary_header, 'EXPTIME'))
hvlevel = sxpar(primary_header, 'HVLEVELR')
; Initialize the MIKE header block
ra_mike_header_init, filename, outfile, primary_header
; Add the wavefile
sxaddpar, primary_header, 'WAVEFILE', file_basename(wavefile), /savecomment
; Add the geometry keywords to the header
ra_mike_add_geometry, primary_header
; Are there are missing packets in this file?
lostpackets = sxpar(primary_header, 'LOSTPKTS')
IF lostpackets GT 0 THEN BEGIN
data_quality = sxpar(primary_header, 'DATAQUAL')
data_quality += 1
sxaddpar, primary_header, 'DATAQUAL', data_quality, /savecomment
ENDIF
; Is the aperture door open? If the door is closed for the whole time,
; the image is a dark, and it's not appropriate to apply dark current,
; effective area and flatfield corrections. Attempt to determine the
; aperture door state from the HK event files.
CSPICE_UTC2ET, sxpar(primary_header, 'STRTSCET'), et_start
CSPICE_UTC2ET, sxpar(primary_header, 'STOPSCET'), et_stop
; Initialize the open_time and closed_time variables
open_time = 0.
closed_time = 0.
hkevind = where(valid_hkevtimes[0,*] LE et_start AND valid_hkevtimes[1,*] GE et_start AND $
valid_hkevtimes[0,*] LE et_stop AND valid_hkevtimes[1,*] GE et_stop, $
hkevfound)
IF hkevfound EQ 1 THEN BEGIN
; The start and stop times of this exposure fall into the time span
; covered by one of the HK event files. Thus, the first aperture door
; motion event prior to the start time should give the initial state
; of the aperture door.
hkev_index = max(where(et_hkev-et_start LT 0))
; The exposure happened before the first door motion event, so the
; door state cannot be determined from the event file.
IF hkev_index EQ -1 THEN GOTO, nodoorinfo
mikeflags.apdoor = apdoor[hkev_index]
; If the door changed during the exposure, or if the error condition
; was set, set the data quality flag and deal with it appropriately.
events = where(et_hkev GE et_start AND et_hkev LT et_stop, nevents)
IF nevents GT 0 THEN BEGIN
apdoor_events = [mikeflags.apdoor, apdoor[events], apdoor[events[nevents-1]]]
event_et = [et_start, et_hkev[events], et_stop]
; Although the door open/close times are of the order of 55 msec., the
; door event sampling rate is usually 1 Hz, so assume if an event
; comes through, that the door position instantly changed at the event
; time.
IF mikeflags.apdoor EQ 1 THEN open = 1 ELSE open = 0
FOR i = 0, n_elements(apdoor_events)-2 DO BEGIN
IF open EQ 1 THEN BEGIN
open_time += event_et[i+1] - event_et[i]
IF apdoor_events[i+1] EQ 0 OR apdoor_events[i+1] GT 1 THEN open = 0
ENDIF ELSE BEGIN
closed_time += event_et[i+1] - event_et[i]
IF apdoor_events[i+1] EQ 1 OR apdoor_events[i+1] EQ 2 THEN open = 1
ENDELSE
ENDFOR
open_time = float(open_time)
closed_time = float(closed_time)
; If there was an aperture door motion error event, all bets are off,
; so don't do any more processing.
IF max(apdoor_events GT 2) THEN mikeflags.apdoor = 3
; Write warnings to the log file and adjust the data quality flag
mike_log, err = 1, strtrim(nevents, 2) + $
' aperture door motion event during exposure!'
data_quality = sxpar(primary_header, 'DATAQUAL')
data_quality += 2
sxaddpar, primary_header, 'DATAQUAL', data_quality, /savecomment
; We have HK event info and there was no aperture door motion event
; during the exposure.
ENDIF ELSE CASE mikeflags.apdoor OF
0: closed_time = exptime
1: open_time = exptime
ELSE:
ENDCASE
sxaddpar, primary_header, 'APDRSRC', 'HK event file', $
' Source of aperture door status', after = 'DATAQUAL'
ENDIF ELSE BEGIN
; We don't have any reliable information from the HK event file, so
; just trust the value in the data header.
nodoorinfo:apdoor_state = strupcase(sxpar(primary_header, 'APDOOR'))
CASE 1 OF
apdoor_state EQ 'CLOSED': BEGIN
mikeflags.apdoor = 0
closed_time = exptime
END
apdoor_state EQ 'OPEN': BEGIN
mikeflags.apdoor = 1
open_time = exptime
END
apdoor_state EQ 'BETWEEN': mikeflags.apdoor = 2
apdoor_state EQ 'ERROR': mikeflags.apdoor = 3
ENDCASE
sxaddpar, primary_header, 'APDRSRC', 'Data header', $
' Source of aperture door status', after = 'DATAQUAL'
ENDELSE
; Add the OPENTIME and CLOSTIME keywords to the primary header
sxaddpar, primary_header, 'OPENTIME', open_time, ' Time aperture door was open', $
after = 'EXPTIME'
sxaddpar, primary_header, 'CLOSTIME', closed_time, ' Time aperture door was closed', $
after = 'OPENTIME'
; If the opened time is less than the total exptime, treat as a science image,
; if the open time is equal to the exptime, it is, by definition, a dark
IF closed_time EQ exptime THEN isdark = 1 ELSE isdark = 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'))
; High voltage set point
hvsetc = sxpar(primary_header, 'hvsetc')
; Split the pipeline here based on the acquisition mode of the
; observation.
CASE mikeflags.mode OF
;=====================================
; HISTOGRAM observations
;=====================================
'HISTOGRAM': BEGIN
image = data.im0
; Create error image. Use the Gaussian approximation where image > 2500
; counts, otherwise use Poisson errors.
error_image = FLTARR(SIZE(image, /DIM))
over2500 = WHERE(image GT 2500, nover2500, COMPLEMENT = under2500, NCOMPLEMENT = nunder2500)
IF nover2500 GT 0 THEN error_image[over2500] = SQRT(image[over2500])
IF nunder2500 GT 0 THEN error_image[under2500] = poisson_errors[image[under2500]]
mkhdr, error_header, error_image, /image
sxaddpar, error_header, 'EXTNAME', 'Errors'
; Create the wavelength image
stim_correction = ra_stim_lya_correction(image, primary_header, stimflag = stimflag, $
lyaflag = lyaflag, /inverse)
IF n_elements(stim_correction) EQ 1024 THEN BEGIN
wavelength_indices = stim_correction
wavelength_indices = rebin(wavelength_indices, 1024, 32)
mike_log, ' Stim pixel correction = YES'
sxaddpar, primary_header, 'WAVEFLAG', 'T', /savecomment
ENDIF ELSE wavelength_indices = rebin(findgen(1024), 1024, 32)
wavelength_indices_no_offset = wavelength_indices
wavelength_indices -= row_offset
wave_image = ra_wavelen(wavesoln, poly2, index = wavelength_indices, wavefile = wavefile)
; 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
; Dead time correction
IF keyword_set(deadflag) THEN BEGIN
countrate = sxpar(primary_header, 'EVENTS')/exptime
; If the (digital) countrate is > 50 kHz, the file is probbaly a test pattern
; or something really weird is going on. Either way, we probbaly shouldn't
; process it.
IF countrate LT 5e4 THEN BEGIN
deadtime_correction = ra_mike_deadtime(countrate)
image *= deadtime_correction
error_image *= deadtime_correction
sxaddpar, primary_header, 'DEADFLAG', 'T', ' Dead-time correction applied?'
ENDIF ELSE BEGIN
mike_log, err = 2, 'File has unrealistically high countrate and will not be processed.'
mike_log, err = 2, 'This may be because the file is a test pattern and not real science data.'
GOTO, error
ENDELSE
ENDIF
; 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
; Correction for dark counts.
IF keyword_set(darkflag) AND isdark EQ 0 THEN BEGIN
ra_darksubtract, calibration_dir, primary_header, image, error_image, wavelength_indices_no_offset
mike_log, ' dark subtraction = YES'
ENDIF
; Subtract off sky background.
IF keyword_set(skyflag) AND isdark EQ 0 THEN BEGIN
ra_mike_skysubtract, primary_header, image, error_image, sky, skyerr, skyfile, $
wavelength_indices_no_offset, skyindex = skyindex
mike_log, ' sky subtraction = YES'
ENDIF
; Divide by the time the aperture door was open, or the exposure time
; if this is a dark.
IF isdark EQ 0 THEN BEGIN
image /= open_time
error_image /= open_time
ENDIF ELSE BEGIN
image /= exptime
error_image /= exptime
ENDELSE
sxaddpar, primary_header, 'BUNIT', 'counts s**-1', /savecomment
sxaddpar, error_header, 'BUNIT', 'counts s**-1', $
'Units of data in the error extension'
; Flatfield Correction
IF keyword_set(flatflag) AND isdark EQ 0 THEN BEGIN
mike_log, ' flatfield correction = YES'
image /= flat
error_image /= flat
sxaddpar, primary_header, 'FLATFLAG', 'T', ' Flatfield correction applied?'
sxaddpar, primary_header, 'FLATFILE', /savecomment, file_basename(flatfile)
ENDIF
; Effective Area Calibration
IF keyword_set(aeffflag) AND isdark EQ 0 THEN BEGIN
aeff_image = ra_aeff_norm(primary_header, wave_image, calibration_dir = calibration_dir)
image /= aeff_image
error_image /= aeff_image
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'
IF n_elements(stimcorrection) EQ 1024 THEN sxaddpar, Hwave, 'COMMENT', $
'Stim pixel wavelength correction: APPLIED' ELSE BEGIN
sxaddpar, hwave, 'COMMENT', 'Stim pixel wavelength correction: NOT APPLIED'
sxaddpar, hwave, 'COMMENT', 'wavelength solution may have innaccuracies of 5-10 Angstoms'
ENDELSE
sxaddpar, Hwave, 'EXTNAME','Wavelengths'
sxaddpar, Hwave, 'BUNIT', 'Angstroms', 'Units for the wavelength array'
sxaddhist, wave_comments, Hwave, /comment
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, data.hdr2, 'EXTNAME','Count Rate'
writefits, OutFile, data.im2, HDRtmp, /app
; Now create the linearized image and write it and its extensions to a FITS file
IF NOT keyword_set(no_linearize) THEN BEGIN
; Divide by the natural wavelength scale, such that the units are:
; photons / s / cm^2 / Angstrom.
; This prevents the integrated flux level from changing when you interpolate.
dwave_image = wave_image - shift(wave_image,-1,0)
image_lin = image / dwave_image
error_image_lin = error_image / dwave_image
lin_wavelength_indices = fltarr(1024, 32)
FOR jj = 0, 31 DO lin_wavelength_indices[*, jj] = interpol(xindex, wave_image[*, jj], lin_wavelen)
lin_wavelength_indices_no_offset = rebin(lin_wavelength_indices[*, 15], 1024, 32)
lin_normalization = total(image_lin[80:950, *])
; Interpolate the image onto the new Linear scale
image_lin = interpolate(image/dwave_image, lin_wavelength_indices, yindex, cubic = -0.5)
error_image_lin = interpolate(error_image/dwave_image, lin_wavelength_indices, yindex, cubic = -0.5)
; Normalize so that we haven't lost any flux due to interpolation errors
lin_normalization /= total(image_lin[80:950, *])
mike_log, ' Image Linearization = YES'
sxaddpar, primary_header, 'FILE_OUT', strupcase(linfile), /savecomment
sxaddpar, primary_header, 'LINFLAG', 'T', /savecomment
sxaddpar, primary_header, 'BUNIT', 'photons cm^-2 s^-1 Angstrom^-1'
sxaddpar, error_header, 'BUNIT', 'photons cm^-2 s^-1 Angstrom^-1'
writefits, linFile, Image_lin, primary_header
writefits, linfile, error_image_lin, error_header, /app
mkhdr, Hwave, lin_wavelen, /image
sxaddpar, Hwave, 'XTENSION', 'IMAGE', 'Image extension: Wavelength Image'
sxaddpar, Hwave, '',''
sxaddpar, Hwave, 'COMMENT', 'R-Alice liearized wavelength solution'
IF n_elements(stimcorrection) EQ 1024 THEN sxaddpar, Hwave, 'COMMENT', $
'Stim pixel wavelength correction: APPLIED' ELSE BEGIN
sxaddpar, hwave, 'COMMENT', 'Stim pixel wavelength correction: NOT APPLIED'
sxaddpar, hwave, 'COMMENT', 'wavelength solution may have innaccuracies of 5-10 Angstoms'
ENDELSE
sxaddpar, Hwave, 'EXTNAME','Wavelength'
sxaddpar, Hwave, 'BUNIT', 'Angstroms', 'Units for the wavelength array'
sxaddhist, wave_comments, Hwave, /comment
writefits, linFile, lin_wavelen, Hwave, /app
HDRtmp = data.hdr1
sxaddpar, HDRtmp, 'EXTNAME','Pulse Height Distribution'
writefits, linfile, data.im1, HDRtmp, /app
HDRtmp = data.hdr2
sxaddpar, HDRtmp, 'EXTNAME','Count Rate'
writefits, linfile, data.im2, HDRtmp, /app
ENDIF
END
;=====================================
'PIXELLIST': BEGIN
;=====================================
image = data.im0
; Create error image. Use the Gaussian approximation where image > 2500
; counts, otherwise use Poisson errors.
error_image = FLTARR(SIZE(image, /DIM))
over2500 = WHERE(image GT 2500, nover2500, COMPLEMENT = under2500, NCOMPLEMENT = nunder2500)
IF nover2500 GT 0 THEN error_image[over2500] = SQRT(image[over2500])
IF nunder2500 GT 0 THEN error_image[under2500] = poisson_errors[image[under2500]]
mkhdr, error_header, error_image, /image
sxaddpar, error_header, 'EXTNAME', 'Errors'
; Create the wavelength image
stim_correction = ra_stim_lya_correction(image, primary_header, stimflag = stimflag, $
lyaflag = lyaflag, /inverse)
IF n_elements(stim_correction) EQ 1024 THEN BEGIN
wavelength_indices = stim_correction
wavelength_indices = rebin(wavelength_indices, 1024, 32)
mike_log, ' Stim pixel correction = YES'
sxaddpar, primary_header, 'WAVEFLAG', 'T', ' Stim pixel correction applied to wavelength solution?'
ENDIF ELSE wavelength_indices = rebin(findgen(1024), 1024, 32)
wavelength_indices_no_offset = wavelength_indices
wavelength_indices -= row_offset
wave_image = ra_wavelen(wavesoln, poly2, index = wavelength_indices, wavefile = wavefile)
; 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 = ra_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 factor = image[notzero]/orig_image[notzero]
sxaddpar, primary_header, 'DEADFLAG', 'T', ' Dead-time correction applied?'
ENDIF
; Correction for dark counts.
IF keyword_set(darkflag) AND mikeflags.apdoor THEN BEGIN
ra_darksubtract, calibration_dir, primary_header, image, error_image, wavelength_indices_no_offset
mike_log, ' dark subtraction = YES'
ENDIF
; Subtract off sky background.
IF keyword_set(skyflag) AND isdark EQ 0 THEN BEGIN
ra_mike_skysubtract, primary_header, image, error_image, sky, skyerr, skyfile, $
wavelength_indices_no_offset, skyindex = skyindex
mike_log, ' sky subtraction = YES'
ENDIF
; Divide by the time the aperture door was open, or the exposure time
; if this is a dark.
IF isdark EQ 0 THEN BEGIN
image /= open_time
error_image /= open_time
ENDIF ELSE BEGIN
image /= exptime
error_image /= exptime
ENDELSE
sxaddpar, primary_header, 'BUNIT', 'counts s**-1'
sxaddpar, error_header, 'BUNIT', 'counts s**-1', $
'Units of data in the error extension'
; 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?'
sxaddpar, primary_header, 'FLATFILE', /savecomment, file_basename(flatfile)
ENDIF
; Effective Area Calibration
IF keyword_set(aeffflag) AND isdark EQ 0 THEN BEGIN
aeff_image = ra_aeff_norm(primary_header, wave_image, calibration_dir = calibration_dir)
image /= aeff_image
error_image /= aeff_image
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'
IF n_elements(stimcorrection) EQ 1024 THEN sxaddpar, Hwave, 'COMMENT', $
'Stim pixel wavelength correction: APPLIED' ELSE BEGIN
sxaddpar, hwave, 'COMMENT', 'Stim pixel wavelength correction: NOT APPLIED'
sxaddpar, hwave, 'COMMENT', 'wavelength solution may have innaccuracies of 5-10 Angstoms'
ENDELSE
sxaddpar, Hwave, 'EXTNAME','Wavelengths'
sxaddpar, Hwave, 'BUNIT', 'Angstroms', 'Units for the wavelength array'
sxaddhist, wave_comments, Hwave, /comment
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
; Now create the linearized image and write it and its extensions to a FITS file
IF NOT keyword_set(no_linearize) THEN BEGIN
; Divide by the natural wavelength scale, such that the units are counts /
; angstrom. This prevents the integrated flux level from changing when you interpolate.
dwave_image = wave_image - shift(wave_image,-1,0)
image_lin = image / dwave_image
error_image_lin = error_image / dwave_image
lin_wavelength_indices = fltarr(1024, 32)
FOR jj = 0, 31 DO lin_wavelength_indices[*, jj] = interpol(xindex, wave_image[*, jj], lin_wavelen)
lin_wavelength_indices_no_offset = rebin(lin_wavelength_indices[*, 15], 1024, 32)
lin_normalization = total(image_lin[80:950, *])
image_lin = interpolate(image/dwave_image, lin_wavelength_indices, yindex, cubic = -0.5)
error_image_lin = interpolate(error_image/dwave_image, lin_wavelength_indices, yindex, cubic = -0.5)
; Normalize so that we haven't lost any flux due to interpolation errors
lin_normalization /= total(image_lin[80:950, *])
mike_log, ' Image Linearization = YES'
sxaddpar, primary_header, 'FILE_OUT', strupcase(linfile), /savecomment
sxaddpar, primary_header, 'LINFLAG', 'T', /savecomment
sxaddpar, primary_header, 'BUNIT', 'photons cm^-2 s^-1 Angstrom^-1'
sxaddpar, error_header, 'BUNIT', 'photons cm^-2 s^-1 Angstrom^-1'
writefits, linfile, Image_lin, primary_header
writefits, linfile, error_image_lin, error_header, /app
mkhdr, Hwave, lin_wavelen, /image
sxaddpar, Hwave, 'XTENSION', 'IMAGE', 'Image extension: Wavelength Image'
sxaddpar, Hwave, '',''
sxaddpar, Hwave, 'COMMENT', 'R-Alice liearized wavelength solution'
IF n_elements(stimcorrection) EQ 1024 THEN sxaddpar, Hwave, 'COMMENT', $
'Stim pixel wavelength correction: APPLIED' ELSE BEGIN
sxaddpar, hwave, 'COMMENT', 'Stim pixel wavelength correction: NOT APPLIED'
sxaddpar, hwave, 'COMMENT', 'wavelength solution may have innaccuracies of 5-10 Angstoms'
ENDELSE
sxaddpar, Hwave, 'EXTNAME','Wavelength'
sxaddpar, Hwave, 'BUNIT', 'Angstroms', 'Units for the wavelength array'
sxaddhist, wave_comments, Hwave, /comment
writefits, linfile, lin_wavelen, Hwave, /app
writefits, linfile, tab, tabhdr, /app
HDRtmp = data.hdr2
sxaddpar, HDRtmp, 'EXTNAME','Count Rate'
writefits, linfile, data.im2, HDRtmp, /app
ENDIF
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 = ra_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 isdark EQ 0 THEN BEGIN
ra_darksubtract, calibration_dir, primary_header, cnt
mike_log, ' dark subtraction = YES'
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, /append
END
ENDCASE
mike_log, 'Processed data written to file: ' + OutFile
IF keyword_set(movefiles) 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 + '/'
outdir_lin = repstr(outdir, 'sci', 'lin')
ENDIF
IF keyword_set(outdir) THEN BEGIN
file_move, outfile, outdir, /overwrite
IF NOT keyword_set(no_linearize) AND mikeflags.mode NE 'COUNTRATE' THEN file_move, linfile, outdir_lin, /overwrite
ENDIF
; 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
; message, /reset_error_state
; catch, /cancel
; ENDIF
error:foo = 0
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