;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.
;
;KEYWORDS:
;	Keywords taken from KeywordStruct.pro file in same directory as this file

pro Dither_Median, rootdir=rootdir, date=date,skipsky=skipsky,skipflat=skipflat, $
		   reduceddir = reduceddir, ObjectName=ObjectName, NReads = NReads,$
		   ThisFilter = ThisFilter

  ;Get defaults for keywords
  @KeywordStruct.pro

  ;Get all of the that match the input parameters
  ObjectFiles = file_search(ObjectName+'*'+strtrim(NReads,2)+'_Reads*.fits')
  ObjectFiles = Return_File_List(ObjectFiles, ThisFilter = ThisFilter)

  dummy=0
  valid_indices = 0

  ;Get a few kewyords
  Fits_Read, ObjectFiles(0), NoData, FitsHeader, /header_only       
  Naxis1     = long(sxpar(FitsHeader, 'NAXIS1'))          
  Naxis2     = long(sxpar(FitsHeader, 'NAXIS2'))         
  TotalReads = sxpar(FitsHeader ,'NAXIS3')       
  ArrSize    = long(Naxis1*Naxis2)

 	
  NumObjectFiles = N_elements(ObjectFiles)
  RAs = Fltarr(NumObjectFiles)
  DECs = fltarr(NumObjectFiles)

  ObjectFiles
        ;Find images for the proper filter
        for i = 0, n_elements(ObjectFiles) - 1 do begin
            print, i
            Fits_Read,  ObjectFiles(i), NoData, FitsHeader, /header_only
            Filter = sxpar(FitsHeader ,'FILTER')
            If stregex(Filter, ThisFilter , /boolean) then $
            valid_indices = [valid_indices, i]
        endfor

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

        FullIms = IntArr(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]

	;Get Sky_Frame for subtraction
	fits_read, reduceddir+'medsky_'+ ThisFilter + '_'+Date+'.fits', Sky_Frame, SkyH
	
	; Form Dithered Image with Desired Filter
	For i= 1, N_elements(valid_indices) - 1 do begin
	    print, valid_indices(i)
            Fits_Read,  ObjectFiles(valid_indices(i)), Data, FitsHeader, /header_only 
	    RAs[valid_indices(i)] = float(sxpar(FitsHeader, 'RA'))  
            DECs[valid_indices(i)] = float(sxpar(FitsHeader, 'DEC'))  

           ;I Filter For Object
           Fits_Read, ObjectFiles(valid_indices(i)), Frame_First, $ 
		First = BlockLength, $
 		Last = long(ArrSize) + BlockLength - 1
           Frame_First = reform(temporary(Frame_First + 32768.), Naxis1, Naxis2)
           Fits_Read, ObjectFiles(valid_indices(i)), Frame_Last, $
		First = (TotalReads - 1) * ArrSize + BlockLength, $
                Last = (TotalReads)*ArrSize +BlockLength - 1
           Frame_Last = reform(temporary(Frame_Last + 32768.), Naxis1, Naxis2)

   	   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, $
		     valid_indices(i)) = long(Frame_Last) - long(Frame_First) - Sky_Frame
	endfor
	
	medim = median(FullIms, dimension = 3)
        writefits,reduceddir+ObjectName+'_'+ThisFilter+'_'+'Dithered.fits',medim

	stop


end