halo_mass_function.py 文件源码

def integrate_inf(fcn, error=1e-3, initial_guess=10):
    """
    Integrate a function *fcn* from zero to infinity, stopping when the answer
    changes by less than *error*. Hopefully someday we can do something
    better than this!
    """
    xvals = np.logspace(0,np.log10(initial_guess), initial_guess+1)-.9
    yvals = fcn(xvals)
    xdiffs = xvals[1:] - xvals[:-1]
    # Trapezoid rule, but with different dxes between values, so np.trapz
    # will not work.
    areas = (yvals[1:] + yvals[:-1]) * xdiffs / 2.0
    area0 = np.sum(areas)
    # Next guess.
    next_guess = 10 * initial_guess
    xvals = np.logspace(0,np.log10(next_guess), 2*initial_guess**2+1)-.99
    yvals = fcn(xvals)
    xdiffs = xvals[1:] - xvals[:-1]
    # Trapezoid rule.
    areas = (yvals[1:] + yvals[:-1]) * xdiffs / 2.0
    area1 = np.sum(areas)
    # Now we refine until the error is smaller than *error*.
    diff = area1 - area0
    area_last = area1
    one_pow = 3
    while diff > error:
        next_guess *= 10
        xvals = np.logspace(0,np.log10(next_guess), one_pow*initial_guess**one_pow+1) - (1 - 0.1**one_pow)
        yvals = fcn(xvals)
        xdiffs = xvals[1:] - xvals[:-1]
        # Trapezoid rule.
        areas = (yvals[1:] + yvals[:-1]) * xdiffs / 2.0
        area_next = np.sum(areas)
        diff = area_next - area_last
        area_last = area_next
        one_pow+=1
    return area_last