##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009 ##Return_Persistence_Quantities ## ##PURPOSE: ## ## To take in a box of pixels and return quantities related to the ## persistence. ## ## ##NOTE: The format of pyfits data is [Frame, Row, Col] !!!!!!!!!!!!!!!!!!! import pyfits import pylab from sg_filter import * import pdb import matplotlib.pylab as mplot import time import scipy.optimize as optimize from FitGaussian import * from ReturnRadialAvg import * from SpatialJacobian import * ##Return_Sat_Quantities ## ##PURPOSE: ## To take in the saturated box of pixels (a slopefitted image) ## and return quantities that define the degree of saturation ## ##INPUTS: ## SubIm: 2-d array ## The YxX box of pixels that ## HxRG: class instance ## An object of the HxRG class that contains all the necessary members ## LeakySubMask: 2-d array ## A mask of the leaky pixels ## DarkSubMask: 2-d array ## A mask of the hot pixels ## DarkIm: 2-d array ## The sub-image in the dark exposure after the saturation ## ##KEYWORDS: ## XCen: int ## The x-coordinate ## YCen: int ## The y-coordinate ## Mag: float ## The magnitude returned by DAOFind ## Sharp: float ## The roundness returned by DAOFind ## TvFlag: int ## 0 -Don't plot ## 1 -Plot and output parameters ## TvText: int ## 0 - Don't output quantities to the figure ## 1 - Output quantities to the plot ## ##OUTPUTS: ## MaxSlopeSec: float ## The maximum slope in units of Int/Second ## MaxSlopeRd: float ## The maximum slope in units of Int/Reads ## MaxCount: float ## The maximum slope times the number of reads ## WellPer: float ## The percentage of the well that was filled assuming that the ## charge accumulates linearly and that HxRG.WellDepth is the maximum ## well depth ## NumSat: int ## The number of pixels saturated ## NumLeaky: int ## The number of leaky pixels in the box ## def Return_Sat_Quantities(SubIm, HxRG, XCen=0, YCen=0, \ Mag=0, Sharp=0, LeakySubMask=[0],\ HotSubMask=[0], TvFlag=0, \ Verbosity=0, DarkIm=[0], SRawIm=[0],\ TvText=0): #Provide dumby arrays if Masks are not provided if len(LeakySubMask) == 1: LeakySubMask = zeros([SubIm.shape[0], SubIm.shape[1]], dtype=byte) if len(HotSubMask) == 1: HotSubMask = zeros([SubIm.shape[0], SubIm.shape[1]], dtype=byte) #Maximum slope in both units of Int/SubIm and Int/Sec MaxSlopeSec = SubIm.max() MaxSlopeSecStr = "%.2f" % MaxSlopeSec MaxSlopeRd = SubIm.max()*HxRG.FrameTime MaxSlopeRdStr = "%.2f" % MaxSlopeRd MaxCount = MaxSlopeSec*HxRG.ITime MaxCountStr = "%.2f" % MaxCount MaxDark = 0 #Fraction of well filled WellPer = MaxCount/HxRG.FullWell WellPerStr = "%.2f" % WellPer #Number of pixels saturated TotalCounts = SubIm*HxRG.ITime SatCri = TotalCounts > HxRG.FullWell NumSat = len((where(SatCri))[0]) #Radius of saturated pixels ImSizeY = SubIm.shape[0] ; MidY = ImSizeY/2 ImSizeX = SubIm.shape[1] ; MidX = ImSizeX/2 if ImSizeY % 2 == 1 and ImSizeX % 2 ==1: ImCoo = mgrid[-MidX:MidX+1, -MidY:MidY+1] else: ImCoo = mgrid[-MidX:MidX, -MidY:MidY]+0.5 Radii = sqrt(ImCoo[0]**2+ImCoo[1]**2) RadMax = Radii.max() if NumSat > 0 : SatRad = Radii[where(SatCri)] SatRadMax = max(SatRad) else: SatRadMax = 0 #Do a radial profile of the slopefit SIntArr = ReturnRadialAvg(SubIm,Mode=0, Verbosity=Verbosity) SRadii = SIntArr[:,0] SIntArr = SIntArr[:,1] SRad1d = SRadii.copy() SRad1d.shape = SRad1d.size SIntArr1d = SIntArr.copy() SIntArr1d.shape = SIntArr1d.size #Estimate the sky values SubIm = SubIm SkyM = mean(SubIm[(where(Radii > .85*RadMax))[0]]) SkyV = std(SubIm[(where(Radii > .85*RadMax))[0]]) SCrit = (where(SubIm < SkyM+2*SkyV)) SkyR = Radii[SCrit].min() HM = .5*SubIm.max() HCrit = (where(SubIm > HM)) FWHM0 = 2*Radii[HCrit].max() Height, MuX, MuY, FWHMX, FWHMY, FWHM = moments(SubIm) #Number of Leaky Pixels in the box NumLeaky = len((where(LeakySubMask != 0))[0]) NumHot = len((where(HotSubMask != 0))[0]) if TvFlag == 1: #Show the star and write some values to the figure mplot.clf() if len(DarkIm) == 1: mplot.imshow(SubIm) MaxDark = 0 else: #Do a radial profile of the slopefit DIntArr = ReturnRadialAvg(DarkIm,Mode=0, Verbosity=Verbosity) DRadii = DIntArr[:,0] DIntArr = DIntArr[:,1] DRad1d = DRadii.copy() DRad1d.shape = DRad1d.size DIntArr1d = DIntArr.copy() DIntArr1d.shape = DIntArr1d.size mplot.subplot(2,2,1) mplot.imshow(SubIm) mplot.subplot(2,2,2) mplot.imshow(DarkIm) mplot.subplot(2,2,3) mplot.plot(SRad1d, SIntArr1d/max(SIntArr1d), 'bo') mplot.plot(DRad1d, DIntArr1d/max(DIntArr1d), 'ro') MaxDark = DIntArr1d.max() if MaxDark > MaxCount: print 'Likely Bad Dark at X: ' +\ str(XCen)+' , '+str(YCen) if TvText == 1: mplot.figtext(0.03,0.8, 'X : '+ str(XCen)) mplot.figtext(0.03,0.75, 'Y : '+ str(YCen)) mplot.figtext(0.03,0.70, 'Mag : '+ str(Mag)) mplot.figtext(0.03,0.65, 'Sharp : '+ str(Sharp)) mplot.figtext(0.03,0.60, 'Max Sls : '+ str(MaxSlopeSecStr)) mplot.figtext(0.03,0.55, 'Max Slr : '+ str(MaxSlopeRdStr)) mplot.figtext(0.03,0.50, 'Max Cnt : '+ str(MaxCountStr)) mplot.figtext(0.03,0.45, 'Well Per: '+ str(WellPerStr)) mplot.figtext(0.03,0.40, 'Num Sat : '+ str(NumSat)) mplot.figtext(0.03,0.35, 'Sat Rad : '+ str(SatRadMax)) mplot.figtext(0.03,0.30, 'Num Leak: '+ str(NumLeaky)) mplot.figtext(0.03,0.25, 'Max Dark: '+ str(MaxDark)) if Verbosity == 1: print 'X : '+ str(XCen) print 'Y : '+ str(YCen) print 'Mag : '+ str(Mag) print 'Sharp : '+ str(Sharp) print 'Max Sls : '+ str(MaxSlopeSecStr) print 'Max Slr : '+ str(MaxSlopeRdStr) print 'Max Cnt : '+ str(MaxCountStr) print 'Well Per: '+ str(WellPerStr) print 'Num Sat : '+ str(NumSat) print 'Sat Rad : '+ str(SatRadMax) print 'Num Leak: '+ str(NumLeaky) print 'Num Hot : '+ str(NumHot) print 'Max Dark: '+ str(MaxDark) print 'Sky : '+ str(SkyM) print 'Sky Rad : '+ str(SkyR) print 'FWHM: : '+ str(FWHM) print '\n\n' #Pause for a moment if the output is being displayed if TvFlag > 0: if MaxDark > MaxCount: mplot.figtext(.55,0.2,\ 'Likely Bad Dark at X: '+str(XCen)+' , '+str(YCen)) time.sleep(2) del SubIm, DarkIm return MaxSlopeSec, MaxSlopeRd, MaxCount, WellPer, NumSat, \ NumLeaky, NumHot, SatRadMax, FWHM, SkyR, SkyM ##Return_Persistence_Quantities ## ##PURPOSE: ##To return the negative persistence quantities from a sub-image ## ##INPUTS: ## SubIm: 2-d array ## The YxX box of pixels that ## HxRG: class instance ## An object of the HxRG class that contains all the necessary members ## LeakySubMask: 2-d array ## A mask of pixels ## AttFit: int ## 0 - Don't attempt to fit the radial profile ## 1 - Attempt to fit the radial profile ## NegPersistence: int ## 0 - The subimage is from a saturation exposure ## 1 - The subimage is from an image following the saturation ## This Keyword must be set to find the zero radii ## Verbosity: ## 0 - Don't output text ## 1 - Output some text ## ReadNum: int ## The read from which the sub-image is taken ## ExpFun1: int ## Choice of the function to fit the radial profile for a given read ## 1 - A*Cos(B*r^C)*Exp(-r^2/D) ## 2 - A*Cos(B*r^C)*Exp(-r^2/D)+E ## 3 - A*Cos(B*r&C)*r**D ## NOTE: For all database entries ExpFunc=1 was used. ## ExpFun2: int ## Function to fit rise in signal due to persistence (IMax) ## 1 - A*(1-Exp(-t/B))+C ## 2 - ## FitFunc: int ## 0 - Fit with a cosine modulated by a gaussian ## 1 - Fit with a cosine modulated by a polynomial ## FitThreshold: int ## The Intensity value that determines whether or not a fit should be applied. ## If all values fall below FitThreshold, a fit will not be attempted. ## The fitting will take an extremely long time if noise is fit. ## FromRadProf: int ## 0 - the input is a 2-d cutout from an image ## 1 - the image is a radial profile (radii must be included in Radii keyword) ## Radii: float array ## The radii used when FromRadProf = 1 ## ##OUTPUTS ## MinInt: float ## Minimum in intensity ## MaxInt: float ## Maximum in intensity ## MinRad: float ## Radius at MinInt ## MaxRad: float ## Radius at MaxInt ## FirstZRad: float ## Radius at first zero (R_core) ## SecondZRad: float ## Radius at second zero (R_halo) ## SumCore: float ## Sum of counts from r=0 - r=FirstZRad ## SumHalo: float ## Sum of counts from r=FirstZRad - r=SecondZRad ## SumOut: float ## Sum of counts from r=SecondZRad- r=R_Max ## SumInt: flota ## Sum of all counts in box def Return_Persistence_Quantities(SubIm, HxRG, AttFit, NegPersistence, \ TvFlag=0, Verbosity=0, ReadNum=0, ExpFunc1=2,\ ExpFunc2 = 2, FitThreshold = 25, FromRadProf=0, Radii=0): ##FITTING FUNCTIONS .................................................... ##Radial profile fitting function if ExpFunc1 == 1: FitFunc = lambda p, x: p[0]*cos(p[1]*x**p[3])*exp(-x**2/p[2]) elif ExpFunc1 == 2: FitFunc = lambda p, x: p[0]*cos(p[1]*x**p[2])*exp(-x**2/p[3])+p[4] elif ExpFunc1 == 3: FitFunc = lambda p, x: p[0]*cos((1./(2*pi))*x**p[1])*exp(-x**2/p[2])+p[3] else: FitFunc = lambda p, x: p[0]*cos(p[1]*x**p[3])*x**(p[2]) ##The error function for the radial profile ErrFunc1 = lambda p, x, y: y-FitFunc(p,x) ##Rise of intensity at center function if ExpFunc2 == 1: RiseFunc = lambda p, x: p[0]*(1-exp(-x/p[1]))+p[2] elif ExpFunc2 == 2: RiseFunc = lambda p, x: p[0]*(1-exp(-x/p[1]))+p[2]*(1-exp(-x/p[3]))+p[4] ##The error function for the rise function ErrFunc2 = lambda p, x, y: y-RiseFunc(p,x) ##..................................................................... if FromRadProf == 0: #Add the radial plot and get the averaged min/max'es IntArr3 = ReturnRadialAvg(SubIm,Mode=0, Verbosity=Verbosity) Radii = IntArr3[:,0] IntArr = IntArr3[:,1] #Get the full intensity array FullRadii, FullIntArr = ReturnRadialAvg(SubIm, Mode=1, Verbosity=Verbosity) #Get the minimum and maximum radii Rad1d = Radii.copy() Rad1d.shape = Rad1d.size IntArr1d = IntArr.copy() IntArr1d.shape = IntArr1d.size FullRad1d = FullRadii.copy() FullRad1d.shape = FullRad1d.size FullInt1d = FullIntArr.copy() FullInt1d.shape = FullInt1d.size if AttFit == 1: FitRad1d = arange(0,FullRad1d.max(),0.1) else: Rad1d = Radii.copy() IntArr1d = SubIm.copy() FullRad1d = Radii.copy() FullInt1d = SubIm.copy() #Try to make sure that maximum value is near center of radius MaxInt = IntArr1d.max() MinInt = IntArr1d.min() MinIntInd = (where(IntArr1d == MinInt))[0] MaxIntInd = (where(IntArr1d == MaxInt))[0] if size(MinIntInd) != 0: MinRad = Rad1d[MinIntInd[0]] else: MinRad = 0 if size(MaxIntInd) != 0: MaxRad = Rad1d[MaxIntInd[0]] else: MaxRad = 0 SumInt = IntArr1d.sum() if NegPersistence == 1: ##Attempt a fit if AttFit = 1 if AttFit == 1 and MaxInt > FitThreshold: #Initial Guess if ExpFunc1 == 1: p0=[float(MaxInt), 0.3, 40, .5] elif ExpFunc1 == 2: p0=[float(MaxInt), 0.15, 1.2, 90, 10] elif ExpFunc1 == 3: p0=[float(MaxInt), 1.0, 80, 10] #Now do the fit p1, Success = optimize.leastsq(ErrFunc1, p0[:], \ args = (FullRad1d, FullInt1d),full_output=0,\ maxfev=100000.,\ Dfun = SpatialJacobian, col_deriv=1,\ xtol = 1.e-12, gtol=1e-15, ftol=1e-15) if ExpFunc1 == 1 or ExpFunc1 == 3: FitC0 = p1[0]; FitC1 = p1[1]; FitC2=p1[2]; FitC3=p1[3]; FitC4=-1 else: FitC0 = p1[0]; FitC1 = p1[1]; FitC2=p1[2]; FitC3=p1[3]; FitC4=p1[4] print '\n' print 'Success: ' + str(Success) #If Verbosity == 1 print the fit coefficients if Success == 1 or Success ==2: if ExpFunc1 == 1: Fit = p1[0]*cos(p1[1]*FitRad1d**p1[2])*exp(-FitRad1d**2/p1[3]) if Verbosity == 1: print 'P1:' +str(p1[0]) + ' , ' + str(p1[1])+ ' , '+ \ str(p1[2]) + ' , ' + str(p1[3]) elif ExpFunc1 == 2: Fit = p1[0]*cos(p1[1]*FitRad1d**p1[2])*exp(-FitRad1d**2/p1[3])+p1[4] if Verbosity == 1: print 'P1:' +str(p1[0]) + ' , ' + str(p1[1])+ ' , ' \ +str(p1[2]) + ' , ' + str(p1[3])+ ' , ' \ +str(p1[4]) elif ExpFunc1 == 3: Fit = p1[0]*cos((1./(2*pi))*FitRad1d**p1[1])*exp(-FitRad1d**2/p1[2])+\ p1[3] if Verbosity == 1: print 'P1:' +str(p1[0]) + ' , ' + str(p1[1])+ ' , '+ \ str(p1[2]) + ' , ' + str(p1[3]) elif ExpFunc1 == 4: Fit = p1[0]*cos(p1[1]*FitRad1d)*FitRad1d**(p1[2]) if len((where(isnan(Fit)))[0]) == 0: NoFit=0 MinFitInt = Fit.min() MaxFitInt = Fit.max() MinFitInd = (where(Fit == MinFitInt))[0] MaxFitInd = (where(Fit == MaxFitInt))[0] if len(MinFitInd) > 0: MinFitRad = FitRad1d[MinFitInd[0]] else: MinFitRad = -1 if len(MaxFitInd) > 0: MaxFitRad = FitRad1d[MaxFitInd[0]] else: MaxFitRad = -1 ##Find the zero crossings FitProd = Fit[0:len(Fit)-1]*Fit[1:len(Fit)] Criterium = (FitProd < 0) | (FitProd == 0) ZeroCrossings = where(Criterium)[0] if len(ZeroCrossings) > 2: FirstFitZRad = FitRad1d[ZeroCrossings[0]] SecondFitZRad = FitRad1d[ZeroCrossings[1]] else: FirstFitZRad = -1 SecondFitZRad = -1 else: NoFit = 1 else: NoFit = 1 else: NoFit = 1 #If the fit was not successful, give -1s if NoFit == 1 and AttFit == 1: MinRad = -1 MaxRad = -1 MinFit = -1 MaxFit = -1 MinFitInt = -1 MaxFitInt = -1 MinFitRad = -1 MaxFitRad = -1 p1=[-1,-1,-1,-1,-1] FitC0 = -1 FitC1 = -1 FitC2 = -1 FitC3 = -1 FitC4 = -1 FirstFitZRad = -1 SecondFitZRad = -1 Fit = arange(0,FullRad1d.max(),0.1) Fit[:] = 0. if AttFit == 0: MinFitInt = -1 MaxFitInt = -1 MinFitRad = -1 MaxFitRad = -1 p1=[-1,-1,-1,-1,-1] FitC0 = -1 FitC1 = -1 FitC2 = -1 FitC3 = -1 FitC4 = -1 FirstFitZRad = -1 SecondFitZRad = -1 Fit = arange(0,FullRad1d.max(),0.1) Fit[:] = 0. ##Find the zero crossings IntProd = IntArr1d[0:len(IntArr1d)-1]*IntArr1d[1:len(IntArr1d)] Criterium = (IntProd < 0) | (IntProd == 0) ZeroCrossings = where(Criterium)[0] if len(ZeroCrossings) > 1: FirstZRad = Rad1d[ZeroCrossings[0]] SecondZRad = Rad1d[ZeroCrossings[1]] ##Get the sum of halo and core FirstCrit = (FullRad1d <= FirstZRad) PixInCore = where(FirstCrit) SumCore = FullInt1d[PixInCore].sum() SecondCrit = (FullRad1d > FirstZRad) & (FullRad1d <= SecondZRad) PixInHalo = where(SecondCrit) SumHalo = FullInt1d[PixInHalo].sum() ThirdCrit = (FullRad1d > SecondZRad) PixOutHalo = where(ThirdCrit) SumOut = FullInt1d[PixOutHalo].sum() else: FirstZRad = -1 SecondZRad = -1 SumCore = -1 SumHalo = -1 SumOut = -1 ##Allow for plotting if TvFlag == 1: if AttFit == 1: mplot.figure(2) mplot.hold(False) mplot.plot(FitRad1d,Fit) mplot.hold(True) mplot.plot(FullRadii, FullIntArr, 'y.') mplot.plot(Rad1d,IntArr1d,'bo') time.sleep(2) else: FirstZRad = -1 SecondZRad = -1 SumCore = -1 SumHalo = -1 SumOut = -1 MinFit = -1 MaxFit = -1 MinFitInt = -1 MaxFitInt = -1 MinFitRad = -1 MaxFitRad = -1 p1=[-1,-1,-1,-1, -1] FitC0 = -1 FitC1 = -1 FitC2 = -1 FitC3 = -1 FitC4 = -1 FirstFitZRad = -1 SecondFitZRad = -1 if Verbosity == 1: print 'Read: ' +str(ReadNum)+' , '+'Sum: '+str(SumInt) print 'MinRad: ' +str(MinRad) +' , '+'Min: '+str(MinInt) print 'MaxRad: ' +str(MaxRad) +' , '+'Max: '+str(MaxInt) if NegPersistence > 0: if Verbosity == 1: print 'First Zero: ' +str(FirstZRad) print 'Second Zero:' +str(SecondZRad) print 'Sum Core: ' + str(SumCore) print 'Sum Halo: ' + str(SumHalo) print 'Sum Out : ' + str(SumOut) return MinInt, MaxInt, MinRad, MaxRad, FirstZRad, SecondZRad, SumCore,\ SumHalo, SumOut, SumInt, FitC0, FitC1, FitC2, FitC3, FitC4,\ MinFitInt, MaxFitInt, MinFitRad, MaxFitRad, FirstFitZRad, SecondFitZRad