;PRO TESTCALIMAGE
;
;PURPOSE:
;	This script deals with the ShackHartmann test-calibration image.  It
;	is called TestCalImage.gs6 and opens with DS9.
;
;	The optimal coordinates of the grid and their offsets with respect to 
;	a 25-pixel-spaced uniform grid centered at the same center are printed
;	to a file
;
;CALLING SEQUENCE:
;	TestCalImage,[[Image_Dir],Stats_Dir]
;
;KEYWORDS:
;	Image_Dir: The directory in which TestCalImage.gs6 is stored
;	Stats_Dir: The directory where the coordinates and offset will be 
;		   written
;	DATE:      The string '10','11',or '12' corresponding to
;                  2005-05-10,2005-05-11, and 2005-05-12
;
;OUTPUS:
;	[stats_dir]/Calibrated_grid.txt
;

pro testcalimage, Date=Date, Image_Dir=Image_Dir, Stats_Dir=Stats_Dir

If (not keyword_set(Date)) then Date='11'
If (not keyword_set(x_center)) then x_center=[542]
If (not keyword_set(y_center)) then y_center=[465]
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'                                                           


get_lun,stat
openw,stat,stat_file


num=25
subx=2*indgen(num)-25           ;The x-reference numbers
suby=indgen(num)-12             ;The y-reference numbers
box=11                          ;The distance used for finding spot

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

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

; XCENTER, YCENTER, XERR, YERR, MAG
readcol,mag_file,x_coord,y_coord,x_err,y_err,$
	format='f,f,f,f'
                                                                               
;Read in the data columns from the output IRAF coo file.  They are
                ; SHARP, ROUNDNESS
readcol,coo_file,sharp,roundness,$
	format='X,X,X,X,f',SKIPLINE=41
                                                                               
;Each of these vectors will have as many elements as the number of points
num_of_iraf_points=N_elements(x_coord)

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")',$
           ID_index,x_coord(x_index),y_coord(y_index),$
           x_ref_col(row_num),y_ref_col(col_num),$
           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)
           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)
                                                                                                                   
valid_x_coords=x_coords(point_indices)
valid_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,valid_x_coords,xerrs,x_center,$
      /noderivative,function_name='linear_function_x',yerror=xerror)
y_fit=curvefit(fully,valid_y_coords,yerrs,y_center,$
      /noderivative,function_name='linear_function_y',yerror=yerror)
                                                                                                                   
print,FORMAT='(%"%f\t%f")',x_center,y_center


;THE CENTER HAS BEEN FOUND, REESTABLISH THE GRID BASED UPON IT
x_coord_cent=fltarr(num*num)
y_coord_cent=fltarr(num*num)
x_off_grid=fltarr(num*num)
y_off_grid=fltarr(num*num)
new_x_ref_col=12.5*subx+x_center
new_y_ref_col=25.0*suby+y_center

for col_num=0,N_elements(new_y_ref_col)-1 do begin
    for row_num=0,N_elements(new_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) eq -1 then begin
          x_coord_cent(ID_index)=new_x_ref_col(row_num)
	  y_coord_cent(ID_index)=new_y_ref_col(col_num)
          x_off_grid(ID_index)=x_coords(id_index)-new_x_ref_col(row_num)
	  y_off_grid(ID_index)=y_coords(id_index)-new_y_ref_col(col_num)
      endif else begin
	  x_coord_cent(ID_index)=-1
	  y_coord_cent(ID_index)=-1
      endelse

    endfor

endfor

partvelvec,x_off_grid(point_indices),y_off_grid(point_indices),$
	x_coord_cent(point_indices),y_coord_cent(point_indices),$
	yrange=[150,800],xrange=[200,850]



stop

end



                                                                                
