;pro LineFit ; ;An older version of the linefit routine used in the slopefitting ;algorithm. Not revised yet. pro LineFitb, x, y1, y2, Mask, m, b, sigm, sigb, Error, TvFlag ;Initialize Accumulators and Results ;Sx = 0 & Sxx = 0 ;Sy = 0 & Sxy = 0 ;m = 0.0 & b = 0.0 ;Sigm = 0.0 & Sigb = 0.0 ;Wt = 1.0 NRej=1 Found=0 Color=0 ThresHold = 300 While NRej eq 1 do begin ;Calclate the slopes from the points ValidX = double(X(Mask)) & ValidY1 = double(Y1(Mask)) ValidY = ValidY1 - double(Y2(Mask)) S = double(N_Elements(Mask)) Sx = Total(ValidX) Sy = Total(ValidY) Sxx = Total(ValidX^2) Sxy = Total(ValidX*ValidY) ;Plot the differences of the object frames RampDiff = y1(1:N_elements(y1)-1) - $ y1(0:N_elements(y1)-2) - $ (y2(1:N_elements(y2)-1)- $ y2(0:N_elements(y2)-2)) TimeDiff = x(1:N_elements(x)-1) - $ x(0:N_elements(x)-2) print, 'Diff :', RampDiff/TimeDiff print, 'Mean :', Mean(RampDiff/TimeDiff) print, 'StdDev:', StdDev(RampDiff/TimeDiff) print, 'NumPoints:', S ;Assume bright objects will have a large difference for all of ;the consecutive reads in the linear A/D range and don't reject ;Get the slope values and the difference of each point from the fit Delta = S*Sxx - Sx*Sx m = (S*Sxy - Sx*Sy)/Delta b = (Sxx*Sy-Sx*Sxy)/Delta Diff = ValidY- (M*ValidX+B) Threshold = 3*StdDev(Diff) print, 'Diff from line:', Diff print, 'StdDev:', StdDev(Diff) ;Find Cosmic Rays where there is a negative to positive spike in the ;difference DiffDiff = Diff(1:S-1) - Diff(0:S-2) print, DiffDiff print, 'Mean:', Mean(DiffDiff) print, 'StdDev:', StdDev(DiffDiff) PossCos = where(DiffDiff gt Threshold) If PossCos(0) ne -1 then begin ;A cosmic ray was found, use slope before it hit If PossCos(0) eq 0 then begin oplot, validx, m*validx+b, color = 255.^(Color+1) NRej = 1 Mask = Mask(1:N_elements(Mask)-1) If N_Elements(Mask) le 2 then begin m = 0 & b = 0 & error = 0 return EndIf EndIf Else Begin Found = 1 Mask = Mask(0:PossCos(0)) NRej = 1 EndElse EndIf Else begin ;No more cosmic rays were found NRej = 0 Found=0 EndElse oplot, validx, m*validx+b, color = fsc_color('red') Color+=1 ;Use Error as Gauge Error = StdDev((ValidY - (M*ValidX + B))) EndWhile end ;**************************************************************************** ;pro WindowStartUp ; ;Start up the windows and widgets pro WindowStartUp Common ExState, State, Masks Common Images, DarkIm, ObjectIm, FlatIm State.WindowBase = widget_base(title='Window', column = 1, uvalue = 'WindowBase') State.MainDrawId = widget_draw(State.WindowBase, xsize=512, ysize = 512, frame = 2, $ uvalue = 'MainWindowGetPix', /motion_events, /button_events,$ keyboard_events=1) Widget_Control, State.WindowBase, Set_Uvalue = info, /No_Copy Widget_Control, /Realize, State.WindowBase Widget_Control, State.MainDrawId, get_value = Tmp_Value State.MainWindowID = Tmp_Value Widget_Control, State.MainDrawID, Event_Pro = 'Examine_Pixel_List_Event' XManager, 'Examine_Pixel_List', State.WindowBase,/no_block End ;**************************************************************************** ;pro Examine_Pixel_List_Event ; ;This routine handles the clicking of the buttons on the image pro Examine_Pixel_List_Event, Event Common ExState case Event.type of 0: begin if Event.PRESS eq 4 then begin device,/close print, 'File Closed' endif else begin print, event.x,event.y State.PixelXSub = floor(event.x * float(State.Width/512.)) State.PixelYSub = floor(event.y * float(State.Height/512.)) State.PixelX = State.PixelXSub + State.XStart State.PixelY = State.PixelYSub + State.YStart print, State.PixelX, State.PixelY PlotPixels, State.PixelXSub, State.PixelYSub, $ State.PixelX, State.PixelY, State.TvFlagLocal endelse end else: begin end endcase end ;**************************************************************************** ;pro PlotPixels, ; ;This routine handles the plotting on the second window pro PlotPixels, PixelXSub, PixelYSub, PixelX, PixelY, TvFlagLocal Common ExState Common Images Common KeyParams, KeyStr Common SlopeParams, SlopeStr, SlopeFitStr If TvFlagLocal ne 2 then !p.multi=[0,2,2] ;Get the true pixel value PixTInd = PixelX + PixelY*State.Naxis2 ;Check to see if it falls in one of the pixel mask categories If (where(Masks.DeadInd eq PixTInd))[0] ne -1 then begin PixelType = 'Dead' EndIf else if (where(Masks.LeakyInd eq PixTInd))[0] ne -1 then begin PixelType = 'Leaky' Endif else if (where(Masks.LeakyNNInd eq PixTInd))[0] ne -1 then begin PixelType = 'Next to Leaky' Endif else if (where(Masks.HotInd eq PixTInd))[0] ne -1 then begin PixelType = 'Hot' Endif else begin PixelType = 'Good' EndElse If KeyStr.TvFlag eq 1 then begin resetplt !p.charsize=1 !p.charthick=1 EndIf Else If KeyStr.TvFlag eq 3 then begin Set_plot,'z' device,set_resolution=[1400,1200] !p.charsize=2 !p.charthick=2.5 !x.thick=1.75 !y.thick=1.75 !p.thick=2.25 !p.noerase=0 loadct, 0 !p.font = 1 device, font = 'Helvetica Bold' EndIf If PixelY gt State.Ref/2 then begin ;Subtract Biases ObjectRamp=long(ObjectIm(PixelXSub, PixelYSub, 1:State.NReads-1)) DarkRamp =long(DarkIm(PixelXSub,PixelYSub,1:State.NReads-1)) If TvFlagLocal ne 2 then begin ;Plot the Dark Image Ramp========================1 Plot, SlopeStr.GroupTimes, DarkIm(PixelXSub, PixelYSub, *), $ title = 'Dark Image for '+strtrim(PixelX,2)+','+$ strtrim(PixelY,2), $ yrange=[min(DarkIm(PixelXSub,PixelYSub,*)), $ max(DarkIm(PixelXSub,PixelYSub,*))], $ color = fsc_color('black'), background = fsc_color('white') ;Plot the Object Image Ramp=======================2 Plot, SlopeStr.GroupTimes, ObjectIm(PixelXSub,PixelYSub,*), $ title = 'Object Image for '+strtrim(PixelX,2)+','+strtrim(PixelY,2), $ yrange= [min(ObjectIm(PixelXSub,PixelYSub,*)),$ max(ObjectIm(PixelXSub,PixelYSub,*))], $ color=fsc_color('black') xyouts, 0.75,0.6, 'Bias Read : ' + Strtrim(ObjectIm(PixelXSub,PixelYSub,0)), $ /normal, $ color = fsc_color('black') xyouts, 0.75,0.65, 'Last Read : ' + $ Strtrim(ObjectIm(PixelXSub,PixelYSub,State.NReads - 1)), /normal, $ color=fsc_color('black') xyouts, 0.75, 0.70, 'Pixel Type : ', /normal, color = fsc_color('black') xyouts, 0.85, 0.70, PixelType, /normal, color = fsc_color('red') ;Plot the Object - Dark===========================3 Plot, FlatIm(PixelXSub, PixelYSub, *), $ title = 'Flat', $ color = fsc_color('black') EndIf ;Plot the Object - Dark - Bias ===================4 Plot, SlopeStr.GroupTimes, ObjectRamp - DarkRamp, $ yrange=[min(ObjectRamp-DarkRamp),max(ObjectRamp-DarkRamp)],$ title = 'Object - Dark : No Bias '+$ strtrim(PixelX,2)+','+strtrim(PixelY,2), $ color = fsc_color('black') ;Get the linear fit to the object -ramp Mask = where(ObjectRamp gt KeyStr.LowerLinLim and $ ObjectRamp lt KeyStr.UpperLinLim) If N_Elements(Mask) gt 2 then begin LineFit, SlopeStr.GroupTimes, ObjectRamp, DarkRamp, Mask, $ m, b, sigm, sigb, Error, TvFlag EndIf Else begin m=0 & b=0 & Error = 0 EndElse ;Output some numbers ==================================== xyouts, 0.75, .20, 'Slope : ' + $ Strtrim(m,2),$ /normal, color=fsc_color('black') xyouts, 0.75, .15, 'Intercept :'+ $ Strtrim(b,2), $ /normal, color=fsc_color('black') xyouts, 0.75, .10, 'Error :'+ $ Strtrim(Error,2), $ /normal, color = fsc_color('black') ;Plot the differences of the object frames Diff = ObjectRamp(1:N_elements(ObjectRamp)-1) - $ ObjectRamp(0:N_elements(ObjectRamp)-2) - $ (DarkRamp(1:N_elements(DarkRamp)-1)- $ DarkRamp(0:N_elements(DarkRamp)-2)) DiffDiff = Diff(1:N_elements(Diff)-1) - $ Diff(0:N_elements(Diff)-2) ;Plot, Diff, $ ; color = fsc_color('black') ;Plot, DiffDiff, $ ; color = fsc_color('black') If KeyStr.TvFlag eq 3 then begin PlotName = 'PixelSlopesFor_'+Strtrim(PixelX,2)+'_'+Strtrim(PixelY,2)+'.png' PngImg = tvrd() tvlct,reds,greens,blues,/get write_png,KeyStr.PlotDir+PlotName, $ PngImg,reds,greens,blues KeyStr.TvFlag = 1 Set_plot,'x' resetplt EndIf EndIf End ;;pro examine_pixel_list ; ;PURPOSE: ; To take in set of pixels and plot their individual ramps for a ; dark and an object file. ; ;INPUTS: ; PIXELS - A one dimensional list of pixel indices into the 2-d array ; (1-d coordinates). The pixel values of those coordinates will ; be plotted ; If pixels is not supplied, the xstart, xstop, ystart, and ystop ; keywords will be used to define a region ;KEYWORDS: ; ; PIXELSKIP - A number of pixels that will be skipped to reach farther back into ; the list ;CALLING SEQUENCE ; Examine_Pixel_List, objectfile='/nfs/slac/g/ki/ki09/lances/07Apr26/M13_Dither9_20ArcSecondEast_20ArcSecondsSouth_H4RG_SIPIN_30_Reads_Apr27_2007_03_11_35.fits*',/gui ; Examine_Pixel_List, objectfile='/nfs/slac/g/ki/ki09/lances/07Apr27/Landolt110956_Dither7_0ArcSecondsEast_20ArcSecondsSouth_H4RG_SIPIN_10_Reads_Apr28_2007_02_53_25.fits*',/gui ; Examine_Pixel_List, objectfile='/nfs/slac/g/ki/ki09/lances/07Apr27/Landolt110956_Dither7_0ArcSecondsEast_20ArcSecondsSouth_H4RG_SIPIN_10_Reads_Apr28_2007_02_40_47.fits',/gui,xstart=1650,ystart=1850 pro Examine_Pixel_List, Pixels, TvFlag = TvFlag, XStart = XStart, $ XStop = XStop, YStart=YStart, YStop = YStop , $ PixelSkip = PixelSkip, DarkFile = DarkFile, $ ObjectFile = ObjectFile, Date = Date, NReads = NReads, $ GUI = GUI, PlotIm = PlotIm, SlopeFit = SlopeFit, $ TvFlagLocal = TvFlagLocal If N_Elements(PixelSkip) eq 0 then PixelSkip = 0 If N_Elements(GUI) eq 0 then GUI = 0 If N_Elements(PlotIm) eq 0 then PlotIm = 0 If N_Elements(SlopeFit) eq 0 then SlopeFit = 0 If N_Elements(TvFlagLocal) eq 0 then TvFlagLocal = 1 ;Coords of bright stars in the default image If N_Elements(XStart) eq 0 then XStart = 1656 If N_Elements(XStop) eq 0 then XStop = XStart+50 If N_Elements(YStart) eq 0 then YStart = 1850 If N_Elements(YStop) eq 0 then YStop = YStart+50 ;Set up some state variables Common ExState Common Images Common KeyParams, KeyStr Common SlopeParams, SlopeStr, SlopeFitStr ;Include all the keywords from the KeywordStruct.pro file @KeywordStruct.pro @PlotSettings.pro State = { WindowBase: 0L, $ MainDrawID: 0L, $ MainWindowID: 0L, $ Naxis1 : Naxis1, $ Naxis2 : Naxis2, $ NumPixels : 0L, $ PixelX : 0L, $ PixelY : 0L, $ PixelXSub : 0L, $ PixelYSub : 0L, $ XStart : XStart, $ YStart : YStart, $ XStop : XStop, $ YStop : YStop, $ Height : 0L, $ Width : 0L, $ Ref : Ref ,$ NReads : NReads , $ SlopeFit : SlopeFit, $ TvFlagLocal : TvFlagLocal } set_plot,'x' Userwindow = !d.window wset, Userwindow ;Initialize the widget WindowStartUp ;A region was supplied to the script or a list of pixels If (Size(Pixels))[0] eq 0 then begin State.Width = long(XStop-XStart+1) & State.Height = long(YStop-YStart+1) State.NumPixels = long(State.Width)*long(State.Height) ;Get the relative coordinates of the pixels in the window PixelsSub = (Indgen(State.NumPixels) mod State.Width) + $ (Indgen(State.NumPixels)/ State.Width)*State.Width PixelsXSub = PixelsSub mod State.Width & PixelsYSub=PixelsSub/State.Width ;Get the absolute value of the pixels on the full frame Pixels = (Indgen(State.NumPixels) mod State.Width)+XStart+$ ((Indgen(State.NumPixels)/State.Width)+YStart)*(Naxis2) PixelsX= Pixels mod Naxis1 & PixelsY = Pixels/Naxis2 EndIf Else begin State.NumPixels = N_Elements(Pixels) XStart = 0 & XStop = Naxis1-1 YStart = 0 & YStop = Naxis2-1 PixelsXSub = Pixels mod Naxis1 & PixelsYSub = Pixels/Naxis2 PixelsX = PixelsXSub & PixelsY = PixelsYSub EndElse ;Open the original object if it's a slope fit If Not keyword_set(ObjectFile) then $ ObjectFile = '/nfs/slac/g/ki/ki09/lances/07Apr28/'+$ 'M13_Dither3_0ArcSecondsEast_20ArcSecondsNorth_H4RG_SIPIN_20_Reads_Apr29_2007_02_47_59.fits' DisplayFile='' ;Set the display plot to the slope fit and the object to be posted If stregex(ObjectFile, 'Slope', /boolean) then begin DisplayFile = ObjectFile ObjectFile = strmid(ObjectFile, 0,stregex(ObjectFile,'_Slope'))+$ '.fits' ObjectFileParts = strsplit(ObjectFile,'/',/extract, count=nparts) ObjectFile = RootDir+Date+PathDelim+ObjectFileParts(NParts-1) EndIf ;Read in the Object File Position = Fits_Open_DataCube(ObjectFile, ObjectFCB, ObjectHeader, 'Read') FileParams = Return_File_Params(ObjectFile) SlopeStr = Return_Slope_Params(ObjectHeader) State.NReads=KeyStr.TotalReads ObjectIm = Fits_Read_DataCube_Region(ObjectFile, ObjectFCB, ObjectHeader, $ XStart, XStop, YStart, YStop, 0, KeyStr.TotalReads-1) ObjectClose = Fits_Close_DataCube(ObjectFCB) Position = Fits_Open_DataCube(ObjectFile, ObjectFCB, ObjectHeader, 'Read') ;Set the display plot to the slope fit and the object to be posted If stregex(DisplayFile, 'Slope', /boolean) then begin Fits_Read, DisplayFile, DisplayIm DisplayIm = DisplayIm(XStart:XStop,YStart:YStop) EndIf Else begin DisplayIm = long(ObjectIm(*,*,State.NReads-1)) - $ long(ObjectIm(*,*,0)) EndElse print, 'Opened Object' ;Get some files for the plotting If Not keyword_set(DarkFile) then DarkFile = KeyStr.DarkName If Not keyword_set(FlatFile) then FlatFile = KeyStr.FlatName ;Read in the files that will be used for the examination Position = Fits_Open_DataCube(DarkFile, DarkFCB, DarkHeader, 'Read') DarkIm = Fits_Read_DataCube_Region(DarkFile, DarkFCB, DarkHeader, $ XStart, XStop, YStart, YStop, 0,KeyStr.TotalReads-1) DarkClose = Fits_Close_DataCube(DarkFCB) print, 'Opened Dark' ;Read in the Flat File Position = Fits_Open_DataCube(FlatFile, FlatFCB, FlatHeader, 'Read') FlatIm = Fits_Read_DataCube_Region(FlatFile, FlatFCB, FlatHeader, $ XStart, XStop, YStart, YStop, 0, 1) FlatClose = Fits_Close_DataCube(FlatFCB) print, 'Opened Flat' wset, state.mainwindowid DisplayIm = BytScl(DisplayIm, $ min = Mean(DisplayIm) - 2*StdDev(DisplayIM), $ max = Mean(DisplayIm) + 2*StdDev(DisplayIm)) tvscl, congrid(DisplayIm, 512,512) ;tvscl, congrid(ObjectIm(*,*,6)/double(FlatIm(*,*,1)),512,512) window,1,xsize=500, ysize = 500 ;Allow for plotting of the image being viewed to a .png If PlotIm eq 1 then begin Set_plot,'z' device,set_resolution=[512,512] DisplayIm = Fltarr(State.width, State.Height,3) DisplayIm = [[[ObjectIm(*,*,State.NReads-1)]], $ [[ObjectIm(*,*,State.NReads-1)]], $ [[ObjectIm(*,*,State.NReads-1)]]] tvimage, Congrid(bytscl(DisplayIm),512,512,3) ,/tv, true = 3 PlotName = 'M13Region.png' PngImg = tvrd() tvlct,reds,greens,blues,/get write_png,KeyStr.PlotDir+PlotName, $ PngImg,reds,greens,blues KeyStr.TvFlag = 1 Set_plot,'x' resetplt EndIf ;Get the bad pixel maps to check what type of pixel we're plotting Masks = Return_Pixel_Masks() ;If a list was supplied then plot the pixels in that list If GUI eq 0 then begin For PixInd = 0L, State.NumPixels-1 do begin ;Allow for skipping of some pixels PixTInd = PixInd + PixelSkip PixelX = PixelsX(PixTInd) & PixelY = PixelsY(PixTInd) PixelXSub = PixelsXSub(PixTInd) & PixelYSub = PixelsYSub(PixTInd) PlotPixels, PixelXSub, PixelYSub, PixelX, PixelY stop EndFor EndIf If State.TvFlagLocal eq 2 then begin set_plot,'ps' device, filename = KeyStr.PlotDir+KeyStr.ObjectName+'Slopes.ps' loadct, 39 device, /color ;allow color on the postscript device,ysize=8.5,/inches ;Height of plot in y device,xsize=11.0,/inches ;Width of plot in x device,yoffset=0.50,/inches ;Y position of lower left corner device,xoffset=0.00,/inches !p.multi=[0,1,1] !p.charsize=1 !p.thick=3 !p.charthick=3 Endif End