pro final_fft_ellip_theta_histogram_plotff,$
	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,$
	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(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

  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)

  ;--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)

  ;--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)

;SIMULATION OFFSETS GATHERED FROM final_phase_screen_angle_histogram_no_text
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

ps_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/final_paper_figs/'

readcol,'/nfs/slac/g/ki/ki08/lsst/CPanalysis/Ellip_Wind_Pointing.txt',$
	basename,uttime,date,frame_e,frame_e1,frame_e2,$
	frame_theta,phi_array,WSpeedArr,WDirArr,$
	PAltArr,PAzArr,xcenters,ycenters,xfcenters,yfcenters,$
	format='a,a,f,f,f,f,f,f,f,f,f,f,f,f,f,f'

;SORT BY TIME SINCE THEY ARE OUT OF ORDER
ro=sort(uttime)
basename=basename(ro)
date=date(ro)
frame_e=frame_e(ro)
frame_e1=frame_e1(ro)
frame_e2=frame_e2(ro)
frame_theta=frame_theta(ro)
phi_array=phi_array(ro)
WSpeedArr=WSpeedArr(ro)
WDirArr=WDirArr(ro)
PAltArr=PAltArr(ro)
PAzArr=PAzArr(ro)
xcenter=xcenters(ro)
ycenters=ycenters(ro)
xfcenters=xfcenters(ro)
yfcenters=yfcenters(ro)

;***********************************************************************
;GET WIND DIRECTION ON CCD AND FIGURE OUT DIFFERENCE BETWEEN WIND AND E1
;MAKE DEGREES RANGE FROM -180 to +180
WDirArrComp=WDirArr		;MAINTAIN COMPASS
PAzArrComp=PAzArr		;MAINTAIN COMPASS
WDirArr(where(WDirArr gt 180))=WDirArr(where(WDirArr gt 180))-360
PAzArr(where(PAzArr gt 180))=PAzArr(where(PAzArr gt 180))-360
;GET DIFFERENCE BETWEEN WIND AND POINTING
XAzArr=PAzArr+90
XAzArr(where(XAzArr gt 180))=XAzArr(where(XAzArr gt 180))-360
phi=(PAzArr+90)-WDirArr
phiComp=phi
phi(where(phi gt 180))=phi(where(phi gt 180))-360
phiplus=where(phi gt 90)
phiminus=where(phi lt -90)
if phiplus(0) ne -1 then phi(phiplus)=phi(phiplus)-180
if phiminus(0) ne -1 then phi(phiminus)=phi(phiminus)+180

;CONVERT TO DEGREES FOR COMPARISON
frame_thetadeg=frame_theta*180/(!pi)

;GET TRACKING SLOPES
num_files=n_elements(PAzArr)
PAzSlopes=PAzArr(1:num_files-1)-PAzArr(0:num_files-2)
PAltSlopes=PAltArr(1:num_files-1)-PAltArr(0:num_files-2)

png=0
if png eq 0 then begin 
;DO THE HISTOGRAMS ***************************************************
;************************************************************ALTITUDE


set_plot,'ps'
device,filename=ps_dir+'Altitude_Theta.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=12.0,/inches        ;Width of plot in x
device,yoffset=0.0,/inches      ;Y position of lower left corner
device,xoffset=0.0,/inches
white='FFFFFF'x
black='000000'x
blue='FF0000'x
red='0000FF'x
!p.charsize=1
!p.thick=3
!p.charthick=3

!p.multi=[0,2,1]
;All Nights
thetamin=-!pi/2				  ;THETA
no_bins=20.
thetabinwidth=!pi/no_bins
thetamax=thetamin+(no_bins-1)*thetabinwidth
thetabinoff=thetabinwidth/2
thetabinoffdeg=thetabinoff*180/!pi

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

;SEPARATE HISTOGRAM INTO TWO BINS
altcut=62.75
highalt=where(PAltArr gt altcut)
lowalt=where(PAltArr le altcut)

print,n_elements(highalt),n_elements(lowalt)
highthetahist=histogram(frame_theta(highalt),Nbins=no_bins,$
		 locations=loc_theta_high,$
                 reverse_indices=theta_ind_high,$
                 max=thetamax,min=thetamin)
lowthetahist=histogram(frame_theta(lowalt),Nbins=no_bins,$
                 locations=loc_theta_low,$
                 reverse_indices=theta_ind_low,$
                 max=thetamax,min=thetamin)

plot,loc_theta_low+thetabinoff,lowthetahist,psym=10,$
	color=black,background=white,position=[.1,.55,.9,.95],$
	title='Ellipticity Angle Distribution at Highest and Lowest Pointing'
oplot,loc_theta_high+thetabinoff,highthetahist,psym=10,$
	color=60
itemst=['1307 Events Alt <= 62.75 degrees','1296 Events Alt > 62.75 degrees']
psymt=[0,0]
colorst=[0,60]
legend,itemst,psym=psymt,colors=colorst

PhiWindETheta=abs(phi-frame_thetadeg)
redPhiWindE=where(PhiWIndETheta ge 90)
PhiWindETheta(redPhiWindE)=-(PhiWindETheta(redPhiWindE)-180)
PhiWindEThetaMax=(thetamax)*180/!pi
PhiWindEThetaMin=(thetamin+!pi/2)*180/!pi
highPhiWindEThetahist=histogram(PhiWindETheta(highalt),Nbins=no_bins,$
                 locations=loc_PhiWindETheta_high,$
                 reverse_indices=PhiWindETheta_ind_high,$
                 max=PhiWindEThetaMax,min=PhiWindEThetaMin)
lowPhiWindEThetahist=histogram(PhiWindETheta(lowalt),Nbins=no_bins,$
                 locations=loc_PhiWindETheta_low,$
                 reverse_indices=PhiWindETheta_ind_high,$
                 max=PhiWindEThetaMax,min=PhiWindEThetaMin)


plot,loc_PhiWindETheta_low+thetabinoffdeg,lowPhiWindEThetahist,psym=10,$
        color=black,background=white,position=[.1,.05,.9,.5],$
	title='Difference Angle between E1 and Wind Direction' 
oplot,loc_PhiWindETheta_high+thetabinoffdeg,highPhiWindEThetahist,psym=10,$
        color=60
itemst=['1307 Events Alt <= 62.75 degrees','1296 Events Alt > 62.75 degrees']
psymt=[0,0]
colorst=[0,60]
legend,itemst,psym=psymt,colors=colorst

!p.multi=[0,2,1]
;SEPARATE HISTOGRAM INTO TWO BINS
highaltcut=65.7
lowaltcut=54.60
highalt=where(PAltArr gt highaltcut)
lowalt=where(PAltArr le  lowaltcut)

print,n_elements(highalt),n_elements(lowalt)
highthetahist=histogram(frame_theta(highalt),Nbins=no_bins,$
                 locations=loc_theta_high,$
                 reverse_indices=theta_ind_high,$
                 max=thetamax,min=thetamin)
lowthetahist=histogram(frame_theta(lowalt),Nbins=no_bins,$
                 locations=loc_theta_low,$
                 reverse_indices=theta_ind_low,$
                 max=thetamax,min=thetamin)

plot,loc_theta_low+thetabinoff,lowthetahist,psym=10,$
        color=black,background=white,position=[.1,.55,.9,.95],$
	title='Ellipticity Angle Distribution at Highest and Lowest Pointing'

oplot,loc_theta_high+thetabinoff,highthetahist,psym=10,$
        color=60
itemst=['426 Events Alt <= 54.6 degrees','418 Events Alt > 65.7 degrees']
psymt=[0,0]
colorst=[0,60]
legend,itemst,psym=psymt,colors=colorst

highPhiWindEThetahist=histogram(PhiWindETheta(highalt),Nbins=no_bins,$
                 locations=loc_PhiWindETheta_high,$
                 reverse_indices=PhiWindETheta_ind_high,$
                 max=PhiWindEThetaMax,min=PhiWindEThetamin)
lowPhiWindEThetahist=histogram(PhiWindETheta(lowalt),Nbins=no_bins,$
                 locations=loc_PhiWindETheta_low,$
                 reverse_indices=PhiWindETheta_ind_high,$
                 max=PhiWindEThetaMax,min=PhiWindEThetaMin)


plot,loc_PhiWindETheta_low+thetabinoffdeg,lowPhiWindEThetahist,psym=10,$
        color=black,background=white,position=[.1,.05,.9,.5],$
        title='Difference Angle between E1 and Wind Direction' 
oplot,loc_PhiWindETheta_high+thetabinoffdeg,highPhiWindEThetahist,psym=10,$
        color=60
itemst=['426 Events Alt <= 54.6 degrees','418 Events Alt > 65.7 degrees']
psymt=[0,0]
colorst=[0,60]
legend,itemst,psym=psymt,colors=colorst


!p.multi=[0,1,2]
plot, PAzArr,psym=1,color=black,background=white,position=[0.05,0.05,.9,.7]
oplot,XAzArr,psym=2,color=60
oplot,WDirArr,psym=4,color=120
oplot,Phi,psym=2,color=180
oplot,(Frame_Theta)*180/!pi,psym=1,color=240
oplot,PAltArr,psym=5,color=30

itemst=['Azimuthal Angle of telescope','Compass Angle of CCD x-axis',$
	'Compass Direction of Wind','Angle between Wind and x-axis',$
	'Ellipticity Angle','Altitude of telescope']
psymt=[1,2,4,2,1,5]
colorst=[0,60,120,180,240,30]
legend,itemst,psym=psymt,colors=colorst,pos=[.5,.65],/norm

plot,PAzSlopes,color=black,psym=2,background=white,$
	yrange=[-2,2],position=[0.05,.75,.9,.95]
oplot,PAltSlopes,color=120,psym=1
itemst=['Az Change','Alt Change']
psymt=[2,1]
colorst=[0,120]
legend,itemst,psym=psymt,colors=colorst,pos=[.05,.8],/norm

!p.multi=[0,1,3]
plot, PAzArr,psym=1,color=black,background=white,position=[0.05,0.05,.9,.45],$
	yrange=[-90,90]
oplot,XAzArr,psym=2,color=60
oplot,WDirArr,psym=4,color=120
oplot,Phi,psym=2,color=180
oplot,(Frame_Theta)*180/!pi,psym=1,color=240
oplot,PAltArr,psym=5,color=30

itemst=['Azimuthal Angle of telescope','Compass Angle of CCD x-axis',$
        'Compass Direction of Wind','Angle between Wind and x-axis',$
        'Ellipticity Angle','Altitude of telescope']
psymt=[1,2,4,2,1,5]
colorst=[0,60,120,180,240,30]
legend,itemst,psym=psymt,colors=colorst,pos=[.5,.55],/norm

!p.thick=2
plot,Xcenters,color=black,psym=1,background=white,$
        yrange=[530,550],position=[0.05,.5,.9,.7]
itemst=['Xcenter']
psymt=[1]
colorst=[0]
legend,itemst,psym=psymt,colors=colorst,pos=[.05,.6],/norm
plot,Ycenters,color=black,psym=1,position=[0.05,.75,.9,.925],$
	yrange=[450,480]
;oplot,XFcenters,color=120,psym=1
itemst=['Ycenter']
psymt=[1]
colorst=[0]
legend,itemst,psym=psymt,colors=colorst,pos=[.05,.8],/norm

!p.multi=[0,1,2]
plot,Xcenters,color=black,psym=1,background=white,$
        yrange=[530,550],$
        ytitle='Pixel Number'
;position=[0.05,.5,.9,.7]
        
plot,Ycenters,color=black,psym=1,$
        yrange=[450,480],$
        ytitle='Pixel Number',$
        xtitle='Frame'


!p.multi=[0,1,1]
plot,Xcenters,color=black,psym=1,background=white,$
        yrange=[530,550],$
	ytitle='Pixel Number'
;position=[0.05,.5,.9,.7]
	
plot,Ycenters,color=black,psym=1,$
        yrange=[450,480],$
	ytitle='Pixel Number'
	xtitle='Frame'
;position=[0.05,.75,.9,.925],$

!p.multi=[0,1,1]
plot,PhiWindETheta,psym=2,color=black,background=white,$
     title='Angle Between E1 and the Wind'
device,/close

endif

;E_THETA_SCATTER**************************************************
set_plot,'z'
erase
device,set_resolution=[1400,1200]
!p.charsize=3
!p.charthick=4.5
!x.thick=3
!y.thick=3
!p.thick=2.4
!p.noerase=0

yrmax=550
yrmin=530
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
!p.multi=[0,1,2]
plot,Xcenters,color=black,psym=1,background=white,$
        yrange=[yrmin,yrmax],$
        ytitle='Pixel Number',$
	xrange=[0,2700],xstyle=1

!p.thick=3.
oplot,[first_night,first_night],[yrmin,yrmax],linestyle=2,color=black
oplot,[second_night,second_night],[yrmin,yrmax],linestyle=2,color=black


yrmin=450
yrmax=480

!p.thick=2.4
plot,Ycenters,color=black,psym=1,$
        yrange=[450,480],$
        ytitle='Pixel Number',$
        xtitle='Frame',$
	xrange=[0,2700],xstyle=1

!p.thick=3.
oplot,[first_night,first_night],[yrmin,yrmax],linestyle=2,color=black
oplot,[second_night,second_night],[yrmin,yrmax],linestyle=2,color=black

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

stop

device,/close

;THE REST IS NOT IMPLEMENTED YET
;*************************************************************************
;*************************************************************************
;SEPARATE HISTOGRAM INTO THREE BINS  ;hialtcut=64.5 ;lowaltcut=59.32 gives equal
hialtcut=64.5
lowaltcut=59.32
highalt=where(PAltArr gt hialtcut)
medalt=where((PAltArr le hialtcut) and (PAltArr gt lowaltcut))
lowalt=where (PAltArr le lowaltcut)

print,n_elements(lowalt),n_elements(medalt),n_elements(highalt)
highthetahist=histogram(frame_theta(highalt),Nbins=no_bins,$
                 locations=loc_theta_high,$
                 reverse_indices=theta_ind_high,$
                 max=thetamax,min=thetamin)
medthetahist=histogram(frame_theta(medalt),Nbins=no_bins,$
                 locations=loc_theta_med,$
                 reverse_indices=theta_ind_med,$
                 max=thetamax,min=thetamin)
lowthetahist=histogram(frame_theta(lowalt),Nbins=no_bins,$
                 locations=loc_theta_low,$
                 reverse_indices=theta_ind_low,$
                 max=thetamax,min=thetamin)

plot,loc_theta_low+thetabinoff,lowthetahist,psym=10
oplot,loc_theta_med+thetabinoff,medthetahist,psym=10,linestyle=2,color=blue
oplot,loc_theta_high+thetabinoff,highthetahist,psym=10,linestyle=3,color=red

;SEPARATE HISTOGRAM INTO FOUR BINS
hialtcut=65.3
midaltcut=62.315
lowaltcut=57.32
highalt=where(PAltArr gt hialtcut)
medhalt=where((PAltArr le hialtcut) and (PAltArr gt midaltcut))
medlalt=where((PAltArr le midaltcut) and (PAltArr gt lowaltcut))
lowalt=where (PAltArr le lowaltcut)

print,n_elements(lowalt),n_elements(medlalt),n_elements(medhalt),n_elements(highalt)
highthetahist=histogram(frame_theta(highalt),Nbins=no_bins,$
                 locations=loc_theta_high,$
                 reverse_indices=theta_ind_high,$
                 max=thetamax,min=thetamin)
medhthetahist=histogram(frame_theta(medhalt),Nbins=no_bins,$
                 locations=loc_theta_medh,$
                 reverse_indices=theta_ind_medh,$
                 max=thetamax,min=thetamin)
medlthetahist=histogram(frame_theta(medlalt),Nbins=no_bins,$
                 locations=loc_theta_medl,$
                 reverse_indices=theta_ind_medl,$
                 max=thetamax,min=thetamin)
lowthetahist=histogram(frame_theta(lowalt),Nbins=no_bins,$
                 locations=loc_theta_low,$
                 reverse_indices=theta_ind_low,$
                 max=thetamax,min=thetamin)

plot,loc_theta_low+thetabinoff,lowthetahist,psym=10
oplot,loc_theta_medl+thetabinoff,medlthetahist,psym=10,linestyle=2,color=blue
oplot,loc_theta_medh+thetabinoff,medhthetahist,psym=10,linestyle=2,color=29999
oplot,loc_theta_high+thetabinoff,highthetahist+thetabinoff,psym=10,linestyle=3,color=red

;AZIMUTH******************************************************************
;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)

;SEPARATE HISTOGRAM INTO TWO BINS
azcut=220
highaz=where(PAzArr gt azcut)
lowaz=where(PAzArr le azcut)

print,n_elements(highaz),n_elements(lowaz)
highthetahist=histogram(frame_theta(highaz),Nbins=no_bins,$
                 locations=loc_theta_high,$
                 reverse_indices=theta_ind_high,$
                 max=thetamax,min=thetamin)
lowthetahist=histogram(frame_theta(lowaz),Nbins=no_bins,$
                 locations=loc_theta_low,$
                 reverse_indices=theta_ind_low,$
                 max=thetamax,min=thetamin)

plot,loc_theta_low+thetabinoff,lowthetahist,psym=10,$
        color=black,background=white
oplot,loc_theta_high+thetabinoff,highthetahist,psym=10,$
        color=60
itemst=['1307 Events Az <= 62.75 degrees','1296 Events Az > 62.75 degrees']
psymt=[0,0]
colorst=[0,60]
legend,itemst,psym=psymt,colors=colorst


					   ;THETA HIGH HIST
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)

;E HIST IN TWO BINS
highehist=histogram(frame_e(highalt),Nbins=no_bins,$
                 locations=loc_e_high,$
                 reverse_indices=e_ind_high,$
                 max=emax,min=emin)
lowehist=histogram(frame_e(lowalt),Nbins=no_bins,$
                 locations=loc_e_low,$
                 reverse_indices=e_ind_low,$
                 max=emax,min=emin)

plot,loc_e_low,lowehist,psym=10
oplot,loc_e_high,highehist,psym=10,linestyle=2

print,n_elements(lowalt),n_elements(medalt),n_elements(highalt)
highehist=histogram(frame_e(highalt),Nbins=no_bins,$
                 locations=loc_e_high,$
                 reverse_indices=e_ind_high,$
                 max=emax,min=emin)
medehist=histogram(frame_e(medalt),Nbins=no_bins,$
                 locations=loc_e_med,$
                 reverse_indices=e_ind_med,$
                 max=emax,min=emin)
lowehist=histogram(frame_e(lowalt),Nbins=no_bins,$
                 locations=loc_e_low,$
                 reverse_indices=e_ind_low,$
                 max=emax,min=emin)

plot,loc_e_low,lowehist,psym=10
oplot,loc_e_med,medehist,psym=10,linestyle=2,color=blue
oplot,loc_e_high,highehist,psym=10,linestyle=3,color=red




;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,'z'
erase
device, set_font='Courier'
device,set_resolution=[800,600]
!p.charsize=1
!p.charthick=1
!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]

;THETA_HISTOGRAM*************************************************
plot,loc_theta+thetabinsize/2,thetahist,psym=10,xtitle='Theta',$
     ytitle='Number of Events',$
     background=white,color=black

oplot,loc_theta_sim+thetabinsize/2,thetasimhist,psym=10,color=black,$
        linestyle=2

  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,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
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+label_string+'.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
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+label_string+'.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=6
!p.charthick=14
!x.thick=2
!y.thick=2
!p.thick=14
!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

