pro final_fft_all_nights_successive_filter,$ stats_directory=stats_directory,$ depth=depth,exam_date=exam_date,$ xwindspeed=xwindspeed,ywindspeed=ywindspeed,$ phase_dim=phase_dim,$ low=low,lcut=lcut If (not keyword_set(exam_date)) then exam_date='1_31' If (not keyword_set(stats_directory)) then stats_directory='stats_'+exam_date+'/' If (not keyword_set(depth)) then depth=0 If (not keyword_set(postscript_dir)) then postscript_dir='distorted_postscripts' If (not keyword_set(xFrameDrift)) then xFrameDrift=25 ;pixels/frame If (not keyword_set(yFrameDrift)) then yFrameDrift=25 ;pixels/frame If (not keyword_set(xWindSpeed)) then xWindSpeed=2;17 ;m/sec If (not keyword_set(yWindSpeed)) then yWindSpeed=0;4 ;m/sec If (not keyword_set(phase_dim)) then phase_dim=long(1024) If (not keyword_set(low)) then low=1.8 If (not keyword_set(lcut)) then lcut=10000.0 If (not keyword_set(MultiWindSpeeds)) then MultiWindSpeeds=1 ;OPTIONS: ********************************************************** ;Avg_sub=1 means subtract the average offset of a given frame ;Avg_sub=2 means set the kx=0,ky=0 componenet of FFT to zero ;Avg_sub=3 means subtract a scaled offset from the distorted grid avg_sub=1 ;fitted_only=1 means use only the fitted points on each frame ;fitted_only=2 means use only the non-fitted points one each frame ;fitted_only=3 means use only points that have non-zero value fitted_only=1 ;UseFid=1 apply the fiducial cut ;UseFid=0 don't apply the fiducial cut UseFid = 0 ;if MultiWindSpeeds eq 1 then begin ; xwinds=[-18.3 , -17.0 , -17.7 , -16.3 , -17.0, -18.1] ; ywinds=[ 6.6 , -8.2 , 1.7 , -5.6 , 5.6, -3.5] ;endif if MultiWindSpeeds eq 1 then begin xwinds=[-18.3 , -17.9 , -17.6, -17.3, -17.0, $ -17.7, -17.3 , -16.9, -16.6, -16.3, $ -17.0, -17.3 , -17.6, -17.9, -18.1] ywinds=[ 6.6 , 2.9, -0.8, -4.5, -8.2, $ 1.7 ,-0.2, -2.0, -3.8, -5.6, $ 5.6, 3.3, 1.0, -1.3, -3.5] endif ;LABELS: *************************************************************** if fitted_only eq 1 then begin label_string='_fitted_only_' endif else if fitted_only eq 2 then begin label_string='_not_fitted_only' endif else if fitted_only eq 3 then begin label_string='_non_zero_only' endif if avg_sub eq 0 then begin label_string=label_string+'no_sub_' endif else if avg_sub eq 1 then begin label_string=label_string+'avg_sub_' endif else if avg_sub eq 2 then begin label_string=label_string+'kx0_ky0_sub_' endif else begin label_string=label_string endelse if UseFid eq 0 then begin label_string = '_NoFiducial_' endif else begin label_string = '_Fiducial_' endelse xwind=strtrim(string(xWindSpeed,format='(%"%10.1f")'),2) ywind=strtrim(string(yWindSpeed,format='(%"%10.1f")'),2) lowinmeters=phase_dim*.17/low highinmeters=phase_dim*.17/lcut lowinmeters=strtrim(string(lowinmeters,format='(%"%10.2f")'),2) highinmeters=strtrim(string(highinmeters,format='(%"%10.4f")'),2) ;FRAMESPEEDS FOR DRAGGING PHASE SCREEN ACROSS CWFS xFrameDrift_m=xFrameDrift*.17 yFrameDrift_m=yFrameDrift*.17 FrameDrifts=[xFrameDrift,yFrameDrift] max_FrameDrift=max(FrameDrifts) ;NUMBER OF CELLS TO USE FOR AVERAGING FROM WINDSPEEDS exp_time=0.030 ; milliseconds numXcells=ceil(xWindSpeed*exp_time/0.17)+1; numYcells=ceil(yWindSpeed*exp_time/0.17)+1; if MultiWindSpeeds eq 0 then begin label_string=label_string+'low_'+lowinmeters+'_high_'+highinmeters+$ '_xwind_'+xwind+'_ywind_'+ywind endif else begin label_string=label_string+'low_'+lowinmeters+'_high_'+highinmeters+$ '6xwinds6ywinds' endelse postscript_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/'+exam_date+'_figs/' postscript_dir='/afs/slac/u/ki/lances/latex/' date=['10','11','12'] white='FFFFFF'x black='000000'x png=0 if png eq 1 then begin set_plot,'z' device, set_font='Courier' device,set_resolution=[2000,1600] !p.multi=[0,1,4] !p.charsize=6 !p.charthick=4.4 !x.thick=6 !y.thick=6 !p.thick=5 !p.noerase=0 plotsymsize=3 endif else begin set_plot,'ps' device,filename=postscript_dir+'EFilteredAllNightsPlusPhaseScreen'+$ label_string+'.ps' !p.multi=[0,1,4] loadct,39 ;load color table 39 device,/color ;allow color on the postscript device,ysize=7.5,/inches ;Height of plot in y device,xsize=7.0,/inches ;Width of plot in x device,yoffset=0.50,/inches ;Y position of lower left corner device,xoffset=0.00,/inches white='FFFFFF'x black='000000'x !P.CHARSIZE=1.4 !P.CHARTHICK=2 !P.THICK=4 plotsymsize=1 endelse nonfiducial=[1,2,3,4,5,6,7,8,9,13,14,17,18,19,20,21,22,23,24,25,$ 26,27,28,29,30,31,32,38,39,45,46,47,48,49,50,$ 51,52,53,54,55,63,64,71,72,73,74,75,$ 76,77,78,79,88,89,98,98,99,100,$ 101,102,103,113,114,123,124,125,$ 126,127,128,138,139,149,150,$ 151,152,163,164,174,175,$ 176,177,188,189,200,$ 201,210,213,214,225,$ 226,238,239,249,$ 251,261,262,263,264,265,266,$ 276,286,287,288,289,290,291,$ 301,302,303,304,305,306,307,308,309,$ 310,311,312,313,314,315,316,$ 317,318,319,320,321,322,323,324,325,$ 326,327,336,337,338,339,340,341,$ 351,361,362,363,364,365,366,$ 376,388,389,$ 401,413,414,425,$ 426,427,438,439,449,450,$ 451,452,454,463,464,473,474,475,$ 476,477,478,488,489,498,499,500,$ 501,502,503,504,513,514,522,523,524,525,$ 526,527,528,529,530,538,539,546,547,548,549,550,$ 551,552,553,554,555,556,563,564,570,571,572,573,574,575,$ 576,577,578,579,580,581,582,588,589,594,595,596,597,598,599,600,$ 601,602,603,604,605,606,607,608,609,610,613,614,617,618,619,620,621,$ 622,623,624,625]-1 ;Extra Index for the scale offset complement,findgen(626),nonfiducial,fiducial fitted=fiducial ;*****************************************OPTIONS***************** ;Avg_sub=1 means subtract the average offset of a given frame ;Avg_sub=2 means set the kx=0,ky=0 componenet of FFT to zero ;Avg_sub=3 means subtract a scaled offset from the distorted grid avg_sub=1 ;fitted_only=1 means use only the fitted points on each frame ;fitted_only=2 means use only the non-fitted points one each frame ;fitted_only=3 means use only points that have non-zero value fitted_only=1 ;CONSTANTS kmax=12 num=25 grid_spots=625 index=findgen(625) ;ARRAY REFERENCE 25x25=625 FOR GRID subx=2*indgen(num)-25 ;The x-reference numbers suby=indgen(num)-12 ;The y-reference numbers Date=['10','11','12'] ;START THE REAL DATA ANALYSIS AND GENERATE 3 PLOTS for i=0,2 do begin stats_directory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date(i)+$ '/ShackHartman/stats_'+exam_date+'/' ;ARRAYS OF FILES offset_files=file_search(stats_directory+'*5132*sub_offsets.txt',count=num_stats_files) ;HUGE DATA CUBE x_offset_hor=complexarr(num_stats_files,num,num) y_offset_hor=complexarr(num_stats_files,num,num) x_offset_ver=complexarr(num_stats_files,num,num) y_offset_ver=complexarr(num_stats_files,num,num) fittedspots_arr=intarr(num_stats_files,num*num) ;ffts xfft_array=complexarr(num_stats_files,num,num) yfft_array=complexarr(num_stats_files,num,num) ;Correlations timee1=fltarr(num_stats_files-depth,num,num) timee2=fltarr(num_stats_files-depth,num,num) timee1num=fltarr(num_stats_files-depth,num,num) timee2num=fltarr(num_stats_files-depth,num,num) timeeden=fltarr(num_stats_files-depth,num,num) timee1_arr=fltarr(num,num) timee2_arr=fltarr(num,num) frame_e1=fltarr(num_stats_files) frame_e2=fltarr(num_stats_files) frame_e=fltarr(num_stats_files) frame_thetap=fltarr(num_stats_files) frame_thetam=fltarr(num_stats_files) kx0_array=fltarr(num_stats_files) ky0_array=fltarr(num_stats_files) phi_array=fltarr(num_stats_files) r_array=fltarr(num_stats_files) ;Average Ellipticity frame_e1_av_h=fltarr(kmax) frame_e2_av_h=fltarr(kmax) frame_e_av_h=fltarr(kmax) frame_e_av_den_h=fltarr(kmax) frame_e1_av_l=fltarr(kmax) frame_e2_av_l=fltarr(kmax) frame_e_av_l=fltarr(kmax) frame_e_av_den_l=fltarr(kmax) ;Grid for partvelvec x_ref_col=12.5*subx+542 y_ref_col=25.0*suby+460 index=findgen(num*num) gridx=fltarr(num*num) gridy=fltarr(num*num) gridx(index)=x_ref_col(index mod num) gridy(index)=y_ref_col(index/num) ;Indices and such no_bins=8 a=findgen(13)+12 b=findgen(12) subn=[a,b] subm=subn readcol,offset_files(0),xcoord,ycoord,format='x,f,x,f,x' ;FILL THE DATA CUBES for file_number=0, num_stats_files-depth-1 do begin ;Take the results from the offset.txt files and fill the cube readcol,offset_files(file_number),x_off,y_off,format='X,X,f,X,f' ;Obtain the difference ordered in COLUMN, ROW between centroid and grid point fittedspots=where(x_off ne -10) fittedspots_arr(file_number,fittedspots)=1 x_off(where(x_off eq -10))=0 y_off(where(y_off eq -10))=0 x_offset=fltarr(num,num) y_offset=fltarr(num,num) x_offset(*,*)=x_off y_offset(*,*)=y_off if avg_sub eq 1 then begin xfit=where(x_offset ne 0) yfit=where(y_offset ne 0) x_offset(xfit)=x_offset(xfit)-mean(x_offset(xfit)) y_offset(yfit)=y_offset(yfit)-mean(y_offset(yfit)) endif else if avg_sub eq 2 then begin x_fft=complexarr(num,num) y_fft=complexarr(num,num) x_fft(*,*)=fft(x_offset) y_fft(*,*)=fft(y_offset) x_fft_re=complexarr(num,num) y_fft_re=complexarr(num,num) for n=0,num-1 do begin for m=0,num-1 do begin x_fft_re(subn(n),subm(m))=x_fft(n,m) y_fft_re(subn(n),subm(m))=y_fft(n,m) endfor endfor kx0_array(file_number)=x_fft_re(12,12) ky0_array(file_number)=y_fft_re(12,12) phi_array(file_number)=atan(y_fft_re(12,12),x_fft_re(12,12)) r_array(file_number)=sqrt(x_fft_re(12,12)^2+y_fft_re(12,12)^2) x_fft_re(12,12)=0 y_fft_re(12,12)=0 for n=0,num-1 do begin for m=0,num-1 do begin x_fft(n,m)=x_fft_re(subn(n),subm(m)) y_fft(n,m)=y_fft_re(subn(n),subm(m)) endfor endfor x_offset=fft(x_fft,/inverse) y_offset=fft(y_fft,/inverse) endif else if avg_sub eq 3 then begin xfit=where(x_offset ne 0) yfit=where(y_offset ne 0) x_offset(xfit)=x_offset(xfit)+xscale*xoff(xfit) y_offset(yfit)=y_offset(yfit)+yscale*yoff(yfit) endif ;Obtain the difference ordered in ROW, COLUMN x_offset_hor(file_number,*,*)=x_offset y_offset_hor(file_number,*,*)=y_offset x_offset_ver(file_number,*,*)=x_offset((index mod 25)*25+index/25) y_offset_ver(file_number,*,*)=y_offset((index mod 25)*25+index/25) xfft_array(file_number,*,*)=fft(x_offset_hor(file_number,*,*)) yfft_array(file_number,*,*)=fft(y_offset_hor(file_number,*,*)) endfor for kval=0,kmax-1 do begin ;Correlations timee1=fltarr(num_stats_files-depth,num,num) timee2=fltarr(num_stats_files-depth,num,num) timee1num_h=fltarr(num_stats_files-depth,num,num) timee2num_h=fltarr(num_stats_files-depth,num,num) timeeden_h=fltarr(num_stats_files-depth,num,num) timee1num_l=fltarr(num_stats_files-depth,num,num) timee2num_l=fltarr(num_stats_files-depth,num,num) timeeden_l=fltarr(num_stats_files-depth,num,num) timee1_arr=fltarr(num,num) timee2_arr=fltarr(num,num) ;for high pass filter frame_e1_h=fltarr(num_stats_files) frame_e2_h=fltarr(num_stats_files) frame_e_h=fltarr(num_stats_files) frame_thetap_h=fltarr(num_stats_files) frame_thetam_h=fltarr(num_stats_files) frame_e_den_h=fltarr(num_stats_files) ;for low pass filter frame_e1_l=fltarr(num_stats_files) frame_e2_l=fltarr(num_stats_files) frame_e_l=fltarr(num_stats_files) frame_thetap_l=fltarr(num_stats_files) frame_thetam_l=fltarr(num_stats_files) frame_e_den_l=fltarr(num_stats_files) x_offset_hor_calc_h=complexarr(num_stats_files,num,num) y_offset_hor_calc_h=complexarr(num_stats_files,num,num) x_offset_hor_calc_l=complexarr(num_stats_files,num,num) y_offset_hor_calc_l=complexarr(num_stats_files,num,num) for file_number=0, num_stats_files-depth-1 do begin ;GET READY TO USE ONLY THE FITTED SPOTS fittedspots=fittedspots_arr(file_number,*) xfit=where(fittedspots eq 1) notfit=where(fittedspots ne 1) ;Average the fft x_fft=complexarr(num,num) y_fft=complexarr(num,num) x_low_fft=complexarr(num,num) x_low_fft_re=complexarr(num,num) y_low_fft=complexarr(num,num) y_low_fft_re=complexarr(num,num) y_high_fft_re=complexarr(num,num) x_high_fft=complexarr(num,num) y_high_fft=complexarr(num,num) x_high_fft_re=complexarr(num,num) x_fft(*,*)=xfft_array(file_number,*,*) y_fft(*,*)=yfft_array(file_number,*,*) x_fft_re=complexarr(num,num) y_fft_re=complexarr(num,num) ;Rearrange indices for n=0,num-1 do begin for m=0,num-1 do begin x_fft_re(subn(n),subm(m))=x_fft(n,m) y_fft_re(subn(n),subm(m))=y_fft(n,m) endfor endfor ;Filter kx=0, ky=0 x_high_fft_re=x_fft_re y_high_fft_re=y_fft_re x_high_fft_re((12-kval):(12+kval),(12-kval):(12+kval))=0 y_high_fft_re((12-kval):(12+kval),(12-kval):(12+kval))=0 x_low_fft_re((12-kval):(12+kval),(12-kval):(12+kval))=$ x_fft_re((12-kval):(12+kval),(12-kval):(12+kval)) y_low_fft_re((12-kval):(12+kval),(12-kval):(12+kval))=$ y_fft_re((12-kval):(12+kval),(12-kval):(12+kval)) ;Rearrange for n=0,num-1 do begin for m=0,num-1 do begin x_low_fft(n,m)=x_low_fft_re(subn(n),subm(m)) y_low_fft(n,m)=y_low_fft_re(subn(n),subm(m)) x_high_fft(n,m)=x_high_fft_re(subn(n),subm(m)) y_high_fft(n,m)=y_high_fft_re(subn(n),subm(m)) endfor endfor x_offset_hor_calc_h(file_number,*,*)=fft(x_high_fft,/inverse) y_offset_hor_calc_h(file_number,*,*)=fft(y_high_fft,/inverse) x_offset_hor_calc_l(file_number,*,*)=fft(x_low_fft,/inverse) y_offset_hor_calc_l(file_number,*,*)=fft(y_low_fft,/inverse) test=0 if test eq 1 then begin set_plot,'x' xoff=fltarr(625) yoff=fltarr(625) xoff(*)=x_offset_hor_calc_h(file_number,*,*) yoff(*)=y_offset_hor_calc_h(file_number,*,*) partvelvec,xoff,yoff,$ gridx,gridy,yrange=[150,800],$ xrange=[200,850],position=[0.1,0.1,0.9,0.9],$ title='High Passed Average Offsets for all images for k='+$ strtrim(kval,2) xyouts,gridx(624)-100,gridy(624)+125,'length='+$ strtrim(xoff(624),2) wait,4 endif ;high pass filter timee1num_h(file_number,*,*)=x_offset_hor_calc_h(file_number,*,*)*$ x_offset_hor_calc_h(file_number+depth,*,*)-$ y_offset_hor_calc_h(file_number,*,*)*$ y_offset_hor_calc_h(file_number+depth,*,*) timee2num_h(file_number,*,*)=(x_offset_hor_calc_h(file_number,*,*)*$ y_offset_hor_calc_h(file_number+depth,*,*)$ +x_offset_hor_calc_h(file_number+depth,*,*)$ *y_offset_hor_calc_h(file_number,*,*)) timeeden_h(file_number,*,*)=(x_offset_hor_calc_h(file_number,*,*)^2+$ y_offset_hor_calc_h(file_number,*,*)^2) timee1temp_h=timee1num_h(file_number,*,*) timee2temp_h=timee2num_h(file_number,*,*) timeedtemp_h=timeeden_h(file_number,*,*) if fitted_only eq 1 then begin frame_e_den_h=mean(timeedtemp_h(xfit)) frame_e1_h(file_number)=mean(timee1temp_h(xfit))/$ mean(timeedtemp_h(xfit)) frame_e2_h(file_number)=mean(timee2temp_h(xfit))/$ mean(timeedtemp_h(xfit)) endif else if fitted_only eq 2 then begin frame_e_den_h=mean(timeedtemp_h(notfit)) frame_e1_h(file_number)=mean(timee1temp_h(notfit))/$ mean(timeedtemp_h(notfit)) frame_e2_h(file_number)=mean(timee2temp_h(notfit))/$ mean(timeedtemp_h(notfit)) endif else begin frame_e_den_h=mean(timeedtemp_h(where(timee1temp_h ne 0))) frame_e1_h(file_number)=mean(timee1temp_h(where(timee1temp_h ne 0)))/$ mean(timeedtemp_h(where(timee1temp_h ne 0))) frame_e2_h(file_number)=mean(timee2temp_h(where(timee2temp_h ne 0)))/$ mean(timeedtemp_h(where(timee2temp_h ne 0))) endelse frame_e_h(file_number)=sqrt(frame_e1_h(file_number)^2+frame_e2_h(file_number)^2) ;Get the position angle associated with the elongation of the object if frame_e2_h(file_number) < 0 then begin frame_thetap_h(file_number)=$ -asin(sqrt(.5*(1-sqrt(1-frame_e2_h(file_number)^2/$ frame_e_h(file_number)^2)))) endif else begin frame_thetap_h(file_number)=$ asin(sqrt(.5*(1-sqrt(1-frame_e2_h(file_number)^2/$ frame_e_h(file_number)^2)))) endelse ;Recalculate e1 and e2 with new cosine angle frame_e1_h(file_number)=frame_e_h(file_number)*$ cos(2*frame_thetap_h(file_number)) frame_e2_h(file_number)=frame_e_h(file_number)*$ sin(2*frame_thetap_h(file_number)) ;low pass filter timee1num_l(file_number,*,*)=x_offset_hor_calc_l(file_number,*,*)*$ x_offset_hor_calc_l(file_number+depth,*,*)-$ y_offset_hor_calc_l(file_number,*,*)*$ y_offset_hor_calc_l(file_number+depth,*,*) timee2num_l(file_number,*,*)=(x_offset_hor_calc_l(file_number,*,*)*$ y_offset_hor_calc_l(file_number+depth,*,*)$ +x_offset_hor_calc_l(file_number+depth,*,*)$ *y_offset_hor_calc_l(file_number,*,*)) timeeden_l(file_number,*,*)=(x_offset_hor_calc_l(file_number,*,*)^2+$ y_offset_hor_calc_l(file_number,*,*)^2) timee1temp_l=timee1num_l(file_number,*,*) timee2temp_l=timee2num_l(file_number,*,*) timeedtemp_l=timeeden_l(file_number,*,*) ;IF K=0, the low pass filter will necessarily give e=1 if kval eq 0 then begin frame_e1_l(file_number)=0 frame_e2_l(file_number)=0 frame_e_l(file_number)=1 endif else if fitted_only eq 1 then begin frame_e_den_l=mean(timeedtemp_l(xfit)) frame_e1_l(file_number)=mean(timee1temp_l(xfit))/$ mean(timeedtemp_l(xfit)) frame_e2_l(file_number)=mean(timee2temp_l(xfit))/$ mean(timeedtemp_l(xfit)) endif else if fitted_only eq 2 then begin frame_e_den_l=mean(timeedtemp_l(notfit)) frame_e1_l(file_number)=mean(timee1temp_l(notfit))/$ mean(timeedtemp_l(notfit)) frame_e2_l(file_number)=mean(timee2temp_l(notfit))/$ mean(timeedtemp_l(notfit)) endif else begin frame_e_den_l=mean(timeedtemp_l(where(timee1temp_l ne 0))) frame_e1_l(file_number)=mean(timee1temp_l(where(timee1temp_l ne 0)))/$ mean(timeedtemp_l(where(timee1temp_l ne 0))) frame_e2_l(file_number)=mean(timee2temp_l(where(timee2temp_l ne 0)))/$ mean(timeedtemp_l(where(timee2temp_l ne 0))) endelse frame_e_l(file_number)=sqrt(frame_e1_l(file_number)^2+frame_e2_l(file_number)^2) ;Get the position angle associated with the elongation of the object if frame_e2_l(file_number) < 0 then begin frame_thetap_l(file_number)=$ -asin(sqrt(.5*(1-sqrt(1-frame_e2_l(file_number)^2/$ frame_e_l(file_number)^2)))) endif else begin frame_thetap_l(file_number)=$ asin(sqrt(.5*(1-sqrt(1-frame_e2_l(file_number)^2/$ frame_e_l(file_number)^2)))) endelse ;Recalculate e1 and e2 with new cosine angle frame_e1_l(file_number)=frame_e_l(file_number)*$ cos(2*frame_thetap_l(file_number)) frame_e2_l(file_number)=frame_e_l(file_number)*$ sin(2*frame_thetap_l(file_number)) endfor ;Take the average e values FOR EACH KVAL frame_e1_av_h(kval)=mean(frame_e1_h) frame_e2_av_h(kval)=mean(frame_e2_h) frame_e_av_h(kval)=mean(frame_e_h) frame_e1_av_l(kval)=mean(frame_e1_l) frame_e2_av_l(kval)=mean(frame_e2_l) frame_e_av_l(kval)=mean(frame_e_l) endfor if fitted_only eq 2 then begin ymax=1 endif else begin ymax=0.6 endelse if i eq 0 then begin erase emaxs=[max(frame_e_av_h),max(frame_e_av_l),max(frame_e_av_den_h),$ max(frame_e_av_den_l)] emax=max(emaxs) emins=[min(frame_e_av_h),min(frame_e_av_l),min(frame_e_av_den_h),$ min(frame_e_av_den_l)] emin=min(emins) plot,frame_e_av_h,psym=2,$ yrange=[0,ymax],color=black,$ background=white,position=[.15,.75,.95,.975],$ xtickname=[replicate(' ',13)],/normal,symsize=plotsymsize,$ ytickname=['0.0','0.2','0.4','0.6'],$ yticks=3 oplot,frame_e_av_l,psym=1,color=black,symsize=plotsymsize if png eq 1 then !p.charsize=4 xyouts,.9,.925,'(a)',color=black,/normal if png eq 1 then !p.charsize=6 endif else if i eq 1 then begin emaxs=[max(frame_e_av_h),max(frame_e_av_l),max(frame_e_av_den_h),$ max(frame_e_av_den_l)] emax=max(emaxs) emins=[min(frame_e_av_h),min(frame_e_av_l),min(frame_e_av_den_h),$ min(frame_e_av_den_l)] emin=min(emins) plot,frame_e_av_h,psym=2,$ ytitle='e',yrange=[0,ymax],color=black,$ background=white,position=[.15,.525,.95,.75],$ xtickname=[replicate(' ',13)],/normal,symsize=plotsymsize,$ ytickname=['0.0','0.2','0.4',' '],$ yticks=3 oplot,frame_e_av_l,psym=1,color=black,symsize=plotsymsize if png eq 1 then !p.charsize=4 xyouts,.9,.7,'(b)',color=black,/normal if png eq 1 then !p.charsize=6 endif else begin emaxs=[max(frame_e_av_h),max(frame_e_av_l),max(frame_e_av_den_h),$ max(frame_e_av_den_l)] emax=max(emaxs) emins=[min(frame_e_av_h),min(frame_e_av_l),min(frame_e_av_den_h),$ min(frame_e_av_den_l)] emin=min(emins) plot,frame_e_av_h,psym=2,$ yrange=[0,ymax],color=black,$ background=white,position=[.15,.3,.95,.525],$ xtickname=[replicate(' ',13)],/normal,symsize=plotsymsize,$ ytickname=['0.0','0.2','0.4',' '],$ yticks=3 oplot,frame_e_av_l,psym=1,color=black,symsize=plotsymsize if png eq 1 then !p.charsize=4 xyouts,.9,.475,'(c)',color=black,/normal if png eq 1 then !p.charsize=6 endelse get_lun,e_file openw,e_file,stats_directory+'e_filtered_values.txt' for j=0,kmax-1 do begin printf,e_file,FORMAT='(%"%f\t%f")',frame_e_av_h(j),frame_e_av_l(j) endfor close,e_file free_lun,e_file endfor ;======================================================================= ;*********************************************************************** ;NOW DO THE PHASE SCREEN lowstr=strtrim(string(low,format='(%"%10.1f")'),2) lcutstr=strtrim(string(lcut,format='(%"%10.1f")'),2) print,xwindspeed print,ywindspeed xwindspeed_m=xwindspeed*.17 ywindspeed_m=ywindspeed*.17 windspeeds=[xwindspeed,ywindspeed] max_windspeed=max(windspeeds) ;DIMENSIONS OF PHASE SCREEN AND SUB_SECTION OF SCREEN sub_sec=512 offset=sub_sec/2-12 grid_spots=long(625) num=25 num_of_cor=(phase_dim-sub_sec)/max_FrameDrift kmax=12 ;SIMULATIONS------------------------------------------------------------- ;ARRAY REFERENCE 25x25=625 FOR GRID subx=2*indgen(num)-25 ;The x-reference numbers suby=indgen(num)-12 ;The y-reference numbers ;Grid for partvelvec x_ref_col=12.5*subx+542 y_ref_col=25.0*suby+460 index=findgen(num*num) gridx=fltarr(num*num) gridy=fltarr(num*num) gridx(index)=x_ref_col(index mod num) gridy(index)=y_ref_col(index/num) ;FILE IO atmosphere_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/lsst_sim/atmosphere/' image_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/ShackHartmann_Cal/' screen_directory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/lsst_sim/atmosphere/' phase_files=file_search(screen_directory+'phase_'+strtrim(lowstr,2)+'_'+$ strtrim(lcutstr,2)+'*.out',$ count=num_phase_files) total_num_of_cor=num_of_cor*num_phase_files ;MAKE SUB ARRAY FOR SLOPES ACROSS SUBAPERTURE x_slopes_sub_arr=fltarr(total_num_of_cor,25,25) y_slopes_sub_arr=fltarr(total_num_of_cor,25,25) magnitude=fltarr(total_num_of_cor,25,25) phi=fltarr(total_num_of_cor,25,25) phi_array=fltarr(total_num_of_cor) r_array=fltarr(total_num_of_cor) ;FFTs xfft_array=complexarr(total_num_of_cor,num,num) yfft_array=complexarr(total_num_of_cor,num,num) ;Angles frame_theta=fltarr(total_num_of_cor) ;************************************************************************ for screen_num=0, num_phase_files-1 do begin ;Open the phase screen binary file and read floating point (data_type=4) openr,phase,phase_files(screen_num),/get_lun screen=read_binary(phase,data_dims=[phase_dim,phase_dim],data_type=4) close,phase free_lun,phase ;Use Central Differencing to get slopes and normalize to 1 xslopes=fltarr(phase_dim,phase_dim) yslopes=fltarr(phase_dim,phase_dim) xslopes(1:1022,1:1022)=(screen(2:1023,1:1022)-screen(0:1021,1:1022))/2 yslopes(1:1022,1:1022)=(screen(1:1022,2:1023)-screen(1:1022,0:1021))/2 yslopes=yslopes/sqrt(mean(yslopes^2)) xslopes=xslopes/sqrt(mean(xslopes^2)) ;WINDSPEED TO MOVE PHASE SCREEN ACROSS SUBAPERTURE ;FILL THE ARRAYS FOR LATER PROCESSING for increment=0,num_of_cor-1 do begin tot_increment=num_of_cor*screen_num+increment If MultiWindSpeeds eq 1 then begin xWindSpeed = xwinds(floor($ float(tot_increment)/(float(total_num_of_cor)/N_elements(xwinds)))) yWindSpeed = ywinds(floor($ float(tot_increment)/(float(total_num_of_cor)/N_elements(ywinds)))) xwindsign=xWindSpeed/abs(xWindSpeed) ywindsign=yWindSpeed/abs(yWindSpeed) endif print,xwindspeed,ywindspeed exp_time=0.030 ; milliseconds numXcells=ceil(abs(xWindSpeed)*exp_time/0.17)+1; numYcells=ceil(abs(yWindSpeed)*exp_time/0.17)+1; ;CORRELATIONS x_segment=increment*xFrameDrift y_segment=increment*yFrameDrift ;SLOPES AVERAGED WITH THE WIND Cells=[numXCells,numYCells] ;Make Array of numcells numCells=max(Cells) ;Total num will be greater of 2 xWindSlopes=fltarr(numCells,num,num) ;Number of x slopes to average yWindSlopes=fltarr(numCells,num,num) ;Number of y slopes to average iniXoffset=phase_dim-offset-x_segment iniYoffset=phase_dim-offset-y_segment for WindInc=0, numCells-1 do begin XWindoff=xWindSign*$ floor((float(numXCells)/float(numCells))*float(WindInc)) yWindoff=yWindSign*$ floor((float(numYCells)/float(numCells))*float(WindInc)) xWindSlopes(WindInc,*,*)=xslopes[$ iniXoffset-xWindOff:iniXoffset-xWindOff+num-1,$ iniYoffset-yWindOff:iniYoffset-yWindOff+num-1] yWindSlopes(WindInc,*,*)=yslopes[$ iniXoffset-xWindOff:iniXoffset-xWindOff+num-1,$ iniYoffset-yWindOff:iniYoffset-yWindOff+num-1] endfor x_slopes_sub_arr(tot_increment,*,*)=avg(xWindSlopes,0) y_slopes_sub_arr(tot_increment,*,*)=avg(yWindSlopes,0) xtemp=fltarr(num*num) ytemp=fltarr(num*num) xtemp=x_slopes_sub_arr(tot_increment,*,*) ytemp=y_slopes_sub_arr(tot_increment,*,*) xtemp(nonfiducial)=0 ytemp(nonfiducial)=0 if avg_sub eq 1 then begin xfit=where(xtemp ne 0) yfit=where(ytemp ne 0) mean_xoff=mean(xtemp(xfit)) mean_yoff=mean(ytemp(yfit)) xtemp(xfit)=xtemp(xfit)-mean_xoff ytemp(yfit)=ytemp(yfit)-mean_yoff phi_array(tot_increment)=atan(mean_yoff,mean_xoff) r_array(tot_increment)=sqrt(mean_xoff^2+mean_yoff^2) endif else if avg_sub eq 2 then begin x_fft=complexarr(num,num) y_fft=complexarr(num,num) x_fft(*,*)=fft(xtemp) y_fft(*,*)=fft(ytemp) x_fft_re=complexarr(num,num) y_fft_re=complexarr(num,num) for n=0,num-1 do begin for m=0,num-1 do begin x_fft_re(subn(n),subm(m))=x_fft(n,m) y_fft_re(subn(n),subm(m))=y_fft(n,m) endfor endfor kx0_array(tot_increment)=x_fft_re(12,12) ky0_array(tot_increment)=y_fft_re(12,12) phi_array(tot_increment)=atan(y_fft_re(12,12),x_fft_re(12,12)) r_array(tot_increment)=sqrt(x_fft_re(12,12)^2+y_fft_re(12,12)^2) x_fft_re(12,12)=0 y_fft_re(12,12)=0 for n=0,num-1 do begin for m=0,num-1 do begin x_fft(n,m)=x_fft_re(subn(n),subm(m)) y_fft(n,m)=y_fft_re(subn(n),subm(m)) endfor endfor xtemp=fft(x_fft,/inverse) ytemp=fft(y_fft,/inverse) endif else if avg_sub eq 3 then begin xfit=where(xtemp ne 0) yfit=where(ytemp ne 0) xtemp(xfit)=xtemp(xfit)+xscale*xoff(xfit) ytemp(yfit)=ytemp(yfit)+yscale*yoff(yfit) endif x_slopes_sub_arr(tot_increment,*,*)=xtemp y_slopes_sub_arr(tot_increment,*,*)=ytemp ;GET FFT OF SLOPES xfft_array(tot_increment,*,*)=fft(x_slopes_sub_arr(tot_increment,*,*)) yfft_array(tot_increment,*,*)=fft(y_slopes_sub_arr(tot_increment,*,*)) endfor endfor ;Correlations total_num_of_cor=num_of_cor*num_phase_files timee1=fltarr(total_num_of_cor-depth,num,num) timee2=fltarr(total_num_of_cor-depth,num,num) timee1num=fltarr(total_num_of_cor-depth,num,num) timee2num=fltarr(total_num_of_cor-depth,num,num) timeeden=fltarr(total_num_of_cor-depth,num,num) timee1_arr=fltarr(num,num) timee2_arr=fltarr(num,num) frame_e1=fltarr(total_num_of_cor) frame_e2=fltarr(total_num_of_cor) frame_e=fltarr(total_num_of_cor) frame_thetap=fltarr(total_num_of_cor) frame_thetam=fltarr(total_num_of_cor) ;Average Ellipticity frame_e1_av_h=fltarr(kmax) frame_e2_av_h=fltarr(kmax) frame_e_av_h=fltarr(kmax) frame_e_av_den_h=fltarr(kmax) frame_e1_av_l=fltarr(kmax) frame_e2_av_l=fltarr(kmax) frame_e_av_l=fltarr(kmax) frame_e_av_den_l=fltarr(kmax) ;Indices and such no_bins=8 a=findgen(13)+12 b=findgen(12) subn=[a,b] subm=subn readcol,offset_files(0),xcoord,ycoord,format='x,f,x,f,x' for kval=0,kmax-1 do begin ;Correlations timee1=fltarr(total_num_of_cor-depth,num,num) timee2=fltarr(total_num_of_cor-depth,num,num) timee1num_h=fltarr(total_num_of_cor-depth,num,num) timee2num_h=fltarr(total_num_of_cor-depth,num,num) timeeden_h=fltarr(total_num_of_cor-depth,num,num) timee1num_l=fltarr(total_num_of_cor-depth,num,num) timee2num_l=fltarr(total_num_of_cor-depth,num,num) timeeden_l=fltarr(total_num_of_cor-depth,num,num) timee1_arr=fltarr(num,num) timee2_arr=fltarr(num,num) ;for high pass filter frame_e1_h=fltarr(total_num_of_cor) frame_e2_h=fltarr(total_num_of_cor) frame_e_h=fltarr(total_num_of_cor) frame_thetap_h=fltarr(total_num_of_cor) frame_thetam_h=fltarr(total_num_of_cor) frame_e_den_h=fltarr(total_num_of_cor) ;for low pass filter frame_e1_l=fltarr(total_num_of_cor) frame_e2_l=fltarr(total_num_of_cor) frame_e_l=fltarr(total_num_of_cor) frame_thetap_l=fltarr(total_num_of_cor) frame_thetam_l=fltarr(total_num_of_cor) frame_e_den_l=fltarr(total_num_of_cor) x_offset_hor_calc_h=complexarr(total_num_of_cor,num,num) y_offset_hor_calc_h=complexarr(total_num_of_cor,num,num) x_offset_hor_calc_l=complexarr(total_num_of_cor,num,num) y_offset_hor_calc_l=complexarr(total_num_of_cor,num,num) for file_number=0,total_num_of_cor-depth-1 do begin ;Average the fft x_fft=complexarr(num,num) y_fft=complexarr(num,num) x_low_fft=complexarr(num,num) x_low_fft_re=complexarr(num,num) y_low_fft=complexarr(num,num) y_low_fft_re=complexarr(num,num) y_high_fft_re=complexarr(num,num) x_high_fft=complexarr(num,num) y_high_fft=complexarr(num,num) x_high_fft_re=complexarr(num,num) x_fft(*,*)=xfft_array(file_number,*,*) y_fft(*,*)=yfft_array(file_number,*,*) x_fft_re=complexarr(num,num) y_fft_re=complexarr(num,num) ;Rearrange indices for n=0,num-1 do begin for m=0,num-1 do begin x_fft_re(subn(n),subm(m))=x_fft(n,m) y_fft_re(subn(n),subm(m))=y_fft(n,m) endfor endfor ;Filter kx=0, ky=0 x_high_fft_re=x_fft_re y_high_fft_re=y_fft_re x_high_fft_re((12-kval):(12+kval),(12-kval):(12+kval))=0 y_high_fft_re((12-kval):(12+kval),(12-kval):(12+kval))=0 x_low_fft_re((12-kval):(12+kval),(12-kval):(12+kval))=$ x_fft_re((12-kval):(12+kval),(12-kval):(12+kval)) y_low_fft_re((12-kval):(12+kval),(12-kval):(12+kval))=$ y_fft_re((12-kval):(12+kval),(12-kval):(12+kval)) ;Rearrange for n=0,num-1 do begin for m=0,num-1 do begin x_low_fft(n,m)=x_low_fft_re(subn(n),subm(m)) y_low_fft(n,m)=y_low_fft_re(subn(n),subm(m)) x_high_fft(n,m)=x_high_fft_re(subn(n),subm(m)) y_high_fft(n,m)=y_high_fft_re(subn(n),subm(m)) endfor endfor x_offset_hor_calc_h(file_number,*,*)=fft(x_high_fft,/inverse) y_offset_hor_calc_h(file_number,*,*)=fft(y_high_fft,/inverse) x_offset_hor_calc_l(file_number,*,*)=fft(x_low_fft,/inverse) y_offset_hor_calc_l(file_number,*,*)=fft(y_low_fft,/inverse) test=0 if test eq 1 then begin set_plot,'x' xoff=fltarr(625) yoff=fltarr(625) xoff(*)=x_offset_hor_calc_h(file_number,*,*) yoff(*)=y_offset_hor_calc_h(file_number,*,*) partvelvec,xoff,yoff,$ gridx,gridy,yrange=[150,800],$ xrange=[200,850],position=[0.1,0.1,0.9,0.9],$ title='High Passed Average Offsets for all images for k='+$ strtrim(kval,2) xyouts,gridx(624)-100,gridy(624)+125,'length='+$ strtrim(xoff(624),2) wait,4 endif ;high pass filter timee1num_h(file_number,*,*)=x_offset_hor_calc_h(file_number,*,*)*$ x_offset_hor_calc_h(file_number+depth,*,*)-$ y_offset_hor_calc_h(file_number,*,*)*$ y_offset_hor_calc_h(file_number+depth,*,*) timee2num_h(file_number,*,*)=(x_offset_hor_calc_h(file_number,*,*)*$ y_offset_hor_calc_h(file_number+depth,*,*)$ +x_offset_hor_calc_h(file_number+depth,*,*)$ *y_offset_hor_calc_h(file_number,*,*)) timeeden_h(file_number,*,*)=(x_offset_hor_calc_h(file_number,*,*)^2+$ y_offset_hor_calc_h(file_number,*,*)^2) timee1temp_h=timee1num_h(file_number,*,*) timee2temp_h=timee2num_h(file_number,*,*) timeedtemp_h=timeeden_h(file_number,*,*) if fitted_only eq 1 then begin frame_e_den_h=mean(timeedtemp_h(xfit)) frame_e1_h(file_number)=mean(timee1temp_h(fitted))/$ mean(timeedtemp_h(fitted)) frame_e2_h(file_number)=mean(timee2temp_h(fitted))/$ mean(timeedtemp_h(fitted)) endif else if fitted_only eq 2 then begin frame_e_den_h=mean(timeedtemp_h(nonfiducial)) frame_e1_h(file_number)=mean(timee1temp_h(nonfiducial))/$ mean(timeedtemp_h(nonfiducial)) frame_e2_h(file_number)=mean(timee2temp_h(nonfiducial))/$ mean(timeedtemp_h(nonfiducial)) endif else begin frame_e_den_h=mean(timeedtemp_h(where(timee1temp_h ne 0))) frame_e1_h(file_number)=mean(timee1temp_h(where(timee1temp_h ne 0)))/$ mean(timeedtemp_h(where(timee1temp_h ne 0))) frame_e2_h(file_number)=mean(timee2temp_h(where(timee2temp_h ne 0)))/$ mean(timeedtemp_h(where(timee2temp_h ne 0))) endelse frame_e_h(file_number)=sqrt(frame_e1_h(file_number)^2+frame_e2_h(file_number)^2) signe1=abs(frame_e1_h(file_number))/frame_e1_h(file_number) signe2=abs(frame_e2_h(file_number))/frame_e2_h(file_number) ;Get the position angle associated with the elongation of the object ;Should be between pi/4 and pi/2 if e1 is negative and between 0 and pi/2 if ;e1 is positive. Should be in second quadrant if e2 is negative frame_thetap_h(file_number)=$ signe2*asin(sqrt(.5*(1-signe1*sqrt(1-frame_e2_h(file_number)^2/$ frame_e_h(file_number)^2)))) ;Angular Distribution we are interested in if kval eq 0 then begin frame_theta(file_number)=frame_thetap_h(file_number) endif ;Recalculate e1 and e2 with new cosine angle frame_e1_h(file_number)=frame_e_h(file_number)*$ cos(2*frame_thetap_h(file_number)) frame_e2_h(file_number)=frame_e_h(file_number)*$ sin(2*frame_thetap_h(file_number)) ;-----------------------------------------------------------LOW PASS BEGIN ;low pass filter timee1num_l(file_number,*,*)=x_offset_hor_calc_l(file_number,*,*)*$ x_offset_hor_calc_l(file_number+depth,*,*)-$ y_offset_hor_calc_l(file_number,*,*)*$ y_offset_hor_calc_l(file_number+depth,*,*) timee2num_l(file_number,*,*)=(x_offset_hor_calc_l(file_number,*,*)*$ y_offset_hor_calc_l(file_number+depth,*,*)$ +x_offset_hor_calc_l(file_number+depth,*,*)$ *y_offset_hor_calc_l(file_number,*,*)) timeeden_l(file_number,*,*)=(x_offset_hor_calc_l(file_number,*,*)^2+$ y_offset_hor_calc_l(file_number,*,*)^2) timee1temp_l=timee1num_l(file_number,*,*) timee2temp_l=timee2num_l(file_number,*,*) timeedtemp_l=timeeden_l(file_number,*,*) ;IF K=0, the low pass filter will necessarily give e=1 if kval eq 0 then begin frame_e1_l(file_number)=0 frame_e2_l(file_number)=0 frame_e_l(file_number)=1 endif else if fitted_only eq 1 then begin frame_e_den_l=mean(timeedtemp_l(fitted)) frame_e1_l(file_number)=mean(timee1temp_l(fitted))/$ mean(timeedtemp_l(fitted)) frame_e2_l(file_number)=mean(timee2temp_l(fitted))/$ mean(timeedtemp_l(fitted)) endif else if fitted_only eq 2 then begin frame_e_den_l=mean(timeedtemp_l(nonfiducial)) frame_e1_l(file_number)=mean(timee1temp_l(nonfiducial))/$ mean(timeedtemp_l(nonfiducial)) frame_e2_l(file_number)=mean(timee2temp_l(nonfiducial))/$ mean(timeedtemp_l(nonfiducial)) endif else begin frame_e_den_l=mean(timeedtemp_l(where(timee1temp_l ne 0))) frame_e1_l(file_number)=mean(timee1temp_l(where(timee1temp_l ne 0)))/$ mean(timeedtemp_l(where(timee1temp_l ne 0))) frame_e2_l(file_number)=mean(timee2temp_l(where(timee2temp_l ne 0)))/$ mean(timeedtemp_l(where(timee2temp_l ne 0))) endelse frame_e_l(file_number)=sqrt(frame_e1_l(file_number)^2+frame_e2_l(file_number)^2) signe1=abs(frame_e1_l(file_number))/frame_e1_l(file_number) signe2=abs(frame_e2_l(file_number))/frame_e2_l(file_number) ;Get the position angle associated with the elongation of the object ;Should be between pi/4 and pi/2 if e1 is negative and between 0 and pi/2 if ;e1 is positive. Should be in second quadrant if e2 is negative frame_thetap_l(file_number)=$ signe2*asin(sqrt(.5*(1+signe1*sqrt(1-frame_e2_l(file_number)^2/$ frame_e_l(file_number)^2)))) ;Recalculate e1 and e2 with new cosine angle frame_e1_l(file_number)=frame_e_l(file_number)*$ cos(2*frame_thetap_l(file_number)) frame_e2_l(file_number)=frame_e_l(file_number)*$ sin(2*frame_thetap_l(file_number)) endfor ;Take the average e values FOR EACH KVAL frame_e1_av_h(kval)=mean(frame_e1_h) frame_e2_av_h(kval)=mean(frame_e2_h) frame_e_av_h(kval)=mean(frame_e_h) frame_e_av_den_h(kval)=mean(timeedtemp_h(where(timee1temp_h ne 0))) if kval eq 0 then begin frame_e1_av_l(kval)=1 frame_e2_av_l(kval)=1 frame_e_av_l(kval)=1 frame_e_av_den_l(kval)=1 endif else begin frame_e1_av_l(kval)=mean(frame_e1_l) frame_e2_av_l(kval)=mean(frame_e2_l) frame_e_av_l(kval)=mean(frame_e_l) frame_e_av_den_l(kval)=mean(timeedtemp_l(where(timee1temp_l ne 0))) endelse endfor emaxs=[max(frame_e_av_h),max(frame_e_av_l),max(frame_e_av_den_h),$ max(frame_e_av_den_l)] emax=max(emaxs) emins=[min(frame_e_av_h),min(frame_e_av_l),min(frame_e_av_den_h),$ min(frame_e_av_den_l)] emin=min(emins) plot,frame_e_av_h,psym=2,$ yrange=[0,ymax],color=black,$ background=white,position=[.15,.075,.95,.3],/normal,$ xtitle='f',symsize=plotsymsize,$ ytickname=['0.0','0.2','0.4',' '],$ yticks=3 oplot,frame_e_av_l,psym=1,color=black,symsize=plotsymsize if png eq 1 then !p.charsize=4 xyouts,.9,.25,'(d)',color=black,/normal if png eq 1 then !p.charsize=6 get_lun,e_file openw,e_file,atmosphere_dir+'e_filtered_values_'+label_string+$ '.txt' for i=0,kmax-1 do begin printf,e_file,FORMAT='(%"%f\t%f")',frame_e_av_h(i),frame_e_av_l(i) endfor close,e_file free_lun,e_file if png eq 1 then begin jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,postscript_dir+'e_filtered_all_nights_plus_phase_screen'+$ label_string+'.png',$ jpgimg,reds,greens,blues device,/close endif else begin device,/close endelse set_plot,'x' stop end