;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_batch,$ xFrameDrift=xFrameDrift,yFrameDrift=yFrameDrift,$ xWindSpeed=xWindSpeed,yWindSpeed=yWindSpeed,$ depth=depth,phase_dim=phase_dim,$ low=low,lcut=lcut 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=5. ;m/sec If (not keyword_set(yWindSpeed)) then yWindSpeed=5. ;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=1000.0 ;LABELS: *************************************************************** xwinds=[0,5,10,15,20] ywinds=[0,5,10,15,20] 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) ;FRAMESPEEDS FOR DRAGGING PHASE SCREEN ACROSS CWFS xFrameDrift_m=xFrameDrift*.17 yFrameDrift_m=yFrameDrift*.17 FrameDrifts=[xFrameDrift,yFrameDrift] max_FrameDrift=max(FrameDrifts) ;KEY TO USE THE FIDUCIAL IN THE CALCULATION UseFid = 0 if UseFid eq 0 then begin label_string = 'NoFiducial_' endif else begin label_string = 'Fiducial_' endelse ;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 ;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 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 ;************************************************************************ ;TAKE THE DATA IN A CUBE TO USE THROUGHOUT phase_files=file_search(screen_directory+'phase_'+strtrim(lowstr,2)+'_'+$ strtrim(lcutstr,2)+'*.out',$ count=num_phase_files) xslopes=fltarr(num_phase_files,phase_dim,phase_dim) yslopes=fltarr(num_phase_files,phase_dim,phase_dim) total_num_of_cor=num_of_cor*num_phase_files 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(screen_num,1:1022,1:1022)=(screen(2:1023,1:1022)-screen(0:1021,1:1022))/2 yslopes(screen_num,1:1022,1:1022)=(screen(1:1022,2:1023)-screen(1:1022,0:1021))/2 yslopes(screen_num,*)=yslopes(screen_num,*)/$ sqrt(mean(yslopes(screen_num,*)^2)) xslopes(screen_num,*)=xslopes(screen_num,*)/$ sqrt(mean(xslopes(screen_num,*)^2)) endfor ;*********************************************************************** ;LOOP THROUGH PLOTS AND WINDS ;png_val=0 POSTSCRIPT ;png_val=1 PNGs png_val=0 if png_val eq 0 then begin set_plot,'ps' label_string=label_string+'low_'+lowinmeters+'_high_'+highinmeters ps_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/final_paper_figs/' device,filename=ps_dir+'Hist_Scat'+label_string+'.ps' !p.multi=[0,5,5] 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 blue='FF0000'x red='0000FF'x endif for plot_val=1, 6 do begin if png_val eq 1 then begin if plot_val gt 0 and plot_val lt 5 then begin set_plot,'z' erase device, set_font='Courier' device,set_resolution=[1400,1200] !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,5,5] endif else begin !p.charsize=1 !p.charthick=1 !x.thick=2 !y.thick=2 !p.thick=1 !p.noerase=0 items_e=['e','e1','e2'] sym_e=[2,1,5] colors_e=[black,blue,red] endelse endif for XWindGlobal=0,N_elements(xwinds)-1 do begin for YWindGlobal=0,N_elements(ywinds)-1 do begin xWindSpeed=xwinds(XWindGlobal) yWindSpeed=ywinds(YWindGlobal) print,xWindSpeed,yWindspeed xwind=strtrim(string(xWindSpeed,format='(%"%10.1f")'),2) ywind=strtrim(string(yWindSpeed,format='(%"%10.1f")'),2) label_string='_xwind_'+xwind+'_ywind_'+ywind ;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; ;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 ;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) 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_theta=fltarr(total_num_of_cor) frame_e_den=fltarr(total_num_of_cor) for screen_num=0,num_phase_files-1 do begin for increment=0,num_of_cor-1 do begin tot_increment=num_of_cor*screen_num+increment ;CORRELATIONS x_segment=increment*xFrameDrift y_segment=increment*yFrameDrift ;SLOPES AVERAGED WITH THE WIND Cells=[numXCells,numYCells] numCells=max(Cells) xWindSlopes=fltarr(numCells,num,num) yWindSlopes=fltarr(numCells,num,num) iniXoffset=phase_dim-offset-x_segment iniYoffset=phase_dim-offset-y_segment for WindInc=0, numCells-1 do begin xWindoff=floor((float(numXCells)/float(numCells))*float(WindInc)) yWindoff=floor((float(numYCells)/float(numCells))*float(WindInc)) xWindSlopes(WindInc,*,*)=xslopes[screen_num,$ iniXoffset-xWindOff:iniXoffset-xWindOff+num-1,$ iniYoffset-yWindOff:iniYoffset-yWindOff+num-1] yWindSlopes(WindInc,*,*)=yslopes[screen_num,$ 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,*,*) ;USE FIDUCIAL OR NOT if UseFid eq 1 then begin xtemp(nonfiducial)=0 ytemp(nonfiducial)=0 endif 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) x_slopes_sub_arr(tot_increment,*,*)=xtemp y_slopes_sub_arr(tot_increment,*,*)=ytemp timee1num(tot_increment,*,*)=x_slopes_sub_arr(tot_increment,*,*)*$ x_slopes_sub_arr(tot_increment+depth,*,*)-$ y_slopes_sub_arr(tot_increment,*,*)*$ y_slopes_sub_arr(tot_increment+depth,*,*) timee2num(tot_increment,*,*)=(x_slopes_sub_arr(tot_increment,*,*)*$ y_slopes_sub_arr(tot_increment+depth,*,*)$ +x_slopes_sub_arr(tot_increment+depth,*,*)$ *y_slopes_sub_arr(tot_increment,*,*)) timeeden(tot_increment,*,*)=(x_slopes_sub_arr(tot_increment,*,*)^2+$ y_slopes_sub_arr(tot_increment,*,*)^2) timee1temp=timee1num(tot_increment,*,*) timee2temp=timee2num(tot_increment,*,*) timeedtemp=timeeden(tot_increment,*,*) frame_e_den=mean(timeedtemp(xfit)) frame_e1(tot_increment)=mean(timee1temp(xfit))/$ mean(timeedtemp(xfit)) frame_e2(tot_increment)=mean(timee2temp(xfit))/$ mean(timeedtemp(xfit)) frame_e(tot_increment)=sqrt(frame_e1(tot_increment)^2+frame_e2(tot_increment)^2) signe1=abs(frame_e1(tot_increment))/frame_e1(tot_increment) signe2=abs(frame_e2(tot_increment))/frame_e2(tot_increment) frame_theta(tot_increment)=$ signe2*asin(sqrt(.5*(1-signe1*sqrt(1-frame_e2(tot_increment)^2/$ frame_e(tot_increment)^2)))) ;Recalculate e1 and e2 with new cosine angle frame_e1(tot_increment)=frame_e(tot_increment)*$ cos(2*frame_theta(tot_increment)) frame_e2(tot_increment)=frame_e(tot_increment)*$ sin(2*frame_theta(tot_increment)) endfor endfor 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 if plot_val eq 1 then begin ;THETA_HISTOGRAM***************************************************** plot,loc_theta+thetabinsize/2,thetahist,psym=10,xtitle='Theta',$ ytitle='Number of Events',xrange=[-1.6,1.6],$ title=label_string,$ background=white,color=black endif else if plot_val eq 2 then begin ;E_HISTOGRAM******************************************************** plot,loc_e,ehist,psym=10,xtitle='e',$ ytitle='Number of Events',$ title=label_string,$ background=white,color=black endif else if plot_val eq 3 then begin ;PHI_HISTOGRAM****************************************************** plot,loc_phi+phibinsize/2,phihist,psym=10,xtitle='Theta',$ ytitle='Number of Events',$ title='Histogram of Centroid Angle for'+label_string,$ background=white,color=black endif else if plot_val eq 4 then begin ;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'+label_string,$ background=white,color=black endif else if plot_val eq 5 then begin ;E_THETA_SCATTER************************************************** 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 oplot,frame_e1,frame_theta,psym=2,color=30 oplot,frame_e2,frame_theta,psym=2,color=60 endif else if plot_val eq 6 then begin ;E_PHI_SCATTER**************************************************** 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 oplot,frame_e1,phi_array,psym=1,color=30 oplot,frame_e2,phi_array,psym=5,color=60 endif endfor endfor if png_val eq 1 then begin if plot_val eq 1 then begin ;THETA_HISTOGRAM***************************************************** jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,postscript_dir+'Theta_Histogram'+label_string+'.png',$ jpgimg,reds,greens,blues endif else if plot_val eq 2 then begin ;E_HISTOGRAM******************************************************** jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,postscript_dir+'E_Histogram_'+label_string+'.png',$ jpgimg,reds,greens,blues endif else if plot_val eq 3 then begin ;PHI_HISTOGRAM****************************************************** jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,postscript_dir+'Phi_Histogram'+label_string+'.png',$ jpgimg,reds,greens,blues endif else if plot_val eq 4 then begin ;PHI_THETA_SCATTER****************************************** jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,postscript_dir+'Theta_Phi_Scatter'+label_string+'.png',$ jpgimg,reds,greens,blues endif else if plot_val eq 5 then begin ;E_THETA_SCATTER************************************************** jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,postscript_dir+'Theta_E_Scatter'+label_string+'.png',$ jpgimg,reds,greens,blues endif else if plot_val eq 6 then begin ;E_PHI_SCATTER**************************************************** jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,postscript_dir+'Phi_E_Scatter'+label_string+'.png',$ jpgimg,reds,greens,blues endif endif endfor device,/close set_plot,'x' stop end