def quadgk_int(f, fmin=-np.inf, fmax=np.inf, difftol=0.1):
    """
    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:
    """
    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