!+
! subroutine LRBBI_SETUP (lat_in, lat_out, particle, i_train, j_car, n_trains_Tot, n_cars, current, rec_taylor, bbiwt, custom_pat)
!
! Subroutine to add lrbbi elements to lat at relevant crossing points 
!
! Input:
!   LAT_IN -- lat_struct: Lat containing the lattice to be modified
!   particle  -- Integer : (+-1) for positrons or electrons
!   i_train   -- Integer : Train of test particle
!   j_car     -- Integer : Car of test particle
!   n_trains_tot -- Integer : Total number of trains
!   n_cars    -- Integer : Total number of cars/train
!   current   -- real : bunch current in ma
!   rec_taylor -- Logical : passed to stash_split_lat. save taylor as digested file
!   bbiwt(:,:)  -- Real, optional : Weights for the current in each 
!                          strong bunch. Range 0 to 1. 
!                  The dimensions should be defined by the  calling 
!                         routine to be (n_trains_Tot,n_cars)
!   custom_pat -- logical,optional : use custom bunch pattern ?
!
! Output:
!    LAT_OUT -- lat_struct: Lat_out includes Beam-beam elements 
!                         at each of crossing points for bunch i_train, j_car
!                         

subroutine lrbbi_setup (lat_in, lat_out, particle, i_train, j_car, n_trains_tot, n_cars,  current, rec_taylor, bbiwt, custom_pat)

use bmad
use bmadz_mod
use cesr_mod
use bunchcross_mod
use constraints_mod

implicit none

interface
 subroutine MARK_LRBBI_ONLY(master_lat, master_lat_oppos, lat, crossings)
  use bmad_struct
  use bmad_interface                    
  implicit none
  type (lat_struct) :: lat
  type (lat_struct) :: master_lat, master_lat_oppos
  real(rp), dimension(:,:) :: crossings
 end subroutine
end interface

type (lat_struct) lat_in, lat_out, lat_1
type (lat_struct), save:: master_lat, master_lat_oppos
type (coord_struct), allocatable, save :: co(:)
type (ele_struct) ele

integer particle, i_train, j_car, n_trains_tot, n_cars
integer i_bunch
integer ierr
integer crossings_tot, bunch_tot, total
integer, dimension(:), allocatable :: ix_lrbbi, ix_master
integer k
integer n
integer i, j
integer ib,it

! Allow adjustment of train and bunch spacing in call to cesr_crossings
! Default spacing
integer :: train_spacing(1:10)= (/140, 140, 147,0,0,0,0,0,0,0/)
integer :: n_car_spacing(1:10)=(/7,0,0,0,0,0,0,0,0,0/)
! Example of 6, 8 ns spacing
!  integer :: n_car_spacing(1:10) /3, 4, 8*0/
! Example of 2, 12 ns spacing
!  integer :: n_car_spacing(1:10) /1, 6, 8*0/

real(rp), dimension(:,:), allocatable :: crossings_1
real(rp), dimension(:), allocatable :: cross_positions
! ADDED TO FIX cesr_crossings
integer, dimension(:), allocatable :: ptrain
integer, dimension(:), allocatable :: pcar
 
real(rp) current

logical rec_taylor

real(rp), dimension(:,:), optional :: bbiwt
logical, optional :: custom_pat

type (pc_struct)  pc
type (constraint_struct)  con

call reallocate_coord( co, lat_in%n_ele_max )

! init

lat_in%param%particle = particle
call twiss_at_start(lat_in)
co(0)%vec = 0.
call closed_orbit_calc(lat_in, co, 4)
call track_all (lat_in, co)
call lat_make_mat6(lat_in,-1,co)
call twiss_at_start(lat_in)
call twiss_propagate_all (lat_in)

! set up for lrbbi
if (present(custom_pat)) then
   if (custom_pat) then
     con%BunchPattern='bunchcross.in'
     call bunchcross (lat_in, con, pc)
     n_trains_tot = con%n_trains
     n_cars = con%n_cars
   endif
endif

if ( i_train*j_car * n_cars * n_trains_tot * current == 0) then
  print *,' LRBBI_SETUP: bunch specifications incomplete '
  stop
endif

i_bunch = (i_train - 1) * n_cars + j_car 

lat_in%param%n_part   = current*0.001 *(lat_in%param%total_length/c_light)/e_charge
master_lat = lat_in
master_lat_oppos = lat_in
lat_1 = lat_in
master_lat_oppos%param%n_part = lat_in%param%n_part
master_lat_oppos%param%particle = -particle

bunch_tot = n_trains_tot * n_cars
crossings_tot = bunch_tot * 2
total = crossings_tot

allocate(cross_positions(1:crossings_tot), stat=ierr)
if(ierr .ne. 0) then
	print*, "CROSS_POSTIONS: ALLOCATION REQUEST DENIED."
	call err_exit
endif

! ADDED TO FIX cesr_crossings
allocate(ptrain(1:crossings_tot), stat=ierr)
if(ierr .ne. 0) then
   print*, "PTRAIN: ALLOCATION REQUEST DENIED."
   call err_exit
endif

allocate(pcar(1:crossings_tot), stat=ierr)
if(ierr .ne. 0) then
   print*, "PCAR: ALLOCATION REQUEST DENIED."
   call err_exit
endif

allocate(crossings_1(1:total, 1:5), stat=ierr)
if(ierr .ne. 0) then
	print*, "CROSSINGS_1: ALLOCATION REQUEST DENIED."
	call err_exit
endif

allocate(ix_lrbbi(1:crossings_tot), stat=ierr)
allocate(ix_master(1:crossings_tot), stat=ierr)


if (present(custom_pat)) then
   ! use custom bunch pattern
   if (custom_pat) then
     do j=1, pc%total_pc
        cross_positions(j)=pc%cross(j)%ele%s / lat_in%param%total_length
     enddo
   else
     call cesr_crossings(i_train, j_car, particle,n_trains_tot,n_cars,cross_positions, ptrain, pcar, &
                          n_car_spacing, train_spacing )
   endif       
else
   call cesr_crossings(i_train, j_car, particle,n_trains_tot,n_cars,cross_positions, ptrain, pcar, &
                          n_car_spacing, train_spacing )
endif

do j = 1, crossings_tot
 crossings_1(j,1) = cross_positions(j)
 crossings_1(j,2) = 1
 crossings_1(j,3) = j  
 crossings_1(j,4) = 0  
 crossings_1(j,5) = 0  
enddo

call MARK_LRBBI_only(master_lat, master_lat_oppos, lat_1, crossings_1)
call stash_split_lat (i_train, j_car, particle, n_trains_tot, n_cars,master_lat, rec_taylor)  
call stash_split_lat (i_train, j_car, particle, n_trains_tot, n_cars,master_lat_oppos, rec_taylor)  
call stash_split_lat (i_train, j_car, particle, n_trains_tot, n_cars,lat_1, rec_taylor)

do i=1, crossings_tot
 find: do k=1,total
  if(int(crossings_1(k,2)) == 1 .and. int(crossings_1(k,3)) == i)ix_lrbbi(i) = crossings_1(k,4)
  if(int(crossings_1(k,2)) == 1 .and. int(crossings_1(k,3)) == i)ix_master(i) = crossings_1(k,5)
 enddo find
end do

call make_lrbbi(master_lat_oppos, lat_1, ix_lrbbi, ix_master)
 
lat_out = lat_1

! If the bunch weight input argument is present
! weight the parasitic crossings with the relative bunch current 
! of the strong beam. For this, we use the independent variable
! ele(charge$). Until now, this has been simply +1 or -1. But it
! is the only element-specific parameter used to calculate the
! dependent variable bbi_const, which gives the magnitude of
! the kick.  18 January 2006  JAC
!     bbi_const$ = param%n_part * m_electron * charge$ * r_e /
!                  (2 * pi * E_TOT$ * (sig_x$ + sig_y$)
!

if ( present(bbiwt) .and. .not. logic_option(.false., (custom_pat))) then
  do i = 1, crossings_tot
    j = ix_lrbbi(i)
    ib = int(pcar(i))
    it = int(ptrain(i))
    lat_out%ele(j)%value(charge$) = bbiwt(it,ib) * lat_out%ele(j)%value(charge$)
  end do
endif

!

deallocate(cross_positions, stat=ierr)
deallocate(crossings_1, stat=ierr)
deallocate(ix_lrbbi, stat=ierr)
deallocate(ix_master, stat=ierr)
deallocate (ptrain, stat=ierr)
deallocate (pcar, stat=ierr)

end subroutine lrbbi_setup