;find_center
;
;PURPOSE:
;	find_center will read in the listed coordinates from the files 
;	of type [filename]_mag.txt and try to find the center of the 
; 	coordinate grid     

pro find_all_center_1letter_diff_weights, Directory=Directory, y_center=y_center, x_center=x_center, $
		 letter=letter,file_number=file_number

If (not keyword_set(Directory)) then Directory='/a/pippin01/Volumes/u08/lsst/CPanalysis/2005-05-10/ShackHartman/'
If (not keyword_set(letter)) then letter='a'
If (not keyword_set(file_number)) then file_number='20'
If (not keyword_set(x_center)) then x_center=[542.5]
If (not keyword_set(y_center)) then y_center=[466]


;Loop through 100 images and obtain the fitted center for each
num=25
num_of_files=100
X_centers=fltarr(num_of_files-1)
Y_centers=fltarr(num_of_files-1)
my_X_centers=fltarr(num_of_files-1)
my_Y_centers=fltarr(num_of_files-1)

for file_number=1,num_of_files-1 do begin

	if (file_number lt 10) then filename=Directory+'stats/hr5132'+letter+'_'+$
        string(file_number,format='(I1)')+'_stats.txt'
        if file_number ge 10 and file_number lt 100 then filename=$
        Directory+'stats/hr5132'+letter+'_'+string(file_number,format='(I2)')+'_stats.txt'
        if file_number ge 100 and file_number lt 1000 then filename=$
        Directory+'stats/hr5132'+letter+'_'+string(file_number,format='(I3)')+'_stats.txt'


	;The x coordinates are symmetric about a point that is not a ShackHartman image
	;The y coordinates are symmetric about a ShackHartman image
	subx=2*indgen(num)-25
	suby=indgen(num)-12

	index=findgen(num*num)

	fullx=fltarr(num*num)
	fully=fltarr(num*num)

	;y is constant for 25 entries; x increases for 25 entries and repeats
	;i.e., scan is along x followed by y
	fullx(index)=subx(index mod num)
	fully=rebin(suby,num*num)

	readcol,filename,x_coord,y_coord,xweights,yweights,format='X,f,f,f,f'

	point_indices=where(x_coord ne -1)

	x_coord=x_coord(point_indices)
	y_coord=y_coord(point_indices)
	xweights=xweights(point_indices)
	yweights=yweights(point_indices)
	fullx=fullx(point_indices)
	fully=fully(point_indices)

	;For some reason, some xerrors and yerrors are zero
	xzeroes=where(xweights eq 0,xcount)
	yzeroes=where(yweights eq 0,ycount)
        If xcount ne 0 then xweights(xzeroes)=0.001
	If ycount ne 0 then yweights(yzeroes)=0.001

	x_fit=curvefit(fullx,x_coord,weights,x_center,/noderivative,function_name='linear_function_x',yerror=xerror)
	y_fit=curvefit(fully,y_coord,weights,y_center,/noderivative,function_name='linear_function_y',yerror=yerror)

	my_X_centers(file_number-1)=total((x_coord-12.5*fullx)/xweights)$
					/total(1/xweights)
	my_Y_centers(file_number-1)=total((y_coord-25*fully)/yweights)$
					/total(1/yweights)
	X_centers(file_number-1)=x_center
	Y_centers(file_number-1)=y_center
	
	print,FORMAT='(%"%s\t%d\t%f\t%f\t%f\t%f")',letter,file_number,x_center,$
		my_X_centers(file_number-1),y_center,my_Y_centers(file_number-1)

	;Errors at each point
	y_err=y_coord-y_fit
	x_err=x_coord-x_fit

endfor

set_plot,'ps'
                                                                                
device,filename=Directory+'postscripts/hr5132'+letter+'center.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]
plot,findgen(99)+1,y_centers,psym=1,title="Y Center Coordinate for files hr5132 "+letter,$
       xtitle="File Number",ytitle="Y Coordinate",background=white,color=black,$
       yrange=[min(y_centers)-10,max(y_centers)+10],POSITION=[0.10,0.10,0.9,0.9]
oplot,findgen(99)+1,my_y_centers,psym=6,color=22
legend,['Manual Fit','IDL Fit'],psym=[6,1],position=[.2,.8],/norm

!P.multi=[0,1,1]
plot,findgen(99)+1,x_centers,psym=1,title="X Center Coordinate for files hr5132 "+letter,$
       xtitle="File Number",ytitle="X Coordinate",background=white,color=black,$
       yrange=[min(x_centers)-10,max(x_centers)+10],POSITION=[0.10,0.10,0.9,0.9]
oplot,findgen(99)+1,my_x_centers,psym=6,color=22
legend,['Manual Fit','IDL Fit'],psym=[6,1],position=[.2,.8],/norm

                                                              
device,/close

set_plot,'X'


end