; LineFit : 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 ; ; y1 = Raw Signal ; y2 = Dark Current ; x = Time at beginning of read ; Mask = Array subscripts - begin as full array - end with indices ; of reads flagged as hit by cosmic rays remvoed from array ; BitMas = An array of 0s or 1s where the 1s are at the indices ; of Mask ; m = slope of fit ; b = y intercept of fit ; sigm = weight of points going into m calc ; ; Error = Error in fit from data points ; TvFlag = Set to 1 in order to plot slopes ; CosmicEDep = If a possible cosmic ray is detected, this will be the ; value of I(read+1)-I(read) where read is the read during ; which the cosmic ray hit ; ; NOTE: If TvFlag = 1 and the user wants to save the plot, he/she should ; enter saveplot=1 after the program has stopped and the enter ; .continue. pro LineFit, x, y1, y2, Mask, m, b, sigm, sigb, Error, CosmicEDep, $ Col, Row, ThresholdFactor , TvFlag, TvPPT ;Local Variables CosmicEDep = 0 TvFlagLocal = 0 NRej=1 NumOrig = N_Elements(Mask) LastInd = Mask(N_Elements(Mask)-1) BitMask = BytArr(NumOrig)+1 ;Assume slope of 0 in case no fit is possible m = 0 & b = 0 & error = 0 SavePlot = 0 NIter = 0 ;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)) ValidY = ValidY1 - double(Y2(Mask)) 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 = 1.*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*SpikeThresh and $ Diff(PossCos(0)-1) lt -1.75*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 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 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)) EndElse EndElse EndIf Else begin ;No more cosmic rays were found NRej = 0 EndElse If TvFlagLocal then begin If NIter eq 0 then begin erase screensize= get_screen_size() window, 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 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') Oplot, ValidX, M*ValidX+B, color =fsc_color('black'),$ linestyle = 2, thick = 6 stop EndIf else begin 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') stop if SavePlot then begin 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 end