Как (правильно) повернуть эти данные, чтобы эклиптика находилась в плоскости XY?

При построении орбит с использованием пакета Skyfield для Python используются барицентрические координаты ICRS . Это помещает центр масс Солнечной системы в начало координат, но$\hat{z}$ось наклонена примерно на 23° по отношению к нормали к эклиптике. Я верю ICRS$\hat{z}$ось совпадает со склонением +90°, но даже не уверен в этом.

Как лучше всего повернуть данные ICRS так, чтобы плоскость эклиптики (что бы это ни значило ) совпадала с плоскостью 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()