program cbpm_anal

 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
 character*100 file_name
 character*6 cnumber
 character*3 cbunch,cbpm
 character*60 directory

 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 ntune_x, ntune_y
 integer system
 integer status

 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(1:5,0:120), amp_y(1:5,0:120)
 real(rp) pie

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

 type bpm_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 (bpm_bunch_struct), allocatable :: bunch(:)
 type (RD_file_struct) RD_files(30)

 logical err, another/.true./
 logical dir_e

 namelist/bpm_data/numbers, tune_x, tune_y
     pie = atan(1.)*4

  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
     do j=1,99
      write(cbpm,'(i3.3)')j
      open(unit=lun, file = 'BPM_'//trim(cbpm)//'.dat',status='REPLACE')
      write(lun,'(2a12,6a12,1x,6a12,1x,6a12,4a6,2a5)')'RD #','Bunch','<x>','<y>', & 
                         'rms_x', 'rms_y', 'tune_x' ,'tune_y', &
                          'amp_x','amp_x(1)','amp_x(2)','amp_x(3)','amp_x(4)','amp_x(5)', 'amp_y', &
                          'amp_y(1)','amp_y(2)','amp_y(3)','amp_y(4)','amp_y(5)', 'I(ma)', 'bunch', "Q'_y", &
                         "Q'_x", 'cpl8', 'cpl9'

      close(unit=lun)
    end do

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

!      call number_to_true_number(number, 'TBT', number,err, .true.) 
      print '(a,i)',' true number = ', RD_files(m)%number
!      print '(a,$)',' file number = '
!       accept *,number
     number = RD_files(m)%number
 
     gain_correct = gain_and_pedestal_correct$
!     print '(a)',' Correct tbt data for gain and pedestal ? '
!     print '(a,$)',' (1=no_correct, 2=pedestal correct, 3=gain and pedestal correct ) '
!     read(*,'(i)')dummy
!     select case(dummy)
!      case(1)
!        gain_correct = no_correct$
!      case(2)
!        gain_correct = pedestal_correct$
!      case(3)
!        gain_correct = gain_and_pedestal_correct$
!      end select
!     print '(a,i)',' gain correct = ', dummy

   call tbt_read_data_file(number, gain_correct, .true., tbt, err)
   print '(a,i7)',' Success reading file ',number

     if(err)then
         print '(a,i7)',' error reading file ', number
         stop
     endif

   call tbt_read_data_file(number+1, gain_correct, .false., tbt, err)
   print '(a,i7)',' Success reading  file ',number+1

     if(err)then
         print '(a,i7)',' error reading file ', number
         stop
     endif

     write(cnumber,'(i5)')number

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

  allocate(bunch(size(tbt%bunch)))

   directory = trim('RD_'//trim(cnumber))
    status= system('mkdir ' // directory)
    if(status == 0)print '(a20,a)',' create directory = ', trim(directory)

  RD_files(m)%current = tbt%total_current

  do i=1,size(tbt%bunch)

    write(cbunch,'(i3.3)')tbt%bunch(i)%ix_bunch
   print *,cbunch, tbt%bunch(i)%ix_bunch
    do j= 1,size(tbt%bunch(i)%bpm)
      if(j > 99) exit
      cbpm(1:3) = '000'
      write(cbpm,'(i3.3)')j
      file_name = trim(directory)//'/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(tbt%bunch(i)%bpm(j)%turn) < 1024)cycle
      allocate(x(size(tbt%bunch(i)%bpm(j)%turn)))
      allocate(y(size(tbt%bunch(i)%bpm(j)%turn)))
      allocate(fftx(size(tbt%bunch(i)%bpm(j)%turn)))
      allocate(ffty(size(tbt%bunch(i)%bpm(j)%turn)))

!      print '(3(a15,i10))','Bunch = ', tbt%bunch(i)%ix_bunch,  ' bpm = ', j, &
!               ' turns = ', size(tbt%bunch(i)%bpm(j)%turn)

      do k=1,size(tbt%bunch(i)%bpm(j)%turn)
        n=n+1
        ave_x = ave_x + tbt%bunch(i)%bpm(j)%turn(k)%x
        ave_y = ave_y + tbt%bunch(i)%bpm(j)%turn(k)%y
      end do
        ave_x = ave_x/n
        ave_y = ave_y/n

      nn=log(float(n))/log(2.)
      n=2**nn
      x(:) = tbt%bunch(i)%bpm(j)%turn(:)%x-ave_x
      y(:) = tbt%bunch(i)%bpm(j)%turn(:)%y-ave_y

      forall(i=1:n)x(i) = x(i) *(2*(sin((pie*i)/n))*(sin((pie*i)/n)))
      forall(i=1:n)y(i) = y(i) *(2*(sin((pie*i)/n))*(sin((pie*i)/n)))

      call realft(x(1:n),1)
      call realft(y(1:n),1)
      fftx(1:n) = x(1:n)
      ffty(1:n) = y(1:n)
      nlow = 100. * n/4000.  !skip bins near zero
      ntune_x = (1.-tune_x)*n*2
      ntune_y = (1.-tune_y)*n*2
      x=0.
      y=0.
      forall(k=ntune_x-100:ntune_x+100)x(k)=fftx(k)
      forall(k=ntune_y-100:ntune_y+100)y(k)=ffty(k)
      do k=1,5 !order 5 biggest peaks in spectrum
         highest_x(k) = maxloc(x**2,1)
         x(highest_x(k))=0.
         highest_y(k) = maxloc(y**2,1)
         y(highest_y(k))=0.
      end do
      tunex(1:5) = 1-float(highest_x(1:5))/n/2
      tuney(1:5) = 1-float(highest_y(1:5))/n/2
      i_qx= minloc(abs(tune_x-tunex(1:5)),1)        
      i_qy= minloc(abs(tune_y-tuney(1:5)),1)       

      write(lun,'(5(a12))')'turn','x','y','fftx','ffty'

      do k=1,size(tbt%bunch(i)%bpm(j)%turn)
        rms_x = rms_x + (tbt%bunch(i)%bpm(j)%turn(k)%x-ave_x)**2
        rms_y = rms_y + (tbt%bunch(i)%bpm(j)%turn(k)%y-ave_y)**2
      end do

      write(lun,'(a1,a15,i10,a15,i10,a15,i10)')'#', 'bunch = ',tbt%bunch(i)%ix_bunch, ' BPM = ',j,' nturns = ',n
      write(lun,'(a1,a15,f6.2,a11,i6,a11,i6,a15,i4,a15,i4)')' #','bunch current = ',RD_files(m)%current,' xqune1 = ',RD_files(m)%chromy, &
               ' xqune2 = ',RD_files(m)%chromx, &
               ' coupling 8 = ',RD_files(m)%coupling8,' coupling 9 = ',RD_files(m)%coupling9
      write(lun,'(a1,a15,es12.4,a15,es12.4)')'#', 'average_x = ',ave_x,' average_y =', ave_y
      write(lun,'(a1,a15,es12.4,a15,es12.4)')'#', 'rms_x = ',sqrt(rms_x/n),' rms_y =', sqrt(rms_y/n)
      write(lun,'(a1,a,5es12.4)')'#','Peaks in horizontal spectrum', 1-float(highest_x(1:5))/n/2
      write(lun,'(a1,a,5es12.4)')'#','Peaks in vertical spectrum',1-float(highest_y(1:5))/n/2
      do k=1,size(tbt%bunch(i)%bpm(j)%turn)
        write(lun,'(i12,4(es12.4))')k,tbt%bunch(i)%bpm(j)%turn(k)%x,tbt%bunch(i)%bpm(j)%turn(k)%y, &
                                 fftx(k), ffty(k)
      end do

      print '(a,a)', trim(file_name),' written'
      close(lun)
     bunch(i)%bpm(j)%ave_x = ave_x
     bunch(i)%bpm(j)%ave_y = ave_y
     bunch(i)%bpm(j)%rms_x = sqrt(rms_x/n)
     bunch(i)%bpm(j)%rms_y = sqrt(rms_y/n)
     amp_x(1:5,j) = fftx(highest_x(1:5))
     amp_y(1:5,j) = ffty(highest_y(1:5))
     deallocate(x,y,fftx,ffty)


     end do !bpm

     open(unit=lun,file = 'RD_'//trim(cnumber)//'/RD_'//trim(cnumber)//'_bunch_'//trim(cbunch)//'.dat')
     write(lun,'(5a12)')'bpm','x','y','rms_x','rms_y'
     do j=1,99
       write(lun,'(i12,4es12.4)')j,bunch(i)%bpm(j)%ave_x,bunch(i)%bpm(j)%ave_y,bunch(i)%bpm(j)%rms_x,bunch(i)%bpm(j)%rms_y
     end do
     close(lun)
     do j=1,99
      write(cbpm,'(i3.3)')j
      open(unit=lun, file = 'BPM_'//trim(cbpm)//'.dat',access = 'append')
       write(lun,'(2i12,6es12.4,1x,6es12.4,1x,6es12.4,f6.2,3i6,2i5)')number,tbt%bunch(i)%ix_bunch,bunch(i)%bpm(j)%ave_x,&
                         bunch(i)%bpm(j)%ave_y, & 
                         bunch(i)%bpm(j)%rms_x, &
                         bunch(i)%bpm(j)%rms_y, tunex(i_qx),tuney(i_qy),amp_x(i_qx,j), amp_x(1:5,j), amp_y(i_qy,j),amp_y(1:5,j), &
                         RD_files(m)%current, RD_files(m)%bunch, RD_files(m)%chromy, &
                         RD_files(m)%chromx, RD_files(m)%coupling8, RD_files(m)%coupling9
      close(lun)
     end do

    end do !bunch
  deallocate(bunch)
 end do !another

end