;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
;pro Return_Persistence_Mask_xRG
;
;PURPOSE:
;
;   To take in a mask of pixels that were saturated in an image and
;   extend the mask radially outward to enclose all of the pixels that
;   have either negative or positive persistence,
;
;INPUTS:
;	SatMask- The mask of saturated pixels.  In most cases it should 
;		       
;KEYWORDS:
;
function Return_Persistence_Mask_xRG, SatMask, TvFlag = TvFlag, RadMax = RadMax, $
	 Verbose=Verbose

If N_Elements(TvFlag)  eq 0 then TvFlag = 0 
If N_Elements(RadMax)  eq 0 then RadMax = 20
If N_Elements(Verbose) eq 0 then Verbose = 0
SatRegSize = 2*RadMax+1

Common KeyParams, KeyStr

;The mask to be returned will have 1's in circles around bright stars
If KeyStr.ElecStr eq 'ASIC' then begin
  PersistenceMask = BytArr(KeyStr.Naxis1-2*KeyStr.Ref/2, KeyStr.Naxis2-2*KeyStr.Ref/2)
EndIf Else Begin
  PersistenceMask = BytArr(KeyStr.Naxis1-2*KeyStr.Ref/2-1, KeyStr.Naxis2-2*KeyStr.Ref/2)
EndElse

;Find the regions that were occupied by bright stars as opposed to 
;very high E cosmic rays or messed up pixels
HMin     = 2.		; minimum counts necessary to see a peak for find.pro
FWHM     = 6.		; FWHM cutoff for find.pro
RoundLim = [-2.,2.]	; roundness paramtetrs for find.pro
SharpLim = [0.1,1.2]	; sharpness parameters for find.pro
PhotAper = [2.]		; the aperture radius used by aper.pro
RdNoise  = [3.]		; the readnoise used by psf.pro
PhPerADU = [1.]		; the number of photons per adu
PSFRad   = 20.		; the radius of the circular area in which psf will be def
FitRad   = 12.
TReads   = KeyStr.TotalReads

Sky,  SatMask, SkyMode, SkySig
Find, SatMask, XStar, YStar, Flux, Sharp, Roundness, HMin, $
	FWHM, RoundLim, SharpLim,/silent 
If N_Elements(XStar) eq 0 then return, PersistenceMask

Aper, SatMask, XStar, YStar, Mag, ErrAp, Sky, SkyErr, 1., PhotAper, -1, $
        [-1., TReads+1], setskyval = SkyMode, /silent

NumStars = N_Elements(XStar)
XSigs = Fltarr(NumStars) & YSigs = Fltarr(NumStars)
print ,'NumStars:' +Strtrim(NumStars,2)
For StarNum = 0, NumStars-1 do begin
  
    XStart = (Round(XStar(StarNum)) - RadMax) > 0 
    XStop  = (Round(XStar(StarNum)) + RadMax) < ((size(SatMask))[1]-1)
    YStart = (Round(YStar(StarNum)) - RadMax) > 0 
    YStop  = (Round(YStar(StarNum)) + RadMax) < ((size(SatMask))[2]-1)

    Centroid = Find_Centroid_Sat_xRG(SatMask, XStart, XStop, YStart, YStop, IsStar, $
			             SatRad, TvFlag = TvFlag, Verbose = KeyStr.Verbose, $
                                     ElecStr=KeyStr.ElecStr, RefPix = KeyStr.Ref/2)
    XSigs(StarNum) = Centroid(2) & YSigs(StarNum) = Centroid(3)

    ;Check if we'll stay in bounds
    If KeyStr.ElecStr eq 'ASIC' then begin
      If Centroid(0) + SatRegSize/2 gt (KeyStr.Naxis1-2*KeyStr.Ref/2) or $
         Centroid(0) - SatRegSize/2 lt 0 or $
         Centroid(1) + SatRegSize/2 gt (KeyStr.Naxis2-2*KeyStr.Ref/2) or $
         Centroid(1) - SatRegSize/2 lt 0 then goto, SkipStar
    Endif Else if KeyStr.ElecStr eq 'LEACH' then begin
      If Centroid(0) + SatRegSize/2 gt (KeyStr.Naxis1-2*KeyStr.Ref/2-1) or $
         Centroid(0) - SatRegSize/2 lt 0 or $
         Centroid(1) + SatRegSize/2 gt (KeyStr.Naxis2-2*KeyStr.Ref/2)  or $
         Centroid(1) - SatRegSize/2 lt 0 then goto, SkipStar

    Endif
    RelXCoords = Indgen(SatRegSize,SatRegSize) mod SatRegSize
    RelXCoords = RelXCoords - Mean(RelXCoords)
    RelYCoords = reverse(Indgen(SatRegSize,SatRegSize),2) / SatRegSize
    RelYCoords = RelYCoords - Mean(RelYCoords)
    RsqCoords  = RelXCoords^2+RelYCoords^2
    If IsStar then begin
       RelPixInRad = where(sqrt(RsqCoords) lt RadMax)
       AbsPixInds  = RelXCoords(RelPixInRad) + Centroid(0) + $
			(RelYCoords(RelPixInRad) + Centroid(1))* KeyStr.Naxis1
       PersistenceMask(AbsPixInds) = 1

    EndIf	
    SkipStar:  ;The star was too far off

EndFor

return, PersistenceMask

End