Kepler Solvers
The heart of this package is being able to take a set of Keplerian elements and output relative positions, velocities, etc.
This normaly requires solving Kepler's equation numerically. This package supports a multitude of solver algorithms that can be passed to orbitsolve:
The last of these RootsMethod, allows one to substitute any algorithm from the Roots.jl package. These include many different classical and modern root finding algorithms.chosen precision, including artibrary precision BigFloats. Using big floats with, for example, Roots.PlanetOrbits.Thukral5B and a tight tolerenace, allows you to solve orbits up to arbitrary precision.
The default choice is Auto, which currently selects Markley for all cases. The Markley algorithm is very fast, reasonably accurate, and always converges, making it a good default choice.
The Markley algorithm is a tweaked version of the algorithm from AstroLib.jl. It is non-iterative and converges with less than 1e-15 relative error across the full range of e between 0 and 1. On my laptop, this solves for a single eccentric anomaly in just 71 ns. Since it is implemented in pure Julia, there is no overhead from calling into a C or Cython compiled function and no need for vectorization.
Examples
using PlanetOrbits, BenchmarkTools
orb = orbit(a=1.2, e=0.1, M=1.0, ω=1.4, τ=0.5)
t = mjd("2025-06-23")
@benchmark orbitsolve(orb, t, PlanetOrbits.Markley())BenchmarkTools.Trial: 10000 samples with 535 evaluations.
Range (min … max): 211.609 ns … 20.987 μs ┊ GC (min … max): 0.00% … 98.96%
Time (median): 215.225 ns ┊ GC (median): 0.00%
Time (mean ± σ): 251.662 ns ± 493.964 ns ┊ GC (mean ± σ): 5.46% ± 2.79%
█▆▃ ▁ ▁ ▂▂▅▅▅▆▅▄▃▂▂ ▂
█████▇▆▅▅▅▄▄▄▅▇██▇▇▇▇▄▄▄▄▃▄▁▃▁▁▅▃▃▃▁▁▁▆██████████████▇██████▇ █
212 ns Histogram: log(frequency) by time 284 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.Goat())BenchmarkTools.Trial: 10000 samples with 154 evaluations.
Range (min … max): 689.734 ns … 50.914 μs ┊ GC (min … max): 0.00% … 98.55%
Time (median): 741.321 ns ┊ GC (median): 0.00%
Time (mean ± σ): 738.287 ns ± 693.161 ns ┊ GC (mean ± σ): 1.33% ± 1.39%
█▅ ▃
▃███▄▃▂▂▂▂▂▂▂▂▂▁▂▁▁▂▃▄▃▂▃██▇▄▃▆█▆▄▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
690 ns Histogram: frequency by time 813 ns <
Memory estimate: 128 bytes, allocs estimate: 1.using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Newton()))BenchmarkTools.Trial: 10000 samples with 323 evaluations.
Range (min … max): 262.347 ns … 25.266 μs ┊ GC (min … max): 0.00% … 98.81%
Time (median): 309.495 ns ┊ GC (median): 0.00%
Time (mean ± σ): 305.894 ns ± 542.737 ns ┊ GC (mean ± σ): 3.96% ± 2.21%
█▄
██▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂▃▅▃▃█▇▄▅▆▄▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
262 ns Histogram: frequency by time 348 ns <
Memory estimate: 128 bytes, allocs estimate: 1.using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Thukral3B()))BenchmarkTools.Trial: 10000 samples with 261 evaluations.
Range (min … max): 303.364 ns … 30.483 μs ┊ GC (min … max): 0.00% … 98.90%
Time (median): 351.770 ns ┊ GC (median): 0.00%
Time (mean ± σ): 346.789 ns ± 519.333 ns ┊ GC (mean ± σ): 2.59% ± 1.71%
▄█▅▃ ▄▅▄▃▆▆▄▄▅▄▃▃▃▁▁▁ ▂
█████▅▆▇▇▆▆▄▄▄▅▅▄▁▃▃▁▁▃▆█▆▅▅▆██████████████████▇▇█▇▆▆▇▇▆▆▆▆▆▇ █
303 ns Histogram: log(frequency) by time 400 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.A42()))BenchmarkTools.Trial: 10000 samples with 110 evaluations.
Range (min … max): 770.627 ns … 73.771 μs ┊ GC (min … max): 0.00% … 98.88%
Time (median): 817.527 ns ┊ GC (median): 0.00%
Time (mean ± σ): 822.376 ns ± 730.320 ns ┊ GC (mean ± σ): 0.89% ± 0.99%
▅█▆▁▁ ▆▇▄▂▁▃▇▆▄▂▁▃▅▄▂▁▁▁▂▁ ▂
█████▆▄▄▅▆▆▄▄▂▃▄▇████████████████████▇▇▇▆▆▆▅▆▅▇▇▆▄▅▇▇▆▄▅▅▆▆▅▅ █
771 ns Histogram: log(frequency) by time 923 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Bisection()))BenchmarkTools.Trial: 10000 samples with 10 evaluations.
Range (min … max): 1.297 μs … 4.779 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.306 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.366 μs ± 122.867 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▇█▆▁ ▄▅▅▄▃▂▂ ▁▁ ▂
█████▇▇▇▆▆▃▄▄▇▇▆▄▁▃▃▁▁▁▁▁█████████▇▇▆▅▆▄▅▆▃▅▃▅▅▅▆▅▅▆▆████▇▇ █
1.3 μs Histogram: log(frequency) by time 1.65 μs <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.SuperHalley()))BenchmarkTools.Trial: 10000 samples with 313 evaluations.
Range (min … max): 263.783 ns … 25.655 μs ┊ GC (min … max): 0.00% … 98.79%
Time (median): 311.989 ns ┊ GC (median): 0.00%
Time (mean ± σ): 305.886 ns ± 504.814 ns ┊ GC (mean ± σ): 3.29% ± 1.98%
▇█▅▂ ▂▂▁▄▅▄▄▆▆▅▅▅▄▂▂▂▁ ▁ ▂
████▇▆█▇▅▄▆▅▄▅▁▃▃▁▄▅▆█▇▆▆▅▄▅▅███████████████████▇██▇▆▇▇▆▇▇▇▇▇ █
264 ns Histogram: log(frequency) by time 353 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Brent()))BenchmarkTools.Trial: 10000 samples with 8 evaluations.
Range (min … max): 3.886 μs … 7.497 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 3.936 μs ┊ GC (median): 0.00%
Time (mean ± σ): 4.029 μs ± 232.658 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▅▆█▅ ▁ ▁▃▅▅▃▁ ▁ ▂
██████▇█▅▃███████▇▇▅▅▄▇██▇▇▆▄▅▃▃▁▄▄▅▆▆▇▇▆▇████▇▅▃▄▅▅▆▆▆▅▇▇█ █
3.89 μs Histogram: log(frequency) by time 5.01 μs <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Order2()))BenchmarkTools.Trial: 10000 samples with 223 evaluations.
Range (min … max): 336.955 ns … 36.186 μs ┊ GC (min … max): 0.00% … 98.90%
Time (median): 385.117 ns ┊ GC (median): 0.00%
Time (mean ± σ): 376.690 ns ± 504.585 ns ┊ GC (mean ± σ): 1.89% ± 1.40%
▆█▅▂ ▅▆▅▃▄▆▆▄▃▄▄▃▂▁ ▂
█████▅▆▇▇▆▅▄▄▄▃▃▄▃▁▃▁▁▄▄▄▁▇███████████████████▆▇▇▆▅▇▆▆▅▆▇▆▆▆▇ █
337 ns Histogram: log(frequency) by time 435 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.AlefeldPotraShi()))BenchmarkTools.Trial: 10000 samples with 146 evaluations.
Range (min … max): 697.055 ns … 56.776 μs ┊ GC (min … max): 0.00% … 98.71%
Time (median): 744.541 ns ┊ GC (median): 0.00%
Time (mean ± σ): 747.190 ns ± 778.221 ns ┊ GC (mean ± σ): 1.47% ± 1.40%
▅█▅▂▂ ▂▆▅▃▁▄▇▆▄▃▃▅▅▄▂▁▂▂▁▁ ▂
█████▆▃▅▇▇▆▇▅▅▅▁▅▃▁████████████████████▇█▇▆▅▆▇▇▆▅▅▆▆▇▇▇▇▇▇▇▇▆ █
697 ns Histogram: log(frequency) by time 825 ns <
Memory estimate: 128 bytes, allocs estimate: 1.High precision
You can solve Kepler's equation in high precision using big floats and tightening the tolerance on the solver.
orb_big = orbit(a=big(1.2), e=big(0.1), M=big(1.0), ω=big(1.4), τ=big(0.5))
sol = orbitsolve(orb_big, big(t), PlanetOrbits.RootsMethod(Roots.Thukral5B(),rtol=1e-30,atol=1e-30,))
radvel(sol)25030.34323194365565831507041947847102607642627168092403020601399008600937133855Comparison
