;function find_centroid_sat
;
;PURPOSE:
;	To find a centroid of a small rectangular region that is known to 
;	contain a star. This program is to be used in conjunction with the 
;       script examine_persistent_pixels.pro.  find_centroid.pro should be used
;       to simply return a centroid. 
;
;INPUTS:
;	Image - The 2d image containing the star
;
;OUTPUTS: A 4 element array
;	  1 - X Centroid Coord
;	  2 - Y Centroid Coord
;	  3 - Second Moment in X
;	  4 - Second Moment in Y
;
;	  SatRad = Radius at which the star becomes unsaturated
;KEYWORDS:
;	TVFLAG
;	  0 - Don't plot or Tv
;	  1 - Tv the star and print out the centroids as iterated
;
;	SATURATIONLIM:
;	  Set this keyword to the saturation value for the star. 
;	  When the mean number of counts per pixel falls below this
;	  for a given aperture radius, this radius will be a good estimate
;	  for the radius that contains the saturated pixels.
;	IsStar:
;	  0 if the object is not a star.  1 if it is a star
;
function find_centroid_sat, Image, XStart, XStop, YStart, YStop, IsStar, $
	 SatRad, TvFlag = TvFlag, $
	 SaturationLim = SaturationLim, Verbose = Verbose

If N_Elements(TvFlag) 	     eq 0 then TvFlag = 0
If N_Elements(SaturationLim) eq 0 then SaturationLim = 1
If N_Elements(Verbose)       eq 0 then Verbose = 0
SatRad = 0
IsStar = 1
Naxis1 = (Size(Image))[1] & Naxis2 = (Size(Image))[2]
X = 0D & Y = 0D
XStartL = XStart & XStopL = XStop
YStartL = YStart & YStopL = YStop
LoopNum = 0
Repeat Begin 
  LoopNum+=1
  SubImage = Image(XStartL:XStopL, YStartL: YStopL)
 
  ;Make the arrays with negative and positive indices for easy weighting 
  SubImSize  = Size(SubImage)
  XCoords    = Indgen(SubImSize(1),SubImSize(2)) mod SubImSize(1) 
  XCoords    = XCoords - Mean(XCoords)
  YCoords    = reverse(Indgen(SubImSize(1),SubImSize(2)),2) / SubImSize(1) 
  YCoords    = YCoords - Mean(YCoords)

  ;Calculate Flux, First and Second Moments
  Flux = Total(SubImage)
  X    = Total(XCoords*SubImage)/Flux
  Y    = Total(YCoords*SubImage)/Flux
  X2   = Total(XCoords^2*SubImage)/Flux
  Y2   = Total(YCoords^2*SubImage)/Flux

  ;Shift the image over and try to recentroid
  XStartL = XStartL + round(X) > 0
  XStopL  = XStopL  + round(X) < Naxis1
  YStartL = YStartL - round(Y) > 0
  YStopL  = YStopL  - round(Y) < Naxis2

  ;Set TvFlag = 1 to check if centroiding is being done correctly
  If TvFlag eq 1 then begin
     If TvFlag eq 1 then window,xsize=768,ysize=768
     tvimage, congrid(bytscl(SubImage),512,512), $
	background= fsc_color('white'),/erase
     print, 'X Centroid: '+Strtrim(X,2)
     print, 'Y Centroid: '+Strtrim(Y,2)
     print, 'X Sec Mom : '+Strtrim(sqrt(X2),2)
     print, 'Y Sec Mom : '+Strtrim(sqrt(Y2),2) 
  EndIf
EndRep Until Abs(X) lt 1 and Abs(Y) lt 1 or (LoopNum gt 40)

;Now find out where the outer edge of the star/saturated mask is
;by checking flux in different apertures
RCoords = XCoords^2 +YCoords^2
MaxRad  = fix(Max(Sqrt(RCoords)))
FluxIn  = Fltarr(MaxRad+1) & FluxOut = Fltarr(MaxRad+1)  
;Keep track of mean count number In, Out, and in 1 pixel Annulus
MeanCountOut = Fltarr(MaxRad+1)
MeanCountAnn = Fltarr(MaxRad+1)
MeanCountIn  = Fltarr(MaxRad+1)

For Rad = 0, fix(MaxRad) do begin
    InRadPix  = where(sqrt(RCoords) le Rad, complement = OutRadPix)
    AnnRadPix = where(sqrt(RCoords) ge Rad $
                  and sqrt(RCoords) lt Rad+1)
    FluxIn(Rad)       = total(SubImage(InRadPix)) 
    FluxOut(Rad)      = total(SubImage(OutRadPix))
    MeanCountIn(Rad)  = Mean(SubImage(InRadPix)) 
    MeanCountOut(Rad) = Mean(SubImage(OutRadPix))
    MeanCountAnn(Rad) = Mean(SubImage(AnnRadPix))    
    If SatRad eq 0 and MeanCountAnn(Rad) lt SaturationLim then $
       SatRad = Rad
   
EndFor

;Determine whether this is really a star
If SatRad eq 0 or Flux lt 15*Max(SubImage) or $
   X2 lt 3 or Y2 lt 3 then IsStar = 0

;The centroid might be different if on the edge
XCen = XStartL + Round(Float(XStopL-XStartL)/2)
YCen = YStartL + Round(Float(YStopL-YStartL)/2)

If Verbose then begin 
  If IsStar then print, 'This is a star'
  If N_Elements(SatRad) gt 0 then print, 'Saturation Radius :'+Strtrim(SatRad,2)
  If IsStar then print, 'XCen : ', XCen, 'YCen : ', YCen
EndIf

If TvFlag then stop

return, [ XCen, YCen, X2, Y2]
end