pro repeat_delta_dark, files1, files2, ofiles
;+
; NAME:
;   REPEAT_DELTA_DARK
; PURPOSE:
;   To subtract a dark offset from repeat frames in an AOR.  The
;   "first-frame" correction for repeats in an IRAC AOR is not optimal as
;   the correction is determined from an extrapolation to short delay (<1s)
;   times.  In addition, the variation in "first-frame" effect is most
;   significant for short delay times.  The "first-frame" effect is largest
;   in channel 3 data.
;
; SYNTAX:
;   repeat_delta_dark, files1, files2, ofiles
;
; INPUTS:
;   files1: Array of filenames of 1st repeats in an AOR.  If a full-frame AOR 
;           uses 2 repeats, this is every 2nd BCD starting with the first BCD 
;           of that frametime.
;   files2: Array of filenames of a particular repeat instance.  There should
;           be an equal number of files in files1 and files2 
;   ofiles: string array containing names of files to write delta-dark
;           corrected images to.
;
; USAGE EXAMPLE:
;   Consider a 100s full frame AOR that consists of a 3 by 3 map with 2
;   cycling dithers and 3 repeats.  For channel 3, there will be 
;   3x3x2x3 + 2 = 56 BCDs.  Exposure numbers 0 and 1 are the short frametime
;   throwaway observations at the beginning of the AOR.  BCDs 2-4 will be
;   the set of 3 repeats at map position 1, first dither position.  BCDs
;   5-7 will be the set of 3 repeats at map position 1, second dither.  The
;   subsequent map/dither positions will each have 3 repeats (8-10, 11-13, etc.)
;   The second and third repeats at each position will need to have a delta
;   dark applied.  To fix the repeats, make an array containing the
;   filenames of the first repeats (exposure numbers 2, 5, 8, 11, etc.) as 
;   files1.
;   IDL> files1 = ['SPITZER_I1_9999999_0002_0000_1_bcd.fits',  $
;   IDL>           'SPITZER_I1_9999999_0005_0000_1_bcd.fits',  $
;   IDL>           'SPITZER_I1_9999999_0008_0000_1_bcd.fits',  $
;   ...
;   IDL>           'SPITZER_I1_9999999_0053_0000_1_bcd.fits']
;   The second repeats should go in array files2:
;   IDL> files2 = ['SPITZER_I1_9999999_0003_0000_1_bcd.fits',  $
;   IDL>           'SPITZER_I1_9999999_0006_0000_1_bcd.fits',  $
;   IDL>           'SPITZER_I1_9999999_0009_0000_1_bcd.fits',  $
;   ...
;   IDL>           'SPITZER_I1_9999999_0054_0000_1_bcd.fits']
;   The corrected file names can have "dd_" appended if you like.  If you
;   use the same filenames, the originals will be overwritten.
;   IDL> ofiles = 'dd_' + files2
;   Now run the procedure
;   IDL> repeat_delta_dark, files1, files2, ofiles
; 
;   To correct the 3rd repeat, make a list of those repeats
;   IDL> files3 = ['SPITZER_I1_9999999_0004_0000_1_bcd.fits',  $
;   IDL>           'SPITZER_I1_9999999_0007_0000_1_bcd.fits',  $
;   IDL>           'SPITZER_I1_9999999_0010_0000_1_bcd.fits',  $
;   ...
;   IDL>           'SPITZER_I1_9999999_0055_0000_1_bcd.fits']
;   IDL> ofiles3 = 'dd_' + files3
;
;   If you want to save the deltadark applied to the 2nd set of repeats,
;   then copy the file 'deltadark.fits' to another filename before rerunning the
;   procedure on the 3rd repeat.
;   IDL> repeat_delta_dark, files1, files3, ofiles3
;
; METHOD:
;   The procedure takes the difference between the 1st repeat and the repeat
;   to be corrected for each pair of repeats as specified in files1 and files2.
;   Each difference image is an estimate of the residual bias offset in the 
;   repeats contained in files2.  The median of the image difference stack 
;   is calculated on a per pixel basis.  This median is the delta-dark
;   correction which is then applied to each image in files2 as a correction. 
;   The corrected images have their headers updated and are written to ofiles. 
;   In addition, the delta-dark used is written to a file deltadark.fits
;
; NOTES:
;   The residual first-frame effect is most significant in channel 3 and this
;   correction should not have to be applied to data in channels 1, 2 and 4.
;
; HISTORY
;   Code written 05 August 2005 SJC
;   Comments added 30 August 2005 and code released SJC
;-

; Check inputs
	if (N_params() ne 3) then begin
		print, 'Syntax -- REPEAT_DELTA_DARK, exposures1, exposures2, $'
		print, '                             out_exposures2'
		return
	endif

; Count number of files in each list
	n1 = n_elements(files1)
	n2 = n_elements(files2)

; If number of first repeats does not match the number of the 2nd, ...
; repeat, then exit
	if (n1 ne n2) then begin
		print, 'REPEAT_DELTA_DARK number of exposure 1 and 2 files must be the same'
		return
	endif

; Get image size
	h0 = headfits(files1[0])
	nx = sxpar(h0, 'NAXIS1')
	ny = sxpar(h0, 'NAXIS2')

; Create image stack for offset calculation
	stack = fltarr(nx, ny, n1)
; Create output image stack
	outimages = fltarr(nx, ny, n1)

; Difference exposures
	for i = 0, n1-1 do begin
		im1 = readfits(files1[i], h)
		im2 = readfits(files2[i], h)
		stack[*, *, i] = im1 - im2
		outimages[*, *, i] = im2
	endfor

; Use pixelwise median of stack as delta dark
; A more robust estimator could be used, but this should 
; suffice for most observations
	delta_dark = median(stack, DIMENSION=3)
	hdark = h0
; Output delta dark used
	sxaddpar, hdark, 'DELTADARK', 'T', 'Image is delta dark'
	writefits, 'deltadark.fits', delta_dark, hdark

; Make output images
	for i = 0, n2-1 do begin
; Apply delta dark to each repeated image
		outimages[*, *, i] = outimages[*, *, i] + delta_dark
		h = headfits(files2[i])
; See if the file has already been corrected. If so, then do not correct it
; further
		ddproc = sxpar(h, 'DDPROC', COUNT=count)
		if (count eq 0) then begin
; Update header and write out file
			sxaddpar, h, 'DDPROC', 'T', ' Delta dark subtracted'
			writefits, ofiles[i], outimages[*, *, i], h
		endif else begin
			print, 'REPEAT_DELTA_DARK: Already corrected ' + files2[i]
		endelse
	endfor
return
end