# Sun and moon's orbit in ecliptic coordinates using skyfield

I am new to using skyfield, is there any doc or help file that can show me on how to get the orbit of Sun and moon in ecliptic coordinate for a particular date and time. This is a follow up question of this question

• I don't really see how this question is different from the one you link to. Nov 28 '19 at 10:21
• @usernumber that was how can be done, this is for how can be done using skyfield Nov 28 '19 at 10:26

I cant help you with skyfield, but i usually use JPL Horizons Web interface. No installation required, you can also print it in a text file if you want:

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

Otherwise i found the documentation for skyfield: https://rhodesmill.org/skyfield/toc.html

And if nothing of that works, i made a little astropy script for you

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("results.npy", pointlist)

• thanks or the reply, it gives me a single point, however i want complete orbit, is that possible, my question may look silly, apologise Nov 29 '19 at 17:44
• create a list of times (in julian days) with >import numpy as np >t = np.linspace(2451545, 2451545+365, 2000) then loop over the above code but create Time object with >t = Time(t1, format="jd") instead Dec 9 '19 at 22:47
• did that help or should i do it for you? Jan 3 '20 at 13:01
• yes, please. that will be really helpful Jan 18 '20 at 19:01
• can you give me your e-mail or something? i cant post all the code here i think... Jan 23 '20 at 11:11