subroutine sync_beta_path(lat, d_amp_x, d_amp_y)
! Turn on RF. Compute closed orbit with i_dim=6. Then with i_dim=4.
! Compute invariant amplitude of 6D and 4D orbits at IP and
! then fractional difference at IP. Divide by total RF voltage.

  use bmad

  implicit none

  type (lat_struct) lat
  type (lat_struct), save :: lat_save
  type (coord_struct) dorb
  type (coord_struct), allocatable, save :: co6_lat(:), co4_lat(:)
  type (ele_struct) ele
  real(rp) total_volts, d_amp_x, d_amp_y, a_0, a_1
  integer i

  call reallocate_coord( co4_lat, lat%n_ele_track )
  call reallocate_coord( co6_lat, lat%n_ele_track )

 total_volts = 0.

  lat_save = lat

  call set_on_off (rfcavity$, lat, on$)

  call  set_z_tune(lat%branch(0), -0.089 * twopi)
  do i =1, lat%n_ele_track
    if(lat%ele(i)%key == rfcavity$)then
      total_volts = total_volts + lat%ele(i)%value(voltage$)/1.e6
    endif
  end do
  
  co4_lat(0)%vec(:) = 0.
  co6_lat(0)%vec(:) = 0.
  call closed_orbit_calc(lat, co6_lat, 6)
  call closed_orbit_calc(lat, co4_lat, 4)

  ele = lat%ele(0)
  a_0 = ele%a%beta * co4_lat(0)%vec(2)**2 + ele%a%alpha * &
       co4_lat(0)%vec(1) * co4_lat(0)%vec(2) + &
       (ele%a%alpha **2 +1)/ele%a%beta * co4_lat(0)%vec(1)**2
  a_1 = ele%a%beta * co6_lat(0)%vec(2)**2 + ele%a%alpha * &
       co6_lat(0)%vec(1) * co6_lat(0)%vec(2) + &
       (ele%a%alpha **2 +1)/ele%a%beta * co6_lat(0)%vec(1)**2
  
  if (total_volts == 0) then
    print *,' SYNCHRO_BETA: total_volts=0, cannot compute d_amp '
    stop
  endif

  d_amp_x = 0.       
  d_amp_x = (sqrt(a_1)-sqrt(a_0))/(total_volts)
       
  a_0 = ele%b%beta * co4_lat(0)%vec(2)**2 + ele%b%alpha * &
       co4_lat(0)%vec(1)*co4_lat(0)%vec(2) + &
       (ele%b%alpha **2 +1)/ele%b%beta * co4_lat(0)%vec(1)**2
  a_1 = ele%b%beta * co6_lat(0)%vec(2)**2 + ele%b%alpha * &
       co6_lat(0)%vec(1)*co6_lat(0)%vec(2) + &
       (ele%b%alpha **2 +1)/ele%b%beta * co6_lat(0)%vec(1)**2

      
  d_amp_y = 0.
  d_amp_y = (sqrt(a_1)-sqrt(a_0))/(total_volts)

  lat = lat_save

  return
end subroutine sync_beta_path

subroutine sync_beta_volt(lat, d_amp_x, d_amp_y)
! Compute closed orbit with i_dim=6 vs rf voltage
! Compute invariant amplitude of 6D orbits at IP and
! then fractional difference at IP. Divide by change in RF voltage.

  use bmad

  implicit none

  type (lat_struct) lat
  type (lat_struct), save :: lat_save
  type (coord_struct) dorb
  type (coord_struct), allocatable, save :: co_high(:), co_low(:)
  type (ele_struct) ele
  real(rp) total_volts_high, total_volts_low, d_amp_x, d_amp_y, a_0, a_1
  integer i

  total_volts_high = 0.

  call reallocate_coord( co_high, lat%n_ele_max )
  call reallocate_coord( co_low, lat%n_ele_max )

  lat_save = lat

  call set_on_off (rfcavity$, lat, on$)
  call  set_z_tune(lat%branch(0), -0.089 * twopi)
  do i =1, lat%n_ele_track
    if(lat%ele(i)%key == rfcavity$)then
      total_volts_high = total_volts_high + lat%ele(i)%value(voltage$)/1.e6
    endif
  end do
  
  co_high(0)%vec(:) = 0.
  call closed_orbit_calc(lat, co_high,6)
  
  total_volts_low = 0.
  call  set_z_tune(lat%branch(0), -0.001 * twopi)
  do i =1, lat%n_ele_track
    if(lat%ele(i)%key == rfcavity$)then
      total_volts_low = total_volts_low + lat%ele(i)%value(voltage$)/1.e6
    endif
  end do
  
  co_low(0)%vec(:) = 0.
  call closed_orbit_calc(lat, co_low,6)
  
  ele = lat%ele(0)
  a_0 = ele%a%beta * co_low(0)%vec(2)**2 + ele%a%alpha * &
       co_low(0)%vec(1)*co_low(0)%vec(2) + &
       (ele%a%alpha **2 +1)/ele%a%beta * co_low(0)%vec(1)**2
  a_1 = ele%a%beta * co_high(0)%vec(2)**2 + ele%a%alpha * &
       co_high(0)%vec(1)*co_high(0)%vec(2) + &
       (ele%a%alpha **2 +1)/ele%a%beta * co_high(0)%vec(1)**2
  
  if(total_volts_high == 0) then
    print *,' SYNCHRO_BETA: total_volts_high=0, cannot compute d_amp '
    stop
  endif
  
  d_amp_x = 0.       
  d_amp_x = (sqrt(a_1)-sqrt(a_0))/(total_volts_high-total_volts_low)
  
  a_0 = ele%b%beta * co_low(0)%vec(2)**2 + ele%b%alpha * &
       co_low(0)%vec(1) * co_low(0)%vec(2) + &
       (ele%b%alpha **2 +1)/ele%b%beta * co_low(0)%vec(1)**2
  a_1 = ele%b%beta * co_high(0)%vec(2)**2 + ele%b%alpha * &
       co_high(0)%vec(1)*co_high(0)%vec(2) + &
       (ele%b%alpha **2 +1)/ele%b%beta * co_high(0)%vec(1)**2
  
  
  d_amp_y = 0.
  d_amp_y = (sqrt(a_1)-sqrt(a_0))/(total_volts_high-total_volts_low)
  
  lat = lat_save
  
  return
end subroutine sync_beta_volt