pro rectify_raw_irac_data, rfiles, out_savefile, MAKEFITS=makefits, VERBOSE=verbose, $
                           BCD_ORIENTATION=bcd_orientation
;+
; NAME:
;    RECTIFY_RAW_IRAC_DATA
;
; PURPOSE:
;    Turns IRAC subarray raw file into data cube with rectified and corrected DN values in
;    floating point format.  The procedure will also rectify full array data but will not 
;    create a useful cube stack.
;
; EXAMPLES:
;    To rectify a set of channel 1 raw files and make corresponding FITS files
;    IDL> files = file_search("SPITZER_I1_*_dce.fits")
;    IDL> rectify_raw_irac_data, files, /MAKEFITS
;
;    To rectify a set of channel 2 raw files, make corresponding FITS files in BCD orientation
;    and write an IDL save file containing the data cubes    
;    IDL> files = file_search("SPITZER_I1_*_dce.fits")
;    IDL> rectify_raw_irac_data, files, 'ch2_cubes.sav', /MAKEFITS, /BCD_ORIENTATION
;
; INPUT:
;    rfiles: string array containing raw data files to rectify
;
; OPTIONAL OUTPUT:
;    ofile: string containing filename to write rectified cubes to in IDL .sav format
;           The save file will consist of three arrays:
;                stack which is a 256x256xN array consisting of all the rectified images in 2d 
;                 format and N is the number of images.  This array is most useful for full 
;                 array data.
;                cube_stack which is a 32x32x64xN array of all the rectified subarray frames, but
;                  a 32x32x64 slice will be zero for any full array images
;                mstack which is an array of dimension N which contains the maximum pixel value in
;                  each image for full array data and each subframe for subarray data.
; 
; OPTIONAL OUTPUT KEYWORDS:
;    MAKEFITS: if set, writes out rectified cubes in input data directory using '_recdce.fits' as
;              filename suffix, instead of '_dce.fits'
;    VERBOSE: Print out some diagnostics
;    BCD_ORIENTATION: If set, raw images will be transformed into BCD orientation.  That is, the
;                     channel 1 and 2 images will be flipped about the x-axis (y -> -y)
;
; PROCEDURES USED:
;	 WRITEFITS, SXADDPAR, SXADDHIST, REMOVE_LEADING_PATH, SJC_RENAME, PRINT_PCT
;
; FUNCTIONS USED:
;    READFITS(), SXPAR(), MAX()
;
; EXTERNAL LIBRARIES NEEDED:
;    astrolib
;   
; METHOD:
;    Following the pipeline description document, each input file is converted from unsigned
;    integers (InSb data are flipped about zero), adjusted for Fowler bias and corrected
;    for wrap around.  If the data is subarray, it is transformed into cubes.  
;
; HISTORY
;    Version 1.0 released to IRSA 15 Aug 2011 SJC
;    Full array bug fix, increased documentation, added BCD_ORIENTATION keyword 15 Aug 2011 SJC
;    Module comments added and some code checking done 11 August 2011 SJC
;-
	if (N_params() lt 1 or N_params() gt 2) then begin
		print, 'Syntax: RECTIFY_RAW_IRAC_DATA, rawfiles, output_cube_file, /MAKEFITS, /VERBOSE, $'
		print, '                                /BCD_ORIENTATION'
		return
	endif

; Figure out number of input files
	n = n_elements(rfiles)

; If output file is given, then prepare output stacks	
	if (N_params() lt 2) then do_stacks = 0 else do_stacks = 1
; Big arrays of all cubes/frames suitable for further processing	
	if (do_stacks) then begin
; stack of cubes	
		cube_stack = fltarr(32, 32, 64, n)
; stack of 2d frames
		stack = fltarr(256, 256, n)
; Maximum pixel in each subframe	
		mstack = fltarr(n*64L)		
	endif 

; Loop over all input files
	if (keyword_set(VERBOSE)) then print, FORMAT='($, "RECTIFY_RAW_IRAC_DATA: Processing")'
	mstack_count = 0L
	for i = 0, n-1 do begin
		if (keyword_set(VERBOSE)) then print_pct, i, n
		im = readfits(rfiles[i], h, /SILENT)
; Get header keywords pertaining to data sampling
; Barrel shift
		bs = sxpar(h, 'A0741D00')
; Fowler Number		
		fn = sxpar(h, 'A0614D00')
; Wait Ticks		
		wt = sxpar(h, 'A0615D00')
; Readmode 1 for subarray, 0 for full		
		readmod = sxpar(h, 'A0617D00')
; Channel number 
		ch = sxpar(h, 'A0610D00')

; Flip sign	for InSb dectectors	
		if (ch lt 3) then nim = 65535 - uint(im)

; Finally convert to floating point
		nim = float(nim)
; Handle bias due to Fowler sampling		
		nim = nim - 0.5 * (1. - 2.^(-bs))
; Treat wraparound correctly
		ptr = where(nim gt 45000., count)
		if (count gt 0) then nim[ptr] = nim[ptr] - 65535.
; Handle barrel shift (this is a NOOP for standard IRAC data, for now anyways)		
		nim = nim * 2.^bs / fn
; Place in 2d stack if warranted
		if (do_stacks) then stack[*, *, i] = nim
		
; turn into cube; If subarray array mode, process file skip this file
		if (readmod eq 1) then begin
			cube = float(reform(nim, 32, 32, 64))
; Concatenate cubes into one 4d array		
			if (do_stacks) then cube_stack[*, *, *, i] = cube
		endif
		
; Write rectified data to FITS files if desired
		if (keyword_set(MAKEFITS)) then begin
			remove_leading_path, rfiles[i], base, LPATH=dir
			sjc_rename, base, '_dce', '_recdce', ofile
			nh = h
; Adjust BITPIX as data is now floating point
			sxaddpar, nh, 'BITPIX', -32
; If subarray cube then change the NAXIS keyowrds appropriately for the 3d format
			if (readmod eq 1) then begin			
				sxaddpar, nh, 'NAXIS', 3
				sxaddpar, nh, 'NAXIS1', 32
				sxaddpar, nh, 'NAXIS2', 32
				sxaddpar, nh, 'NAXIS3', 64
; Place in BCD orientation if desired, only Channels 1 and 2 need to be flipped
				if (keyword_set(BCD_ORIENTATION) and ch lt 3) then begin
					for k = 0, 63 do begin
						sub = rotate(cube[*, *, k], 7)
						cube[*, *, k] = sub
					endfor
					sxaddhist, 'Flipped y-axis so image in BCD coordinates', nh
				endif
				writefits, dir + ofile, cube, nh
			endif else begin
; Place full array image in BCD orientation if desired			
				if (keyword_set(BCD_ORIENTATION) and ch lt 3) then begin
					nim = rotate(nim, 7)
					sxaddhist, 'Flipped y-axis so image in BCD coordinates', nh
				endif
				writefits, dir + ofile, nim, nh
			endelse 
		endif
; Put max value in output array	for each image (full array) or each subimage (subarray)	
		if (do_stacks) then begin
			if (readmod eq 1) then begin
; Is subarray so loop over the 64 subframes in the cube
				for k = 0L, 63L do begin
					mstack[mstack_count] = max(cube[*, *, k], /NAN)
; Increment to next subframe
					mstack_count = mstack_count + 1L
				endfor
; For full array just find maximum of image
			endif else begin
				mstack[mstack_count] = max(nim, /NAN)
; Increment counter to next DCE
				mstack_count = mstack_count + 1L
			endelse
		endif
	endfor
	

	if (do_stacks) then begin
; Truncate the maximum stack for full array mode, should be a NOOP for subarray	
		mstack = mstack[0L:mstack_count-1L]
; Output stack.sav file
		save, FILENAME=out_savefile, stack, cube_stack, mstack
	endif
return
end

pro remove_leading_path, istr, ostr, LPATH=lpath
;+
; NAME:
;    REMOVE_LEADING_PATH
;
; PURPOSE:
;    Removes the leading directories from an input string
; 
; INPUT:
;    istr: string containing input pathname
;
; OUTPUT:
;    ostr: string containing output filename
; 
; OUTPUT Keyword
;    lpath: if set, returns the leading path of the input string
;
; METHOD:
;    Find the last '/' in the string and remove all characters from the
;    first to that position.  For lpath, returns from the first character to
;    the last '/'.  If there are no '/', then lpath is a null string.
;
; HISTORY
;    Initial coding SJC 29 Dec 2006
;    Module checked SJC 29 Dec 2006
;    Added leading path option 23 April 2008
;-

; Check input parameters
	if (N_params() ne 2) then begin
		print, 'Syntax -- remove_leading_path, in_string, out_string, $'
		print, '              LPATH=leading_path'
		return
	endif

; Set slash string
    slash_str = '/'

; Copy the input string to the output string
    ostr = istr

; Find the position of the last slash in the input string
    pos = strpos(ostr, slash_str, /REVERSE_SEARCH)

; If the lash slash is not the last character then add a slash to the end of 
; the string
    if (pos ne -1) then begin
        lpath = strmid(ostr, 0, pos+1)
    	ostr = strmid(ostr, pos+1, strlen(ostr)-pos-1) 
    endif else lpath = ''
    
return
end

pro sjc_rename, istr, oldsub, newsub, nstr
;+
; NAME:
;    SJC_RENAME
; 
; PURPOSE:
;    To change a filename by replacing a substring with a new substring
;
; INPUT:
;    istr: string containing input filename
;    oldsub: string containing substring to change
;    newsub: string containing substring to change to
;
; OUTPUT: 
;    nstr: string containing output filename
;    
; HISTORY
;    Initial coding SJC 15 Nov 2006
;    Module checked SJC 27 Dec 2006
;-
    nosub = strlen(oldsub)
    nstr = istr

    pos = strpos(nstr, oldsub)
    while (pos gt -1) do begin
        nnstr = strlen(nstr)
        nstr = strmid(nstr, 0, pos) + newsub + $
                        strmid(nstr, pos+nosub, nnstr-pos+nosub)
        pos = strpos(nstr, oldsub)
    endwhile

return
end

pro print_pct, count, nmax, pcts=pcts
;+
; NAME:
;    PRINT_PCT
; 
; PURPOSE:
;    To print progress stepping through a set of values
;
; INPUT:
;    count: integer containing current step
;    nmax: maximum number of iterations
;
; KEYWORD: 
;    PCTS: if set, determines the percentages at which to print out a value
;
; NOTES:
;    Algorithm does not produce great results for small maximum number of iterations.  Also
;    procedure is designed for steps that end on the maximum value.
;    
; HISTORY
;    Comments added SJC 11 Aug 2011
;-
; Default to printing a step ever 10%
	if (not keyword_set(PCTS)) then pcts = 10
; Determine indices at which to print a value	
	indices = (indgen(fix(nmax/pcts)>1)+1) * fix(pcts)
	vals = long(float(indices) / 100.0 * nmax)
; See if current value is one of those indices
	ptr = where(count eq vals, n)
; If so, then print a step
	if (n eq 1) then print,FORMAT='($,"...",I3,"%")', indices[ptr[0]]
; At end, print completion
	if (count ge nmax-1) then begin 
		print, FORMAT='($,"...100%")'
		print, ''
	endif

return
end