code: def _forceInt(x,y,z,a2,b2,c2,n,i):
      
      def integrand(tau):
          return (x*(i==0) + y*(i==1) + z*(i==2))/(a2*(i==0) + b2*(i==1) + c2*(i==2) + tau) * \
              _FracInt(x, y, z, a2, b2, c2, tau, n)
      return integrate.quad(integrand, lowerlim(x**2, y**2, z**2, a2, b2, c2), np.inf, epsabs=1e-12)[0]
docs: Integral involved in the force at (x,y,z)
      integrates 1/A B^n (x_i/(tau+a_i)) where
      A = sqrt((tau+a)(tau+b)(tau+c)) and B = (1-x^2/(tau+a)-y^2/(tau+b)-z^2/(tau+c))
      from lambda to infty with respect to tau.
      The lower limit lambda is given by lowerlim function.