;xy_offsets
;
;PURPOSE:
;	This script will return the average offset over the whole night for a given
;	point on the ShackHartman grid.  The covariance spectrum will then be obtained.
;

pro new_xy_offsets_average,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


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

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+'/'

;ARRAYS OF FILES
offset_files=file_search(stats_directory+'*5132*offsets.txt',count=num_stats_files)
xcor_name=stats_directory+'average_xcorrelations.txt'
ycor_name=stats_directory+'average_ycorrelations.txt'
off_name=stats_directory+'average_offsets.txt'


get_lun,xcor
openw,xcor,xcor_name
get_lun,ycor
openw,ycor,ycor_name
get_lun,off
openw,off,off_name

;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)
x_coord_arr=fltarr(num_stats_files,grid_spots)
y_coord_arr=fltarr(num_stats_files,grid_spots)
;OFFSET AVERAGES
y_off_avg=fltarr(grid_spots)
x_off_avg=fltarr(grid_spots)
x_coord_avg=fltarr(grid_spots)
y_coord_avg=fltarr(grid_spots)

;FILL THE DATA CUBES WITH THE OFFSETS
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_coord,x_offset,y_coord,y_offset,$
	format='X,f,f,f,f'
	
	;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
	x_coord_arr(file_number,*)=x_coord
	y_coord_arr(file_number,*)=y_coord

endfor
		
;DATA CUBEs are now in place. Now obtain the mean offset for each grid point
for spot_number=0,grid_spots-1 do begin

  x_valid=where(x_offset_hor(*,spot_number) ne -10)
  y_valid=where(y_offset_hor(*,spot_number) ne -10)

  if where(x_valid ne y_valid) ne -1 then stop
  
  if x_valid(0) ne -1 then begin 
  ;obtain average offsets for fiducial points
  x_off_avg(spot_number)=mean(x_offset_hor(x_valid,spot_number))  
  y_off_avg(spot_number)=mean(y_offset_hor(x_valid,spot_number))
  ;obtain average coordinates as well
  x_coord_avg(spot_number)=mean(x_coord_arr(x_valid,spot_number))
  y_coord_avg(spot_number)=mean(y_coord_arr(x_valid,spot_number))
  endif else begin
  x_off_avg(spot_number)=0
  y_off_avg(spot_number)=0
  endelse

  printf,off,FORMAT='(%"%i\t%f\t%f\t%f\t%f")',$
             spot_number,x_coord_avg(spot_number),x_off_avg(spot_number),$
		y_coord_avg(spot_number),y_off_avg(spot_number)


endfor

free_lun,off
close,off

;OFFSETS are arranged by column, then row.  Also obtain arrangement by row, then column
x_off_avg_hor=x_off_avg
y_off_avg_hor=y_off_avg
x_off_avg_ver=x_off_avg((index mod 25)*25+index/25)
y_off_avg_ver=y_off_avg((index mod 25)*25+index/25)




  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_off_avg_hor(spot_index) ne 0) and $
		    (x_off_avg_hor(spot_index+spacing) ne 0) $
		    then begin
		    ;X_off and Y_off contain offsets to tweek the plots
		    Cxl_arr(spot_index)=(x_off_avg_hor(spot_index)+x_off)*$
		    (x_off_avg_hor(spot_index+spacing)+x_off)
		    Cxt_arr(spot_index)=(y_off_avg_hor(spot_index)+y_off)*$
                    (y_off_avg_hor(spot_index+spacing)+y_off)
		    ;print,FORMAT='(%"%i\t%i\t")',spot_index,spot_index+spacing
	
		 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,0),Cxt(spacing,0),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_off_avg_ver(spot_index) ne 0) and $
                   (y_off_avg_ver(spot_index+spacing) ne 0) $
		   then begin

                   Cyl_arr(spot_index)=(y_off_avg_ver(spot_index)+y_off)*$
                   (y_off_avg_ver(spot_index+spacing)+y_off)
                   Cyt_arr(spot_index)=(x_off_avg_ver(spot_index)+x_off)*$
                   (x_off_avg_ver(spot_index+spacing)+x_off)
                   ;print,FORMAT='(%"%i\t%i\t")',spot_index,spot_index+spacing

          	   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,0),Cyt(spacing,0),N_elements(correlated)


     endfor
	
endfor	
	
free_lun,xcor
free_lun,ycor
close,xcor
close,ycor


set_plot,'ps'
device,filename=Postscript_dir+'average_cn.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=.7
!P.THICK=1.2

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
!P.multi=[0,1,1]

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


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


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

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

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

!p.multi=[0,1,1]
 fitted=where(x_off_avg ne 0)
 partvelvec,x_off_avg(fitted),y_off_avg(fitted),$
       x_coord_avg(fitted),y_coord_avg(fitted),$
       yrange=[100,820],xrange=[180,880],$
       background=white,color=black,$
       POSITION=[0.10,0.10,0.9,0.9],$
       xtitle='x coordinate',ytitle='y coordinate',$
       title='Average Offsets'


device,/close
set_plot,'X'


end


