La órbita del sol y la luna en coordenadas eclípticas usando skyfield

No puedo ayudarte con Skyfield, pero normalmente uso la interfaz web de JPL Horizons. No requiere instalación, también puede imprimirlo en un archivo de texto si lo desea:

https://ssd.jpl.nasa.gov/horizons.cgi

De lo contrario, encontré la documentación para skyfield: https://rhodesmill.org/skyfield/toc.html

Y si nada de eso funciona, hice un pequeño guión de astrología para ti.

from astropy import units as u
from astropy.coordinates import SkyCoord, EarthLocation, AltAz, get_body
from astropy.time import Time
import numpy as np

# Create 1000 Timepoints between Time 1 and Time 2 (one year later)
t = np.linspace(2451545, 2451545+365, 1000)


pointlist = []

#Loop through this times
for tn in t:
    # For every timepoint, create an astropy_time object
    astropy_time = Time(tn, format="jd")
    # Get Planet (as string, "earth", "moon", "mercury" etc. in aequatorial coordinates
    planet_aequatorial = get_body("moon", time = astropy_time)
    #Transform to Barycentric True Ecliptic (relative to the center of mass of the solar system).
    planet_ecliptic =  planet_aequatorial.transform_to("barycentrictrueecliptic")

    # Add a point to the orbit. Every point is described as (longitude [deg], latitude (ecliptic coords), distance (km))
    pointlist.append([planet_ecliptic.lon.deg, planet_ecliptic.lat.deg, planet_ecliptic.distance.km])
    print(planet_ecliptic.distance.km)
    # So pointslist is a 2D array. The rows are all the 1000 points of the orbit

# In every point there is 3 columns for [Long, Lat, Distance]
print(pointlist)


# You can also save the result with
pointlist = np.array(pointlist)

np.save("resultados.npy", lista de puntos)