pro ellip_offset_ind,Date=Date,stats_directory=stats_directory,$
        x_center_file=x_center_file,y_center_file=y_center_file,$
        depth=depth,x_off=x_off,y_off=y_off


If (not keyword_set(Date)) then Date='10'
If (not keyword_set(stats_directory)) then stats_directory='stats'
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(depth)) then depth=1
If (not keyword_set(postscript_dir)) then postscript_dir='distorted_postscripts'
If (not keyword_set(x_off)) then x_off=0
If (not keyword_set(y_off)) then y_off=0

stats_directory='/a/pippin01/Volumes/u08/lsst/CPanalysis/2005-05-'+Date+$
          '/ShackHartman/'+stats_directory+'/'
postscript_dir='/a/pippin01/Volumes/u08/lsst/CPanalysis/fft_ellip/'


nonfiducial=[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

complement,findgen(625),nonfiducial,fiducial
fitted=fiducial

;ARRAYS OF FILES
offset_files=file_search(stats_directory+'*5132*newest_offsets.txt',count=num_stats_files)

;CONSTANTS
num=25
grid_spots=625
index=findgen(625)

;ARRAY REFERENCE 25x25=625 FOR GRID
subx=2*indgen(num)-25           ;The x-reference numbers
suby=indgen(num)-12             ;The y-reference numbers

;HUGE DATA CUBE
x_offset_hor=fltarr(num_stats_files,grid_spots)
y_offset_hor=fltarr(num_stats_files,grid_spots)
x_offset_ver=fltarr(num_stats_files,grid_spots)
y_offset_ver=fltarr(num_stats_files,grid_spots)

;Correlations
timee1=fltarr(num_stats_files-depth,grid_spots)
timee2=fltarr(num_stats_files-depth,grid_spots)

timee1_arr=fltarr(grid_spots)
timee2_arr=fltarr(grid_spots)

;FILL THE DATA CUBES
for file_number=0,num_stats_files-1 do begin

        ;Take the results from the offset.txt files and fill the cube
        readcol,offset_files(file_number),x_offset,y_offset,format='X,X,f,X,f'

        ;Obtain the difference ordered in COLUMN, ROW between centroid and grid         point
	
	x_offset(where(x_offset eq -10))=0
        y_offset(where(y_offset eq -10))=0

        x_offset_hor(file_number,*)=x_offset
        y_offset_hor(file_number,*)=y_offset

        ;Obtain the difference ordered in ROW, COLUMN
        x_offset_ver(file_number,*)=x_offset((index mod 25)*25+index/25)
        y_offset_ver(file_number,*)=y_offset((index mod 25)*25+index/25)

endfor

readcol,offset_files(0),xcoord,ycoord,format='x,f,x,f,x'
totalxoff=total(x_offset_hor,1)
totalyoff=total(y_offset_hor,1)

alphaxsq=x_offset_hor^2
alphaysq=y_offset_hor^2
alphaxsqmean=avg(alphaxsq,0)
alphaysqmean=avg(alphaysq,0)
psue1=(alphaxsq-alphaysq)/(alphaxsq+alphaysq)
psue2=(2*x_offset_hor*y_offset_hor)/(alphaxsq+alphaysq)
psue1(where(finite(psue1) eq 0))=0
psue2(where(finite(psue2) eq 0))=0

;Compute the value of <alphaxsq>-<alphaysq>/<alphaxsq>+<alphaysq>
e1ensavg=(alphaxsqmean-alphaysqmean)/(alphaxsqmean+alphaysqmean)
e2ensavg=(2*alphaxsqmean*alphaysqmean)/(alphaxsqmean+alphaysqmean)

;Compute the value of <alphaxsq-alphaysq>/<alphasq+alphasq>


;Compute the value of <(alphaxsq-alphaysq)/(alphaxsq+alphaysq)
avgpsue1=avg(psue1,0)
avgpsue2=avg(psue2,0)

e1mean=mean(avgpsue1)
e1var=variance(avgpsue1)
e2mean=mean(avgpsue2)
e2var=variance(avgpsue2)

set_plot,'ps'
device,filename=Postscript_dir+'E1_E2_stick'+Date+'.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.0
!P.THICK=4.
                                                                        
!p.multi=[0,1,1]
                                                                        
partvelvec,avgpsue1(fitted),avgpsue2(fitted),xcoord(fitted),ycoord(fitted),$
        yrange=[100,820],xrange=[180,880],xstyle=17,ystyle=17,$
        position=[.05,.05,.7,.9],color=black,background=white
xyouts,.75,.9,'Mean e1='+strtrim(e1mean,2),/normal
xyouts,.75,.87,'Variance e1='+strtrim(e1var,2),/normal
xyouts,.75,.45,'Mean e2='+strtrim(e2mean,2),/normal
xyouts,.75,.42,'Variance e2='+strtrim(e2var,2),/normal
xyouts,.75,.83,'Max abs(e1)='+strtrim(max(abs(avgpsue1)),2),/normal
xyouts,.75,.39,'Max abs(e2)='+strtrim(max(abs(avgpsue2)),2),/normal

depthmax=4

;Do time delayed psuedoellipticity
for depth=0, depthmax do begin

;Correlations
timee1=fltarr(num_stats_files-depth,grid_spots)
timee2=fltarr(num_stats_files-depth,grid_spots)

timee1num=fltarr(num_stats_files-depth,grid_spots)
timee2num=fltarr(num_stats_files-depth,grid_spots)
timeden=fltarr(num_stats_files-depth,grid_spots)
                                                                 
timee1_arr=fltarr(grid_spots)
timee2_arr=fltarr(grid_spots)

  for file_number=0,num_stats_files-depth-1 do begin
                                                                       
   timee1(file_number,*)=(x_offset_hor(file_number,*)*x_offset_hor(file_number+depth,*)-$
	 y_offset_hor(file_number,*)*y_offset_hor(file_number+depth,*))$
       /(x_offset_hor(file_number,*)*x_offset_hor(file_number+depth,*)+$
         y_offset_hor(file_number,*)*y_offset_hor(file_number+depth,*))
   
   timee2(file_number,*)=(2*x_offset_hor(file_number,*)*y_offset_hor(file_number+depth,*)$
       )/(x_offset_hor(file_number,*)*x_offset_hor(file_number+depth,*)+$
         y_offset_hor(file_number,*)*y_offset_hor(file_number+depth,*))

   timee1num(file_number,*)=x_offset_hor(file_number,*)*x_offset_hor(file_number+depth,*)-$
         y_offset_hor(file_number,*)*y_offset_hor(file_number+depth,*)

   timee2num(file_number,*)=2*x_offset_hor(file_number,*)*y_offset_hor(file_number+depth,*)

   timeden(file_number,*)=x_offset_hor(file_number,*)*x_offset_hor(file_number+depth,*)+$
         y_offset_hor(file_number,*)*y_offset_hor(file_number+depth,*)                                                                        
 endfor

;Get average of numerator and denominator
timee1numavg=avg(timee1num,0)
timee2numavg=avg(timee2num,0)
timedenavg=avg(timeden,0)

;Get total of numeratro and denominator
timee1totnum=total(timee1num,1)
timee2totnum=total(timee2num,1)
timetotden=total(timeden,1)

tote1=timee1totnum/timetotden
tote2=timee2totnum/timetotden

tote1(where(finite(tote1) eq 0))=0
tote2(where(finite(tote2) eq 0))=0

tote1mean=mean(tote1)
tote1var=variance(tote1)
tote2mean=mean(tote2)
tote2var=variance(tote2)

timee1avg=timee1numavg/timedenavg
timee2avg=timee2numavg/timedenavg
timee1avg(where(finite(timee1avg) eq 0))=0
timee2avg(where(finite(timee2avg) eq 0))=0

timee1avgmean=mean(timee1avg)
timee1avgvar=variance(timee1avg)
timee2avgmean=mean(timee2avg)
timee2avgvar=variance(timee2avg)

timee1(where(finite(timee1) eq 0))=0
timee2(where(finite(timee2) eq 0))=0

timee1arr=avg(timee1,0)
timee2arr=avg(timee2,0)

timee1mean=mean(timee1arr)
timee1var=variance(timee1arr)
timee2mean=mean(timee2arr)
timee2var=variance(timee2arr)

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

                                                                        
partvelvec,timee1arr(fitted),timee2arr(fitted),xcoord(fitted),ycoord(fitted),$
        yrange=[100,820],xrange=[180,880],xstyle=17,ystyle=17,$
        position=[.05,.05,.7,.9],color=black,background=white

xyouts,.75,.9,'Mean e1='+strtrim(timee1mean,2),/normal
xyouts,.75,.87,'Variance e1='+strtrim(timee1var,2),/normal
xyouts,.75,.45,'Mean e2='+strtrim(timee2mean,2),/normal
xyouts,.75,.42,'Variance e2='+strtrim(timee2var,2),/normal
xyouts,.75,.83,'Max abs(e1)='+strtrim(max(abs(timee1arr)),2),/normal
xyouts,.75,.39,'Max abs(e2)='+strtrim(max(abs(timee2arr)),2),/normal

!p.charthick=1.2
xyouts,.05,.95,'<(!4a!6!lx!N!6(x,y,t)!4a!6!lx!N!6(x,y,t+'+strtrim(7*depth,2)+$
        ')-!4a!6!ly!N!6(x,y,t)!4a!6!ly!N!6(x,y,t+'+strtrim(7*depth,2)+'))/'$
        +'(!4a!6!lx!N!6(x,y,t)!4a!6!lx!N!6(x,y,t+'+strtrim(7*depth,2)+$
    ')-!4a!6!ly!N!6(x,y,t)!4a!6!ly!N!6(x,y,t+'+strtrim(7*depth,2)+'))>',$
	/normal

partvelvec,timee1avg(fitted),timee2avg(fitted),xcoord(fitted),ycoord(fitted),$
        yrange=[100,820],xrange=[180,880],xstyle=17,ystyle=17,$
        position=[.05,.05,.7,.9],color=black,background=white
                                                                           
                                                                           
xyouts,.75,.9,'Mean e1='+strtrim(timee1avgmean,2),/normal
xyouts,.75,.87,'Variance e1='+strtrim(timee1avgvar,2),/normal
xyouts,.75,.45,'Mean e2='+strtrim(timee2avgmean,2),/normal
xyouts,.75,.42,'Variance e2='+strtrim(timee2avgvar,2),/normal
xyouts,.75,.83,'Max abs(e1)='+strtrim(max(abs(timee1avg)),2),/normal
xyouts,.75,.39,'Max abs(e2)='+strtrim(max(abs(timee2avg)),2),/normal
                                                                           
!p.charthick=1.2
xyouts,.05,.95,'<!4a!6!lx!N!6(x,y,t)!4a!6!lx!N!6(x,y,t+'+strtrim(7*depth,2)+$
        ')-!4a!6!ly!N!6(x,y,t)!4a!6!ly!N!6(x,y,t+'+strtrim(7*depth,2)+')>/'$
        +'<!4a!6!lx!N!6(x,y,t)!4a!6!lx!N!6(x,y,t+'+strtrim(7*depth,2)+$
    ')-!4a!6!ly!N!6(x,y,t)!4a!6!ly!N!6(x,y,t+'+strtrim(7*depth,2)+')>',$
	/normal
                                                                           

partvelvec,tote1(fitted),tote2(fitted),xcoord(fitted),ycoord(fitted),$
        yrange=[100,820],xrange=[180,880],xstyle=17,ystyle=17,$
        position=[.05,.05,.7,.9],color=black,background=white
                                                                               
                                                                               
xyouts,.75,.9,'Mean e1='+strtrim(tote1mean,2),/normal
xyouts,.75,.87,'Variance e1='+strtrim(tote1var,2),/normal
xyouts,.75,.45,'Mean e2='+strtrim(tote2mean,2),/normal
xyouts,.75,.42,'Variance e2='+strtrim(tote2var,2),/normal
xyouts,.75,.83,'Max abs(e1)='+strtrim(max(abs(tote1)),2),/normal
xyouts,.75,.39,'Max abs(e2)='+strtrim(max(abs(tote2)),2),/normal
                                                                               
!p.charthick=1.2
xyouts,.05,.95,'!MI!N(!4a!6!lx!N!6(x,y,t)!4a!6!lx!N!6(x,y,t+'+strtrim(7*depth,2)+$
        ')-!4a!6!ly!N!6(x,y,t)!4a!6!ly!N!6(x,y,t+'+strtrim(7*depth,2)+'))/'$
        +'!MI!N(!4a!6!lx!N!6(x,y,t)!4a!6!lx!N!6(x,y,t+'+strtrim(7*depth,2)+$
    ')-!4a!6!ly!N!6(x,y,t)!4a!6!ly!N!6(x,y,t+'+strtrim(7*depth,2)+'))',$
        /normal
                                                                               

endfor

device,/close
set_plot,'x'

stop

end