;function Median_Sky
;
;PURPOSE:
;	To produce a median value for a pixel that consists of sky background
;	in a given filter
;
function Sky_Median, Rootdir = Rootdir, Date=Date, ReducedDir = ReducedDir, $
		skipsky=skipsky,skipflat=skipflat, $
		BlockLength = BlockLength, DarkSub = DarkSub, $
		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 = 'Dark'
 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(SkyName)) then SkyName = 'Cuyo'
 If not(keyword_set(OverRide)) then OverRide = 0

 ;***************************************************************************************
 ;Check to see if the Median Dark already exists and only write a fits file if it doesn't
 DarkSubStr = ''
 If DarkSub then DarkSubStr = 'DarkSub' 
 SkyMedianFile = file_search(ReducedDir+'MedianSky*'+DarkSubStr+'*'+$
		   ThisFilter+'*'+strtrim(NReads,2)+'Reads*.fits')
 
 If DivideByFlat then begin
    FlatFile = file_search(ReducedDir+'MedianFlat*'+DarkSubStr+'*'+ThisFilter+'*.fits')
    Fits_Read, FlatFile, FlatIm, FlatHeader
    FlatIm = FlatIm/(max(FlatIm) - min(FlatIm))
 EndIf

 If SkyMedianFile eq '' or OverRide then begin
 
   print, 'Forming Median Sky Image for ' + strtrim(NReads,2) + ' Read Ramp from ' + SkyName

   ;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
   ;Set up Indices for the images in the filter we want
   If SkyName eq 'Cuyo' and ThisFilter eq 'I' then NReads = 20
   Files=file_search(SkyName+'*Dither*'+strtrim(NReads,2)+'_Reads*.fits')
   dummy=0
   valid_indices = 0

   ;Find Images with the right filters
   for i = 0, n_elements(Files) - 1 do begin
      Fits_Read,  Files(i), NoData, FitsHeader, /header_only
      Filter = sxpar(FitsHeader ,'FILTER')
      If stregex(Filter, ThisFilter , /boolean) then $
	     valid_indices = [valid_indices, i]
   endfor

   ;Get a Few Keywords
   Fits_Read,  Files(valid_indices(0)), NoData, FitsHeader, /header_only
   Naxis1     = long(sxpar(FitsHeader, 'NAXIS1'))
   Naxis2     = long(sxpar(FitsHeader, 'NAXIS2'))
   TotalReads = sxpar(FitsHeader ,'NAXIS3')
   ArrSize = long(Naxis1*Naxis2)
	
   ;Form a datacube big enough to hold the array
   Im=lonarr(Naxis1,Naxis2,n_elements(valid_indices-1))

   for i=1, n_elements(valid_indices)-1 do begin
       index=valid_indices(i)
       Fits_Read,  Files(index), NoData, FitsHeader, /header_only
       Filter = sxpar(FitsHeader ,'FILTER')
       print, i,' ', Files(index), ' ' , Filter

       ;Get differnece
       If SkyName eq 'Cuyo' and ThisFilter eq 'I' then TotalReads = 15
       Im(*,*,i-1)= return_cube_diff(Files(index), ArrSize, Naxis1, Naxis2, 1, TotalReads) 
       If DarkSub then im(*,*,i-1) = im(*,*,i-1) - long(MedianDark)
       If DivideByFlat then begin
          im(*,*,i-1)= im(*,*,i-1)/float(FlatIm)
       endif
       delvarx, dummy

   endfor

   medim=median(Im,dimension=3)
   writefits, Reduceddir+'MedianSky_'+DarkSubStr+'_'+ThisFilter+'_'+date+'_'+$
	strtrim(TotalReads, 2)+'Reads.fits', medim
   print, 'Wrote Median Sky Image: ' + 'MedianSky_'+DarkSubStr+'_'+$
	ThisFilter+'_'+date+'_'+strtrim(TotalReads, 2)+'Reads.fits'
   return, 1

 EndIf else begin

    return, 2

 Endelse
 
stop

end