; NAME:
;   latent_analysis.pro (version 2.7)
;
; PURPOSE:
;   This procedure reads in reduced FITS images generated by latent.pro and
;   produces plots of the results.
;
; CALLING SEQUENCE:
;   latent_analysis, testdir, [outdir, imgfile, region]
;
; INPUTS:
;   testdir     - directory where input files and images are located
;
; KEYWORD PARAMETERS:
;   outdir      - directory to which output is written (assumed same as testdir
;                 if not supplied)
;   infile      - name of input file containing names of reduced images
;                 (default name is 'reduced_images.txt')
;   region      - 1-d integer array with 4 elements specifying rectangular
;                 analysis region in the form[x1, x2, y1, y2]
;
; EXAMPLE
;   Create latency plots for a data set in current directory
;   IDL> latent_analysis, '.', region = [4, 511, 4, 2043]
;
; REFERENCE:
;
; MODIFICATION HISTORY
;   Written by:
;       Ernie Morse, IDTL, December 2002
;   Modified by:
;       Ernie Morse, IDTL, February 15, 2003
;           Added pix keyword for specification of analysis region
;       Ernie Morse, IDTL, March 5-6, 2003 (v2.2)
;           Added plot of cumulative persistence
;           Added annotation for analysis region and reference pixel version &
;           type
;           Added name of plot below x axis
;           Added production of jpeg plots as well as PostScript
;       Ernie Morse, IDTL, March 24, 2003 (v2.3)
;           Changed method of writing JPEGs from using Ghostview to convert PS
;           files to writing to the 'Z' device and then using tvrd and
;           write_jpeg
;       Ernie Morse, IDTL, April 7, 2003 (v2.4)
;           Added SB304 to supported detector types in the default pix selection
;           block
;           Modified JPG writing section to make higher-resolution plots
;       Ernie Morse, IDTL, April 15, 2003 (v2.5)
;           Removed zero-bin from consideration in histogram fitting
;       Ernie Morse, IDTL, April 16, 2003 (v2.5)
;           zero-bin exclusion didn't work, so it was removed
;           NaN values are written into reduced files by latent_reduce.  These
;           values won't be included in histogram
;       Ernie Morse, IDTL, May 28, 2003(v2.6)
;           Changed pix keyword to region to match convention of other
;           reduction scripts
;           Changed existing region variable to sigma_region
;       Ernie Morse, IDTL, February 10, 2004 (v2.7)
;           Edited to ensure compliance with IDTL coding standards.
;           Removed special region setting for H2RG-003-4.8mu, as that device
;           has been destroyed.
;           Set nan keywords on calls to stddev instead of using
;           where(finite()) to select finite portion of input
;       Ernie Morse, IDTL, May 18, 2004 (v2.8)
;           Added NaN flag to histogram command to enable handling of those values
;
;-
Pro latent_analysis, testdir, outdir = outdir, infile = infile, region = region, $
        rowinc = rowinc

    delim =path_sep()

    ;set defaults for unspecified keyword values
    if (NOT(keyword_set(outdir))) then begin
        outdir = testdir
    endif
    if (NOT(keyword_set(infile))) then begin
        infile = 'reduced_files.txt'
    endif

    if (not(keyword_set(rowinc))) then begin
        rowinc = 1
    endif

    ; save current device name so it can be reset later
    dname = !D.NAME

    ; read images and sequence numbers to analyze from input file
    readcol, testdir + delim + infile, files, format = '(A)'

    ; get detector type from header of 1st image
    fits_read, testdir + delim + files[0], 0, h, /header_only
    detname = sxpar(h, 'detname')
    det_type = strmid(detname, 0, strpos(detname, '-'))
    ref_ver = strtrim(sxpar(h, 'ref_ver'), 2)
    ref_typ = strtrim(sxpar(h, 'ref_typ'), 2)

    ; assign default analysis regions based on detector type if the calling
    ; context has not set the region keyword
    if (not keyword_set(region)) then begin
        case det_type of
            'H1RG': begin
                region = [4, 1019, 4, 1019]
            end
            'H2RG': begin
                region = [4, 2043, 4, 2043]
            end
            'SB304': begin
                region = [4, 2051, 0, 2047]
            end
            else : begin
                print, 'Detector type not supported.  Exiting procedure.'
                return
            end
        endcase
    endif

    ; make string representation of the region for use in annotating plots
    region_str = strarr(4)
    for i = 0, 3 do begin
        region_str[i] = strtrim(string(region[i]), 2)
    endfor
    region_out = '[' + region_str[0] + ':' + region_str[1] + ', ' $
                 + region_str[2] + ':' + region_str[3] + ']'

    ; processing loop for input reduced images
    for i = 0, n_elements(files)-1 do begin

        ; read header of 1st image and retrieve selected values from it
        fits_read, testdir + delim + files[i], 0, hd, /header_only
        temp = sxpar(hd, 'dettemp')
        detname = sxpar(hd, 'detname')
        exptime = float(sxpar(hd, 'exp_time'))
        cols = fix(sxpar(hd, 'naxis1'))
        rows = fix(sxpar(hd, 'naxis2'))
        nreads = fix(sxpar(hd, 'naxis3'))
        fluence = float(sxpar(hd, 'fluence'))
        level = float(sxpar(hd, 'flu_lvl'))
        wheel1 = strtrim(sxpar(hd, 'prefilt1'), 2)
        wheel2 = strtrim(sxpar(hd, 'prefilt2'), 2)
        flu_time = strtrim(string(float(sxpar(hd, 'flu_time')), $
                format = '(F6.1)'), 2)
        t_prior = strtrim(string(float(sxpar(hd, 't_prior')), format = '(F4.1)'), 2)
        n_resets = strtrim(string(sxpar(hd, 'n_resets') + 1, $
                format = '(I3)'), 2)

        ; since IDTL headers are non-standard, the following two function calls
        ; are necessary to get the date and time
        gettime, hd, hour
        getdate, hd, day

        ; determine size of one read
        readsize = long(cols) * rows

        ; read 1st frame of image and reform into 2-D array
        fits_read, testdir + delim + files[i], darksub, first = 0, $
                last = readsize - 1
        darksub = reform(darksub, cols, rows, /overwrite)

        if (n_resets eq '0') then n_resets = 'autoflush'

        ; create arrays for exposure times and persistent charge at each read
        time = (findgen(nreads) + 1) * (exptime / (nreads))
        subcounts = fltarr(nreads)
        cumulative_counts = fltarr(nreads + 1)


        ; find median dark-subtracted counts at each UTR interval
        ; set up histogram parameters for finding the mode
        ; find standard deviation and median of difference array
        ; standard deviation doesn't like the NaN values, so set the
        ; nan keyword to get stddev function to ignore them
        sigma_region = darksub[region[0]:region[1], region[2]:region[3]:rowinc]
        sigma = stddev(sigma_region, /nan)
        med = median(darksub[region[0]:region[1], region[2]:region[3]:rowinc], /even)

        ; use sigma to set histogram limits.  Sigma*5 should be inclusive enough.
        ; Set binsize to 1.0 for now.
;        _binsize = 1.0
        _min = long(med - (sigma * 5))
        _max = long(med + (sigma * 5))
        _nbins = 128

        theHist = histogram(darksub[region[0]:region[1], $
                region[2]:region[3]:rowinc], min = _min, max = _max, $
                nbins = _nbins, locations = theVals, /nan)


        ; adjust values to center of bins
        theVals = theVals + ((_max - _min) / float(_nbins - 1)) / 2.0

        theFit = gaussfit(theVals, theHist, C)
        subcounts[0] = (cumulative_counts[1] = C[1])
        prevdarksub = darksub
        for k = 1, nreads-1 do begin

            ; read one frame of reduced image
            fits_read, testdir + delim + files[i], darksub, $
                    first = readsize * k, last = readsize * (k + 1) - 1
            darksub = reform(darksub, cols, rows, /overwrite)

            ; subtract the previous read
            intervalcounts = darksub - prevdarksub

            ; save this read as the previous for the next time through
            prevdarksub = darksub

            ; find median dark-subtracted counts at each UTR interval
            ; set up histogram parameters for finding the mode
            ; find standard deviation and median of difference array
            ; standard deviation doesn't like the NaN values, so set the nan
            ; keyword for the stddev function
            sigma_region $
                    = intervalcounts[region[0]:region[1], $
                    region[2]:region[3]:rowinc]
            sigma = stddev(sigma_region, /nan)
            med = median(intervalcounts[region[0]:region[1], $
                    region[2]:region[3]:rowinc], /even)

            ; use sigma to set histogram limits.  Sigma*5 should be inclusive
            ; enough.
            ; Set binsize to 1.0 for now.
            ;_binsize = 1.0
            _min = long(med - (sigma * 5))
            _max = long(med + (sigma * 5))

            theHist = histogram(intervalcounts[region[0]:region[1], $
                    region[2]:region[3]:rowinc], min = _min, max = _max, $
                    nbins = _nbins, locations = theVals, /nan)

            ; adjust values to center of bins
            theVals = theVals + ((_max - _min) / float(_nbins - 1)) / 2.0

            theFit = gaussfit(theVals, theHist, C)
            subcounts[k] = C[1]
            cumulative_counts[k + 1] = cumulative_counts[k] + subcounts[k]
        endfor

        ; divide all interval values by fluence to get percentage of
        ; persistence for interval
        percents = (subcounts / fluence) * 100
        cumulative_percents = (cumulative_counts / fluence) * 100

        ; set up plot limits so that plot data takes up 80% of plot window
        pmax = max([cumulative_percents,percents])
        pmin = min([percents,cumulative_percents])
        tmax = max(time)
        tmin = min(time)
        yrng = (pmax - pmin) / 0.8
        ylo = (0.0 - yrng * 0.2) < (pmin - yrng * 0.2)
        yhi = pmax + yrng * 0.1
        xrng = (tmax - tmin) / 0.8
        xlo = tmin - xrng * 0.1
        xhi = tmax + xrng * 0.1

        ; set plot name
        psname = outdir + delim + 'Latent_' + wheel1 + '_' + wheel2 + '_' + $
                strtrim(string(round(level * 100)), 2) + '_' $
                + strtrim(string(round(temp)), 2) + 'K.ps'
        jpgname = outdir + delim + 'Latent_' + wheel1 + '_' + wheel2 + '_' + $
                strtrim(string(round(level * 100)), 2) + '_' $
                + strtrim(string(round(temp)), 2) + 'K.jpg'

        ; set plot colors
        black = 0
        white = 255

        ; make two plots, one PS and one JPEG
        for j = 1, 1 do begin
            if (j eq 0) then begin

                ; set up for landscape Postscript plotting
                set_plot, 'ps'
                !p.font = -1
                device, filename = psname, /landscape

                ; fsub is filename to be written at bottom of plot
                fsub = psname

                !x.minor = 10
                !y.minor = 5
                !y.margin = [6, 2]
                !x.margin = [14, 4]

                ; set some size values
                asize = 1.12        ; charsize of annotations
                athick = 4.0        ; charthick of annotations
                lthick = 8.0        ; line thickness
            endif else begin

                ; set up to plot to memory buffer.  This plot will later be
                ; written from memory to a JPEG file.
                set_plot, 'z'
                device, set_resolution = [8000, 6000]
                device, set_font = 'Courier'

                ; fsub is filename to be written at bottom of plot
                fsub = jpgname

                ; set plot margin & tick marks
                !Y.MINOR = 5
                !x.minor = 10
                !y.margin = [6, 4]
                !x.margin = [14, 4]

                ; set some size values
                asize = 8.0         ; charsize of annotations
                athick = 15.0       ; charthick of annotations
                lthick = 20.0       ; line thickness
            endelse

            fval = fluence + 500.0
            fval = long(fval - (fval mod 1000))
            ptitle = wheel1 + ' + ' + wheel2 + ' persistence, ' $
                    + strtrim(string(fval), 2) + ' ADU fluence'

            ; plot differential persistence
            plot, time, percents, psym = -5, xtitle = 'Time (seconds)', $
                    ytitle = '% of fluence liberated', $
                    title = ptitle, yrange = [ylo, yhi], xrange = [xlo, xhi], $
                    xstyle = 1, ystyle = 1, $
                    BACKGROUND = white, $
                    charsize = asize * 1.25, charthick = athick, $
                    xthick = lthick, ythick = lthick, $
                    COLOR = black, symsize = asize, thick = lthick

            ; plot cumulative persistence
            oplot, [0, time], cumulative_percents, psym = -2, linestyle = 2, $
                    color = black, symsize = asize, thick = lthick

            ; plot the line through y = 0
            oplot, [xlo, xhi], [0.0, 0.0], thick = lthick * 1.5, color = black,$
                    symsize = asize

            ; plot legend symbols
            oplot, [xlo + xrng * 0.6, xlo + xrng * 0.65], $
                    [yhi-yrng * 0.075, yhi-yrng * 0.075], linestyle = 2, $
                    color = black, thick = lthick
            oplot, [xlo + xrng * 0.6, xlo + xrng * 0.65], $
                    [yhi-yrng * 0.045, yhi-yrng * 0.045], color = black, $
                    thick = lthick

            ; put annotations on the lower left portion of the plot
            xyouts, xlo + (xrng * 0.05), ylo + (yrng * 0.145), $
                    'Fluence Time: ' + flu_time + ' seconds', $
                    charsize = asize, color = black, charthick = athick
            xyouts, xlo + (xrng * 0.05), ylo + (yrng * 0.115), $
                    'Delay: ' + t_prior + ' seconds, Resets:  ' + n_resets, $
                    charsize = asize, color = black, charthick = athick
            xyouts, xlo + (xrng * 0.05), ylo + (yrng * 0.055), $
                    'Source:  ' + files[i], charsize = asize, color = black, $
                    charthick = athick
            xyouts, xlo + (xrng * 0.05), ylo + (yrng * 0.085), $
                    'Date, Time (UT): ' + day + ', ' + hour, charsize = asize, $
                    color = black, charthick = athick
            xyouts, xlo + (xrng * 0.05), ylo + (yrng * 0.025), $
                    '1 ADU = ' + strtrim(string(100.0 / fluence), 2) + '%', $
                    charsize = asize, color = black, charthick = athick

            ; put annotations on the upper left part of the plot
            xyouts, xlo + (xrng * 0.05), yhi - (yrng * 0.075), $
                    'Detector Temperature:  ' $
                    + strtrim(string(temp, format = '(F7.3)'), 2) + ' K', $
                    charsize = asize, color = black, charthick = athick
            xyouts, xlo + (xrng * 0.05), yhi - (yrng * 0.045), $
                    'Detector:  ' + detname, charsize = asize, color = black, $
                    charthick = athick

            ; put annotations next to the legend symbols at upper right
            xyouts, xlo + (xrng * 0.67), yhi - (yrng * 0.045), $
                    '= Differential', charsize = asize, color = black, $
                    charthick = athick
            xyouts, xlo + (xrng * 0.67), yhi - (yrng * 0.075), $
                    '= Cumulative', charsize = asize, color = black, $
                    charthick = athick

            ; put annotations on the lower right part of the plot
            xyouts, xlo + (xrng * 0.95), ylo + (yrng * 0.145), $
            'SCA Region:  ' + region_out, charsize = asize, alignment = 1.0, $
                    color = black, charthick = athick
            if (ref_ver ne '0') then begin
                xyouts, xlo + (xrng * 0.95), ylo + (yrng * 0.115), $
                        'Refsub Version: ' + ref_ver, charsize = asize, $
                        alignment = 1.0, color = black, charthick = athick
            endif
            if (ref_typ ne '0') then begin
                xyouts, xlo + (xrng * 0.95), ylo + (yrng * 0.085), $
                        'Refsub Type: ' + ref_typ , charsize = asize, $
                        alignment = 1.0, color = black, charthick = athick
            endif

            ; write name of file outside of plot margins at bottom of page
            xyouts, 1.0, 40.0, fsub, charsize = asize, charthick = athick, $
                    color = black, /device

            ; For the JPEG plot, the plot image must be read in from memory and
            ; written out to a file
            if (j eq 1) then begin
                xyouts, 1.0, 0.0, '_', charsize = 1.1, charthick = 1.8, /device
                jpgimg = tvrd()
                write_jpeg, jpgname, congrid(jpgimg, 1200, 960, /interp, $
                        /center), quality = 100
            endif

           ; finish up plot by closing device
           ; device must be open for the above JPEG code to work, so close must
           ; happen here at the end
           device, /close
        endfor      ; j loop
    endfor          ; i loop

    ; set plot back to the default device
    set_plot, dname
End