;spot_ellipticity

;PURPOSE:
;	This script will analyze the spots in each of the CWFS raw/*.fits files
;	for the [date] directory and try to fit an ellipse to them.  
;
;	An overall net ellipticity will point to a wind effect that is described 
;	in the Image motion paper by Martin
;

pro one_spot_ellipticity

If (not keyword_set(Date)) then Date='10'
If (not keyword_set(stats_directory)) then stats_directory='stats'
If (not keyword_set(HalfSpot)) then HalfSpot = 5

RawDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$
	'/ShackHartman/raw/'
StatsDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$
	'/ShackHartman/stats_1_31/'

;ARRAY OF RAW FILES
RawFiles=file_search(RawDir+'h*.fits',count=NumRawFiles)
StatsFiles=file_search(StatsDir+'*5132*sub_offsets.txt',count=NumStatsFiles)
SpotNo=614
for FileNum = 0, NumStatsFiles -1 do begin

  fits_read,RawFiles(FileNum),Im
  readcol,StatsFiles(FileNum),Index,Xcoord,Xoffset,Ycoord,Yoffset,$
          format='f,f,f,f,f'
  tvscl,congrid(Im(Xcoord(SpotNo)-HalfSpot:Xcoord(SpotNo)+HalfSpot,$
                   Ycoord(SpotNo)-HalfSpot:Ycoord(SpotNo)+HalfSpot),250,250,$
        /center)
  Spot=Im(Xcoord(SpotNo)-HalfSpot:Xcoord(SpotNo)+HalfSpot,$
           Ycoord(SpotNo)-HalfSpot:Ycoord(SpotNo)+HalfSpot)
  Spot=Spot/max(Spot)
  Flux=total(Spot)
  ;XSpot=rebin(findgen(2*HalfSpot+1)+Xcoord(SpotNo)-HalfSpot,2*HalfSpot+1,2*HalfSpot+1)
  ;YSpot=rebin(findgen(2*HalfSpot+1)+Ycoord(SpotNo)-HalfSpot,2*HalfSpot+1,2*HalfSpot+1)
  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-----------
  e1=(I11-I22)/(I11+I22)
  e2=2*I12/(I11+I22)
  e=sqrt(e1^2+e2^2)
  pa=0.5*atan((I12),(I11-I22))
  print,I11/Flux,I22/Flux,I12,e1,e2,e,pa
  stop

endfor

stop
end