Listing 1: FORTRAN code to calculate [6]

         rk = dlog(k1/k2)
         dk = k1-k2
         eu = dexp(-u)
         if (dabs(dk)/k2 .le. 1.d-4) then
            if (u .le. 0.0036d0) then
               wv = u/6.d0
            else if (u .ge. 36.d0) then
               wv = (u-2.d0)/u
            else
               wv = (u-2.d0+(u+2.d0)*eu)/(u*(1.d0-eu))
            end if
        else
            if (u .le. 0.0036d0) then
               wv = (dk*(1.d0+0.5d0*rk)-k1*rk)*u/(-rk*(dk-k2*rk))
            else if (u .ge. 36.d0) then
               wv = 1.d0 - k2*rk*rk/(u*(k1-k2*(1.d0+rk)))
            else if (dabs(u-rk)/u .le. 1.d-4) then
               wv = (dk*dk-k1*k2*rk*rk)/(dk*(dk-k2*rk))
            else 
               wv = ( (1.d0-eu)*dk*(u-rk) +(k1*eu-k2)*rk*u) /
     +               ( (1.d0-eu)*((dk-k2*rk)*u-dk*rk+k2*rk*rk) )
            end if
        end if