subroutine fourier_comp (lat, sum_c, sum_s, res_amp)

! compute the fourier component of 2Qx-Qs resonance with m=2Qx

  use bmad
!  use bmadz_interface
  
  implicit none
  
  type (lat_struct) lat
  type (ele_struct) ave, ele
  
  integer i, j, m, jmax
  
  real(rp) circ, s, k
  real(rp) theta, sint, cost
  real(rp) Qx
  real(rp) beta_x, alpha_x, rbz, eta
  real(rp) real_contrib, imag_contrib, sum_c, sum_s, res_amp
  logical ok

  sum_c=0.
  sum_s=0.
  
  Qx = lat%ele(lat%n_ele_track)%a%phi
  circ = lat%ele(lat%n_ele_track)%s
  m= 21 !(2*Qx-Qs)
  
  do i = 1, lat%n_ele_track
    ele = lat%ele(i)
    if (lat%ele(i)%key /= sextupole$ .and. lat%ele(i)%key /= quadrupole$ &
         .and. lat%ele(i)%key /= sol_quad$ .and. &
         lat%ele(i)%key /= wiggler$) cycle
    S=0.
    K=0.
    if (lat%ele(i)%key == sextupole$) &
         S = lat%ele(i)%value(k2$) * lat%ele(i)%value(l$) 
    if (lat%ele(i)%key == quadrupole$ .or. lat%ele(i)%key == sol_quad$) &
         k = lat%ele(i)%value(k1$) * lat%ele(i)%value(l$) 
    if (ele%key == wiggler$  .and. ele%name(1:5) /= 'WIG_W' &
         .and. ele%name(1:5) /= 'WIG_E') then

      ! Needed to add the following if statement after totalview showed
      ! examples where the term was unassociated and had a size of one,
      ! thus going through the loop once and dying when trying to access
      ! the data.         2006.01.24 mjf
      if (associated(ele%taylor(2)%term)) then
        jmax = size(ele%taylor(2)%term)
!        print *,' jmax ', jmax, ' element ',ele%name
        do j=1,jmax
          if (ele%taylor(2)%term(j)%expn(2) /= 0 .or. &
               ele%taylor(2)%term(j)%expn(3) /= 0) cycle
          if (ele%taylor(2)%term(j)%expn(4) /= 0 .or. &
               ele%taylor(2)%term(j)%expn(5) /= 0) cycle
          if (ele%taylor(2)%term(j)%expn(1) == 1 .and. &
               ele%taylor(2)%term(j)%expn(6) == 0) &
               k = -ele%taylor(2)%term(j)%coef
          if (ele%taylor(2)%term(j)%expn(1) == 2 .and. &
               ele%taylor(2)%term(j)%expn(6) == 0) &
               s = -2*ele%taylor(2)%term(j)%coef
        end do
      endif
    endif

!      print '(a12,4(a7,e12.4))',ele%name,'    s =',s,'    k =',k,'  s_t =',s_t,'  k_t  ',k_t
    call twiss_at_element(lat%ele(i), average = ave)
    beta_x = ave%a%beta 
    eta = ave%a%eta
    theta = m * lat%ele(i)%a%phi / Qx * twopi
    sint = sin(theta)
    cost = cos(theta)

!      real_contrib = (s*eta - k)* (beta_x)**1.5 * cost  
!      imag_contrib = (s*eta - k)* (beta_x)**1.5 * sint
    real_contrib = (s*eta - k)* (beta_x)**1 * cost  
    imag_contrib = (s*eta - k)* (beta_x)**1 * sint

    sum_c = sum_c + real_contrib
    sum_s = sum_s + imag_contrib
 
    !  print '(1x,a12,7f12.4)',lat%ele(i)%name,s,k,beta_x,eta,theta, real_contrib, imag_contrib

!    j = j+1
!      if(mod(j,20) == 0)pause



  end do
!     print '(a9, f12.4, a9, f12.4)','   A_m = ', sum_c,'   B_m = ', sum_s
  res_amp = sqrt(sum_c**2 + sum_s**2)
  return

end subroutine fourier_comp