;pro dither_median
;
;PURPOSE:
;	To form a single image from a 3x3 set of dithered images.  
;	For each image, the median sky at a given pixel is removed.
;
;INPUTS: 
;	ThisFilter - The letter of the filter that should be sought : 'I', 'G', 'Y'

function Object_Median, Rootdir = Rootdir, Date=Date, ReducedDir = ReducedDir, $
                skipsky=skipsky,skipflat=skipflat, $
                BlockLength = BlockLength, DarkSub = DarkSub, SkySub = SkySub, $
                ObjectName = ObjectName, SkyName = SkyName, DivideByFlat = DivideByFlat, $
                NReads = NReads, ThisFilter = ThisFilter, OverRide = OverRide

 If not(keyword_set(date)) then date='07Apr29'
 If not(keyword_set(RootDir)) then Rootdir='/Volumes/SEA_DISC/'
 If not(keyword_set(ReducedDir)) then ReducedDir = '~/Desktop/LSST/KPNO/Images/Reduced/'+Date+'/'
 ;cd,rootdir+date

 ;Block Length for Offset
 If not(keyword_set(BlockLength)) then BlockLength = 0
 If not(keyword_set(skipsky)) then skipsky=0
 If not(keyword_set(skipflat)) then skipflat=0
 If not(keyword_set(ThisFilter)) then ThisFilter = 'I'
 If not(keyword_set(ObjectName)) then ObjectName = 'M16'
 If not(keyword_set(NReads)) then NReads = 20
 If not(keyword_set(DarkSub)) then DarkSub = 1
 If not(keyword_set(DivideByFlat)) then DivideByFlat = 0
 If not(keyword_set(SkySub)) then SkySub = 1
 If not(keyword_set(SkyName)) then SkyName = 'Cuyo'
 If not(keyword_set(OverRide)) then OverRide = 0

 ;******************************************************************
 ; DARK Subtraction, FLAT fielding and SKY Subtraction
 DarkSubStr = ''
 If DarkSub then DarkSubStr = 'DarkSub'
 ;Get the Median Dark for this number of Reads
 If DarkSub then begin
     DarkMedianFile = file_search(ReducedDir+'MedianDark*'+strtrim(NReads,2)+'Reads*.fits')
     Fits_Read, DarkMedianFile(0), MedianDark , DarkHeader
 Endif

 ;Get the Median Flat for this filter and this number of reads
 If DivideByFlat then begin
    FlatFile = file_search(ReducedDir+'MedianFlat*'+DarkSubStr+'*'+ThisFilter+'*.fits')
    Fits_Read, FlatFile, FlatIm, FlatHeader
    FlatIm = FlatIm/(max(FlatIm) - min(FlatIm))
 EndIf

 ;Get the Median Sky for this filter and this number of Reads
 If SkySub then begin
   SkyFile = File_Search(ReducedDir+'MedianSky_'+DarkSubStr+'*'+ThisFilter+'*'+strtrim(NReads, 2) +'Reads.fits')
   fits_read, SkyFile, Median_Sky, SkyH
 EndIf
 ;*******************************************************************
 ;Get all of the Files for this object and number of reads
 Object_Files = file_search(ObjectName+'*'+strtrim(NReads,2)+'_Reads*.fits')
 dummy=0
 valid_indices = 0

 ;Get a few kewyords
 Fits_Read,  Object_Files(0), NoData, FitsHeader, /header_only       
 Naxis1     = long(sxpar(FitsHeader, 'NAXIS1'))          
 Naxis2     = long(sxpar(FitsHeader, 'NAXIS2'))         
 TotalReads = sxpar(FitsHeader ,'NAXIS3')       
 ArrSize = long(Naxis1*Naxis2)
	
 NumObject_Files = N_elements(Object_Files)
 RAs = Fltarr(NumObject_Files)
 DECs = fltarr(NumObject_Files)
 ;Find images for the proper filter
 for i = 0, n_elements(Object_Files) - 1 do begin
     print, i
     Fits_Read,  Object_Files(i), NoData, FitsHeader, /header_only
     Filter = sxpar(FitsHeader ,'FILTER')
     If stregex(Filter, ThisFilter , /boolean) and $
	not(stregex(Object_Files(i), 'Test', /boolean)) then $
        valid_indices = [valid_indices, i]
 endfor

 ;Get all offsets
 Offset = 167
 Edge = Naxis1 + Offset
 FullDim = Naxis1+2*Offset

 FullIms = lonArr(FullDim, FullDim, N_elements(valid_indices)-1 )
 RAPix = [0, Offset, Offset, 0, -Offset, -Offset, -Offset, 0, Offset]
 DecPix = [0, 0, -Offset, -Offset, -Offset, 0, Offset, Offset, Offset]

 ;*******************************************************************
 ; Form Dithered Image with Desired Filter
 For i= 1, N_elements(valid_indices) - 1 do begin
     index = valid_indices(i)
     Fits_Read,  Object_Files(index), Data, FitsHeader, /header_only 
     RAs[index] = float(sxpar(FitsHeader, 'RA'))  
     DECs[index] = float(sxpar(FitsHeader, 'DEC'))  
     Filter = sxpar(FitsHeader, 'FILTER')
     print, index, ' ', Object_Files(index), ' ' ,Filter

     ;Get the difference Image
     DiffIm = return_cube_diff(Object_Files(index), ArrSize, NAxis1, Naxis2, 1, TotalReads-1, $
	BlockLength = BlockLength)

     ;Do the Dark Subtraction, Flat Field and then SkySubtraction if specified in keywords
     If DarkSub then DiffIm = DiffIm - MedianDark
     If DivideByFlat then DiffIm = DiffIm/FlatIm
     If SkySub then DiffIm = DiffIm - SkySub

     ;Put the Image in the right spot
     FullIms(Offset + RaPix(valid_indices(i)) : Offset + RaPix(valid_indices(i)) + Naxis1 - 1, $
		    Offset + DecPix(valid_indices(i)): Offset + DecPix(valid_indices(i)) + Naxis2 - 1, $
		    i - 1) = DiffIm

  endfor

  medim = median(FullIms, dimension = 3)
  tvscl,congrid(medim,512,512)
  writefits, ReducedDir + ObjectName+'_'+ThisFilter+'_'+'Dithered.fits',medim
  print, 'Wrote Dithered Image to :' + ReducedDir+ObjectName+'_'+ThisFilter+'_'+'Dithered.fits'
  stop


end