;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 phase_screen_wind_rms_power_spectrum,$
		   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=40   ;pixels/frame
If (not keyword_set(yFrameDrift)) then yFrameDrift=40   ;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 = 1

;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)

 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))

endfor
;
;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,'ps'
label_string=label_string+'low_'+lowinmeters+'_high_'+highinmeters
device,filename=Postscript_dir+'rms_power_spect'+label_string+'.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=10.0,/inches        ;Width of plot in x
device,yoffset=1.0,/inches      ;Y position of lower left corner
!p.multi=[0,2,3]
white='FFFFFF'x
black='000000'x         

symx=[1,5]
colorsx=[black,31]
itemsx=['X Power Spectrum','Y Power Spectrum']
xwinds=[0,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75]
NoXwinds=N_elements(xwinds)
XoffrmsArr=fltarr(NoXWinds)
YoffrmsArr=fltarr(NoXWinds)

plot_no=0

;**********************************************************
for XWindGlobal= 0, NoXWinds-1 do begin
 xWindSpeed=xwinds(XWindGlobal)
 yWindSpeed=0 
 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

 xpowerarr=fltarr(total_num_of_cor,13)
 ypowerarr=fltarr(total_num_of_cor,13)

 if plot_no mod 3 eq 0 then begin
   plot,[-4,4],[-4,4],/xstyle,/ystyle,xrange=[-4,4],yrange=[-4,4],/nodata,$
	color=black,background=white,position=[0.05,0.70,0.30,0.95]
 endif else if plot_no mod 3 eq 1 then begin
   plot,[-4,4],[-4,4],/xstyle,/ystyle,xrange=[-4,4],yrange=[-4,4],/nodata,$
        color=black,background=white,position=[0.05,0.40,0.30,0.65]
 endif else if plot_no mod 3 eq 2 then begin
   plot,[-4,4],[-4,4],/xstyle,/ystyle,xrange=[-4,4],yrange=[-4,4],/nodata,$
        color=black,background=white,position=[0.05,0.10,0.30,0.35]
 endif
 for screen_num=0,num_phase_files-1 do begin

  ;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
    if tot_increment lt 15 then begin
    oplot, x_slopes_sub_arr(tot_increment,*,*),$
          y_slopes_sub_arr(tot_increment,*,*),psym=2
    endif
    ;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,*,*))
  
    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=fltarr(13,1)
    ypower=fltarr(13,1)
    xpower(*)=abs(x_fft_re(12:24,12)^2)
    ypower(*)=abs(y_fft_re(12,12:24)^2)
    xpowerarr(tot_increment,*)=xpower
    ypowerarr(tot_increment,*)=ypower

  endfor

 endfor

    xoffrms=sqrt(mean(x_slopes_sub_arr(where(x_slopes_sub_arr ne 0))^2))
    yoffrms=sqrt(mean(y_slopes_sub_arr(where(y_slopes_sub_arr ne 0))^2))
    XoffrmsArr(XWindGlobal)=xoffrms
    YoffrmsArr(XWindGlobal)=yoffrms
    if plot_no mod 3 eq 0 then begin
     plot,avg(xpowerarr,0),psym=1,color=black,background=white,$
	position=[0.4,0.7,0.95,.95],$
	title=xwind 
     oplot,avg(ypowerarr,0),psym=5,color=31
     xyouts,0.8,0.87,'xrms='+strtrim(xoffrms,2),/norm
     xyouts,0.8,0.84,'y rms ='+strtrim(yoffrms,2),/norm
    endif else if plot_no mod 3 eq 1 then begin
     plot,avg(xpowerarr,0),psym=1,color=black,background=white,$
        position=[0.4,0.4,0.95,.65],$
	title=xwind
     oplot,avg(ypowerarr,0),psym=5,color=31
     xyouts,0.8,0.60,'xrms='+strtrim(xoffrms,2),/norm
     xyouts,0.8,0.55,'y rms ='+strtrim(yoffrms,2),/norm
    endif else if plot_no mod 3 eq 2 then begin
     plot,avg(xpowerarr,0),psym=1,color=black,background=white,$
        position=[0.4,0.1,0.95,.35],$
	title=xwind
     oplot,avg(ypowerarr,0),psym=5,color=31
     xyouts,0.8,0.30,'xrms='+strtrim(xoffrms,2),/norm
     xyouts,0.8,0.25,'y rms ='+strtrim(yoffrms,2),/norm
     legend,itemsx,psym=symx,colors=colorsx,textcolors=black,$
	position=[0.05,0.05],/normal
   endif
    plot_no=plot_no+1
    
endfor

!p.multi=[0,1,1]
plot,XoffrmsArr^2,psym=2,title='Mean Squares of the offsets with wind',$
	xtitle='Exposure Time * Wind speed / Aperture Size',$
	ytitle='Mean Square'
oplot,YoffrmsArr^2,psym=1
itemsrms=['X offsets','Y offsets']
psymrms=[2,1]
legend,itemsrms,psym=psymrms,textcolors=black

device,/close
set_plot,'x'
                                                                           
stop
                                                                                end