!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Subroutine sim_orm_data(ring, bpm, beam, orm_data)
!
! Simulate a full set of Orbit Response Matrix (ORM analysis) measurements.
! Lattice and BPM errors are taken into account.
!
! Input:
!   ring    -- lat_struct
!   bpm     -- det_struct, holds information about a single BPM
!   beam    -- meas_beam_struct; stores beam current, lifetime
!   
!
! Output:
!   orm_data(:) -- orm_data_struct(:).
!--------------------------------------------------------------------------
  subroutine sim_orm_data(ring, bpm, beam, orm_data)

    use bmad
    use sim_decay_mod
    use sim_bpm_mod

    implicit none


    type orm_data_struct
       type(det_struct) :: bpm(1:99)
       real(rp) :: mag_hkick = 0, mag_vkick = 0
       integer :: hkick_ix = 981, vkick_ix = 984
       character(20) :: ele_name = ''
    end type orm_data_struct

    type(lat_struct) :: ring
    type(det_struct) :: bpm(:)
    type(meas_beam_struct) :: beam
    type(coord_struct), allocatable :: orb(:)
    type(orm_data_struct) :: orm_data(:)
    real(rp) :: harvest
    type (all_pointer_struct), allocatable, save :: hkick_pointer(:), vkick_pointer(:)
    integer i, j, n_bpms, n_steering_tot
    logical err
    character(40) lat_type

    lat_type = ring%use_name

    if (.not. allocated(orb)) allocate(orb(0:ring%n_ele_track))
  
    n_steering_tot = ubound(orm_data,1)/2
    n_bpms = ubound(bpm,1)

    ! Reference orbits - generated for both output types
    
    do i=1, n_steering_tot
       
       call closed_orbit_calc(ring, orb, 4)
       
       bpm(:)%vec(1) = orb(bpm(:)%ix_lat)%vec(1) + bpm(:)%x_offset
       bpm(:)%vec(3) = orb(bpm(:)%ix_lat)%vec(3) + bpm(:)%y_offset
       
       do j=1, n_bpms
          call apply_bpm_errors(bpm(j),beam%current, lat_type)
       enddo
       
       orm_data(2*i-1)%bpm(1:n_bpms)%vec(1) = bpm(:)%vec(1)
       orm_data(2*i-1)%bpm(1:n_bpms)%vec(3) = bpm(:)%vec(3)
       
    enddo  ! closes reference orbit loop
  
    do i=1, n_steering_tot

     call pointers_to_attribute(ring, ring%ele(orm_data(2*i)%hkick_ix)%name, "HKICK", .true., hkick_pointer, err)
     hkick_pointer(1)%r = hkick_pointer(1)%r + orm_data(2*i)%mag_hkick ! mag_{h,v}kick is user-defined in the input file
     call pointers_to_attribute(ring, ring%ele(orm_data(2*i)%vkick_ix)%name, "VKICK", .true., vkick_pointer, err)
     vkick_pointer(1)%r = vkick_pointer(1)%r + orm_data(2*i)%mag_vkick ! mag_{h,v}kick is user-defined in the input file

     call closed_orbit_calc(ring, orb, 4)
       
     bpm(:)%vec(1) = orb(bpm(:)%ix_lat)%vec(1) + bpm(:)%x_offset
     bpm(:)%vec(3) = orb(bpm(:)%ix_lat)%vec(3) + bpm(:)%y_offset

     do j=1, n_bpms
        call rotate_meas(bpm(j)%vec(1), bpm(j)%vec(3), bpm(j)%tilt)
        call apply_bpm_errors(bpm(j), beam%current, lat_type)
     enddo

     orm_data(2*i)%bpm(1:n_bpms)%vec(1) = bpm(:)%vec(1)
     orm_data(2*i)%bpm(1:n_bpms)%vec(3) = bpm(:)%vec(3)

     hkick_pointer(1)%r = hkick_pointer(1)%r - orm_data(2*i)%mag_hkick ! reset to zero change
     vkick_pointer(1)%r = hkick_pointer(1)%r - orm_data(2*i)%mag_vkick

  enddo ! closes kick loop
  
end subroutine sim_orm_data