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

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

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

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

• That z-axis was Earth's north axis on January 1, 2000. Earth's axis has been moving about 20 arcseconds per year since then. The ecliptic is the plane of Earth's orbit around the Sun (or the solar system barycenter), usually as of January 1, 2000 as well. That plane rotates about 0.5 arcseconds per year. May 6 '16 at 14:22

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