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 442 evaluations per sample.
Range (min … max): 232.885 ns … 38.264 μs ┊ GC (min … max): 0.00% … 99.13%
Time (median): 235.514 ns ┊ GC (median): 0.00%
Time (mean ± σ): 254.744 ns ± 493.181 ns ┊ GC (mean ± σ): 3.97% ± 4.16%
▄█▇▃ ▂▂▂▁ ▁▂ ▁ ▂▁ ▁
████▇▅▇▇▇▆▄▅▃▃▇██████▇▇▇▇▇▆▅▆▄▄▃▃▄▃▃▆▆█████▇██▇██▆▄▄▄▅▆▅▅▆▆▆▆ █
233 ns Histogram: log(frequency) by time 316 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.Goat())BenchmarkTools.Trial: 10000 samples with 133 evaluations per sample.
Range (min … max): 714.158 ns … 186.846 μs ┊ GC (min … max): 0.00% … 99.39%
Time (median): 769.353 ns ┊ GC (median): 0.00%
Time (mean ± σ): 784.730 ns ± 2.000 μs ┊ GC (mean ± σ): 3.30% ± 1.40%
▃██▄▁ ▆▇▆▄▃▅▄▄▂▁▃▄▃▁ ▁▁▁ ▁▂▂ ▁▁ ▂
██████▇▅▅▆▆▇▆▆▅▃▃▁▄▇██████████████████▇▆▅▅▆▆█████████▇█▇▇▆▅▆▇ █
714 ns Histogram: log(frequency) by time 875 ns <
Memory estimate: 128 bytes, allocs estimate: 1.using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Newton()))BenchmarkTools.Trial: 10000 samples with 352 evaluations per sample.
Range (min … max): 254.389 ns … 80.854 μs ┊ GC (min … max): 0.00% … 99.46%
Time (median): 262.128 ns ┊ GC (median): 0.00%
Time (mean ± σ): 303.309 ns ± 946.813 ns ┊ GC (mean ± σ): 5.58% ± 2.33%
▂▆██▆▂ ▁▂▃▂▂▁ ▂▅▆▅▆▄▃▄▅▄▅▃▁ ▁▁▂▁▁ ▁▁ ▃
██████▇▆▇█▇▅▅▄▁▃▇███████▇▇▆▅▆▆▅▆▅███████████████▆▆███████████ █
254 ns Histogram: log(frequency) by time 350 ns <
Memory estimate: 128 bytes, allocs estimate: 1.using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Thukral3B()))BenchmarkTools.Trial: 10000 samples with 247 evaluations per sample.
Range (min … max): 305.725 ns … 135.279 μs ┊ GC (min … max): 0.00% … 99.62%
Time (median): 337.105 ns ┊ GC (median): 0.00%
Time (mean ± σ): 362.294 ns ± 1.512 μs ┊ GC (mean ± σ): 6.36% ± 1.72%
▇█▇▄ ▁▂▂▁▁▁▂▄▃▄▅▅▆▅▄▂▂▃▃▄▄▃ ▁▁▁▁ ▂
█████▆▆▆▇▇▆▆▄▃▃▃▁▁▆███████████████████████▇▆▅▆▆▇██████▇█▇▇▇▇▇ █
306 ns Histogram: log(frequency) by time 418 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.A42()))BenchmarkTools.Trial: 10000 samples with 197 evaluations per sample.
Range (min … max): 468.269 ns … 169.754 μs ┊ GC (min … max): 0.00% … 99.53%
Time (median): 531.051 ns ┊ GC (median): 0.00%
Time (mean ± σ): 544.819 ns ± 1.805 μs ┊ GC (mean ± σ): 4.25% ± 1.40%
▄█▇▃ ▁▂▂▁▅▇▇▆▅▅▄▄▅▅▃▃▂▂ ▂▃▂▁▁▁▁▁▁ ▁ ▃
▅▅█████▆▇▇▇▆▄▄▄▁▁▃▆███████████████████▆▆████████████▆▇█▆▇▆▆▃▅ █
468 ns Histogram: log(frequency) by time 623 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Bisection()))BenchmarkTools.Trial: 10000 samples with 9 evaluations per sample.
Range (min … max): 2.226 μs … 5.979 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 2.242 μs ┊ GC (median): 0.00%
Time (mean ± σ): 2.330 μs ± 193.933 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▅ ▄▆▄▂▁▁ ▁
███▅▇▅▄▆▅▅▇███████▆▅▆▄▅▄▄▄▄▅▃▅▄▄▅▇█▆▆▄▃▃▄▁▁▃▁▁▁▁▁▁▆▆▇▆▆▇▆▆▅ █
2.23 μs Histogram: log(frequency) by time 3.28 μs <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.SuperHalley()))BenchmarkTools.Trial: 10000 samples with 369 evaluations per sample.
Range (min … max): 247.499 ns … 91.257 μs ┊ GC (min … max): 0.00% … 99.52%
Time (median): 303.737 ns ┊ GC (median): 0.00%
Time (mean ± σ): 310.742 ns ± 1.152 μs ┊ GC (mean ± σ): 7.12% ± 2.21%
▂▆█▇▄ ▁▂▂▂▁ ▁▄▅▆▆▆▆▆▅▅▅▅▅▄▁ ▁▂▁▂▂▂▁▂▂▂▁ ▃
██████▆▅▇▆▆▅▁▃▄███████▇▆▆▆▆▆▅▅▅███████████████▆█████████████ █
247 ns Histogram: log(frequency) by time 345 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Brent()))BenchmarkTools.Trial: 10000 samples with 208 evaluations per sample.
Range (min … max): 380.284 ns … 161.126 μs ┊ GC (min … max): 0.00% … 99.58%
Time (median): 436.786 ns ┊ GC (median): 0.00%
Time (mean ± σ): 449.925 ns ± 1.720 μs ┊ GC (mean ± σ): 4.92% ± 1.40%
█▃
▃▇██▄▂▂▂▂▂▂▂▂▂▂▂▁▁▁▂▂▂▃▂▃▇█▇▆███▅▄▄▃▃▄▄▄▃▂▂▂▂▂▂▂▂▃▂▂▃▂▂▂▂▂▂▂▂ ▃
380 ns Histogram: frequency by time 506 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Order2()))BenchmarkTools.Trial: 10000 samples with 239 evaluations per sample.
Range (min … max): 313.904 ns … 140.463 μs ┊ GC (min … max): 0.00% … 99.59%
Time (median): 369.385 ns ┊ GC (median): 0.00%
Time (mean ± σ): 380.070 ns ± 1.509 μs ┊ GC (mean ± σ): 5.15% ± 1.40%
▆█▆▁ ▁▂▁▁ ▃▆▆▇▆▆▆▇▅▄▄▄▄▄▄▄▂▁ ▁▂▁▁▂▂▁▁▁ ▃
████▇▅▆▄▇▇▅▄▅▁▁▅▅▅▄▇█████▆███████████████████▆███████████████ █
314 ns Histogram: log(frequency) by time 426 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.AlefeldPotraShi()))BenchmarkTools.Trial: 10000 samples with 197 evaluations per sample.
Range (min … max): 468.822 ns … 169.545 μs ┊ GC (min … max): 0.00% … 99.55%
Time (median): 525.812 ns ┊ GC (median): 0.00%
Time (mean ± σ): 539.752 ns ± 1.815 μs ┊ GC (mean ± σ): 4.35% ± 1.40%
█▂ ▅▃
▂███▃▂▂▂▂▂▂▂▂▂▂▂▁▁▁▂▂▂▃▃▂▅██▇▄▅▅▄▃▄▅▄▃▃▃▂▂▂▂▂▂▂▃▃▃▂▂▂▂▂▂▂▂▂▂▂ ▃
469 ns Histogram: frequency by time 599 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)24650.95447361446997022436542900109713321175234563222583018047585031821959181256Comparison
