(Un-)expected inaccuracies when propagating heliocentric orbits (?)

[s-m-e]

I am exploring poliastro's propagators with Minor Planet Center "classical" orbital elements and poliastro's bodies vs. JPL Horizons' ephemerides data. I am not entirely sure if the results are within the expected accuracy limits or if there is a problem in poliastro or if I have not read the manual properly. If I had to guess I'd says it's the first given poliastro's focus on two-body propagators.

I have attached a set of three short notebooks illustrating my experiments: poliastro_debug.tar.gz

🐞 Problem

In experiment 1, I compared Eros' orbit based on MPC data and JPL data:

I know MPC's data is not as good as JPL's computations, but this seems a little too far off ...

The relevant section(s) from my code look about as follows:

```python eros_mpc = { "Epoch": 2458700.5, "M": 103.1926, "Peri": 178.84531, "Node": 304.30277, "i": 10.82872, "e": 0.2227412, "n": 0.55971324, "a": 1.4582288, } def orbit_from_mpc(data): return Orbit.from_classical( Sun, a = (data['a'] * u.AU).to(u.km), ecc = data['e'] * u.one, inc = (data['i'] * u.deg).to(u.rad), raan = (data['Node'] * u.deg).to(u.rad), argp = (data['Peri'] * u.deg).to(u.rad), nu = (data['M'] * u.deg).to(u.rad), epoch = Time(data["Epoch"], format = 'jd'), plane = Planes.EARTH_ECLIPTIC, ) eros_mpc_orb = orbit_from_mpc(eros_mpc) astro_times = Time([Time(Time('2020-01-01 00:00') + timedelta(days = day)) for day in range(0, 31, 1)]) relative_times = astro_times - eros_mpc_orb.epoch eros_mpc_eph = propagate(eros_mpc_orb, relative_times).xyz.to(u.AU).value # plotted eros_jpl_eph = Ephem.from_horizons( name = 'Eros', epochs = astro_times, attractor = Sun, plane = Planes.EARTH_ECLIPTIC, ).sample().xyz.to(u.AU).value # plotted ```

EDIT: MPC provides mean anomalies, Orbit.from_classical expects true anomalies. A quick & dirty solution looks as follows:

```python from orbital.utilities import true_anomaly_from_mean def orbit_from_mpc(data): return Orbit.from_classical( Sun, a = (data['a'] * u.AU).to(u.km), ecc = data['e'] * u.one, inc = (data['i'] * u.deg).to(u.rad), raan = (data['Node'] * u.deg).to(u.rad), argp = (data['Peri'] * u.deg).to(u.rad), nu = true_anomaly_from_mean(e = data['e'], M = float((data['M'] * u.deg).to(u.rad).value)) * u.rad, # FIX epoch = Time(data["Epoch"], format = 'jd'), plane = Planes.EARTH_ECLIPTIC, ) ```

In experiments 2 and 3, I looked at Earth's orbit around the Sun.

Computations are based on an example in poliastro's manual where an Orbit object can be generated via an intermediate step through an Ephem object from the Earth object. In the above plot, I essentially varied the epoch when running Ephem.from_body.

The relevant section(s) from my code look about as follows:

```python astro_times = Time([Time(Time('2020-01-01 00:00').datetime + timedelta(days = day)) for day in range(0, 120, 30)]) EARTH_EPOCH = Time("2010-01-01 00:00") # some "random" time(s) earth_orb = Orbit.from_ephem( attractor = Sun, ephem = Ephem.from_body( body = Earth, epochs = EARTH_EPOCH.tdb, plane = Planes.EARTH_ECLIPTIC, attractor = Sun, ), epoch = EARTH_EPOCH.tdb, ) relative_times = astro_times - earth_orb.epoch earth_poliastro_eph = propagate(earth_orb, relative_times).xyz.to(u.AU).value # plotted earth_jpl_eph = Ephem.from_horizons( name = '399', epochs = astro_times, attractor = Sun, plane = Planes.EARTH_ECLIPTIC, id_type = 'majorbody', ).sample().xyz.to(u.AU).value # plotted ```

Choosing different values for EARTH_EPOCH yields vastly different results. Interestingly, it's not so much the difference between astro_times and EARTH_EPOCH that appears to dominate the inaccuracies. At least within a decade of difference, it appears to be the relative position of Earth on its orbit (i.e. "how far into the year" EARTH_EPOCH is) that dominates the variations.

🖥 Please paste the output of following commands

• pip freeze

``` ads==0.12.3 appdirs==1.4.4 argon2-cffi==20.1.0 asciitree==0.3.3 astroid==2.4.2 astropy==4.0.1.post1 astropy-helpers==4.0.1 astroquery==0.4.1 attrs==19.3.0 autopep8==1.5.4 backcall==0.2.0 beautifulsoup4==4.9.1 -e git+github:pleiszenburg/bewegung.git@7d0cc4e23c348576deace98a2ba76afeac39ef93#egg=bewegung black==19.10b0 bleach==3.1.5 bokeh==2.1.1 certifi==2020.6.20 cffi==1.14.2 chardet==3.0.4 click==7.1.2 cloudpickle==1.5.0 colorama==0.4.3 colorcet==2.0.2 cryptography==3.1 cycler==0.10.0 Cython==0.29.21 dask==2.22.0 datashader==0.11.0 datashape==0.5.2 decorator==4.4.2 defusedxml==0.6.0 distributed==2.22.0 docutils==0.16 entrypoints==0.3 fasteners==0.15 flake8==3.8.3 fsspec==0.8.0 HeapDict==1.0.1 html5lib==1.1 httpretty==0.8.10 idna==2.10 iniconfig==1.0.1 ipykernel==5.3.4 ipympl==0.5.7 ipython==7.17.0 ipython-genutils==0.2.0 ipywidgets==7.5.1 isort==4.3.21 jedi==0.17.2 jeepney==0.4.3 Jinja2==2.11.2 jplephem==2.15 json5==0.9.5 jsonschema==3.2.0 jupyter-client==6.1.6 jupyter-core==4.6.3 jupyterlab==2.2.6 jupyterlab-server==1.2.0 keyring==21.4.0 kiwisolver==1.2.0 lazy-object-proxy==1.4.3 llvmlite==0.34.0 locket==0.2.0 MarkupSafe==1.1.1 matplotlib==3.3.1 mccabe==0.6.1 mimeparse==0.1.3 mistune==0.8.4 mock==4.0.2 monotonic==1.5 more-itertools==8.4.0 msgpack==1.0.0 multipledispatch==0.6.0 nbconvert==5.6.1 nbformat==5.0.7 nodeenv==1.4.0 notebook==6.1.1 numba==0.51.2 numcodecs==0.6.4 numpy==1.19.1 packaging==20.4 pandas==1.1.1 pandocfilters==1.4.2 param==1.9.3 parso==0.7.1 partd==1.1.0 pathspec==0.8.0 pexpect==4.8.0 pickleshare==0.7.5 Pillow==7.2.0 pkginfo==1.5.0.1 plotly==4.9.0 pluggy==0.13.1 poliastro==0.14.0 prometheus-client==0.8.0 prompt-toolkit==3.0.6 psutil==5.7.2 ptyprocess==0.6.0 py==1.9.0 pycairo==1.19.1 pycodestyle==2.6.0 pycparser==2.20 pyct==0.4.6 pydocstyle==5.0.2 pyflakes==2.2.0 Pygments==2.6.1 PyGObject==3.36.1 pylint==2.5.3 pymongo==3.11.0 pyparsing==2.4.7 pyrsistent==0.16.0 pytest==6.0.1 python-dateutil==2.8.1 python-jsonrpc-server==0.3.4 python-language-server==0.34.1 pytz==2020.1 pyvo==1.1 PyYAML==5.3.1 pyzmq==19.0.2 readme-renderer==26.0 regex==2020.7.14 requests==2.24.0 requests-toolbelt==0.9.1 retrying==1.3.3 rfc3986==1.4.0 rope==0.17.0 sbpy==0.2.2 scipy==1.5.2 seaborn==0.10.1 SecretStorage==3.1.2 Send2Trash==1.5.0 six==1.15.0 snowballstemmer==2.0.0 sortedcontainers==2.2.2 soupsieve==2.0.1 synphot==0.1.3 tblib==1.7.0 terminado==0.8.3 testpath==0.4.4 toml==0.10.1 toolz==0.10.0 tornado==6.0.4 tqdm==4.48.2 traitlets==4.3.3 twine==3.2.0 typed-ast==1.4.1 typeguard==2.9.1 typing-extensions==3.7.4.2 ujson==1.35 urllib3==1.25.10 wcwidth==0.2.5 webencodings==0.5.1 Werkzeug==1.0.1 widgetsnbextension==3.5.1 wrapt==1.12.1 xarray==0.16.0 yapf==0.30.0 zarr==2.4.0 zict==2.0.0 ```

• python -c "import poliastro.testing; poliastro.testing.test()"

ModuleNotFoundError: No module named 'poliastro.testing'

[astrojuanlu]

Hi @s-m-e! Just a quick question while I find some spare time to closely look at this: the blue and green line in your first plot overlap, correct?

[s-m-e]

Correct, they do - at least at this zoom level.

[s-m-e]

Take your time ... whenever I end up writing bug reports, I find (at least part of) the problem at my end. MPC provides mean anomalies, Orbit.from_classical expects true anomalies. At least this part of the result looks much better now.

[no-response[bot]]

This issue has been automatically closed because there has been no response to our request for more information from the original author. With only the information that is currently in the issue, we don't have enough information to take action. Moreover, the problems discussed in the issue may not be directly from poliastro. If you still face this problem, reopen the issue or comment on it. We would be happy to help again!

[astrojuanlu]

We haven't yet determined if this is a bug or not, naughty bot.