;pro ShowStar
;	a function to plot the star and some quantities if requested
pro ShowStar, FileName, Centroid, StarIm, XStart, YStart, FileNum, IterNum

Common KeyParams, KeyStr
Common StarKeys, ImageSize, Masks

	 If FileNum lt KeyStr.NumDithers then begin
	   FileParts = StrSplit(FileName, '/', /extract, count = NumParts)
           FileKey   = FileParts(NumParts-1)
	 EndIf else begin
	   FileKey = 'Median of Dithered Star'
	 EndElse

         StarImTmp = StarIm
         BadPix = where(StarImTmp lt 0, NumBadPix)
         If NumBadPix gt 0 then StarImTmp(BadPix) = 0
         plot, [0],[0], /nodata, background = fsc_color('white'), $
	       xstyle=4,ystyle=4
         !p.noerase = 1
         xyouts, 0.05,0.9, FileKey,$
         /normal, color=fsc_color('black')
         xyouts, 0.7,0.80, 'Centroid : ' + Strtrim(Centroid(0),2)+$
         ', '+Strtrim(Centroid(1),2),$
         /normal, color=fsc_color('black')
         xyouts, 0.7,0.75, 'Second Moments : ' + Strtrim(sqrt(Centroid(2)),2)+$
	 ', '+Strtrim(sqrt(Centroid(3)),2),$
         /normal, color=fsc_color('black')
         xyouts, 0.7,0.70, 'Flux : ' + Strtrim(Total(StarImTmp),2), $
         /normal, color=fsc_color('black')
         xyouts, 0.7,0.65, 'Min : ' +Strtrim(min(StarIm),2), $
         /normal, color=fsc_color('black')
         xyouts, 0.7,0.60, 'Max : ' + Strtrim(max(StarIm),2), $
         /normal, color=fsc_color('black')
         MinIm = Mean(StarIm)-2*StdDev(StarIm)
         MaxIm = Mean(StarIm)+2*StdDev(StarIm)
         StarImDisplay = bytscl(StarIm, min = MinIm, max = MaxIm)
	 StarImDisplay = congrid(StarImDisplay,512,512)
         StarImDisplay = [[[StarImDisplay]],[[StarImDisplay]],[[StarImDisplay]]] 
         tvimage, StarImDisplay, $
             position = [0.025,0.175,0.675,0.825], $
	     /normal,true=3

	 ;Map bad pixels for raw slopes
	 If FileNum lt KeyStr.NumDithers then begin
           ;Identify the bad pixels in the image
           ZeroedPixInd = where(StarIm eq 0, NumZero)
           xyouts, 0.7,0.55, 'Number of Zeroed Pixels : ' + $
           Strtrim(NumZero,2), $
           /normal, color=fsc_color('black')
           If NumZero gt 0 then begin
             BadPixelTypes = StrArr(NumZero)
             For i = 0, NumZero - 1 do begin
               ZeroedPixTInd = $
                long(XStart) + KeyStr.RAs(FileNum) + $
                ZeroedPixInd(i) mod ImageSize + $
                (long(YStart) + KeyStr.DECs(FileNum) + $
                ZeroedPixInd(i)/ImageSize)*$
                KeyStr.Naxis1
                BadPixelTypes(i) = Return_Pixel_Type(ZeroedPixTInd, Masks)
               EndFor
           
              HotInds   = where(BadPixelTypes eq 'Hot'  , NumHot)
              LeakyInds = where(BadPixelTypes eq 'Leaky', NumOpen)
              DeadInds  = where(BadPixelTypes eq 'Dead' , NumDead)
	    EndIf else begin
              NumHot  = 0
              NumOpen = 0
	      NumDead = 0
	    EndElse 
            xyouts,0.7,0.50, 'Number of Hot Pixels :' + Strtrim(NumHot,2),$
            /normal, color=fsc_color('black')
            xyouts,0.7,0.45, 'Number of Open Pixels:' + Strtrim(NumOpen,2),$
            /normal, color=fsc_color('black')
            xyouts,0.7,0.40, 'Number of Dead Pixels:' + Strtrim(NumDead,2),$
            /normal, color=fsc_color('black')
	  EndIf Else begin
            xyouts,0.7,0.55, 'Number of Iterations to find centroid : ' + $
		Strtrim(IterNum,2), $
		/normal, color=fsc_color('black')
	  EndElse
          !p.noerase = 0
          If KeyStr.TvFlag eq 1 then stop
	  erase
End

;pro return_dither_offsets
;
;PURPOSE:
;	To find the best offset coordinates for the 9 dithers taken in 
;	a dither sequence for one field in one filter.
;
;	The function operates by finding a bright star, calculating the 
;	centroid, and returning the coordinates of that centroid. 
;	The first guess at the offsets should already be instantiated 
;	in KeyStr.RAs and KeyStr.DECs
;
;KEYWORDS:

pro Return_Dither_Offsets, XStart, YStart, $
		   RootDir=RootDir, Date=Date, $
		   ReducedDir = ReducedDir, $
		   SkySubtract= SkySubtract, ObjectName=ObjectName, $
		   NReads = NReads, ThisFilter = ThisFilter, $
		   FlatFielded = FlatFielded, SkipFileNum = SkipFileNum, $
		   StarRegionSize= StarRegionSize, TvFlag = TvFlag, $
		   WriteFits = WriteFits, SkipStar = SkipStar

  ;Make some common blocks
  Common KeyParams, KeyStr
  Common StarKeys, ImageSize, Masks

  If N_Elements(StarRegionSize) eq 0 then StarRegionSize = 30
  If N_Elements(WriteFits) eq 0 then WriteFits = 0
  If N_Elements(SkipStar) eq 0 then SkipStar = 0
  ImageSize = StarRegionSize

  ;Form a filename for the output postscript if requested
  If TvFlag eq 2 then $
        device, FileName = PlotDir+ObjectName+'DitheredStar_'+$
	  Strtrim(KeyStr.ThisFilter,2)+'_'+$
	  Strtrim(XStart,2)+'_'+Strtrim(YStart,2)+'.ps'

  ;Get the pixel masks in case an identification is required
  Masks = Return_Pixel_Masks()

  ;Form arrays to hold the star in each of the 9 images
  StarArr = Fltarr(StarRegionSize, StarRegionSize, KeyStr.NumDithers)

  ; *****************************************************************
  ; Adjust Offsets until all of the star centroids line up
  ; CentroidsFound = 1 after all offsets yield centered images
  CentroidsFound = 0 & TwiceThrough = 0
  If TvFlag eq 2 then  TwiceThrough = 1
  IterNum = 0
  NumStars = N_Elements(XStart)

  ; Loop through a few stars to get the best offsets
  For StarNum = 0, NumStars-1 do begin

   ;Form arrays to hold the star in each of the 9 images
   StarArr = Fltarr(StarRegionSize, StarRegionSize, KeyStr.NumDithers)

   For PlotSeq = 0, TwiceThrough do begin
    CentroidsFound = 0
    While CentroidsFound ne 1 do begin
     ;Each time check to see if all of the centroids are at the center
     CentroidsFound = 1
     If PlotSeq eq 0 then IterNum += 1

     ;Determine the number of files that are being used
     NumRawFiles = KeyStr.NumDithers

     For FileNum= 0, NumRawFiles - 1 do begin
       print, 'Dithering with File : ' + KeyStr.DitherFiles(FileNum)
       ;Read in the sub region containing the star
       Fits_Read_DataCube, KeyStr.DitherFiles(FileNum), StarIm, StarHeader, $
 	 xstart = XStart(StarNum)+KeyStr.RAs(FileNum, StarNum), $
	 xstop  = XStart(StarNum)+KeyStr.RAs(FileNum, StarNum)+StarRegionSize-1, $
	 ystart = YStart(StarNum)+KeyStr.Decs(FileNum, StarNum), $
	 ystop  = YStart(StarNum)+KeyStr.Decs(FileNum, StarNum)+StarRegionSize-1

       ;Put the Star in the Array
       StarArr(*,*,FileNum) = StarIm
       
       ;Do a simple centroid detection without sigma clipping or 
       ;background subtraction
       Centroid = Return_Centroid(StarIm,/RemoveBad)

       ;Test to see if the centroid is more than a pixel off
       CentroidEval = where(abs(Centroid(0:1)) gt 1, NumOff)
       If NumOff gt 0 then CentroidsFound = 0
       print, 'Star number ' +Strtrim(StarNum, 2)
       print, NumOff, CentroidsFound 
       ;Allow user to view on the screen to validate correct centroid
       If KeyStr.TvFlag eq 1 and (StarNum ge SkipStar) $
	 or (KeyStr.TvFlag eq 2 and PlotSeq eq 1) then begin
	 print,'plotting'
         ShowStar, KeyStr.DitherFiles(FileNum), Centroid, StarIm, $
		XStart, YStart, FileNum
       EndIf
       
       ;Adjust the offsets for another round
       KeyStr.RAs(FileNum, StarNum)=KeyStr.RAs(FileNum,StarNum)+$
		Round(Centroid(0))
       KeyStr.DECs(FileNum, StarNum)=KeyStr.DECs(FileNum,StarNum)-$
		Round(Centroid(1))
       print, 'XOffsets for Star '+Strtrim(StarNum,2)+ ': ',keystr.Ras(FileNum)
       print, 'YOffsets for Star '+Strtrim(StarNum,2)+ ': ',keystr.Decs(FileNum)

      Endfor
     EndWhile
    EndFor
    print, 'Offsets in X:' , KeyStr.RAs(*,StarNum)
    print, 'Offsets in Y:' , KeyStr.DECs(*,StarNum)
    ; **********************************************************
    ; Take the Median of the stack
    MedIm  = median(StarArr(*,*,0:NumRawFiles-1), dimension = 3)
    Centroid =Return_Centroid(MedIm)
    If KeyStr.TvFlag eq 1 or (TvFlag eq 2) then $
         ShowStar, '', Centroid, MedIm, 0, 0, 9, IterNum

    ; Adjust the offsets so that Dither 0's at the center for this star
    KeyStr.RAs(*,StarNum)  = KeyStr.RAs(*,StarNum)-KeyStr.RAs(0,StarNum)
    KeyStr.DECs(*,StarNum) = KeyStr.DECs(*,StarNum)-KeyStr.DECs(0,StarNum)

    ; Give these to the next star so it has a better start
    IF StarNum lt NumStars -1 then begin
       KeyStr.RAs(*,StarNum+1) = KeyStr.RAs(*,StarNum)+StarRegionSize/2 
       KeyStr.DEcs(*,StarNum+1) = KeyStr.DECs(*,StarNum)+StarRegionSize/2
    EndIf

   EndFor
   
   ;Average the offsets from the various stars that were used
   RAOffs = KeyStr.RAs(*,0:NumStars-1)
   DECOffs = KeyStr.DECs(*,0:NumStars-1)
   If NumStars eq 1 then begin
     KeyStr.XOffs = RAOffs
     KeyStr.YOffs = DECOffs
   EndIf else begin 
     KeyStr.XOffs = Round(Avg(RAOffs,1))
     KeyStr.YOffs = Round(Avg(DECOffs,1))
   EndElse
   ; **********************************************************
   ; Take the Median of the stack
   MedIm  = median(StarArr(*,*,0:NumRawFiles-1), dimension = 3)
   Centroid =Return_Centroid(MedIm)
   If TvFlag eq 1 or (TvFlag eq 2) then $
         ShowStar, '', Centroid, MedIm, 0, 0, KeyStr.NumDithers, IterNum

   If WriteFits then begin
    MkHdr, MedHeader, MedIm
    Writefits, ReducedDir+ObjectName+'_'+Strtrim(KeyStr.ThisFilter,2)+$
	'_'+'Slope'+KeyStr.FlatFieldStr+'Dithered.fits', MedIm, MedHeader
   EndIF
   If TvFlag eq 2 then device, /close


end
