;pro median_difference_of_reads
;
;PURPOSE:
;	To take a collection of images of a given type, i.e. a combination of 
;	a given number of reads and object type (Flat, Dark, M16,...) and get 
;	the median value of the pixels and then write them to a fits file.
;
;	The resulting file will be an Naxis1xNaxis2 image since there is 
;	only a difference of two reads and the median of this difference is 
;	taken across multiple files.
;
;INPUTS:
;
;KEYWORDS:
;	ROOTDIR: 
;	  The Directory holding the directories of Raw files
;	DATE:
;	  The Date during which they were taken.  Also the subdirectory that 
;	  the files are in.
;	REDUCEDDIR:
;	  The directory where the processed images and plots will be placed in.
;	BLOCKLENGTH:
;	  An offset for files that were written improperly.  The readu 
;	  function reads the values with this offset in place.
;	OBJECTNAME:
;	  The string that indicates what files are being sought.  The filenames
;	  should begin with this string ('Dark', 'Flat', etc.).
;	NREADS:
;	  The number of reads (NAXIS3) that are present in the file.  This
;	  is also a string in the filename that will be used as a regular
;	  expression in searching for matched files.
;	FIRSTREAD:
;	  The first read to be used in the difference 
;	      Diff = LastRead - FirstRead
;	LASTREAD: 
;	  The last read to be used in the difference
;	THISFILTER:
;	  The filter through which the files were taken. ('I', 'G', or 'Y')
;OUTPUTS:
;
function Median_Difference_Of_Reads, RootDir=RootDir, Date=Date, $
	 ReducedDir = ReducedDir, BlockLength = BlockLength, $ 
	 ObjectName=ObjectName, NReads = NReads, ThisFilter = ThisFilter,$
	 OverRide = OverRide, FirstRead = FirstRead, LastRead = LastRead

  ;Directory and FileName Keywords
  If Not keyword_set(RootDir) then RootDir = '/nfs/slac/g/ki/ki09/lances/'
  If Not keyword_set(Date) then Date='07Apr30'
  If Not keyword_set(ReducedDir) then ReducedDir= RootDir+'/Reduced/'+Date+'/'
  If Not keyword_set(BlockLength) then BlockLength = 0
  If Not keyword_set(ObjectName) then ObjectName = 'Dark'
  If Not keyword_set(NReads) then NReads = 20
  If Not keyword_set(FirstRead) then FirstRead = 0
  If Not keyword_set(LastRead) then LastRead = NReads - 1
  If Not keyword_set(ThisFilter) then ThisFilter = 'Blank'
  If Not keyword_set(OverRide) then OverRide = 0
  cd, Rootdir+Date
	
  ;****************************************************************************
  ;Check to see if the Median Dark already exists and only write if it doesn't
  MedianFile = file_search(ReducedDir+'Median'+ObjectName+'_'+$
			   ThisFilter+'*'+strtrim(NReads,2)+'Reads*.fits')

  ;If Dark does not exist, form the median and write the image
  If MedianFile(0) eq '' or OverRide then begin
    print, 'Forming Median Image for Combination of:' 
    print, 'Object: ' + ObjectName
    print, 'Filter: ' + ThisFilter
    print, '# Reads:' +  strtrim(NReads,2)  

    ;Get all of the Dark Files for this Number of Reads
    RawFiles = File_Search(ObjectName+'*'+strtrim(NReads,2)+'_Reads*.fits')
    Dummy = 0 
    ValidIndices = 0

    ;Get a Few Kewyords
    Fits_Read,  RawFiles(0), NoData, FitsHeader, /header_only       
    Naxis1     = long(sxpar(FitsHeader, 'NAXIS1'))          
    Naxis2     = long(sxpar(FitsHeader, 'NAXIS2'))         
    TotalReads = sxpar(FitsHeader,'NAXIS3')      
    ArrSize = long(Naxis1*Naxis2)
	
    ;Only use the images that are certain to have Blank as the filter
    for i = 0, n_elements(RawFiles) - 1 do begin
        Fits_Read, RawFiles(i), NoData, FitsHeader, /header_only
        Filter = SxPar(FitsHeader ,'FILTER')
	print, i, ' ' , Filter
        If stregex(Filter, ThisFilter , /boolean) and $
	   not stregex(RawFiles(i), '42_44', /boolean) then $
           ValidIndices = [ValidIndices, i]
        endfor

	;Form an Array big enough for all of the darks of a given Naxis3
	NumRawFiles = N_elements(ValidIndices)
	RawIms = LonArr(Naxis1, Naxis2, NumRawFiles-1)
	
	;Form Cube of Darks and then Take Median at the End 
	for FileNum = 1, NumRawFiles - 1 do begin
	   Index = ValidIndices(FileNum)
	   Fits_Read,  RawFiles(Index), NoData, FitsHeader, /header_only
           Filter = SxPar(FitsHeader ,'FILTER')
	   print, Index, ' ', Filter, ' ', RawFiles(Index) 

           ;Get Difference of Dark Frames
	   RawIms(*,*, FileNum - 1) = Return_Cube_Diff(RawFiles(Index), $
		  ArrSize, Naxis1, Naxis2, FirstRead, LastRead)

	endfor
	
	;Write a fits file with the Median Dark
	MedIm = median(RawIms, dimension = 3)
        WriteFits, ReducedDir+'Median'+ObjectName+'_'+ThisFilter+$
		Strtrim(NReads,2)+'Reads_'+strtrim(LastRead,2)+'_Minus'+$
                strtrim(FirstRead,2)+'.fits', Medim
	print, 'Median Dark Image written with filename: '
	print,  ReducedDir+'Median'+ObjectName+'_'+ThisFilter+'_'+$
                Strtrim(NReads,2)+'Reads_'+strtrim(LastRead,2)+'_Minus'+$
		strtrim(FirstRead,2)+'.fits'

	return, 1

  EndIf else begin
        
	;The File already exists
	print, 'Median Dark Image: ' +  ReducedDir+'Median'+ObjectName+'_'+$
		Strtrim(NReads,2)+'Reads_'+strtrim(LastRead,2)+'_Minus'+$
                strtrim(FirstRead,2)+'.fits'+$
		' already exists...continuing'
	return, 2
 
  EndElse

stop

end