CODEWITHSUNDEEP
pythonastronomyskyfield
Fabio
Fabio

Reputation: 85

How to calculate distance to comets using Skyfields?

I'm using skyfield to compute the relative distance of planets to Earth as a function of time (as described on the skyfield home page). It works great and now I'm trying to implement Earth=>comet distance (e.g. 67P/ Tchouri).

I've found at NASA JPL, a way to create Spice SPK files for comets (here) but it produces xsp files that I cannot seem to read with the load command from skyfield.

Another possibility I considered is to use orbital information as suggested for pyephem (see here) but I don't know how to read them in Skyfield.

I also saw that comets were on the roadmap for skyfield coding sprint so that's maybe my answer but if you know a way to make it work with the current version that would be very helpful.

Upvotes: 1

Views: 323

Answers (1)

Brandon Rhodes
Brandon Rhodes

Reputation: 89415

Skyfield did gain support for comets! You can find the details here:

https://rhodesmill.org/skyfield/kepler-orbits.html#comets

Adapting the code from the documentation, here is the distance to a comet from the Minor Planet Center database:

from skyfield.api import load
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN
from skyfield.data import mpc

ts = load.timescale()

eph = load('de421.bsp')
sun, earth = eph['sun'], eph['earth']

with load.open(mpc.COMET_URL) as f:
    comets = mpc.load_comets_dataframe(f)
comets = comets.set_index('designation', drop=False)
row = comets.loc['1P/Halley']
comet = sun + mpc.comet_orbit(row, ts, GM_SUN)

t = ts.utc(2020, 10, 17)
ra, dec, distance = earth.at(t).observe(comet).radec()

print('Distance in AU:', distance.au)

I see the result:

Distance in AU: 35.22790002485247

Upvotes: 1

Related Questions