program cbpm_anal

 use precision_def
 use tbt_mod
 use cesr_read_data_mod
 use nr
 use cbpm_interface

 implicit none

 type (tbt_all_data_struct) tbt
 character*140 line, input_file, junk
 character*100 file_name, file_name1,file_name2
 character*6 cnumber, cnumber1, cnumber2
 character*3 cbunch,cbpm
 character*60 directory
 character*50 lattice/'No name'/

 integer number/-99999/
 integer ios, i,j,k,lun, lun_avg, lun1, lun3, lun2
 integer gain_correct
 integer dummy
 integer n
 integer highest_x(5), highest_y(5)
 integer numbers(100)/100*0/,m, bunches(100)/100*0/
 integer i_qx, i_qy
 integer nargs
 integer ix
 integer nn
 integer nlow
 integer ntune_x, ntune_y
 integer system
 integer status
 integer n_bpm, fft_window/100/
 integer save_set/0/
 integer b1(0:120), b2(0:120)

 real(rp) rms_x, rms_y, ave_x, ave_y, ave_but(1:4)
 real(rp), allocatable :: x(:), y(:), fftx(:), ffty(:)
 real(rp) tune_x, tune_y
 real(rp) tunex(5), tuney(5)
 real(rp) fit_xtune(5), fit_ytune(5)
 real(rp) amp_x(1:5,0:120), amp_y(1:5,0:120)
 real(rp) pie
 real(rp) bunch_currents(100)/100*0/
 real(rp) x1(0:120), x2(0:120), y1(0:120), y2(0:120)
 real(rp) but11(0:120),but12(0:120), but13(0:120),but14(0:120), but21(0:120),but22(0:120), but23(0:120),but24(0:120)
 real(rp) zero/0.0/

 type bpm_struct 
   real(rp) ave_x,ave_y,rms_x,rms_y, tune_x, tune_y,ave_but(1:4)
 end type

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

 type (bpm_bunch_struct), allocatable :: bunch(:)
 type (RD_file_struct) RD_files(100)

 logical err, another/.true./
 logical dir_e

 namelist/bpm_data/numbers, bunches, bunch_currents, tune_x, tune_y, lattice, save_set
     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
     lun_avg = lunget()
    open(unit=lun_avg, file='average_tunes.dat') 
    write(lun_avg,'(2a12,2a14,a12,a6)')'RD_file','bunch','tune_x_avg', 'tune_y_avg', 'current', 'bunch'

    do m=1,size(RD_files)
  
      RD_files(m)%number = numbers(m)
      RD_files(m)%bunch = bunches(m)
      RD_files(m)%bunch_current = bunch_currents(m)
      if(RD_files(m)%number == 0)cycle

!      call number_to_true_number(number, 'TBT', number,err, .true.) 
      print '(a,2i,es12.4)',' true number = ', RD_files(m)%number, RD_files(m)%bunch, RD_files(m)%bunch_current
      write(21, '(2i,es12.4)') RD_files(m)%number, RD_files(m)%bunch, RD_files(m)%bunch_current
     number = RD_files(m)%number
 
     gain_correct = gain_and_pedestal_correct$
     gain_correct = no_correct$

   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

     write(cnumber,'(i5)')number

  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+1)/2
   print *,cbunch, (tbt%bunch(i)%ix_bunch+1)/2
   if(tbt%bunch(i)%ix_bunch+1/2 /= RD_files(m)%bunch )cycle 
    do j= 0,size(tbt%bunch(i)%bpm)-1
      if(tbt%bunch(i)%bpm(j)%data_ok)then

      if(i == 1)then
         open(unit=lun, file = 'BPM_'//trim(tbt%bpm(j)%name)//'.dat',access = 'append')
         if(m == 1)write(lun,'(2a12,a12,1x,6a12,4a6,2a5)')'RD #','Bunch','<x>','<y>', & 
                         'rms_x', 'rms_y', 'tune_x' ,'tune_y', &
                          'Bunch I', 'I(ma)', 'bunch', "Q'_y", &
                         "Q'_x", 'cpl8', 'cpl9'
         if(m > 1)write(lun,'(/,a1,/)' ) '#'
        close(lun)
      endif

      file_name = trim(directory)//'/RD_'//trim(cnumber)//'_bunch_'//trim(cbunch)//'_BPM_'//trim(tbt%bpm(j)%name)
      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)))

      n=size(tbt%bunch(i)%bpm(j)%turn)
      nn=log(float(n))/log(2.)
      n=2**nn

      fft_window = abs(tune_x-tune_y)*n !tune split * number of fft points
   print '(i10,a14,1x,a,i10)', n,' points in fft',' search window = ',fft_window

      do k=1,n
        ave_x = ave_x + tbt%bunch(i)%bpm(j)%turn(k)%x
        ave_y = ave_y + tbt%bunch(i)%bpm(j)%turn(k)%y
        ave_but(1:4) = ave_but(1:4) + tbt%bunch(i)%bpm(j)%turn(k)%button(1:4)
      end do

        ave_x = ave_x/n
        ave_y = ave_y/n
        ave_but(1:4) = ave_but(1:4)/n


      do k=1,n 
       x(k) = tbt%bunch(i)%bpm(j)%turn(k)%x!-ave_x
       y(k) = tbt%bunch(i)%bpm(j)%turn(k)%y!-ave_y
       fftx(k) = x(k) *(2*(sin((pie*k)/n))*(sin((pie*k)/n)))
       ffty(k) = y(k) *(2*(sin((pie*k)/n))*(sin((pie*k)/n)))
      end do

      call realft(fftx(1:n),1)
      call realft(ffty(1:n),1)

      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-fft_window:ntune_x+fft_window)x(k)=fftx(k)
      forall(k=ntune_y-fft_window:ntune_y+fft_window)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

! fit each of the five peaks
      do k=1,5
       call fit_peak(fftx, highest_x(k),n, fit_xtune(k))
       call fit_peak(ffty, highest_y(k),n, fit_ytune(k))
     end do

      tunex(1:5) = 1-float(highest_x(1:5))/n/2
      tuney(1:5) = 1-float(highest_y(1:5))/n/2
      fit_xtune(1:5) = 1-fit_xtune(1:5)/2
      fit_ytune(1:5) = 1-fit_ytune(1:5)/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_peak','ffty_peak'

      do k=1,n !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+1)/2, ' 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)')'#','Fitted Peaks in horizontal spectrum', fit_xtune(1:5)
      write(lun,'(a1,a,5es12.4)')'#','Peaks in vertical spectrum',1-float(highest_y(1:5))/n/2
      write(lun,'(a1,a,5es12.4)')'#','Fitted Peaks in vertical spectrum',fit_ytune(1:5)
      do k=1,n !size(tbt%bunch(i)%bpm(j)%turn)
        write(lun,'(i12,4(es14.6))')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)%ave_but(1:4) = ave_but(1:4)
     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))
     bunch(i)%bpm(j)%tune_x = tunex(1)
     bunch(i)%bpm(j)%tune_y = tuney(1)
     bunch(i)%bpm(j)%tune_x = fit_xtune(1)
     bunch(i)%bpm(j)%tune_y = fit_ytune(1)
     deallocate(x,y,fftx,ffty)

     endif
     end do !bpm

     open(unit=lun,file = 'RD_'//trim(cnumber)//'/RD_'//trim(cnumber)//'_bunch_'//trim(cbunch)//'.dat')
     lun1 = lunget()
     open(unit=lun1,file = 'RD_'//trim(cnumber)//'/orbit'//'_bunch_'//trim(cbunch)//'.dat')
     write(lun,'(5a12)')'bpm','x','y','rms_x','rms_y'
      write(lun1, '(a16)')'&DATA_PARAMETERS'
      write(lun1, '(a,a)')'file_type = ',"'ORBIT DATA'"
      write(lun1, '(a,a)')'lattice = ',"'"//trim(lattice)//"'"
      write(lun1, '(a,a)')'data_date = ',"'"//tbt%time_stamp//"'"
      write(lun1, '(a,i10)')'save_set = ',save_set
      write(lun1,'(a4)')'/END'
      write(lun1, '(a,es16.8)')'total_current = ',tbt%total_current
      write(lun1,'(a)')'        DET      |X_orbit [mm]|      Y_orbit[mm]    |        but 1               but 2              but 3              but 4     '
     n_bpm=0
     bunch(i)%tune_x_avg = 0
     bunch(i)%tune_y_avg = 0
     do j=1,99
       if(.not. tbt%bunch(1)%bpm(j)%data_ok)write(lun1,'(i10,6es20.6)')j,zero, zero, zero, zero, zero, zero
       if(tbt%bunch(i)%bpm(j)%data_ok)then
         write(lun1,'(i10,6es20.6)')j,bunch(i)%bpm(j)%ave_x*1000, bunch(i)%bpm(j)%ave_y*1000, bunch(i)%bpm(j)%ave_but(1:4)
         n_bpm = n_bpm+1
         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
         bunch(i)%tune_x_avg = bunch(i)%tune_x_avg + bunch(i)%bpm(j)%tune_x
         bunch(i)%tune_y_avg = bunch(i)%tune_y_avg + bunch(i)%bpm(j)%tune_y
       endif
     end do
     close(lun)
     close(lun1)
     do j=1,99
      if(tbt%bunch(i)%bpm(j)%data_ok )then
      open(unit=lun, file = 'BPM_'//trim(tbt%bpm(j)%name)//'.dat',access = 'append')
       write(lun,'(2i12,6es12.4,1x,es12.4,f6.2,3i6,2i5)')number,(tbt%bunch(i)%ix_bunch+1)/2,bunch(i)%bpm(j)%ave_x,&
                         bunch(i)%bpm(j)%ave_y, & 
                         bunch(i)%bpm(j)%rms_x, &
                         bunch(i)%bpm(j)%rms_y, bunch(i)%bpm(j)%tune_x,bunch(i)%bpm(j)%tune_y,RD_files(m)%bunch_current, &
                         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)
      endif
     end do

     write(lun_avg,'(2i12,2es14.6,f12.3,i6)')number,(tbt%bunch(i)%ix_bunch+1)/2, bunch(i)%tune_x_avg/n_bpm, bunch(i)%tune_y_avg/n_bpm,&
                                            RD_files(m)%current, RD_files(m)%bunch
    end do !bunch
  deallocate(bunch)
  write(lun_avg,'(/,a1,/)')'#'
 end do !another RD file

! combine orbit files for first and second half
     
    do m=1,size(RD_files),2
      if(RD_files(m)%number == 0)cycle
      write(cnumber1,'(i5)')RD_files(m)%number
      write(cnumber2,'(i5)')RD_files(m+1)%number
      do i=1,1 !size(tbt%bunch)
      write(cbunch,'(i3.3)')(tbt%bunch(i)%ix_bunch+1)/2
      if(tbt%bunch(i)%ix_bunch+1/2 /= RD_files(m)%bunch )cycle 
      lun1=lunget()
      file_name1 = 'RD_'//trim(cnumber1)//'/orbit'//'_bunch_'//trim(cbunch)//'.dat'
      open(unit=lun1,file = file_name1)
      lun2=lunget()
      file_name2 = 'RD_'//trim(cnumber2)//'/orbit'//'_bunch_'//trim(cbunch)//'.dat'
      open(unit=lun2,file = file_name2)

      lun3=lunget()
      file_name = 'RD_'//trim(cnumber1)//'_orbit'//'_bunch_'//trim(cbunch)//'.dat'
      open(unit=lun3,file = file_name)
      do j=1,8
       read(lun1,'(a)')junk
       write(lun3,'(a)')junk
       read(lun2, '(a)')junk
      end do
      do j=1,99
       x1(j)=0.
       x2(j)=0.
       y1(j)=0
       y2(j)=0
       read(lun1,'(i10,6es20.6)')b1(j),x1(j),y1(j), but11(j), but12(j), but13(j), but14(j)       
       read(lun2,'(i10,6es20.6)')b2(j),x2(j),y2(j), but21(j), but22(j), but23(j), but24(j)       
       write(lun3,'(i10,6es20.6)')b1(j),x1(j)+x2(j), y1(j)+y2(j), but11(j)+but21(j), but12(j)+but22(j), but13(j)+but23(j), but14(j)+but24(j)
      end do       
      close(unit=lun1)
      close(unit=lun2)
      close(unit=lun3)
     end do
   end do

end