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

subroutine dx_dpretz1_calc ( lat, prz1_diff)

  use bmadz_mod
  use bmadz_interface, dummy => dx_dpretz1_calc
  use bmad_interface

  implicit none

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

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

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


!
! find horizontal separators

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


  kick_8w = lat%ele(ix_8w)%value(hkick$)
  kick_8e = lat%ele(ix_8e)%value(hkick$)
  kick_45w = lat%ele(ix_45w)%value(hkick$)
  kick_45e = lat%ele(ix_45e)%value(hkick$)
  dk = 0.1

  lat%ele(ix_8w)%value(hkick$) = kick_8w*(1+dk)
  lat%ele(ix_8e)%value(hkick$) = kick_8e*(1+dk)
  lat%ele(ix_45w)%value(hkick$) = kick_45w*(1+dk)
  lat%ele(ix_45e)%value(hkick$) = kick_45e*(1+dk)

  call make_mat6(lat%ele(ix_8w), lat%param)
  call make_mat6(lat%ele(ix_8e), lat%param)
  call make_mat6(lat%ele(ix_45w), lat%param)
  call make_mat6(lat%ele(ix_45e), lat%param)

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

  lat%ele(ix_8w)%value(hkick$) = kick_8w*(1-dk)
  lat%ele(ix_8e)%value(hkick$) = kick_8e*(1-dk)
  lat%ele(ix_45w)%value(hkick$) = kick_45w*(1-dk)
  lat%ele(ix_45e)%value(hkick$) = kick_45e*(1-dk)

  call make_mat6(lat%ele(ix_8w), lat%param)
  call make_mat6(lat%ele(ix_8e), lat%param)
  call make_mat6(lat%ele(ix_45w), lat%param)
  call make_mat6(lat%ele(ix_45e), 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
  lat%ele(ix_8e)%value(hkick$) = kick_8e
  lat%ele(ix_45w)%value(hkick$) = kick_45w
  lat%ele(ix_45e)%value(hkick$) = kick_45e

  call make_mat6(lat%ele(ix_8w), lat%param)
  call make_mat6(lat%ele(ix_8e), lat%param)
  call make_mat6(lat%ele(ix_45w), lat%param)
  call make_mat6(lat%ele(ix_45e), lat%param)

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

end