;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009 ;pro tv_1_window_wcentroid ; ;PURPOSE: ; To calculate the centroid ;INPUTS: ; FitsFileName - The path to the file that contains the naxis1 x naxis2 x ; number of reads x number of windows data cube ; pro Tv_EpsLyra_Part_file, FitsFileName, TvFlag = TvFlag, $ PlotCentroid = PlotCentroid If N_Elements(PlotCentroid) eq 0 then PlotCentroid =0 If N_Elements(TvFlag) eq 0 then TvFlag = 0 Common FitsCube, FCB, unit @KeywordStruct.pro @PlotSettings.pro ;Open the FitsFile to the first valid data byte after the header. ;The common file control block FCB holds the structure to the fits file Position = fits_open_datacube(FitsFileName, FCB, Header, 'Read') naxis1 = fcb.axis(0) naxis2 = fcb.axis(1) naxis3 = fcb.axis(2) naxis4 = fcb.axis(3) bitpix = fcb.bitpix(0) bzero = fcb.bzero(0) bytes_per_word = (abs(bitpix)/8) ;Zero the bad pixels in each window Masks = Return_Pixel_Masks() SubMask = BytArr(KeyStr.Naxis1,KeyStr.Naxis2) SubMasks = BytArr(Naxis4, Naxis1, Naxis2) SubFlats = Fltarr(Naxis4, Naxis1, Naxis2) WindowLimits = LonArr(Naxis4, 4) ;Get the Flat and the subsections FileParams= Return_File_Params(FitsFileName) Fits_Read_DataCube, KeyStr.FlatName, Flat, FHeader, $ ZStart=1, ZStop =1 For WindowNum=0, Naxis4-1 do begin WindowLimits(WindowNum, 0) = (XMin =$ Long(SxPar(Header, 'WNDW'+Strtrim(WindowNum,2)+'XMI'))) WindowLimits(WindowNum, 1) = (XMax =$ Long(SxPar(Header, 'WNDW'+Strtrim(WindowNum,2)+'XMA'))) WindowLimits(WindowNum, 2) = (YMin = $ Long(SxPar(Header, 'WNDW'+Strtrim(WindowNum,2)+'YMI'))) WindowLimits(WindowNum, 3) = (YMax =$ Long(SxPar(Header, 'WNDW'+Strtrim(WindowNum,2)+'YMA'))) SubMask(XMin:XMax,YMin:YMax)=1 SubMask = SubMask and Masks.BadMask SubMasks(WindowNum,*,*) = ~SubMask(XMin:XMax-1, YMin:YMax-1) SubFlats(WindowNum,*,*) = Flat(XMin:XMax-1, YMin:YMax-1) EndFor ;Arrays for holding NumWindows = 4000 Star1Centroids = Fltarr(NumWindows,2) Star2Centroids = Fltarr(NumWindows,2) Star12Centroids = Fltarr(NumWindows,2) Star3Centroids = Fltarr(NumWindows,2) Star4Centroids = Fltarr(NumWindows,2) Time1Arr = DblArr(NumWindows) TIme2Arr = DblArr(NumWindows) ;Image stacking Arrays StackS = 50 Stack1Centroids = DblArr(NumWindows/StackS,2) Stack2Centroids = DblArr(NumWindows/StackS,2) TimeStackArr = DblArr(NumWindows/StackS,2) Stack1 = Dblarr(Naxis1,Naxis2) NumStack1 = Dblarr(Naxis1,Naxis2) Stack1Arr = Dblarr(Naxis1,Naxis2,NumWindows/StackS) Stack2 = Dblarr(Naxis1,Naxis2) NumStack2 = Dblarr(Naxis1,Naxis2) Stack2Arr = Dblarr(Naxis1,Naxis2,NumWindows/StackS) ;Window Time Tp = 10.e-6 ;Pixel time in seconds Tr = (Naxis1*Tp+18.e-6) ;Row time in seconds Tw = Naxis2*Tr ;Window time in seconds Tc = 2*3.*Tw+2*2*18.e-6+150.e-6 ;Cycle time in seconds ;WINDOW COORDINATES SS = 200 YRat = SS/Naxis1 & XRat = SS/Naxis2 X4Min = 18 & X4Max = 42 & X4Cen = (X4Min+X4Max)/2 Y4Min = 8 & Y4Max = 24 & Y4Cen = (Y4Min+Y4Max)/2 X1Min = Double(05) & X1Max = Double(23) & X1Cen = (X1Min+X1Max)/2 Y1Min = Double(20) & Y1Max = Double(45) & Y1Cen = (Y1Min+Y1Max)/2 X2Min = Double(30) & X2Max = Double(48) & X2Cen = (X2Min+X2Max)/2 Y2Min = Double(20) & Y2Max = Double(45) & Y2Cen = (Y2Min+Y2Max)/2 X3Min = Double(14) & X3Max = Double(42) & X3Cen = (X3Min+X3Max)/2 Y3Min = Double(26) & Y3Max = Double(48) & Y3Cen = (Y3Min+Y3Max)/2 X1Min0 = X1Min & Y1Min0 = Y1Min & X1Max0 = X1Max & Y1Max0 = Y1Max X2Min0 = X2Min & Y2Min0 = Y2Min & X2Max0 = X2Max & Y2Max0 = Y2Max X3Min0 = X3Min & Y3Min0 = Y3Min & X3Max0 = X3Max & Y3Max0 = Y3Max If KeyStr.TvFlag eq 1 then begin If PlotCentroid then begin window,0, xsize = 512, ysize = 512 window,1, xsize = 512, ysize = 512 Plot,[0, NumWindows],[-2,2],/nodata EndIf Else Begin window, xsize = 512, ysize =768 EndElse EndIf StartFrame =200. FirstFrame = 1. for FrameNum = StartFrame , NumWindows-1 do begin ; Naxis3 - 1 do begin ;Look at image stacks If FrameNum mod StackS eq 0 and FrameNum gt StartFrame+1 then begin X1Min = X1Min0 & X1Max = Y1Min0 & X1Cen = (X1Min+X1Max)/2 Y1Min = Y1Min0 & Y1Max = Y1Max0 & Y1Cen = (Y1Min+Y1Max)/2 Stack1Cen = Find_Centroid(Stack1, X1Min, X1Max, Y1Min, Y1Max, /global) Stack1Centroids(FrameNum/StackS,0) = Stack1Cen(0) Stack1Centroids(FrameNum/StackS,1) = Stack1Cen(1) Stack1Arr(*,*,FrameNum/StackS) = Stack1 X3Min = X3Min0 & X3Max = X3Max0 Y3Min = Y3Max0 & Y3Max = Y3Max0 Stack2Cen = Find_Centroid(Stack2, X2Min, X2Max, Y2Min, Y2Max, /global) Stack2Centroids(FrameNum/StackS,0) = Stack2Cen(0) Stack2Centroids(FrameNum/StackS,1) = Stack2Cen(1) Stack2Arr(*,*,FrameNum/StackS) = Stack2 TimeStackArr(FrameNum/StackS)=Tc*(FrameNum+WindowNum) Stack1 = Dblarr(Naxis1,Naxis2) NumStack1 = Dblarr(Naxis1,Naxis2) Stack2 = Dblarr(Naxis1,Naxis2) NumStack2 = Dblarr(Naxis1,Naxis2) EndIf for windownum = 0, Naxis4 - 1 do begin Fits_Read_Datacube,FitsFileName, Data, Header, $ ZStart = Naxis4*FrameNum+WindowNum, $ ZStop = Naxis4*FrameNum+WindowNum ThisWindow = Data MinWindow = Min(ThisWindow) MeanWin = Mean(ThisWindow) StdDevWin = StdDev(ThisWindow) ThisWindow = Data * SubMasks(WindowNum,*,*) If KeyStr.TvFlag eq 1 then begin wset,0 tv, bytscl(congrid(thisWindow,SS,SS), $ min = MeanWin - 2*StdDevWin, $ max = MeanWin + 2*StdDevWin), $ (WindowNum mod 2)*(SS), $ (WindowNum/2)*(SS) EndIf position = position + long(naxis1)*long(naxis2)*bytes_per_word centroid=1 if centroid then begin If WindowNum eq 0 then begin ;Left Double Stars Stack1(*,*)= Stack1(*,*)+Double(ThisWindow) ;Attempt a recursive search for the centroid Star1Cen = [4,4] X1Min = X1Min0 & X1Max = X1Max0 & X1Cen = (X1Min+X1Max)/2 Y1Min = Y1Min0 & Y1Max = Y1Max0 & Y1Cen = (Y1Min+Y1Max)/2 Star1Cen = Find_Centroid(ThisWindow, X1Min, X1Max, $ Y1Min, Y1Max, /global) Star1 = ThisWindow(X1Min:X1Max, Y1Min:Y1Max) Star2Cen = [4,4] X2Min = X2Min0 & X2Max = X2Max0 & X2Cen = (X2Min+X2Max)/2 Y2Min = Y2Min0 & Y2Max = Y2Max0 & Y2Cen = (Y2Min+Y2Max)/2 Star2Cen = Find_Centroid(ThisWindow, X2Min, X2Max, $ Y2Min, Y2Max, /global) Star2 = ThisWindow(X2Min:X2Max, Y2Min:Y2Max) X1Circ = Star1Cen(0)*XRat Y1Circ = Star1Cen(1)*YRat X2Circ = Star2Cen(0)*XRat Y2Circ = Star2Cen(1)*YRat DStar1 = Star1 DStar2 = Star2 If FirstFrame eq 1 then begin X1StartCirc = X1Circ Y1StartCirc = Y1Circ X2StartCirc = X2Circ Y2StartCirc = Y2Circ ;X12StartCirc = X12Circ ;Y12StartCirc = Y12Circ X3StartCirc = 0 Y3StartCirc = 0 X4StartCirc = 0 Y4StartCirc = 0 EndIf Time1Arr(FrameNum) = Tc*(FrameNum+WindowNum) Star1Centroids(FrameNum,0) = Star1Cen(0) Star1Centroids(FrameNum,1) = Star1Cen(1) Star2Centroids(FrameNum,0) = Star2Cen(0) Star2Centroids(FrameNum,1) = Star2Cen(1) EndIf Else Begin Stack2(*,*) = Stack2(*,*)+Double(ThisWindow) ;Right Double Stars X3Min = X3Min0 & X3Max = X3Max0 Y3Min = Y3Min0 & Y3Max = Y3Max0 Star3Cen = Find_Centroid(ThisWindow, X3Min, X3Max, $ Y3Min, Y3Max, /global) Star3 = ThisWindow(X3Min:X3Max, Y3Min:Y3Max) Star4 = (DStar2 = ThisWindow(X4Min:X4Max,Y4Min:Y4Max)) Star4Min = Min(Star4(where(Star4 ne 0))) Star4Cen = Return_Centroid(Star4-Star4Min) Star3Centroids(FrameNum,0) = Star3Cen(0) Star3Centroids(FrameNum,1) = Star3Cen(1) Star4Centroids(FrameNum,0) = Star4Cen(0) Star4Centroids(FrameNum,1) = Star4Cen(1) X1Circ = (Star3Cen(0))*XRat Y1Circ = (Star3Cen(1))*YRat X2Circ = (X4Cen+Star4Cen(0))*XRat Y2Circ = (Y4Cen+Star4Cen(1))*YRat DStar1 = Star3 Time2Arr(FrameNum) = Tc*(FrameNum+WindowNum) If FirstFrame eq 1 then begin X3StartCirc = X1Circ Y3StartCirc = Y1Circ X4StartCirc = X2Circ Y4StartCirc = Y2Circ FirstFrame = 0 EndIf EndElse Endif If KeyStr.TvFlag eq 1 then begin wset,0 tv, bytscl(congrid(DStar1,200,200), $ min = MeanWin - 2*StdDevWin, $ max = MeanWin + 2*StdDevWin), $ 0, $ 200 tv, bytscl(congrid(DStar2,200,200), $ min = MeanWin - 2*StdDevWin, $ max = MeanWin + 2*StdDevWin), $ 200, $ 200 If FrameNum mod StackS gt 1 then begin tv, bytscl(congrid(Stack1,200,200), $ min = Min(Stack1(where(Stack1 ne 0)))),$ 0, $ 600 tv, bytscl(congrid(Stack2,200,200), $ min = Min(Stack2(where(Stack2 ne 0)))),$ 200, $ 600 EndIf wait,0.01 tvcircle, 3, X1Circ+(WindowNum mod 2)*SS,$ Y1Circ, fsc_color('green'),/fill tvcircle, 3, X2Circ+(WindowNum mod 2)*SS,$ Y2Circ, fsc_color('green'),/fill tvcircle, 3, X1StartCirc, Y1StartCirc, fsc_color('red'),/fill tvcircle, 3, X2StartCirc, Y2StartCirc, fsc_color('red'),/fill tvcircle, 3, X3StartCirc+ SS, Y3StartCirc, fsc_color('red'),/fill tvcircle, 3, X4StartCirc+ SS, Y4StartCirc ,fsc_color('red'),/fill print, FrameNum If PlotCentroid then begin wset,1 oplot, [FrameNum], [Star1Cen(0)], psym=2 EndIf EndIf endfor endfor ;Get the relative motions with respect to the previous measurements RelCen1X = Star1Centroids(1:NumWindows-1,0)-Star1Centroids(0:NumWindows-2,0) RelCen1Y = Star1Centroids(1:NumWindows-1,1)-Star1Centroids(0:NumWindows-2,1) RelCen2X = Star2Centroids(1:NumWindows-1,0)-Star2Centroids(0:NumWindows-2,0) RelCen2Y = Star2Centroids(1:NumWindows-1,1)-Star2Centroids(0:NumWindows-2,1) RelCen12X = Star12Centroids(1:NumWindows-1,0)-Star12Centroids(0:NumWindows-2,0) RelCen12Y = Star12Centroids(1:NumWindows-1,1)-Star12Centroids(0:NumWindows-2,1) RelCen3X = Star3Centroids(1:NumWindows-1,0)-Star3Centroids(0:NumWindows-2,0) RelCen3Y = Star3Centroids(1:NumWindows-1,1)-Star3Centroids(0:NumWindows-2,1) RelCen4X = Star4Centroids(1:NumWindows-1,0)-Star4Centroids(0:NumWindows-2,0) RelCen4Y = Star4Centroids(1:NumWindows-1,1)-Star4Centroids(0:NumWindows-2,1) ;Get the relative motions with respect to the first measurement AbsCen1X = Star1Centroids(1:NumWindows-1,0)-Star1Centroids(0,0) AbsCen1Y = Star1Centroids(1:NumWindows-1,1)-Star1Centroids(0,1) AbsCen2X = Star2Centroids(1:NumWindows-1,0)-Star2Centroids(0,0) AbsCen2Y = Star2Centroids(1:NumWindows-1,1)-Star2Centroids(0,1) AbsCen12X = Star12Centroids(1:NumWindows-1,0)-Star12Centroids(0,0) AbsCen12Y = Star12Centroids(1:NumWindows-1,1)-Star12Centroids(0,1) AbsCen3X = Star3Centroids(1:NumWindows-1,0)-Star3Centroids(0,0) AbsCen3Y = Star3Centroids(1:NumWindows-1,1)-Star3Centroids(0,1) AbsCen4X = Star4Centroids(1:NumWindows-1,0)-Star4Centroids(0,0) AbsCen4Y = Star4Centroids(1:NumWindows-1,1)-Star4Centroids(0,1) If KeyStr.TvFlag eq 2 then begin device, filename = KeyStr.PlotDir+'EpsilonLyraCentroids.eps',/encapsul EndIf DiffMotion12X = AbsCen1X-AbsCen2X DiffMotion13X = AbsCen1X-AbsCen3X DiffMotion12Y = AbsCen1Y-AbsCen2Y DiffMotion13Y = AbsCen1Y-AbsCen3Y !P.Noerase=0 Erase ArcAx = 0.125*(findgen(15)-7) Plot,Time1Arr(200:900)-Time1Arr(200), $ abscen1x(200:900)-mean(abscen1x(200:900)), $ xrange=[Time1Arr(0),Time1Arr(700)], xstyle=1,$ yrange=[-7,7], ystyle = 9, $ background = fsc_color('white'), color=fsc_color('black'), $ xtitle = 'Time (s)', ytitle = 'Pixels', position = [0.075,0.1,0.55,.9],$ /normal axis,yaxis=1,color=fsc_color('black'), yrange=[min(ArcAx),max(ArcAx)], $ ytitle='Arcseconds' oplot, Time2Arr(200:900)-Time1Arr(200), $ abscen3x(200:900)-mean(abscen3x(200:900)), $ color = fsc_color('red') items = ['Eps02 Lyr C', 'Eps01 Lyr A'] pcolors= [fsc_color('black'), fsc_color('red')] ptextcolors = [fsc_color('black'), fsc_color('black')] psyms = [0,0] pthick = [4,4] pcharsize = [1.5,1.5] legend, items, colors = pcolors, textcolors=ptextcolors, psym=psyms, $ position = [0.275,0.255],/normal, thick = pthick, charsize=2 stack12=Fltarr(Naxis1,Naxis2,3) stack12(*,*,0) = stack2arr(*,*,5) stack12(*,*,1) = stack2arr(*,*,5) stack12(*,*,2) = stack2arr(*,*,5) stack11=Fltarr(Naxis1,Naxis2,3) stack11(*,*,0) = stack1arr(*,*,5) stack11(*,*,1) = stack1arr(*,*,5) stack11(*,*,2) = stack1arr(*,*,5) stack22=Fltarr(Naxis1,Naxis2,3) stack22(*,*,0) = stack2arr(*,*,10) stack22(*,*,1) = stack2arr(*,*,10) stack22(*,*,2) = stack2arr(*,*,10) stack21=Fltarr(Naxis1,Naxis2,3) stack21(*,*,0) = stack1arr(*,*,10) stack21(*,*,1) = stack1arr(*,*,10) stack21(*,*,2) = stack1arr(*,*,10) stack32=Fltarr(Naxis1,Naxis2,3) stack32(*,*,0) = stack2arr(*,*,17) stack32(*,*,1) = stack2arr(*,*,17) stack32(*,*,2) = stack2arr(*,*,17) stack31=Fltarr(Naxis1,Naxis2,3) stack31(*,*,0) = stack1arr(*,*,17) stack31(*,*,1) = stack1arr(*,*,17) stack31(*,*,2) = stack1arr(*,*,17) tvimage, bytscl(stack11,min=min(stack11(where(stack11 ne 0)))),$ position = [0.67,0.65,0.82,0.9], /normal tvimage, bytscl(stack12,min=min(stack12(where(stack12 ne 0)))),$ position = [0.845,0.65,0.995,0.9], /normal xyouts, 0.85,0.675, $ 't = '+Strtrim(string(TimeStackArr(5)-TimeStackArr(5),format='(F10.4)'),2)+$ ' s',$ color=fsc_color('red'),/normal,charsize=1.75 tvimage, bytscl(stack21,min=min(stack21(where(stack21 ne 0)))),$ position = [0.67,0.375,0.82,0.625], /normal tvimage, bytscl(stack22,min=min(stack22(where(stack22 ne 0)))),$ position = [0.845,0.375,0.995,.625], /normal xyouts, 0.85,0.40, $ 't = '+Strtrim(string(TimeStackArr(10)-TimeStackArr(5),format='(F10.4)'),2)+$ ' s',$ color=fsc_color('red'),/normal,charsize=1.75 tvimage, bytscl(stack31,min=min(stack31(where(stack31 ne 0)))),$ position = [0.67,0.1,0.82,0.35], /normal tvimage, bytscl(stack32,min=min(stack32(where(stack32 ne 0)))),$ position = [0.845,0.1,0.995,.35], /normal xyouts, 0.85,0.125, $ 't = '+Strtrim(string(TimeStackArr(17)-TimeStackArr(5),format='(F10.4)'),2)+$ ' s',$ color=fsc_color('red'),/normal, charsize=1.75 If KeyStr.TvFlag eq 2 then device,/close stop If KeyStr.TvFlag ne 2 then stop Plot, DiffMotion13Y, background=fsc_color('white'), $ color=fsc_color('black'), title = 'Differential Centroid Motion', $ xtitle = 'Seconds', ytitle = 'Centroid Difference' OPlot, DiffMotion12Y, color=fsc_color('red') If KeyStr.TvFlag ne 2 then stop Plot, DiffMotion13X, background=fsc_color('white'), $ color=fsc_color('black'), title = 'Differential Centroid Motion', $ xtitle = 'Seconds', ytitle = 'Centroid Difference' OPlot, DiffMotion12X, color=fsc_color('red') stop end