¿Cómo rotar (correctamente) estos datos para que la eclíptica esté en el plano XY?

Cuando se trazan órbitas con el paquete Skyfield de Python , las coordenadas son ICRS baricéntricas . Esto pone el centro de masa del sistema solar en el origen, pero el$\hat{z}$El eje está inclinado unos 23° con respecto a la normal a la eclíptica. Creo que el ICRS$\hat{z}$el eje coincide con la declinación de +90°, pero ni siquiera estoy seguro de eso.

¿Cuál sería la forma más apropiada de rotar los datos del ICRS para que el plano de la eclíptica (lo que sea que eso signifique exactamente ) coincida con el plano xy?

from skyfield.api import load
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

eph = load('de421.bsp')

venus = eph['venus']
earth = eph['earth']
mars  = eph['mars barycenter']

ts = load.timescale()

t  = ts.utc(2016, 1, range(700))

vp = venus.at(t).position.km
ep = earth.at(t).position.km
mp = mars.at(t).position.km

ps = [vp, ep, mp]
hw = max([np.abs(p).max() for p in ps])
size = [-hw, hw]

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection ='3d')
for p in ps:
    x, y, z = p
    ax.plot(x, y, z)
ax.set_xlim(size), ax.set_ylim(size), ax.set_zlim(size)
plt.show()