def rij_by_rij3(ri, rj):
    return (ri - rj) / np.linalg.norm(ri - rj) ** 3


def three_body_system(X, t, M, k1, k2):
    x_v_vectors = X.reshape(
        2, 3, 3
    )  # dim-0 contains position vecs and dim-1 contains velocity vecs
    x_v_deriv = np.empty(x_v_vectors.shape)

    n_objects = np.shape(M)[0]
    obj_indices = range(n_objects)

    # loop over each object and compute velocity and acceleration
    for i, m_i in enumerate(M):
        # acceleration
        rest_of = np.delete(obj_indices, i)
        x_v_deriv[1, i, :] = k1 * np.sum(
            np.array(
                [
                    [M[i] * rij_by_rij3(x_v_vectors[0, i, :] - x_v_vectors[0, j, :])]
                    for j in rest_of
                ]
            ),
            axis=0,
        )

        # velocity
        x_v_deriv[0, i, :] = k2 * x_v_vectors[1, i, :]

        return np.ravel(x_v_deriv)