;pro 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

        if not(keyword_set(date)) then date='07Apr29'
        rootdir='/Volumes/SEA_DISC/'
        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 = 'G'
	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

                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)
	
		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)

                   dummy=0

                endfor

                medim=median(im,dimension=3)
                writefits, Reduceddir+'medsky_'+ThisFilter+'_'+date+'.fits',medim

        endif else begin

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

	return, medim

stop

end