program doublebunch
 
  use bmad
  use beam_mod
  use mode3_mod
  implicit none

  !================================================!
  ! Data structure to hold elements of the lattice !
  ! that use custom field calculations             !
  !================================================!

  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
  integer multiuse_counter

  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

  type ListElem
    type(ele_struct) :: value;
    type(ListElem), pointer :: next;
  end type ListElem

  type(ListElem), pointer :: head
  type(ListElem), pointer :: temp_pointer
  nullify(head)
  nullify(temp_pointer)

  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"
    !===============================!
    ! Now remove custom field calcs !
    !===============================!
    multiuse_counter = 0
    do i=0, lat%n_ele_track-1
      print *, i, lat%ele(i)%field_calc
      if (lat%ele(i)%field_calc == custom$) then
        if (associated(head)) then
          allocate(temp_pointer)
          temp_pointer%value = lat%ele(i)
          temp_pointer%next => head
          head => temp_pointer
        else
          allocate(head)
          head%value = lat%ele(i)
        end if
        lat%ele(i)%field_calc = bmad_standard$
        multiuse_counter = multiuse_counter + 1
      endif
    end do
    print *, "Removed cloud from", multiuse_counter, "elements"
    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
    !========================!
    ! Reenable custom fields !
    !========================!
    multiuse_counter = 0
    temp_pointer => head
    do while (associated(temp_pointer))
      temp_pointer%value%field_calc = custom$
      temp_pointer => temp_pointer%next
      multiuse_counter = multiuse_counter + 1
    end do
    print *, "Added cloud to", multiuse_counter, "elements"
  else
    print *, "[+] OPEN GEOMETRY: Using pre-assigned Twiss values"
    !values from backup_periodic
    beam_init%a_emit = 6.4246E-09
    beam_init%b_emit = 1.1612E-14
    beam_init%sig_z = 3.1346E-03
    beam_init%sig_e = 1.6544E-04
  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 *, " Tune:",lat%a%tune, lat%b%tune
  open(unit=13, file="tunes.data")
  write(13,*) "A tune:", lat%a%tune/twopi
  write(13,*) "B tune:", lat%b%tune/twopi
  close(13)
  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 *, '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(3) = beam%bunch(1)%particle(:)%vec(3) + leader_y_off
  beam%bunch(1)%particle(:)%vec(5) = beam%bunch(1)%particle(:)%vec(5) + leader_z_off

  write(12,*) "Elements:"

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

  do i=1,size(beam%bunch)
      beam%bunch(i)%particle(:)%vec(3) = 0.001
  end do

  write(11,*) 'Starting element-by-element tracking'
  do i=0,lat%n_ele_track-1
    sel_ele => lat%ele(i+1)
    write(11,*) '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