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 550 evaluations.
Range (min … max): 209.211 ns … 20.581 μs ┊ GC (min … max): 0.00% … 98.85%
Time (median): 215.276 ns ┊ GC (median): 0.00%
Time (mean ± σ): 248.313 ns ± 454.402 ns ┊ GC (mean ± σ): 5.08% ± 2.78%
██▃ ▁ ▃▅▅▆▆▆▅▄▃▂▁ ▂
███▇██▆▆▇▅▄▄▄▄▇██▇▇▆▄▄▃▄▅▄▅▅▅▄▁▃▃▃▃▃▁▁▃▃████████████▇▇██▇▇███ █
209 ns Histogram: log(frequency) by time 281 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.Goat())BenchmarkTools.Trial: 10000 samples with 159 evaluations.
Range (min … max): 662.811 ns … 50.226 μs ┊ GC (min … max): 0.00% … 98.29%
Time (median): 726.836 ns ┊ GC (median): 0.00%
Time (mean ± σ): 720.678 ns ± 690.826 ns ┊ GC (mean ± σ): 1.35% ± 1.39%
█ ▁▅ ▇▂ ▂
▄█▃▃▂▂▂██▃▃▂▂▂▂▂▂▂▂▁▂▂▁▂▂▂▂▃██▅▄██▅▄▃▄▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
663 ns Histogram: frequency by time 800 ns <
Memory estimate: 128 bytes, allocs estimate: 1.using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Newton()))BenchmarkTools.Trial: 10000 samples with 334 evaluations.
Range (min … max): 263.159 ns … 28.713 μs ┊ GC (min … max): 0.00% … 98.54%
Time (median): 309.189 ns ┊ GC (median): 0.00%
Time (mean ± σ): 306.051 ns ± 534.370 ns ┊ GC (mean ± σ): 3.88% ± 2.21%
█▆▁ ▃▃▂▅▅▄▆▅▄▄▃▂▁ ▂
███▇▆█▄▅▃▄▁▁▃▄▄▃▇█▇▆▅▆▆▃▃▁█████████████▆█▇▇▇▇▆█▇▇▇▇▆▆▄▆▆▅▄▄▃▄ █
263 ns Histogram: log(frequency) by time 365 ns <
Memory estimate: 128 bytes, allocs estimate: 1.using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Thukral3B()))BenchmarkTools.Trial: 10000 samples with 251 evaluations.
Range (min … max): 296.131 ns … 31.054 μs ┊ GC (min … max): 0.00% … 99.01%
Time (median): 345.709 ns ┊ GC (median): 0.00%
Time (mean ± σ): 339.626 ns ± 529.790 ns ┊ GC (mean ± σ): 2.70% ± 1.71%
▁█▇▃▁ ▂▂ ▁▆▆▄▄▆▆▄▄▄▄▂▁▁ ▂
█████▆▅▇▇▅▁▁▄▁▃▁▁▅▁▄▁▁▁▁▇████████████████████▆▆▆▆▆▅▇▇▇▆▇▆▇▅▆▇ █
296 ns Histogram: log(frequency) by time 392 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.A42()))BenchmarkTools.Trial: 10000 samples with 116 evaluations.
Range (min … max): 758.060 ns … 1.173 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 812.897 ns ┊ GC (median): 0.00%
Time (mean ± σ): 803.607 ns ± 34.451 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▅█▄▂▁ ▁▆▆▃▁▂▇▇▄▃▃▆▆▃▂▂▂▂ ▂
█████▅▆▆▆▇▆▃▃▁███████████████████▇▇▇▆▆▇▅▆▆▅▆▇█▆▇▆▇█▆▇▇▆▇▇▆▇▇ █
758 ns Histogram: log(frequency) by time 919 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 … 3.846 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.304 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.366 μs ± 123.671 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
██ ▂▅▅▄▃▂▁ ▁▁ ▂
███▇▆▆▆▆▆▅▅▄▃▇▅▁▁▁▁▁▁▁▁▁▄████████▆▆▆▅▄▅▅▆▅▄▄▄▄▄▄▃▅▆▇███▇▇▆▆ █
1.3 μs Histogram: log(frequency) by time 1.66 μs <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.SuperHalley()))BenchmarkTools.Trial: 10000 samples with 328 evaluations.
Range (min … max): 263.360 ns … 24.337 μs ┊ GC (min … max): 0.00% … 98.74%
Time (median): 311.649 ns ┊ GC (median): 0.00%
Time (mean ± σ): 307.826 ns ± 531.856 ns ┊ GC (mean ± σ): 3.86% ± 2.21%
▆█▅▂ ▄▄▂▅▆▅▅▆▅▄▄▃▂▁▁ ▂
████▇▇█▆▄▄▄▁▃▁▄▃▁▁▁▅█▇▅▆▃▃▄▃▄▃▆███████████████▇▇▇▇▇▇▇▇▆▆▇▇▆▇▇ █
263 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.832 μs … 6.243 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 3.865 μs ┊ GC (median): 0.00%
Time (mean ± σ): 3.938 μs ± 162.489 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▅█▆▄▄▄▃▂▂▁ ▂▅▅▃▃▂▂▁ ▁ ▂
██████████▆▆▆██████████▇▆▅▇████▇▇▆▆▄▄▃▄▅▁▃▃▃▄▄▃▃▃▄▄▅▁▁▄▁▁▁▄ █
3.83 μs Histogram: log(frequency) by time 4.73 μs <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Order2()))BenchmarkTools.Trial: 10000 samples with 221 evaluations.
Range (min … max): 332.072 ns … 36.422 μs ┊ GC (min … max): 0.00% … 98.94%
Time (median): 380.443 ns ┊ GC (median): 0.00%
Time (mean ± σ): 375.588 ns ± 619.999 ns ┊ GC (mean ± σ): 2.85% ± 1.71%
██
██▃▃▂▂▂▂▂▂▂▂▁▁▁▁▂▁▁▁▂▁▂▂▂▂▂▂█▇▄▃▇█▄▃▄▄▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
332 ns Histogram: 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 147 evaluations.
Range (min … max): 697.565 ns … 54.075 μs ┊ GC (min … max): 0.00% … 98.64%
Time (median): 746.231 ns ┊ GC (median): 0.00%
Time (mean ± σ): 742.369 ns ± 534.281 ns ┊ GC (mean ± σ): 0.72% ± 0.99%
▆█▆▃▁ ▄▆▄▁▁▅▇▆▄▂▄▆▅▃▂▂▂▂▁ ▂
█████▅▆▆▇▇▅▅▃▁▃▃▃▁▁▁█████████████████████▇▇▆▅▆▇▆▆▆▇▇▇▆▅▇█▇▇▆▇ █
698 ns Histogram: log(frequency) by time 824 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
