sample_one_sphere

def shift_bit_length(n):
    return 1 << (n - 1).bit_length()

def sample_one_sphere(coords, radius, n, sobol=True):
    if sobol:
        m = int(np.log2(shift_bit_length(n)))
        x = norm.ppf(qmc.Sobol(d=3, scramble=True).random_base2(m))
    else:
        x = np.random.normal(size=(n, 3))
    ssq = np.sum(x**2, axis=1)
    fr = radius * gammainc(3 / 2, ssq / 2) ** (1 / 3) / np.sqrt(ssq)
    frtiled = np.tile(fr.reshape(n, 1), (1, 3))
    p = coords + np.multiply(x, frtiled)
    return p