pro final_fft_ellip_theta_histogram,$
	Date=Date,stats_directory=stats_directory,$
        depth=depth,x_off=x_off,y_off=y_off,exam_date=exam_date,$
	xscale=xscale,yscale=yscale,low=low,lcut=lcut,$
        xWindSpeed=xWindSpeed,yWindSpeed=yWindSpeed,$
	phase_dim=phase_dim

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(depth)) then depth=0
If (not keyword_set(x_off)) then x_off=0
If (not keyword_set(y_off)) then y_off=0
If (not keyword_set(xscale)) then xscale=0
If (not keyword_set(yscale)) then yscale=0
If (not keyword_set(low)) then low=1.8
If (not keyword_set(lcut)) then lcut=10000.0
If (not keyword_set(xWindSpeed))  then xWindSpeed=17.    ;m/sec
If (not keyword_set(yWindSpeed))  then yWindSpeed=4.    ;m/sec
If (not keyword_set(phase_dim)) then phase_dim=long(1024)
If (not keyword_set(MultiWindSpeeds)) then MultiWindSpeeds=1

postscript_dir='/afs/slac/u/ki/lances/latex/'

;OPTIONS: **********************************************************
;Avg_sub=0 means do not subract any offsets
;Avg_sub=1 means subtract the average offset of a given frame
;Avg_sub=2 means set the kx=0,ky=0 componenet of FFT to zero
;Avg_sub=3 means subtract a scaled offset from the distorted grid
avg_sub=1
;fitted_only=1 means use only the fitted points on each frame
;fitted_only=2 means use only the non-fitted points one each frame
;fitted_only=3 means use only points that have non-zero value
fitted_only=1
;UseFid=1      apply the fiducial cut
;UseFid=0      don't apply the fiducial cut
UseFid = 0

;LABELS: ***************************************************************
if fitted_only eq 1 then begin
   label_string='_fitted_only_'
endif else if fitted_only eq 2 then begin
   label_string='_not_fitted_only_'
endif else if fitted_only eq 3 then begin
   label_string='_non_zero_only_'
endif

if avg_sub eq 0 then begin
   label_string=label_string+'no_sub_'
endif else if avg_sub eq 1 then begin
   label_string=label_string+'avg_sub_'
endif else if avg_sub eq 2 then begin
   label_string=label_string+'kx0_ky0_sub_'
endif else begin
   label_string=label_string
endelse

if UseFid eq 0 then begin
   label_string = '_NoFiducial_'
endif else begin
   label_string = '_Fiducial_'
endelse

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

;SIMULATION OFFSETS GATHERED FROM final_phase_screen_angle_histogram_no_text
xwind=strtrim(string(xWindSpeed,format='(%"%10.1f")'),2)
ywind=strtrim(string(yWindSpeed,format='(%"%10.1f")'),2)
lowinmeters=phase_dim*.17/low
highinmeters=phase_dim*.17/lcut
lowinmeters=strtrim(string(lowinmeters,format='(%"%10.2f")'),2)
highinmeters=strtrim(string(highinmeters,format='(%"%10.4f")'),2)

if MultiWindSpeeds eq 0 then begin
label_string=label_string+'low_'+lowinmeters+'_high_'+highinmeters+$
             '_xwind_'+xwind+'_ywind_'+ywind
endif else begin
label_string=label_string+'low_'+lowinmeters+'_high_'+highinmeters+$
             '6xwinds6ywinds'
endelse

atm_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/lsst_sim/atmosphere/'
atm_file=atm_dir+'ellipticity_low_'+lowinmeters+$
        '_high_'+highinmeters+'.txt'

if avg_sub eq 1 then begin 
   atm_file=atm_dir+'ellipticity_avg_sub_low_'+label_string+'.txt'
endif else if avg_sub eq 2 then begin
   atm_file=atm_dir+'ellipticity_kx0_ky0_sub_low_'+label_string+'.txt'
endif else begin
   atm_file=atm_dir+'ellipticity_low_'+lowinmeters+$
        '_high_'+highinmeters+'.txt'
endelse

readcol,atm_file,theta_sim,e_sim,format='f,f'
;Distorted Offsets

;ARRAYS OF FILES
all_nights=1
if all_nights eq 1 then begin
  ;--10th
  stats_directory1='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-10'+$
          '/ShackHartman/'+stats_directory+'/'
  offset_files1=file_search(stats_directory1+'*5132*sub_offsets.txt',$
          count=num_stats_files1)
  readcol,stats_directory1+'Average_offsets.txt',xoff1,yoff1,$
          format='X,X,f,X,f'
  raw_files1=file_search('/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-10'+$
	'/ShackHartman/raw/'+'*raw.fits',count=num_raw_files1)
  readcol,stats_directory1+'newest_x_centers.txt',XCenter1,XFCenter1,$
	format='X,f,f'
  readcol,stats_directory1+'newest_y_centers.txt',YCenter1,YFCenter1,$
	format='X,f,f'

  ;--11th
  stats_directory2='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-11'+$
          '/ShackHartman/'+stats_directory+'/'
  offset_files2=file_search(stats_directory2+'*5132*sub_offsets.txt',$
          count=num_stats_files2)
  readcol,stats_directory2+'Average_offsets.txt',xoff2,yoff2,$
          format='X,X,f,X,f'
  raw_files2=file_search('/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-11'+$
        '/ShackHartman/raw/'+'*raw.fits',count=num_raw_files2)
  readcol,stats_directory2+'newest_x_centers.txt',XCenter2,XFCenter2,$
	format='X,f'
  readcol,stats_directory2+'newest_y_centers.txt',YCenter2,YFCenter2,$
	format='X,f'

  ;--12th
  stats_directory3='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-12'+$
          '/ShackHartman/'+stats_directory+'/'
  offset_files3=file_search(stats_directory3+'*5132*sub_offsets.txt',$
          count=num_stats_files3)
  readcol,stats_directory3+'Average_offsets.txt',xoff3,yoff3,$
          format='X,X,f,X,f'
  raw_files3=file_search('/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-12'+$
        '/ShackHartman/raw/'+'*raw.fits',count=num_raw_files3)
  readcol,stats_directory3+'newest_x_centers.txt',XCenter3,XFCenter3,$
	format='X,f,f'
  readcol,stats_directory3+'newest_y_centers.txt',YCenter3,YFCenter3,$
	format='X,f,f'

  num_stats_files=num_stats_files1+num_stats_files2+num_stats_files3
  offset_files=[offset_files1,offset_files2,offset_files3]  
  raw_files=[raw_files1(0:num_stats_files1-1),$
	     raw_files2(0:num_stats_files2-1),$
	     raw_files3(0:num_stats_files3-1)]
  xcenters=[XCenter1(0:num_stats_files1-1),$
	    XCenter2(0:num_stats_files2-1),$
	    XCenter3(0:num_stats_files3-1)]
  ycenters=[YCenter1(0:num_stats_files1-1),$
            YCenter2(0:num_stats_files2-1),$
            YCenter3(0:num_stats_files3-1)]
  xfcenters=[XFCenter1(0:num_stats_files1-1),$
             XFCenter2(0:num_stats_files2-1),$
             XFCenter3(0:num_stats_files3-1)]
  yfcenters=[YFCenter1(0:num_stats_files1-1),$
             YFCenter2(0:num_stats_files2-1),$
             YFCenter3(0:num_stats_files3-1)]

  Date='All_Nights'

  ;These are the offsets of the distorted grid
  xoff=(xoff1+xoff2+xoff3)/3
  yoff=(yoff1+yoff2+yoff3)/3
endif else begin
  stats_directory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$
          '/ShackHartman/'+stats_directory+'/'
  offset_files=file_search(stats_directory+'*5132*sub_offsets.txt',$
          count=num_stats_files)
  Date='5_'+Date+'_'
endelse

;WINDS******************************************************************
wind_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/winds/'
readcol,WIND_dir+'WIND_numbers.txt',Wyear,Wmonth,Wdate,Whr,Wminute,$
	Wsec,Wspeed,Wdir,format='(d,d,f,f,f,f,f,f)'
WHrsDec=WDate+Whr/24+Wminute/(24*60)+Wsec/(24*3600)
WMatchRange=2.75/(24*60)

;POINTING***************************************************************
pointing_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/'
readcol,pointing_dir+'SOAR_tracking.txt',PDateArr,PTime,PAz,PAlt,$
	format='a,a,f,f'
NoPlines=N_elements(PDateArr)
PDate=fltarr(NoPlines)
PMatchRange=2.55/(24*60)

for i=0, N_elements(PDate)-1 do begin
  date_arr=strsplit(PDateArr(i),'-',/extract)
  time_arr=strsplit(PTime(i),':',/extract)
  
  date_val=valid_num(date_arr(2),PDateInd)
  sec_val=valid_num(time_arr(2),PSecInd)
  min_val=valid_num(time_arr(1),PMinInd)
  hr_val=valid_num(time_arr(0),PHrInd)
 
  PHrInd=PHrInd-4   ;Account for 4 hour offset
  PDate(i)=PDateInd+PHrInd/24+PMinInd/(24*60)+PSecInd/(24*3600) 

endfor
;***********************************************************
;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

;Grid for partvelvec
x_ref_col=12.5*subx+542
y_ref_col=25.0*suby+460
index=findgen(num*num)
gridx=fltarr(num*num)
gridy=fltarr(num*num)
gridx(index)=x_ref_col(index mod num)
gridy(index)=y_ref_col(index/num)
                                                                                  
;Indices and such
no_bins=8
a=findgen(13)+12
b=findgen(12)
subn=[a,b]
subm=subn

;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)
fittedspots_arr=intarr(num_stats_files,num*num)

;WINDDIRECTIONS AND POINTINGS
WSpeedArr=fltarr(num_stats_files)
WDirArr=fltarr(num_stats_files)
PAzArr=fltarr(num_stats_files)
PAltArr=fltarr(num_stats_files)

;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)
timeeden=fltarr(num_stats_files-depth,num,num)

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

;for high pass filter
frame_e1=fltarr(num_stats_files)
frame_e2=fltarr(num_stats_files)
frame_e=fltarr(num_stats_files)
frame_theta=fltarr(num_stats_files)
frame_e_den=fltarr(num_stats_files)

kx0_array=fltarr(num_stats_files)
ky0_array=fltarr(num_stats_files)
phi_array=fltarr(num_stats_files)
r_array=fltarr(num_stats_files)

;Setup file for output
get_lun,e_file
openw,e_file,'/nfs/slac/g/ki/ki08/lsst/CPanalysis/Ellip_Wind_Pointing.txt'

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_off,y_off,format='X,X,f,X,f'

	;FIND OUT WHAT WIND WAS AND POINTING
        unshifted_filename=Raw_files(file_number)
        fits_read,unshifted_filename,unshifted_image,header
        uttime=SXPAR(header,'DATE-OBS')

  	first_arr=strsplit(uttime,'-',/extract)
	second_arr=strsplit(first_arr(2),'T',/extract)
	third_arr=strsplit(second_arr(1),':',/extract)

	yr_val=valid_num(first_arr(0),year)
	month_val=valid_num(first_arr(1),month)
	date_val=valid_num(second_arr(0),DecDate)

	sec_val=valid_num(third_arr(2),secs)
	min_val=valid_num(third_arr(1),mins)
	hr_val=valid_num(third_arr(0),hrs)
	mins=mins+secs/60
	hrs=hrs+mins/60
	DecDate=DecDate+hrs/24
	;Match To Wind Direction and Speed
        wmatch=where(abs(WHrsDec - DecDate) le WMatchRange)
        WSpeedArr(file_number)=WSpeed(wmatch(0))
        WDirArr(file_number)=WDir(wmatch(0))

	;Match to Pointings
	pmatch=where(abs(PDate-DecDate) le PMatchRange)
        PAzArr(file_number)=PAz(pmatch(0))
	PAltArr(file_number)=PAlt(pmatch(0))
	;Obtain the difference ordered in COLUMN, ROW between centroid and grid         point
	fittedspots=where(x_off ne -10)
        fittedspots_arr(file_number,fittedspots)=1
	x_off(where(x_off eq -10))=0
        y_off(where(y_off eq -10))=0
          
	x_offset=fltarr(num,num)
	y_offset=fltarr(num,num)
	x_offset(*,*)=x_off
	y_offset(*,*)=y_off

        ;Avg_sub=1 means subtract the average offset of a given frame
	;Avg_sub=2 means set the kx=0,ky=0 componenet of FFT to zero
	;Avg_sub=3 means subtract a scaled offset from the distorted grid
        if avg_sub eq 1 then begin
           xfit=where(x_offset ne 0)
           yfit=where(y_offset ne 0)
           mean_xoff=mean(x_offset(xfit))
	   mean_yoff=mean(y_offset(yfit))
           x_offset(xfit)=x_offset(xfit)-mean_xoff
           y_offset(yfit)=y_offset(yfit)-mean_yoff
           phi_array(file_number)=atan(mean_yoff,mean_xoff)
           r_array(file_number)=sqrt(mean_xoff^2+mean_yoff^2)
 
        endif else if avg_sub eq 2 then begin
          x_fft=complexarr(num,num)
          y_fft=complexarr(num,num)
          x_fft(*,*)=fft(x_offset)
          y_fft(*,*)=fft(y_offset)
          x_fft_re=complexarr(num,num)
	  y_fft_re=complexarr(num,num)
       	  for n=0,num-1 do begin
            for m=0,num-1 do begin
                x_fft_re(subn(n),subm(m))=x_fft(n,m)
                y_fft_re(subn(n),subm(m))=y_fft(n,m)
            endfor
          endfor
          kx0_array(file_number)=x_fft_re(12,12)
	  ky0_array(file_number)=y_fft_re(12,12)
          phi_array(file_number)=atan(y_fft_re(12,12),x_fft_re(12,12))
          r_array(file_number)=sqrt(x_fft_re(12,12)^2+y_fft_re(12,12)^2)

          x_fft_re(12,12)=0
	  y_fft_re(12,12)=0

          for n=0,num-1 do begin
            for m=0,num-1 do begin
                x_fft(n,m)=x_fft_re(subn(n),subm(m))
                y_fft(n,m)=y_fft_re(subn(n),subm(m))
            endfor
          endfor

          x_offset=fft(x_fft,/inverse)
          y_offset=fft(y_fft,/inverse)
        
        
        endif else if avg_sub eq 3 then begin
           xfit=where(x_offset ne 0)
           yfit=where(y_offset ne 0)
           x_offset(xfit)=x_offset(xfit)+xscale*xoff(xfit)
           y_offset(yfit)=y_offset(yfit)+yscale*yoff(yfit)
	
        endif

        
        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)


   ;GET READY TO USE ONLY THE FITTED SPOTS
   fittedspots=fittedspots_arr(file_number,*)
   xfit=where(fittedspots eq 1)
   notfit=where(fittedspots ne 1)

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

   timeeden(file_number,*,*)=(x_offset_hor(file_number,*,*)^2+$
                              y_offset_hor(file_number,*,*)^2)

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

   if fitted_only eq 1 then begin
        frame_e_den=mean(timeedtemp(xfit))
        frame_e1(file_number)=mean(timee1temp(xfit))/$
                         mean(timeedtemp(xfit))
        frame_e2(file_number)=mean(timee2temp(xfit))/$
                         mean(timeedtemp(xfit))
   endif else if fitted_only eq 2 then begin
        frame_e_den=mean(timeedtemp(notfit))
        frame_e1(file_number)=mean(timee1temp(notfit))/$
                         mean(timeedtemp(notfit))
        frame_e2(file_number)=mean(timee2temp(notfit))/$
                         mean(timeedtemp(notfit))
   endif else begin
        frame_e_den=mean(timeedtemp(where(timee1temp ne 0)))  
        frame_e1(file_number)=mean(timee1temp(where(timee1temp ne 0)))/$
                         mean(timeedtemp(where(timee1temp ne 0)))
        frame_e2(file_number)=mean(timee2temp(where(timee2temp ne 0)))/$
                         mean(timeedtemp(where(timee2temp ne 0)))
   endelse
   frame_e(file_number)=sqrt(frame_e1(file_number)^2+frame_e2(file_number)^2)
                                                                                
   signe1=abs(frame_e1(file_number))/frame_e1(file_number)
   signe2=abs(frame_e2(file_number))/frame_e2(file_number)
                      
   ;Get the position angle associated with the elongation of the object
   ;Should be between pi/4 and pi/2 if e1 is negative and between 0 and pi/2 if 
   ;e1 is positive.  Should be in second quadrant if e2 is negative
   frame_theta(file_number)=$
   signe2*asin(sqrt(.5*(1-signe1*sqrt(1-frame_e2(file_number)^2/$
                                   frame_e(file_number)^2))))
                                                                                
   ;Recalculate e1 and e2 with new cosine angle
   frame_e1(file_number)=frame_e1(file_number)*cos(frame_theta(file_number))+$
                         frame_e2(file_number)*sin(frame_theta(file_number))
   frame_e2(file_number)=frame_e1(file_number)*(-sin(frame_theta(file_number)))+$
                         frame_e2(file_number)*cos(frame_theta(file_number))
   
   ;PRINT OUT TO THE FILE
   basename=file_basename(offset_files(file_number))
   printf,e_file,$
	 FORMAT=$
	'(%"%s\t%s\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f")',$
	 basename,uttime,DecDate,$
	 frame_e(file_number),frame_e1(file_number),frame_e2(file_number),$
	 frame_theta(file_number),phi_array(file_number),$
	 WSpeedArr(file_number),WDirArr(file_number),$
	 PAltArr(file_number),PAzArr(file_number),$
	 xcenters(file_number),ycenters(file_number),$
	 XFcenters(file_number),yFCenters(file_number)
endfor

close,e_file
free_lun,e_file


;DO THE HISTOGRAMS ***************************************************
;GET THE RATIO OF SIMULATED EVENTS TO DATA EVENTS
SDratio=float(n_elements(frame_theta))/float(n_elements(theta_sim))
;All Nights
thetamin=-!pi/2				  ;THETA
no_bins=20.
thetabinwidth=!pi/no_bins
thetamax=thetamin+(no_bins-1)*thetabinwidth

thetahist=histogram(frame_theta,Nbins=no_bins,locations=loc_theta,$
		 reverse_indices=theta_ind,$
		 max=thetamax,min=thetamin)
thetabinsize=(thetamax-thetamin)/(no_bins-1)
thetasimhist=histogram(theta_sim,Nbins=no_bins,locations=loc_theta_sim,$
                 reverse_indices=theta_sim_ind,$
                 max=thetamax,min=thetamin)


emin=0					   ;ELLIPTICITY
noe_bins=20.
ebinwidth=(1-emin)/noe_bins
emax=emin+(noe_bins)*ebinwidth
ehist=histogram(frame_e,Nbins=noe_bins,locations=loc_e,$
                 reverse_indices=e_ind,$
		 max=emax,min=emin)
ebinsize=(emax-emin)/(noe_bins-1)

esimhist=histogram(e_sim,Nbins=noe_bins,locations=loc_e_sim,$
                 reverse_indices=e_sim_ind,$
                 max=emax,min=emin)


;Separate Nights
if all_nights eq 1 then begin

first_night=num_stats_files1-1
second_night=num_stats_files1+num_stats_files2-1
third_night=num_stats_files1+num_stats_files2+num_stats_files3-1

thetahist1=histogram(frame_theta(0:first_night),$
                 Nbins=no_bins,$
                 locations=loc_theta1,$
                 reverse_indices=theta_ind1,$
                 max=thetamax,min=thetamin)
thetahist2=histogram(frame_theta(first_night+1:second_night),$
                 Nbins=no_bins,$
                 locations=loc_theta2,$
                 reverse_indices=theta_ind2,$
                 max=thetamax,min=thetamin)
thetahist3=histogram(frame_theta(second_night+1:third_night),$
                 Nbins=no_bins,$
                 locations=loc_theta3,$
                 reverse_indices=theta_ind3,$
                 max=thetamax,min=thetamin)

ehist1=histogram(frame_e(0:first_night),$
                 Nbins=no_bins,$
                 locations=loc_e1,$
                 reverse_indices=e_ind1,$
                 max=emax,min=emin)
ehist2=histogram(frame_e(first_night+1:second_night),$
                 Nbins=no_bins,$
                 locations=loc_e2,$
                 reverse_indices=e_ind2,$
                 max=emax,min=emin)
ehist3=histogram(frame_e(second_night+1:third_night),$
                 Nbins=no_bins,$
                 locations=loc_e3,$
                 reverse_indices=e_ind3,$
                 max=emax,min=emin)

endif


phimin=-!pi                                ;CENTROID OFFSET ANGLE
no_bins=20.
phibinwidth=2*!pi/no_bins
phimax=phimin+(no_bins-1)*phibinwidth


phihist=histogram(phi_array,Nbins=no_bins,locations=loc_phi,$
                 reverse_indices=phi_ind,$
                 max=phimax,min=phimin)
phibinsize=(phimax-phimin)/(no_bins-1)

;PLOTS: ***************************************************************

set_plot,'ps'
device,filename=postscript_dir+'ThetaHistogram_'+Date+'plus_phase_screen'+$
      label_string+'.ps'

!p.multi=[0,1,1]
loadct,39                       ;load color table 39
device,/color                   ;allow color on the postscript
device,ysize=7.5,/inches        ;Height of plot in y
device,xsize=7.0,/inches        ;Width of plot in x
device,yoffset=0.50,/inches      ;Y position of lower left corner
device,xoffset=0.00,/inches
white='FFFFFF'x
black='000000'x
!P.CHARSIZE=.9
!P.CHARTHICK=1.8
!P.THICK=4

;THETA_HISTOGRAM*************************************************
plot,loc_theta_sim+thetabinsize/2,SDratio*thetasimhist,psym=10,$
        linestyle=2,$
        xtitle='Theta',$
        ytitle='Number of Events',$
        background=white,color=black,$
        yrange=[0,320]
oplot,loc_theta+thetabinsize/2,thetahist,psym=10,color=black

device,/close

device,filename=postscript_dir+'EHistogram_'+Date+'plus_phase_screen'+$
      label_string+'.ps'
plot,loc_e,ehist,psym=10,xtitle='e',$
     ytitle='Number of Events',$
     background=white,color=black
oplot,loc_e_sim,SDratio*esimhist,color=black,psym=10,linestyle=2

device,/close

set_plot,'z'
erase
device, set_font='Courier'
device,set_resolution=[800,600]
!p.charsize=1
!p.charthick=3
!x.thick=2
!y.thick=2
!p.thick=3
!p.noerase=0
loadct,39

white='FFFFFF'x
black='000000'x
blue='FF0000'x
red='0000FF'x

!p.multi=[0,1,1]
erase
;THETA_HISTOGRAM*************************************************
plot,loc_theta_sim+thetabinsize/2,SDratio*thetasimhist,psym=10,$
        linestyle=2,$
	xtitle='Theta',$
        ytitle='Number of Events',$
        background=white,color=black,$
	yrange=[0,320]
oplot,loc_theta+thetabinsize/2,thetahist,psym=10,color=black
  jpgimg=tvrd()
  tvlct,reds,greens,blues,/get
write_png,postscript_dir+'Theta_Histogram_'+Date+'plus_phase_screen'+$
      '_low_'+lowinmeters+'_high_'+highinmeters+label_string+'.png',$
      jpgimg,reds,greens,blues
 
erase

;E_HISTOGRAM*****************************************************
plot,loc_e,ehist,psym=10,xtitle='e',$
     ytitle='Number of Events',$
     background=white,color=black
oplot,loc_e_sim,SDratio*esimhist,color=black,psym=10,linestyle=2

jpgimg=tvrd()
tvlct,reds,greens,blues,/get
     write_png,postscript_dir+'E_Histogram_'+Date+'plus_phase_screen'+$
	'_low_'+lowinmeters+'_high_'+highinmeters+label_string+'.png',$
        jpgimg,reds,greens,blues


if all_nights eq 1 then begin

erase
items_e=['10','11','12']
linestyle_e=[1,2,3]
sym_e=[10,10,10]
textcolors_e=[black,black,black]
colors_e=[black,20,40]

plot,loc_e1,ehist1,psym=10,xtitle='e',$
     ytitle='Number of Events',$
     background=white,color=black,linestyle=1
oplot,loc_e2,ehist2,psym=10,color=20,linestyle=2
oplot,loc_e3,ehist3,psym=10,color=40,linestyle=3
legend,items_e,linestyle=linestyle_e,textcolors=textcolors_e,$
	       colors=colors_e
jpgimg=tvrd()
tvlct,reds,greens,blues,/get
     write_png,postscript_dir+'E_Histogram_'+Date+'.png',$
        jpgimg,reds,greens,blues
erase


plot,loc_theta3,thetahist3,psym=10,xtitle='e',linestyle=3,$
     ytitle='Number of Events',$
     background=white,color=40
oplot,loc_theta2,thetahist2,psym=10,color=20,linestyle=2
oplot,loc_theta1,thetahist1,psym=10,color=black,linestyle=1
legend,items_e,linestyle=linestyle_e,textcolors=textcolors_e,$
               colors=colors_e
jpgimg=tvrd()
tvlct,reds,greens,blues,/get
     write_png,postscript_dir+'Theta_Histogram_'+Date+'.png',$
        jpgimg,reds,greens,blues
erase

endif

jpgimg=tvrd()
tvlct,reds,greens,blues,/get
     write_png,postscript_dir+'E_Histogram_'+Date+'plus_phase_screen'+$
        '_low_'+lowinmeters+'_high_'+highinmeters+label_string+'.png',$
        jpgimg,reds,greens,blues

erase

;PHI_HISTOGRAM***************************************************
plot,loc_phi+phibinsize/2,phihist,psym=10,xtitle='Theta',$
     ytitle='Number of Events',$
     title='Histogram of Centroid Angle for '+Date,$
     background=white,color=black

jpgimg=tvrd()
tvlct,reds,greens,blues,/get
write_png,postscript_dir+'Phi_Histogram_'+Date+$
         label_string+$
          '.png',$
      jpgimg,reds,greens,blues

erase
;THETA_PHI_SCATTER***********************************************
plot,phi_array,frame_theta,psym=2,$
     xtitle='Azimuthal angle of centroid',$
     ytitle='Ellipticity angle that maximizes e1 (Theta)',$
     title='Scatter Plot of Theta vs. Phi for'+Date,$
     background=white,color=black

jpgimg=tvrd()
tvlct,reds,greens,blues,/get
write_png,postscript_dir+'Theta_Phi_Scatter_'+Date+$
        label_string+$
        '.png',$
        jpgimg,reds,greens,blues

erase
;E_THETA_SCATTER**************************************************
device,set_resolution=[1400,1200]
!p.charsize=1
!p.charthick=1
!x.thick=2
!y.thick=2
!p.thick=1
!p.noerase=0

!p.multi=[0,2,2]
items_e=['e','e1','e2']
sym_e=[2,1,5]
colors_e=[black,blue,red]

plot,frame_e,frame_theta,psym=2,$
     xtitle='e',$
     ytitle='Ellipticity angle that maximizes e1 (Theta)',$
     title='Scatter Plot of e vs. Theta for '+Date,$
     background=white,color=black
plot,frame_e1,frame_theta,psym=2,$
     xtitle='e1',$
     ytitle='Ellipticity angle that maximizes e1 (Theta)',$
     title='Scatter Plot of e1 vs. Theta for '+Date,$
     background=white,color=black
plot,frame_e2,frame_theta,psym=2,$
     xtitle='e',$
     ytitle='Ellipticity angle that maximizes e1 (Theta)',$
     title='Scatter Plot of e2 vs. Theta for '+Date,$
     background=white,color=black

jpgimg=tvrd()
tvlct,reds,greens,blues,/get
write_png,postscript_dir+'Theta_E_Scatter_'+Date+$
        label_string+$
        '.png',$
        jpgimg,reds,greens,blues

erase
;E_PHI_SCATTER****************************************************
!p.multi=[0,2,2]
plot,frame_e,phi_array,psym=2,$
     xtitle='e',$
     ytitle='Azimuthal angle of centroid',$
     title='Scatter of e vs. Phi for '+Date,$
     background=white,color=black

plot,frame_e1,phi_array,psym=1,$
     xtitle='e1',$
     ytitle='Azimuthal angle of centroid',$
     title='Scatter of e1 vs. Phi for '+Date,$
     background=white,color=black

plot,frame_e2,phi_array,psym=5,$
     xtitle='e2',$
     ytitle='Azimuthal angle of centroid',$
     title='Scatter of e2 vs. Phi for '+Date,$     
     background=white,color=black

jpgimg=tvrd()
tvlct,reds,greens,blues,/get
write_png,postscript_dir+'Phi_E_Scatter_'+Date+$
        label_string+$
        '.png',$
        jpgimg,reds,greens,blues




device,/close
set_plot,'x'

stop

end

