; NAME:
;	darkcurrent_perpix_analysis.pro (v1.0)
;
;
; PURPOSE:
;	To produce composite per-pixel dark current and noise images from a
;	collection of dark current slope images at a given detector temperature.
;	The input slope images are produced by darkcurrent_reduce.pro
;
;
; CALLING SEQUENCE:
;	darkcurrent_perpix_analysis, testdir [, infile[, outdir]]
;
;
; INPUTS:
;	testdir		the full path of the directory where the input data is stored.
;
;
;
; KEYWORD PARAMETERS:
;	infile		the file name of the list of input images.  The file is assumed
;				to consist of two columns, with the first column containing
;				names of input darkcurrent slope images produced by
;				darkcurrent_reduce and the second column containing the detector
;				temperature for the raw data exposure from which the slope
;				image was generated.  It is further assumed that the images
;				are sorted in order of increasing detector temperature.  The
;				default value for this variable is 'slopeimgs.txt'
;
;	outdir		the name of the directory where output images are written. The
;				default is equal to the value of the testdir input variable.
;
;
; OUTPUTS:
;	Produces a composite dark current image and a composite noise image for
;	each unique detector temperature given in the input file.
;
;
; EXAMPLE:
;	IDL > darkcurrent_perpix_analysis, '.'
;
;
; MODIFICATION HISTORY:
;	Created by:
;		Ernie Morse, IDTL, October 10, 2003
;1234567890123456789012345678901234567890123456789012345678901234567890123456789
Pro darkcurrent_perpix_analysis, testdir, infile=infile, outdir=outdir
    ; set directory path delimiter & 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
		pathdelim = '/'
		endofline = string(10b)
	endif else begin
		pathdelim = '\'
		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 begin
	    testdir = testdir + pathdelim
	endif

	; set outdir to testdir if not set by user
	if (keyword_set(outdir) eq 0) then begin
		outdir = testdir
	endif

	; set defaults for keyword values not provided by the calling context
	if (not keyword_set(infile)) then begin
		infile = 'slopeimgs.txt'
	endif

	; read the names and temperatures from the input file
	readcol, testdir + infile, slopeimgs, dettemps, format = '(A, F)'

	; find the last index where each unique detector temperature occurs.
	lastindex = uniq(dettemps)

	; loop through input images, creating composite image for each detector
	; temperature.
	firstindex = 0
	for tempcount = 0, n_elements(lastindex) - 1 do begin

		; get the detector temperature and find out how many images there are
		; for it
		temp = dettemps[lastindex[tempcount]]
	    numimgs = lastindex[tempcount] - firstindex + 1

		; read header of first image to find out the frame size and
		fits_read, testdir + slopeimgs[firstindex], 0, hdr, /header_only
		xsize = sxpar(hdr, 'naxis1')
		ysize = sxpar(hdr, 'naxis2')
		detname = sxpar(hdr, 'detname')

		; construct a cube to hold all of the images at this temperature
		tempcube = ptr_new(fltarr(xsize, ysize, numimgs))

		; read all of the images at this temperature into the big cube
		for imgcount = firstindex, lastindex[tempcount] do begin
		    fits_read, testdir + slopeimgs[imgcount], scaframe
		    (*tempcube)[*,*, imgcount - firstindex] = scaframe
		endfor

		; make the average and standard deviation images, and a frame
		; that shows how many pixels were used to create them
		; start out with all pixels = NaN, and replace the good pixels
		; with the correct values
		meanimg = ptr_new(make_array(xsize, ysize, /float, val = !values.f_nan))
		stdimg = ptr_new(make_array(xsize, ysize, /float, val = !values.f_nan))
		countimg = ptr_new(intarr(xsize, ysize))

		for row = 0, ysize - 1 do begin
		    for col = 0, xsize - 1 do begin
		        goodpix = where(finite((*tempcube)[col, row, *]) gt 0)
		        if (goodpix[0] eq -1) then begin
		            (*countimg)[col, row] = 0
		        endif else begin
		            (*countimg)[col, row] = n_elements(goodpix)
		        endelse
		        if ((*countimg)[col, row] ge 2) then begin
		            (*meanimg)[col, row] = mean((*tempcube)[col, row, *], /nan)
		            (*stdimg)[col, row] = stddev((*tempcube)[col, row, *], /nan)
		        endif
		    endfor
		endfor

		; define base image header and file name for output FITS files
		mkhdr, hdr, 4, [xsize, ysize]
		basename = 'dark_' + $
		    strtrim(string(dettemps[firstindex], format = '(F5.1)'), 2) + 'K_'

		; write out the mean image
		sxaddpar, hdr, 'IMGTYPE', 'Mean dark current'
		fits_write, outdir + basename + 'mean.fits', *meanimg, hdr

		; write out the standard deviation image
		sxaddpar, hdr, 'IMGTYPE', 'Variance of dark current'
		fits_write, outdir + basename + 'variance.fits', *stdimg, hdr

		; write out the good pixel count image
		sxaddpar, hdr, 'IMGTYPE', 'Good pixel count'
		sxaddpar, hdr, 'NUMIMGS', numimgs
		sxaddpar, hdr, 'BITPIX', 16
		fits_write, outdir + basename + 'pixcount.fits', *countimg, hdr

		; now make histograms for the mean and variance images
		labels = strarr(1)
		labels[0] = 'Detector:  ' + detname
		makehistogram_single, *meanimg, 'Histogram of Mean Dark Current', $
		    'Dark Current (ADU/sec)', labels, $
		    outdir + basename + '_mean_hist.ps', $
		    outdir + basename + '_mean_hist.jpg'

		; free up the memory used for the mean and std images so that
		; more won't have to be reallocated the next time through the
		; loop.
		ptr_free, meanimg, stdimg, countimg, tempcube

		firstindex = lastindex[tempcount] + 1
	endfor

End