def METHOD_NAME(sl1, sl2, distance, color):
    # get distance in terms of unit cells
    d_cell = (distance + site_offsets[sl1] - site_offsets[sl2]) @ np.linalg.inv(
        basis_vectors
    )
    if not np.all(is_approx_int(d_cell, atol=atol)):
        # error out
        msg = f"{distance} is invalid distance vector between sublattices {sl1}->{sl2}"
        # see if the user flipped the vector accidentally
        d_cell = (distance + site_offsets[sl2] - site_offsets[sl1]) @ np.linalg.inv(
            basis_vectors
        )
        if np.all(is_approx_int(d_cell, atol=atol)):
            msg += f" (but valid {sl2}->{sl1})"
        raise ValueError(msg)
    d_cell = np.asarray(np.rint(d_cell), dtype=int)
    # catches self-referential and other unrealisable long edges
    if not np.all(d_cell < extent):
        raise ValueError(
            f"Distance vector {distance} does not fit into the lattice"
        )
    # Unit cells of starting points
    start_min = np.where(pbc, 0, np.maximum(0, -d_cell))
    start_max = np.where(pbc, extent, extent - np.maximum(0, d_cell))
    start_ranges = [slice(lo, hi) for lo, hi in zip(start_min, start_max)]
    start = np.mgrid[start_ranges].reshape(len(extent), -1).T
    end = (start + d_cell) % extent
    # Convert to site indices
    start = site_to_idx((start, sl1), extent, site_offsets)
    end = site_to_idx((end, sl2), extent, site_offsets)
    return [(*edge, color) for edge in zip(start, end)]