Pack a number of non-intersecting spheres into a system.

    Can specify packing by number of spheres or packing fraction.

    This implementation uses a naive uniform distribution of spheres,
    and the Tabula Rasa rule (start from scratch if an intersection occurs).

    This is likely to be very slow for high packing fractions

    Parameters
    ----------
    R: float
        Sphere radius.
    L: float array, shape (d,)
        System lengths.
    pf: float or None
        Packing fraction
    n: integer or None
        Number of spheres.
    rng: RandomState or None
        Random number generator. If None, use inbuilt numpy state.
    periodic: bool
        Whether or not the system is periodic.

    Returns
    -------
    r: float array, shape (n, d)
        Coordinates of the centres of the spheres for a valid configuration.
    R_actual: float
        Actual sphere radius used in the packing.
        In this implementation this will always be equal to `R`;
        it is returned only to provide a uniform interface with the
        Metropolis-Hastings implementation.

def pack_simple(R, L, pf=None, n=None, rng=None, periodic=False):
    """Pack a number of non-intersecting spheres into a system.

    Can specify packing by number of spheres or packing fraction.

    This implementation uses a naive uniform distribution of spheres,
    and the Tabula Rasa rule (start from scratch if an intersection occurs).

    This is likely to be very slow for high packing fractions

    Parameters
    ----------
    R: float
        Sphere radius.
    L: float array, shape (d,)
        System lengths.
    pf: float or None
        Packing fraction
    n: integer or None
        Number of spheres.
    rng: RandomState or None
        Random number generator. If None, use inbuilt numpy state.
    periodic: bool
        Whether or not the system is periodic.

    Returns
    -------
    r: float array, shape (n, d)
        Coordinates of the centres of the spheres for a valid configuration.
    R_actual: float
        Actual sphere radius used in the packing.
        In this implementation this will always be equal to `R`;
        it is returned only to provide a uniform interface with the
        Metropolis-Hastings implementation.
    """
    if rng is None:
        rng = np.random
    if pf is not None:
        if pf == 0.0:
            return np.array([]), R
        # If packing fraction is specified, find required number of spheres
        # and the actual packing fraction this will produce
        n, pf_actual = pf_to_n(L, pf, R)
    elif n is not None:
        if n == 0:
            return np.array([]), R
    if periodic:
        return _pack_simple_periodic(R, L, n, rng)
    else:
        return _pack_simple(R, L, n, rng)