program doublebunch
 
  use bmad
  use beam_mod
  use mode3_mod
  implicit none

  type (lat_struct) lat

  type (ele_struct) ele 
  type (ele_struct), pointer :: sel_ele

  type (coord_struct), allocatable :: orbit(:)

  type (beam_init_struct) beam_init
  type (beam_struct) beam
!  type (bunch_params_struct) bunch_params

  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 n_bunch, n_particle, observer_index
  integer ios

  real(rp) bunch_spacing, bunch_charge, leader_x_off, leader_y_off, leader_z_off, leader_charge
!  real(rp) sigma_s(6,6), s(6,6), charge(1:beam_init%n_particle)

  character*19 file_name
  character*3 num

  character*140 lat_file
  character*120 line, last_line
  logical error/.false./
  namelist /doublebunch_nml/n_bunch, bunch_spacing, n_particle, bunch_charge, leader_x_off, leader_y_off, leader_z_off, leader_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='doublebunch_nml.in', STATUS ='old')
  read(lun, nml=doublebunch_nml, IOSTAT=ios)
  close(unit=lun)
  write(6,nml=doublebunch_nml)

  call bmad_parser (lat_file, lat)
  call reallocate_coord (orbit, lat%n_ele_track)

  if (lat%param%geometry==2) then
    print *, "[+] CLOSED GEOMETRY: Calculating radiation integrals"
    call lat_make_mat6(lat, -1)
    call twiss_at_start(lat)
    call twiss_propagate_all(lat)
    call closed_orbit_calc(lat,orbit,6)
    call track_all(lat,orbit)
    call radiation_integrals (lat, orbit, mode, ix_cache, 0, rad_int)
    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
  else
    print *, "[+] OPEN GEOMETRY: Using pre-assigned Twiss values"
    beam_init%a_emit = 2.e-9
    beam_init%b_emit = 10.e-12
    beam_init%sig_z = 0.01
    beam_init%sig_e = 0.0008
  end if

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

  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 *, 'Initializating beam distribution'
  call init_beam_distribution(lat%ele(0),lat%param, beam_init, beam)
  
  ! Change leader charge
  beam%bunch(1)%charge_tot = leader_charge
  ! Offset the leader
  beam%bunch(1)%particle(:)%vec(1) = beam%bunch(1)%particle(:)%vec(1) + leader_x_off
  beam%bunch(1)%particle(:)%vec(2) = beam%bunch(1)%particle(:)%vec(2) + leader_y_off
  beam%bunch(1)%particle(:)%vec(3) = beam%bunch(1)%particle(:)%vec(3) + leader_z_off

  print *, "Elements:"

  print *, ' Ix  Name              Ele_type                                      S          Beta_a'
  do i=0, lat%n_ele_track
    sel_ele => lat%ele(i)
    print '(i4,2x,a16,2x,a,2f12.4)', i, sel_ele%name, key_name(sel_ele%key), sel_ele%s, sel_ele%a%beta
  end do

  print *, 'Starting element-by-element tracking'
  do i=0,lat%n_ele_track-1
    sel_ele => lat%ele(i+1)
    print *, 'Tracking ele:',sel_ele%name, key_name(sel_ele%key)!,i,'-->',(i+1)
    call track_beam(lat,beam,ele1=lat%ele(i),ele2=lat%ele(i+1), err=error)
    call write_turn_data(beam_init,beam,i)
  end do
end program