!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Subroutine sim_gain_map_data(this_bpm, beam, nx, ny, dx, dy, gain_map_data)
!
! Simulate a gain map measurement at one BPM. Each measurement is a grid of 
! orbit measurements spanning the cross-section of a BPM. Some BPM errors are 
! taken into account, but not lattice errors.
!
! NOTE: Only use this for CESR-type lattices!
!
! Input:
!   this_bpm -- det_struct, holds information about a single BPM
!   beam     -- meas_beam_struct; holds current, I0, and lifetime information.
!   nx, ny   -- integers, number of horizontal and vertical points in the grid
!   dx, dy   -- real(rp), distance between data points in x and y (meters)
!
! Output:
!   gain_map_data(:) -- gain_map_data_struct(nx*ny). Contains (x,y) orbit and 4 
!                       button amplitudes for each measurement.
!--------------------------------------------------------------------------
subroutine sim_gain_map_data(this_bpm, beam, nx, ny, dx, dy, gain_map_data)
  
  use bmad
  use sim_bpm_mod
  use sim_utils
  use nonlin_bpm_mod
  use sim_decay_mod
 
  
  implicit none



  type gain_map_data_struct
     real(rp) :: vec(6) = 0
     real(rp) :: amp(4) = 0
  end type gain_map_data_struct


  type(det_struct) :: this_bpm
  type(meas_beam_struct) :: beam
  type(gain_map_data_struct), allocatable :: gain_map_data(:)
  integer, intent(in) :: nx, ny
  real(rp), intent(in) :: dx, dy
  real(rp) :: harvest
  real(rp) I0,t ! initial current, total time
  
  integer :: ix, jx, det_ix, j = 0 ! j = orbit index
  real(rp) x0, y0, d(4,0:2,0:2)
  character(40) :: lat_type = 'CESR'
  
  if (.not. allocated(gain_map_data)) allocate(gain_map_data(nx*ny))
  
  
  j = 0 ! need to reinitialize on every call
  x0 = -((nx-1) * dx)/2
  y0 = -((ny-1) * dy)/2
  h_loop: do ix = 0, nx-1 ! loop over horizontal kicks
     v_loop: do jx = 0, ny-1 ! vert kicks
        
        j=j+1 
        
        this_bpm%vec(1) = x0 + ix * dx
        this_bpm%vec(3) = y0 + jx * dy
        
        t = (30.*60/(nx*ny))*j
        call sim_decay(beam, t) ! assume 30mins for measurement.
        
        call nonlin_bpm_set_pointers(this_bpm%ix_db)
        call nonlin_bpm_interpolate((/this_bpm%vec(1),this_bpm%vec(3),beam%current/),d)
        this_bpm%amp = d(:,0,0)
        
        call apply_bpm_errors(this_bpm, beam%current, lat_type)
        
        
        ! now format such that bpm_gain_anal can read it in:
        
        gain_map_data(j)%vec(1) = this_bpm%vec(1)
        gain_map_data(j)%vec(3) = this_bpm%vec(3)
        gain_map_data(j)%amp   = this_bpm%amp
        
        
     enddo v_loop ! loop over y-kicks
     
  enddo h_loop ! loop over x-kicks
  
  
end subroutine sim_gain_map_data