Pro get_reference_grid,Date=Date,stats_dir=stats_dir,image_dir=image_dir 

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

If (not keyword_set(Date)) then Date='11'
If (not keyword_set(stats_dir)) then stats_dir='/a/pippin01/Volumes/u08/lsst/CPanalysis/2005-05-'+date+'/ShackHartman/stats/'
If (not keyword_set(image_dir)) then image_dir='/a/pippin01/Volumes/u08/lsst/CPanalysis/ShackHartmann_Cal/'
                                                                                
                                                                                
raw_file=image_dir+'TestCalImage.fits'
coo_file=image_dir+'TestCalImage.fits.coo.1'
mag_file=image_dir+'TestCalImage.mag.1.txt'
center_file=image_dir+'TestCalImage_centers.txt'
stat_file=image_dir+'TestCalImage_stats.txt' 
offset_file=stats_dir+'average_offsets.txt'

index=indgen(625)

readcol,stat_file,xcoords,ycoords,format='x,f,f'
readcol,offset_file,xavgcoord,yavgcoord,format='X,f,X,f,X'

xcoord_x=xcoords
ycoord_x=ycoords

cum_invalid=where(ycoord_x eq -1)
cum_valid=where(ycoord_x ne -1)

;Get the Reference Grid Center
ycoord_arr=fltarr(25,25)
ycoord_arr(index)=ycoords
avg_col=fltarr(25)
for i=0, 24 do begin
   valid=where(ycoord_arr(*,i) ne -1)
   if valid(0) ne -1 then begin
	avg_col(i)=avg(ycoord_arr(where(ycoord_arr(*,i) ne -1),i))
   endif
endfor


xcoord_y=xcoords((index mod 25)*25+index/25)
ycoord_y=ycoords((index mod 25)*25+index/25)
xcoord_arr=fltarr(25,25)
xcoord_arr(index)=xcoord_y
avg_row=fltarr(25)
for i=0, 24 do begin
   valid=where(xcoord_arr(*,i) ne -1)
   if valid(0) ne -1 then begin     
    avg_row(i)=avg(xcoord_arr(where(xcoord_arr(*,i) ne -1),i))
   endif
endfor

x_cal_center=avg(avg_row)
y_cal_center=avg(avg_col(where(avg_col ne 0)))

;Get the Night's Center
ycoord_arr=fltarr(25,25)
ycoord_arr(index)=yavgcoord
avg_col=fltarr(25)
for i=0, 24 do begin
   valid=where(ycoord_arr(*,i) ne 0)
   if valid(0) ne -1 then begin
        avg_col(i)=avg(ycoord_arr(where(ycoord_arr(*,i) ne 0),i))
   endif
endfor
                                                                                               
                                                                                               
xcoord_y=xavgcoord((index mod 25)*25+index/25)

xcoord_arr=fltarr(25,25)
xcoord_arr(index)=xcoord_y
avg_row=fltarr(25)
for i=0, 24 do begin
   valid=where(xcoord_arr(*,i) ne 0)
   if valid(0) ne -1 then begin
    avg_row(i)=avg(xcoord_arr(where(xcoord_arr(*,i) ne 0),i))
   endif
endfor
                                                                                               
x_data_center=avg(avg_row)
y_data_center=avg(avg_col(where(avg_col ne 0)))

x_data_off_center=xavgcoord-x_data_center
y_data_off_center=yavgcoord-y_data_center
x_data_off_center(cum_invalid)=-1
y_data_off_center(cum_invalid)=-1

xcoord_dif=xcoord_x(1:624)-xcoord_x(0:623)
ycoord_dif=ycoord_y(1:624)-ycoord_y(0:623)

get_lun,delta_file
openw,delta_file,stats_dir+'Delta_offsets.txt'

for i=0, 624 do begin

  printf,delta_file,format='(%"%i\t%f\t%f\t%f\t%f")',$
         i,xavgcoord(i),x_data_off_center(i),yavgcoord(i),y_data_off_center(i)

endfor

close,delta_file
free_lun,delta_file


!p.multi=[0,1,2]
loadct,39

ycoord_dif_max=28
ycoord_dif_min=22
n_bins=40
ycoord_dif_hist=histogram($
	ycoord_dif(where((ycoord_dif gt 10) and (ycoord_dif lt 40))),$
	nbins=n_bins,max=ycoord_dif_max,min=ycoord_dif_min,$
	locations=ycoord_dif_locations)
plot,ycoord_dif_locations,ycoord_dif_hist,psym=10

xcoord_dif_max=28
xcoord_dif_min=22
n_bins=40
xcoord_dif_hist=histogram($ 
        xcoord_dif(where((xcoord_dif gt 10) and (xcoord_dif lt 40))),$
        nbins=n_bins,max=xcoord_dif_max,min=xcoord_dif_min,$ 
        locations=xcoord_dif_locations) 
oplot,xcoord_dif_locations,xcoord_dif_hist,psym=10

plot,xcoord_x(where(xcoord_x gt 0)),ycoord_x(where(ycoord_x gt 0)),psym=3
oplot,xavgcoord(cum_valid),yavgcoord(cum_valid),color=200,psym=3
stop

end