Return separation vectors between each pair of the two sets of points.

    Parameters
    ----------
    ra, rb: float array-like, shape (n, d) and (m, d) in d dimensions.
        Two sets of points.
    L: float array, shape (d,)
        System lengths.

    Returns
    -------
    csep: float array-like, shape (n, m, d)
        csep[i, j] is the separation vector from point j to point i.
        Note the un-intuitive vector direction.

def csep_periodic(ra, rb, L):
    """Return separation vectors between each pair of the two sets of points.

    Parameters
    ----------
    ra, rb: float array-like, shape (n, d) and (m, d) in d dimensions.
        Two sets of points.
    L: float array, shape (d,)
        System lengths.

    Returns
    -------
    csep: float array-like, shape (n, m, d)
        csep[i, j] is the separation vector from point j to point i.
        Note the un-intuitive vector direction.
    """
    seps = ra[:, np.newaxis, :] - rb[np.newaxis, :, :]
    for i_dim in range(ra.shape[1]):
        seps_dim = seps[:, :, i_dim]
        seps_dim[seps_dim > L[i_dim] / 2.0] -= L[i_dim]
        seps_dim[seps_dim < -L[i_dim] / 2.0] += L[i_dim]
    return seps