;PHASE_Screen_ellipticity_filter ; ;PURPOSE: ; This script will take in a phase screen generated by turb2d.c (written ; by Garret Jernigan) and obtain the Fourier Spectrum for a frame ; identical to the ones obtained from the SOAR CWFS wavefront sensor. ; ; With that Fourier Spectrum, it will do high and low pass filtering ; and calculate the Ellipticity for each frame and then average them. ;CALLING SEQUENCE: ; phase_screen_ellipticity_filter,[date,[FrameDrift,[depth]]] ; ;KEYWORDS: ; Date: The date (10,11, or 12) of the directory to which the stats will ; be written ; Windspeed: The speed at which the frozen screen will be dragged across ; the screen. The units pro final_phase_screen_ellip_theta_histogram,$ xFrameDrift=xFrameDrift,yFrameDrift=yFrameDrift,$ xWindSpeed=xWindSpeed,yWindSpeed=yWindSpeed,$ depth=depth,phase_dim=phase_dim,$ low=low,lcut=lcut,$ MultiWindSpeeds=MultiWindSpeeds 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=17. ;m/sec If (not keyword_set(yWindSpeed)) then yWindSpeed=4. ;m/sec If (not keyword_set(depth)) then depth=0 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=2 If (not keyword_set(x_off)) then x_off=0 If (not keyword_set(y_off)) then y_off=0 ;OPTIONS: ********************************************************** ;Avg_sub=0 means don't subtract any offset ;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 ;Covariance=0 don't do covariance analysis ;Covariance=1 do covariance analaysis Covariance=1 ;MultiWindSpeeds=1 Use multiple wind speeds ;MultiWindSpeeds=2 Use only one wind speed ;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 else if MultiWindSpeeds eq 2 then begin xwinds=[0] ywinds=[0] 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 =label_string+'NoFiducial_' endif else begin label_string =label_string+'Fiducial_' endelse xwind=strtrim(string(xWindSpeed,format='(%"%10.1f")'),2) ywind=strtrim(string(yWindSpeed,format='(%"%10.1f")'),2) lowstr=strtrim(string(low,format='(%"%10.1f")'),2) lcutstr=strtrim(string(lcut,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) 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 ;WINDSPEEDS 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(abs(xWindSpeed)*exp_time/0.17)+1; numYcells=ceil(abs(yWindSpeed)*exp_time/0.17)+1; ;DIMENSIONS OF PHASE SCREEN AND SUB_SECTION OF SCREEN sub_sec=512 offset=sub_sec/2-12 grid_spots=long(625) index=findgen(625) num=25 num_of_cor=(phase_dim-sub_sec)/max_FrameDrift kmax=12 ;ARRAY REFERENCE 25x25=625 FOR GRID subx=2*indgen(num)-25 ;The x-reference numbers suby=indgen(num)-12 ;The y-reference numbers ;FILE IO atmosphere_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/lsst_sim/atmosphere/' Postscript_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/phase_screen_figs/' screen_directory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/lsst_sim/atmosphere/' ;Indices and such no_bins=8 a=findgen(13)+12 b=findgen(12) subn=[a,b] subm=subn ;POSSIBLY EXIT IF THE CELLS ARE BAD if numXcells gt offset or numYcells gt offset then begin print,"Wind Speed will bring frame out of bounds" stop endif if avg_sub eq 1 then begin evals_file_name=screen_directory+'ellipticity_avg_sub_low_'+$ label_string+'.txt' endif else if avg_sub eq 2 then begin evals_file_name=screen_directory+'ellipticity_kx0_ky0_sub_low_'+$ label_string+'.txt' endif else if avg_sub eq 3 then begin evals_file_name=screen_directory+'ellipticity_low_'+lowinmeters+$ '_high_'+highinmeters+'.txt' endif else begin evals_file_name=screen_directory+'ellipticity_low_'+lowinmeters+$ '_high_'+highinmeters+'.txt' endelse ;PRINT OUT THE ELLIPTICITIES FOR EACH FRAME TO A FILE get_lun, evals openw,evals,evals_file_name 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 xfit=fitted notfit=nonfiducial 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) x_slopes_sub_arr_hor=fltarr(total_num_of_cor,num*num) y_slopes_sub_arr_hor=fltarr(total_num_of_cor,num*num) x_slopes_sub_arr_ver=fltarr(total_num_of_cor,num*num) y_slopes_sub_arr_ver=fltarr(total_num_of_cor,num*num) magnitude=fltarr(total_num_of_cor,25,25) phi=fltarr(total_num_of_cor,25,25) kx0_array=fltarr(total_num_of_cor) ky0_array=fltarr(total_num_of_cor) 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) ;PREPARE TO GET THETA AND E AND PRINT OUT FOR EACH FRAME TO A FILE get_lun, evals openw,evals,evals_file_name 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)) ;DRIFT FRAME 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 else if MultiWindSpeeds eq 2 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,*,*) if UseFid eq 1 then begin xtemp(nonfiducial)=0 ytemp(nonfiducial)=0 endif 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 ;Horizontally arranged x_slopes_sub_arr(tot_increment,*,*)=xtemp y_slopes_sub_arr(tot_increment,*,*)=ytemp ;Veritcally arranged for Covariance x_slopes_sub_arr_hor(tot_increment,*)=xtemp y_slopes_sub_arr_hor(tot_increment,*)=ytemp x_slopes_sub_arr_ver(tot_increment,*)=xtemp((index mod 25)*25+index/25) y_slopes_sub_arr_ver(tot_increment,*)=ytemp((index mod 25)*25+index/25) ;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 ;******************************************************* ;DATA IS IN; NOW PROCESS 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) ;for high pass filter frame_e1=fltarr(total_num_of_cor) frame_e2=fltarr(total_num_of_cor) frame_e=fltarr(total_num_of_cor) frame_theta=fltarr(total_num_of_cor) frame_e_den=fltarr(total_num_of_cor) ;COVARIANCE Cxt_full=fltarr(total_num_of_cor,num,depth+1) Cxl_full=fltarr(total_num_of_cor,num,depth+1) Cyt_full=fltarr(total_num_of_cor,num,depth+1) Cyl_full=fltarr(total_num_of_cor,num,depth+1) for file_number=0,total_num_of_cor-depth-1 do begin timee1num(file_number,*,*)=x_slopes_sub_arr(file_number,*,*)*$ x_slopes_sub_arr(file_number+depth,*,*)-$ y_slopes_sub_arr(file_number,*,*)*$ y_slopes_sub_arr(file_number+depth,*,*) timee2num(file_number,*,*)=(x_slopes_sub_arr(file_number,*,*)*$ y_slopes_sub_arr(file_number+depth,*,*)$ +x_slopes_sub_arr(file_number+depth,*,*)$ *y_slopes_sub_arr(file_number,*,*)) timeeden(file_number,*,*)=(x_slopes_sub_arr(file_number,*,*)^2+$ y_slopes_sub_arr(file_number,*,*)^2) timee1temp=timee1num(file_number,*,*) timee2temp=timee2num(file_number,*,*) timeedtemp=timeeden(file_number,*,*) if fitted_only eq 1 then begin frame_e_den=mean(timeedtemp(xfit)) frame_e1(file_number)=mean(timee1temp(xfit))/$ mean(timeedtemp(xfit)) frame_e2(file_number)=mean(timee2temp(xfit))/$ mean(timeedtemp(xfit)) endif else if fitted_only eq 2 then begin frame_e_den=mean(timeedtemp(notfit)) frame_e1(file_number)=mean(timee1temp(notfit))/$ mean(timeedtemp(notfit)) frame_e2(file_number)=mean(timee2temp(notfit))/$ mean(timeedtemp(notfit)) endif else begin frame_e_den=mean(timeedtemp(where(timee1temp ne 0))) frame_e1(file_number)=mean(timee1temp(where(timee1temp ne 0)))/$ mean(timeedtemp(where(timee1temp ne 0))) frame_e2(file_number)=mean(timee2temp(where(timee2temp ne 0)))/$ mean(timeedtemp(where(timee2temp ne 0))) endelse frame_e(file_number)=sqrt(frame_e1(file_number)^2+frame_e2(file_number)^2) signe1=abs(frame_e1(file_number))/frame_e1(file_number) signe2=abs(frame_e2(file_number))/frame_e2(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_theta(file_number)=$ signe2*asin(sqrt(.5*(1-signe1*sqrt(1-frame_e2(file_number)^2/$ frame_e(file_number)^2)))) ;Recalculate e1 and e2 with new cosine angle frame_e1(file_number)=frame_e(file_number)*$ cos(2*frame_theta(file_number)) frame_e2(file_number)=frame_e(file_number)*$ sin(2*frame_theta(file_number)) printf,evals,FORMAT='(%"%f\t%f")',$ frame_theta(file_number),frame_e(file_number) ;*************************************************COVARIANCE if covariance eq 1 then begin Cxl=fltarr(num,depth+1) Cxt=fltarr(num,depth+1) Cyl=fltarr(num,depth+1) Cyt=fltarr(num,depth+1) for time=0, depth do begin for spacing=0, num-1 do begin Cxl_arr=fltarr(num*num) Cxt_arr=fltarr(num*num) Cyl_arr=fltarr(num*num) Cyt_arr=fltarr(num*num) ;SCAN FIRST ALONG THE ROWS for column_index=0, num-1 do begin for row_index=0, num-1 do begin spot_index=column_index*num+row_index ;If we have reached the edge, start on the next row If (row_index+spacing) ge num then break ;Check to see if both have valid values If (x_slopes_sub_arr_hor(file_number,spot_index) ne -10) and $ (x_slopes_sub_arr_hor(file_number+time,spot_index+spacing) ne -10) $ then begin ;add offset in to tweek the plots Cxl_arr(spot_index)=$ (x_slopes_sub_arr_hor(file_number,spot_index)+x_off)*$ (x_slopes_sub_arr_hor(file_number+time,spot_index+spacing)+x_off) Cxt_arr(spot_index)=$ (y_slopes_sub_arr_hor(file_number,spot_index)+y_off)*$ (y_slopes_sub_arr_hor(file_number+time,spot_index+spacing)+y_off) endif endfor endfor correlated=where(Cxl_arr ne 0) If correlated(0) ne -1 then begin Cxl(spacing,time)=mean(Cxl_arr(correlated)) Cxt(spacing,time)=mean(Cxt_arr(correlated)) endif else begin print,'no correlations',file_number Cxl(spacing,time)=-10 Cxt(spacing,time)=-10 endelse ;SCAN SECOND ALONG Y for column_index=0, num-1 do begin for row_index=0, num -1 do begin spot_index=column_index*num+row_index ;If we have reached the edge, start on the next row If (row_index+spacing) ge num then break ;Check to see if both have valid values If (y_slopes_sub_arr_ver(file_number,spot_index) ne -10) and $ (y_slopes_sub_arr_ver(file_number+time,spot_index+spacing) ne -10) $ then begin Cyl_arr(spot_index)=$ (y_slopes_sub_arr_ver(file_number,spot_index)+y_off)*$ (y_slopes_sub_arr_ver(file_number+time,spot_index+spacing)+y_off) Cyt_arr(spot_index)=$ (x_slopes_sub_arr_ver(file_number,spot_index)+x_off)*$ (x_slopes_sub_arr_ver(file_number+time,spot_index+spacing)+x_off) endif endfor endfor correlated=where(Cyl_arr ne 0) If correlated(0) ne -1 then begin Cyl(spacing,time)=mean(Cyl_arr(correlated)) Cyt(spacing,time)=mean(Cyt_arr(correlated)) endif else begin print,'no correlations',file_number Cyl(spacing,time)=-10 Cyt(spacing,time)=-10 endelse endfor endfor Cxt_full(file_number,*,*)=Cxt Cxl_full(file_number,*,*)=Cxl Cyt_full(file_number,*,*)=Cyt Cyl_full(file_number,*,*)=Cyl endif ;**********************************END COVARIANCE endfor ;*********************************************************** close, evals free_lun, evals ;*******************************************CALCULATE Bn if covariance eq 1 then begin Cxt_means=fltarr(num,depth+1) Cxl_means=fltarr(num,depth+1) Cyt_means=fltarr(num,depth+1) Cyl_means=fltarr(num,depth+1) for m=0, depth do begin get_lun,correlations openw,correlations,atmosphere_dir+'PhaseScreen'+label_string+$ 'Depth'+string(strtrim(m,2))+'Correlations.txt' printf,correlations,$ FORMAT='(%"%s\t%s\t%s\t%s")','Long. X Cor.','Trans. X Cor.','Long. Y Cor.','Trans. Y Cor' for n=0, num-1 do begin cxt_temp=cxt_full(*,n,m) inc_cxt=where((cxt_temp ne 0) and (cxt_temp ne -10)) Cxt_means(n,m)=mean(Cxt_temp(inc_cxt)) cxl_temp=cxl_full(*,n,m) inc_cxl=where((cxl_temp ne 0) and (cxl_temp ne -10)) Cxl_means(n,m)=mean(Cxl_temp(inc_cxl)) cyt_temp=cyt_full(*,n,m) inc_cyt=where((cyt_temp ne 0) and (cyt_temp ne -10)) Cyt_means(n,m)=mean(Cyt_temp(inc_cyt)) cyl_temp=cyl_full(*,n,m) inc_cyl=where((cyl_temp ne 0) and (cyl_temp ne -10)) Cyl_means(n,m)=mean(Cyl_temp(inc_cyl)) printf,correlations,FORMAT='(%"%f\t%f\t%f\t%f")',$ Cxl_means(n,m),Cxt_means(n,m),Cyl_means(n,m),Cyt_means(n,m) endfor close,correlations free_lun,correlations endfor endif ;*******************************************Begin HISTOGRAMS thetamin=-!pi/2 ;THETA no_bins=20. thetabinwidth=!pi/no_bins thetamax=thetamin+(no_bins-1)*thetabinwidth thetahist=histogram(frame_theta,Nbins=no_bins,locations=loc_theta,$ reverse_indices=theta_ind,$ max=thetamax,min=thetamin) thetabinsize=(thetamax-thetamin)/(no_bins-1) emin=0 ;ELLIPTICITY noe_bins=20. ebinwidth=(1-emin)/noe_bins emax=emin+(noe_bins)*ebinwidth ehist=histogram(frame_e,Nbins=noe_bins,locations=loc_e,$ reverse_indices=e_ind,$ max=emax,min=emin) ebinsize=(emax-emin)/(noe_bins-1) phimin=-!pi ;CENTROID OFFSET ANGLE no_bins=20. phibinwidth=2*!pi/no_bins phimax=phimin+(no_bins-1)*phibinwidth phihist=histogram(phi_array,Nbins=no_bins,locations=loc_phi,$ reverse_indices=phi_ind,$ max=phimax,min=phimin) phibinsize=(phimax-phimin)/(no_bins-1) ;**************************************************************** ;PLOT set_plot,'z' erase device, set_font='Courier' device,set_resolution=[800,600] !p.charsize=1 !p.charthick=1.2 !x.thick=2 !y.thick=2 !p.thick=4 !p.noerase=0 loadct,39 white='FFFFFF'x black='000000'x blue='FF0000'x red='0000FF'x !p.multi=[0,1,1] ;THETA_HISTOGRAM***************************************************** plot,loc_theta+thetabinsize/2,thetahist,psym=10,xtitle='Theta',$ ytitle='Number of Events',xrange=[-1.6,1.6],$ background=white,color=black jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,postscript_dir+'Theta_Histogram'+label_string+'.png',$ jpgimg,reds,greens,blues erase ;E_HISTOGRAM******************************************************** plot,loc_e,ehist,psym=10,xtitle='e',$ ytitle='Number of Events',$ background=white,color=black jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,postscript_dir+'E_Histogram_'+label_string+'.png',$ jpgimg,reds,greens,blues erase ;PHI_HISTOGRAM****************************************************** plot,loc_phi+phibinsize/2,phihist,psym=10,xtitle='Theta',$ ytitle='Number of Events',$ title='Histogram of Centroid Angle for low_'$ +lowinmeters+'_high_'+highinmeters,$ background=white,color=black jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,postscript_dir+'Phi_Histogram'+label_string+'.png',$ jpgimg,reds,greens,blues erase ;PHI_THETA_SCATTER************************************************* plot,phi_array,frame_theta,psym=2,$ xtitle='Azimuthal angle of centroid',$ ytitle='Ellipticity angle that maximizes e1 (Theta)',$ title='Scatter Plot of Theta vs. Phi for low_'$ +lowinmeters+'_high_'+highinmeters,$ background=white,color=black jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,postscript_dir+'Theta_Phi_Scatter'+label_string+'.png',$ jpgimg,reds,greens,blues ;E_THETA_SCATTER************************************************** device,set_resolution=[1400,1200] !p.charsize=1 !p.charthick=1 !x.thick=2 !y.thick=2 !p.thick=1 !p.noerase=0 !p.multi=[0,2,2] items_e=['e','e1','e2'] sym_e=[2,1,5] colors_e=[black,blue,red] plot,frame_e,frame_theta,psym=2,$ xtitle='e',$ ytitle='Ellipticity angle that maximizes e1 (Theta)',$ title='Scatter Plot of e vs. Theta for '+label_string,$ background=white,color=black plot,frame_e1,frame_theta,psym=2,$ xtitle='e1',$ ytitle='Ellipticity angle that maximizes e1 (Theta)',$ title='Scatter Plot of e1 vs. Theta for '+label_string,$ background=white,color=black plot,frame_e2,frame_theta,psym=2,$ xtitle='e',$ ytitle='Ellipticity angle that maximizes e1 (Theta)',$ title='Scatter Plot of e2 vs. Theta for '+label_string,$ background=white,color=black jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,postscript_dir+'Theta_E_Scatter'+label_string+'.png',$ jpgimg,reds,greens,blues erase ;E_PHI_SCATTER**************************************************** !p.multi=[0,2,2] plot,frame_e,phi_array,psym=2,$ xtitle='e',$ ytitle='Azimuthal angle of centroid',$ title='Scatter of e vs. Phi for '+label_string,$ background=white,color=black plot,frame_e1,phi_array,psym=1,$ xtitle='e1',$ ytitle='Azimuthal angle of centroid',$ title='Scatter of e1 vs. Phi for '+label_string,$ background=white,color=black plot,frame_e2,phi_array,psym=5,$ xtitle='e2',$ ytitle='Azimuthal angle of centroid',$ title='Scatter of e2 vs. Phi for '+label_string,$ background=white,color=black jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,postscript_dir+'Phi_E_Scatter'+label_string+'.png',$ jpgimg,reds,greens,blues device,/close if covariance eq 1 then begin set_plot,'ps' If MultiWindSpeeds eq 1 then begin device,filename=Postscript_dir+'PhaseScreenMultiWinds'+Label_String+$ 'Log_depth_'+string(strtrim(depth,2))+'.ps' endif else if MultiWindSpeeds eq 2 then begin device,filename=Postscript_dir+'PhaseScreenNoWind'+Label_String+$ 'Log_depth'+string(strtrim(depth,2))+'.ps' endif loadct,39 ;load color table 39 device,/color ;allow color on the postscript device,ysize=8.5,/inches ;Height of plot in y device,xsize=10.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='0000FF'x !P.CHARSIZE=.7 !P.THICK=4. items_time=['t=0'] sym_time=[1] items_k_tran=['Data Points','n!e-7/32!n','n!e-13/64!n','n!e-27/128!n'] sym_k_tran=[1,0,0,0] line_k_tran=[0,1,2,3] items_k_lon=['Data Points','n!e-1/3!n'] sym_k_long=[0,1] colors_time=[0] ;*************FIRST PLOT !P.multi=[0,1,1] plot,Cxt_means(*,0),psym=1,title="Transverse correlation along X direction",$ xtitle="n",ytitle="Covariance",$ background=white,color=black,$ yrange=[-1,1],$ POSITION=[0.10,0.10,0.9,0.9] oplot,Cxl_means(*,0),psym=2,color=80 oplot,Cyt_means(*,0),psym=4,color=160 oplot,Cyl_means(*,0),psym=5,color=240 itemsC=['Transverse X', 'Longitudinal X', 'Transverse Y', 'Longitudinal Y'] psymC=[1,2,4,5] colorsC=[0,80,160,240] legend,itemsC,psym=psymC,colors=colorsC,pos=[0.8,0.8],/norm xyouts,.2,.2,'xwinds='+strtrim(xwinds,2),/norm xyouts,.2,.3,'ywinds='+strtrim(ywinds,2),/norm device,/close endif set_plot,'x' stop end