program singleturn
  use bmad
  use beam_mod
  use mode3_mod
  use ecloudmod
  implicit none

  type (lat_struct) lat
  type (ele_struct) ele 
  type (ele_struct), pointer :: sel_ele
  type (track_struct) tracker
  integer bunch_num, part_num, k
  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
  integer nturns, n_bunch, n_particle
  integer ios
  integer multiuse_counter

  real(rp) bunch_spacing, bunch_charge
  real(rp)x_off, y_off, z_off
  real(rp) feedback_val/0/

  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./
  logical cloud_centered/.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/.true./
  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, cloud_centered

  print *, "Singleturn"

!  print *, "[!] Adjusting default error settings for runge_kutta"
!  bmad_com%rel_tol_adaptive_tracking = 1
!  bmad_com%abs_tol_adaptive_tracking = 1
  
  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)

  do i=0,lat%n_ele_track
    write(13,'(i,6es12.4)') i, orbit(i)%vec
  end do


  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))', 'horizontal emittance = ',mode%a%emittance, 'vertical emittance = ', mode%b%emittance 
  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)

  ! 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)
    do j=1,6
      beam%bunch(i)%particle(:)%vec(j) = orbit(0)%vec(j) + beam%bunch(i)%particle(:)%vec(j)
    end do
    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)
  ! === ================================ ===

  !=====================!
  ! Write out log files !
  !=====================!
  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

  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)
    ! Single particle tracking
    do bunch_num = 1, size(beam%bunch)
      write(11,*) "Bunch",bunch_num, "Charge:", beam%bunch(bunch_num)%charge_live
      do part_num = 1, size(beam%bunch(bunch_num)%particle)
        call track1(beam%bunch(bunch_num)%particle(part_num), sel_ele, lat%param, &
                    beam%bunch(bunch_num)%particle(part_num), tracker)
      end do
    end do
    if (cloud_centered) then
      call adjustcloudoffsets(beam)
    end if

    call write_turn_data(beam_init,beam,i+1)
  end do
end program