How to (correctly) rotate this data so the Ecliptic is in the XY plane?

Question

When plotting orbits using the Python package Skyfield, the coordinates are barycentric ICRS. This puts the center of mass of the solar system at the origin, but the $\hat{z}$ axis is tilted about 23° with respect to the normal to the ecliptic. I believe the ICRS $\hat{z}$ axis is coincident with declination of +90°, but not even sure of that.

What would be the most appropriate way to rotate the ICRS data so that the plane of the ecliptic (whatever that exactly means) coincides with the xy plane?

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

Answer 1

You should be able to fetch ecliptic coordinates with:

x, y, z = venus.at(t).ecliptic_position().km

Here is the current documentation — it is a bit brief, so I will try adding an example to them this weekend. You will see that there is also a method, in case you ever need it, for ecliptic latitude and longitude:

http://rhodesmill.org/skyfield/api-position.html#skyfield.positionlib.ICRF.ecliptic_position