Implement unbound parabolic and hyperbolic orbits

[astrojuanlu]

I'm trying the library for the first time and noticed a couple of strange things. It's my very first time running Julia, so bear with me.

I'm using the following semimajor axis and attractor mass (which I understand it's expressed in solar masses):

a = 10000e3
M = 3e-6

And the following code for the orbit definition and plotting:

using PlanetOrbits
using Plots

elems = KeplerianElementsDeg(
       a = a,
       M = M,
       i = 0,
       e = e,
       ω = 0,
       Ω = 0,
       plx = 100,
       τ = 0
       )

plot(elems)

(what's the parallax for, by the way? I typed 100 because I didn't know what to write but also it didn't seem to affect the results)

For near-parabolic orbits e = 0.99999, I noticed that the representation rounding is a bit too agressive, appearing as 1.0. In any case, the correct value is stored:

julia> elems = KeplerianElementsDeg(
       a = a,
       M = M,
       i = 0,
       e = 0.99999,
       ω = 0,
       Ω = 0,
       plx = 100,
       τ = 0
       )
KeplerianElements{Float64}
─────────────────────────
a   [au ] = 1.0e7
e         = 1.0
i   [°  ] = 0.0
ω   [°  ] = 0.0
Ω   [°  ] = 0.0
τ         = 0.0
M   [M⊙ ] = 3.0e-6 
plx [mas] = 100.0 
──────────────────────────
period      [yrs ] : 1.82468092739242e13 
distance    [pc  ] : 10.0 
mean motion [°/yr] : 1.97e-11 
──────────────────────────


julia> elems.e
0.99999

In this case, a straight line is produced as a plot, probably because the x range is not properly computed:

In the hyperbolic case, the eccentricity is set to 1.0 silently:

julia> elems = KeplerianElementsDeg(
       a = a,
       M = M,
       i = 0,
       e = 1.2,
       ω = 0,
       Ω = 0,
       plx = 100,
       τ = 0
       )
KeplerianElements{Float64}
─────────────────────────
a   [au ] = 1.0e7
e         = 1.0
i   [°  ] = 0.0
ω   [°  ] = 0.0
Ω   [°  ] = 0.0
τ         = 0.0
M   [M⊙ ] = 3.0e-6 
plx [mas] = 100.0 
──────────────────────────
period      [yrs ] : 1.82468092739242e13 
distance    [pc  ] : 10.0 
mean motion [°/yr] : 1.97e-11 
──────────────────────────


julia> elems.e
1.0

And a seemingly empty plot is produced:

[astrojuanlu]

Just noticed that I was setting an absurd SMA since the units are au and not m, but it's only a scale effect. Double checked that with a = 10000e3 / 150e9 the results are the same.

[sefffal]

Hi @astrojuanlu thanks for the bug report and for trying the package out (and Julia)!

There's a few things going on.

First, the parameter plx is the parallax distance to the system in milliarcseconds. The main type KeplerianElements represents the visual orbit of a secondary in the plane of the sky. We should probably rename this and create another type to represent a general orbit in 3D space.

Next, the plot recipe for the orbital elements has aspectratio=1 by default, but we can turn that off

plot(elems, aspectratio=:auto)

Next, I see that the heuristic for distributing points along the orbit is not working well for this case:

scatter(elems, aspectratio=:auto, title="e=0.995")

I used scatter instead of plot to see the points instead of a line. Right now the plot recipe solves the orbit at a set of points that are equidistant in true anomaly.

I think a good solution would be to plot points that are equidistant in edit:~mean~eccentric anomaly instead:

Here's solutions for a wide range of eccentricities:

elems = [
    KeplerianElementsDeg(
              a = 10000e3/150e9,
              M = 3e-6,
              i = 65,
              e = e,
              ω = 0,
              Ω = 0,
              plx = 100,
              τ = 0
      )
     for e in 1 .- 10.0 .^(-9:0)
]
plot(elems, aspectratio=:auto)

I'll go ahead and tag a new version with that fix.

Finally, hyperbolic orbits are not yet supported by the kepler solver we are using, but I would welcome either a pull request implementing it or a pointer towards a good algorithm for hyperbolic cases.

[sefffal]

Version 0.1.1 has the plotting fix. Once merged in the registry you can update by entering package mode and running update (type ] update).

[astrojuanlu]

I think a good solution would be to plot points that are equidistant in mean anomaly instead:

You might want to check https://apps.dtic.mil/dtic/tr/fulltext/u2/a605040.pdf out:

In poliastro we are doing uniform sampling of the eccentric anomaly for closed orbits:

https://github.com/poliastro/poliastro/blob/321ca443f47e64eac40b462c31ca72a80cefa3b5/src/poliastro/twobody/sampling.py#L8-L41

[sefffal]

Ok, PlanetOrbits.jl now does the same uniform sampling of eccentric anomaly. Thanks for your help @astrojuanlu!

[astrojuanlu]

Thanks for the swift response by the way! Answering one last thing:

Finally, hyperbolic orbits are not yet supported by the kepler solver we are using, but I would welcome either a pull request implementing it or a pointer towards a good algorithm for hyperbolic cases.

I see you implemented Markley's algorithm. In poliastro we are using a method by Farnocchia et al http://dx.doi.org/10.1007/s10569-013-9476-9 that has the property of also converging in the near-parabolic region (both elliptic and hyperbolic), you can read the source code here:

https://github.com/poliastro/poliastro/blob/0.16.x/src/poliastro/core/propagation/farnocchia.py

It's MIT licensed, so the only thing we ask is that the attribution is retained. And you can also cite our upcoming paper :)

[astrojuanlu]

(Our own implementation of Markley is to be found here https://github.com/poliastro/poliastro/blob/3c504db/src/poliastro/core/propagation/markley.py)

[sefffal]

Okay nice, I'll take a look. One thing I really liked about the Markley algorithm is that it's non iterative. This lets it work on the GPU and works nicely with some kinds of auto-diff backends like Enzyme.jl. Which means we can write a model or cost function that computes orbits and automatically get gradients to accelerate optimization etc.

I think the best solution would be to support multiple algorithms by passing them through to orbitsolve, e.g. orbitsolve(elem, Markley()) or orbitsolve(elem, Farnocchia()).

[sefffal]

I addded support for choosing different kepler solver algorithms including options with arbitrary precision: https://sefffal.github.io/PlanetOrbits.jl/dev/kepler/

Everything is now in place to allow solving unbound orbits, except that some assumptions in the code will have to be relaxed. For example, we pre-calculate some factors like √((1 + e)/(1 - e)) which are no longer valid for e > 1.

[sefffal]

@astrojuanlu I'm working to add support for hyperbolic orbits here.

Do you happen to have a reference for determining a hyperbolic orbit from an (x,y,z,vx,vy,vz) point?

I haven't been able to find anything in the literature.

[astrojuanlu]

I don't, and I'm not actively working on poliastro or Astrodynamics anymore 😢

[sefffal]

Oh, sorry for the notification spam then!

[astrojuanlu]

No worries! Best luck ✨