Offer an alternative to Orbit.propagate for the Moon and planets

[astrojuanlu]

🐞 Problem

Currently Orbit.from_body_ephem is a convenient way to get the osculating orbit of a planet or the Moon for plotting, to retrieve some classical element, or other things. However, there are some inconsistencies:

• It gives the false impression that these objects can be accurately propagated or sampled, while the correct way to retrieve planetary trajectories would be directly querying the ephemerides or using the undocumented build_ephem_interpolant.
• For planets, it sets the attractor to the Sun, while in all seriousness it should be the Solar System Barycenter, to match the original data.

The places where from_body_ephem is used are:

1. Get approximate orbital elements from the planets to compute spheres of influence in poliastro.threebody.soi
2. Plot approximate (osculating) orbits of planets
3. Get precise position vectors of planets to compute Lambert arcs

🎯 Goal

We should provide an alternative that works equally well for all use cases but that prevents the user from shooting themself in the foot.

💡 Possible solutions

We can have a special "attractor", the SSB, that does not have a k value and therefore prevents the orbit from being propagated and sampled. However, this still does not fix the Moon case, unless we have another object which is EMB (Earth-Moon Barycenter).

📋 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 🎉

[astrojuanlu]

Idea:

We have Orbit objects, which represent an osculating orbit. Of course, they can be propagated with any propagator, including one that includes perturbations.

Perhaps we should have another type of object that follows a similar API than Orbit, but whose propagate method does a totally different thing and, for example, looks for external ephemerides.

[Sedictious]

@Juanlu001 I'd be interested in tackling the problem, but I fear I'm a bit lost. If I understand correctly, the main problem is that while the attractor is defined as the sun (and the moon respectively). This isn't completely the case as the motion is affected by other bodies, so we want to define the orbit in terms of the barycenter of the system. So I guess the solution would be to translate the position into the actual coordinate reference system (rather than GCRS or ICRS). Is my intuition correct?

[astrojuanlu]

Hi @Sedictious, sorry for the delay!

You more or less got it. The problem in general is that Orbit.propagate works great for spacecraft, satellites and anything small that fulfills the assumptions of the two body problem plus having a negligible mass with respect to the attractor. However, for planetary bodies, things are more complicated and we don't want to Orbit.propagate the Earth or Moon, but to use the ephemerides data that is readily available. Using Orbit objects for both things has worked well so far but I think it's somehow broken.

[Sedictious]

@Juanlu001 Thanks for the reply! The only part that still eludes me, is when you refer to the ephemerides. Like you mentioned in your OP, the correct way to retrieve the ephemerides would be to retrieve the data from an external source:

1) Is there an existing function (or an interest in one) for retrieving and parsing ephemerides data from an external, reputable source?

2) Wouldn't it help to recalculate the attractor of the system if the mass of the orbiting body isn't negligible (like in the case of satellites and spacecrafts you mentioned)? I may be neglecting the nuances of the problem, but what keeps us from introducing another function which recalculates the barycentric coordinates depending on the mass of the planetary objects?

[astrojuanlu]

Is there an existing function (or an interest in one) for retrieving and parsing ephemerides data from an external, reputable source?

Yes: Orbit.from_body_ephem uses Astropy and jplephem internally, which retrieves the data directly from JPL, and recently we implemented Orbit.from_horizons and Orbit.from_sbdb as well, which also use reputable sources.

Regarding your second question, I didn't get it entirely so I will try to rephrase: there is a mismatch when saying that Orbit.from_body_ephem(Earth) orbits the Sun and that its reference frame is ICRS. This gets worse when propagating, because that kind of implies that the barycenter of the Solar System matches the center of the Sun, which it doesn't. If we want to solve this problem we either have to

a) Keep ICRS, but change the attractor to something like SSB (Solar System Barycenter), or b) Forget about .propagate() ing for bodies altogether.

Let's think of a extreme case: the Sun. The first option would allow us for things like Orbit.from_body_ephem(Sun).propagate(), which would give totally wrong results, because the Sun wobbles around the SSB in a totally non-keplerian orbit:

(Source)

For a spacecraft in the Solar System, we can choose to study its theoretical orbit around the Sun, assuming no perturbations from other bodies (trivial with poliastro) or actually try to model these perturbations as well (not trivial, but easy thanks to poliastro.core.perturbations.third_body).

However, even if we try to "fix" this issue, I argue that giving the user the option of naïvely Orbit.propagate()-ing a body, when we actually know with very high precision its past and future positions, it's too misleading and can cause confusion. This last point is the one that is open to discussion.

[Sedictious]

@Juanlu001 Ah thank you, I think I now get the issue. So propagating non-Keplerian orbits gives us an error margin so large, its rendered more or less useless or even directly misleading.

While offering the naive version of propagate() might prove to be a source of confusion, looking back at the docs, I think the "preferred" way is thoroughly explained. Maybe, it would help to expand the docs to explain your reasoning and why propagate() should be avoided when dealing with planets.

[astrojuanlu]

Maybe, it would help to expand the docs to explain your reasoning and why propagate() should be avoided when dealing with planets.

True, but that would be the "easy" option for the developer, and perhaps more confusing to the user :) If we provide an API, but there are "red flags" they have to be aware of, then it's not so convenient. Perhaps it's not that problematic, but at the moment I don't have the time to pause and reflect about this in detail.

[jorgepiloto]

By creating the attractor SSB (as said in possible solutions' section) and automatically assigning it when using from_body_ephem(), it would be possible to use it as a flag when calling the method propagate. Something like:

# Inside propagate function

if orbit.attractor == SSB:
    # We propagate by build_ephem_interpolant
else:
    # We are under the two-body assumption

Do you think this implementation is a good approach to the problem?

[Sedictious]

@jorgepiloto Couldn't we define a more generic system barycenter which could default to SSB, but allow for alternatives (e.g. Earth-Moon barycenter)

[jorgepiloto]

We could have something like _Barycenter class and follow the same idea _Body does. The problem comes when defining k = 0. * u.m ** 3 / u.s ** 2. Since many internal methods make use of this value, different errors arise. How could we deal with this? What about adding all solar system bodies masses for the SSB case? Is it a 'brute force' solution? :rofl:

[astrojuanlu]

Every time I get back to this I have to carefully think about it again :innocent: My favourite explanation on barycenters and orbits is this Space.SE answer.

First, from a physical point of view, if you experiment a bit with things like https://hermann.is/gravity/, you will notice that a small object very close to one attractor in a binary system does not orbit the barycenter of the system, but its main attractor.

Taking into account that the Solar System is much more complex than that, I'd say that within the usual assumptions that hold in poliastro, "orbiting around the SSB" does not make much sense. Remember that poliastro is mostly concerned with spacecraft (i.e. the mass of the object is negligible with respect to the mass of its attractor) and the two-body problem (extending it using the patched conics framework and also allowing some applications of the restricted three-body problem).

To recap, the only interface that gives us ICRS coordinates (therefore, SSB-centered) is from_horizons and from_body_ephem, which use JPL HORIZONS in both cases. Regarding what @Sedictious said:

Couldn't we define a more generic system barycenter which could default to SSB, but allow for alternatives (e.g. Earth-Moon barycenter)

If we had an object far enough from the Earth that could be considered to orbit the EMB (Earth-Moon barycenter), it would, in fact, be in a heliocentric orbit. Therefore, we probably don't need to define other barycenters.

The places where from_body_ephem is used are:

1. Get approximate orbital elements from the planets to compute spheres of influence in poliastro.threebody.soi.
2. Plot approximate (osculating) orbits of planets
3. Get precise position vectors of planets to compute Lambert arcs

I still think we should serve these two use cases (approximate, osculating orbits and precise position vectors for targetting) with two different APIs and without giving an Orbit object (which has methods that don't make any sense for planets, like apply_maneuver! apart from all the subtleties listed before).

It took me a lot of time to write the comment above and discarded several ideas along the way, so I will let this sink and postpone it to the next release.

[astrojuanlu]

One other thing I forgot after rereading the above:

I think from_horizons is a very good example of the dangers of putting real-world data in an Orbit object: think of comets whose orbits get deflected by nearby passes around the planets, or even deep space exploration spacecraft that perform trajectory correction maneuvers and flybys! If we naïvely .propagate an object coming from Orbit.from_horizons, we automatically neglect all these effects. For me, it's another case in favor of being careful of what's an Orbit and what's not.

[astrojuanlu]

sbpy solved this by having an Orbit object and an Ephem object:

https://sbpy.readthedocs.io/en/latest/sbpy/data.html#what-are-ephem-orbit-and-phys-and-how-to-use-them

(this came up when I reviewed their submission to the Journal of Open Source Software https://github.com/openjournals/joss-reviews/issues/1426)

[jorgepiloto]

Interesting, I think that Ephem class could be a good option. We may have something like:

ss = Ephem.from_body(args)
ss_f = ss.propagate(1 * u.year)

The principle behind this propagate would be the one defined at the very beginning of the issue: use the build_ephem_interpolant or get_body_barycentric_posvel.

Should this Ephem be inherited from the Orbit class and the propagate method overridden?

[astrojuanlu]

I'm not sure Ephem.propagate would make sense (just for consistency, but it would be a bit of a stretch in my view). Ephem.sample, however, would make total sense.

Let's talk about this more in detail for the next release. Before talking about inheritance, perhaps we would need to tackle https://github.com/poliastro/poliastro/issues/659 first. Or, in general, the fact that the Orbit class is growing quite a lot.

[hpdumm]

Orekit uses a CelestialBody object and provides a method to retrieve the current position and velocity. That would be less generic than Ephem.

[astrojuanlu]

Thanks @hpdumm for the pointer! This is the relevant documentation:

https://www.orekit.org/site-orekit-9.3.1/apidocs/org/orekit/bodies/CelestialBody.html

And indeed it provides a getPVCoordinates method that receives a date and a frame, and not a propagate method.

In my view, an Ephem object would be a vectorized version of that (i.e. allowing an array of times for simplicity), but the idea is the same.

[astrojuanlu]

I will tackle this and #790 together.

[astrojuanlu]

These are my notes: https://github.com/poliastro/poliastro/wiki/Body-orbits

[astrojuanlu]

Roadmap:

1. [x] Make plane part of BaseState objects because it's where it belongs: we shouldn't have elements sets without a reference.
2. [x] Use ClassicalState to provide a static, referenced, representative (mean) orbit for bodies (see this comment) hence fix #790.
3. [x] Address #824.
4. [ ] Introduce Ephem objects holding time histories and states that should not be propagated (see this comment)
5. [ ] Replace all internal uses of Orbit.from_body_ephem with the mean Body orbit.
6. [ ] Write documentation that explains the transition and emphasize the Orbit.from_body_ephem deprecation.
7. [ ] Release
8. [ ] Remove Orbit.from_body_ephem

[astrojuanlu]

Proposal: https://github.com/poliastro/poliastro/wiki/Ephem-objects

[astrojuanlu]

A new Ephem object has been introduced that fixes this dichotomy. I will update the User Guide soon, stay tuned! In the meanwhile, take a look at the code and especially the notebooks, which are updated already.