;spot_ellipticity ;PURPOSE: ; This script will analyze the spots in each of the CWFS raw/*.fits files ; for the [date] directory and try to fit an ellipse to them. ; ; An overall net ellipticity will point to a wind effect that is described ; in the Image motion paper by Martin ; pro spot_ellipticity,$ Date=Date,File100Analysis=File100Analysis,$ MinFileNumInt=MinFileNumInt,$ MaxFileNumInt=MaxFileNumInt If (not keyword_set(Date)) then Date='10' If (not keyword_set(stats_directory)) then stats_directory='stats' If (not keyword_set(HalfSpot)) then HalfSpot = 5 If (not keyword_set(NumGS)) then NumGS=625 If (not keyword_set(ChopSize)) then ChopSize = 100 If (not keyword_set(PaMinAng)) then PaMinAng = -!pi/2 If (not keyword_set(File100Analysis)) then File100Analysis=0 If (not keyword_set(MinFileNumInt)) then MinFileNumInt=0 If (not keyword_set(MaxFileNumInt)) then MaxfileNumInt=895 RawDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$ '/ShackHartman/raw/' StatsDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$ '/ShackHartman/stats_1_31/' PostDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/final_paper_figs/' readcol,'/nfs/slac/g/ki/ki08/lsst/CPanalysis/Ellip_Wind_Pointing.txt',$ basename,uttime,dateEntry,frame_e,frame_e1,frame_e2,$ frame_theta,phi_array,WSpeedArr,WDirArr,$ PAltArr,PAzArr,xcenters,ycenters,xfcenters,yfcenters,$ format='a,a,f,f,f,f,f,f,f,f,f,f,f,f,f,f' if Date eq '10' then begin StartF=0 EndF=897 frame_e=frame_e(StartF:EndF) frame_theta=frame_theta(StartF:EndF) frame_e=frame_e(MinFileNumInt:MaxFileNumInt) frame_theta=frame_theta(MinFileNumInt:MaxFileNumInt) endif else if Date eq '11' then begin StartF=898 EndF=1706 frame_e=frame_e(StartF:EndF) frame_theta=frame_theta(StartF:EndF) frame_e=frame_e(MinFileNumInt:MaxFileNumInt) frame_theta=frame_theta(MinFileNumInt:MaxFileNumInt) endif else if Date eq '12' then begin StartF=1707 EndF=2601 frame_e=frame_e(StartF:EndF) frame_theta=frame_theta(StartF:EndF) frame_e=frame_e(MinFileNumInt:MaxFileNumInt) frame_theta=frame_theta(MinFileNumInt:MaxFileNumInt) endif nonfiducial=return_fiducial() !p.multi=[0,1,2] ;ARRAY OF RAW FILES StatsFiles=file_search(StatsDir+'*5132*sub_offsets.txt',count=NumStatsFiles) if Date eq '11' then begin StatsTFiles=strarr(NumStatsFiles) StatsTFiles(NumStatsFiles-300:NumStatsFiles-1)=StatsFiles(0:299) StatsTFiles(0:NumStatsFiles-301)=StatsFiles(300:NumStatsFiles-1) StatsFiles=StatsTFiles endif if File100Analysis eq 0 then begin StatsFiles=StatsFiles(MinFileNumInt:MaxFileNumInt) PostScriptFileName=PostDir+'Spots'+Date+'Start'+strtrim(MinFileNumInt,2)+$ 'End'+strtrim(MaxFileNumInt,2)+'.ps' NumStatsFiles=N_elements(StatsFiles) endif else begin PostScriptFileName=PostDir+'Spots_100_Average'+Date+'.ps' endelse EArr=fltarr(NumStatsFiles,NumGS) PaArr=fltarr(NumStatsFiles,NumGS) I11Arr=fltarr(NumStatsFiles,NumGS) I22Arr=fltarr(NumStatsFiles,NumGS) E1Arr=fltarr(NumStatsFiles,NumGS) E2Arr=fltarr(NumStatsFiles,NumGS) E1mean=fltarr(NumStatsFiles) E2mean=fltarr(NumStatsFiles) Emean=fltarr(NumStatsFiles) PaMean=fltarr(NumStatsFiles) Evmean=fltarr(NumStatsFiles) TotalFlux=fltarr(NumStatsFiles) set_plot,'ps' device,filename=PostScriptFileName loadct,39 ;load color table 39 device,/color ;allow color on the postscript device,ysize=10.5,/inches ;Height of plot in y device,xsize=12.0,/inches ;Width of plot in x device,yoffset=1.0,/inches ;Y position of lower left corner white='FFFFFF'x black='000000'x red='FF0000'x green='00FF00'x blue='FF0000'x !P.CHARSIZE=1.2 !P.THICK=4. !P.multi=[0,1,2] !P.multi=[0,1,3] !p.noerase=1 for FileNum = 0, NumStatsFiles -1 do begin readcol,StatsFiles(FileNum),$ Index,Xcoord,Xoffset,Ycoord,Yoffset,$ I11,I22,I12,d,e1,e2,e,pa,flux,$ format='i,f,f,f,f,x,x,x,f,f,f,f,f,f,f,f,f' I11(nonFiducial)=-10 I22(nonFiducial)=-10 e(nonFiducial)=-10 pa(nonFiducial)=-10 e1(nonFiducial)=-10 e2(nonFiducial)=-10 flux(nonFiducial)=-10 Index(nonFiducial)=-10 Xcoord(nonFiducial)=-10 Ycoord(nonFiducial)=-10 Spots=where(I11 ne -10) EArr(FileNum,*)=e PaArr(FileNum,*)=pa I11Arr(FileNum,*)=I11 I22Arr(FileNum,*)=I22 E1Arr(FileNum,*)=e1 E2Arr(FileNum,*)=e2 E1Mean(FileNum)=mean(E1(Spots)) E2Mean(FileNum)=mean(E2(Spots)) EMean(FileNum)=mean(e(Spots)) PaMean(FileNum)=.5*atan(E2Mean(FileNum),E1Mean(FileNum)) EvMean(FileNum)=1./N_elements(Spots)*$ sqrt(total(E1(Spots)^2)+total(E2(Spots)^2)) TotalFlux(FileNum)=total(flux(Spots)) If File100Analysis eq 0 then begin E100=EArr(FileNum,*) Pa100=PaArr(FileNum,*) I11100=I11Arr(FileNum,*) I22100=I22Arr(FileNum,*) E1100=E1Arr(FileNum,*) E2100=E2Arr(FileNum,*) Spots100=where(E100 ne -10) E1Mean100=fltarr(NumGS) E2Mean100=fltarr(NumGS) I11Mean100=fltarr(NumGS) I22Mean100=fltarr(NumGS) EMean100=fltarr(NumGS) PaMean100=fltarr(NumGS) for SpotNo = 0, NumGS -1 do begin I11Spot=I11100(*,SpotNo) I22Spot=I22100(*,SpotNo) E1Spot=E1100(*,SpotNo) E2Spot=E2100(*,SpotNo) ESpot=E100(*,SpotNo) SpotOnF=where(I11Spot ne -10) if SpotOnF(0) ne -1 then begin I11Mean100(SpotNo)=mean(I11Spot(SpotOnF)) I22Mean100(SpotNo)=mean(I22Spot(SpotOnF)) EMean100(SpotNo)=mean(ESpot(SpotOnF)) PaMean100(SpotNo)=0.5*atan(mean(E2Spot(SpotOnF)),$ mean(E1Spot(SpotOnF))) endif else begin I11Mean100(SpotNo)=-10 I22Mean100(SpotNo)=-10 EMean100(SpotNo)=-10 PaMean100(SpotNo)=-10 E1Mean100(SpotNo)=-10 E2Mean100(SpotNo)=-10 endelse endfor I11Mean100=[I11Mean100,10] I22Mean100=[I22Mean100,10] EMean100=[EMean100,0.4] PaMean100=[PaMean100,!pi/4] PlSpots=[Spots,625] XCoord=[XCoord,760] YCoord=[YCoord,760] !p.multi=[0,1,3] Partvelvec,I11Mean100(PlSpots),fltarr(N_elements(PlSpots)),$ XCoord(PlSpots),YCoord(PlSpots),$ yrange=[100,820],xrange=[180,880],$ xstyle=5,ystyle=5,position=[0.025,0.4,0.625,1.0],$ color=80 !p.noerase=1 Partvelvec,fltarr(N_elements(PlSpots)),I22Mean100(PlSpots)/5,$ XCoord(PlSpots),YCoord(PlSpots),$ yrange=[100,820],xrange=[180,880],$ xstyle=5,ystyle=5,position=[0.025,0.4,0.625,1.0],$ color=160 Partvelvec,EMean100(PlSpots)*cos(PaMean100(PlSpots)),$ EMean100(PlSpots)*sin(PaMean100(PlSpots)),$ XCoord(PlSpots),YCoord(PlSpots),$ yrange=[100,820],xrange=[180,880],$ xstyle=5,ystyle=5,position=[0.025,0.4,0.625,1.0],$ color=240 PaMin=PaMinAng ;THETA NoBins=20. PaBinWidth=!pi/NoBins Pamax=PaMin+(NoBins-1)*PaBinWidth PaHist=histogram(Pa100(Spots100),Nbins=NoBins,locations=LocPa,$ reverse_indices=PaInd,$ max=PaMax,min=PaMin) PaBinSize=(PaMax-PaMin)/(NoBins-1) Plot,LocPa+PaBinSize/2,PaHist,psym=10,$ title='Ellipticity Magnitude and Angle for Frame '+$ strtrim(StatsFiles(FileNum)),$ xtitle='Radians',ytitle='Number of Spots',$ position=[0.025,0.025,0.45,0.375] oplot,[Frame_theta(FileNum),Frame_theta(FileNum)],$ [0,1000],color=240 EMin=0 ;ELLIPTICITY NoEBins=20. EBinWidth=(1-Emin)/NoEBins EMax=EMin+(NoEBins)*EBinWidth EHist=histogram(E100(Spots100),Nbins=NoEBins,locations=LocE,$ reverse_indices=EInd,$ max=EMax,min=EMin) EBinsize=(EMax-EMin)/(NoEBins-1) Plot,LocE,EHist,psym=10,$ xtitle='E',ytitle='Number of Spots',$ position=[0.5,0.025,0.925,0.375] oplot,[Frame_E(FileNum),Frame_E(FileNum)],$ [0,1000],color=240 Plot,Pa100(Spots100),E100(Spots100),psym=2,$ xtitle='Theta',ytitle='e',$ position=[0.65,0.4,0.925,0.675] oplot,[Frame_theta(FileNum),Frame_theta(FileNum)],$ [0,1000],color=240 xyouts,.7,.95,'Total I(x,y)x^2 10 pixels^2',color=80,/norm xyouts,.7,.91,'Total I(x,y)y^2 10 pixels^2',color=160,/norm xyouts,.7,.87,'Ellipticity = 0.4',color=240,/norm xyouts,.7,.83,'Mean Spot E = '+strtrim(Mean(E100(Spots100)),2),$ color=0,/norm xyouts,.7,.79,'Mean WvFt E = '+strtrim(Frame_e(FileNum),2),$ color=240,/norm xyouts,.7,.75,'Mean Spot Ang = '+strtrim(Mean(Pa100(Spots100)),2),$ color=0,/norm xyouts,.7,.71,'Mean WvFt Ang = '+strtrim(Frame_theta(FileNum),2),$ color=240,/norm !p.noerase=0 endif else if File100Analysis eq 1 then begin if (((FileNum+1) mod ChopSize) eq 0) or $ FileNum eq (NumStatsFiles-1) then begin StartInd=FileNum-99 EndInd=FileNum E100=EArr(StartInd:EndInd,*) Pa100=PaArr(StartInd:EndInd,*) I11100=I11Arr(StartInd:EndInd,*) I22100=I22Arr(StartInd:EndInd,*) E1100=E1Arr(StartInd:EndInd,*) E2100=E2Arr(StartInd:EndInd,*) Spots100=where(E100 ne -10) E1Mean100=fltarr(NumGS) E2Mean100=fltarr(NumGS) I11Mean100=fltarr(NumGS) I22Mean100=fltarr(NumGS) EMean100=fltarr(NumGS) PaMean100=fltarr(NumGS) for SpotNo = 0, NumGS -1 do begin I11Spot=I11100(*,SpotNo) I22Spot=I22100(*,SpotNo) E1Spot=E1100(*,SpotNo) E2Spot=E2100(*,SpotNo) ESpot=E100(*,SpotNo) SpotOnF=where(I11Spot ne -10) if SpotOnF(0) ne -1 then begin I11Mean100(SpotNo)=mean(I11Spot(SpotOnF)) I22Mean100(SpotNo)=mean(I22Spot(SpotOnF)) EMean100(SpotNo)=mean(ESpot(SpotOnF)) PaMean100(SpotNo)=0.5*atan(mean(E2Spot(SpotOnF)),$ mean(E1Spot(SpotOnF))) endif else begin I11Mean100(SpotNo)=-10 I22Mean100(SpotNo)=-10 EMean100(SpotNo)=-10 PaMean100(SpotNo)=-10 E1Mean100(SpotNo)=-10 E2Mean100(SpotNo)=-10 endelse endfor I11Mean100=[I11Mean100,10] I22Mean100=[I22Mean100,10] EMean100=[EMean100,0.4] PaMean100=[PaMean100,!pi/4] PlSpots=[Spots,625] XCoord=[XCoord,760] YCoord=[YCoord,760] !p.multi=[0,1,3] Partvelvec,I11Mean100(PlSpots),fltarr(N_elements(PlSpots)),$ XCoord(PlSpots),YCoord(PlSpots),$ yrange=[100,820],xrange=[180,880],$ xstyle=5,ystyle=5,position=[0.025,0.4,0.625,1.0],$ color=80 !p.noerase=1 Partvelvec,fltarr(N_elements(PlSpots)),I22Mean100(PlSpots)/5,$ XCoord(PlSpots),YCoord(PlSpots),$ yrange=[100,820],xrange=[180,880],$ xstyle=5,ystyle=5,position=[0.025,0.4,0.625,1.0],$ color=160 Partvelvec,EMean100(PlSpots)*cos(PaMean100(PlSpots)),$ EMean100(PlSpots)*sin(PaMean100(PlSpots)),$ XCoord(PlSpots),YCoord(PlSpots),$ yrange=[100,820],xrange=[180,880],$ xstyle=5,ystyle=5,position=[0.025,0.4,0.625,1.0],$ color=240 PaMin=PaMinAng ;THETA NoBins=20. PaBinWidth=!pi/NoBins Pamax=PaMin+(NoBins-1)*PaBinWidth PaHist=histogram(Pa100(Spots100),Nbins=NoBins,locations=LocPa,$ reverse_indices=PaInd,$ max=PaMax,min=PaMin) PaBinSize=(PaMax-PaMin)/(NoBins-1) Plot,LocPa+PaBinSize/2,PaHist,psym=10,$ title='Ellipticity Magnitude and Angle for Frames '+$ strtrim(StartInd,2)+'-'+strtrim(EndInd,2),$ xtitle='Radians',ytitle='Number of Spots',$ position=[0.025,0.025,0.45,0.375] EMin=0 ;ELLIPTICITY NoEBins=20. EBinWidth=(1-Emin)/NoEBins EMax=EMin+(NoEBins)*EBinWidth EHist=histogram(E100(Spots100),Nbins=NoEBins,locations=LocE,$ reverse_indices=EInd,$ max=EMax,min=EMin) EBinsize=(EMax-EMin)/(NoEBins-1) Plot,LocE,EHist,psym=10,$ xtitle='E',ytitle='Number of Spots',$ position=[0.5,0.025,0.925,0.375] Plot,Pa100(Spots100),E100(Spots100),psym=2,$ xtitle='Theta',ytitle='e',$ position=[0.65,0.4,0.925,0.675] xyouts,.7,.9,'Total I(x,y)x^2 10 pixels^2',color=80,/norm xyouts,.7,.85,'Total I(x,y)y^2 10 pixels^2',color=160,/norm xyouts,.7,.8,'Ellipticity = 0.4',color=240,/norm xyouts,.7,.75,'Mean E = '+strtrim(Mean(E100(Spots100)),2),color=0,/norm xyouts,.7,.7,'Mean E Angle = '+strtrim(Mean(Pa100(Spots100)),2),color=0,$ /norm !p.noerase=0 endif endif endfor !p.multi=[0,1,2] plot,EMean,psym=1,title='Ellipticity',xtitle='File Number',$ yrange=[0,1] oplot,frame_e,psym=2,color=240 items_e=['Spots','Wvft'] psym_e=[1,2] colors_e=[0,240] legend,items_e,psym=psym_e,colors=colors_e plot,180*PaMean/!pi,psym=1,title='Ellipticity Angle',xtitle='File Number',$ ytitle='Degrees',yrange=[-90,90] oplot,180*Frame_theta/!pi,psym=2,color=240 !p.multi=[0,1,1] plot,frame_theta,frame_e,psym=2 plot,frame_e,TotalFlux,psym=2,xtitle='Ellipticity',ytitle='Total Flux' plot,frame_theta,TotalFlux,psym=2,xtitle='E Angle',ytitle='Total Flux' plot,TotalFlux,psym=4,title='Total Flux',ytitle='Sum of pixels',$ xtitle='Frame Number' plot,EvMean,psym=4,title='Vector Ellipticity' device,/close set_plot,'X' stop end