Wie kann man diese Daten (richtig) drehen, damit sich die Ekliptik in der XY-Ebene befindet?

Beim Zeichnen von Umlaufbahnen mit dem Python-Paket Skyfield sind die Koordinaten baryzentrische ICRS . Damit liegt der Massenschwerpunkt des Sonnensystems im Ursprung, aber der$\hat{z}$Die Achse ist gegenüber der Ekliptiknormalen um 23° geneigt. Ich glaube dem ICRS$\hat{z}$Achse fällt mit der Deklination von +90° zusammen, aber nicht einmal dessen ist sicher.

Was wäre der geeignetste Weg, um die ICRS-Daten so zu drehen, dass die Ebene der Ekliptik (was auch immer das genau bedeutet) mit der xy-Ebene zusammenfällt?

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()