;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_all_letter_ws, Directory=Directory, y_center=y_center,$
	x_center=x_center, stats_directory=stats_directory,$
	letter=letter,file_number=file_number

If (not keyword_set(Directory)) then Directory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-10/ShackHartman/'
If (not keyword_set(letter)) then letter=['a','b','c','d','e','f','g','h','i']
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]
If (not keyword_set(stats_directory)) then stats_directory='stats'

;Loop through 100 images and obtain the fitted center for each
num=25
num_of_files=100
;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)


;Allocate Arrays to hold data
letters=N_elements(letter)
X_centers_ws=fltarr(letters*(num_of_files-1))
Y_centers_ws=fltarr(letters*(num_of_files-1))
X_centers_wos=fltarr(letters*(num_of_files-1))
Y_centers_wos=fltarr(letters*(num_of_files-1))
Point_indices_arr_ws=intarr(letters*(num_of_files-1)) 
Point_indices_arr_wos=intarr(letters*(num_of_files-1))
Num_of_iraf_points_arr=intarr(letters*(num_of_files-1))
mag_arr_ws=[0]
sharp_arr_ws=[0]
mag_arr_wos=[0]
sharp_arr_wos=[0]

for letter_ind=0,letters-1 do begin

	for file_number=1,num_of_files-1 do begin

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

	;allocate space for full set of basis numbers
	fullx=fltarr(num*num)
	fully=fltarr(num*num)

	;For fits
	x_center_ws=[542.5]
	x_center_wos=[542.5]
	y_center_ws=[466]
	y_center_wos=[466]

	;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)

	;Take in the number of IRAF points, and the x_coordinates and y_coordinates
	readcol,filename,num_of_iraf_points,format='d',numline=1
	readcol,filename,x_coord,y_coord,mag,sharp,format='X,f,f,X,X,f,f'
		
	point_indices=where(x_coord ne -1)
	
	;Temporary Vectors for spots with spider included
	x_coord_ws=x_coord(point_indices)
	y_coord_ws=y_coord(point_indices)
	fullx_ws=fullx(point_indices)
	fully_ws=fully(point_indices)
	
	;Temporary Vectors for spots without spider
	x_coord_wos=[x_coord(0:299),x_coord(325:624)]
	y_coord_wos=[y_coord(0:299),y_coord(325:624)]
	fullx_wos=[fullx(0:299),fullx(325:624)]
	fully_wos=[fully(0:299),fully(325:624)]
	
	point_indices_wos=where(x_coord_wos ne -1)
	fullx_wos=fullx_wos(point_indices_wos)
	fully_wos=fully_wos(point_indices_wos)
	x_coord_wos=x_coord_wos(point_indices_wos)
	y_coord_wos=y_coord_wos(point_indices_wos)
	
	mag_arr_ws=[mag_arr_ws,mean(mag(point_indices))]
	sharp_arr_ws=[sharp_arr_ws,mean(sharp(point_indices))]
	
	mag_arr_wos=[mag_arr_wos,mean(mag(point_indices_wos))]
	sharp_arr_wos=[sharp_arr_wos,mean(sharp(point_indices_wos))]

	x_fit_ws=curvefit(fullx_ws,x_coord_ws,weights,x_center_ws,/noderivative,function_name='linear_function_x',yerror=xerror)
	y_fit_ws=curvefit(fully_ws,y_coord_ws,weights,y_center_ws,/noderivative,function_name='linear_function_y',yerror=yerror)

	x_fit_wos=curvefit(fullx_wos,x_coord_wos,weights,x_center_wos,/noderivative,function_name='linear_function_x',yerror=xerror)
        y_fit_wos=curvefit(fully_wos,y_coord_wos,weights,y_center_wos,/noderivative,function_name='linear_function_y',yerror=yerror)

	;Cumulative Vectors
	X_centers_ws(letter_ind*(num_of_files-1)+file_number-1)=x_center_ws
	Y_centers_ws(letter_ind*(num_of_files-1)+file_number-1)=y_center_ws
	X_centers_wos(letter_ind*(num_of_files-1)+file_number-1)=x_center_wos
        Y_centers_wos(letter_ind*(num_of_files-1)+file_number-1)=y_center_wos

	Point_indices_arr_ws(letter_ind*(num_of_files-1)+file_number-1)=N_elements(point_indices)
	Point_indices_arr_wos(letter_ind*(num_of_files-1)+file_number-1)=N_elements(point_indices_wos)
	Num_of_iraf_points_arr(letter_ind*(num_of_files-1)+file_number-1)=$
					num_of_iraf_points
	
	;Output to show processing
	print,FORMAT='(%"%s\t%d\t\t%f\t\t%f")',letter(letter_ind),file_number,x_center_wos,y_center_wos

	endfor

endfor

;These arrays contain zero as their initial element.  Get rid of it
mag_arr_ws=mag_arr_ws(where(mag_arr_ws ne 0))
sharp_arr_ws=sharp_arr_ws(where(sharp_arr_ws ne 0))
mag_arr_wos=mag_arr_wos(where(mag_arr_wos ne 0))
sharp_arr_wos=sharp_arr_wos(where(sharp_arr_wos ne 0))

set_plot,'ps'
                                                                                
device,filename=Directory+'postscripts/hr5132_unweighted_'+stats_directory+'_'+letter(0)+'_'+letter(letters-1)+'_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.NOERASE=0 
;!P.POSITION=[0.10,0.10,0.9,0.9]

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

for letter_ind=0,letters-1 do begin
	
	plot,findgen(99)+1,$
	y_centers_ws(letter_ind*(num_of_files-1):(letter_ind+1)*(num_of_files-1)-1),$
	psym=1,title="Y Center Coordinate for files hr5132"+letter(letter_ind),$
        xtitle="File Number",ytitle="Y Coordinate",background=white,color=black, $
        yrange=[min(y_centers_ws)-10,max(y_centers_ws)+10]

	oplot,findgen(99)+1,y_centers_wos(letter_ind*(num_of_files-1):(letter_ind+1)*(num_of_files-1)-1),psym=6
                                    
endfor

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

for letter_ind=0,letters-1 do begin

        plot,findgen(99)+1,$
        x_centers_ws(letter_ind*(num_of_files-1):(letter_ind+1)*(num_of_files-1)-1),$
        psym=1,title="X Center Coordinate for files hr5132"+letter(letter_ind),$
        xtitle="File Number",ytitle="X Coordinate",background=white,color=black, $
        yrange=[min(x_centers_ws)-10,max(x_centers_ws)+10]
                   
	oplot,findgen(99)+1,x_centers_wos(letter_ind*(num_of_files-1):(letter_ind+1)*(num_of_files-1)-1),psym=6                 
endfor

!p.multi=[0,1,1]
plot,findgen(letters*(num_of_files-1))+1,y_centers_ws,psym=1,title="Y Center Coordinate for files hr5132 "+letter(0)+'-'+letter(letters-1),$
        xtitle="File Number",ytitle="Y Coordinate",background=white,color=black, $
	yrange=[min(y_centers_ws)-10,max(y_centers_ws)+10],$
	POSITION=[0.10,0.10,0.9,0.9]

oplot,findgen(letters*(num_of_files-1))+1,y_centers_wos,psym=6

!p.multi=[0,1,1]
plot,findgen(letters*(num_of_files-1))+1,x_centers_ws,psym=1,title="X Center Coordinate for files hr5132 "+letter(0)+'-'+letter(letters-1),$
        xtitle="File Number",ytitle="X Coordinate",background=white,color=black, $
        yrange=[min(x_centers_ws)-10,max(x_centers_ws)+10],$
	POSITION=[0.10,0.10,0.9,0.9]
                
oplot,findgen(letters*(num_of_files-1))+1,x_centers_wos,psym=6                                                     

!p.multi=[0,1,1]
plot,findgen(letters*(num_of_files-1))+1,mag_arr_ws,psym=1,title="Average Magnitude for files hr5132 "+letter(0)+'-'+letter(letters-1),$
        xtitle="File Number",ytitle="Average Magnitude",background=white,$
	color=black,yrange=[min(mag_arr_ws)-.2,max(mag_arr_ws)+.2],$
        POSITION=[0.10,0.10,0.9,0.9]

oplot,findgen(letters*(num_of_files-1))+1,mag_arr_wos,psym=6

!p.multi=[0,1,1]
plot,findgen(letters*(num_of_files-1))+1,sharp_arr_ws,psym=1,$
	title="Average Sharpness for files hr5132 "+letter(0)+'-'+letter(letters-1),$
        xtitle="File Number",ytitle="Average Sharpness",background=white,$
        color=black,yrange=[min(sharp_arr_ws)-.2,max(sharp_arr_ws)+.2],$
        POSITION=[0.10,0.10,0.9,0.9]

oplot,findgen(letters*(num_of_files-1))+1,sharp_arr_wos,psym=6

!p.multi=[0,1,1]
plot,findgen(letters*(num_of_files-1))+1,Point_indices_arr_ws,psym=1,$
        title="Number of Points Fitted for files hr5132 "+letter(0)+$
	'-'+letter(letters-1),$
        xtitle="File Number",ytitle="Number of Points",background=white,$
        color=black,yrange=[min(Point_indices_arr_ws)-10,max(Point_indices_arr_ws)+10],$
        POSITION=[0.10,0.10,0.9,0.9]

oplot,findgen(letters*(num_of_files-1))+1,Point_indices_arr_wos,psym=1
!p.multi=[0,1,1]
plot,findgen(letters*(num_of_files-1))+1,Num_of_iraf_points_arr,psym=1,$
        title="Number of Points Found by IRAF for files hr5132 "+letter(0)+$
        '-'+letter(letters-1),$
        xtitle="File Number",ytitle="Number of Points",background=white,$
        color=black,yrange=[min(Point_indices_arr_ws)-10,max(Point_indices_arr_ws)+10],$
        POSITION=[0.10,0.10,0.9,0.9]

device,/close

set_plot,'X'

end