subroutine intq(T, k0, q, thmax, xstot)

      implicit none

      include 'const.inc'

      integer i, n
      parameter(n=200)
      double precision T, k0, q, xstot, eq, w
      double precision rc, en0, cs, thmin, thmax, dth, th, sum
      double precision ep, kp, dxs, ds0dw, hedens, denat
*
      rc  = (-4.7092+7.1047*T-1.9944*T**2+0.21769*T**3)/sqrt(3.)
      en0 = (h_par*k0)**2/(2.*m_n)
      cs  = 2.4170              ! (Ang/psec)
*
      call dsg0(T,ds0dw)
      call wq(q,w)
*
      eq = (h_par*q)**2/(2.*m_n)
      ep = (en0-(h_par*w))
      kp = sqrt(2*m_n*ep)/h_par
*
      sum   = 0.
      thmin = 0.
c      thmax = dacos((k0**2 + kp**2 - q**2)/(2.*k0*kp))
      dth   = (thmax-thmin)/dble(n)
      do 100 i = 1, n
         th  = thmin + dth*i
         dxs = ds0dw*dsin(th) /(1.+(rc*q)**2)
         sum = sum + dxs
 100  continue
*
      call dsg0(T,ds0dw)
      call heden(T,hedens)
      denat = hedens*avog/he4mol
*
      xstot = 2.*pi*sum*dth/denat
      return
      end