! Subroutine dx_dpretz13_CALC (LAT, prz13_DIFF)
!
! Subroutine to calculate the dependence of orbit at detectors 0E and 0W
! on pretzing 13.
!
!
! Input:
!     LAT      -- lat_struct: Lat
!
! Output:
!-    prz1_diff(0:n_ele_lat)  -- Coord_struct
!                        Difference between closed orbits with horizontal
!                          separators at 8 +- dk symmetrically

subroutine dx_dpretz13_calc ( lat, prz13_diff)

  use bmadz_mod
  use bmadz_interface, dummy => dx_dpretz13_calc
  use bmad_interface

  implicit none

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

  real(rp) :: kick_8w,kick_8e, dk, vec0(6) = 0
  integer i, n_track
  integer :: ix_8w =0, ix_8e = 0

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

!
! find horizontal separators

  i=1
   do while(ix_8w == 0 .or. ix_8e == 0)
   if(lat%ele(i)%name == 'H_SEP_08W')ix_8w=i
   if(lat%ele(i)%name == 'H_SEP_08E')ix_8e=i
   i=i+1
   if( i > n_track) then
    print *,' DX_DPRETZ13_CALC: cannot find horizontal separators'
    stop
   endif
  enddo


  kick_8w = lat%ele(ix_8w)%value(hkick$)
  kick_8e = lat%ele(ix_8e)%value(hkick$)

  dk = 0.1

  lat%ele(ix_8w)%value(hkick$) = kick_8w*(1+dk)
  call attribute_bookkeeper(lat%ele(ix_8w), .true.)
  call make_mat6(lat%ele(ix_8w), lat%param)

  lat%ele(ix_8e)%value(hkick$) = kick_8e*(1-dk)
  call attribute_bookkeeper(lat%ele(ix_8e), .true.)
  call make_mat6(lat%ele(ix_8e), lat%param)

  call closed_orbit_calc( lat, coplus, 4, +1 )

  lat%ele(ix_8w)%value(hkick$) = kick_8w*(1-dk)
  call attribute_bookkeeper(lat%ele(ix_8w), .true.)
  call make_mat6(lat%ele(ix_8w), lat%param)

  lat%ele(ix_8e)%value(hkick$) = kick_8e*(1+dk)
  call attribute_bookkeeper(lat%ele(ix_8e), .true.)
  call make_mat6(lat%ele(ix_8e), 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_8w)%value(hkick$) = kick_8w
  call attribute_bookkeeper(lat%ele(ix_8w), .true.)
  call make_mat6(lat%ele(ix_8w), lat%param)

  lat%ele(ix_8e)%value(hkick$) = kick_8e
  call attribute_bookkeeper(lat%ele(ix_8e), .true.)
  call make_mat6(lat%ele(ix_8e), lat%param)

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

end