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