;pro coo_mag_to_stats_fid
;
;PURPOSE:
;	This script will take the input from the .mag.1 files returned by DIGIPHOT>DAOPHOT>PHOT
;	(and processed with testscript.cl to extract the necessary columns) and the .coo.1 files
;	returned by DIGIPHOT>DAOPHOT>DAOFIND and fit the fitted coordinate points to a grid.
;
;	The grid is set up to have [num]x[num] points.  The center of the initial grid is eyeballed
;	by examining hr5132a_1_raw.fits and the spacing is assumed to be 600 microns/24 microns per 
;	pixel=25 pixels.  A fit is performed on this grid, and a new center is
;	determined.  This center is then fed to the next iteration and a new grid is layed out. 
;	The process repeats through all [num_of_files]*[letters] files.
;
;	In this way, the grid should shift so that two centroids are not fit to the same reference
;	point.
;
;	The file will output statistics from the .mag.1 and .coo.1 files in [stats_directory]
;
;CALLING SEQUENCE:
;	new_coo_mag_to_stat_fid,[Date=Date,[Directory=Directory,[stats_directory=stats_directory,
;		[xcenter=xcenter,[ycenter=ycenter,[center_file=center_file]]]]]]]]
;
;INPUTS
;       DATE: the string '10','11',or '12' corresponding to 
;       2005-05-10,2005-05-11, and 2005-05-12
;       
;       LETTER: an array of strings corresponding to the letters of the 
;       files you want to analyze
;       
;       STATS_DIRECTORY: the directory to which the statistics will be output
;
;       XCENTER: a scalar indicating the approximate x coordinate for the 
;       center of the grid for the first image 
;
;       YCENTER: a scalar indicating the approximate y coordinate for the 
;       center of the grid for the first image
;
;       CENTER_FILE: The name of the file to which the fitted centers will
;       be written
;
;OUTPUTS
;
;       FILES:
;
;       <Directory/[stats_directory]/hr5132[letter]_[number_stats].txt>
;               ---this file contains the following entries
;               ID_index=A number 1-625 indicating where on the grid the point is located
;               x_coord=The x coordinate of the fitted centroid
;               y_coord=The y coordinate of the fitted centroid
;               x_error=The error from the IRAF fit
;               y_err=The error from the IRAF fit
;               mag=The magnitude assessed by IRAF
;               sharp=IRAF parameter specificying sharpness
;               roundness=IRAF parameter specifiying how "round" the fit was 
;               x_ref_dif,y_ref_dif
;
;       <Directory/[stats_directory]/centers.txt
;
;               ---this file contains the fitted center coordinates
;               x_center=the x coordinate of the center
;               y_center=the y coordinate of the center
;  
;	
;
;
pro new_coo_mag_to_stats_fid, Date=Date,Directory=Directory,$
	stats_directory=stats_directory,xcenter=xcenter,ycenter=ycenter,$
	center_file=center_file
	
If (not keyword_set(Date)) then Date='11'
If (not keyword_set(Directory)) then Directory='/a/pippin01/Volumes/u08/lsst/CPanalysis/2005-05-'$
	+date+'/ShackHartman/'
If (not keyword_set(stats_directory)) then stats_directory='stats_temp'
If (not keyword_set(x_center)) then x_center=[542]
If (not keyword_set(y_center)) then y_center=[465]
If (not keyword_set(center_file)) then center_file=Directory+stats_directory+'/centers.txt'

raw_files=FILE_SEARCH(Directory+'raw/h*.fits',count=num_raw_files)
coo_files=FILE_SEARCH(Directory+'coo/h*.coo.1',count=num_coo_files)
mag_files=FILE_SEARCH(Directory+'mag/h*.mag.1.txt',count=num_mag_files)


;FILE UNITS
get_lun,cent
openw,cent,center_file
printf,cent,FORMAT='(%"%s\t%s\t%s\t%s")','Y CENTER','X CENTER','Y ERROR','X ERROR'

;CONSTANTS
num=25				;The number of rows and columns in reference grid
num_of_files=100		;The number of files for each letter
subx=2*indgen(num)-25		;The x-reference numbers
suby=indgen(num)-12		;The y-reference numbers
box=11				;The distance used for assigning coordinates to the ref grid

;FIDUCIAL-These points will not be included in the analysis

fiducial=[1,2,3,4,5,6,7,8,18,19,20,21,22,23,24,25,$
          26,27,28,29,30,31,45,46,47,48,49,50,$
          51,52,53,54,55,72,73,74,75,$
          76,77,78,79,98,99,100,$
          101,102,103,123,124,125,$
          126,127,149,150,$
          151,152,174,175,$
          176,200,$
          201,225,$
          238,239,$
          261,262,263,264,265,266,$
          286,287,288,289,290,291,$
          301,302,303,304,305,306,307,308,309,$
          310,311,312,313,314,315,316,$
          317,318,319,320,321,322,323,324,325,$
          336,337,338,339,340,341,$
          361,362,363,364,365,366,$
          388,389,$
          401,425,$
          426,427,449,450,$
          451,452,452,473,474,475,$
          476,477,478,498,499,500,$
          501,502,503,504,522,523,524,525,$
          526,527,528,529,530,546,547,548,549,550,$
          551,552,553,554,555,556,570,571,572,573,574,575,$
          576,577,578,579,580,581,582,594,595,596,597,598,599,600,$
          601,602,603,604,605,606,607,608,609,617,618,619,620,621,$
          622,623,624,625]-1


;SUB ARRAY ALLOCATION
X_coords=fltarr(num*num)
Y_coords=fltarr(num*num)
Mags=fltarr(num*num)
Sharps=fltarr(num*num)
Xerrs=fltarr(num*num)
Yerrs=fltarr(num*num)

;FULL ARRAY ALLOCATION
X_centers=fltarr(num_raw_files)				;Computed centers
Y_centers=fltarr(num_raw_files)				;Computed centers
Point_indices_arr=intarr(num_raw_files)			;Number of points matched to grid 
Num_of_iraf_points_arr=intarr(num_raw_files)		;Number of points found by iraf
mag_arr=[0]						;Array to hold average magnitudes
sharp_arr=[0]						;Array to hold average sharpness


;MAIN LOOP
for file_number=0, num_raw_files-1 do begin
                                 
		;Check to see if the raw, mag, and coo files match up
		raw_path_parts=strsplit(raw_files(file_number),'/',/extract)
		raw_file_ref=strsplit(raw_path_parts(N_elements(raw_path_parts)-1),'raw.fits',$
		/extract,/regex)
		mag_path_parts=strsplit(mag_files(file_number),'/',/extract)
                mag_file_ref=strsplit(mag_path_parts(N_elements(mag_path_parts)-1),'raw.fits.mag.1.txt',$
		/extract,/regex)
		coo_path_parts=strsplit(coo_files(file_number),'/',/extract)
                coo_file_ref=strsplit(coo_path_parts(N_elements(coo_path_parts)-1),'raw.fits.coo.1',$
                /extract,/regex)
		stats_file=Directory+stats_directory+'/'+raw_file_ref+'stats.txt'

		if (mag_file_ref ne coo_file_ref) or (raw_file_ref ne coo_file_ref) then stop 

		;Figure out the new reference grid, applied around x_center and y_center
		x_ref_col=12.5*subx+x_center[0]
		y_ref_col=25.0*suby+y_center[0]
		
		;prepare arrays to be used in the linear fit
		x_coords=fltarr(num*num)
		y_coords=fltarr(num*num)
		Mags=fltarr(num*num)
                Sharps=fltarr(num*num)
                Xerrs=fltarr(num*num)
                Yerrs=fltarr(num*num)
		x_coords(*)=-1
		y_coords(*)=-1
		Mags(*)=-1
		Sharps(*)=-1
		Yerrs(*)=-1
		Xerrs(*)=-1

		index=findgen(num*num)
		fullx=fltarr(num*num)
	        fully=fltarr(num*num)
        	fullx(index)=subx(index mod num)
        	fully=rebin(suby,num*num)

		print,file_number
		;Read in the data columns from the output IRAF mag file.  They are
		; XCENTER, YCENTER, XERR, YERR, MAG
		readcol,mag_files(file_number),x_coord,y_coord,x_err,y_err,mag,format='f,f,f,f,f'

		;Read in the data columns from the output IRAF coo file.  They are 
		; SHARP, ROUNDNESS
		readcol,coo_files(file_number),sharp,roundness,format='X,X,X,f,f',SKIPLINE=41
	
		;Each of these vectors will have as many elements as the number of points
		;that IRAF fitted
		num_of_iraf_points=N_elements(x_coord)
		
		get_lun,stat
		openw,stat,stats_file
		;This loop will go through all of the possible places where a point could be
		;and check to see if there is an image point near it
		for col_num=0,N_elements(y_ref_col)-1 do begin
			for row_num=0,N_elements(x_ref_col)-1 do begin

				;This index will provide the location on the grid
				ID_index=col_num*N_elements(x_ref_col)+row_num

				;Check to see if this point is in the fiducial region
				check=where(fiducial eq ID_index)
				If check(0) ne -1 then begin
					found_point=0
				endif else begin

				;Pick points where x might be near and those where y might 
				;be near
				x_index_pos=where(abs(x_coord-x_ref_col(row_num)) lt box)
				y_index_pos=where(abs(y_coord-y_ref_col(col_num)) lt box)
				
				;Do some indexing gymnastics to determine if they match
				for y_sub_index=0,N_elements(y_index_pos)-1 do begin
					x_sub_index=where(x_index_pos eq $
							  y_index_pos(y_sub_index))
					
					;The grid coordinates have been matched to a spot
					if ((x_sub_index(0) ne -1) and (y_index_pos(0) ne -1)$
					and (x_index_pos(0) ne -1)) then begin
						x_index=x_index_pos(x_sub_index)
						y_index=y_index_pos(y_sub_index)
						x_ref_dif=x_coord(x_index)-x_ref_col(row_num)
						y_ref_dif=y_coord(y_index)-y_ref_col(col_num)
						
						;Check to see if it is 
						;Print this to the stats file
						
						printf,stat,FORMAT=$
		'(%"%3.0d\t%0.3f\t%0.3f\t%0.3f\t%0.3f\t%0.3f\t%0.3f\t%0.3f\t%0.3f\t%0.3f")',$
						ID_index,x_coord(x_index),y_coord(y_index),$
						x_err(x_index),y_err(y_index),mag(x_index),$
						sharp(x_index),roundness(x_index),$
						x_ref_dif,y_ref_dif
						;found_point=1
						
						;Add it to the full grid array
						X_coords(ID_index)=x_coord(x_index)
						Y_coords(ID_index)=y_coord(y_index)
						Mags(ID_index)=mag(x_index)
						Sharps(ID_index)=sharp(x_index)
						Xerrs(ID_index)=x_err(x_index)
						Yerrs(ID_index)=y_err(x_index)
						found_point=1
						break
					endif else begin
						found_point=0
					endelse
				
				endfor
				
				endelse

				;The grid coordinates had no spot within a box		
				if found_point eq 0 then begin
				printf,stat,$
				FORMAT='(%"%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d")',$
				ID_index,-1,-1,-1,-1,-1,-1,-1,-1,-1
				X_coords(ID_index)=-1
                                Y_coords(ID_index)=-1
                                Mags(ID_index)=-1
                                Sharps(ID_index)=-1
                                Xerrs(ID_index)=-1
                                Yerrs(ID_index)=-1
				endif
			endfor
		endfor 
		close,stat
		free_lun,stat

		;Now we have an array that mimics the contents of the _stats.txt files
		;Fill the Arrays with the processed grids, removing the -1s
		point_indices=where(x_coords ne -1)
                                                                                
	        x_coords=x_coords(point_indices)
       		y_coords=y_coords(point_indices)
       		xerrs=xerrs(point_indices)
        	yerrs=yerrs(point_indices)
        	fullx=fullx(point_indices)
        	fully=fully(point_indices)
			
		;For some reason, some xerrors and yerrors are zero
	        xzeroes=where(xerrs eq 0,xcount)
       		yzeroes=where(yerrs eq 0,ycount)
        	If xcount ne 0 then xerrs(xzeroes)=0.001
       		If ycount ne 0 then yerrs(yzeroes)=0.001

		x_fit=curvefit(fullx,x_coords,xerrs,x_center,$
		      /noderivative,function_name='linear_function_x',yerror=xerror)
        	y_fit=curvefit(fully,y_coords,yerrs,y_center,$
		      /noderivative,function_name='linear_function_y',yerror=yerror)
                                                                                
		;Fill the Cumulative Arrays
		X_centers(file_number)=x_center
	        Y_centers(file_number)=y_center
       		Point_indices_arr(file_number)=N_elements(point_indices)
        	Num_of_iraf_points_arr(file_number)=num_of_iraf_points

		printf,cent,FORMAT='(%"%f\t%f\t%f\t%f\t%s")',$
		x_center,y_center,xerror,yerror,raw_file_ref

		print,FORMAT='(%"%s\t\t\%f\t%f")',mag_file_ref,x_center,y_center
	

endfor

close,cent
free_lun,cent

end