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
*
c   Using model from Tsipenyuk, S(q) =S(0)/(1+(rc*q)**2)
c
      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,en0,k0)
*
      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.
cc     thmax = dacos((k0**2 + kp**2 - q**2)/(2.*k0*kp))   bec
c      dth   = (thmax-thmin)/dble(n)
c      do 100 i = 1, n
c         th  = thmin + dth*i
c         dxs = ds0dw*dsin(th) /(1.+(rc*q)**2)
c         sum = sum + dxs
c 100  continue
*
cc      call dsg0(T,ds0dw)
      call heden(T,hedens)
      denat = hedens*avog/he4mol
*
c      xstot = 2.*pi*sum*dth/denat
c		We want probabilty of scattering, found from interaction length (dreact) found below.  
c         This depends on the total cross section.  The angular dependence is taken care of 
c		in choosing q from S(q).  To try and include some angular dependence here as well, 
c		would be redundant.  (see dsg0.f for same modification)
c		Another way to see this is to realize that integrating theta from 0 to theta_sc 
c		biases larger angles (larger integration), which is no sensical.   --bec
c
       xstot = 4.*pi*ds0dw/(1.+(rc*q)**2)/denat  
      return
      end