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) bunch_a_emit/0/
  real(rp) bunch_b_emit/0/
  real(rp) bunch_sig_z/0/
  real(rp) bunch_sig_e/0/

  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/.false./
  logical cloud_motion/.false./
  logical cloud_growth/.false./
  logical cloud_pinch/.false./
  logical cloud_centered/.false./
  ! cloud_start is the total charge per unit length of the central gaussian
  ! this can be converted to electron density of the central gaussian with:
  ! - cloud_start / (cloud_start_sxx * cloud_start_syy * e_charge)
  real(rp) cloud_start/0/ 
  real(rp) cloud_pf/0/
  real(rp) cloud_start_sxx/1e-4/
  real(rp) cloud_start_syy/1e-3/

  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, bunch_a_emit, bunch_b_emit, bunch_sig_z, bunch_sig_e
  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, cloud_centered

  print *, "Multibunch"
  
  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)
  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)

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

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

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

  ! Turn on cloud if needed 
  call setstate(.false.)


  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

  if (bunch_a_emit .eq. 0) then
    beam_init%a_emit = mode%a%emittance !beam_init%a_emit = 2.e-9
  else
    beam_init%a_emit = bunch_a_emit
  endif
  if (bunch_b_emit .eq. 0) then
    beam_init%b_emit = mode%b%emittance !beam_init%b_emit = 10.e-12
  else
    beam_init%b_emit = bunch_b_emit
  endif
  if (bunch_sig_z .eq. 0) then
    beam_init%sig_z = mode%sig_z        !beam_init%sig_z = 0.01
  else
    beam_init%sig_z = bunch_sig_z
  endif
  if (bunch_sig_e .eq. 0) then
    beam_init%sig_e = mode%sigE_E       !beam_init%sig_e = 0.0008
  else
    beam_init%sig_E = bunch_sig_e
  endif


  !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/twopi, lat%b%tune/twopi, lat%z%tune/twopi
  
  print '(2(a26,es12.4))', 'horizontal emittance = ',beam_init%a_emit, 'vertical emittance = ', beam_init%b_emit
  print '(2(a26,es12.4))', 'bunch length = ',beam_init%sig_z,'sig E/E = ', beam_init%sig_e
  print '(a26,es12.4)','bunch charge =', beam_init%bunch_charge
  print '(1(a26,i10))','number of bunches =', beam_init%n_bunch
  print '(1(a26,i10))', 'particles/bunch =', beam_init%n_particle
  print '(1(a26,i10))', 'number of turns =', nturns

  print *, 'Cloud Settings:'
  print '(a26,es12.4)', 'Linear Charge Density =', cloud_start
  print '(a26,2es12.4)', 'Gaussian shape =', cloud_start_sxx, cloud_start_syy
  print '(a26,es12.4)', 'Electron Density at beam =', - cloud_start / (cloud_start_sxx * cloud_start_syy * e_charge)

  write(6,nml=beam_def)
  write(6,nml=beam_offset)
  ! write(6,nml=slope_def)
  write(6,nml=cloud_settings)

  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 setdebug(all_debug)
  !call setSlopes(slope,cslope,slopeslope,cslopeslope)
  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)
  !=== ================================ ===

  if (.not. cloud_centered) then
    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
  else
    do i=1,nturns
      if (last_debug .and. i==nturns-9) then
        call setdebug(.true.)
      end if
      do j=0,lat%n_ele_track-1
        call track_beam(lat,beam,ele1=lat%ele(j),ele2=lat%ele(j+1), err=error)  
        call adjustcloudoffsets(beam)
      end do
      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 if

end program