;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,[windspeed,[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_angle_histogram_no_text, Date=Date,$ xwindspeed=xwindspeed,ywindspeed=ywindspeed,$ depth=depth,phase_dim=phase_dim,$ low=low,lcut=lcut If (not keyword_set(Date)) then Date='11' If (not keyword_set(xwindspeed)) then xwindspeed=25 If (not keyword_set(ywindspeed)) then ywindspeed=25 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 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_windspeed 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 ;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) outer_scale=float(phase_dim*.17/low) print,outer_scale 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) ;FILE IO stats_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$ '/ShackHartman/coo/' atmosphere_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/lsst_sim/atmosphere/' image_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/ShackHartmann_Cal/' Postscript_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/phase_screen_figs/' stats_directory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$ '/ShackHartman/stats_1_31/' screen_directory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/lsst_sim/atmosphere/' ;HOW TO DEAL WITH AVERAGE SUBTRACTION ;Indices and such no_bins=8 a=findgen(13)+12 b=findgen(12) subn=[a,b] subm=subn avg_sub=0 if avg_sub eq 1 then begin evals_file_name=screen_directory+'ellipticity_avg_sub_low_'+lowinmeters+$ '_high_'+highinmeters+'.txt' endif else if avg_sub eq 2 then begin evals_file_name=screen_directory+'ellipticity_kx0_ky0_sub_low_'+lowinmeters+$ '_high_'+highinmeters+'.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 ;ARRAYS OF FILES offset_files=file_search(stats_directory+'*5132*sub_offsets.txt',$ count=num_stats_files) get_lun,test_out openw,test_out,'./test_out.txt' 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 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) 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) ;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) yslopes(1:1022,1:1022)=(screen(2:1023,1:1022)-screen(0:1021,1:1022))/2 xslopes(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 ;CORRELATIONS x_segment=increment*xwindspeed y_segment=increment*ywindspeed ;SLOPES x_slopes_sub_arr(tot_increment,*,*)=xslopes[$ phase_dim-offset-x_segment:phase_dim-x_segment-1-offset+num,$ phase_dim-offset-y_segment:phase_dim-y_segment-1-offset+num] y_slopes_sub_arr(tot_increment,*,*)=yslopes[$ phase_dim-offset-x_segment:phase_dim-x_segment-1-offset+num,$ phase_dim-offset-y_segment:phase_dim-y_segment-1-offset+num] ;25x25 grid, make one-dimensional and apply fiducial 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 ;Avg_sub=1 means subtract the average offset of a given frame ;Avg_sub=2 means subtract the kx=0,ky=0 componenet of FFT ;Avg_sub=3 means subtract a scaled offset from the distorted grid if avg_sub eq 1 then begin xfit=where(xtemp ne 0) yfit=where(ytemp ne 0) xtemp(xfit)=xtemp(xfit)-mean(xtemp(xfit)) ytemp(yfit)=ytemp(yfit)-mean(ytemp(yfit)) 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) readcol,offset_files(0),xcoord,ycoord,format='x,f,x,f,x' set_plot,'ps' device,filename=Postscript_dir+'Phase_Screen_Multi_Successive_Filter_Individual.ps' loadct,0 ;load color table 39 device,/color ;allow color on the postscript device,ysize=8.5,/inches ;Height of plot in y device,xsize=11.0,/inches ;Width of plot in x device,yoffset=1.0,/inches ;Y position of lower left corner 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_k0=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,*,*) 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))) 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) frame_e_k0(file_number)=frame_e_h(file_number) printf,evals,FORMAT='(%"%f\t%f")',$ frame_theta(file_number),frame_e_k0(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 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 begin frame_e_den_l=mean(timeedtemp_h(where(timee1temp_h 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))) frame_e_l(file_number)=sqrt(frame_e1_l(file_number)^2+frame_e2_l(file_number)^2) endelse 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 ;---------------------------------------------------------- ;Plot e1 for each kval items_e1=['e1 high pass','e1 low pass'] sym_e1=[2,1] plot,frame_e1_h,psym=2,title='low passed e1 and e2 for kmin='+strtrim(kval,2) oplot,frame_e1_l,psym=1 legend,items_e1,psym=sym_e1,position=[.8,.8],/normal items_e2=['e2 high pass','e2 low pass'] sym_e2=[2,1] plot,frame_e2_h,psym=2,$ color=10,$ title='high and low passed e2 for kmin='+strtrim(kval,2) oplot,frame_e2_l,psym=1,$ color=20 legend,items_e2,psym=sym_e2,position=[.8,.8],/normal ;Plot e for each KVAL items_e=['e high pass','e low pass'] sym_e=[2,1] plot,frame_e_h,psym=2,$ color=10,$ title='e for kmin='+strtrim(kval,2) oplot,frame_e_l,psym=1,$ color=20 legend,items_e,psym=sym_e,position=[.8,.8],/normal ;Plot samples of the low and high passed offset x_offset_h=fltarr(num,num) y_offset_h=fltarr(num,num) x_offset_l=fltarr(num,num) y_offset_l=fltarr(num,num) for row_i=0,24 do begin for col_j=0,24 do begin x_offset_h(row_i,col_j)=avg(x_offset_hor_calc_h(*,row_i,col_j)) y_offset_h(row_i,col_j)=avg(y_offset_hor_calc_h(*,row_i,col_j)) x_offset_l(row_i,col_j)=avg(x_offset_hor_calc_l(*,row_i,col_j)) y_offset_l(row_i,col_j)=avg(y_offset_hor_calc_l(*,row_i,col_j)) endfor endfor x_offset_h(624)=max(x_offset_h) y_offset_h(624)=0 partvelvec,x_offset_h(fitted),y_offset_h(fitted),$ gridx(fitted),gridy(fitted),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)+25,'length='+$ strtrim(x_offset_h(624),2) x_offset_l(624)=max(x_offset_l) y_offset_l(624)=0 partvelvec,x_offset_l(fitted),y_offset_l(fitted),$ gridx(fitted),gridy(fitted),yrange=[150,800],$ xrange=[200,850],position=[0.1,0.1,0.9,0.9],$ title='Low Passed Average Offsets for all images for k='+$ strtrim(kval,2) xyouts,gridx(624)-100,gridy(624)+25,'length='+$ strtrim(x_offset_l(624),2) ;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))) 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))) endfor device,/close ;----------------------CLOSE EVALS FILE close, evals free_lun, evals ;------------------------------------------Begin Multi Plot device,filename=Postscript_dir+$ 'Ellipticity_Filtered_Cumulative_outer_'$ +strtrim(lowinmeters,2)+'_m_inner_'+$ strtrim(highinmeters,2)+'_m.ps' ;Plot averages after all is said and done !p.multi=[0,1,1] items_all=['high pass e1','low pass e1','denominator high pass',$ 'denominator low pass'] sym_all=[2,1,4,5] e1maxs=[max(frame_e1_av_h),max(frame_e1_av_l),max(frame_e_av_den_h),$ max(frame_e_av_den_l)] e1max=max(e1maxs) e1mins=[min(frame_e1_av_h),min(frame_e1_av_l),min(frame_e_av_den_h),$ min(frame_e_av_den_l)] e1min=min(e1mins) plot,frame_e1_av_h,psym=2,title='e1 for filter up to kval for',$ xtitle='k',ytitle='e1',yrange=[e1min,e1max] oplot,frame_e1_av_l,psym=1 oplot,frame_e_av_den_h,psym=4 oplot,frame_e_av_den_l,psym=5 xyouts,.8,.6,'Outer (m)='+strtrim(lowinmeters,2),/normal,color=black xyouts,.8,.65,'Inner (m)='+strtrim(highinmeters,2),/normal,color=black xyouts,.8,.55,'Total frames='+strtrim(total_num_of_cor,2),/normal,color=black legend,items_all,psym=sym_all,position=[.7,.8],/normal items_all=['high pass e2','low pass e2','denominator high pass',$ 'denominator low pass'] sym_all=[2,1,4,5] e2maxs=[max(frame_e2_av_h),max(frame_e2_av_l),max(frame_e_av_den_h),$ max(frame_e_av_den_l)] e2max=max(e2maxs) e2mins=[min(frame_e2_av_h),min(frame_e2_av_l),min(frame_e_av_den_h),$ min(frame_e_av_den_l)] e2min=min(e2mins) !p.multi=[0,1,1] plot,frame_e2_av_h,psym=2,title='e2 for filter up to kval',$ xtitle='k',ytitle='e2',yrange=[e2min,e2max] oplot,frame_e2_av_l,psym=1 oplot,frame_e_av_den_h,psym=4 oplot,frame_e_av_den_l,psym=5 xyouts,.8,.6,'Outer (m)='+strtrim(lowinmeters,2),/normal,color=black xyouts,.8,.65,'Inner (m)='+strtrim(highinmeters,2),/normal,color=black xyouts,.8,.55,'Total frames='+strtrim(total_num_of_cor,2),/normal,color=black legend,items_all,psym=sym_all,position=[.7,.8],/normal items_all=['high pass e','low pass e','denominator high pass',$ 'denominator low pass'] sym_all=[2,1,4,5] 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) !p.multi=[0,1,1] plot,frame_e_av_h,psym=2,title='e for filter up to kval',$ xtitle='k',ytitle='e',yrange=[emin,emax] oplot,frame_e_av_l,psym=1 oplot,frame_e_av_den_h,psym=4 oplot,frame_e_av_den_l,psym=5 xyouts,.8,.6,'Outer (m)='+strtrim(lowinmeters,2),/normal,color=black xyouts,.8,.65,'Inner (m)='+strtrim(highinmeters,2),/normal,color=black xyouts,.8,.55,'Total frames='+strtrim(total_num_of_cor,2),/normal,color=black legend,items_all,psym=sym_all,position=[.7,.8],/normal thetamin=-!pi/2 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) png=1 if png eq 1 then begin set_plot,'z' erase device, set_font='Courier' device,set_resolution=[800,600] !p.charsize=.8 !p.charthick=1.2 !x.thick=2 !y.thick=2 !p.thick=.8 !p.noerase=0 white='FFFFFF'x black='000000'x !P.CHARSIZE=1.0 !P.THICK=4. endif white='FFFFFF'x black='000000'x plot,loc_theta+thetabinsize/2,thetahist,psym=10,xtitle='Theta',$ ytitle='Number of Events',xrange=[-1.6,1.6],$ background=white,color=black if png eq 1 then begin jpgimg=tvrd() tvlct,reds,greens,blues,/get if avg_sub eq 0 then begin write_png,postscript_dir+'Theta_Histogram_low_'+lowinmeters+$ '_high_'+highinmeters+'.png',$ jpgimg,reds,greens,blues endif else if avg_sub eq 1 then begin write_png,postscript_dir+'Theta_Histogram_Avg_sub'+lowinmeters+$ '_high_'+highinmeters+'.png',$ jpgimg,reds,greens,blues endif else begin write_png,postscript_dir+'Theta_Histogram_kx0_ky0_sub'+lowinmeters+$ '_high_'+highinmeters+'.png',$ jpgimg,reds,greens,blues endelse endif device,/close set_plot,'x' stop end