! Subroutine dx_dvcros5_CALC (LAT, VC_DIFF)
!
! Subroutine to calculate the dependence of orbit at detectors 0E and 0W
! on vcrosing 5
!
!
! Input:
!     LAT      -- lat_struct: Lat
!
! Output:
!-    vc_diff(0:lat%n_ele_track)  -- Coord_struct
!                        Difference between closed orbits with vertical
!                          separators +- dk antisymmetrically

subroutine dx_dvcros5_calc ( lat, vc_diff)

  use bmadz_mod
  use bmadz_interface, dummy => dx_dvcros5_calc
  use bmad_interface

  implicit none

  type (lat_struct)  lat
  type (coord_struct), allocatable ::  vc_diff(:)
  type (coord_struct), allocatable, save ::  coplus(:)
  type (coord_struct), allocatable, save ::  cominus(:)

  real(rp) :: kick_west, kick_east, dk, vec0(6) = 0
  integer :: i, ix_west = 0, ix_east = 0, n_track 

  n_track = lat%n_ele_track
  call reallocate_coord(  coplus, lat)
  call reallocate_coord(  cominus, lat)


!
! find vertical separators

  i=0
  do while(ix_east == 0 .or. ix_west == 0)
   if(lat%ele(i)%name == 'V_SEP_48W')ix_west = i
   if(lat%ele(i)%name == 'V_SEP_48E')ix_east = i
   i=i+1
   if( i > n_track) then
    print *,' DX_DVCROS5_CALC: cannot find vertical separators'
    stop
   endif
  enddo


  kick_west = lat%ele(ix_west)%value(vkick$)
  kick_east = lat%ele(ix_east)%value(vkick$)
  dk = 0.1

  lat%ele(ix_west)%value(vkick$) = kick_west*(1+dk)
    call make_mat6(lat%ele(ix_west), lat%param)
  lat%ele(ix_west)%value(vkick$) = kick_east*(1-dk)
    call make_mat6(lat%ele(ix_east), lat%param)
  call closed_orbit_calc( lat, coplus, 4, +1 )

  lat%ele(ix_west)%value(vkick$) = kick_west*(1-dk)
    call make_mat6(lat%ele(ix_west), lat%param)
  lat%ele(ix_west)%value(vkick$) = kick_east*(1+dk)
    call make_mat6(lat%ele(ix_east), lat%param)

  call init_coord (cominus(n_track), vec0, lat%ele(n_track), downstream_end$, antiparticle(lat%param%particle), -1)
  call closed_orbit_calc( lat, cominus, 4, -1 )

  lat%ele(ix_west)%value(vkick$) = kick_west
  lat%ele(ix_east)%value(vkick$) = kick_east
  call make_mat6(lat%ele(ix_west), lat%param)
  call make_mat6(lat%ele(ix_east), lat%param)
 

  do i = 1,n_track
   vc_diff(i)%vec(:) = (coplus(i)%vec(:) - cominus(i)%vec(:))/(2*dk)
  end do

  return
end