program multibunch
  use bmad
  use beam_mod
  use mode3_mod
  implicit none

  type (lat_struct) lat
  type (ele_struct) ele 
  type (coord_struct), allocatable :: orbit(:)
  type (beam_init_struct) beam_init
  type (beam_struct) beam
  type (normal_modes_struct) mode
  type (rad_int_all_ele_struct) rad_int

  integer unit, number, m, i, ix, nargs,j
  integer ix_cache/0/
  integer lun
  integer nturns, n_bunch, n_particle
  integer ios
  integer multiuse_counter

  real(rp) bunch_spacing, bunch_charge

  character*19 file_name
  character*3 num

  character*140 lat_file
  character*120 line, last_line
  logical error/.false./
  namelist /beam_def/nturns, n_bunch, bunch_spacing, n_particle, bunch_charge

  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)
  close(unit=lun)
  write(6,nml=beam_def)

  call bmad_parser (lat_file, lat)

  call reallocate_coord (orbit, lat%n_ele_track)
  call lat_make_mat6(lat, -1)
  call twiss_at_start(lat)
  call twiss_propagate_all(lat)
  call closed_orbit_calc(lat,orbit,6)
  orbit(0)%vec(:) = 0
  call track_all(lat,orbit)
  call radiation_integrals (lat, orbit, mode, ix_cache, 0, rad_int)


  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%a_emit = 2.e-9
  beam_init%b_emit = mode%b%emittance !beam_init%b_emit = 10.e-12
  beam_init%sig_z = mode%sig_z        !beam_init%sig_z = 0.01
  beam_init%sig_e = mode%sigE_E       !beam_init%sig_e = 0.0008


  print *, "Tune: ", lat%a%tune, lat%b%tune
  print *, "Orbit:", orbit(0)%vec

  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

  call init_beam_distribution(lat%ele(0),lat%param, beam_init, beam)
  print *, "Beam information:"
  print *, "     " ,"bunch"," ", "t_center", "# particles"
  
  do i=1, size(beam%bunch)
    beam%bunch(i)%particle(:)%vec(3) = 0.001
  end do

  do i=1,nturns
    call write_turn_data(beam_init,beam,i)
    call track_beam(lat,beam,ele1=lat%ele(0),ele2=lat%ele(lat%n_ele_track), err=error)  
  end do

end program