subroutine sext_detune (lat, axx,axy,ayy)

! compute the detuning coefficients axx, axy, ayy from sextupoles
!  from S.Y. Lee page 191

  use bmad
!  use bmad_struct
!  use bmad_interface

  implicit none

  type (lat_struct) lat
  type (ele_struct) ele, ele2

  real(rp) axx, axy, ayy
  real(rp) pi_nux, pi_nuy, den_3, den_1, den_2p, den_2m
  real(rp) rbxi, rbyi, rbxj, rbyj
  real(rp) S_i, S_j
  real(rp) dps_x, dps_y
  real(rp) cos3x, cos1x, cos2xy_p, cos2xy_m
  real(rp) sum_xx, sum_xy, sum_yy
 
  integer i,j


  pi_nux = 0.5*lat%a%tune
  pi_nuy = 0.5*lat%b%tune
  den_3 = 1./sin(3*pi_nux)
  den_1 = 1./sin(pi_nux)
  den_2p = 1./sin(2*pi_nuy + pi_nux)
  den_2m = 1./sin(2*pi_nuy - pi_nux)

  sum_xx=0.
  sum_xy=0.
  sum_yy=0.

  do i = 1, lat%n_ele_track
  if(lat%ele(i)%key == sextupole$)then
     call twiss_at_element (lat%ele(i), average = ele)
     rbxi = sqrt(ele%a%beta) 
     rbyi = sqrt(ele%b%beta)
     S_i = ele%value(k2$) * lat%ele(i)%value(l$)

     do j = 1, lat%n_ele_track
     if(lat%ele(j)%key == sextupole$)then
        call twiss_at_element (lat%ele(j), average = ele2)
        
        dps_x = ele2%a%phi - ele%a%phi
        if(dps_x < 0)dps_x = dps_x + 2*pi_nux
        dps_y = ele2%b%phi - ele%b%phi
        if(dps_y < 0)dps_y = dps_y + 2*pi_nuy
        
        cos3x = cos(3*(pi_nux - dps_x))
        cos1x = cos((pi_nux - dps_x))
        cos2xy_p = cos(2*(pi_nuy - dps_y) + pi_nux - dps_x)
        cos2xy_m = cos(2*(pi_nuy - dps_y) - pi_nux + dps_x)

        rbxj = sqrt(ele2%a%beta) 
        rbyj = sqrt(ele2%b%beta) 
        S_j = ele2%value(k2$)*ele2%value(l$)     


        sum_xx = sum_xx + S_i*S_j*(rbxi*rbxj)**3 &
                       * (cos3x * den_3 + 3* cos1x * den_1)
        sum_xy = sum_xy + S_i*S_j*(rbxi*rbxj)*(rbyi*rbyj)**2 &
                       * (cos2xy_p * den_2p +cos2xy_m * den_2m) &
                       -2 * S_i*S_j*rbxi**3*rbxj*rbyj**2 * cos1x * den_1
        sum_yy = sum_yy + S_i*S_j*rbxi*rbxj*(rbyi*rbyj)**2 &
           * (cos2xy_p * den_2p + cos2xy_m * den_2m + 3* cos1x * den_1)

      endif
      end do

    endif
    end do

    axx=1/32./twopi * sum_xx
    axy = 1/16./twopi * sum_xy
    ayy = 1./32./twopi * sum_yy

  return

end subroutine