;xy_offsets
;
;PURPOSE:
;	This script will return the difference in coordinates between a 
;	ShackHartman image centroid and the corresponding point on the 
;	reference grid.  I.e. it will give the offsets between the 
;	Image centroid and the supposed true center point of a ShackHartman
;	spot
;	
;	It will also calculate the magnitude and angle of the deflection 
;	vector. 
;
;OUTPUT:
;  	The program will write to the file
;	[stats_directory]/[raw_file_name].offsets.txt
;
;	The following columns will be written
;
;	index|xcoord|ycoord|xoffset|yoffset|magnitude|angle


pro tv_spot_ellipticity,Date=Date,stats_directory=stats_directory,$
	x_center_file=x_center_file,y_center_file=y_center_file,$	
	exam_date=exam_date,HalfSpot=HalfSpot
If (not keyword_set(exam_date)) then exam_Date='1_31'
If (not keyword_set(Date)) then Date='10'
If (not keyword_set(stats_directory)) then stats_directory='stats_'+exam_date+'/'
If (not keyword_set(x_center_file)) then x_center_file='newest_x_centers.txt'
If (not keyword_set(y_center_file)) then y_center_file='newest_y_centers.txt'
If (not keyword_set(HalfSpot)) then HalfSpot = 5


;-----------------------------------------STATS AND RAW .fits IMAGES
stats_directory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$
	  '/ShackHartman/'+stats_directory+'/'
RawDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$
        '/ShackHartman/raw/'
Stats_files=file_search(stats_Directory+'*5132*_sub_offsets.txt',count=num_of_stats_files)
RawFiles=file_search(RawDir+'h*.fits',count=NumRawFiles)

;OBTAIN THE FITTED CENTERS FOR ALL IMAGES
readcol,stats_Directory+x_center_file,x_fit_center,format='X,X,f'
readcol,stats_Directory+y_center_file,y_fit_center,format='X,X,f'

;-----------------------------------------------------------CONSTANTS
num_of_files=N_elements(x_fit_center)
num=25					;
subx=2*indgen(num)-25           	;The x-reference numbers
suby=indgen(num)-12             	;The y-reference numbers

print,num_of_stats_files
print,num_of_files

!p.charsize=3
!p.charthick=2
;***********************************************************************
;FOURTH LOOP
;Run through again 1) subtract offsets
;		   2) Get Spot Characteristics
for file_number=2,num_of_files-1 do begin

        readcol,Stats_files(file_number),grid_index,XCoord,YCoord,$
                                XErr, YErr, bolo_mags, Sharp, $
        		        format='i,f,X,f,X,f,f,f,f'

        stats_path_parts=strsplit(stats_files(file_number),'/',/extract)
        stats_file_ref=strsplit(stats_path_parts(N_elements(stats_path_parts)-1)$
        ,'stats.txt',/extract,/regex)
        stats_file_ref=stats_file_ref(0)

        ;Take the centroid values from the stats.txt files
        stats_file_name=stats_files(file_number)
        offsets_sub_file_name=stats_directory+stats_file_ref+'sub_offsets.txt'
	fits_read,RawFiles(file_number),Im

	get_lun,offsets
        openw,offsets,offsets_sub_file_name

	Spots=where(XCoord ne -10,complement=NonSpots)

        for SpotNo=0, num*num-1 do begin

	  if XCoord(SpotNo) ne -10 then begin
	    ;USE FITS FILE TO GET SPOTS
	    Spot=Im(floor(Xcoord(SpotNo))-HalfSpot:$
		    floor(Xcoord(SpotNo))+HalfSpot,$
                    floor(Ycoord(SpotNo))-HalfSpot:$
		    floor(Ycoord(SpotNo))+HalfSpot)
            window,xsize=512,ysize=600
            !p.noerase=1
	    whiteback=fltarr(88,512)
            whiteback(*)=1000
            tv,congrid(whiteback,512,88),0,513
            tvscl,congrid(Spot,512,512,/center)
	    Spot=Spot/max(Spot)
	    Flux=total(Spot)
	    XSpot=rebin(findgen(2*HalfSpot+1)-HalfSpot,$
              2*HalfSpot+1,2*HalfSpot+1)
  	    YSpot=rebin(transpose(findgen(2*HalfSpot+1)-HalfSpot),$
              2*HalfSpot+1,2*HalfSpot+1)
	    ;QUADROPOLE MOMENTS---------------
	    I11=total(XSpot*XSpot*Spot)
	    I22=total(YSpot*YSpot*Spot)
	    I12=total(XSpot*YSpot*Spot)
	    ;ELLIPTICITY COMPONENTS-----------
	    d=sqrt(.5*(I11+I22)/Flux)
	    e1=(I11-I22)/(I11+I22)
	    e2=2*I12/(I11+I22)
	    e=sqrt(e1^2+e2^2)
	    pa=0.5*atan((I12),(I11-I22))
            plot,[e,e],[pa-!pi,pa],/polar,color='00FF00'x,$
	         xrange=[-.2,.2],yrange=[-.2,.2],$
	         xstyle=5,ystyle=5,$
		 position=[0,0,1,.8]
            xyouts,.03,.95,'e='+strtrim(string(FORMAT='(%"%8.5f")',e),2),/norm,$
		 color='000000'x
	    xyouts,.03,.875,'Angle='+$
		 strtrim(string(FORMAT='(%"%10.2f")',180*pa/!pi),2),/norm,$
	         color='000000'x
	    xyouts,.5,.95,'e1='+$
	         strtrim(string(FORMAT='(%"%8.5f")',e1),2),/norm,$
                 color='000000'x
            xyouts,.5,.875,'e2='+$
		 strtrim(string(FORMAT='(%"%8.5f")',e2),2),/norm,$
                 color='000000'x

	    stop
	    I11=I11/flux
	    I12=I12/flux
	    I22=I22/flux
	  endif else begin
            I11=-10
            I22=-10
            I12=-10
	    d=-10
            e1=-10
            e2=-10
            e=-10
            pa=-10
	  endelse
 
        endfor


endfor

stop

end