; NAME:
;   pro darkcurrent_reduce.pro (v2.16)
;
; PURPOSE:
;   This program takes a images listed in filename and produces estimates of
;   dark current from them.  A per-pixel dark current frame is produced, and
;   individual dark current slope plots can optionally be created.
;
; CALLING SEQUENCE:
;   darkcurrent_reduce, testdir, [outdir [, infile [, valoutfile $
;       [, noref [, ma_refcor [, region [,imgoutfile] [, groupfactor]]]]]]
;
; INPUTS:
;   testdir             string containing directory name wherein input
;                           images and a list of the image file names
;                           can be found
;
; KEYWORD PARAMETERS:
;   outdir              directory where output will be written; defaults to
;                           testdir + '\SlopePlots'
;   infile              input file name.  File must be in directory specified by
;                           testdir.  Default name is 'dark.lst'
;   valoutfile          name for output file containing read-by-read mode values.
;                           Will be written to outdir.  Default is
;                           'slopevals.txt'.
;   imgoutfile          name for output file containing names of slope images.
;                           Will be written to outdir.  Default is
;                           'slopeimgs.txt'.
;   region              rectangular region of SCA to consider for
;                           analysis.  Should be entered as 4-element
;                           1-d integer array, i.e. [x1,x2,y1,y2] to
;                           consider pixels in rows x1 to x2 and
;                           columns y1 to y2.  Default is specified as
;                           all non-reference pixels
;   ma_refcor           set for maximal averaging reference pixel correction
;   noref               set to disable reference pixel correction
;   gdmask              set equal to the name of a good pixel mask.
;                           1=good & 0=bad
;   groupfactor         sets number of reads placed together in a group for
;                           processing; defaults to 1 (which means
;                           mode of each read is determined).
;
;   NOTE:  Default behavior is to do spatial average reference pixel correction.
;   If both ma_refcor and noref are set, the noref keyword over-rides the
;   ma_refcor keyword and no reference pixel correction will be done.
;
; EXAMPLE
;
;   Reduce data from specified directory, using default file names,
;   generating slope plots and using maximal averaging reference pixel
;   correction, using default region.
;
;   IDL> dkdir = '\\rabbit\rockwell3\H2RG-006-5.0mu\cold1\dktest.3Apr03'
;   IDL> darkcurrent_reduce, dkdir, /ma_refcor
;
; REFERENCE:
;
; MODIFICATION HISTORY:
;   Written by:
;     Don Figer, IDTL, May 3, 2003
;   Modified by:
;     Don Figer, IDTL, May 5, 2002
;           Added output switches
;       Bernie Rauscher, IDTL, August 5, 2002
;           Added reference pixel
;           correction, glow correction, and spiffy line-fit plots.
;       Ernie Morse, IDTL, October 21, 2002 (v1.1)
;           Modified code so that choosing bjrmods option does ONLY the bjrmods
;           section, and not the rest of the procedure
;           Added code to create output file for use by darkjumbo2.pro
;           Added readtime and month2days functions for use in creating
;           darkjumbo output file
;           Increased min and max histogram bin limits--they were too small for
;           higher-signal images
;       Ernie Morse, IDTL, October 31, 2002 (v1.2)
;           Modified bjrmods section so that PostScript output files are put in
;           outdir, not the current directory.
;       Ernie Morse, IDTL, Nov. 7, 2002 (v1.2.1)
;           Now gets SCA type from header, instead of assuming H1RG
;           Calls refsub with SCA type instead of refsub_h1rg directly
;       Ernie Morse, IDTL, January 7, 2003 (v1.2.2)
;           Changed BJRMODS histogram setup so that histogram bins are
;           centered around median signal in data instead of around 0.
;           In the previous method, the actual distribution curve tended to
;           move off the edge of the histogram range as signal increased and
;           there were a lot of empty bins with negative values.
;       Ernie Morse, IDTL, February 10-12, 2003 (v1.3)
;           Changed reading of FITS files to frame-at-a-time to allow for
;           reading arbitrarily large files to prevent limitations imposed by
;           the amount of RAM available.
;           Added region keyword to allow specification of custom region.
;           Removed hard-coded assumptions that array size is 1024 x 1024
;           (array size now taken from header keywords).
;           Removed hard-coded assumptions of analysis region.
;           Changed plot parameters (increased line & char size & thickness,
;           only print det. temp. to 1 mK accuracy).
;       Ernie Morse, IDTL, March 5, 2003 (v1.3)
;           Added keywords to allow specification of type of reference pixel
;           correction
;           Added information about analysis region and reference pixel
;           correction to plots.
;       Ernie Morse, IDTL, March 13, 2003 (v1.3)
;           Added name of file to bottom margin area of slopefit plots in
;           bjrmods section.
;       Ernie Morse, IDTL, April 4, 2003 (v2.0)
;           Started all over from scratch.
;       Ernie Morse, IDTL, April 8, 2003 (v2.1)
;           Cast values read in from headers to int and float data types to
;           insure that the slope array
;           written out at the end will be a float array in order to decrease
;           size of the resulting FITS file.
;       B.J. Rauscher, IDTL, April 10, 2003 (v2.2)
;           Added ability to use good pixel masks
;       Ernie Morse, IDTL, April 11, 2003 (v2.3)
;           Changed calls to refsub--must now pass SCA name instead of SCA
;           type
;       Ernie Morse, IDTL, April 14, 2003 (v2.4)
;           Removed the bin representing 0 from the gaussfit to protect against
;           spurious zeros caused by railed pixels.
;       Ernie Morse, IDTL, April 16, 2003 (v2.5)
;           Corrected bug found in implementation of zero-bin rejection.
;       Ernie Morse, IDTL, April 30, 2003 (v2.6)
;           Added the rd_time keyword and code to allow the processing of the
;           Fowler glow images.
;       Ernie Morse, IDTL, May 9, 2003 (v2.7)
;           Eliminated sa_refcor keyword--spatial average reference pixel
;           correction is now default
;           Added noref keyword--can be set to explicitly disable reference
;           pixel correction
;           Made changes in code to reflect new default case (spatial averaging).
;           Changed default output directory to testdir+'SlopePlots'
;       Ernie Morse, IDTL, May 15, 2003 (v2.8)
;           Added write of reference pixel correction version number and type
;           to header of slope images
;       Ernie Morse, IDTL, October 10, 2003 (v2.9)
;           Added cosmic ray detection to creation of per-pixel images
;       Ernie Morse, IDTL, October 17, 2003 (v2.10)
;           Added groupfactor keyword and modifications to allow grouping
;           of reads for processing.  Split outfile keyword into valoutfile and
;           imgoutfile.  Made some changes to format of valoutfile (it now
;           must include read times as they can not be calculated in the
;           analysis script anymore due to the addition of grouping).
;       Ernie Morse, IDTL, December 2, 2003 (v2.11)
;           Corrected bugs associated with operation when groupfactor = 1.
;           Added graceful exit for case when all pixels in input image are
;           railed.
;       Ernie Morse, IDTL, July 15, 2004(v2.12)
;           Changed from using binsize = 1.0 to starting out with 128 bins with
;           size determined by data distribution.
;           Added code to catch math errors (from curvefit)
;           Corrected problems that arose when first reads of input image are
;           at bottom rail.
;       Ernie Morse, IDTL, July 22, 2004 (v2.13)
;           Added rowinc and colinc keywords.
;       Ernie Morse, IDTL, September 29, 2004 (v2.14)
;           Modified error handler so that math errors are distinguished from
;           non-math errors.
;       Ernie Morse, IDTL, October 8, 2004 (v2.15)
;           Made changes to prevent blank lines being written to imgoutfile for
;           input images that are railed in every read.  Moved opening of
;           imgoutfile to end of procedure, conditional on there being at least
;           one non-railed image and trimmed imginfo array to only include
;           entries with non-null file names.
;       Ernie Morse, IDTL, September 2, 2005 (v2.16)
;           Now uses djsiterstat to determine median and standard
;           deviation, which is used to calculate min and max
;           for histogram procedure.  Number of bins has been altered so
;           that bin size = 1.0 ADU.
;
Pro darkcurrent_reduce, testdir, outdir = outdir, infile = infile, $
        valoutfile = valoutfile, imgoutfile = imgoutfile, $
        ma_refcor = ma_refcor, noref=noref, region = region, $
        gdmask = gdmask, groupfactor = groupfactor, $
        rowinc = rowinc, colinc = colinc

    ; set directory path delimiter
    pathdelim=path_sep()

    ; set end-of-line character based on operating system in use (non-unix is assumed to be Windows).
    if (!version.os_family eq 'unix') then begin
        endofline = string(10b)
    endif else begin
        endofline = string(13b) + string(10b)
    endelse

    ; check to see if input directory name ends in path delimiter--if not,
    ; add it to the end
    if (strmid(testdir, strlen(testdir)-1, 1) ne pathdelim) then testdir = testdir + pathdelim

    ; set outdir to testdir + 'SlopePlots' if not set by user
    if (keyword_set(outdir) eq 0) then outdir = testdir + 'SlopePlots' + pathdelim

    ; check to see if outdirectory name ends in path delimiter--if not, add it to the end
    if (strmid(outdir, strlen(outdir)-1, 1) ne pathdelim) then outdir = outdir + pathdelim

    ; set defaults for keyword values not provided by the calling context
    if (not keyword_set(infile)) then infile = 'dark.lst'
    if (not keyword_set(valoutfile)) then valoutfile = 'slopevals.txt'
    if (not keyword_set(imgoutfile)) then imgoutfile = 'slopeimgs.txt'
    if (not keyword_set(groupfactor)) then groupfactor = 1

    ; If reference pixel correction is to be done, query the refsub function for
    ; the version number, otherwise, just set a flag indicating it was not done
    if (not keyword_set(noref)) then begin
        ref_ver = ' '
        refsub, 0, sca=' ', version = ref_ver
        if (not keyword_set(ma_refcor)) then begin
            ref_type = 'Spatial averaging'
        endif else begin
            ref_type = 'Maximal averaging'
        endelse
    endif else begin
        ref_ver = 'Not used'
        ref_type = 'N/A'
    endelse

    if (not keyword_set(rowinc)) then begin
        rowinc = 1
    endif
    if (not keyword_set(colinc)) then begin
        colinc = 1
    endif

    ; open input file and read all filenames into an array
    readcol, testdir + infile, imgfiles, format='(A)'

    ; make an array of structs to hold filename and detector temperature for
    ; output slope images
    imginfo = replicate({name:'', dettemp:''}, n_elements(imgfiles))

    ; open output file for slope values
    file_mkdir,outdir
    openw, lun_val, outdir + valoutfile, /get_lun

    ; start by reading just the header of the first image in the list to get
    ; some required information
    fits_read_datacube, testdir + imgfiles[0], dummy, firsthdr, first=0, last=10

    xsize = fix(sxpar(firsthdr, 'naxis1'))
    ysize = fix(sxpar(firsthdr, 'naxis2'))
    scaName = sxpar(firsthdr, 'detname')

    ; define the number of pixels to read in for one group of sample reads
    numpixrd = long(xsize) * ysize

    ; extract SCA type from SCA name
    scaType = strmid(scaName, 0, strpos(scaName, '-'))

    ; if the user didn't set the region, define default regions based on the SCA type
	if (not keyword_set(region)) then begin
		case scaType OF
			'HRG' : region = [4, 4091, 4, 4091]
			'H4RG' : region = [4, 4091, 4, 4091]
			'H1RG' : region = [4, 1019, 4, 1019]
			'H2RG' : region = [4, 2043, 4, 2043]
			'SB304': region = [4, 2051, 0, 2047]
		ELSE   : begin
			print, 'Detector type ' + scaType + $
			' not supported.  Exiting procedure.'
			return
			end
		endcase
	endif

    ; construct a string representation of the region array
    region_str = strtrim(string(region), 2)
    region_out = '[' + region_str[0] + ':' + region_str[1] + ':' $
            + strtrim(string(colinc), 2) + ', ' + region_str[2] + ':' $
            + region_str[3] + ':' + strtrim(string(rowinc), 2) + ']'

    ; write values common to all the images to the header of the output file
    printf, lun_val, ref_ver
    printf, lun_val, ref_type
    printf, lun_val, region_out

    ; construct arrays for cosmic ray detection
    mask10sig = intarr(xsize, ysize)

    ; read mask array if present, otherwise construct mask of
    ; all 1s.
    if (keyword_set(gdmask)) then begin
        fits_read, gdmask, gdmask_data, gdmask_header, /data_only
    endif else begin
        gdmask_data = make_array(xsize, ysize, /integer, value=1)
    endelse

    region_mask = gdmask_data[region[0]:region[1]:colinc, region[2]:region[3]:rowinc]

    ; loop through images listed in input file and process them one frame at
    ; a time
    for imgcount = 0, n_elements(imgfiles) - 1 do begin

        ; To keep it simple, the standard formula from Numerical
        ; Recipes will be used to calculate slope per pixel in an incremental
        ; fashion.
        ;
        ; slope = N * S[x_i * y_i] - S[x_i] * S[y_i] / ( N * S[x_i^2] -
        ;         S[x_i]^2 )

        ; First initialize all sums that we will need. To keep it
        ; simple, I'll keep some redundant sums (e.g., the
        ; information already is available to the program)

        ; declare these running-sum arrays for each file so they'll be zeroed
        ; at the start each time
        Sxy = dblarr(xsize, ysize)
        Sx = dblarr(xsize, ysize)
        Sy = dblarr(xsize, ysize)
        Sxx = dblarr(xsize, ysize)

        ; read the 1st frame header
        fits_read_datacube, testdir + imgfiles[imgcount], 0, imghdr, first=0, last=10

        ; extract information from header
        if (scaType ne 'HRG') and (scaType ne 'H4RG') then begin
	        exposuretime = double(sxpar(imghdr, 'exp_time'))
	        numreads = fix(sxpar(imghdr, 'numutr'))
	        samples_per_read = fix(sxpar(imghdr, 'numsamrd'))
	        dettemp = float(sxpar(imghdr, 'dettemp'))
	        ; Evidently, the last few reads in the dark current experiment
	        ; for H2RG-003-SIPIN are corrupted. This can be seen by displaying the
	        ; corresponding images. In order to address this problem, I (Don) added
	        ;  the following lines of code that throw away the last 10 reads.
	        ; if (scaname eq 'H2RG-003-SIPIN' and numreads gt 12) then numreads = numreads - 10
		endif else begin
	        numreads = fix(sxpar(imghdr, 'NGROUPS'))
	        samples_per_read = fix(sxpar(imghdr, 'NREADS'))
	        if (samples_per_read gt numreads) then begin
	        	numreads = fix(sxpar(imghdr, 'NREADS'))
	        	samples_per_read = fix(sxpar(imghdr, 'NGROUPS'))
	        endif
	        ndrops = fix(sxpar(imghdr, 'NDROPS'))
	        dettemp = float(sxpar(imghdr, 'MOLYTABL'))
	        exposuretime = 4096.*4096./32.*10.e-6*(ndrops*numreads+samples_per_read*numreads)
		endelse

        ; find out how many frame groups there are
        numframes = numreads * samples_per_read
        numgroups = ceil(float(numframes) / groupfactor)
        groupends = indgen(numgroups) * groupfactor + (groupfactor - 1)
        groupends[numgroups - 1] = groupends[numgroups - 1] < (numframes - 1)
        groupstart = 0

        ; create an array containing the time for each group of reads
        readtimes = dindgen(groupfactor, numgroups) / (numreads - 1) * exposuretime

        ; if groupfactor > 1 and the last group isn't completely filled
        ; (i.e., the number of reads in the exposure is not evenly
        ; divisible by the groupfactor), the readtimes array will contain
        ; exposure times for reads that aren't present because it assumes that
        ; all groups contain the same number of reads.  Mark these invalid
        ; exposure times as NaN so they can be ignored in processing.

        overtimes = where(readtimes gt exposuretime)
        if (overtimes[0] ne -1) then readtimes[overtimes] = !values.f_nan

        ; make an array to hold the exposure times for each group.  Each
        ; element is the mean value of the exposure times of the reads
        ; contained in the group.
        grouptimes = fltarr(numgroups)
        for groupnum = 0, numgroups - 1 do begin
        	grouptimes[groupnum] = mean(readtimes[*, groupnum], /nan)
        endfor

        ; read the frame groups in until the 1st non-railed group is encountered
        groupnum = -1
        repeat begin
            groupnum = groupnum + 1
            fits_read_datacube, testdir + imgfiles[imgcount], firstgroup, first = (groupends[groupnum]) * numpixrd, last = (groupends[groupnum + 1]) * numpixrd - 1

            ; transform the first group of reads into an average
            framecount = groupends[groupnum + 1] - groupends[groupnum]
            firstgroup = reform(float(firstgroup), xsize, ysize, framecount)
            firstgroup = total(firstgroup, 3) / framecount

            regionpix = firstgroup[region[0]:region[1]:colinc,region[2]:region[3]:rowinc]

            lowrailpix = where(regionpix lt 0.,numlowrailpix)
            hirailpix = where(regionpix gt 65535.,numhighrailpix)

            lowfraction = numlowrailpix / float(n_elements(regionpix))
            hifraction = numhighrailpix / float(n_elements(regionpix))

            hi_or_low = hifraction ge lowfraction ? 65535. : 0.0

        endrep until ((lowfraction le 0.5 and hifraction le 0.5) $
                OR (groupnum eq (numgroups - 2)))

        ; skip to next image if all reads are railed low
        if (groupnum eq numgroups - 2) then begin

            ; make an array to hold the output combination of group time and
            ; corresponding mode value and raw signal
            outcombo = fltarr(3, numgroups)
            outcombo[0, *] = grouptimes
            outcombo[2, *] = hi_or_low

            ; write out time and signal values for each read to the file,
            ; darkcurrent_analyze can be called later to make slope plots.
            printf, lun_val, testdir+imgfiles[imgcount] + ' ' + strtrim(string(numgroups), 2) + ' ' + strtrim(string(groupfactor), 2)
            printf, lun_val, outcombo, format='(F14.3, F14.3, F14.3)'
            printf, lun_val     ; print blank space to separate data for different images
            continue
        endif

        ; the first group is the zero-level for the dark current, so adjust
        ; the grouptimes to be relative to the time of the first group
        firstgroupnum = groupnum
        starttime = grouptimes[firstgroupnum]
        grouptimes = grouptimes - starttime

        ; find out where the railed pixels are
        railedpix = where(firstgroup gt 65535. OR firstgroup lt 0.,numrailedpix)

        ; exit the procedure if every pixel in the read is railed
        if (numrailedpix eq n_elements(firstgroup)) then begin
            stop, 'Processing cannot continue because every pixel in input ' + $
                    'image is railed at the top or bottom of the A-D range.'
        endif

        ; declare array to hold median and mode values
        theMode=fltarr(numgroups)
        theMedian=fltarr(numgroups)

        ; declare array to hold raw median signal for each read group
        rawsig = fltarr(numgroups)
        rawsig[firstgroupnum] = median(firstgroup[region[0]:region[1]:colinc, region[2]:region[3]:rowinc])

        if (firstgroupnum ne 0) then rawsig[0:firstgroupnum - 1] = 0

        ; set up array used to hold previous read for cosmic ray detection
        previousgroup = fltarr(xsize,ysize)
		; define mask for outliers
        mask10sig[*] = 0

        ; loop through the reads of an individual input file
        for groupcount = firstgroupnum + 1, numgroups - 1 do begin
			print,'Reading group number ',groupcount
            ; read all the frames composing one sample group
            fits_read_datacube, testdir + imgfiles[imgcount], scagroup, first = (groupends[groupcount - 1] + 1) * numpixrd, $
                    last = (groupends[groupcount] + 1) * numpixrd - 1

            ; transform the group of reads into an average
            framecount = groupends[groupcount] - groupends[groupcount - 1]
            scagroup = reform(float(scagroup), xsize, ysize, framecount)
            scagroup = total(scagroup, 3) / framecount

            ; subtract off the first group as bias correction
            scagroup = scagroup - firstgroup

            ; if it was requested, perform reference pixel correction on all reads
            if (not keyword_set(noref)) then begin
                print, 'Applying reference pixel correction to ' + imgfiles[imgcount] + ', read # ' + strtrim(string(groupcount),2)
                if (not keyword_set(ma_refcor)) then begin
                    refsub, scagroup, scaName=scaName
                endif else begin
                    refsub, scagroup, scaName=scaName, /simple
                endelse
            endif else begin
                print, 'No ref. pixel correction for ' + imgfiles[imgcount] + ', read #' + strtrim(string(groupcount),2)
            endelse

            ; record the median signal. Note that firstgroup was already subtracted from scagroup in order to perform reference pixel correction, so we have
            ; to add it back in.
            rawsig[groupcount] = median(scagroup[region[0]:region[1]:colinc, region[2]:region[3]:rowinc]+firstgroup[region[0]:region[1]:colinc, region[2]:region[3]:rowinc])

            ; do cosmic ray detection
            diffFrame = scagroup - previousgroup
            crsdev = stddev(diffFrame[4:xsize-5, 4:ysize-5])
            crmed = median(diffFrame[4:xsize-5, 4:ysize-5])
            outliers = where(diffFrame gt (crmed + crsdev*10))
            if (outliers[0] ne -1) then mask10sig[outliers] = mask10sig[outliers] + 1

			; store the present group as the previous group so that it is available for cosmic ray rejection algorithm when the next group is read
            previousgroup = scagroup

            ; do the incremental slope calculation operations
            Sxy = Sxy + grouptimes[groupcount] * scagroup
            Sx = Sx + grouptimes[groupcount]
            Sy = Sy + scagroup
            Sxx = Sxx + grouptimes[groupcount] ^ 2

            ; extract the region of interest from the frame and apply the mask
            ; (which has already been cut down to the size of the region)
            sca_roi = scagroup[region[0]:region[1]:colinc, region[2]:region[3]:rowinc]
            sca_mask = sca_roi[where(region_mask eq 1)]

            ; begin the process of finding signal values for individual reads
            ; begin by getting the median and stddev using iterative sigma
            ; clipping algorithm.
            djs_iterstat, sca_mask, sigrej=4.0, sigma=djsig, median=djmed
            theMedian[groupcount] = djmed

            ; We want to find the mode. The easiest way to do this is to fit a
            ; gaussian to a histogram.
            ; set min and max histogram bins, centered around the median of the
            ; data
            _min = floor(djmed - (djsig * 4)) + 0.5
            _max = floor(djmed + (djsig * 4)) - 0.5
            _nbins = _max - _min

            ; set up error handler to catch and deal with any errors
            ; caused by gaussfit procedure.
            math_err_stat = check_math(mask = 208)
            math_errcount = 0
            catch, errval
            if (errval ne 0) then begin
                ; if there is a math error which has been thrown manually
                ; by a call to the message procedure, reduce the number of bins
                ; and retry
                if (strpos(!error_state.msg, 'math error') ne -1) then begin
                    _nbins = _nbins / 2
                    message, /reset
                    math_errcount = math_errcount + 1
                endif else begin
                    ; if this is an unexpected non-math error, cancel the
                    ; error handler and retry, allowing the error to be
                    ; handled by the default IDL handler (which should exit the
                    ; procedure).
                    catch, /cancel
                endelse
            endif

            ; make a histogram from the masked data in the region of interest.
            theHist = histogram(sca_mask, nbins=_nbins, min=_min, max=_max, /nan, locations = theVals)

            ; adjust values to center of bins
            theVals = theVals + ((_max - _min) / float(_nbins - 1)) / 2.0

            ; fit the histogram counts and values to a gaussian
            theResult = gaussfit(theVals, theHist, A)
            math_err_stat = check_math(mask = 208)
            if (math_err_stat ne 0 and math_errcount lt 5) then message, 'math error'

            ; we want the mode, which is element 1 of the A array returned from gaussfit
            ;print, math_errcount, a[1]
            if (groupcount lt n_elements(themode)) then begin
                theMode[groupcount] = A[1]
            	print,groupcount,rawsig(groupcount),themedian(groupcount),themode(groupcount)
            endif
            wait,0.5
        endfor

        ; make an array to hold the output combination of group time and
        ; corresponding median value and raw signal
        outcombo = fltarr(4, n_elements(theMode))
        outcombo[0, *] = grouptimes + starttime
        outcombo[1, *] = theMode
        outcombo[2, *] = rawsig
        outcombo[3, *] = theMedian

        ; write out time and signal values for each read,
        ; darkcurrent_analyze can be called later to make slope plots.
        printf, lun_val, testdir+imgfiles[imgcount] + ' ' + strtrim(string(numgroups), 2) + ' ' + strtrim(string(groupfactor), 2)
        printf, lun_val, outcombo, format='(F14.3, F14.3, F14.3, F14.3)'
        printf, lun_val     ; print blank space to separate data for different images

        ; calculate the final per-pixel slope values
        outarray = (numgroups * Sxy - Sx * Sy) / (numgroups * Sxx - Sx * Sx)

        ; change output array to float to make output FITS file smaller
        outarray = float(outarray)

        ; set the detected cosmic rays/hot pixels to NaN
        if (mask10sig[0] ne -1) then begin
            crpix = where(mask10sig ne 0)
            if (crpix[0] ne -1) then outarray[crpix] = !values.f_nan
        endif

        ; correct values in the original header to indicate 1 frame of float
        ; values
        sxaddpar, imghdr, 'naxis3', 1
        sxaddpar, imghdr, 'bitpix', -32
        sxaddpar, imghdr, 'refver', ref_ver
        sxaddpar, imghdr, 'reftype', ref_type

        ; write out slope array with filename equal to original filename
        ; with '_slope' inserted prior to '.fits'
        nameparts = strsplit(imgfiles[imgcount], '.', /extract)
        outname = strjoin(nameparts[0:n_elements(nameparts)-2], '.') + '_slope.fits'
        print,'Writing file ',outdir+outname
        fits_write, outdir + outname, outarray, imghdr

        ; enter the output name and detector temperature (rounded to nearest
        ; 0.1 K) into imginfo array
        imginfo[imgcount].name = outname
        imginfo[imgcount].dettemp = string(dettemp, format='(F5.1)')
    endfor

    ; get rid of null file names (representing images that are railed in every read)
    validnames = where(imginfo.name ne '', numvalid)
    if (numvalid gt 0) then begin
        imginfo = imginfo[validnames]

        ; open output file for names of slope images
        openw, lun_img, outdir + imgoutfile, /get_lun

        ; sort the imginfo struct array by the detector temperature field
        imginfo = imginfo[sort(imginfo.dettemp)]

        ; find the lengths of all the filenames and set the length of the field for
        ; all names to the maximum value + 5 characters (for output formatting)
        namelen = strlen(imginfo.name)
        fieldlen = max(namelen) + 5
        for i=0, n_elements(imginfo) - 1 do begin
            imginfo[i].name = imginfo[i].name + string(make_array(fieldlen - strlen(imginfo[i].name), /byte, val=32b))
        endfor

        ; add EOL character(s) to output dettemp so that output will not all be on
        ; the same line
        imginfo.dettemp = imginfo.dettemp + endofline

        ; write the imginfo struct array out to the image list file
        writeu, lun_img, imginfo.name + endofline
        free_lun, lun_img
    endif

    ; close the output files
    free_lun, lun_val

End
