program EnergyDispersion

use bmad
implicit none

real(rp) eta/8.5/,Trev/149.0e-9/, Q/0.90387/,Qp/-0.13854/
real(rp) Delta0, omega
real(rp) t, dt/1.e-9/
real(rp) tn
real(rp), allocatable :: x(:), den(:), xod(:), intensity(:)

integer nmax, Turns/100/
integer n
integer i
integer lun, lun1
integer lun2
integer lun3
Delta0 = 0.0012
omega = twopi/Trev
lun=lunget()
open(unit=lun, file='xaverage_energy_1ns.dat')
lun1=lunget()
open(unit=lun1, file='xaverage_energy_10ns.dat')
lun2=lunget()
open(unit=lun2, file='xaverage_energy_100ns.dat')
print '(a,$)',' Turns ? '
read(5,*)turns

nmax=Turns * Trev/dt
allocate(x(0:nmax+1), den(0:nmax+1), xod(0:nmax+1), intensity(0:nmax+1))
write(lun,'(a,i10)'),' Number of turns = ',Turns
write(lun1,'(a,i10)'),' Number of turns = ',Turns
write(lun2,'(a,i10)'),' Number of turns = ',Turns
t=0.
i=0
x(i) = 0
do while (t < nmax * Trev .and. i+1 < nmax)
 i=i+1
 t=i*dt
 x(i)=0.
 den(i)=0.
 intensity(i)=0.
 do n=1,Turns * 2
  tn = t/float(n)/Trev - 1.0_rp
!  x(i) = x(i) + eta*tn*(1.0_rp-cos((Q+(Qp-Q)*tn)*omega*t))*exp(-tn**2/2/Delta0**2)
  if(abs(tn) < 0.02)then
   x(i) = x(i) + 1/sqrt(twopi)/Delta0 * eta*tn*(1.0_rp-cos((Q*n*Trev+Qp*(t-n*Trev))*omega))*exp(-tn**2/2/Delta0**2)
   den(i) = den(i) + 1/sqrt(twopi)/Delta0 *exp(-tn**2/2/Delta0**2)
   intensity(i) = intensity(i)+ exp(-(t/n/Trev -1)**2/2./Delta0**2)/sqrt(2*pi)/Delta0/n/Trev
  endif
!  print '(4es12.4,2i10,4es16.8)',Q,Qp,omega,Delta0,i,n,t,tn,x(i),exp(-tn**2/2/Delta0**2)
 end do
 if(abs(x(i)) < 1.e-12)x(i)=0
! print '(es12.4,i10,es16.8)', t,i,x(i)
 xod(i) = x(i)/den(i)
 write(lun,'(es12.4,i10,3e20.12)')t,i,x(i), x(i)/den(i), intensity(i)
 if(10*((i+5)/10) == i+5) write(lun1,'(es12.4,i10,3e20.12)')t,i,sum(x(i-4:i+5))/10., sum(xod(i-4:i+5))/10., sum(intensity(i-4:i+5))/10.
 if(100*((i+50)/100) == i+50) write(lun2,'(es12.4,i10,3e20.12)')t,i,sum(x(i-49:i+50))/100.,sum(xod(i-49:i+50))/100., sum(intensity(i-49:i+50))/100.
end do
close(unit=lun)
close(unit=lun1)
print '(a)', 'write to'
print  '(a)', 'xaverage_energy_1ns.dat' 
print '(a)', 'xaverage_energy_10ns.dat' 
print '(a)', 'xaverage_energy_100ns.dat'
end