program multibunch
  use bmad
  use beam_mod
  use mode3_mod
  use ecloudmod
  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, lun2, lun3, lun4, lun5
  integer nturns, n_bunch, n_particle
  real(rp) feedback_val/0/
  integer ios

  real(rp) bunch_spacing, bunch_charge

  real(rp)x_off, y_off, z_off

  real(rp)slope, cslope, slopeslope, cslopeslope

  character*19 file_name
  character*3 num

  character*140 lat_file
  character*120 line, last_line
  logical error/.false./
  character*5 active_string
  logical cloud_active/.true./
  logical cloud_motion/.false./
  logical cloud_growth/.false./
  logical cloud_pinch/.false./
  real(rp) cloud_start/0/ 
  real(rp) cloud_start_sxx/1e-4/
  real(rp) cloud_start_syy/1e-3/
  real(rp) cloud_pf/1/

  real(rp) total_charge/0/

  logical all_debug/.false./
  logical last_debug/.true./

  namelist /beam_def/nturns, n_bunch, bunch_spacing, n_particle, bunch_charge, feedback_val
  namelist /beam_offset/x_off, y_off, z_off
  namelist /slope_def/slope, cslope, slopeslope, cslopeslope
  namelist /cloud_settings/ cloud_active, cloud_motion, cloud_growth, cloud_pinch, cloud_start,cloud_start_sxx,cloud_start_syy,cloud_pf

  print *, "Multibunch v2.0"
  
  nargs = cesr_iargc()
  if(nargs == 1 .or. nargs==2)then
    call cesr_getarg(1,lat_file)
    
    call cesr_getarg(2,active_string)
    if (active_string == "false") then
      cloud_active = .false.
    endif
    
    print *, 'Using ', trim(lat_file)
    print *, 'Cloud: ', cloud_active
  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)

  lun2 = lunget()
  open(unit=lun2, file='beam_offsets.in', STATUS ='old')
  read(lun2, nml=beam_offset, IOSTAT=ios)
  close(unit=lun2)
  write(6,nml=beam_offset)

  lun3 = lunget()
  open(unit=lun3, file='snapshot/slope_def.in',STATUS='old')
  read(lun3,nml=slope_def,IOSTAT=ios)
  close(unit=lun3)
  write(6,nml=slope_def)

  lun4 = lunget()
  open(unit=lun4, file='snapshot/cloud_settings.in',STATUS='old')
  read(lun4,nml=cloud_settings,IOSTAT=ios)
  close(unit=lun4)
  write(6,nml=cloud_settings)

  ! === SET SLOPES ===
  call setdebug(all_debug)
  call setSlopes(slope,cslope,slopeslope,cslopeslope)
  call setstate(.false.)
  call setmotion(cloud_motion)
  call setgrowth(cloud_growth)
  call setpinch(cloud_pinch)
  call setstartcharge(cloud_start)
  call setstartsigx(cloud_start_sxx)
  call setstartsigy(cloud_start_syy)
  call setpinchfactor(cloud_pf)
  ! === ========== ===

  print *, "Parsing Lattice:"
  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)
  call track_all(lat,orbit)
  call radiation_integrals (lat, orbit, mode, ix_cache, 0, rad_int)
  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%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

  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

  print *, "Tune: ", lat%a%tune, lat%b%tune, lat%z%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)

  ! For pinch effect calculations
  call setbunchlength(beam_init%sig_z)

  ! Takes ~20,000 turns to have an effect
  bmad_com%radiation_damping_on = .true.
  bmad_com%radiation_fluctuations_on = .true.

  do i=1, size(beam%bunch)
    beam%bunch(i)%particle(:)%vec(1) = x_off + beam%bunch(i)%particle(:)%vec(1)
    beam%bunch(i)%particle(:)%vec(3) = y_off + beam%bunch(i)%particle(:)%vec(3)
    beam%bunch(i)%particle(:)%vec(5) = z_off + beam%bunch(i)%particle(:)%vec(5)
  end do

  ! === Turn to whatever value is needed ===
  call setstate(cloud_active)
  ! === ================================ ===

  call write_turn_data(beam_init,beam,0)
  do i=1,nturns
    if (last_debug .and. i==nturns-9) then
      call setdebug(.true.)
    end if
    call track_beam(lat,beam,ele1=lat%ele(0),ele2=lat%ele(lat%n_ele_track), err=error)  
    call xyfeedback(beam, feedback_val)
    call write_turn_data(beam_init,beam,i)

    total_charge = 0
    do j=1, size(beam%bunch)
      total_charge = total_charge + beam%bunch(j)%charge_live
    end do
    if (total_charge .eq. 0) then
      exit
    end if

  end do

end program