; NAME:
;   detector_expose
;
; PURPOSE:
;   Driver program that creates h2rg detector object, loads input count-rate and
;   dark current images and produces a simulated detector readout with the
;   same exposure time and read mode as the dark current image.  Output is in
;   the form of a FITS file for each read of the detector.
; 
; INPUTS:
;   none
; 
; KEYWORDS:
;   im_countrate    full path name of input count-rate FITS image.  Must be
;                   a 2048 x 2048 floating-point array.  User will be prompted
;                   with a graphic dialog pickfile box if this keyword is not
;                   set.
;
;   im_dark         full path name of input dark current image.  User is
;                   prompted if it is not provided.
;
;   outdir          directory where output files should be written.  User is
;                   prompted if it is not provided.
;
;   addcr           set to add simulated cosmic rays to the output images.
;
;   crout           set to write out FITS images containing only the cosmic rays
;                   added to each read of the image
;
; EXTRA KEYWORDS (optional, passed on to other procedures called by this one)
;   rate            cosmic ray incidence rate in number of hits per cm^2 per
;                   second--has no effect if addcr is not set.  If not set,
;                   the default value of 5/cm^2/second is used.
;
;   mask            set to enable charge diffusion to neighboring pixels from
;                   cosmic ray hits--has no effect if addcr is not set.
;
; OUTPUTS:
;   none
;
Pro detector_expose, im_countrate = im_countrate, im_dark = im_dark, $
        outdir = outdir, addcr = addcr, crout = crout, _ref_extra = crkeys

    ; define an error handler for all of the code that follows (including
    ; nested procedure or function calls).
    ; when an error occurs, control is transferred here.
    catch, errval
    if (errval ne 0) then begin

        ; free up the memory used by objects and pointers
        if (n_elements(zeromed) ne 0) then begin
            ptr_free, zeromed
        endif
        if (n_elements(H2RG_SCA) ne 0) then begin
            obj_destroy, H2RG_SCA
        endif

        ; print message and stop execution
        print, !error_state.msg
        return
    endif

    ; determine the proper path delimiter to use based on OS
    if (!version.os_family eq 'Windows') then begin
        delim = '\'
    endif else begin
        delim = '/'
    endelse

    ; get location of output directory (if not provided on command line)
    if (not(keyword_set(outdir))) then begin
        outdir = dialog_pickfile(title = $
            'Please select directory for output.', /dir)

        ; throw an error if the user hit the Cancel button in the dialog
        if (outdir eq '') then begin
            message, 'Output directory not selected.  Exiting.'
        endif
        
    endif else begin
        
        ; make sure the directory name provided on the command line ends with
        ; the proper path delimiter.
        if (strmid(outdir, strlen(outdir) - 1) ne delim) then begin
            outdir = outdir + delim
        endif
    endelse

    ; wait a couple of seconds so user doesn't think that the same dialog is
    ; popping up again
    wait, 2

    ; get location of input count-rate image (if not provided on command line)
    if (not(keyword_set(im_countrate))) then begin
        im_countrate = dialog_pickfile(title = $
            'Please select input count-rate image.', path = outdir, /must_exist)

        if (im_countrate eq '') then begin
            message, 'Input countrate image not selected.  Exiting.'
        endif

    endif


    wait, 2

    ; get location of dark current image (if not provided on command line)
    if (not(keyword_set(im_dark))) then begin
        im_dark = dialog_pickfile(title = $
            'Please select input dark current image.', path = outdir, $
            /must_exist)

        if (im_dark eq '') then begin
            message, 'Input dark current image not selected.  Exiting.'
        endif

    endif

    ; create a text string from the current date and time.  This will serve as a
    ; a base for file names of output images.
    caldat, systime(/julian), month, day, year, hour, minute

    timestrs = [$
        ptr_new(strtrim(string(year), 2)), $
        ptr_new(string(month, format='(I2)')), $
        ptr_new(string(day, format='(I2)')), $
        ptr_new(string(hour, format='(I2)')), $
        ptr_new(string(minute, format='(I2)'))]

    for tc = 1, 4 do begin
        if (strpos(*timestrs[tc], ' ') ne -1) then begin
            strput, *timestrs[tc], '0', 0
        endif
    endfor

    basename = outdir + *timestrs[0] + '_' + *timestrs[1] + '_' + $
        *timestrs[2] + '_' + *timestrs[3] + '_' + *timestrs[4] + '_'

    ptr_free, timestrs

    ; read the count-rate image 
    fits_read, im_countrate, countimg

    ; get some parameters from the dark current image header
    fits_read, im_dark, 0, dkhdr, /header_only
    numutr = sxpar(dkhdr, 'numutr')
    numsamrd = sxpar(dkhdr, 'numsamrd')
    exp_time = sxpar(dkhdr, 'exp_time')
    xsize = sxpar(dkhdr, 'naxis1')
    ysize = sxpar(dkhdr, 'naxis2')
    dettemp = sxpar(dkhdr, 'dettemp')

    ; get the name of the median frame that goes along with the dark image.
    ; This frame contains the median zero-signal value of each
    ; detector pixel.
    medname = sxpar(dkhdr, 'medimg')

    ; get the name of the flat field frame to use with this dark current image.
    ; this frame will be used to recreate the pixel sensitivity differences of
    ; the array
    flatname = 'flat_' + strtrim(string(round(dettemp)), 2) + 'K.fits'

    ; extract the directory path from the name of the dark current image
    ; and prepend it to the names of the median and flat images
    indir = strmid(im_dark, 0, strpos(im_dark, delim, /reverse_search) + 1)
    medname = indir + medname
    flatname = indir + flatname

    ; prepare an array on the heap to hold the zero median image and then
    ; read the image into it.  this is a pointer so that either an array or a 
    ; single value that applies to all pixels can be used by the same code.
    zeromed = ptr_new(fltarr(xsize, ysize))
    fits_read, medname, *zeromed

    ; get the flat-field image also. this one doesn't need to be a pointer.
    fits_read, flatname, flatimg

    ; create H2RG detector object
    H2RG_SCA = obj_new('h2rg', 'd1', gain = 1.3, welldepth = 1.35e5)

    ; reset detector pixel array
    H2RG_SCA->reset_pix

    ; get time to read out detector array
    t_frm = H2RG_SCA->get_frametime()

    ; calculate timing of individual reads in the exposure
    rd_times = H2RG_SCA->calc_readtimes(numsamrd, numutr, exp_time, t_frm)
    
    ; perform the required number of reads of the array and write the results
    ; to individual FITS files for each read.
    for readcount = 0, numutr * numsamrd - 1 do begin
    
        ; calculate interval time for this read
        if (readcount eq 0) then begin
            t_int = t_frm
        endif else begin
            t_int = rd_times[readcount] - rd_times[readcount - 1]
        endelse

        ; calculate accumulated signal since last read
        sigaccum = countimg * t_int

        ; add signal accumulated since last read
        H2RG_SCA->add, sigaccum

        ; calculate and add shot noise
        H2RG_SCA->add, calc_shot_noise(sigaccum)

        ; add cosmic rays if the addcr keyword is set
        if (keyword_set(addcr)) then begin
            crplane = add_cosmic_rays(xsize, ysize, t_int, $
                    H2RG_SCA->get_pixdim(), _extra = crkeys)            
            H2RG_SCA->add, crplane
        endif
        
        ; read out array
        outimg = H2RG_SCA->readout(readcount, t_int, zeromed, $
            darkimg = im_dark, flatimg = flatimg)

        ; write out the FITS image
        imname = basename + H2RG_SCA->get_name() + '_' + $
                strtrim(string(readcount), 2) + '.fits' 
        hdr = make_header(outimg, rd_times[readcount], numutr, numsamrd, $
                H2RG_SCA->get_name())
        fits_write, imname, outimg, hdr

        ; write out cosmic ray image if it was generated and keyword for
        ; cosmic ray output was set.
        if (n_elements(crplane) ne 0 and keyword_set(crout)) then begin
            
            crimname = basename + H2RG_SCA->get_name() + '_' + $
                       strtrim(string(readcount), 2) + '_cr.fits'

            ; change the header to reflect the fact that the cr image array
            ; is floating point
            sxaddpar, hdr, 'BITPIX', -32 
            fits_write, crimname, crplane, hdr
        endif

    endfor

    ; free up the memory used by objects and pointers
    ptr_free, zeromed
    obj_destroy, H2RG_SCA

end