;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 
;
;	WAITPERSIST:
;	  The number of dithers to wait until the pixels that were saturated
;	  are used in the mean or median of the mosaic.
;	REJECTPERSIST:
;	  If this keyword is set, the pixels that have been saturated in 
;	  a frame will be set to zero for WAITPERSIST frames afterwards.
;	RADMAX:
;	  The radius (in pixels) from the center of a star for which the 
;	  pixels will be zeroed if REJECTPERSIST is set.

;Star for M57: 3100, 3560
;CALLING SEQUENCE:
;
; dither_slope_median,'/nfs/slac/g/ki/ki09/lances/Reduced/07Apr29/M16/M16GSlopes.lst',/rejectpersist
pro Dither_Slope_Median, FileList, XStart = XStart, YStart=YStart, $
		   RootDir=RootDir, Date=Date, $
		   ReducedDir = ReducedDir, $
		   SkySubtract= SkySubtract, ObjectName=ObjectName, $
		   NReads = NReads, ThisFilter = ThisFilter, $
		   FlatFielded = FlatFielded, SkipFileNum = SkipFileNum,$
		   TvFlag = TvFlag, DitherCheck = DitherCheck, $
		   SkipStar = SkipStar, Rem60Hz = Rem60Hz, $
		   RejectPersist = RejectPersist, WaitPersist = WaitPersist, $
		   Verbose = Verbose, RadMax = RadMax

  If N_Elements(WaitPersist) eq 0 then WaitPersist = 3
  ;Make some common blocks
  Common KeyParams, KeyStr
  Common ColorStr, Color, ColorNum, ColorTri, ColorNumTri

  ;Get defaults for keywords
  @KeywordStruct.pro
  @PlotSettings
  cd, ReducedDir

  If N_Elements(SkipFileNum) eq 0 then SkipFileNum = 0
  If FlatFielded then FlatFieldStr = 'FlatFielded' else FlatFieldStr = ''
  If N_Elements(SkipStar)    eq 0 then SkipStar = 0
  If N_Elements(DitherCheck) eq 0 then DitherCheck = 0
  If N_Elements(RadMax)      eq 0 then RadMax = 50

  ;Check to see whether persistence rejection is going to happen
  If KeyStr.RejectPersist then begin
     KeyStr.RejectPersistStr = 'PersistReject'+Strtrim(WaitPersist,2)+'_'
  EndIf

  ;Get the files from a list in a file and the offsets
  ReadCol, FileList, RawFiles, XOff, YOff, Format = 'A,I,I'
  ;Put the info in the structure
  KeyStr.DitherFiles(0:N_Elements(RawFiles)-1) = RawFiles
  KeyStr.NumDithers = N_Elements(RawFiles)

  ;Get the bright star coordinates from the same file
  ReadCol, FileList, XStart, YStart, Format = 'I,I', $
    SkipLine = N_Elements(RawFiles)
  For i = 0, N_Elements(XStart)-1 do begin 
    KeyStr.RAs(0:N_Elements(RawFiles)-1,i)  = XOff
    KeyStr.DECs(0:N_Elements(RawFiles)-1,i) = YOff
  EndFor

  If KeyStr.TvFlag then window, xsize=700,ysize=700
  ;Get a few kewyords, pointings, and offsets
  FileCheck = Return_File_Params(RawFiles(0), /slope)	

  ;Try out a first guess at the offsets from our 20 arcsecond dithers
  Return_Dither_Offsets, XStart, YStart, TvFlag=TvFlag, SkipStar = SkipStar
  NumRawFiles = N_elements(RawFiles)

  ;Find the Maximums for the dimensions of the array
  MinXOff  = Min(KeyStr.XOffs)  & MaxXOff  = Max(KeyStr.XOffs)
  MinYOff  = Min(KeyStr.YOffs)  & MaxYOff = Max(KeyStr.YOffs)
  NewArrXSize = KeyStr.Naxis1 + (MaxXOff - MinXOff)
  NewArrYSize = KeyStr.Naxis2 + (MaxYOff - MinYOff)

  ;Form the Array for Airmasses and UTs
  AirMasses = [0.]
  UTs	    = [0.]

  ;Form an array big enough for all the dithers
  FullIms    = DblArr(NewArrXSize, NewArrYSize, NumRawFiles)
  ColorArr   = LonArr(NewArrXSize, NewArrYSize, 3)
  PersistArr = BytArr(KeyStr.Naxis1, KeyStr.Naxis2)
  ;Start arrays and header that will contain the names of the files used in the
  ;dither. Form Array for Median and Mean of Non-Outliers of the stack
  MedIm  = FltArr(NewArrXSize, NewArrYSize)
  MeanIm = FltArr(NewArrXSize, NewArrYSize)
  MkHdr, MeanHeader, MeanIm
  MeanHeader = Update_Header(MeanHeader)
  MkHdr, MedHeader, MedIm
  MedHeader = Update_Header(MedHeader)

  ; ****************************************************DITHER
  ; Form Dithered Image with Desired Filter
  For FileNum= 0, NumRawFiles - 1 do begin
      print, 'Dithering with File : ' + RawFiles(FileNum)

      ;Read in the frame and the saturation mask
      Fits_Read,  RawFiles(FileNum), Frame, FitsHeader, exten_no = 0
      Fits_Read,  RawFiles(FileNum), SatMask, exten_no = 1
      ;Get the Airmass and UTs
      KeyStr.AirMasses(FileNum) = Double(SxPar(FitsHeader,'AIRMASS'))
      KeyStr.UTs(FileNum) = Double(SxPar(FitsHeader,'UT'))
      print,min(frame), max(frame)
      SxAddPar, MedHeader, 'FILE'+Strtrim(FileNum,2), $
	  File_BaseName(RawFiles(FileNum))
      SxAddPar, MeanHeader, 'FILE'+Strtrim(FileNum,2), $
	  File_BaseName(RawFiles(FileNum))

      ;Zero the pixels that were previously saturated and
      If KeyStr.RejectPersist then begin
        PersistPix = where(PersistArr gt 0, NumPersist)
        If NumPersist gt 0 then begin 
	  Frame(PersistPix) = 0
	  PersistArr(PersistPix) -= 1
        EndIf
      EndIf
      ;Allow subtraction of commond mode noise before dithering
      If KeyStr.Rem60Hz then $
	Frame = Subtract_CM(Frame)
 
      ;Shift the images in the opposite direction to stack
      XOff = -(KeyStr.XOffs(FileNum)-MaxXOff)
      YOff = -(KeyStr.YOffs(FileNum)-MaxYOff)
 
      If DitherCheck eq 1 then begin
         ColorArr(XOff:XOff+KeyStr.Naxis1-1, $
		  YOff : YOff +KeyStr.Naxis2-1, $
		  *) = $
	          [[[replicate(ColNumStr.ColNums(0,FileNum),$
		   KeyStr.Naxis1,KeyStr.Naxis2)]],$ 
		   [[replicate(ColNumStr.ColNums(1,FileNum),$
		   KeyStr.Naxis1,KeyStr.Naxis2)]],$
		   [[replicate(ColNumStr.ColNums(2,FileNum),$
		   KeyStr.Naxis1,KeyStr.Naxis2)]]]
	 !p.noerase = 1
         tvimage, congrid(colorarr,700,700,3), true=3,/tv
      EndIf
      ;Put it in the right spot in the dithered large cube
      FullIms(XOff : XOff +KeyStr.Naxis1-1, $
	      YOff : YOff +KeyStr.Naxis2-1, $
	      FileNum) = Frame
 
      ;Tag the saturated pixels in the previous array and set them to zero
      ;Wait at least 2 dithers to use these pixels
      If KeyStr.RejectPersist then begin
	PersistenceMask = Return_Persistence_Mask(SatMask, RadMax= RadMax, $
		TvFlag=KeyStr.TvFlag)
	SatMask = SatMask + PersistenceMask
        PersistPix = where(SatMask gt 0,NumPersist)
        If NumPersist gt 0 then PersistArr(PersistPix) += WaitPersist
      EndIf 

  Endfor

  ; **********************************************************
  ; Run through the stacks 
  For Row = 0, NewArrYSize - 1 do begin 
    If Row mod (KeyStr.Naxis2/4) eq 0 then print, 'Got to Row: '+Strtrim(Row,2)

    For Col = 0, NewArrXSize - 1 do begin 
        PixelStack = FullIms(Col,Row,*)
        ;First off: remove Nans
        OutRange = where(finite(PixelStack) eq 0, NumOut,Complement =InRange )
	If NumOut gt 0 then PixelStack = PixelStack(InRange)

	;Next off: only deal with positive count values
        GoodPix = where(PixelStack gt 0, NumGood)

	If NumGood eq 0 then begin
           ;No good values - set to zero
           MedIm(Col,Row)  = 0
           MeanIm(Col,Row) = 0
	EndIf Else If NumGood eq 1 then begin
           ;Only one good value - go with it
           MedIm(Col,Row)  = PixelStack(GoodPix)
           MeanIm(Col,Row) = PixelStack(GoodPix)
        EndIf Else Begin
	   ;More than one value - do some cases
           PixelStackKeep = PixelStack(GoodPix)
           MeanPixStack = Mean(PixelStackKeep)
           StdDevPixStack = StdDev(PixelStackKeep)
           If NumGood eq 2 then begin
	     ;Two Good Values, take mean for both since median chooses higher
             ;Check to see if one is very high - most likely semi-hot pixel
             ;so choose lower value
             If Max(PixelStackKeep) gt 3*Min(PixelStackKeep) then begin
		MedIm(Col,Row)  = Min(PixelStackKeep)
		MeanIm(Col,Row) = Min(PixelStackKeep)
             EndIf Else Begin 
	        MedIm(Col,Row)  = Mean(PixelStack(GoodPix))
                MeanIm(Col,Row) = Mean(PixelStack(GoodPix))
	     EndElse
           EndIf Else Begin
             PixelStackKeep = PixelStack(GoodPix)
             ;Take the median - should reject outliers by definition
             MedIm(Col,Row) = Median(PixelStackKeep)

	     ;Do sigma clipping for the mean value
     	     Rej = 1
             While Rej eq 1 do begin

	       ;Look for an obvious outlier hot pixel 
               If Max(PixelStackKeep) gt 3* Min(PixelStackKeep) then $
		  PixelStackKeep = PixelStackKeep($
			where(PixelStackKeep ne max(PixelStackKeep)))
	       If N_Elements(PixelStackKeep) lt 2 then begin
		  MeanIm(Col,Row) = Mean(PixelStackKeep)
		  break
	       EndIf
	       ;Look for the mean - Reject values outside the mean
	       MedianPixStack = Median(PixelStackKeep)
	       MeanPixStack = Mean(PixelStackKeep)
               StdDevPixStack = StdDev(PixelStackKeep)
	     
	       SigThresh = 1.75
	       ;Do a few simple clips
	       If MeanPixStack gt 2*MedianPixStack then SigThresh = 1
	     
	       ClipInds=where($
	       abs(PixelStackKeep-MeanPixStack) gt SigThresh*StdDevPixStack,$
		  NumClip,complement = KeepInds)
               If NumClip gt 0 then begin
	          PixelStackKeep = PixelStackKeep(KeepInds)
                  Rej = 1
	       EndIf Else begin
                  MeanIm(Col,Row) = Mean(PixelStackKeep)
		  Rej = 0
               EndElse
	     EndWhile

	   EndElse
	EndElse
     EndFor
  EndFor

  ;**************************************************************************
  ;Update the file that holds the list of filenames and offsets with the 
  ;current ones. Don't change the star coordinates
  get_lun, OffsetsFile
  openw, OffsetsFile, FileList
  For FileNum = 0, NumRawFiles-1 do begin
     printf, OffsetsFile, Format= '(%"%s\t%I\t%I")', $
	RawFiles(FileNum), KeyStr.XOffs(FileNum), KeyStr.YOffs(FileNum)
  EndFor
  printf, OffsetsFile, Format = '(%"%s")', $
	';Bright Star Coordinates'
  For StarNum = 0, N_Elements(XStart)-1 do begin
     printf, OffsetsFile, Format= '(%"%I\t%I")', $
        XStart(StarNum), YStart(StarNum)
  EndFor
  close, OffsetsFile
  free_lun, OffsetsFile 

  ;Write the full fits files and find the offset of the most populated 
  ;area of the image
  XCenOff = -(KeyStr.XOffs(0) - MaxXOff)
  YCenOff = -(KeyStr.YOffs(0) - MaxYOff)
  FxAddPar, MedHeader, 'DITHXOFF', XCenOff , 'Offset in X of Central Dither'
  FxAddPar, MedHeader, 'DITHYOFF', YCenOff , 'Offset in Y of Central Dither'
  FxAddPar, MeanHeader, 'DITHXOFF', XCenOff , 'Offset in X of Central Dither'
  FxAddPar, MeanHeader, 'DITHYOFF', YCenOff , 'Offset in Y of Central Dither'

  ;Compute the Average AirMass using Stetson's formula with Simpson avg
  ValidAirMasses = KeyStr.AirMasses(0:NumRawFiles-1)
  ValidUTs       = KeyStr.UTs(0:NumRawFiles-1)
  ValidUTsSort   = ValidUTs(sort(ValidUTs))
  ValidAirMasses = ValidAirMasses(sort(ValidUTs))
  MiddleUT	 = ValidUTs(NumRawFiles/2)

  ;AirMass = (AM_beg+4*AM_middle+AM_last)/6
  AvgAirMass = (ValidAirMasses(0) + 4*ValidAirMasses(NumRawFiles/2) +$
	       ValidAirMasses(NumRawFiles-1))/6
  FxAddPar, MedHeader,  'AIRMASS', AvgAirMass, 'Average Airmass for mosaic'
  FxAddPar, MeanHeader, 'AIRMASS', AvgAirMass, 'Average Airmass for mosaic'
  FxAddPar, MedHeader,  'ITIME',  KeyStr.ITIME, 'Integration time (slope = 1 sec.)'
  FxAddPar, MeanHeader, 'ITIME', KeyStr.ITIME, 'Integration time (slope = 1 sec.)'
  FxAddPar, MedHeader,  'UTMID', MiddleUT, 'Universal time at middle of mosaic'
  FxAddPar, MeanHeader, 'UTMID', MiddleUT, 'Universal time at middle of mosaic'
 
  Writefits, KeyStr.ObjectDir+KeyStr.ObjectName+'_'+$
	Strtrim(KeyStr.ThisFilter,2)+'_'+$
	'Slope'+FlatFieldStr+KeyStr.Rem60HzStr+$
	KeyStr.RejectPersistStr+$
	'DitheredMedian.fits', MedIm, MedHeader
  Writefits, KeyStr.ObjectDir+KeyStr.ObjectName+'_'+$
	Strtrim(KeyStr.ThisFilter,2)+'_'+$
	'Slope'+FlatFieldStr+KeyStr.Rem60HzStr+$
	KeyStr.RejectPersistStr+$
	'DitheredMean.fits', MeanIm, MeanHeader

  ;Form the image of the center file that has Naxis1 columns and Naxis2 rows
  XCenOff = -(KeyStr.XOffs(0) - MaxXOff) 
  YCenOff = -(KeyStr.YOffs(0) - MaxYOff)
  MeanCenIm = MeanIm(XCenOff : XCenOff + KeyStr.Naxis1-1, $
		     YCenOff : YCenOff + KeyStr.Naxis2-1 )
  MedCenIm  =  MedIm(XCenOff : XCenOff + KeyStr.Naxis1-1, $
                     YCenOff : YCenOff + KeyStr.Naxis2-1 )

  ;Modify the headers
  FxAddPar, MedHeader,  'NAXIS1', KeyStr.Naxis1
  FxAddPar, MedHeader,  'NAXIS2', KeyStr.Naxis2
  FxAddPar, MeanHeader, 'NAXIS1', KeyStr.Naxis1
  FxAddPar, MeanHeader, 'NAXIS2', KeyStr.Naxis2
  Writefits, KeyStr.ObjectDir+KeyStr.ObjectName+'_'+$
	Strtrim(KeyStr.ThisFilter,2)+'_'+$
        'Slope'+FlatFieldStr+KeyStr.Rem60HzStr+$
	KeyStr.RejectPersistStr+$
	'CenteredDitheredMedian.fits', MedCenIm, MedHeader
  Writefits, KeyStr.ObjectDir+KeyStr.ObjectName+'_'+$
	Strtrim(KeyStr.ThisFilter,2)+'_'+$
        'Slope'+FlatFieldStr+KeyStr.Rem60HzStr+$
	KeyStr.RejectPersistStr+$
	'CenteredDitheredMean.fits', MeanCenIm, MeanHeader

  stop


end