program space_charge_mc
  use bmad
  use bmadz_interface
  use bsim_cesr_interface
  use scan_parameters
  use sc_track_mod
  use sc_track_mod_mpi
  use ele_noise_mod

  !use mode3_mod
 
  type (lat_struct) ring
  type (ele_struct), pointer :: ele
  type (coord_struct), allocatable :: orb(:)

  type(ele_noise_struct) :: ele_noise(n_ele_noise_max) ! may need to change the bounds on this
  
  integer, parameter :: max_obs_points = 10
  integer i, ix, jx
  integer ios
  integer n_turn
  integer :: seed = -1, seed_used
  integer mpierr

  real(rp) current
  integer :: n_part_track = 1000, n_turns_recalc = 500
  real(rp) :: target_tunes(3) = 0.
  real(rp) :: coupling = 0.005 ! needed for vertical in ideal lattice

  character(200) :: lat_file, file_name

  character(200) :: obs_point(max_obs_points) = ""
  integer :: obs_point_idx(max_obs_points) = -1


  logical ok
  logical :: space_charge_on = .false.
  logical :: auto_bookkeeper = .false.
  logical :: radiation = .true.

  logical spokesnode

  integer readstatus
  integer :: arg_num,iargc

  namelist / parameters / target_tunes, n_part_track, lat_file, n_turn, radiation, &
       current, space_charge_on, n_turns_recalc, coupling, &
       auto_bookkeeper, ele_noise, seed, obs_point, obs_point_idx
  
  !------------------------------Initialize MPI
  CALL mympi_initialize(spokesnode)
  
  ! USE COMMAND LINE ARG TO INPUT FILENAME
  arg_num=iargc()
  if(arg_num==0) then
     file_name='mc_space_charge.in'
  else
     call getarg(1, file_name)
  end if
  call string_trim (file_name, file_name, ix)
  open(unit= 1, file=file_name, action='read', status = 'old', iostat = ios)
  
  if(spokesnode) then
    if(ios.ne.0) then
       print *
       print *, 'ERROR: CANNOT OPEN FILE: ', trim(file_name)
    endif
  endif
    
  rewind(1)
  read(1, nml = parameters,iostat=readstatus)
  if(readstatus > 0) then
     if(spokesnode) then
        print *,"CAN NOT READ IN FROM ",file_name
     endif
     stop
  endif
  
  bmad_com%auto_bookkeeper = auto_bookkeeper
  bmad_com%space_charge_on = space_charge_on
  
  !#############################
  ! Initialize random-number generator from time or with specific value.
  ! if no user-provided seed, only determine random seed from spokesnode
  ! and propagate to other nodes.
  if (seed < 1) then
     seed = 0 ! seed randomizer with CPU time
     if (spokesnode)  then
        call ran_seed_put(seed)
        call ran_seed_get(seed_used)
     endif
  else ! user provided an input seed; use this on all nodes
     seed_used = seed
  endif

  call MPI_Bcast(seed_used, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr)
  call ran_seed_put(seed_used)

  if (spokesnode) write(*,*) "Random seed used: ", seed_used
  !#############################


  !----------------------------Parse lattice
  CALL mympi_bmad_parser(lat_file,ring)

  !=====================================================
  ! Observation point for emittance / beamsize measurements:
  if (.not. (len(obs_point(1)) .eq. 0) .and. .not. (all(obs_point_idx .eq. -1))) then
     write(*,*) 'Both obs_point (string) and obs_point_idx (integer) provided. Choose one! Stopping here...'
     stop
  endif

  i = 1 ! index in obs_point_idx
  ! find ele indices corresponding to requested names. This will also
  ! ensure obs_point_idx(:) is in numerical order, to make life easier.
  ele_loop: do jx = 0, ring%n_ele_max
     ele => ring%ele(jx)
     obs_loop: do ix = 1, size(obs_point) 
        if (len(trim(obs_point(ix))) .eq. 0) cycle obs_loop
        if ( match_reg(upcase(ele%name), upcase(obs_point(ix))) ) then
           obs_point_idx(i) = jx
           i = i+1
           cycle ele_loop
        endif
     enddo obs_loop
  enddo ele_loop

  ! if user didn't specify anything, just use first element in lattice:
  if (all(obs_point_idx .lt. 0)) obs_point_idx(1) = 0
  !=====================================================

  ! set apertures and check wiggler tracking method
  do i = 1, ring.n_ele_max
     if(ring.ele(i).value(x1_limit$) == 0.)ring.ele(i).value(x1_limit$) = 0.05
     if(ring.ele(i).value(x2_limit$) == 0.)ring.ele(i).value(x2_limit$) = 0.05
     if(ring.ele(i).value(y1_limit$) == 0.)ring.ele(i).value(y1_limit$) = 0.05
     if(ring.ele(i).value(y2_limit$) == 0.)ring.ele(i).value(y2_limit$) = 0.05
  enddo
  
  bmad_com%aperture_limit_on = .true.
  
  call set_on_off (rfcavity$, ring, off$)
  call closed_orbit_calc(ring, orb, 4)
  call lat_make_mat6(ring, -1, orb)
  call twiss_at_start(ring)
  call twiss_propagate_all(ring)
  call set_on_off (rfcavity$, ring, on$)

  ! damping / excitation:
  if (radiation .eqv. .true.) then
     bmad_com%radiation_damping_on      = .true.
     bmad_com%radiation_fluctuations_on = .true.
  else
     bmad_com%radiation_damping_on      = .false.
     bmad_com%radiation_fluctuations_on = .false.
  endif

  if (any(target_tunes .gt. 1.e-12)) then
     ok = set_tune3(ring, target_tunes, .true.)
     call closed_orbit_calc(ring, orb, 4)
     call lat_make_mat6(ring, -1, orb)
     call twiss_at_start(ring)
     call twiss_propagate_all(ring)
  endif

  ! prep ele_noise_struct:
  call init_ele_noise(ring, ele_noise)

  ! track the particles, etc.
  call sc_track(ring, n_part_track, n_turn, n_turns_recalc, space_charge_on, current, coupling, ele_noise, obs_point_idx, file_name)

  call mympi_shutdown()
    
end program space_charge_mc