Body orbits - poliastro/poliastro GitHub Wiki

Objectives

  1. poliastro should prevent users from propagating the osculating orbit of a planet
  2. There should be an easy way to retrieve mean/representative orbital elements of planets

Current problems

  • Orbit.heliosynchronous has this hardcoded:
n_sunsync = 1.991063853e-7 * u.one / u.s  # 2 * np.pi * u.rad / (86400 * 365.2421897 * u.s)

therefore, it's only valid for the Earth.

  • Orbit.geostationary asks for the period (or the angular velocity) of the body:
def geostationary(cls, attractor, angular_velocity=None, period=None, hill_radius=None):
    ...
    geo_radius = np.cbrt(attractor.k / np.square(angular_velocity.to(1 / u.s)))
  • threebody.soi does this:
if a is None:
    # https://github.com/poliastro/poliastro/pull/679#issuecomment-503902597
    from poliastro.twobody.orbit import Orbit

    a = Orbit.from_body_ephem(body, J2000_TDB).a

which internally does this:

r, v = get_body_barycentric_posvel(body.name, J2000_TDB)
...
# TODO: The attractor is not really the Sun, but the Solar System Barycenter
ss = cls.from_vectors(Sun, r.xyz.to(u.km), v.xyz.to(u.km / u.day), J2000_TDB)
ss._frame = ICRS()  # Hack!

Use cases

When I want "the orbit of a body", I might mean several things:

  • The osculating orbit in a particular instant (for targetting):
# Get position and velocities for departure and arrival
rr_dpt_body, vv_dpt_body = _get_state(departure_body, t_launch)
rr_arr_body, vv_arr_body = _get_state(target_body, t_arrival)

# Transform into Orbit objects
attractor = departure_body.parent
ss_dpt = Orbit.from_vectors(attractor, rr_dpt_body, vv_dpt_body, epoch=t_launch)
ss_arr = Orbit.from_vectors(attractor, rr_arr_body, vv_arr_body, epoch=t_arrival)

# Define time of flight
tof = ss_arr.epoch - ss_dpt.epoch

...
man_lambert = Maneuver.lambert(ss_dpt, ss_arr)

Although in reality, Maneuver.lambert retrieves position and velocity again, so do I need a "true osculating" orbit in this case? It's an artifact of the Maneuver API and therefore Maneuver.lambert.

  • Mean/representative orbital elements for design calculations.

Apart from the problems above... I also need the period for plotting! Otherwise, how do I know when to stop? Retrieving the true ephemerides won't give a closed orbit, and without the period I also don't know what time span should I use.

  • Remember, "what's the period of a planet?" At some point at pyastro I went ahead and looked how the "365 days" was measured.
  • Related: does the Earth mean orbit from JPL have this 365 days period?
  • A closed orbit appealing for plots:
def plot_solar_system(...):
    ...
    for body in bodies:
        orb = Orbit.from_body_ephem(body, epoch)
        orbit_plotter.plot(orb, label=str(body))

Why is a closed orbit appealing? Because if I retrieve true positions there is no such thing as a "complete revolution": there will either be a gap or some overlap. Also, two planet orbits from different epochs will be different but close, giving an ugly impression on plots and making them very difficult to distinguish.

  • An ephemerides (i.e. not an orbit, but a trajectory) (for what? well, for everything else!)

Possible solutions

Objective 1:

I still think we should remove Orbit.from_body_ephem. The catch is that we should cover all its use cases to prevent users from asking for it again (although it should definitely go through a deprecation cycle before we remove it).

Ideas:

  • Forbidding .propagate if frame is ICRS
    • Setting the frame to ICRS is a hack in itself and we should not rely on it
  • Add an Ephemerides/Trajectory object
    • Does't really replace .from_body_ephem, unless people use it with a single time value
    • Perhaps just a wrapper for get_body_barycentric_posvel
  • Create a new type of orbit without propagate

Objective 2:

We definitely need to add the mean elements to Body, and all the "problems" above would be solved.

Ideas:

  • Add raw elements only, and methods or functions on top of that
    • ...
  • Add an Orbit around the parent
    • Circular dependencies in two_body_ephem (about to disappear), default of heliosynchronous (could go), some constructors that assume the Solar System (from_horizons, from_sbdb), and to_icrs (which should probably disappear too!)
    • Also, probably too much?

If the Orbit has an attractor (a Body) and the attractor can itself have an Orbit around another attractor... Either this is a smell, or an excellent opportunity for composability.

The attractor of an Orbit can either have a parent (therefore, an orbit) or not (for example, an orbit around the Sun).

The end result would be:

n_sunsync = body.orbit.n
a = body.orbit.a

for body in bodies:
    orbit_plotter.plot(body.orbit, label=str(body))

Retort: what should body.orbit.propagate do? Perhaps I do need a new type of orbit.

Other notes

  • Do I even need to propagate when doing .sample()? It is just changing the anomaly. Also, retrieving time from true anomaly does not require solving for E in M = E - ecc \sin{E}, but the other way around (there's a closed-form formula from true to eccentric).

    • The problem is that the implementation would be even probably uglier because we don't have KeplerianRepresentation.

  • Should I plot the spheres of influence of the planets?

  • About the orbit.frame: it is only used in represent_as (unnecessary!), to_icrs (should go, see above), __str__ if available, and propagate if available.

    • I see this lineage: OsculatingOrbit (no propagate, no frame), Orbit (no frame), SolarSystemOrbit (the current Orbit, with frame and also maybe constructors depending on the Sun), and later EarthOrbit!