Help propagating Asteroid perturbed orbit

[Errabundo66]

I'm trying to get a perturbed orbit result for Asteroids using Poliastro.

I didn't find info about using Poliastro for minor planets bodies, so I'm not sure the code is correct.

For example, I take 1234 Elyna asteroid with osculating elements from MPC, dated 2023-02-25.

I just want to get the x,y,z coordinates (ecliptic frame) two years later, 2025-02-25, for:

• Non perturbated orbit
• Perturbated orbit with main planets
• JPL Horizons real position

The results are:

x, y, z, twobody keplerian Poliastro: [ 2.3806107 -1.65426043 0.15480046] AU
x, y, z, perturbed Poliastro: [ 2.37968736 -1.65574948 0.15441254] AU
x, y, z, JPL Horizons query: [ 2.37437146 -1.66096256 0.15495679] AU

Edit:
Just for comparing, I got the result using Skyfield for the unperturbed orbit with the same data:
x, y, z, Skyfield unperturbed:[ 2.39064251 -1.63976216 0.1546959 ] AU

What can I do to improve the result?
The perturbed orbit is not very different from the original one.

Another difference question is about Frames.
I'm not sure if by default Poliastro is using the same coordinates origine in JPL and in Orbit/Ephem system.

Any advise is welcome!
Thx

I attach the code:

from matplotlib import pyplot as plt
import numpy as np

from astropy.coordinates import solar_system_ephemeris
from astropy import time
from astropy import units as u

from poliastro.bodies import Earth, Jupiter, Sun, Mars, Venus, Saturn

from poliastro.core.elements import rv2coe
from poliastro.core.perturbations import (
    atmospheric_drag_exponential,
    third_body,
    J2_perturbation,
)
from poliastro.core.propagation import func_twobody
from poliastro.ephem import build_ephem_interpolant, Ephem
from poliastro.plotting import OrbitPlotter3D
from poliastro.twobody import Orbit
from poliastro.twobody.propagation import CowellPropagator
from poliastro.twobody.sampling import EpochsArray
from poliastro.util import norm
from poliastro.util import time_range
from poliastro.core import angles
from poliastro.frames import Planes

# More info: https://plotly.com/python/renderers/
import plotly.io as pio

pio.renderers.default = "plotly_mimetype+notebook_connected"


# database keeping positions of bodies in Solar system over time
solar_system_ephemeris.set("de432s")

def f(t0, state, k):
    du_kep = func_twobody(t0, state, k)

    ax, ay, az = third_body(
        t0,
        state,
        k,
        k_third = Jupiter.k.to(u.km**3 / u.s**2).value,
        perturbation_body=body_jupiter,
    )
    du_ju = np.array([0, 0, 0, ax, ay, az])

    ax, ay, az = third_body(
        t0,
        state,
        k,
        k_third = Earth.k.to(u.km**3 / u.s**2).value,
        perturbation_body=body_earth,
    )
    du_ea = np.array([0, 0, 0, ax, ay, az])

    ax, ay, az = third_body(
        t0,
        state,
        k,
        k_third = Mars.k.to(u.km**3 / u.s**2).value,
        perturbation_body=body_mars
    )
    du_ma = np.array([0, 0, 0, ax, ay, az])

    ax, ay, az = third_body(
        t0,
        state,
        k,
        k_third = Venus.k.to(u.km**3 / u.s**2).value,
        perturbation_body=body_venus
    )
    du_ve = np.array([0, 0, 0, ax, ay, az])

    ax, ay, az = third_body(
        t0,
        state,
        k,
        k_third = Saturn.k.to(u.km**3 / u.s**2).value,
        perturbation_body=body_saturn
    )
    du_sa = np.array([0, 0, 0, ax, ay, az])

    return du_kep + du_ju + du_ea + du_ma + du_ve + du_sa


# Elementos keplerianos del asteroide 1234 Elyna
inc = 8.52628 << u.deg # Inclinación en grados
raan = 304.33904 << u.deg # Ascensión recta del nodo ascendente en grados
ecc = 0.0869957 << u.one # Excentricidad
argp = 90.41267 << u.deg # Argumento del pericentro en grados
a = 3.0150261 << u.AU # Semieje mayor en UA
mu= 162.15207 << u.deg # Anomalía media en grados

# Calculate True Anomaly nu
Ecc_anomaly= angles.M_to_E(np.radians(mu), ecc)
nu= np.degrees(angles.E_to_nu(Ecc_anomaly, ecc)) * u.deg

#MPCORB date
date_mpc = time.Time("2023-02-25 00:00:00", scale="utc")

#Oher dates... date for checking and date_end final date period calculations
date = time.Time("2025-02-25 00:00:00", scale="utc")
date_end = time.Time("2026-02-25 23:00:00", scale="utc")

#Create Orbit from COE in Ecliptic frame
asteroid_orbit = Orbit.from_classical(Sun, a, ecc, inc, raan, argp, nu, date_mpc, plane=Planes.EARTH_ECLIPTIC)

#Check x,y,z at date time without perturbations
ephem2 = Ephem.from_orbit(asteroid_orbit, date, plane=Planes.EARTH_ECLIPTIC)
rv1=ephem2.rv(date)
print("x, y, z, twobody keplerian Poliastro:",rv1[0].to(u.AU))


# create interpolant of 3rd body coordinates (calling in on every iteration will be just too slow)
body_jupiter = build_ephem_interpolant(
    Jupiter,
    4332.615 * u.day,
    (date_mpc, date_end),
    rtol=1e-2,
    attractor=Sun
)

body_earth = build_ephem_interpolant(
    Earth,
    365 * u.day,
    (date_mpc, date_end),
    rtol=1e-2,
    attractor=Sun
)

body_mars = build_ephem_interpolant(
    Mars,
    687 * u.day,
    (date_mpc, date_end),
    rtol=1e-2,
    attractor=Sun
)

body_venus = build_ephem_interpolant(
    Venus,
    225 * u.day,
    (date_mpc, date_end),
    rtol=1e-2,
    attractor=Sun
)

body_saturn = build_ephem_interpolant(
    Saturn,
    10759 * u.day,
    (date_mpc, date_end),
    rtol=1e-2,
    attractor=Sun
)

tofs = time.TimeDelta(np.linspace(0, 800 * u.day, num=1000))


# Calculate the perturbated Orbit between date_mpc and date_end
ephem = asteroid_orbit.to_ephem(
    EpochsArray(asteroid_orbit.epoch + tofs, method=CowellPropagator(rtol=1e-6, f=f)),
)

#Check x,y,z at date time with Perturbations 
rv=ephem.rv(date)
print("x, y, z, perturbated Poliastro:",rv[0].to(u.AU))

#Check with real JPL Horizons coordinates
Elyna= Orbit.from_sbdb("Elyna")
date = time.Time("2025-02-25 00:00:00", scale="utc")
Elyna_eph = Ephem.from_horizons("Elyna", date, plane=Planes.EARTH_ECLIPTIC)

rv2=Elyna_eph.rv(date)
print("x, y, z, JPL Horizons query:",rv2[0].to(u.AU))

-->

🖥 Please paste the output of following commands

• conda info -a (only if you have conda)
• conda list (only if you have conda)
• pip freeze

# Paste your output here:

# Paste your output here:

💡 Possible solutions

📋 Steps to solve the problem

• Comment below about what you've started working on.
• Add, commit, push your changes
• Submit a pull request and add this in comments - Addresses #<put issue number here>
• Ask for a review in comments section of pull request
• Celebrate your contribution to this project 🎉