;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 new_xy_offsets_10_cn_depth_ind_letter,Date=Date,stats_directory=stats_directory,$
	x_center_file=x_center_file,y_center_file=y_center_file,$
	depth=depth


If (not keyword_set(Date)) then Date='10'
If (not keyword_set(stats_directory)) then stats_directory='stats_temp'
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=2
If (not keyword_set(postscript_dir)) then postscript_dir='postscripts'

letter=['hr5132a','hr5132b','hr5132c','hr5132d','hr5132e','hr5132f','hr5132g','hr5132h','hr5132i']

stats_directory='/a/pippin01/Volumes/u08/lsst/CPanalysis/2005-05-'+Date+$
	  '/ShackHartman/'+stats_directory+'/'
postscript_dir='/a/pippin01/Volumes/u08/lsst/CPanalysis/2005-05-'+Date+$
          '/ShackHartman/'+postscript_dir+'/'


;RUN THIS SCRIPT FOR EACH 100 FILES
for j=0, N_elements(letter)-1 do begin
                                                                                                                                                   
offset_files=file_search(stats_directory+letter(j)+'*offsets.txt',count=num_stats_files)


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

;ARRAY REFERENCE 25x25=625 FOR GRID
subx=2*indgen(num)-25           ;The x-reference numbers
suby=indgen(num)-12             ;The y-reference numbers

;HUGE DATA CUBE
x_offset_hor=fltarr(num_stats_files,grid_spots)
y_offset_hor=fltarr(num_stats_files,grid_spots)
x_offset_ver=fltarr(num_stats_files,grid_spots)
y_offset_ver=fltarr(num_stats_files,grid_spots)

;Correlations
Cxt_full=fltarr(num_stats_files,num,depth)
Cxl_full=fltarr(num_stats_files,num,depth)
Cyt_full=fltarr(num_stats_files,num,depth)
Cyl_full=fltarr(num_stats_files,num,depth)


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

	;Take the results from the offset.txt files and fill the cube
	readcol,offset_files(file_number),x_offset,y_offset,format='X,X,f,X,f'
	

	not_fitted=where(x_offset eq -1,complement=fitted)
	
	;Obtain the difference ordered in COLUMN, ROW between centroid and grid 	point
	x_offset_hor(file_number,*)=x_offset
	y_offset_hor(file_number,*)=y_offset

	;Obtain the difference ordered in ROW, COLUMN
	x_offset_ver(file_number,*)=x_offset((index mod 25)*25+index/25)
	y_offset_ver(file_number,*)=y_offset((index mod 25)*25+index/25)
	
   endfor
	
   ;DATA CUBEs are now in place. Begin the covariance

   ;Loop through all of the files in the 
   for file_number=0,num_stats_files-depth do begin

	;CORRELATIONS
	Cxl=fltarr(num,depth)
	Cxt=fltarr(num,depth)
	Cyl=fltarr(num,depth)
	Cyt=fltarr(num,depth)

	for time=0, depth-1 do begin

		for spacing=0, num-1 do begin
	
			Cxl_arr=fltarr(num*num)
			Cxt_arr=fltarr(num*num)
			Cyl_arr=fltarr(num*num)
			Cyt_arr=fltarr(num*num)

			;SCAN FIRST ALONG THE ROWS
	
			for column_index=0, num-1 do begin
				for row_index=0, num -1 do begin
				spot_index=column_index*num+row_index
			
				;If we have reached the edge, start on the next row
				If (row_index+spacing) ge num then break
				;Check to see if both have valid values
				If (x_offset_hor(file_number,spot_index) ne -1) and $
				   (x_offset_hor(file_number+time,spot_index+spacing) ne -1) then begin

					Cxl_arr(spot_index)=x_offset_hor(file_number,spot_index)*$
					  x_offset_hor(file_number+time,spot_index+spacing)
					Cxt_arr(spot_index)=y_offset_hor(file_number,spot_index)*$
                                          y_offset_hor(file_number+time,spot_index+spacing)
				;print,FORMAT='(%"%i\t%i\t")',spot_index,spot_index+spacing
	
			;	C_num=C_num+1
				endif

				endfor
			endfor
		
			correlated=where(Cxl_arr ne 0)
			If correlated(0) ne -1 then begin
				Cxl(spacing,time)=mean(Cxl_arr(correlated))
				Cxt(spacing,time)=mean(Cxt_arr(correlated))
			endif else begin
				Cxl(spacing,time)=0
				Cxt(spacing,time)=0
			endelse
		
			;printf,xcor,FORMAT='(%"%f\t%f\t%i")',Cxl(spacing-1),Cxt(spacing-1),N_elements(correlated)
		
	
			;SCAN SECOND ALONG Y

               		 for column_index=0, num-1 do begin
                       		 for row_index=0, num -1 do begin
                       			 spot_index=column_index*num+row_index
	
        	               		 ;If we have reached the edge, start on the next row
                		        If (row_index+spacing) ge num then break
                     			;Check to see if both have valid values
                       			If (y_offset_ver(file_number,spot_index) ne -1) and $
                           		(y_offset_ver(file_number+time,spot_index+spacing) ne -1) then begin

                                		Cyl_arr(spot_index)=y_offset_ver(file_number,spot_index)*$
                                          	y_offset_ver(file_number+time,spot_index+spacing)
                                		Cyt_arr(spot_index)=x_offset_ver(file_number,spot_index)*$
                                          	x_offset_ver(file_number+time,spot_index+spacing)
                        ;print,FORMAT='(%"%i\t%i\t")',spot_index,spot_index+spacing

         ;               C_num=C_num+1
          	                         endif

                       		 endfor
               		 endfor

               		 correlated=where(Cyl_arr ne 0)
              		 If correlated(0) ne -1 then begin
                       		 Cyl(spacing,time)=mean(Cyl_arr(correlated))
                        	 Cyt(spacing,time)=mean(Cyt_arr(correlated))
                         endif else begin
                        	 Cyl(spacing,time)=0
                       		 Cyt(spacing,time)=0
                	endelse

                ;printf,ycor,FORMAT='(%"%f\t%f\t%i")',Cyl(spacing-1),Cyt(spacing-1),N_elements(correlated)

     	   endfor
	
	endfor	
	;partvelvec,x_offset(fitted),y_offset(fitted),gridx(fitted),$
	;	gridy(fitted),yrange=[150,800],xrange=[200,850]

	;close,xcor
	;free_lun,xcor		
	;close,ycor
	;free_lun,ycor
	
	Cxt_full(file_number,*,*)=Cxt
	Cxl_full(file_number,*,*)=Cxl
	Cyt_full(file_number,*,*)=Cyt
	Cyl_full(file_number,*,*)=Cyl


   endfor	
		

Cxt_means=fltarr(num,depth)
Cxl_means=fltarr(num,depth)
Cyt_means=fltarr(num,depth)
Cyl_means=fltarr(num,depth)

for n=0, num-1 do begin
	for m=0, depth-1 do begin

	Cxt_means(n,m)=mean(Cxt_full(*,n,m))
	Cxl_means(n,m)=mean(Cxl_full(*,n,m))
	Cyl_means(n,m)=mean(Cyl_full(*,n,m))
	Cyt_means(n,m)=mean(Cyt_full(*,n,m))

	endfor
endfor



set_plot,'ps'
device,filename=Postscript_dir+letter(j)+'cn_time'+string(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=1.2
!P.THICK=4.

!p.multi=[0,1,1]

items_time=['t=0']
sym_time=[1]
items_k_tran=['Data Points','n!e-7/32!n','n!e-13/64!n','n!e-27/128!n']
sym_k_tran=[1,0,0,0]
line_k_tran=[0,1,2,3]
items_k_lon=['Data Points','n!e-1/3!n']
sym_k_long=[0,1]
colors_time=[0]

;*************FIRST PLOT
plot,Cxt_means(*,0),psym=1,title="Transverse correlation along X direction",$
               xtitle="n",ytitle="Covariance",$
               background=white,color=black,$
               yrange=[min(Cxt_means(*,0))-.5,max(Cxt_means(*,0))+.5],$
               POSITION=[0.10,0.10,0.9,0.9]

;Over Plot the various separations in time and also generate legends

for i=1,depth-1 do begin

	oplot,Cxt_means(*,i),psym=(((i+1) mod 5)+1),color=20*((i+1)/5)
	sym_time=[sym_time,((i+1) mod 5)+1]
	items_time=[items_time,string(strtrim(i*7,2))+' seconds']
	colors_time=[colors_time,(20*((i+1)/5))]

endfor

legend,items_time,psym=sym_time,position=[.7,.8],colors=colors_time,/normal

;*************SECOND PLOT
plot,(Cxt_means(*,0)-Cxt_means(*,depth-1)),psym=1,$
	       title="Transverse correlation along X direction for t=0-t="+$
	       strtrim(string(depth*7),2)+' seconds',$
               xtitle="n",ytitle="Covariance",$
               background=white,color=black,$
               yrange=[0,2],$
               POSITION=[0.10,0.10,0.9,0.9]

oplot,findgen(num-1)+1,(findgen(num-1)+1)^(-7./32.),psym=0,color=black,$
        linestyle=1
oplot,findgen(num-1)+1,(findgen(num-1)+1)^(-13./64.),psym=0,color=black,$
        linestyle=2
oplot,findgen(num-1)+1,(findgen(num-1)+1)^(-27./128.),psym=0,color=black,$
	linestyle=3

legend,items_k_tran,linestyle=line_k_tran,position=[.6,.8],/normal

;**************THIRD PLOT
!p.multi=[0,5,5]

for n=0, num-1 do begin

	Cxt_max=floor(max(Cxt_full(*,n,0)-Cxt_full(*,n,depth-1)))+.5
	Cxt_min=floor(min(Cxt_full(*,n,0)-Cxt_full(*,n,depth-1)))-.5
	n_bins=40
	Cxt_hist=histogram((Cxt_full(*,n,0)-Cxt_means(n,depth-1)),nbins=n_bins,$
	max=Cxt_max,min=Cxt_min,locations=Cxt_locations)

	plot,Cxt_locations,Cxt_hist,psym=10,title='Spacing= '+string(strtrim(n,2)),$
	xrange=[Cxt_min,Cxt_max],xtitle="Correlation",ytitle="Number",$
	background=white,color=black
        
endfor

!p.multi=[0,1,1]

;***************FOURTH PLOT
plot,Cxl_means(*,0),psym=1,title="Longitudinal correlation along X direction",$
               xtitle="n",ytitle="Covariance (Pixels Squared)",$
               background=white,color=black,$
               yrange=[min(Cxl_means(*,0))-.5,max(Cxl_means(*,0))+.5],$
               POSITION=[0.10,0.10,0.9,0.9]

for i=1,depth-1 do begin
	
        oplot,Cxl_means(*,i),psym=(((i+1) mod 5)+1),color=20*((i+1)/5)
	
endfor

legend,items_time,psym=sym_time,position=[.2,.4],/normal,colors=colors_time

;**************FIFTH PLOT
plot,(Cxl_means(*,0)-Cxl_means(*,depth-1)),psym=1,$
               title="Longitudinal correlation along X direction for t=0-t="+$
               strtrim(string(depth*7),2)+' seconds',$
               xtitle="n",ytitle="Covariance",$
               background=white,color=black,$
               yrange=[0,max(Cxl_means(*,0)-Cxl_means(*,depth-1))+.5],$
               POSITION=[0.10,0.10,0.9,0.9]

oplot,findgen(num-1)+1,(findgen(num-1)+1)^(-1./3.),psym=0,color=black

legend,items_k_lon,linestyle=line_k_lon,position=[.6,.8],/normal

;**************SIXTH PLOT
!p.multi=[0,5,5]
                                                                                                                                                         
for n=0, num-1 do begin
                                                                                                                                                         
        Cxl_max=floor(max(Cxl_full(*,n,0)-Cxl_full(*,n,depth-1)))+.5
        Cxl_min=floor(min(Cxl_full(*,n,0)-Cxl_full(*,n,depth-1)))-.5
        n_bins=40
        Cxl_hist=histogram((Cxl_full(*,n,0)-Cxl_full(*,n,depth-1)),nbins=n_bins,$
        max=Cxl_max,min=Cxl_min,locations=Cxl_locations)
                                                                                                                                                         
        plot,Cxl_locations,Cxl_hist,psym=10,title='Spacing= '+string(strtrim(n,2)),$
        xrange=[Cxl_min,Cxl_max],xtitle="Correlation",ytitle="Number",$
        background=white,color=black
                                                                                                                                                         
endfor

!P.multi=[0,1,1]
;**************SEVENTH PLOT
plot,Cyt_means(*,0),psym=1,title="Transverse correlation along Y direction",$
               xtitle="n",ytitle="Covariance (Pixels Squared)",$
               background=white,color=black,$
               yrange=[min(Cyt_means(*,0))-.5,max(Cyt_means(*,0))+.5],$
               POSITION=[0.10,0.10,0.9,0.9]

for i=1,depth-1 do begin

        oplot,Cyt_means(*,i),psym=(((i+1) mod 5)+1),color=20*((i+1)/5)

endfor

legend,items_time,psym=sym_time,position=[.6,.8],/normal,colors=colors_time

;***************EIGHT PLOT
plot,(Cyt_means(*,0)-Cyt_means(*,depth-1)),psym=1,$
               title="Transverse correlation along Y direction for t=0-t="+$
               strtrim(string(depth*7),2)+' seconds',$
               xtitle="n",ytitle="Covariance",$
               background=white,color=black,$
               yrange=[min(Cyt_means(*,0)-Cyt_means(*,depth-1))-.5,$
                       max(Cyt_means(*,0)-Cyt_means(*,depth-1))+.5],$
               POSITION=[0.10,0.10,0.9,0.9]

oplot,findgen(num-1)+1,(findgen(num-1)+1)^(-7./32.),psym=0,color=black,$
        linestyle=1
oplot,findgen(num-1)+1,(findgen(num-1)+1)^(-13./64.),psym=0,color=black,$
        linestyle=2
oplot,findgen(num-1)+1,(findgen(num-1)+1)^(-27./128.),psym=0,color=black,$
        linestyle=3                                                                        
legend,items_k_tran,linestyle=line_k_tran,position=[.6,.8],/normal

;***************NINTH PLOT
!p.multi=[0,5,5]
                                                                                                                                                         
for n=0, num-1 do begin
                                                                                                                                                         
        Cyt_max=floor(max(Cyt_full(*,n,0)-Cyt_full(*,n,depth-1)))+.5
        Cyt_min=floor(min(Cyt_full(*,n,0)-Cyt_full(*,n,depth-1)))-.5
        n_bins=40
        Cyt_hist=histogram((Cyt_full(*,n,0)-Cyt_full(*,n,depth-1)),nbins=n_bins,$
        max=Cyt_max,min=Cyt_min,locations=Cyt_locations)
                                                                                                                                                         
        plot,Cyt_locations,Cyt_hist,psym=10,title='Spacing= '+string(strtrim(n,2)),$
        xrange=[Cyt_min,Cyt_max],xtitle="Correlation",ytitle="Number",$
        background=white,color=black
                                                                                                                                                         
endfor

;***************TENTH PLOT
!p.multi=[0,1,1]

plot,Cyl_means(*,0),psym=1,title="Longitudinal correlation along Y direction",$
               xtitle="n",ytitle="Covariance (Pixels Squared)",$
               background=white,color=black,$
               yrange=[min(Cyl_means(*,0))-.5,max(Cyl_means(*,0))+.5],$
               POSITION=[0.10,0.10,0.9,0.9]

for i=1,depth-1 do begin

        oplot,Cyl_means(*,i),psym=(((i+1) mod 5)+1),color=20*((i+1)/5)

endfor

legend,items_time,psym=sym_time,position=[.6,.8],/normal,colors=colors_time

;***************ELEVENTH PLOT
plot,(Cyl_means(*,0)-Cyl_means(*,depth-1)),psym=1,$
               title="Longitudinal correlation along Y direction for t=0-t="+$
               strtrim(string(depth*7),2)+' seconds',$
               xtitle="n",ytitle="Covariance",$
               background=white,color=black,$
               yrange=[0,max(Cyl_means(*,0)-Cyl_means(*,depth-1))+.5],$
               POSITION=[0.10,0.10,0.9,0.9]

oplot,findgen(num-1)+1,(findgen(num-1)+1)^(-1./3.),psym=0,color=black

legend,items_k_lon,linestyle=line_k_lon,position=[.6,.8],/normal

;***************TWELTH PLOT
!p.multi=[0,5,5]
                                                                                                                                                         
for n=0, num-1 do begin
                                                                                                                                                         
        Cyl_max=floor(max(Cyl_full(*,n,0)-Cyl_full(*,n,depth-1)))+.5
        Cyl_min=floor(min(Cyl_full(*,n,0)-Cyl_full(*,n,depth-1)))-.5
        n_bins=40
        Cyl_hist=histogram((Cyl_full(*,n,0)-Cyl_full(*,n,depth-1)),nbins=n_bins,$
        max=Cyl_max,min=Cyl_min,locations=Cyl_locations)
                                                                                                                                                         
        plot,Cyl_locations,Cyl_hist,psym=10,title='Spacing= '+string(strtrim(n,2)),$
        xrange=[Cyl_min,Cyl_max],xtitle="Correlation",ytitle="Number",$
        background=white,color=black
                                                                                                                                                         
endfor

device,/close
set_plot,'X'


endfor

end

