program multiparticle

   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
      type (beam_struct) 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 k,l


       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
       
       character*3 num

      character*140 lat_file
      character*120 line, last_line
      logical error/.false./, excitation_and_damping/.false./,quad_offsets/.false./, absolute_feedback/.false./
      namelist /beam_def/nturns, n_bunch, bunch_spacing, n_particle, bunch_charge, &
                         excitation_and_damping, initial_offset, quad_offsets, y_off,cutoff,seed
      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='mparticle_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)
       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)
       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

       print '(a)',lat_file
       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

       call twiss_at_start(lat)
       call calc_z_tune(lat)

       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

       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)

        forall(i=1:6)beam%bunch(1)%particle(:)%vec(i) = orbit(0)%vec(i)+beam%bunch(1)%particle(:)%vec(i) + initial_offset(i)

        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
              call track_beam(lat,beam,ele1=lat%ele(j),ele2=lat%ele(j+1), err=error)  
             end do

! to write out the position of every particle in bunch k at the end of  turn i
             do k=1, n_bunch
               write(21, '(/,a,i5,a,i5)')' Turn = ', i, '          Bunch = ',k
               do m=1, n_particle
                write(21,'(6es12.4)')beam%bunch(k)%particle(m)%vec(1:6)
               end do
              end do

        end do
end