subroutine intq(T, k0, qmx, xstot)

      implicit none

      include 'const.inc'

      integer i, n
      parameter(n=200)
      double precision T, k0, qmx, xstot
      double precision rc, en0, cs, qmn, dq, sum, q, eq, w
      double precision ep, kp, th, 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)
*
      qmn = 0.
      dq  = (qmx - qmn)/dble(n)
      sum = 0.
      do 100 i = 1, n
         q  = qmn + dq*i
         eq = (h_par*q)**2/(2.*m_n)
         w  = q*cs
         ep = (en0-(h_par*w))
         kp = sqrt(2*m_n*ep)/h_par
         th = dacos((k0**2 + kp**2 - q**2)/(2.*k0*kp))
         if(th.lt.-1..or.th.gt.1..or.
     &        (k0**2 + kp**2 - q**2).gt.(2.*k0*kp) )
     &        write(*,*) '--nnnn-'
         dxs = ds0dw*dsin(th) /(1.+(rc*qmx)**2)
c         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*dq/denat
      return
      end