program multibunch

   use bmad
   use beam_mod
   use mode3_mod
   use multibunch_mod
   use multibunch_interface
   use bmadz_other_code_interface   

     implicit none

      type (lat_struct) lat
      type (ele_struct) ele 
      type (coord_struct), allocatable :: orbit(:), co(:)
      type (beam_init_struct) beam_init, fft_beam_init
      type (beam_struct) beam, fft_beam
      type (normal_modes_struct) mode
      type (rad_int_all_ele_struct) rad_int
      type (moment_struct), allocatable :: bunch_moments(:)
      
       integer unit, number, m, i, ix, nargs,j
       integer ix_cache/0/
       integer lun
       integer nturns, n_bunch, n_particle
       integer ios
       integer mturn/1/
       integer seed/1/
       integer bunch
       integer nb
       integer k,l
       integer, allocatable:: n(:)
       integer n_fft/2048/

       real(rp) bunch_spacing, bunch_charge
       real(rp) density/1.e12/, sigx/0.01/, sigy/0.002/, sigz/0.01/
       real(rp) initial_offset(1:6), y_off, cutoff
       real(rp), allocatable:: sum(:,:),ave(:,:) 
       real(rp) Q(3),phase(3)
       real(rp) qx_before, qy_before, qx_after, qy_after
       real(rp) bbi_length, d_bbi
       real(rp) emitx0/0./, emity0/0./
       
       character*29 beamfield_file
       character*3 num
       character*3 ecloud_type

      character*140 lat_file
      character*120 line, last_line
      logical error/.false./, ecloud/.false./,excitation_and_damping/.false./,quad_offsets/.false./, absolute_feedback/.false./
      namelist /beam_def/nturns, n_bunch, bunch_spacing, n_particle, bunch_charge, &
                         ecloud, density, sigx, sigy, sigz, excitation_and_damping, mturn, &
                         initial_offset, quad_offsets, y_off,cutoff,seed, absolute_feedback, bunch, ecloud_type, beamfield_file, &
                         bbi_length, d_bbi, emitx0, emity0

      bmad_com%radiation_damping_on = .false.
      bmad_com%radiation_fluctuations_on = .false.

      nargs = cesr_iargc()
      if(nargs == 1)then
         call cesr_getarg(1,lat_file)
         print *, 'Using ', trim(lat_file)
       else
        lat_file = 'bmad.'
        print '(a,$)',' Lattice file name ? (default= bmad.) '
        read(5,'(a)') line
        call string_trim(line, line, ix)
        lat_file = line
        if(ix == 0) lat_file = 'bmad.'
        print *, ' lat_file = ', lat_file
      endif

     lun = lunget()
     open(unit=lun, file='beam_def.in', STATUS ='old')
     read(lun, nml=beam_def, IOSTAT=ios)
     rewind(lun)
     read(lun, nml=beam_def, IOSTAT=ios)
     close(unit=lun)
     write(6,nml=beam_def)

       open(unit=13,file='bunch.dat')

       call bmad_parser (lat_file, lat)
       print *,' bmad_parser done'
       call reallocate_coord (orbit, lat%n_ele_track)
       call reallocate_coord (co, n_fft)
       if(quad_offsets)call implement_quad_offsets(lat,y_off,cutoff,seed)
       call closed_orbit_calc(lat,orbit,6)
       call track_all(lat,orbit)
       call lat_make_mat6(lat, -1, ref_orb = orbit)
       do i=1,lat%n_ele_track
          write(14,'(i10,1x,7es12.4)')i,lat%ele(i)%s, orbit(i)%vec
       end do
       call twiss_at_start(lat)
       call twiss_propagate_all(lat)
       call track_all(lat,orbit)
       call radiation_integrals (lat, orbit, mode, ix_cache, 0, rad_int)
       if(emitx0 /= 0)mode%a%emittance = emitx0
       if(emity0 /= 0)mode%b%emittance = emity0
       call calc_z_tune(lat)

       write(15,'(8a12)')'s','x','y','de','betax','betay','etax','etay'
       do i=1,lat%n_ele_track
        write(15,'(8es12.4)')lat%ele(i)%s,orbit(i)%vec(1),orbit(i)%vec(3),orbit(i)%vec(6),lat%ele(i)%a%beta,&
                     lat%ele(i)%b%beta, lat%ele(i)%x%eta, lat%ele(i)%y%eta
       end do

       write(13,'(a)')lat_file
       print '(a)',lat_file
       print *,' No ecloud'
       print '(3(a,es12.4))','horizontal tune = ', lat%a%tune/twopi,' vertical tune = ', lat%b%tune/twopi,' synch tune = ',lat%z%tune/twopi
       print '(4(a,es12.4))',' beta_x = ',lat%ele(lat%n_ele_track)%a%beta,' alpha_x = ',lat%ele(lat%n_ele_track)%a%alpha, &
               ' beta_y = ',lat%ele(lat%n_ele_track)%b%beta,' alpha_y = ',lat%ele(lat%n_ele_track)%b%alpha
       write(13,'(a)')' No ecloud'
       qx_before = lat%a%tune/twopi
       qy_before = lat%b%tune/twopi
       write(13, '(3(a,es12.4))')'horizontal tune = ', lat%a%tune/twopi,' vertical tune = ', lat%b%tune/twopi,' synch tune = ',lat%z%tune/twopi
       if(ecloud)call implement_ecloud(lat, sigx,sigy,sigz, density, bunch, ecloud_type, beamfield_file, bbi_length, d_bbi)
       call twiss_at_start(lat)
       call calc_z_tune(lat)

       print *
       print *,' With ecloud'
       print '(3(a,es12.4))','horizontal tune = ', lat%a%tune/twopi,' vertical tune = ', lat%b%tune/twopi,' synch tune = ',lat%z%tune/twopi
       write(13,'(a)')' With ecloud'
       write(13, '(3(a,es12.4))')'horizontal tune = ', lat%a%tune/twopi,' vertical tune = ', lat%b%tune/twopi,' synch tune = ',lat%z%tune/twopi
       qx_after = lat%a%tune/twopi
       qy_after = lat%b%tune/twopi
       write(13,'(2a,a,i10)')'ecloud_type =',ecloud_type,' bunch ',bunch 

       beam_init%n_bunch = n_bunch
       beam_init%dt_bunch = bunch_spacing
       beam_init%n_particle=n_particle
       beam_init%bunch_charge = bunch_charge

       beam_init%a_emit = mode%a%emittance
       beam_init%b_emit = mode%b%emittance
       beam_init%sig_z = mode%sig_z
       beam_init%sig_e = mode%sigE_E

       fft_beam_init%n_bunch = n_bunch
       fft_beam_init%dt_bunch = bunch_spacing
       fft_beam_init%n_particle = n_fft
       fft_beam_init%bunch_charge = bunch_charge

       fft_beam_init%a_emit = mode%a%emittance
       fft_beam_init%b_emit = mode%b%emittance
       fft_beam_init%sig_z = mode%sig_z
       fft_beam_init%sig_e = mode%sigE_E

       write(13, '(2(a23,es12.4))') 'horizontal emittance = ',beam_init%a_emit,'vertical emittance = ', beam_init%b_emit
       write(13, '(2(a23,es12.4))') 'bunch length = ',beam_init%sig_z,'sig E/E = ', beam_init%sig_e
       write(13, '(a23,es12.4)')'bunch charge =', beam_init%bunch_charge
       write(13, '(1(a23,i10))')'number of bunches =', beam_init%n_bunch
       write(13, '(1(a23,i10))') 'particles/bunch =', beam_init%n_particle
       write(13, '(1(a23,i10))') 'number of turns =', nturns
       write(13, '(a,6es12.4)')'initial_offset =',initial_offset(1:6)

       print '(2(a23,es12.4))', 'horizontal emittance = ',beam_init%a_emit,'vertical emittance = ', beam_init%b_emit
       print '(2(a23,es12.4))', 'bunch length = ',beam_init%sig_z,'sig E/E = ', beam_init%sig_e
       print '(a23,es12.4)','bunch charge =', beam_init%bunch_charge
       print '(1(a23,i10))','number of bunches =', beam_init%n_bunch
       print '(1(a23,i10))', 'particles/bunch =', beam_init%n_particle
       print '(1(a23,i10))', 'number of turns =', nturns
       print '(a,6es12.4)','initial_offset =',initial_offset(1:6)

       call init_beam_distribution(lat%ele(0),lat%param, beam_init, beam)
       call init_beam_distribution(lat%ele(0),lat%param, fft_beam_init, fft_beam)
       allocate(bunch_moments(0:beam_init%n_bunch))
       allocate(ave(1:6,1:beam_init%n_bunch), n(1:beam_init%n_bunch))
       allocate(sum(1:6,1:beam_init%n_bunch))
        forall(i=1:6)beam%bunch(1)%particle(:)%vec(i) = orbit(0)%vec(i)+beam%bunch(1)%particle(:)%vec(i) + initial_offset(i)

        !        write(13, '(a10,10a12)')'Turn','Bunch01 x','Bunch02 x','Bunch03 x','Bunch04 x','Bunch05 x','Bunch06 x','Bunch07 x','Bunch08 x','Bunch09 x','Bunch10 x'
        write(13, '(2a10,6a12)')'Turn','Bunch','ave_x','ave_y','ave_z','sigx','sigy','sigz'

        if(excitation_and_damping)then
           bmad_com%radiation_damping_on = .true.
           bmad_com%radiation_fluctuations_on = .true.
           print '(a,l)',' radiation_damping_on = ',bmad_com%radiation_damping_on
           print '(a,l)',' radiation_fluctuations_on = ',bmad_com%radiation_fluctuations_on 
        endif
        ix=0
        do i=1,nturns
           do j = 0,lat%n_ele_track-1
             if(index(lat%ele(j+1)%name, 'ECLOUD')/= 0)then !compute average vec(1:6)
               do nb=1,n_bunch
                n(nb)=0
                sum(1:6,nb)=0
                 do l = 1, size(beam%bunch(nb)%particle) !count number of particles
                  if(beam%bunch(nb)%particle(l)%state /= alive$)cycle             
                   n(nb)=n(nb)+1
                   forall(k=1:6)sum(k,nb) =  beam%bunch(nb)%particle(l)%vec(k) + sum(k,nb) !average
                 end do
                 ave(1:6,nb) = sum(1:6,nb)/n(nb)
                 forall(k=1:6)beam%bunch(nb)%particle(:)%vec(k) = beam%bunch(nb)%particle(:)%vec(k) - ave(k,nb)
                 write(18,'(i10,es12.4,1x,i10,1x,6es12.4)')nb,(i-1)*lat%ele(lat%n_ele_track)%s+lat%ele(j+1)%s, n(nb), ave(1:6,nb)
               end do

              call track_beam(lat,beam,ele1=lat%ele(j),ele2=lat%ele(j+1), err=error)  
              do nb=1,n_bunch
               forall(k=1:6)beam%bunch(nb)%particle(:)%vec(k) = beam%bunch(nb)%particle(:)%vec(k) + ave(k,nb)
              end do
             else
              call track_beam(lat,beam,ele1=lat%ele(j),ele2=lat%ele(j+1), err=error)  
             endif
           end do
           if(absolute_feedback)call feedback(beam, bunch_moments)
           call beam_moments(beam, bunch_moments)
           if(mturn *(i/mturn) == i)call write_phase_space(i,beam)
!           print '(a,i10)',' Turn = ', i
!            write(13, '(i10,10es12.4)')i,(beam%bunch(j)%particle(1)%vec(1), j=1,n_bunch)
           write(13, '(2i10,3es12.4,3es12.4,1x,i10)')(i,j,bunch_moments(j)%ave(1),bunch_moments(j)%ave(3),bunch_moments(j)%ave(5), bunch_moments(j)%sigma(1,1), &
                       bunch_moments(j)%sigma(3,3),bunch_moments(j)%sigma(5,5),bunch_moments(j)%n_alive,j=1,n_bunch)
           write(6,  '(2i10,3es12.4,3es12.4,1x,i10)')(i,j,bunch_moments(j)%ave(1),bunch_moments(j)%ave(3),bunch_moments(j)%ave(5), bunch_moments(j)%sigma(1,1), &
                bunch_moments(j)%sigma(3,3),bunch_moments(j)%sigma(5,5),bunch_moments(j)%n_alive,j=1,n_bunch)
           if(nturns - i < n_fft) then
            ix=ix+1
            forall(j=1:n_bunch)fft_beam%bunch(j)%particle(ix)%vec = bunch_moments(j)%ave
           endif
        end do

        do j=1,n_bunch
         forall(ix=1:n_fft)co(ix)%vec = fft_beam%bunch(j)%particle(ix)%vec
         call fft_plane(n_fft,co,Q,phase)
         write(13,'(i10, a,3es12.4,a,3es12.4)')j, ' Q = ', Q, ' phase = ', phase
         print '(i10, a,3es12.4,a,3es12.4)', j, ' Q = ', Q, ' phase = ', phase
        end do
        write(13,'(a)')lat_file
        write(13,'(a)')' No ecloud'
        write(13, '(3(a,es12.4))')'horizontal tune = ', qx_before,' vertical tune = ', qx_before,' synch tune = ',lat%z%tune/twopi

        write(13,'(a)')' With ecloud'
        write(13, '(3(a,es12.4))')'horizontal tune = ', qx_after,' vertical tune = ', qx_after,' synch tune = ',lat%z%tune/twopi
        write(13,'(2a,a,i10)')'ecloud_type =',ecloud_type,' bunch ',bunch 

        close(13)
end