pro LineFit2, 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 ;If there's no dark ramp, make some zeros If (Size(y2))[0] eq 0 then y2 = Fltarr(N_Elements(y1)) 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) ;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) ;Set the signma value to 3 sigma away from the line Threshold = 6*StdDev(Diff) ;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) 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 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 Mask = Mask(0:PossCos(0)) NRej = 1 EndElse EndIf Else begin ;No more cosmic rays were found NRej = 0 EndElse If TvFlag then begin plot, ValidX,ValidY, yrange=[min(ValidY), max(ValidY)] Oplot, ValidX, M*ValidX+B xyouts, 0.2,0.80,'M : '+Strtrim(M,2),/normal xyouts, 0.2,0.75,'B : '+Strtrim(B,2),/normal stop End ;Use Error as Gauge Error = StdDev((ValidY - (M*ValidX + B))) EndWhile end ;pro Plot_AllNights - When the input to Examine_Ref_Pixels is a list of ; files, the median of Ch0 will be plotted for all of the nights on ; one plot pro Plot_AllNights, FitsFiles, ReadNum = ReadNum Common KeyParams, KeyStr Common ColorStr, Color, ColorNum, ColorTri, ColorNumTri If KeyStr.TvFlag eq 2 then device,filename=KeyStr.PlotDir+'Ch0Noise.eps',$ /encapsul If N_Elements(ReadNum) eq 0 then ReadNum =0 Readcol, FitsFiles, RawFiles, format='A' !p.charsize=3 If KeyStr.TvFlag eq 1 then window,xsize=1000,ysize=1000 plot, [0], [0], /nodata, background = fsc_color('white'),$ color=fsc_color('black'),$ yrange=[-200,200], xrange=[0, KeyStr.Naxis2] ;!p.noerase=1 !p.multi=[0,2,2] For FileNum = 0, N_Elements(RawFiles)-1 do begin Fits_Read_DataCube, RawFiles(FileNum), Ch0, RawHeader, $ XStart = 0, XStop = KeyStr.Ch0W-1, YStart=0, $ YStop = KeyStr.Naxis2-1, ZStart = ReadNum, ZStop = ReadNum Ch0Med = Median(Ch0, dimension = 1) Ch0MedNoOffset = Ch0Med - Mean(Ch0) ;KPNO Data are different for different gain. Individual nights: If FileNum eq 0 or FileNum eq 2 then begin plot, Ch0MedNoOffset, color = fsc_color('black'), $ background = fsc_color('white'), $ xrange=[0,4100], xstyle=1, xtitle='Row', ytitle='ADU' ColorInd=1 EndIf else if FileNum eq 4 then begin plot, Ch0MedNoOffset, color = fsc_color('black'), $ background = fsc_color('white'),$ xtitle='Row', ytitle = 'ADU',$ xrange=[0,4100], xstyle=1 ColorInd=1 EndIf else begin ColorInd +=1 oplot, Ch0MedNoOffset, color = ColorNum.ColorNums(ColorInd) EndElse EndFor xyouts,0.55,0.425,'Upper Left', $ /normal, color=fsc_color('black'), charsize=3 legend,['07Apr24','07Apr25'], psym=[0,0],$ colors=[fsc_color('black'),ColorNum.ColorNums(2)], $ position = [0.7,0.45], thick = [6,6], /normal, $ textcolors=[fsc_color('black'),fsc_color('black')], $ charsize=3 xyouts,0.55,0.30,'Upper Right', /normal, color=fsc_color('black'), $ charsize=3 legend,['07Apr26','07Apr27'], psym=[0,0],$ colors=[fsc_color('black'),ColorNum.ColorNums(2)], $ position = [0.7,0.32], thick = [6,6], /normal, $ textcolors=[fsc_color('black'),fsc_color('black')], $ charsize=3 xyouts,0.55,0.18,'Lower Left', /normal, color=fsc_color('black'), $ charsize=3 legend,['07Apr28','07Apr29','07Apr30'], psym=[0,0,0],$ colors=[fsc_color('black'),ColorNum.ColorNums(2), $ ColorNum.ColorNums(3)], $ position = [0.7,0.20], thick=[6,6,6],/normal, $ textcolors=[fsc_color('black'),fsc_color('black'),fsc_color('black')],$ charsize=3 If KeyStr.TvFlag eq 3 then begin PngImg = tvrd() tvlct,reds,greens,blues,/get write_png,KeyStr.PlotDir+'Ch0NoiseAllNights.png', $ PngImg,reds,greens,blues EndIf !p.noerase=0 set_plot,'x' stop end ;pro Examine_Ch0 ; ;PURPOSE: ; To take an image and produce quantitative measure- ; ments of the reference pixels and channel 0 using the fast fourier ; transform. ; ;INPUTS: ; FITSFILENAME: ; The File Name for which the analysis should be done ; ;KEYWORDS: ; ROOTDIR: ; The Directory holding the directories of Raw files ; DATE: ; The Date during which they were taken. Also the subdirectory that ; the files are in. ; REDUCEDDIR: ; The directory where the processed images and plots will be placed in. ; BLOCKLENGTH: ; An offset for files that were written improperly. The readu ; function reads the values with this offset in place. ; OBJECTNAME: ; The string that indicates what files are being sought. The filenames ; should begin with this string ('Dark', 'Flat', etc.). ; NREADS: ; The number of reads (NAXIS3) that are present in the file. This ; is also a string in the filename that will be used as a regular ; expression in searching for matched files. ; FIRSTREAD: ; The first read to be used in the difference ; Diff = LastRead - FirstRead ; LASTREAD: ; The last read to be used in the difference ; THISFILTER: ; The filter through which the files were taken. ('I', 'G', or 'Y') ; TVFLAG ; 0 - Don't Plot anything ; 1 - Plot to the X Window ; 2 - Plot a postscript file that was opened by calling script ; 3 - Plot a png file tha can be converted for a paper ; AVGENTIREFRAME ; 0 - Don't do a fourier transform on entire frame ; 1 - Do a fourier transform on entire frame ; REMOVALMODE ; 0 - Try to remove noise by applying filter to fourier spectra ; 1 - Remove noise by subtracting median of Channel 0 from the rest ; of the array ; SKIPROW ; The number of rows to skip before plotting a slope if TVFLAG is ; set to 1. ; MODE: ; 1 - just read in the reference areas of the images ; 2 - read in areas and perform subtractions/slopefits, etc. ;OUTPUTS: ; ; TvFlag = 2, ; A postscript file will be written with all of the Fourier Plots ; TvFlag = 3 ; A .png file will be written with the fourier spectrum for channel 0 ;CALLING SEQUENCE: ; ; pro Examine_Ch0, FitsFiles, RootDir = RootDir, Date=Date, $ ReducedDir = ReducedDir, $ ObjectName=ObjectName, NReads = NReads, $ ReadNum = ReadNum, ThisFilter = ThisFilter,$ OverRide = OverRide, TvFlag = TvFlag, AvgEntireFrame = AvgEntireFrame, $ ExamineEachChannel = ExamineEachChannel, SubtractDark = SubtractDark, $ BiasSub = BiasSub, SkipRow = SkipRow, Mode = Mode If N_Elements(SkipRow) eq 0 then SkipRow = 0 If N_Elements(Mode) eq 0 then Mode = 1 Common KeyParams, KeyStr Common ColorStr, Color, ColorNum, ColorTri, ColorNumTri ;Get all of the defaults for the keywords from the KeywordStruct file @KeywordStruct.pro @PlotSettings.pro File_Mkdir, FourierDir If stregex(FitsFiles, '.lst', /boolean) then $ Plot_AllNights, FitsFiles ;POSTSCRIPT SET-UP********************************************************** If TvFlag eq 2 then begin device,filename=FourierDir+$ KeyStr.RawKey+'_FourierSpectra.ps' EndIf Case Mode of 1 : Begin Fits_Read_DataCube, FitsFiles, ImCh0, ImH, $ XStart=0, XStop = KeyStr.Ch0W-1 If Subtractdark then begin Fits_Read_DataCube, KeyStr.DarkName, DarkImCh0, DarkH, $ XStart=0, XStop = KeyStr.Ch0W-1 EndIf End 2 : Begin Fits_Read, FitsFiles, NoData, FitsHeader, /header_only FileParams = Return_File_Params(FitsFiles) SlopeStr = Return_Slope_Params(FitsHeader) ;Read the image that will be worked on ReadNum = 0 Im = long(Fits_Read_DataCube_Region(FitsFiles, FCB, $ FitsHeader, 0, KeyStr.Naxis1-1, 0, $ KeyStr.Naxis2-1, $ 0, KeyStr.TotalReads-1)) If Subtractdark then $ Dark = long(Fits_Read_DataCube_Region(KeyStr.DarkName, DarkFCB, $ DarkHeader, 0, KeyStr.Naxis1-1, 0, $ KeyStr.Naxis2-1, $ 0, KeyStr.TotalReads-1)) ;Subtract Bias Bias = long(Im(*,*,0)) If KeyStr.BiasSub then begin for ReadNum = 0, KeyStr.Naxis3-1 do begin Im(*,*,ReadNum) = long(Im(*,*,ReadNum))- Bias endfor EndIf If Not SubtractDark then begin RightRefMeds = Median(Im(KeyStr.Naxis1-4:KeyStr.Naxis1-1,*,*),dimension=1) Ch0RowMeds = Median(Im(0:KeyStr.Ch0W-1,*,*),dimension = 1) EndIf else begin ;right reference RightRefMeds = Median(Im(KeyStr.Naxis1-4:KeyStr.Naxis1-1,*,*),$ dimension=1)-$ Median(Dark(KeyStr.Naxis1-4:KeyStr.Naxis1-1,*,*),$ dimension=1) ;dark subtracted channel 0 Ch0RowDSMeds = Fltarr(KeyStr.Naxis2, KeyStr.Naxis3) Ch0RowDSMeds(*,1:KeyStr.Naxis3-1) = Median(($ Im(0:KeyStr.Ch0W-1,*,1:KeyStr.Naxis3-1)-$ Dark(0:KeyStr.Ch0W-1,*,1:KeyStr.Naxis3-1)),$ dimension=1) ;channel 0 no dark subtraction Ch0RowMeds = Median(Im(0:KeyStr.Ch0W-1,*,*),dimension = 1) ;the dark Dark(*,*,0) = 0 DarkRowMeds = Median(Dark(0:KeyStr.Ch0W-1,*,*),dimension = 1) EndElse RightRefMean = Mean(RightRefMeds) Ch0RowMean = Mean(Ch0RowMeds) TimePoints = SlopeStr.GroupTimes(1:KeyStr.Naxis3-1) RefRowSlopes = Fltarr(KeyStr.Naxis2) RefRowIntercepts = Fltarr(KeyStr.Naxis2) Ch0RowSlopes = Fltarr(KeyStr.Naxis2) Ch0RowIntercepts = Fltarr(KeyStr.Naxis2) Ch0RowDSSlopes = Fltarr(KeyStr.Naxis2) Ch0RowDSIntercepts = Fltarr(KeyStr.Naxis2) DarkRowSlopes = Fltarr(KeyStr.Naxis2) DarkRowIntercepts = Fltarr(KeyStr.Naxis2) TvFlag = 0 ;Slope Fit the Reference Rows - Doesn't really work too well. for Row = 0, KeyStr.Naxis2-1 do begin Mask = Indgen(KeyStr.Naxis3-1) ;slopefit reference rows LineFit2, TimePoints, RightRefMeds(Row,1:KeyStr.Naxis3-1), 0, $ Mask, m, b, sigm, sigb, error, 0 RefRowSlopes(Row)=m RefRowIntercepts(Row)=b If Row gt SkipRow and KeyStr.TvFlag then TvFlag = 1 If TvFlag eq 1 then print, Row ;slopefit channel 0 with no dark subtraction LineFit2, TimePoints, Ch0RowMeds(Row,1:KeyStr.Naxis3-1),0, $ Mask,m,b, sigm, sigb, error, TvFlag Ch0RowSlopes(Row)=m Ch0RowIntercepts(Row)=b ;slopefit channel 0 with dark subtraction LineFit2, TimePoints, Ch0RowDSMeds(Row,1:KeyStr.Naxis3-1),0, $ Mask,m,b, sigm, sigb, error, TvFlag Ch0RowDSSlopes(Row)=m Ch0RowDSIntercepts(Row)=b ;slopefit the dark channel 0 LineFit2, TimePoints, DarkRowMeds(Row,1:KeyStr.Naxis3-1),0, $ Mask,m,b, sigm, sigb, error, TvFlag DarkRowSlopes(Row)=m DarkRowIntercepts(Row)=b endfor fwindow = 64 order = 0 degree = 3 savgol_filter = savgol(fwindow, fwindow, order, degree) RefRowSmoothSlopes = convol(RefRowSlopes[4:KeyStr.Naxis2-5], $ savgol_filter, /EDGE_TRUNCATE) plot, RefRowSmoothSlopes, title = 'Reference Row Slopes' stop !p.multi=[0,1,2] Ch0RowSmoothSlopes = convol(Ch0RowSlopes[4:KeyStr.Naxis2-5], $ savgol_filter, /EDGE_TRUNCATE) plot, Ch0RowSmoothSlopes, title = 'Channel 0 Slopes' plot, Ch0RowSlopes, yrange=[-1,1] stop !p.multi=[0,1,2] Ch0RowDSSmoothSlopes = convol(Ch0RowDSSlopes[4:KeyStr.Naxis2-5], $ savgol_filter, /EDGE_TRUNCATE) plot, Ch0RowDSSmoothSlopes, title = 'Channel 0 Dark Subtracted Slopes' plot, Ch0RowDSSlopes, yrange=[-1,1] stop !p.multi=[0,1,1] Ch0RowSlopeDiff = Ch0RowDSSlopes-Ch0RowSlopes Ch0RowSlopeDiffSmoothSlopes = convol(Ch0RowSlopeDiff[4:KeyStr.Naxis2-5], $ savgol_filter, /EDGE_TRUNCATE) plot, Ch0RowSlopeDiffSmoothSlopes, $ title = 'Channel 0 Difference Smoothed Subtracted Slopes' stop plot, Ch0RowDSSlopes-Ch0RowSlopes, $ title = 'Difference of Ch0 with/without Dark' stop plot, Ch0RowDSSmoothSlopes - Ch0RowSmoothSlopes, $ title = 'Difference of Smooth Ch0 with/without Dark' stop DarkRowSmoothSlopes = convol(DarkRowSlopes[4:KeyStr.Naxis2-5], $ savgol_filter, /EDGE_TRUNCATE) plot, DarkRowSmoothSlopes, title = 'Dark Smoothed Slopes' stop plot, DarkRowSlopes, yrange=[-1,1], title= 'Dark Slopes' stop !p.multi=[0,1,1] ;Look at the effect of filtering each read for ReadNum = 0, KeyStr.NAxis3 - 1 do begin djs_iterstat, RightRefMeds(*,ReadNum), Mask = Good plot, RightRefMeds(*,ReadNum) - Mean(RightRefMeds(*,ReadNum)),$ title = 'Right Reference Rows for Read Number '+$ strtrim(ReadNum,2) RightRefMedsSmooth = convol(RightRefMeds(*,ReadNum), savgol_filter, $ /EDGE_TRUNCATE) stop plot, RightRefMedsSmooth-RightRefMean, $ title = 'Right Ref Rows for Read Number ' + $ strtrim(ReadNum,2) stop plot, Ch0RowDSMeds(*,ReadNum) - Ch0RowMean,$ title = 'Ch0 Rows for Read Number ' + $ strtrim(ReadNum,2) stop Ch0RowDSReadSmooth = convol(Ch0RowDSMeds(*,ReadNum), savgol_filter, $ /EDGE_TRUNCATE) plot,Ch0RowDSReadSmooth, title = 'Ch0 Rows for Read Number ' + $ strtrim(ReadNum,2) stop endfor ;Default back to the X Windowi set_plot,'x' End EndCase end