program cbpm_anal_ac

 use precision_def
 use tbt_mod
 use cesr_read_data_mod
 use nr

 implicit none

 type (tbt_all_data_struct) tbt
 character*120 line, input_file, string
 character*100 file_name
 character*5 cnumber
 character*3 cbunch,cbpm

 integer number/-99999/
 integer ios, i,j,k,lun
 integer gain_correct
 integer dummy
 integer n
 integer highest_x(5), highest_y(5)
 integer numbers(30)/30*0/,m
 integer i_qx, i_qy
 integer nargs
 integer ix
 integer nn
 integer nlow
 integer size_ac
 integer lun1

 real(rp) rms_x, rms_y, ave_x, ave_y
 real(rp), allocatable :: x(:), y(:), fftx(:), ffty(:)
 real(rp) tune_x, tune_y
 real(rp) tunex(5), tuney(5)
 real(rp) amp_x(0:120), amp_y(0:120)
 real(rp) y_pos(4096,30), y_pos_bunch(30)

 type bpm_struct 
   real(rp) ave_x,ave_y,rms_x,rms_y
 end type

 type bunch_struct
   type(bpm_struct) bpm(0:120)
 end type
 
 type RD_file_struct
   integer number, bunch, coupling8, coupling9, chromx, chromy
   real(rp) current
 end type

 type (bunch_struct), allocatable :: bunch(:)
 type (RD_file_struct) RD_files(30)

 logical err, another/.true./

 namelist/bpm_data/numbers, tune_x, tune_y


  nargs = cesr_iargc()
  if (nargs == 1) then
    call cesr_getarg(1, input_file)
  else
    input_file = 'input.dat'
    print '(a,$)','  Input file name?  (Default = input.dat)  '
    read(5,'(a)') line
    call string_trim(line, line, ix)
    input_file = line
    if (ix==0) input_file = 'input.dat'
  endif

    open (unit=5, file= input_file, Status='old', IOSTAT=ios)
      print *,' iostat = ',ios
     read(5,nml = bpm_data, IOSTAT=ios)
      print *,' iostat = ',ios
      write(6,nml=bpm_data)
    close(5)

     lun=13
     j=8
      write(cbpm,'(i3.3)')j
      open(unit=lun, file = 'BPM_'//trim(cbpm)//'.dat',status='REPLACE')
      write(lun,'(2a12,6a12,a)')'RD','Bunch','<y>', & 
                          'rms_y', 'tune_y', 'amp_y','highest 5 fft y'

      close(unit=lun)

    do m=1,size(RD_files)
      RD_files(m)%number = numbers(m)
      if(RD_files(m)%number == 0)cycle

     number = RD_files(m)%number

     write(cnumber,'(i5)')number
     lun1 = lunget()
     open(unit=lun1, file ='RD_'//trim(cnumber)//'_ffti.dat')
     print *,' file = ','RD_'//trim(cnumber)//'_ffti.dat'

     i=0

!     read(lun1,'(a)')junk
!     print *,junk
     ios=0
     do while(ios >= 0)
      i=i+1
      read(lun1,*,IOSTAT=ios)y_pos_bunch
      y_pos(i,1:30) = y_pos_bunch(1:30)
!      print '(i10,1x,10es12.4)',ios,y_pos_bunch(1:10)
     end do
     size_ac=i
     print '(a,1x,i10,1x,a)','RD_'//trim(cnumber)//'_ffti.dat' ,i,'rows'
     close(unit=lun1)

      j=8
      write(cbpm,'(i3.3)')j
      open(unit=lun, file = 'BPM_'//trim(cbpm)//'.dat',access = 'append')
       write(lun,'(/,a1,/)')'#'
      close(lun)

  do i=1,30

    write(cbunch,'(i3.3)')i
    print *,cbunch, i
      j=8
      cbpm(1:3) = '000'
      write(cbpm,'(i3.3)')j
      file_name = 'RD_'//trim(cnumber)//'_bunch_'//trim(cbunch)//'_BPM_'//trim(cbpm)
!      lun = lunget()
      lun=13
      open(lun,file = file_name)

      ave_x=0
      ave_y=0
      rms_x=0
      rms_y=0
       n=0

      if(size_ac<1024)cycle
      allocate(y(size_ac))
      allocate(ffty(size_ac))


      do k=1,size_ac ! cycle over turns for file cnumber bunch i
        n=n+1
        ave_y = ave_y + y_pos(k,i)
      end do
      ave_y = ave_y/n

      nn=log(float(n))/log(2.)
      n=2**nn
      y(:) = y_pos(:,i)
      call realft(y(1:n),1)
      ffty(1:n) = y(1:n)
      nlow = 100. * n/4000.  !skip bins near zero
      forall(k=1:nlow)y(k)=0.
      do k=1,5 !order 5 biggest peaks in spectrum
         highest_y(k) = maxloc(y**2,1)
         y(highest_y(k))=0.
      end do
      tuney(1:5) = 1-float(highest_y(1:5))/n/2
      i_qy= minloc(abs(tune_y-tuney(1:5)),1)       


      do k=1,size_ac
        rms_y = rms_y + (y_pos(k,i)-ave_y)**2
      end do

      write(lun,'(a1,a15,i10,a15,i10,a15,i10)')'#', 'bunch = ',i, ' BPM = ',j,' nturns = ',n
      write(lun,'(a1,a15,es12.4,a15,es12.4)')'#', ' average_y =', ave_y
      write(lun,'(a1,a15,es12.4,a15,es12.4)')'#', ' rms_y =', sqrt(rms_y/n)
      write(lun,'(a1,a,5es12.4)')'#','Peaks in vertical spectrum',1-float(highest_y(1:5))/n/2
      write(lun,'(a1,a,5i10)')'#','Peaks in vertical spectrum index',highest_y(1:5)
      write(lun,'(5(a12))')'turn','y','ffty'
      do k=1,size_ac
        write(lun,'(i12,4(es12.4))')k,y_pos(k,i), ffty(k)
      end do

      print '(a,a)', trim(file_name),' written'
      close(lun)
     amp_y(j) = ffty(highest_y(i_qy))



     open(unit=lun,file = 'RD_'//trim(cnumber)//'_bunch_'//trim(cbunch)//'.dat')
     write(lun,'(5a12)')'bpm','y','rms_y'
     j=8
     write(lun,'(i12,4es12.4)')j,ave_y,rms_y
     close(lun)
     write(cbpm,'(i3.3)')j
     open(unit=lun, file = 'BPM_'//trim(cbpm)//'.dat',access = 'append')
!       write(lun,*)number,i,ave_y, rms_y,tuney(i_qy),amp_y(j), &
!               ffty(highest_y(1)),ffty(highest_y(2)),ffty(highest_y(3)),ffty(highest_y(4)),ffty(highest_y(5))
       write(lun,'(2i12,4es12.4,1x,5es12.4)')number,i,ave_y, rms_y,tuney(i_qy),amp_y(j), &
               ffty(highest_y(1)),ffty(highest_y(2)),ffty(highest_y(3)),ffty(highest_y(4)),ffty(highest_y(5))
      close(lun)
     deallocate(y,ffty)
    end do !bunch
 end do !another file

end