pro routine_ellip_wo_denominator,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,exam_date=exam_date


If (not keyword_set(Date)) then Date='10'
If (not keyword_set(exam_date)) then exam_date='1_31'
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(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='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$
          '/ShackHartman/'+stats_directory+'/'
postscript_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/'+exam_date+'_figs/'



nonfiducial=[1,2,3,4,5,6,7,8,9,13,14,17,18,19,20,21,22,23,24,25,$
          26,27,28,29,30,31,32,38,39,45,46,47,48,49,50,$
          51,52,53,54,55,63,64,71,72,73,74,75,$
          76,77,78,79,88,89,98,98,99,100,$
          101,102,103,113,114,123,124,125,$
          126,127,128,138,139,149,150,$
          151,152,163,164,174,175,$
          176,177,188,189,200,$
          201,210,213,214,225,$
          226,238,239,249,$
          251,261,262,263,264,265,266,$
          276,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,$
          326,327,336,337,338,339,340,341,$
          351,361,362,363,364,365,366,$
          376,388,389,$
          401,413,414,425,$
          426,427,438,439,449,450,$
          451,452,454,463,464,473,474,475,$
          476,477,478,488,489,498,499,500,$
          501,502,503,504,513,514,522,523,524,525,$
          526,527,528,529,530,538,539,546,547,548,549,550,$
          551,552,553,554,555,556,563,564,570,571,572,573,574,575,$
          576,577,578,579,580,581,582,588,589,594,595,596,597,598,599,600,$
          601,602,603,604,605,606,607,608,609,610,613,614,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*sub_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+'Psuedo_Time_Ellip_wo_denom'+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]

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

   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,*)=(x_offset_hor(file_number,*)*y_offset_hor(file_number+depth,*)+x_offset_hor(file_number+depth)*y_offset_hor(file_number,*))

   timeden(file_number,*)=$
     sqrt(x_offset_hor(file_number,*)^2+x_offset_hor(file_number+depth,*)^2)*$
     sqrt(y_offset_hor(file_number,*)^2+y_offset_hor(file_number+depth,*)^2)                                                                        
   timee1(file_number,*)=timee1num(file_number,*)/timeden(file_number,*)

   timee2(file_number,*)=timee2num(file_number,*)/timeden(file_number,*)

  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(where(tote1 ne 0)))
tote1var=variance(tote1(where(tote1 ne 0)))
tote2mean=mean(tote2(where(tote2 ne 0)))
tote2var=variance(tote2(where(tote2 ne 0)))

timee1avg=timee1numavg
timee2avg=timee2numavg

timee1avgmean=mean(timee1avg(where(timee1avg ne 0)))
timee1avgvar=variance(timee1avg(where(timee1avg ne 0)))
timee2avgmean=mean(timee2avg(where(timee2avg ne 0)))
timee2avgvar=variance(timee2avg(where(timee2avg ne 0)))

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(where(timee1arr ne 0)))
timee1var=variance(timee1arr(where(timee1arr ne 0)))
timee2mean=mean(timee2arr(where(timee1arr ne 0)))
timee2var=variance(timee2arr(where(timee1arr ne 0)))

!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,.87,'Mean e1='+strtrim(timee1mean,2),/normal
xyouts,.75,.83,'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,.79,'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,.92,'!4e!61=<(!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

xyouts,.05,.97,'!4e!62=<(!4a!6!lx!N!6(x,y,t)'+$
	'!4a!6!ly!N!6(x,y,t+'+strtrim(7*depth,2)+$
        ')+!4a!6!ly!N!6(x,y,t)!4a!6!lx!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,.87,'Mean e1='+strtrim(timee1avgmean,2),/normal
xyouts,.75,.83,'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,.79,'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,.92,'!4e!61=<(!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 

xyouts,.05,.97,'!4e!62=<(!4a!6!lx!N!6(x,y,t)'+$
	'!4a!6!ly!N!6(x,y,t+'+strtrim(7*depth,2)+$
        ')+!4a!6!ly!N!6(x,y,t)!4a!6!lx!N!6(x,y,t+'+strtrim(7*depth,2)+'))>/'$
        ,/normal

endfor

device,/close
set_plot,'x'

stop

end