subroutine sext_energy (lat, co, con, res_cos, res_sin, res_amp, max_x)

! Output
! res_amp :  is the fourier amplitude for m=21 of driving term for 2Qx-Qs resonance line.
! max_x :  is the max amplitude of particle tracked on resonance for con%qs_2qx_turns

  use bmad
  use bmadz_interface
  use constraints_mod

  implicit none

  type (lat_struct) lat
  type (lat_struct), save :: lat_1
  type (coord_struct) start, init, co
  type (coord_struct), allocatable, save ::  orbit(:), co_lat(:)
  type (constraint_struct) con
  
  real(rp) res_cos, res_sin, res_amp
  real(rp) qy/0.58/, qz/-0.089/, phi_x, phi_y
  real(rp), allocatable :: dk1(:)
  real(rp) max_x, max_arc
  real(rp) qx

  integer i
  integer int_Q_x, int_Q_y

  logical ok
  logical, save :: first = .true.
  
!  track

  if (any(con%c(1:con%n_constraint)%variable == Qs_2Qx_track$)) then

    lat_1 = lat

    if (first) then
      print '(a31,i5,a6)',' SEXT_ENERGY: tracking set for ', con%qs_2qx_turns,' turns'
      first=.false.
    endif

! set synch tune
    call set_on_off (rfcavity$, lat_1, on$)
    call  set_z_tune(lat_1%branch(0), qz*twopi)

!set betatron tunes
    qx = (1.-qz)/2   !2Qz+Qx resonance
    qx=0.542
    
    int_Q_x = int(lat_1%ele(lat_1%n_ele_track)%a%phi / twopi)
    int_Q_y = int(lat_1%ele(lat_1%n_ele_track)%b%phi / twopi)
    phi_x = (int_Q_x + qx)*twopi
    phi_y = (int_Q_y + qy)*twopi
    
    call reallocate_coord(orbit,lat_1%n_ele_max)
    allocate(dk1(lat_1%n_ele_max))
    call choose_quads(lat_1, dk1)
    forall(i=0:lat%n_ele_max)orbit(i)%vec(1:6)=0.
    call custom_set_tune (phi_x, phi_y, dk1, lat_1, orbit, ok)
!    deallocate(orbit)
    deallocate(dk1)
    if (.not. ok) print *,' SEXT_ENERGY:  Qtune failed'
!    print '(1x, a27, 3f12.4)',' SEXT_ENERGY, Qx, Qy, Qz = ', &
!                    lat_1%a%tune/twopi, lat_1%b%tune/twopi, lat_1%z%tune/twopi


    init%vec(1:6) = 0.
    init%vec(1) = 0.0025
    init%vec(3) = 0.000025
    init%vec(6) = 0.002
    co_lat(0)%vec(:) = co%vec(:)

    call closed_orbit_calc(lat_1, co_lat, 6)
    
    orbit(0)%vec = co_lat(0)%vec + init%vec 
!   print '(1x,a23,6e12.4)',' closed orbit calc at ele 0 ', co_lat(0)%vec     
    max_x = abs(orbit(0)%vec(1))
    max_arc = 0.
    do i =1, con%qs_2qx_turns
      call track_all(lat_1, orbit)
      orbit(0)%vec = orbit(lat_1%n_ele_track)%vec
      max_arc = max(maxval(abs(orbit(1:lat_1%n_ele_track)%vec(1))),max_arc)
      max_x = max(abs(orbit(0)%vec(1)), max_x)
    end do

!    print '(1x,a22,f12.4,a12, f12.4)',' SEXT_ENERGY: max x = ',max_x,'  max_arc = ', max_arc
    max_x = max_x/init%vec(1)
    call fourier_comp(lat_1, res_cos, res_sin, res_amp)
    
  elseif( any(con%c(1:con%n_constraint)%variable == Qs_2Qx$))then
    call fourier_comp(lat, res_cos, res_sin, res_amp)
  endif
  
  return
  
end subroutine sext_energy