;function Median_Sky
;
;PURPOSE:
;	To produce a median value for a pixel that consists of sky background
;	in a given filter
;

function Median_Sky, date=date,skipsky=skipsky,skipflat=skipflat, Thisfilter = Thisfilter, $
		BlockLength = BlockLength, DarkSub = DarkSub

        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/'
        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

	;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

	skipflat = 1
	skipsky = 0
        if skipsky eq 0 then begin

		;Set up Indices for the images in the filter we want
                files=file_search('Cuyo*Dither*20_Reads*.fits')
		dummy=0
		valid_indices = 0
		
		;Find I filter images
		for i = 0, n_elements(files) - 1 do begin
		   print, i
                   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=intarr(Naxis1,Naxis2,n_elements(valid_indices-1))

                for i=1, n_elements(valid_indices)-1 do begin
		   print, i
		   index=valid_indices(i)
		   Fits_Read, Files(index), Dummy_Last, First = (TotalReads - 2) * ArrSize, $
                                       Last = (TotalReads-1)*ArrSize - 1
                   Dummy_Last = reform(temporary(Dummy_Last), Naxis1, Naxis2)

       		   Fits_Read, Files(index), Dummy_First, First = 0, $
					Last = long(ArrSize) - 1
		   Dummy_First = reform(temporary(Dummy_First), Naxis1, Naxis2)
 
                   im(*,*,i)= temporary(Dummy_Last - Dummy_First)
		   if DarkSub then im(*,*,i) = im(*,*,i) - MedianDark
                   dummy=0

                endfor

                medim=median(im,dimension=3)
		if DarkSub then begin
                   writefits, Reduceddir+'MedianSky_DarkSub_'+ThisFilter+'_'+date+'.fits',medim
		Endif Else begin
                   writefits, Reduceddir+'MedianSky_'+ThisFilter+'_'+date+'.fits',medim
		EndElse

        endif else begin

                medim=readfits('medsky'+date+'.fits',hsky)

        endelse
 
stop

end