def f1d2p(f, h, f0 = None):
    """
    A two-point finite difference stencil.
    This function does either two computations or one,
    depending on whether the 'center' value is supplied.
    This is done in order to avoid recomputing the center
    value many times when we repeat this function for each
    index of the gradient.

    How to use: use fdwrap or something similar to generate
    a one-variable function from the (usually) much more complicated
    function that we wish to differentate.  Then pass it to this function.

    Inputs:
    f  = The one-variable function f(x) that we're differentiating
    h  = The finite difference step size, usually a small number

    Outputs:
    fp = The finite difference derivative of the function f(x) around x=0.
    """
    if f0 is None:
        f0, f1 = [f(i*h) for i in [0, 1]]
    else:
        f1 = f(h)
    fp = (f1-f0)/h
    return fp