;xy_offsets
;
;PURPOSE:
;	This script will return the difference in coordinates between a 
;	ShackHartman image centroid and the corresponding point on the 
;	reference grid
;
;	Since the center of the grid seems to be drifting 

pro cn_ind_new,Date=Date,stats_directory=stats_directory,$
	x_center_file=x_center_file,y_center_file=y_center_file,$
	depth=depth,x_off=x_off,y_off=y_off,num_stats_files=num_stats_files


If (not keyword_set(Date)) then Date='10'
If (not keyword_set(stats_directory)) then stats_directory='stats'
If (not keyword_set(x_center_file)) then x_center_file='x_centers.txt'
If (not keyword_set(y_center_file)) then y_center_file='y_centers.txt'
If (not keyword_set(depth)) then depth=1
If (not keyword_set(postscript_dir)) then postscript_dir='postscripts'
If (not keyword_set(x_off)) then x_off=0
If (not keyword_set(y_off)) then y_off=0

if Date eq '10' then begin
  letter=['hr5132a','hr5132b','hr5132c','hr5132d','hr5132e','hr5132f',$
          'hr5132g','hr5132h','hr5132i']
endif else if Date eq '11' then begin
  letter=['h5132g','h5132h','h5132i','hr5132a','hr5132b','hr5132c','hr5132d',$
          'hr5132e','hr5132f','hr5132g','hr5132h']
endif else if Date eq '12' then begin
  letter=['hr5132a','hr5132b','hr5132c','hr5132d','hr5132e','hr5132f',$
           'hr5132g','hr5132h','hr5132i']
endif

stats_dir='/a/pippin01/Volumes/u08/lsst/CPanalysis/2005-05-'+Date+$
          '/ShackHartman/stats/'
coo_dir='/a/pippin01/Volumes/u08/lsst/CPanalysis/2005-05-'+Date+$
	  '/ShackHartman/coo/'
postscript_dir='/a/pippin01/Volumes/u08/lsst/CPanalysis/2005-05-'+Date+$
          '/ShackHartman/'+postscript_dir+'/'
Cor_directory='/a/pippin01/Volumes/u08/lsst/CPanalysis/2005-05-'+Date+$
          '/ShackHartman/correlations/'

readcol,coo_dir+'fiducial_reference.txt',spots,format='x,x,f'
spot_arr=intarr(25,25)
spot_arr(spots)=130
complement,findgen(625),spots,non_spots


set_plot,'ps'
device,filename=Postscript_dir+'all_ind_both_angles_half_'+date+'.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.

for j=0, N_elements(letter)-1 do begin

;ARRAYS OF FILES-Averages have already been subtracted
offset_files=file_search(stats_dir+letter(j)+'*offsets.txt',$
	count=num_off_files)
;If (not keyword_set(num_stats_files)) then num_stats_files=num_off_files
num_stats_files=num_off_files

;CONSTANTS
num=25
grid_spots=625
index=findgen(625)

;DATA CUBE CONTAINING ALL ANGLES AND MAGNITUDES
phi_arr=fltarr(num_stats_files,grid_spots)
magnitude_arr=fltarr(num_stats_files,grid_spots)

;FILL THE DATA CUBES
for file_number=0, num_stats_files-depth do begin

	;Take the magnitude and angles from the offset.txt files and fill the cube
	readcol,offset_files(file_number),magnitude,phi,format='X,X,X,X,X,f,f'
	magnitude(non_spots)=-10
	phi(non_spots)=-10
	magnitude_arr(file_number,*)=magnitude
	phi_arr(file_number,*)=phi

endfor

;DATA CUBEs are now in place. Begin the covariance
;CORRELATION ARRAYS
Cl_arr=fltarr(num_stats_files,(2*num-1)*(num),depth)
Ct_arr=fltarr(num_stats_files,(2*num-1)*(num),depth)
mags=fltarr((2*num-1)*(num))
angles=fltarr((2*num-1)*(num))

;Loop through all of the files in the stack
for file_number=0, num_stats_files-depth do begin
  print,file_number,offset_files(file_number)

  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)
        
	;Get the grid angles and spacings
	 if file_number eq 0 then begin
		angles(spacing_index)=grid_phi
		mags(spacing_index)=grid_sep
	 endif
         
	;Cl=fltarr(grid_spots,depth)
        Ct=fltarr(grid_spots-abs(delta_x)*abs(delta_y))
	Cl=fltarr(grid_spots-abs(delta_x)*abs(delta_y))

	; Cl=fltarr(grid_spots,depth)
	;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(file_number*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((file_number+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((file_number+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_arr(first_spot) ne -10) and $
                   (phi_arr(second_spot) ne -10))$
                then begin
                ;Longitudinal will be cosine and transverse will be sine
                                                                                
                 Cl(c_index)=magnitude_arr(first_spot)*$
                         cos(phi_arr(first_spot)-grid_phi)*$
                         magnitude_arr(second_spot)*$
                         cos(phi_arr(second_spot)-grid_phi)
                 Ct(c_index)=magnitude_arr(first_spot)*$
                         sin(phi_arr(first_spot)-grid_phi)*$
                         magnitude_arr(second_spot)*$
                         sin(phi_arr(second_spot)-grid_phi)
                c_index=c_index+1
                                                                                
                endif
             endif
           endif                                                                                 
          endfor
                                                                                
                                                                                
        endfor
	
	 ;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(file_number,spacing_index,time)=$
	    mean(Cl(correlated))
	    Ct_arr(file_number,spacing_index,time)=$
	    mean(Ct(correlated))
	
	endif else begin
	 ;If there were no pairs, assign a system no value
	    Cl_arr(file_number,spacing_index,time)=!values.f_nan
	    Ct_arr(file_number,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

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

;Print this to a file
get_lun,cor_angles
openw,cor_angles,stats_directory+'cor_angles'+letter(j)+'.txt'

for i=0, N_elements(sort_index)-1 do begin

printf,cor_angles,FORMAT='(%"%f\t%f\t%f\t%f\t%f")',sorted_mags(i),sorted_angles(i),$
	sorted_Cl_means(i),sorted_Ct_means(i)

endfor

close,cor_angles
free_lun,cor_angles


color_arr=[0]
angle_arr=['0']
sym_arr=[0]

;Set up plot 
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=Letter(j)+'Average Longitudinal Correlations for 5/'+$
	Date+'/2005',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))
     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))
     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



!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=Letter(j)+':Average Transverse Correlations for 5/'+$
        Date+'/2005',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))
     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))
     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


