;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009 ;NAME: ; pro LineFit_xrg ; ;PURPOSE: ; Compute a linear fit of the form y = mx+b to a set of data points ; x (Time Points), y1 (Counts from raw image), and y2(Counts from dark image). ; Code is copied from CNA_LINFIT from CALNICA software suite, which uses a ; slope calculation from Numerical Recipes. Code is adapted to use IDL ; functions ; ;INPUTS: (And OUTPUTS) ; x = Time at beginning of read ; y1 = Raw Signal ; m = slope of fit ; b = y intercept of fit ; sigm = weight of points going into m calc ; Error = Error in fit from data points ; CosmicEDep = If a possible cosmic ray is detected, this will be the ; value of I(read)-I(read-1) where read is the read during ; which the cosmic ray hit ; CosmicRead = The read in which the cosmic ray hit ; CosmicType = 0 - No cosmic was found ; 1 - An electrical spike likely ; 2 - A real cosmic ray likely ;KEYWORDS: ; ; THRESHOLDFACTOR - the sigma away from the line at which a point will ; be flagged as a cosmic ray ; Y2 = Dark Current or another signal that will be subtracted from y1 ; MASK = Array subscripts - begin as full array - end with indices ; of reads flagged as hit by cosmic rays remvoed from array ; TVFLAG = Set to 1 in order to plot slopes ; TVPPT 0 = Don't plot ; 1 = Plot Every Slope besides Reference Pixels ; 2 = Plot the cosmic rays ; 3 = Plot the electrical spikes ; WAITROW: ; If TVFLAG is set, the slope will not be plotted until Row >= WAITROW ; WAITCOL: ; If TVFLAG is set, the slope will not be plotted until Col >= WAITCOL ; ; NOTE: If TvFlag = 1 and the user wants to save the plot, he/she should ; enter saveplot=1 when the program has stopped and then enter ; .continue ; ;CALLING SEQUENCE - ; linefit_xrg, x, y1, m, b, sigm, sigb, Error, CosmicEDep, CosmicRead ; pro LineFit_xRG, x, y1, m, b, sigm, sigb, Error, $ CosmicEDep, CosmicRead, CosmicType, y2=y2, Mask = Mask, $ Col=Col, Row=Row, ThresholdFactor=ThresholdFactor , $ TvFlag = TvFlag, TvPPT = TvPPT, $ WaitRow = WaitRow, WaitCol = WaitCol If N_Elements(Col) eq 0 then Col = 0 If N_Elements(Row) eq 0 then Row = 0 If N_Elements(y2) eq 0 then y2 = Dblarr(N_Elements(y1)) If N_Elements(Mask) eq 0 then Mask= BytArr(N_Elements(y1))+1 If N_Elements(TvFlag) eq 0 then TvFlag = 0 If N_Elements(TvPPT) eq 0 then TvPPT = 0 If N_Elements(ThresholdFactor) eq 0 then ThresholdFactor = 3 If N_Elements(WaitRow) eq 0 then WaitRow = 0 If N_Elements(WaitCol) eq 0 then WaitCol = 0 ;Local Variables CosmicEDep = 0 CosmicRead = 0 CosmicType = 2 TvFlagLocal = 0 NRej=1 NumOrig = N_Elements(Mask) LastInd = Mask(N_Elements(Mask)-1) BitMask = BytArr(NumOrig)+1 XOrig = X ;Assume slope of 0 in case no fit is possible m = 0 & b = 0 & error = 0 SavePlot = 0 NIter = 0 ;Check to see if we are looking for plotting reference rows and cols if TvPPT eq 1 then begin if Row lt 4 or Col lt 4 then TvPPT = 0 endif ;Iterate until there are no further rejections While NRej eq 1 do begin ;Each Iteration updates BitMask. It stays same size. ;Mask contains the valid indices. It shrinks in size Mask = where(BitMask eq 1) ;Calclate the slopes from the points ValidX = double(X(Mask)) & ValidY1 = double(Y1(Mask)) ValidY2 = double(Y2(Mask)) ValidY = ValidY1 - ValidY2 S = double(N_Elements(Mask)) Sx = Total(ValidX) Sy = Total(ValidY) Sxx = Total(ValidX^2) Sxy = Total(ValidX*ValidY) ;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) ;Very bright stars might only have two valid points; treat differently ;than cosmics If NumOrig le 2 then begin PossCos = -1 EndIf Else Begin ;Find Cosmic Rays where there is a negative to positive spike in the ;difference. Form Difference of Distance from the line from the line DiffDiff = Diff(1:S-1) - Diff(0:S-2) ;Set the sigma value to ThresholdFactor sigma away from the line Threshold = ThresholdFactor*StdDev(DiffDiff) SpikeThresh = 3.*StdDev(Diff) PossCos = where(DiffDiff gt Threshold) EndElse ;Base slope on where in the ramp the cosmic ray hit If PossCos(0) ne -1 then begin ;If a large spike is present at the first read, use following reads If PossCos(0) eq 0 or PossCos(0) eq 1 or PossCos(0) eq 2 then begin NRej = 1 BitMask(0:Mask(PossCos(0))) = 0 If N_Elements(where(BitMask) eq 1) le 2 then begin m = m & b = b & error = error return EndIf EndIf Else Begin ;Look for large negative going spikes that show up in ;columns 3500-3700 or so If Diff(PossCos(0)) lt -2.5*SpikeThresh and $ Diff(PossCos(0)-1) lt -2.0*SpikeThresh and $ Diff(PossCos(0)-2) lt -1.5*SpikeThresh and $ Mask(PossCos(0)) lt LastInd then begin NegSpots = 1 & NegSpotInds = 0 ;Move backward in the ramp until a positive difference is found CosmicType = 1 While NegSpots eq 1 do begin If PossCos(0) - NegSpotInds eq 1 then begin return Endif else begin RampDiff = ValidY(PossCos(0) - NegSpotInds) - $ ValidY(PossCos(0) - NegSpotInds -1) if RampDiff lt 0 then begin NegSpots = 1 BitMask(Mask(PossCos(0))-NegSpotInds)=0 NegSpotInds += 1 endIf else begin NegSpots = 0 BitMask(Mask(PossCos(0))) = 0 endElse Endelse If (TvPPT eq 3 or TvPPT eq 1) and TvFlag then $ TvFlagLocal = 1 Endwhile EndIf Else Begin ;Now check for cosmic rays BitMask(Mask(PossCos(0))+1:N_Elements(BitMask)-1) = 0 CosmicRead = Mask(PossCos(0))+1 NRej = 1 If (TvPPT eq 2 or TvPPT eq 1) and TvFlag then $ TvFlagLocal=1 ; Put the difference in an image made up of cosmic rays If PossCos(0) lt N_Elements(ValidY)-1 then $ CosmicEDep = ValidY(PossCos(0)+1)-ValidY(PossCos(0)) CosmicRead = PossCos(0)+1 EndElse EndElse EndIf Else begin ;No more cosmic rays were found NRej = 0 EndElse ;If TvFlag = TVPPT = 1, plot everytime If TvFlag eq 1 and TvPPT eq 1 then TvFlagLocal = 1 ;Plotting If TvFlagLocal and Col ge WaitCol and Row ge WaitRow then begin If NIter eq 0 then begin OutputPlot = 0 If OutputPlot eq 1 then begin screensize= get_screen_size() window, 0, xsize=screensize(0)-100,$ ysize=screensize(1)-100 !p.charsize=3.5 !p.charthick=3.5 !p.thick=3 !x.thick=3.25 !y.thick=3.25 EndIf window, 0 erase plot, ValidX, ValidY, yrange=[min(ValidY), max(ValidY)],$ title = 'Slope for Row : '+Strtrim(uint(Row),2) + $ ', Col : '+Strtrim(uint(Col),2),$ background = fsc_color('white'), color=fsc_color('black'), $ xtitle = 'Seconds', ytitle = 'ADU',$ position = [0.175,0.1,0.95,0.95],/normal xyouts,0.55,0.4,'!3Slope = '+Strtrim(M,2), /normal, $ color=fsc_color('black') xyouts,0.55,0.3,' !4r!l!3DD!3!N =' +$ Strtrim(StdDev(DiffDiff),2),/normal,$ color=fsc_color('black') xyouts,0.55,0.2, 'Max DD = ' +Strtrim(Max(DiffDiff),2), $ /normal, color=fsc_color('black') Oplot, ValidX, M*ValidX+B, color =fsc_color('black'),$ linestyle = 2, thick = 6 stop EndIf else begin wset, 0 oplot, ValidX, M*ValidX+B, color=fsc_color('red'), $ linestyle = 2, thick =6 xyouts,0.8,0.4,Strtrim(M,2),/normal,$ color=fsc_color('red') xyouts,0.8,0.3,Strtrim(StdDev(DiffDiff),2),/normal,$ color=fsc_color('red') !p.noerase=0 If (where(y2 ne 0))[0] ne -1 then begin window, 1 erase plot, X, Y2, $ color=fsc_color('black'), background=fsc_color('white'), $ xtitle = 'Seconds', ytitle = 'Dark Current (ADU)' EndIf stop if SavePlot then begin wset, 0 PngImg = tvrd(true=1) tvlct,reds,greens,blues,/get write_png, 'CosmicAt_Row_'+ $ Strtrim(Row,2) +'_Col_'+Strtrim(Col,2)+'.png',$ PngImg,reds,greens,blues endif EndElse End ;Use Error as Gauge Error = StdDev((ValidY - (M*ValidX + B))) NIter+=1 EndWhile ;If no deviations were found, no cosmic or spike was present If NIter eq 1 then CosmicType = 0 end