;pro slopefit_region ; ;PURPOSE: ; ; To fit a slope to each of the pixels in a small region. This ; serves as a test of the slope fitting algorithm being used on ; the full frame. ; ; This routine should properly handle saturated pixels in a multi-read ; exposure with the exception of pixels that got saturated on the first read. ; ;INPUTS: ; RawFileName - The name of the file to be used ; ; ; XCenter - The center x coordinate ; XHWidth - The width of the box (Length = 2xHWidth + 1) ; YCenter - The center y coordinate ; YHWidth - The height of the box(Length = 2xHWidth + 1) ; ;OUTPUTS: ; SlopeFitStr - A structure containing several arrays described below ; all the arrays are (2 * XHWdith+1) x (2 * YHWidth+1) ; SlopeFitStr.Slopes - Slopes calculated (m in y=mx+b) ; SlopeFitStr.Intercepts - Intercepts calculated (b in y = mx+b) ; SlopeFitStr.Error - Error in the fit to the pixel values ; SlopeFItStr.SN - Read Number of sample points used to calculate ; SlopeFitStr.SatPix - A mask of pixels that were saturated ; SlopeFitStr.DeadPix - A mask of pixels that did not integrate ; SlopeFitStr.Mask10Sig - A mask of pixels that were 10Sigma outliers ; SlopeFitStr.CosmicMask - A mask of pixels that were never fit ; SlopeFitStr.CosmicRays - Difference values between two reads for pixels ; that were flagged as cosmic rays ; ;REFRENCES: ; Based on code written by Figer et. al at STSCI ;**************************************************************************** function Return_Slopefit_Region, FileName_or_Cube, $ XCenter = XCenter, XHWidth = XHWidth, $ YCenter = YCenter, YHWidth = YHWidth, $ DarkCube = DarkCube, $ SubtractDark = SubtractDark, FilterStr = FilterStr, $ NReads = NReads, $ TvFlag = TvFlag, Rem60Hz = Rem60Hz, Verbose = Verbose, $ RawHeader = RawHeader If not keyword_set(Verbose) then Verbose = 0 If not keyword_set(SubtractDark) then SubtractDark = 0 Common KeyParams, SlopeStr, KeyStr, SlopeFitStr ;Check to see if input is a filename or a datacube S = Size(FileName_or_Cube) & DataType = S[S[0]+1] ;Either 7 - A string, or an array of some type (float, double, etc.) Case DataType of 7: begin ;Open the DataCube to be read out later Position = Fits_Open_DataCube(RawFile, RawFCB, $ RawHeader, 'Read') ;Get axes and keywords related to the ramp sequence Width = 2*XHWidth+1 Height = 2*YHWidth+1 XMin = XCenter - XHWidth XMax = XCenter + XHWidth YMin = YCenter - YHWidth YMax = YCenter + YHWidth Naxis1 = RAWFcb.axis(0) Naxis2 = RAWFcb.axis(1) TotalReads = RAWFcb.axis(2) ;Number of Reads ;Check if File is a dark itself, if it is, don't subtract If stregex(RawFile, 'Dark', /boolean) then SubtractDark =0 If SubtractDark then begin DarkFileName = File_Search(DarkDir+'MedianDark*ForEachOf*'+$ Strtrim(SlopeStr.NumGroups,2)+'Reads'+FilterStr+'.fits') If DarkFileName(0) eq '' then $ DarkFileName = File_Search(DarkDir+'MedianBias*ForEachOf*'+$ Strtrim(SlopeStr.NumGroups,2)+'Reads'+FilterStr+'.fits') If DarkFileName(0) eq '' then stop Position2 = Fits_Open_DataCube(DarkFileName, DarkFCB, DarkHeader, 'Read') EndIf ;Subtract 60 hz if requested If Rem60Hz then begin RawCube = lonarr(Width,Height,TotalReads) RawTempCube = long(Fits_Read_DataCube_Region( RawFile, RawFCB,$ RawHeader, $ 0, Naxis1-1, YMin, YMax, 0, TotalReads-1)) For ReadNum = 0, TotalReads -1 do begin RawTempCube(*,*,ReadNum) = Subtract_60Hz(RawTempCube(*,*,ReadNum)) RawCube(*,*,ReadNum) = RawTempCube(XMin:XMax, *, ReadNum) EndFor Endif Else begin RawCube = long(Fits_Read_DataCube_Region( RawFile, RawFCB,$ RawHeader, $ XMin, XMax, YMin, YMax, 0, TotalReads-1)) Frame=RawCube EndElse RawBiasFrame = RawCube(*,*,0) ;Get bias frame for dark exposure If SubtractDark then begin DarkCube = long(Fits_Read_DataCube_Region(DarkFileName, DarkFCB, $ DarkHeader, XMin, XMax, YMin, YMax, 0, TotalReads-1)) DarkBiasFrame = DarkCube(*,*,0) EndIf End Else: begin Width = S[1] Height = S[2] TotalReads = S[3] RawCube = FileName_Or_Cube RawBiasFrame = RawCube(*,*,0) If SubtractDark then DarkBiasFrame = DarkCube(*,*,0) End EndCase ;Get the slope parameters from the header SlopeStr = Return_Slope_Params(RawHeader) ;Masks and arrays for interesting things SatPix = bytArr(Width,Height) ;A mask that maps saturated pixels DeadPix = bytarr(Width,Height) ;A mask that maps dead pixels Mask10Sig = intarr(Width, Height) ;Values outside the 10 sigma range CosmicRays = UIntarr(Width, Height) ;Values of the cosmic ray hits CosmicMask = bytarr(Width, Height) ;A mask to eliminate bright objects ;Running Sums, Cross-Terms and Values according to 15.2 in Numerical Recipes ;Delta = S*Sxx - Sx^2 -- a = (SxxSy - SxSxy)/Delta -- b = (SSxy - SxSy)/Delta ;Number of valid points is stored in SN and the indices of these points is Sxy = dblarr(Width, Height) Sx = dblarr(Width, Height) Sy = dblarr(Width, Height) Sxx = dblarr(Width, Height) SN = bytarr(Width, Height, TotalReads) ;Make an array that will hold the previous frame for loop subtraction PreviousFrame = fltarr(Width, Height) For ReadNum = 1, TotalReads-1 do begin If Verbose then print, 'Working on Read ' + Strtrim(ReadNum,2) ;Get the time of this read in the exposure Time = float(SlopeStr.GroupTimes(ReadNum)) ;Get the current Raw frame RawFrame = RawCube(*,*,ReadNum) ;Find the saturated pixels and stop calculating slopes for these NewSatPix = where(RawFrame lt KeyStr.LowerLinLim or $ RawFrame gt KeyStr.UpperLinLim or $ CosmicMask gt 0, Complement= NotSatPix) If NotSatPix(0) eq -1 then break ;Subtract the bias Frame = RawFrame-RawBiasFrame ;Print some values ValidPixVals = Frame If Verbose then begin print, 'Number of Pixels in Nonlinear regime: ' + $ strtrim(N_elements(NewSatPix),2) print, 'Number of Pixels below zero after bias sub: '+ $ strtrim(N_elements(where(ValidPixVals lt 0)),2) EndIf ;Increment the dead pixel mask where the difference between the ;current read and the zeroth read is 0 DeadPixInd = where(Frame lt 0) If DeadPixInd(0) ne -1 then DeadPix(DeadPixInd) += 1 If SubtractDark then begin ;Get the Dark Frame DarkFrame = DarkCube(*,*,ReadNum) ;Subtract bias from dark (is this right?) DarkFrame = DarkFrame-DarkBiasFrame Frame = long(Frame) - long(DarkFrame) EndIf ;Get difference between this and last frame DiffFrame = Frame - PreviousFrame CRSDev = stddev(DiffFrame) CRMed = median(DiffFrame) Outliers = where (DiffFrame gt (CRMed + CRSDev*10)) If (Outliers[0] ne -1) then begin Mask10Sig[Outliers] = Mask10sig[Outliers] + 1 CosmicMask[Outliers] += 1 CosmicRays[Outliers] = DiffFrame[Outliers] Endif ;Increment the running sums Sx(NotSatPix) = Sx(NotSatPix) + Time Sy(NotSatPix) = Sy(NotSatPix) + Frame(NotSatPix) Sxy(NotSatPix) = Sxy(NotSatPix) + Time * Frame(NotSatPix) Sxx(NotSatPix) = Sxx(NotSatPix) + Time^2 SNThisFrame = BytArr(Width, Height) SNThisFrame(NotSatPix) = 1 SN(*,*,ReadNum)= SNThisFrame ;Update previous frame PreviousFrame = Frame EndFor If DataType eq 7 then begin ;Close the datacubes If SubtractDark then DarkClose=Fits_Close_DataCube(DarkFCB) RawClose=Fits_Close_DataCube(RawFCB) EndIf ;Calculate slope according to p. 662 in Numerical Recipes Delta = SlopeStr.NumGroups * Sxx - Sx * Sx SlopeArray = (SlopeStr.NumGroups * Sxy - Sx * Sy) / Delta InterceptArray = (Sxx*Sy-Sx*Sxy)/ Delta ;Remove objects that repeatedly had big difference value CosmicMaskInds = where(CosmicMask gt 1) If CosmicMaskInds(0) ne -1 then CosmicRays(CosmicMaskInds) = 0 ;************************************************DEAD PIXELS ;Write the pixels that should be working, but aren't to a mask DeadPix(CosmicMask) = 0 ;Saturated Pixels ;Calculate the errors ErrorArray = DblArr(Width, Height) For Row=0, Height - 1 do begin For Col=0, Width -1 do begin ErrorPts = Frame - SlopeStr.NumGroups*SlopeArray(Col,Row) +$ InterceptArray(Col,Row) ErrorArray(Col, Row) = StdDev(ErrorPts) EndFor EndFor ;Delete the old file control blocks Delvarx, RawFCB Delvarx, DarkFCB ;Form the Return Structure SlopeFitStr = { $ Slopes: SlopeArray, $ Intercepts: InterceptArray, $ Errors: ErrorArray, $ SN : SN, $ SatPix: SatPix,$ DeadPix: DeadPix, $ Mask10Sig: Mask10Sig, $ CosmicMask: CosmicMask, $ CosmicRays: CosmicRays $ } return, SlopeFitStr End