program synch_bpm_anal
use nr
use precision_def
use multibunch_interface
implicit none
character*60 file_name
character*7 string
character*100 long_string
character*50 directory/'/nfs/cesr/online/instr/cbpm_data/tsharc/raw'/
character*2 sd
character*100 dir_and_file
integer lun, c(4),num, lunget
integer n
integer i,j,n_fft, start_turn, end_turn
integer ios
integer lun1, lun2, lun3
integer ix
integer nt/20000/
real(rp),allocatable:: x(:), y(:), ftx(:), fty(:)
real(rp) pie
real(rp) tunex(10), tuney(10)

namelist/sbpm_init/file_name,start_turn,n_fft, end_turn
     pie = atan(1.)*4
     lun = lunget()
     open(unit=lun, file='sbpm_init.dat', STATUS ='old')
     read(lun, nml=sbpm_init, IOSTAT=ios)
     rewind(lun)
     read(lun, nml=sbpm_init, IOSTAT=ios)
     close(unit=lun)
     write(6,nml=sbpm_init)
     sd = file_name(5:6)
     dir_and_file = trim(directory)//'/'//sd//'/'//trim(file_name)
     ix = index(file_name,'mac')
     if(ix /=0)dir_and_file = '/Users/dlr10/development/bmad_dist/darkphoton/machine_data/'//trim(file_name(1:ix-1))
     print *,dir_and_file
     open(unit=lun, file = dir_and_file)

    allocate(x(nt),y(nt),ftx(nt),fty(nt))

  lun1 = lunget()
  open (unit=lun1,file = 'but_position.dat')
  lun2 = lunget()
  open(unit=lun2, file= 'fft.dat')
  lun3 = lunget()
  open(unit=lun3, file= 'fft_peaks.dat')
  write(lun1,'(a10,4a10,4a12)')'turn','top','inside','bottom','outside','x','y','xraw','yraw'
  write(lun2,'(a10,5a12)')'turn','fraction', 'x','y','fftx','ffty'
  write(lun3, '(a6,a1,a5,a3,5a12,1x,5a12)')'turn n','-','n+1','tunex 1','2','3','4','5','tuney 1','2','3','4','5'
n=0
 do while (index(long_string,'BEGIN')==0)
  read(lun,*)long_string
 end do
 do while(.true.)
  read(lun,*, end=99)num,string,c(1:4)
  n=n+1
  if(sum(abs(c)) /=0)then
   x(n) = ((c(4))-(c(2)))/float(sum((c)))
   y(n) = ((c(1))-(c(3)))/float(sum((c)))
  else
   x(n)=0
   y(n)=0
  endif
  print '(i10, 4i10,2es12.4,1x,i10,1x,2es12.4)',n, c, x(n), y(n), sum(c), x(n)*sum(c),y(n)*sum(c)
  write(lun1,'(i10,4i10,4es12.4)')n,c,x(n),y(n), x(n)*sum(c),y(n)*sum(c)
end do
99 continue
print '(i10,1x,a)',n,'data points'

      j=start_turn
     print '(a10,a60,a60)','turns','horizontal (1:5)','vertical (1:5)'
      do while(j+n_fft < end_turn)
        forall(i=1:n_fft)ftx(i) = x(j+i-1) *(2*(sin((pie*i)/n_fft))*(sin((pie*i)/n_fft)))
        forall(i=1:n_fft)fty(i) = y(j+i-1) *(2*(sin((pie*i)/n_fft))*(sin((pie*i)/n_fft)))
        call realft(ftx(1:n_fft),1)     
        call realft(fty(1:n_fft),1)
        write(lun2,'(/,a1,2i10,/)')'#',j,j+n_fft
        do i=1,n_fft
         write(lun2,'(i10,es12.4,4es12.4)')i,i/float(n_fft)/2,x(j+i-1),y(j+i-1),ftx(i),fty(i)
        end do
        call find_peaks(ftx(1:n_fft), tunex)
        call find_peaks(fty(1:n_fft), tuney)
        print '(i5,a1,i5,1x,5es12.4,1x,5es12.4)',j,'-',j+n_fft,tunex(1:5)/2,tuney(1:5)/2
        write(lun3, '(i5,a1,i5,1x,5es12.4,1x,5es12.4)')j,'-',j+n_fft,tunex(1:5)/2,tuney(1:5)/2

        j=j+n_fft/2     
     end do
 end