;tv_phase_screen_wind ; ;PURPOSE: ; This script will take in a phase screen generated by turb2d.c (written ; by Garret Jernigan) and obtain the Fourier Spectrum for each frame with ; a wind drifting the phase screen across the aperture. ; ; ;CALLING SEQUENCE: ; tv_phase_screen_wind ; ;KEYWORDS: ; xFrameDrift, yFrameDrift: ; The speed at which the entire phase screen is shifted before ; a new sample is taken ; xWindspeed, yWindspeed: ; The speed at which the frozen screen will be dragged across ; the aperture. The higher the speed, the more angles of arrival ; averaged for the calculation pro tv_phase_screen_wind,$ 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=80. ;m/sec If (not keyword_set(yWindSpeed)) then yWindSpeed=0. ;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 ;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 ;LABELS: *************************************************************** if UseFid eq 0 then begin label_string = '_NoFiducial_' endif else begin label_string = '_Fiducial_' endelse 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 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) label_string=label_string+'low_'+lowinmeters+'_high_'+highinmeters+$ '_xwind_'+xwind+'_ywind_'+ywind ;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(xWindSpeed*exp_time/0.17)+1; numYcells=ceil(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) 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) 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) set_plot,'x' !p.multi=[0,1,2] !p.noerase=1 window,xsize=1000,ysize=1000 ;********************************************************** 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 ;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=floor((float(numXCells)/float(numCells))*float(WindInc)) yWindoff=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) ;PLOT WHAT THE INTENSITY DISTRIBUTION WOULD LOOK LIKE 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 x_slopes_sub_arr(tot_increment,*,*)=xtemp y_slopes_sub_arr(tot_increment,*,*)=ytemp ;PLOT WHAT THE INTENSITY DISTRIBUTION WOULD LOOK LIKE plot, x_slopes_sub_arr(tot_increment,*,*),$ y_slopes_sub_arr(tot_increment,*,*),$ psym=2,xrange=[-3,3],yrange=[-3,3],$ position=[0.20,0.35,0.80,.95] ;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,*,*)) ;Rearrange to plot 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 xpower=abs(x_fft_re)^2 ypower=abs(y_fft_re)^2 plot,xpower(12:24,12),psym=1,color=white,$ position=[0.05,0.05,0.95,.3] oplot,ypower(12,12:24),psym=5,color=25555 wait,2 endfor endfor device,/close set_plot,'x' stop end