;spot_drift
;
;PURPOSE:
;	This script will read in all of the shackhartman data from a series 
;	of images with path given by
;	
;	[Directory]hr5132[letter]_file_number_stats.txt
;
;	and record into a data cube the centroids of the spots
;
;CALLING SEQUENCE
;
;	spot_drift,[Directory=Directory,[letter=letter],[id_array]]]
;
;INPUTS
;
;	Directory	The path to the directory where the Shackhartman data
;			are located
;	letter		The letter included in the files. This ranges from 
;			a-j, corresonding to early-late times
;	id_array	The id numbers (stored in a one-dimensional array)
;			of the points that should be plotted

pro spot_drift, Directory=Directory,letter=letter,id_array=id_array

;set directory path delimiter
    if (!version.os_family eq 'unix') then begin
       delim='/'
    endif else begin
       delim='\'
    endelse

If (not keyword_set(Directory)) then Directory='/a/pippin01/Volumes/u08/lsst/CPanalysis/2005-05-10/ShackHartman/'
If (not keyword_set(id_array)) then id_array=[300,301,302]
If (not keyword_set(letter)) then letter='a'

;create an array that has the write size in 1_dimension
d_size=N_elements(id_array)
x_offset_array=findgen(d_size)
stack_x=findgen(729,99)
stack_y=findgen(729,99)
for file_number=1,99 do begin                                                                              
	if (file_number lt 10) then filename=Directory+'hr5132'+letter+'_'+$
	string(file_number,format='(I1)')+'_stats.txt'
	if file_number ge 10 and file_number lt 100 then filename=$
	Directory+'hr5132'+letter+'_'+string(file_number,format='(I2)')+'_stats.txt'
	if file_number ge 100 and file_number lt 1000 then filename=$
	Directory+'hr5132'+letter+'_'+string(file_number,format='(I3)')+'_stats.txt'

 
	;readcol,filename,id,x_coord,y_coord,x_err,y_err,mag,sharp,roundness,$
	;	x_offset,y_offset,format='X,X,X,X,X,X,X,X,f,f'
	readcol,filename,x_coord,y_coord,x_offset,y_offset,format='X,f,f,X,X,X,X,X,f,f'

	;Obtain sample offset array for plotting
	x_offset_subarray=x_offset(id_array)
	x_offset_array=[[x_offset_array],[x_offset_subarray]]
	
	;create data cube of (x,y) coordinate pairs
	stack_x(*,file_number-1)=x_coord
	stack_y(*,file_number-1)=y_coord

endfor

weighted_x_arr=findgen(729)
weighted_y_arr=findgen(729)

;For each spot, go through and calculate its average (x,y) coordinate
for spot_number=0,728 do begin
	
	;Any spots where the x_coord and y_coord are -1 will not be 
	;included.  valid_point_indices are the array indices where this
	;is the case
	valid_point_indices=where(stack_x(spot_number,*) gt 0)
	
	;For some spots, there will be no valid fits
	if valid_point_indices(0) ne -1 then begin

		weighted_x=mean(stack_x(spot_number,valid_point_indices))
		weighted_y=mean(stack_y(spot_number,valid_point_indices))
	endif else begin	
		weighted_x=0
		weighted_y=0
	endelse
	
	weighted_x_arr(spot_number)=weighted_x
	weighted_y_arr(spot_number)=weighted_y
	
endfor

weighted_x_arr=weighted_x_arr(where(weighted_x_arr gt 0))
weighted_y_arr=weighted_y_arr(where(weighted_y_arr gt 0))

;!P.multi=[0,1,2]
;plot,x_offset_array(0,*)
;oplot,x_offset_array(1,*)
;oplot,x_offset_array(2,*)

range_max=max(weighted_y_arr)
range_min=min(weighted_y_arr)
domain_max=max(weighted_x_arr)
domain_min=min(weighted_x_arr)

set_plot,'ps'
                                                                                
device,filename='hr5132'+letter+'_reference'+'.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=7.0,/inches        ;Width of plot in x
device,yoffset=1.0,/inches      ;Y position of lower left corner

!P.CHARSIZE=2
!P.THICK=10.

plot,[200,900],[100,900],/NODATA,symsize=2.5,title='Plot of ShackHartman Focal Plane',xtitle='X coordinate',ytitle='Y coordinate',background=white,color=000,/xstyle,/ystyle,position=[.1,.1,.9,.9]

oplot,weighted_x_arr,weighted_y_arr,psym=3,symsize=2.5

device,/close

end