;PHASE_OUT_PLOT
;
;PURPOSE:
;	This script will take in a phase screen generated by turb2d.c (written
;	by Garret Jernigan) and generate a color image of that screen.  It also
;	takes the slope of the wavefronts in the x and y direction as listed
;	in the see.x and see.y files and plots them as a field.
;
;	It will also compute the auto-correlation for the available angles on 
;	the ShackHartman grid from the see.x and see.y wavefront slopes.  
;	A plot will be generated of the correlations.
;
;CALLING SEQUENCE:
;	phase_out_plot,[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 see_depth_scroll_backwards,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=0
If (not keyword_set(ywindspeed)) then ywindspeed=0
If (not keyword_set(depth)) then depth=1
If (not keyword_set(k)) then k=100
If (not keyword_set(phase_dim)) then phase_dim=long(2048)
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)

outer_scale=float(phase_dim*.17/knum)
print,outer_scale

;FILE IO
stats_dir='/a/pippin01/Volumes/u08/lsst/CPanalysis/2005-05-'+Date+$
	'/ShackHartman/coo/'
atmosphere_dir='/afs/slac.stanford.edu/u/ki/lances/lsst_sim/atmosphere/'
image_dir='/a/pippin01/Volumes/u08/lsst/CPanalysis/ShackHartmann_Cal/'
Postscript_dir='/a/pippin01/Volumes/u08/lsst/CPanalysis/2005-05-'+Date+$
        '/ShackHartman/postscripts/'
get_lun,test_out
openw,test_out,'./test_out.txt'

;GET THE WAVEFRONT SLOPES. 
;They are stored as Columns in see.x and see.y
readcol,atmosphere_dir+'see_'+strtrim(phase_dim,2)+'_'+strtrim(k,2)+'.x',x_slopes,format='f'
readcol,atmosphere_dir+'see_'+strtrim(phase_dim,2)+'_'+strtrim(k,2)+'.y',y_slopes,format='f'
readcol,image_dir+'TestCalImage_stats.txt',x_grid_coor,y_grid_coor,$
				        format='x,f,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)

;PULL OUT THE FIDUCIAL SPOTS ACROSS THE SOAR APERTURE
readcol,stats_dir+'fiducial_reference.txt',spots,format='x,x,f'
spot_arr=intarr(25,25)
spot_arr(spots)=130
complement,findgen(625),spots,non_spots


;MAKE SUB ARRAY FOR SLOPES ACROSS SUBAPERTURE
x_slopes_sub_arr=fltarr(25,25,num_of_cor)
y_slopes_sub_arr=fltarr(25,25,num_of_cor)
magnitude=fltarr(25,25,num_of_cor)
phi=fltarr(25,25,num_of_cor)

;PREPARE CORRELATION ARRAYS
Cl_arr=fltarr(num_of_cor,(2*num-1)*(num),depth)
Ct_arr=fltarr(num_of_cor,(2*num-1)*(num),depth)
mags=fltarr((2*num-1)*(num))
angles=fltarr((2*num-1)*(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]

  magnitude(*,*,increment)=sqrt(x_slopes_sub_arr(*,*,increment)^2+$
	y_slopes_sub_arr(*,*,increment)^2)
  phi(*,*,increment)=atan(y_slopes_sub_arr(*,*,increment),$
	x_slopes_sub_arr(*,*,increment))
  magnitude(non_spots+increment*grid_spots)=-10
  phi(non_spots+increment*grid_spots)=-10
  
endfor

for increment=long(0),num_of_cor-1-depth do begin

  for time=long(0), depth-1 do begin
    
    spacing_index=0
    for delta_x=-(num-1),(num-1) do begin
        
       for delta_y=0,(num-1) do begin
         
         ;THE ANGLE AND SEPARTATION FOR A GIVEN SET OF DELTAS
         grid_phi=atan(delta_y,delta_x)
         grid_sep=sqrt(delta_x^2+delta_y^2)

 	 if increment eq 0 then begin
           angles(spacing_index)=grid_phi
           mags(spacing_index)=grid_sep
         endif
         Ct=fltarr(grid_spots-abs(delta_x)*abs(delta_y))
         Cl=fltarr(grid_spots-abs(delta_x)*abs(delta_y))
         ;This index represents the number of points going into the sum for
         ;a given separation
	 c_index=0
                                                                                
         for sub_col_index=0,num-1 do begin
                                 
           for sub_row_index=0, num-1 do begin
             ;Check to see if a valid spot is located at the spot we're interrogating
             if (((sub_row_index+delta_x) lt num) and $
                ((sub_row_index+delta_x) ge 0)) then begin
               
	       first_check=num*sub_col_index+sub_row_index
	       first_spot=long(num*sub_col_index+sub_row_index+$
                          long(increment*grid_spots))
               
             if delta_x ge 0 then begin
		second_check=num*sub_col_index+sub_row_index+num*delta_y+delta_x
                second_spot=long(num*sub_col_index+sub_row_index+$
			  long((increment+time)*grid_spots)+num*delta_y+delta_x)
             endif else begin
 		second_check=num*sub_col_index+sub_row_index-num*delta_y+delta_x
		second_spot=long(num*sub_col_index+sub_row_index+$
                          long((increment+time)*grid_spots)-num*delta_y+delta_x)
	     endelse
    		 if ((second_check ge 0) and (second_check lt 625) and $
               (first_check ge 0) and (first_check lt 625)) then begin
                                                                       
                ;If it is on the grid, check to see if it has a partner
                if ((phi(first_spot) ne -10) and $
                   (phi(second_spot) ne -10))$
                then begin
		;Longitudinal will be cosine and transverse will be sine

	         Cl(c_index)=magnitude(first_spot)*$
                         cos(phi(first_spot)-grid_phi)*$
                         magnitude(second_spot)*$
                         cos(phi(second_spot)-grid_phi)
                 Ct(c_index)=magnitude(first_spot)*$
                         sin(phi(first_spot)-grid_phi)*$
                         magnitude(second_spot)*$
                         sin(phi(second_spot)-grid_phi)
                c_index=c_index+1
		
		endif
             endif
	   endif                                                                            
          endfor
                                                                                 
        endfor
	;if grid_phi eq 0 then print,c_index,spacing_index
        ;Correlations for a given delta_x and delta_y are now in place
	 correlated=where(Cl ne 0)
         ;If there were any pairs for a covariance, get there mean
         If correlated(0) ne -1 then begin
            Cl_arr(increment,spacing_index,time)=$
            mean(Cl(correlated))
            Ct_arr(increment,spacing_index,time)=$
            mean(Ct(correlated))
        endif else begin
        ;If there were no pairs, assign a system no value
            Cl_arr(increment,spacing_index,time)=!values.f_nan
            Ct_arr(increment,spacing_index,time)=!values.f_nan
        endelse
        ;Now onto the next set of spacings.  Increment delta_y
        spacing_index=spacing_index+1
                                                                                 
       endfor
  
    endfor

  endfor
  
endfor

close,test_out
free_lun,test_out

Cl_means=fltarr((2*num-1)*(num),depth)
Ct_means=fltarr((2*num-1)*(num),depth)
                                                                                 
;Get the average for each of these guys

for m=0, depth-1 do begin
  for n=0, ((2*num-1)*(num))-1 do begin
                    
    finite_val=where(finite(cl_arr(*,n,m)) eq 1)
    if finite_val(0) ne -1 then begin
    Cl_means(n,m)=mean(cl_arr(where(finite(cl_arr(*,n,m)) eq 1),n,m))
    Ct_means(n,m)=mean(ct_arr(where(finite(ct_arr(*,n,m)) eq 1),n,m))
    endif else begin
    Cl_means(n,m)=!values.f_nan
    Ct_means(n,m)=!values.f_nan
    endelse
  
  endfor                                                                               
endfor

sort_index=sort(angles)
                                                                                 
;Sort things by angle
sorted_mags=mags(sort_index)
sorted_angles=angles(sort_index)
sorted_Cl_means=Cl_means(sort_index,*)
sorted_Ct_means=Ct_means(sort_index,*)
                                                                    
uniq_indices=uniq(sorted_angles)

set_plot,'ps'
device,filename=Postscript_dir+'Simulated_covariance_back_2048'+$
		'_X_w='+string(strtrim(xwindspeed_m,2))+$
	        '_Y_w='+string(strtrim(ywindspeed_m,2))+$
		'_OS='+strtrim(outer_scale,2)+$
		'_Depth='+strtrim(depth,2)+'.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
                                                                                 
white='FFFFFF'x
black='000000'x
!P.CHARSIZE=.6
!P.THICK=4.
!P.NOERASE=0
                                                                         
for time=0, depth-1 do begin
        
  color_arr=[0]
  angle_arr=['0']
  sym_arr=[0]
                                                                                 
  ;Set up plot
  !p.multi=[0,1,1]
  plot,[0,34],[-1.5,1.5],/xstyle,/ystyle,/nodata,xr=[0,26],yr=[-0.5,2],background=white,$
        position=[0.10,0.10,0.8,0.9],title='Average Longitudinal Correlatiions for Xwindspeed='+string(strtrim(xwindspeed_m,2))+' m/s Ywindspeed='$
        +string(strtrim(ywindspeed_m,2))+' m/s Outer Scale='+string(strtrim(outer_scale,2))+$
        ' m Depth= '+string(strtrim(time,2)),$
        xtitle='Separation in n',ytitle='Covariance'
                                                                                 
color_index=0

for ind=0, N_elements(uniq_indices)-2 do begin
                                                                                 
  if ind eq 0 then begin
                                                                                 
     temp_mags=sorted_mags(0:uniq_indices(ind))
     temp_angles=sorted_angles(0:uniq_indices(ind))
     temp_cl_means=sorted_cl_means(0:uniq_indices(ind),time)
     temp_sorted_indices=sort(temp_mags)
     oplot,temp_mags,temp_cl_means,color=10*color_index,psym=1
     sym_arr=[sym_arr,1]
     color_arr=[color_arr,(10*color_index)]
     angle_arr=[angle_arr,string(strtrim(180*temp_angles(0)/!pi,2))]
     color_index=color_index+1
                                                                                 
  endif
                                                                                 
  if ((uniq_indices(ind+1)-uniq_indices(ind)) gt 4) then begin
     temp_mags=sorted_mags(uniq_indices(ind)+1:uniq_indices(ind+1))
     temp_angles=sorted_angles(uniq_indices(ind)+1:uniq_indices(ind+1))
     temp_cl_means=sorted_cl_means(uniq_indices(ind)+1:uniq_indices(ind+1),time)
     temp_sorted_indices=sort(temp_mags)
     print,temp_mags(temp_sorted_indices)
     print,temp_angles(temp_sorted_indices)
     print,temp_cl_means(temp_sorted_indices)
     oplot,temp_mags,temp_cl_means,color=10*color_index,psym=1
     sym_arr=[sym_arr,1]
     color_arr=[color_arr,(10*color_index)]
     angle_arr=[angle_arr,string(strtrim(180*temp_angles(0)/!pi,2))]
     color_index=color_index+1
                                                                                 
    endif
endfor

legend,angle_arr,psym=sym_arr,colors=color_arr,position=[.82,.9],/normal

endfor                                              

for time=0, depth-1 do begin
                                   
!p.multi=[0,1,1]
plot,[0,30],[-1.5,1.5],/xstyle,/ystyle,/nodata,xr=[0,26],yr=[-.5,2],background=white,$
        position=[0.10,0.10,0.8,0.9],title='Average Transverse Correlations for Xwindspeed='+string(strtrim(xwindspeed_m,2))+' m/s Ywindspeed='$
        +string(strtrim(ywindspeed_m,2))+' m/s Outer Scale='+string(strtrim(outer_scale,2))+$
        ' m Depth= '+string(strtrim(time,2)),$
        xtitle='Separation in n',ytitle='Covariance'
                                                                                 
color_arr=[0]
angle_arr=['0']
sym_arr=[0]
                                                                                 
color_index=0
                                                                                 
for ind=0, N_elements(uniq_indices)-2 do begin
                                                                                 
  if ind eq 0 then begin
                                                                                 
     temp_mags=sorted_mags(0:uniq_indices(ind))
     temp_angles=sorted_angles(0:uniq_indices(ind))
     temp_ct_means=sorted_ct_means(0:uniq_indices(ind),time)
     temp_sorted_indices=sort(temp_mags)
     oplot,temp_mags,temp_ct_means,color=10*color_index,psym=1
     sym_arr=[sym_arr,1]
     color_arr=[color_arr,(10*color_index)]
     angle_arr=[angle_arr,string(strtrim(180*temp_angles(0)/!pi,2))]
     color_index=color_index+1
                                                                                 
  endif
  if ((uniq_indices(ind+1)-uniq_indices(ind)) gt 4) then begin
     temp_mags=sorted_mags(uniq_indices(ind)+1:uniq_indices(ind+1))
     temp_angles=sorted_angles(uniq_indices(ind)+1:uniq_indices(ind+1))
     temp_ct_means=sorted_ct_means(uniq_indices(ind)+1:uniq_indices(ind+1),time)
     temp_sorted_indices=sort(temp_mags)
     ;print,temp_mags(temp_sorted_indices)
     ;print,temp_angles(temp_sorted_indices)
     ;print,temp_ct_means(temp_sorted_indices)
     oplot,temp_mags,temp_ct_means,color=10*color_index,psym=1
     sym_arr=[sym_arr,1]
     color_arr=[color_arr,(10*color_index)]
     angle_arr=[angle_arr,string(strtrim(180*temp_angles(0)/!pi,2))]
     color_index=color_index+1
                                                                                 
                                                                                 
   endif
endfor
                                                                                 
legend,angle_arr,psym=sym_arr,colors=color_arr,position=[.82,.9],/normal
                                                                                
endfor

device,/close
set_plot,'X'
                                                                                 
stop
                                                                                 
end