!+
! subroutine beambeam_separation (lat, delta_ip, i_dim, diff_co)
!
! Subroutine to compute differential displacement of strong and weak beam at IP
!
! Input:
!   LAT -- lat_struct: Lat containing the lattice to be modified
! Output:
!    LAT -- lat_struct: Lat includes Beam-beam element at IP
!    delta_ip -- Coord_struct: Relative displacement of two beams at IP                     
!
!
! Added diff_co as mandatory argument    16 Dec 2005  JAC

subroutine beambeam_separation (lat, delta_ip, i_dim, diff_co)

  use bmad_interface
  use bmadz_mod

  implicit none

  type ( lat_struct ) lat
  type (lat_struct), save :: lat_oppos
  type (coord_struct), allocatable, save ::  co(:), co_oppos(:)
  type (coord_struct) delta_ip
!   type (coord_struct), optional, allocatable :: diff_co(:)
   type (coord_struct), optional :: diff_co(0:)
! Discovered that optional, intent(out) results in shifting contents by one position 
! in array    16 Dec 2005
  type (ele_struct) ele

  character(*), parameter :: r_name='beambeam_separation'

  integer i, i_dim
  logical err

  call reallocate_coord( co, lat%n_ele_max )
  call reallocate_coord( co_oppos, lat%n_ele_max )

  if(lat%ele(0)%a%beta == 0.)call twiss_at_start(lat)

  lat_oppos = lat
  lat_oppos%param%n_part = 0.

   do i=0,lat_oppos%n_ele_max
    co_oppos(i)%vec = 0.
   enddo

!
  co_oppos(:)%species = antiparticle(lat_oppos%param%particle)
  lat_oppos%param%particle = co_oppos(0)%species
  call closed_orbit_calc (lat_oppos, co_oppos, i_dim, 1, err_flag = err)

! update offsets for lat
  do i = 0, lat_oppos%n_ele_track
   ele = lat%ele(i)
   if( ele%key == beambeam$ )then

     lat%ele(i)%value(x_offset$) = co_oppos(i)%vec(1)
     lat%ele(i)%value(y_offset$) = co_oppos(i)%vec(3)


    ! DLR suggested remaking BBI matrices due to offset definition
    ! in BBI elements requilat a call to bookkeeper  16 dec 2005

     call lat_make_mat6(lat, i )

   endif

  end do

!

  co(0)%vec(:) = 0.
  co(0)%species = lat%param%particle
  call closed_orbit_calc (lat, co, i_dim )

  delta_ip%vec(1:i_dim) = co_oppos(0)%vec(1:i_dim) - co(0)%vec(1:i_dim)

  if ( present (diff_co) ) then
   do i = 0, lat_oppos%n_ele_track
    diff_co(i)%vec = co_oppos(i)%vec - co(i)%vec
   enddo
  endif

  call out_io(s_info$,r_name,(/ 'closed orbit strong beam = \4e12.4\ ', 'closed orbit   weak beam = \4e12.4\ ' /) ,&
       r_array = (/ co_oppos(0)%vec(1:4), co(0)%vec(1:4)/) )

end subroutine beambeam_separation