;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 phase_screen_reconstruction, Date=Date,$ xwindspeed=xwindspeed,ywindspeed=ywindspeed,$ depth=depth,k=k,phase_dim=phase_dim If (not keyword_set(Date)) then Date='11' If (not keyword_set(xwindspeed)) then xwindspeed=10 If (not keyword_set(ywindspeed)) then ywindspeed=10 If (not keyword_set(depth)) then depth=0 If (not keyword_set(k)) then k=100 If (not keyword_set(phase_dim)) then phase_dim=long(1024) 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 kstr=valid_num(k,knum) 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) red='FF0000'x outer_scale=float(phase_dim*.17/knum) print,outer_scale ;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/1_31_figs/' stats_directory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$ '/ShackHartman/stats_1_31/' ;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 readcol,image_dir+'TestCalImage_stats.txt',x_grid_coor,y_grid_coor,$ format='x,f,f' ;GET THE WAVEFRONT SLOPES. ;They are stored as Columns in see.x and see.y readcol,atmosphere_dir+'see.x',x_slopes,format='f' readcol,atmosphere_dir+'see.y',y_slopes,format='f' ;TURN THESE SLOPES INTO 2-d arrays see_length=n_elements(x_slopes) see_index=lindgen(see_length) see_sub_dim=sqrt(see_length) x_slopes_arr=fltarr(see_sub_dim,see_sub_dim) y_slopes_arr=fltarr(see_sub_dim,see_sub_dim) x_slopes_arr(see_index mod see_sub_dim, see_index/see_sub_dim)=$ x_slopes(see_index) y_slopes_arr(see_index mod see_sub_dim, see_index/see_sub_dim)=$ y_slopes(see_index) ;spots is the fidicual, non_spots is non_fiducial ;MAKE SUB ARRAY FOR SLOPES ACROSS SUBAPERTURE x_slopes_sub_arr=fltarr(num_of_cor,25,25) y_slopes_sub_arr=fltarr(num_of_cor,25,25) x_slopes_nofid_sub_arr=fltarr(num_of_cor,25,25) y_slopes_nofid_sub_arr=fltarr(num_of_cor,25,25) magnitude=fltarr(num_of_cor,25,25) phi=fltarr(num_of_cor,25,25) ;FFTs xfft_array=complexarr(num_of_cor,num,num) yfft_array=complexarr(num_of_cor,num,num) xfft_nofid_array=complexarr(num_of_cor,num,num) yfft_nofid_array=complexarr(num_of_cor,num,num) ;WINDSPEED TO MOVE PHASE SCREEN ACROSS SUBAPERTURE ;FILL THE ARRAYS FOR LATER PROCESSING for increment=0,num_of_cor-1 do begin ;CORRELATIONS x_segment=increment*xwindspeed y_segment=increment*ywindspeed ;SLOPES x_slopes_sub_arr(increment,*,*)=x_slopes_arr[$ 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(increment,*,*)=y_slopes_arr[$ phase_dim-offset-x_segment:phase_dim-x_segment-1-offset+num,$ phase_dim-offset-y_segment:phase_dim-y_segment-1-offset+num] x_slopes_nofid_sub_arr(increment,*,*)=x_slopes_sub_arr(increment,*,*) y_slopes_nofid_sub_arr(increment,*,*)=y_slopes_sub_arr(increment,*,*) x_slopes_sub_arr(nonfiducial+increment*grid_spots)=0 y_slopes_sub_arr(nonfiducial+increment*grid_spots)=0 xfft_array(increment,*,*)=fft(x_slopes_sub_arr(increment,*,*)) yfft_array(increment,*,*)=fft(y_slopes_sub_arr(increment,*,*)) xfft_nofid_array(increment,*,*)=fft(x_slopes_nofid_sub_arr(increment,*,*)) yfft_nofid_array(increment,*,*)=fft(y_slopes_nofid_sub_arr(increment,*,*)) endfor ;Correlations timee1=fltarr(num_of_cor-depth,num,num) timee2=fltarr(num_of_cor-depth,num,num) timee1num=fltarr(num_of_cor-depth,num,num) timee2num=fltarr(num_of_cor-depth,num,num) timeeden=fltarr(num_of_cor-depth,num,num) timee1_arr=fltarr(num,num) timee2_arr=fltarr(num,num) frame_e1=fltarr(num_of_cor) frame_e2=fltarr(num_of_cor) frame_e=fltarr(num_of_cor) frame_thetap=fltarr(num_of_cor) frame_thetam=fltarr(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' set_plot,'ps' device,filename=Postscript_dir+'Test_Fourier_Reconstruction.ps' 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=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(num_of_cor-depth,num,num) timee2=fltarr(num_of_cor-depth,num,num) timee1num_h=fltarr(num_of_cor-depth,num,num) timee2num_h=fltarr(num_of_cor-depth,num,num) timeeden_h=fltarr(num_of_cor-depth,num,num) timee1num_l=fltarr(num_of_cor-depth,num,num) timee2num_l=fltarr(num_of_cor-depth,num,num) timeeden_l=fltarr(num_of_cor-depth,num,num) timee1_arr=fltarr(num,num) timee2_arr=fltarr(num,num) ;for high pass filter frame_e1_h=fltarr(num_of_cor) frame_e2_h=fltarr(num_of_cor) frame_e_h=fltarr(num_of_cor) frame_thetap_h=fltarr(num_of_cor) frame_thetam_h=fltarr(num_of_cor) frame_e_den_h=fltarr(num_of_cor) ;for low pass filter frame_e1_l=fltarr(num_of_cor) frame_e2_l=fltarr(num_of_cor) frame_e_l=fltarr(num_of_cor) frame_thetap_l=fltarr(num_of_cor) frame_thetam_l=fltarr(num_of_cor) frame_e_den_l=fltarr(num_of_cor) x_offset_hor_calc_h=complexarr(num_of_cor,num,num) y_offset_hor_calc_h=complexarr(num_of_cor,num,num) x_offset_hor_calc_l=complexarr(num_of_cor,num,num) y_offset_hor_calc_l=complexarr(num_of_cor,num,num) x_nofidoffset_hor_calc_h=complexarr(num_of_cor,num,num) y_nofidoffset_hor_calc_h=complexarr(num_of_cor,num,num) x_nofidoffset_hor_calc_l=complexarr(num_of_cor,num,num) y_nofidoffset_hor_calc_l=complexarr(num_of_cor,num,num) for file_number=0,num_of_cor-depth-1 do begin ;Average the fft fiducial in place 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) ;Average the fft no fiducial in place x_nofidfft=complexarr(num,num) y_nofidfft=complexarr(num,num) x_nofidlow_fft=complexarr(num,num) x_nofidlow_fft_re=complexarr(num,num) y_nofidlow_fft=complexarr(num,num) y_nofidlow_fft_re=complexarr(num,num) y_nofidhigh_fft_re=complexarr(num,num) x_nofidhigh_fft=complexarr(num,num) y_nofidhigh_fft=complexarr(num,num) x_nofidhigh_fft_re=complexarr(num,num) x_nofidfft(*,*)=xfft_nofid_array(file_number,*,*) y_nofidfft(*,*)=yfft_nofid_array(file_number,*,*) x_nofidfft_re=complexarr(num,num) y_nofidfft_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) x_nofidfft_re(subn(n),subm(m))=x_nofidfft(n,m) y_nofidfft_re(subn(n),subm(m))=y_nofidfft(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)) x_nofidhigh_fft_re=x_nofidfft_re y_nofidhigh_fft_re=y_nofidfft_re x_nofidhigh_fft_re((12-kval):(12+kval),(12-kval):(12+kval))=0 y_nofidhigh_fft_re((12-kval):(12+kval),(12-kval):(12+kval))=0 x_nofidlow_fft_re((12-kval):(12+kval),(12-kval):(12+kval))=$ x_nofidfft_re((12-kval):(12+kval),(12-kval):(12+kval)) y_nofidlow_fft_re((12-kval):(12+kval),(12-kval):(12+kval))=$ y_nofidfft_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)) x_nofidlow_fft(n,m)=x_nofidlow_fft_re(subn(n),subm(m)) y_nofidlow_fft(n,m)=y_nofidlow_fft_re(subn(n),subm(m)) x_nofidhigh_fft(n,m)=x_nofidhigh_fft_re(subn(n),subm(m)) y_nofidhigh_fft(n,m)=y_nofidhigh_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) x_nofidoffset_hor_calc_h(file_number,*,*)=fft(x_nofidhigh_fft,/inverse) y_nofidoffset_hor_calc_h(file_number,*,*)=fft(y_nofidhigh_fft,/inverse) x_nofidoffset_hor_calc_l(file_number,*,*)=fft(x_nofidlow_fft,/inverse) y_nofidoffset_hor_calc_l(file_number,*,*)=fft(y_nofidlow_fft,/inverse) test=1 if (test eq 1 and file_number eq 0) then begin erase !p.noerase=1 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),color=60 xyouts,gridx(624)-100,gridy(624)+125,'length='+$ strtrim(xoff(624),2) xoff=fltarr(625) yoff=fltarr(625) xoff(*)=x_nofidoffset_hor_calc_h(file_number,*,*) yoff(*)=y_nofidoffset_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),color=0 xyouts,.1,.95,'Purple is with fiducial; Black is without',$ color=0,/normal endif endfor endfor xpower=avg(abs(xfft_array),0) ypower=avg(abs(yfft_array),0) xpower_nofid=avg(abs(xfft_nofid_array),0) ypower_nofid=avg(abs(yfft_nofid_array),0) erase plot,xpower_nofid,$ title="Power in x offsets with/without fiducial region",xtitle='k values',$ ytitle='Power',color=60,background=white oplot,xpower,color=0 xyouts,.1,.95,'Purple is with fiducial; Black is without',color=0,/normal erase plot,ypower_nofid,$ title="Power in x offsets with/without fiducial region",xtitle='k values',$ ytitle='Power',color=60,background=white oplot,ypower,color=0 xyouts,.1,.95,'Purple is with fiducial; Black is without',color=0,/normal loadct,0 ypower_re=complexarr(num,num) xpower_re=complexarr(num,num) ypower_nofid_re=complexarr(num,num) xpower_nofid_re=complexarr(num,num) ;Rearrange for n=0,num-1 do begin for m=0,num-1 do begin xpower_re(n,m)=xpower(subn(n),subm(m)) ypower_re(n,m)=ypower(subn(n),subm(m)) xpower_nofid_re(n,m)=xpower_nofid(subn(n),subm(m)) ypower_nofid_re(n,m)=ypower_nofid(subn(n),subm(m)) endfor endfor erase shade_surf,congrid(xpower_re,600,600),background=white,$ color=black,xtitle='Kx',ytitle='Ky',$ position=[.1,.1,.9,.9],/normal xyouts,.1,.95,'X offsets Power with fiducial',color=black,/normal erase shade_surf,congrid(xpower_nofid_re,600,600),background=white,$ color=black,xtitle='Kx',ytitle='Ky',$ position=[.1,.1,.9,.9],/normal xyouts,.1,.95,'X offsets Power without fiducial',color=black,/normal erase shade_surf,congrid(xpower_re-xpower_nofid_re,600,600),background=white,$ color=black,xtitle='Kx',ytitle='Ky',$ position=[.1,.1,.9,.9],/normal xyouts,.1,.95,'Difference in X offsets Power with-without fiducial',color=black,/normal erase shade_surf,congrid(ypower_re,600,600),background=white,$ color=black,xtitle='Kx',ytitle='Ky',$ position=[.1,.1,.9,.9],/normal xyouts,.1,.95,'Y offsets Power with fiducial',color=black,/normal erase shade_surf,congrid(ypower_nofid_re,600,600),background=white,$ color=black,xtitle='Kx',ytitle='Ky',$ position=[.1,.1,.9,.9],/normal xyouts,.1,.95,'Y offsets Power without fiducial',color=black,/normal erase shade_surf,congrid(ypower_re-ypower_nofid_re,600,600),background=white,$ color=black,xtitle='Kx',ytitle='Ky',$ position=[.1,.1,.9,.9],/normal xyouts,.1,.95,'Difference in Y offsets Power with-without fiducial',color=black,/normal erase device,/close set_plot,'x' stop end