!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Subroutine sim_tbt_data(ring, bpm, n_turns, track6x6, n_TTs, tt_params, n_damping, &
!                                     init_vec, bpm_save, tt_tunes, ele_noise)
!
! Simulate a turn-by-turn measurement. Lattice and BPM errors 
! are taken into account. Note: be sure to propagate Twiss parameters, etc.
! before this routine.
!
! Input:
!   ring  -- lat_struct
!   bpm(:) -- det_struct(:), each ele holds information about a single BPM
!   beam   -- meas_beam_struct; hold info about beam current, lifetime
!   n_turns -- integer, total number of turns to sample in the measurement
!   track6x6 -- logical. If .true., use full 6x6 damping with resonant 
!              excitation. If .false., offset particle and track as 4x4.
!   n_TTs    -- integer, optional. number of tune trackers. if zero, no damping/excitation used
!   tt_params -- type(tt_param_struct), optional. for tune trackers.
!   id        -- integer (max_tt), optional. for tune trackers.
!
!   n_damping -- integer, optional. If track_6x6 = .true., how many turns to damp 
!                before recording orbit data.
!   init_vec -- real(rp)(6), in meters. How far to offset particles before tracking.
!   ele_noise(:) -- ele_noise_struct(:), optional.  Noise to apply to elements when tracking.
!                   Assumes init_ele_noise has already been called with ele_noise and ring as parameters.
!
! Output:
!   bpm_save(:,:)  -- det_struct(n_turns,n_bpms). BPM data from all turns.
!   tt_tunes  -- optional, (n_turns,n_TTs) array containing all tunetracker tunes from all turns collected
!--------------------------------------------------------------------------
  subroutine sim_tbt_data(ring, bpm, n_turns, track6x6, n_TTs, tt_params, id, n_damping, init_vec, bpm_save, tt_tunes, ele_noise)
  
    use bmad
    use sim_bpm_mod
    use tune_tracker_mod
    use tt_tools_mod
    use ele_noise_mod

    implicit none
    type(lat_struct) :: ring
    type(lat_struct), save :: ring_ideal
    type(coord_struct), allocatable :: co(:), orb(:)
    type(det_data_struct), allocatable :: bpm_save(:,:)
    type(det_struct) :: bpm(:)
    type(tt_param_struct), optional :: tt_params(:)
    type(ele_noise_struct), optional :: ele_noise(:)
    type(all_pointer_struct) a_ptr

    integer, intent(in), optional :: id(:)
    integer, optional :: n_TTs

    logical, intent(in) :: track6x6
    integer, intent(in) :: n_turns
    real(rp), intent(in) :: init_vec(6)
    integer, optional :: n_damping
    integer :: n_damping_use = 0, n_TTs_use = 0
    real(rp) :: harvest

    integer :: i_turn = 0, ix, jx, n_bpms, loc
    character(40) plane
    logical :: err, everything_ok
    logical, save :: first_pass = .true.
    real(rp), allocatable, optional :: tt_tunes(:,:)
    real(rp), allocatable :: tt_tunes_temp(:)  
    real(rp), allocatable :: on_off_save(:)

    real(rp) :: kick0(1:max_tt) = 0.
    real(rp) :: current = 1000000
    character(40) lat_type

    lat_type = ring%use_name  

    n_bpms = ubound(bpm,1)
    if (present(n_damping)) n_damping_use = n_damping
    if (present(n_TTs)) then
       n_TTs_use = n_TTs
    else
       n_TTs_use = 0
    endif
    
    if (.not. allocated(bpm_save)) allocate(bpm_save(n_turns,n_bpms))

    if (n_TTs_use .gt. 0) then
       if (present(tt_tunes)) then
          if (.not. allocated(tt_tunes)) allocate(tt_tunes(n_turns,n_TTs_use))
       end if
       do ix = 1, n_TTs_use
          loc = tt_params(ix)%kck_loc
          if (tt_params(ix)%orientation .eq. 'h') then
             plane = 'HKICK'
          elseif (tt_params(ix)%orientation .eq. 'v') then
             plane = 'VKICK'
          elseif (tt_params(ix)%orientation .eq. 'l') then
             plane = 'PHI0'
          endif
          kick0(ix) = value_of_attribute(ring%ele(loc),plane,err)
       enddo
    endif

    call set_on_off(rfcavity$, ring, save_state$, saved_values = on_off_save)
    if (track6x6 .eqv. .true.) then
       call set_on_off(rfcavity$, ring, on$)
       call closed_orbit_calc(ring, co, 6)
    else
       call set_on_off(rfcavity$, ring, off$)
       call closed_orbit_calc(ring, co, 4)
    endif
    
    ! Initialize ele_noise if present.
    if (present(ele_noise)) then
       if (any(ele_noise(:)%initialized) .eqv. .false.) then
          call init_ele_noise(ring, ele_noise)
          ring_ideal = ring  ! Save the ideal no-noise lattice
       endif
    end if

    if (.not. allocated(orb)) allocate(orb(0:ring%n_ele_track))
    
    orb(0)%vec(:) = co(0)%vec(:) + init_vec(:) ! displace particle by some amount, to start oscillations
    
    i_turn = 0 ! reinitialize, in case of sequential calls

    do ix=1,n_damping_use + n_turns ! n_damping = # of turns to wait before measuring

       ! useful progress report output
       if ((ix .eq. 1) .and. (n_damping_use .gt. 0)) write(*,'(a,i0,a)') "Start particle tracking. Damp for ", n_damping_use, " turns"
       if (ix .eq. (n_damping_use + 1)) write(*,'(a,i0,a)') "***Begin saving turns data for ", n_turns, " turns***"
       if ((ix .lt. n_damping_use+1) .and. (mod(ix,5000) .eq. 0) .or. (ix .eq. n_damping_use)) write(*,'(a,i0,a,i0)') "Damping turns: ", ix, "/", n_damping_use
       if ((ix .gt. n_damping_use) .and. (mod((ix-n_damping_use),500) .eq. 0) .or. (ix .eq. n_damping_use + n_turns)) write(*,'(a,i0,a,i0)') "Recording turns: ", ix-n_damping_use, "/", n_turns

       if (present(ele_noise)) then
          if (.not. all(len_trim(ele_noise(:)%attrib_name) .eq. 0)) then
             call apply_ele_noise(ring, ring_ideal, ele_noise, ix)
             if (track6x6) then
                call twiss3_at_start(ring, everything_ok) 
                call twiss3_propagate_all(ring)
             else
                call twiss_at_start(ring)
                call twiss_propagate_all(ring)
             endif
          endif
       end if
       call track_all(ring, orb)

       if ((n_TTs_use .gt. 0) .and. (track6x6 .eqv. .true.)) then ! damping + tune trackers on 
          call track_tunes(ring, orb, tt_params, n_TTs_use, id, tt_tunes_temp, kick0)
       endif
       

       orb(0) = orb(ring%n_ele_track) ! also saves time info, if necessary

       if (ix .gt. n_damping_use) then ! save this data

          i_turn = ix - n_damping_use

          if (n_TTs_use .gt. 0 .and. present(tt_tunes)) then
             tt_tunes(i_turn,:) = tt_tunes_temp(:) 
          endif

          do jx = 1, 6
          bpm(:)%vec(jx) = orb(bpm(:)%ix_lat)%vec(jx)
       enddo

          do jx = 1, n_bpms  
             if (bmad_com%absolute_time_tracking) then
                call apply_bpm_errors(bpm(jx), current, lat_type,.true.)  ! note: also transfers vec(6)-->amp(4) button signals
             else
                call apply_bpm_errors(bpm(jx), current, lat_type)  ! note: also transfers vec(6)-->amp(4) button signals
             endif
  

             bpm_save(i_turn,jx)%vec(:) = bpm(jx)%vec(:)
             bpm_save(i_turn,jx)%amp(:) = bpm(jx)%amp(:)
             bpm_save(i_turn,jx)%ix_db  = bpm(jx)%ix_db
             bpm_save(i_turn,jx)%ix_lat = bpm(jx)%ix_lat
          enddo
       endif

    enddo

    ! Return the ring to its ideal state. (No noise)
    if (present(ele_noise)) ring = ring_ideal

    ! in case tune trackers were used, reset kicks to nominal in the lattice
    orb(0)%vec(:) = init_vec(:)
       do ix = 1, n_TTs_use
          loc = tt_params(ix)%kck_loc
          if (tt_params(ix)%orientation .eq. 'h') then
             plane = 'HKICK'
          elseif (tt_params(ix)%orientation .eq. 'v') then
             plane = 'VKICK'
          elseif (tt_params(ix)%orientation .eq. 'l') then
             plane = 'PHI0'
          endif
          call pointer_to_attribute(ring%ele(loc),plane,.true.,a_ptr,err)
          a_ptr%r = kick0(ix)
          call set_flags_for_changed_attribute(ring%ele(loc), a_ptr%r)
       enddo ! n_TTs_use

    call set_on_off(rfcavity$, ring, restore_state$, saved_values = on_off_save)
    call lattice_bookkeeper(ring)
           
  end subroutine sim_tbt_data