;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: ; 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) ;REFRENCES: ; Based on code written by Figer et. al at STSCI ;**************************************************************************** function Slopefit_Region, RawFileName, XCenter, XHWidth, YCenter, YHWidth, $ RootDir = RootDir, Date = Date, $ ReducedDir = ReducedDir, ObjectName = ObjectName, $ SubtractDark = SubtractDark, FilterStr = FilterStr, NReads = NReads, $ TvFlag = TvFlag, Rem60Hz = Rem60Hz ;Include all of the keywords from the KeywordStruct file @KeywordStruct.pro ;Move to the directory cd, RootDir + Date If Rem60Hz then FilterStr = '_60HzRemoved' If stregex(RawFileName(0), 'Dark', /boolean) then begin ObjectName = 'Dark' EndIf FileNum = 0 ;Get the current file that's being operated on RawFile = RawFileName(FileNum) ;Get the base of the filename RawFileBase = (strsplit(RawFile,'.',/extract))[0] ;Get the unique keys in the filename for new filenames RawFileKey = strmid(RawFile, stregex(RawFile,'Reads_')+6) ;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 ;Get the slope parameters from the header SlopeStr = Return_Slope_Params(RawHeader) ;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 ;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 Sxy = dblarr(Width, Height) Sx = dblarr(Width, Height) Sy = dblarr(Width, Height) Sxx = dblarr(Width, Height) ;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)) 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 ;Make an array that will hold the previous frame for loop subtraction PreviousFrame = fltarr(Width, Height) For ReadNum = 1, TotalReads-1 do begin 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.LowerLim or $ RawFrame gt KeyStr.UpperLim or $ CosmicMask gt 0, Complement= NotSatPix) ;Subtract the bias Frame = RawFrame-RawBiasFrame ;Print some values ValidPixVals = Frame 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) ;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 ;Display histogram if keyword is set If TvFlag then begin FrameHist= histogram(Frame, locations=PixLoc) plot,PixLoc, FrameHist,psym=10 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 ;Update previous frame PreviousFrame = Frame EndFor ;Close the datacubes If SubtractDark then DarkClose=Fits_Close_DataCube(DarkFCB) RawClose=Fits_Close_DataCube(RawFCB) ;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 ;Delete the old file control blocks Delvarx, RawFCB Delvarx, DarkFCB ;Plotting if TvFlag = 1 If TvFlag eq 1 then begin !p.multi = [0,2,2] NegSlopes = where(SlopeArray lt 1000, NumNeg) For Pix = 0L, NumNeg - 1 do begin ;Get the X and Y Coordinates for the pixel w/r to the cube XPix = NegSlopes(Pix) mod Width & YPix = NegSlopes(Pix)/Height !p.multi = [0,2,2] ;Plot the Raw Image by default plot, RawCube(XPix, YPix,*), title = 'Object' If SubtractDark then begin ;Plot the uncorrected versions ObjMinDarkBiasCorrected = $ (RawCube(XPix, YPix,*) - RawBiasFrame(XPix,YPix)) - $ (DarkCube(XPix,YPix,*) - DarkBiasFrame(XPix,YPix)) plot, DarkCube(XPix, YPix, *), title = 'Dark' plot, RawCube(XPix, YPix, *) - DarkCube(XPix, YPix, *), title = 'Diff' plot, ObjMinDarkBiasCorrected, $ title = 'Diff - Bias' ErrorPts = ObjMinDarkBiasCorrected - $ SlopeStr.GroupTimes*SlopeArray(XPix,YPix)+$ InterceptArray(XPix,YPix) EndIf Else Begin ObjMinBias = RawCube(XPix,YPix,*) - RawBiasFrame(XPix,YPix) plot, ObjMinBias, title = 'Object - Bias' ErrorPts = ObjMinBias - $ SlopeStr.GroupTimes*SlopeArray(XPix,YPix)+$ InterceptArray(XPix,YPix) EndElse Error = StdDev(ErrorPts) oplot, GroupTimes*SlopeArray(XPix,YPix)+InterceptArray(XPix,YPix),$ linestyle = 2 xyouts, 0.2, .5, 'Slope : ' + strtrim(slopearray(XPix,YPix),2), $ /normal xyouts, 0.5, .5, 'Intercept: '+strtrim(InterceptArray(XPix,YPix),2), $ /normal xyouts, 0.8, .5, 'Error: ' + strtrim(Error, 2), $ /normal stop erase EndFor EndIf stop return, 1 End