pro routine_fft_ellip_gs,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=0
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

;Indices for Rearranging FFT
a=findgen(num/2+1)
b=-reverse(findgen(num/2)+1)
ns=[a,b]
subn=sort(ns)
n_sort=ns(subn)
ms=ns
subm=subn
m_sort=n_sort

;Low pass indices
pos_ind=findgen(25,25)
neg_ind=624-findgen(25,25)

low_pos_ind=pos_ind(1:5,1:5)
low_neg_ind=neg_ind(0:4,0:4)


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

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

timee1num=fltarr(num_stats_files-depth,num,num)
timee2num=fltarr(num_stats_files-depth,num,num)
timeden=fltarr(num_stats_files-depth,num,num)

timee1_arr=fltarr(num,num)
timee2_arr=fltarr(num,num)

frame_e1=fltarr(num_stats_files)
frame_e2=fltarr(num_stats_files)
;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'


  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,*,*))

   timee1temp=timee1num(file_number,*,*)
   timee2temp=timee2num(file_number,*,*)

   frame_e1(file_number)=mean(timee1temp(where(timee1temp ne 0)))
   frame_e2(file_number)=mean(timee2temp(where(timee2temp ne 0)))


endfor

no_bins=8
mean_xfft_abs=complexarr(no_bins,num,num)
mean_yfft_abs=complexarr(no_bins,num,num)
mean_xfft_pha=complexarr(no_bins,num,num)
mean_yfft_pha=complexarr(no_bins,num,num)

if date eq '11' or date eq '10' then begin

e1max=2.4
e1min=-2.4

endif else if date eq '12' then begin

e1max=1.8
e1min=-1.8

endif

e1hist=histogram(frame_e1,Nbins=no_bins,locations=loc_e1,$
		 reverse_indices=e1_ind,$
		 max=e1max,min=e1min)

binsize=(e1max-e1min)/(no_bins-1)

set_plot,'ps'
device,filename=Postscript_dir+'Ellip_vs_Filenumber_gs_'+Date+'.ps'                  
loadct,0                      ;load color table 39
device,/color                   ;allow color on the postscript
device,ysize=8.0,/inches        ;Height of plot in y
device,xsize=12.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,2]

print,e1hist

for i=0, no_bins-1 do begin
	
	rel_indices=e1_ind(e1_ind(i):e1_ind(i+1)-1)
	xfft_array=complexarr(N_elements(rel_indices),num,num)
	yfft_array=complexarr(N_elements(rel_indices),num,num)

	for j=0, N_elements(rel_indices)-1 do begin
		xfft_array(j,*,*)=fft(x_offset_hor(rel_indices(j),*,*))	
		yfft_array(j,*,*)=fft(y_offset_hor(rel_indices(j),*,*))
		
	endfor
	mean_xfft_abs(i,*,*)=avg(abs(xfft_array),0)
	mean_yfft_abs(i,*,*)=avg(abs(yfft_array),0)
	mean_xfft_pha(i,*,*)=avg(atan(imaginary(xfft_array),$
					  real_part(xfft_array)),0)
	mean_yfft_pha(i,*,*)=avg(atan(imaginary(yfft_array),$
				          real_part(yfft_array)),0)
	loadct,11
	!p.multi=[0,1,2]
	tvmean_xfft_abs=fltarr(25,25)
	tvmean_xfft_abs(*,*)=mean_xfft_abs(i,*,*)
	tvmean_xfft_abs=tvmean_xfft_abs/1.2
	tvmean_xfft_re=fltarr(25,25)
	print,max(tvmean_xfft_abs)

        ;Rearrange indices
        for n=0,num-1 do begin
            for m=0,num-1 do begin
                tvmean_xfft_re(subn(n),subm(m))=tvmean_xfft_abs(n,m)
            endfor
        endfor


	tvmean_yfft_abs=fltarr(25,25)
        tvmean_yfft_abs(*,*)=mean_yfft_abs(i,*,*)
        tvmean_yfft_abs=tvmean_yfft_abs/1.2
        tvmean_yfft_re=fltarr(25,25)
	print,max(tvmean_yfft_abs)

	;Rearrange Arrays
        for n=0,num-1 do begin
            for m=0,num-1 do begin
                tvmean_yfft_re(subn(n),subm(m))=tvmean_yfft_abs(n,m)

            endfor
        endfor

	s1=sort(tvmean_xfft_abs)
	s2=sort(tvmean_yfft_abs)
	
	max1=tvmean_xfft_abs(s1(N_elements(s1)-1))
	max2=tvmean_yfft_abs(s2(N_elements(s2)-1))
	bothmax=max([max1,max2])

	tv,(255*congrid(tvmean_xfft_re,625,625)),.5,1,$
	   xsize=5,ysize=5,/inches
	tv,(255*congrid(tvmean_yfft_re,625,625)),6,1,$
	   xsize=5,ysize=5,/inches
	colorbar,range=[0.0,0.6],format='(F5.3)',position=[0.83,.1,.9,.5]

	xyouts,.6,.83,strtrim(loc_e1(i),2)+' < e1 < ' +$
		       strtrim(loc_e1(i)+binsize,2),/normal
	xyouts,.6,.93,'Fourier Components (Magnitude) of X Offsets (Left)',$
	/normal
	xyouts,.6,.88,'Fourier Components (Magnitude) of Y Offsets (Right)',$
	/normal	
	xyouts,.6,.78,'Number of events= '+strtrim(e1hist(i),2),/normal
	xyouts,.3,.08,'Kx',/normal
	xyouts,.6,.08,'Kx',/normal
	xyouts,.93,.4,'Ky',/normal
	erase
	
	loadct,0

	shade_surf,mean_xfft_pha(i,*,*),background=white,$
                color=black,xtitle='Kx',ytitle='Ky',$
                position=[.1,.1,.9,.9],/normal
        xyouts,.1,.93,strtrim(loc_e1(i),2)+' < e1 < ' +$
                       strtrim(loc_e1(i)+binsize,2),$
		       /normal
        xyouts,.1,.98,'Fourier Components (Phase) of X Offsets',/normal        
        xyouts,.5,.93,'Number of events= '+strtrim(e1hist(i),2),/normal

	erase

	shade_surf,mean_yfft_pha(i,*,*),background=white,$
                color=black,xtitle='Kx',ytitle='Ky',$
                position=[.1,.1,.9,.9],/normal
        xyouts,.1,.93,strtrim(loc_e1(i),2)+' < e1 < ' +$
                       strtrim(loc_e1(i)+binsize,2),$
	 	       /normal
        xyouts,.1,.98,'Fourier Components (Phase) of Y Offsets',/normal 
        xyouts,.5,.93,'Number of events= '+strtrim(e1hist(i),2),/normal

	erase	
endfor


plot,frame_e1
plot,frame_e2

!p.multi=[0,1,1]
plot,frame_e1,frame_e2,psym=2

device,/close
set_plot,'x'

stop

end

