; NAME:
;   latent_reduce.pro (version 4.4)
;
; PURPOSE:
;   This procedure reads creates reduced latent charge images from the latency
;   and dark exposures generated by the latency experiment script,
;   latent_test.pro.  Lists of file names for latency and dark exposures are
;   read in from the input files.  Both dark and latency images are processed to
;   remove residual bias and bias drift effects, darks are subtracted from
;   latency images, and the resulting processed images are written to FITS
;   files.
; CALLING SEQUENCE:
;   latent, testdir, [outdir, imgfile, darkfile, noref, ma_refcor]
;
; INPUTS:
;   testdir     - directory where input files and images are located
;
; KEYWORD PARAMETERS:
;   outdir      - directory to which output is written (default is a
;                 subdirectory of testdir named 'Results')
;   imgfile     - name of input file containing names of latency images
;                 (default name is 'filelist.txt')
;   darkfile    - name of input file containing names of dark images
;                 (default name is 'darklist.txt')
;   noref       - set to disable reference pixel correction.
;   ma_refcor   - flag to be set if maximal averaging reference pixel correction
;                 is to be done.
;
;   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
;
;   Create reduced latency images from data in the current directory using the
;   default input files. Write data to subdirectory names "Results"
;
;   latent, '.', outdir = 'Results'
;
; REFERENCE:
;
; MODIFICATION HISTORY:
;   Written by:
;       Ernie Morse, IDTL, November, 2002
;   Modified:
;       Ernie Morse, IDTL, Dec. 10, 2002
;           Added explanatory header
;       Ernie Morse, IDTL, March 7, 2003
;           Added keywords to allow specification of type of reference pixel
;           correction
;           Added call of refsub to check version # of procedure
;           Added writing of refsub version and type into FITS headers of
;           reduced images
;       Ernie Morse, IDTL, April 11, 2003 (v4.1)
;           Changed call to refsub to input detector name instead of detector
;           type
;       Ernie Morse, IDTL, April 16, 2003 (v4.2)
;           Made modifications to write NaN in the output image wherever railed
;           pixels appear
;       Ernie Morse, IDTL, May 9, 2003 (v4.3)
;           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 + 'Results'
;       Ernie Morse, IDTL, February 9, 2004 (v4.4)
;           Edited to ensure compliance with IDTL coding standards.
;
pro latent_reduce, testdir, outdir = outdir, imgfile = imgfile, $
        darkfile = darkfile, noref = noref, ma_refcor = ma_refcor

    pathdelim = path_sep()

    ; 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 begin
        testdir = testdir + pathdelim
    endif

    ; set outdir to testdir + 'Results' if not set by user
    if (keyword_set(outdir) eq 0) then begin
        outdir = testdir + 'Results' + pathdelim
    endif
    file_mkdir,outdir

    ; 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 begin
        outdir = outdir + pathdelim
    endif

    ;if keyword parameters are not provided, set them to defaults
    if (NOT(keyword_set(imgfile))) then begin
        imgfile = 'filelist.txt'
    endif
    if (NOT(keyword_set(darkfile))) then begin
        darkfile = 'darklist.txt'
    endif

    ; 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

        ; set up non-empty string required when calling refsub to get version
        ; number.  The version number will be returned to this string by the
        ; call to refsub with the version keyword set.
        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

    ; read names of latency and dark images
    readcol, testdir + imgfile, input, format = 'A'
    readcol, testdir + darkfile, dkimgs, format = 'A'

    ; make array of structs to hold image file names and temperatures
    images = replicate({img, name:'', temp:0}, n_elements(input))
    images.name = input
    for m = 0, n_elements(images)-1 do begin
        fits_read, testdir + images[m].name, 0, hdr, /header_only
        images[m].temp = round(sxpar(hdr, 'dettemp'))
    endfor

    ; make array of structs to hold dark image file names and temperatures
    darks = replicate({dk, name:'', temp:0}, n_elements(dkimgs))
    darks.name = dkimgs
    for n = 0, n_elements(dkimgs)-1 do begin
        fits_read, testdir + darks[n].name, 0, hdr, /header_only
        darks[n].temp = round(sxpar(hdr, 'dettemp'))
    endfor

    ; sort images by temperature
    sortorder = sort(images.temp)
    images = images[sortorder]

    ; get image dimensions and SCA type from FITS header of last-read image
    cols = fix(sxpar(hdr, 'naxis1'))
    rows = fix(sxpar(hdr, 'naxis2'))
    nreads = fix(sxpar(hdr, 'naxis3'))
    scaName = sxpar(hdr, 'detname')

	; open file for results
    openw, redu, outdir + 'reduced_files.txt', /get_lun

    ;loop through all files and create reduced images

    current_temp = 0
    framesize = long(cols) * rows

    ; create array to hold all darks of current temperature
    darksubs = ptr_new(fltarr(cols, rows, nreads - 1))
    darkmed = darksubs

    for i = 0, n_elements(images) - 1 do begin

        ; make combined dark current image for each new detector temperature
        ; encountered
        if (images[i].temp ne current_temp) then begin
            current_temp = images[i].temp
            print, 'Reading dark current images for T = ' + $
                    strtrim(string(current_temp), 2) + 'K'

            ; find all dark files for current temperature
            current_dks = where(darks.temp eq current_temp)

            alldarks = fltarr(cols, rows, n_elements(current_dks))
            zeroreads = intarr(framesize, n_elements(current_dks))

            ; read and process all current darks
            for rc = 0, nreads - 1 do begin
                for dc = 0, n_elements(current_dks) - 1 do begin
                    fits_read, testdir + darks[current_dks[dc]].name, $
                            darkcurrentimg, first = rc * framesize, $
                            last = (rc + 1) * framesize - 1
                    if (rc eq 0) then begin
                        zeroreads[*, dc] = darkcurrentimg
                        alldarks[*, *, dc] = darkcurrentimg
                    endif else begin

                        ; residual bias subtraction
                        alldarks[*, *, dc] = darkcurrentimg - zeroreads[*, dc]

                        ; refpixel correction
                        if (not keyword_set(noref)) then begin
                            if (not keyword_set(ma_refcor)) then begin
                                simopt = 0
                            endif else begin
                                simopt = 1
                            endelse
                            refcorimg = alldarks[*, *, dc]
                            refsub, refcorimg, sca = scaName, simple = simopt
                            alldarks[*, *, dc] = refcorimg
                        endif
                    endelse
                endfor

                ; make a median of this read across the three darks for all
                ; reads after the 0th (if there are more than 1 darks)
                if (n_elements(current_dks) ne 1) then begin
                    if (rc ne 0) then begin
                        (*darkmed)[*, *, rc - 1] = $
                                median(alldarks, dimension = 3, /even)
                    endif
                endif else begin
                    if (rc ne 0) then begin
                        (*darkmed)[*, *, rc - 1] = reform(alldarks)
                    endif
                endelse
            endfor

            ; write median dark to a FITS file
            fits_write_ptr, outdir + 'darkmedian.fits', darkmed
        endif else begin

            ; read the previously created median dark image in from the file and
            ; put it in memory
            for md = 0, nreads - 2 do begin
                fits_read, outdir + 'darkmedian.fits', refcorimg, $
                        first = md * framesize, last = (md + 1) * framesize - 1
                (*darkmed)[*, *, md] = refcorimg
            endfor
        endelse

        ; zero point correction
        for zpc = 0, nreads - 1 do begin
            fits_read, testdir + images[i].name, image, h, $
                    first = zpc * framesize, last = (zpc + 1) * framesize - 1
            if (zpc eq 0) then begin
                zeroreads[*, 0] = image
                railedpix = where(zeroreads[*, 0] gt 32760 $
                                  OR zeroreads[*, 0] lt -32760)
            endif else begin
                image = temporary(image) - zeroreads[*, 0]
                image = reform(image, cols, rows, /overwrite)

                ; reference pixel correction
                if (not keyword_set(noref)) then begin
                    if (not keyword_set(ma_refcor)) then begin
                        simopt = 0
                    endif else begin
                        simopt= 1
                    endelse
                    refsub, image, sca = scaName, simple = simopt
                endif
                (*darksubs)[*, *, zpc - 1] = image - (*darkmed)[*, *, zpc - 1]
                (*darksubs)[railedpix + framesize*(zpc - 1)] = !values.f_nan
            endelse
        endfor      ; zpc loop

        ; get components of output file name and construct file name from them
        pf1 = strtrim(sxpar(h, 'prefilt1'), 2)
        pf2 = strtrim(sxpar(h, 'prefilt2'), 2)
        fval = long(float(sxpar(h, 'fluence')) + 500.0)
        fluence = strtrim(string(fval - (fval mod 1000)), 2)
        reducedfile = 'Latent_reduced_' + pf1 + '_' + pf2 + '_' + fluence $
                + '_' + strtrim(string(current_temp), 2) + '.fits'

        ; make some corrections to the header to stop writefits from complaining
        sxaddpar, h, 'NAXIS3', nreads - 1
        sxaddpar, h, 'BITPIX', -32
        sxaddpar, h, 'REF_VER', ref_ver
        sxaddpar, h, 'REF_TYP', ref_type

        ; write out reduced image + header and make entry into output text file
        fits_write_ptr, outdir + reducedfile, darksubs, h
        printf, redu, reducedfile, format = '(A)'
    endfor

    ptr_free, darksubs

    free_lun, redu

    ; run analysis code
    ; region is hardwired for H1RG-022-SIPIN
    latent_analysis, outdir, region=[100,900,400,900]

end