def METHOD_NAME(npt, ndim, rstate=None):
    """ Simulate points in ndim ball """
    # use Barthe2005
    x = rstate.standard_normal(size=(npt, ndim))
    y = rstate.exponential(0.5, size=npt)
    x1 = x / np.sqrt((y + (x**2).sum(axis=1)))[:, None]
    return x1