code: def quadgk_int(f, fmin=-np.inf, fmax=np.inf, difftol=0.1):
      
      difftol = 1e-4
      def trans_func(r):
          r2 = np.square(r)
          x = r / (1-r2)
          dx_dr = (1 + r2)/(1-r2)**2
          return f(x)*dx_dr
  
      integrand = quadvgk(trans_func, -1., 1., difftol, difftol)
      return integrand
docs: Integrate f from fmin to fmax,
      do integration by substitution
      x = r / (1-r**2)
      when r goes from -1 to 1 , x goes from -inf to inf.
      the interval for quadgk function is from -1 to +1, so we transform the space from (-inf,inf) to (-1,1)
      :param f:
      :param fmin:
      :param fmax:
      :param difftol:
      :return: