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 310 evaluations per sample.
Range (min … max): 270.342 ns … 58.133 μs ┊ GC (min … max): 0.00% … 99.18%
Time (median): 274.445 ns ┊ GC (median): 0.00%
Time (mean ± σ): 303.065 ns ± 758.717 ns ┊ GC (mean ± σ): 4.78% ± 3.23%
▃██▅▁ ▂▂▂▁ ▁▁▁▃▃▃▂▁ ▁▂▁▁ ▁
█████▇▇█▇▅▄▂▄▂▃▃▅██████▇▇▆▆▄▄▄▄▆██████████████▆▅▄▄▅▅▆▆▆▆▅▅▅▅▅ █
270 ns Histogram: log(frequency) by time 378 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.Goat())BenchmarkTools.Trial: 10000 samples with 73 evaluations per sample.
Range (min … max): 832.630 ns … 200.124 μs ┊ GC (min … max): 0.00% … 99.16%
Time (median): 894.630 ns ┊ GC (median): 0.00%
Time (mean ± σ): 931.789 ns ± 1.997 μs ┊ GC (mean ± σ): 2.13% ± 0.99%
█▂ ▂▇▁▅▃▂▃▂▂ ▂▁▂▁ ▁
██▅▆█████████▇▇██████▇▆▆▇▇▇▅▅▆▅▆▅▄▅▄▅▅▄▅▅▅▆▅▅▃▄▅▅▅▄▅▅▃▄▄▃▄▄▃▄ █
833 ns Histogram: log(frequency) by time 1.56 μs <
Memory estimate: 128 bytes, allocs estimate: 1.using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Newton()))BenchmarkTools.Trial: 10000 samples with 255 evaluations per sample.
Range (min … max): 294.992 ns … 92.685 μs ┊ GC (min … max): 0.00% … 99.41%
Time (median): 300.020 ns ┊ GC (median): 0.00%
Time (mean ± σ): 339.610 ns ± 1.131 μs ┊ GC (mean ± σ): 5.61% ± 1.96%
▂██▅▁ ▂▂▂▁▁ ▄▄▅▄▂ ▁▂▃▃▂ ▂
█████▇▇█▇▆▄▅▁▁▁▁▁▃▇██████▇▆▅▅▅█████████████▇▇▆▅▅▅▆██▇█▇▆▅▅▇▇ █
295 ns Histogram: log(frequency) by time 416 ns <
Memory estimate: 128 bytes, allocs estimate: 1.using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Thukral3B()))BenchmarkTools.Trial: 10000 samples with 211 evaluations per sample.
Range (min … max): 354.090 ns … 154.839 μs ┊ GC (min … max): 0.00% … 99.55%
Time (median): 417.360 ns ┊ GC (median): 0.00%
Time (mean ± σ): 433.363 ns ± 1.782 μs ┊ GC (mean ± σ): 6.42% ± 1.72%
█▇ ▁▁▁
▃██▅▂▂▂▂▂▂▂▂▂▂▁▂▁▁▂▃▃▂▂▂▆▇▆▄███▅▃▃▃▃▄▄▄▂▂▂▂▂▂▂▂▃▃▂▂▂▂▂▂▂▂▂▂▂▂ ▃
354 ns Histogram: frequency by time 501 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.A42()))BenchmarkTools.Trial: 10000 samples with 191 evaluations per sample.
Range (min … max): 525.293 ns … 200.663 μs ┊ GC (min … max): 0.00% … 99.56%
Time (median): 589.681 ns ┊ GC (median): 0.00%
Time (mean ± σ): 603.567 ns ± 2.143 μs ┊ GC (mean ± σ): 4.57% ± 1.40%
▄█▇▄ ▃▃▂ ▅█▇▅▃▁ ▄▅▄▃ ▁▃▃▂▁ ▁▁ ▂
█████▆▆▇█▇▁▁▄▄▁▁▁▁▁▁████████████████████▆▆▅▆▆▅▆██████▇▆▇███▇▇ █
525 ns Histogram: log(frequency) by time 669 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.542 μs … 13.372 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 2.556 μs ┊ GC (median): 0.00%
Time (mean ± σ): 2.589 μs ± 288.147 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▆▂ ▁
███▇▅▅▃▄▅▃▅▄▄▄▃▄▄▆▄▄▁▅▁▃▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅▆███ █
2.54 μs Histogram: log(frequency) by time 3.59 μs <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.SuperHalley()))BenchmarkTools.Trial: 10000 samples with 247 evaluations per sample.
Range (min … max): 303.939 ns … 155.940 μs ┊ GC (min … max): 0.00% … 99.65%
Time (median): 369.462 ns ┊ GC (median): 0.00%
Time (mean ± σ): 382.412 ns ± 1.779 μs ┊ GC (mean ± σ): 7.20% ± 1.72%
▂██▅▁ ▂▂▂ ▃▅▅▄▇▇▇▅▄▃▂▅▅▅▄▁ ▁ ▁▂▂▂▁ ▁▁▁ ▃
█████▆▇▇▇▅▄▄▁▅▃▃▁█████▇▆▆████████████████▇▆█████████▇█████▇▆▅ █
304 ns Histogram: log(frequency) by time 440 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Brent()))BenchmarkTools.Trial: 10000 samples with 200 evaluations per sample.
Range (min … max): 417.325 ns … 192.093 μs ┊ GC (min … max): 0.00% … 99.61%
Time (median): 479.575 ns ┊ GC (median): 0.00%
Time (mean ± σ): 495.571 ns ± 2.057 μs ┊ GC (mean ± σ): 5.36% ± 1.40%
▂█ ▄
▃██▇▃▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▃▂▂▃███▅▄▆▅▄▃▄▄▃▃▃▃▃▂▂▂▂▂▃▃▃▂▂▂▂▂▂▂▂▂▂▂ ▃
417 ns Histogram: frequency by time 556 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Order2()))BenchmarkTools.Trial: 10000 samples with 202 evaluations per sample.
Range (min … max): 387.757 ns … 192.368 μs ┊ GC (min … max): 0.00% … 99.61%
Time (median): 452.812 ns ┊ GC (median): 0.00%
Time (mean ± σ): 467.272 ns ± 2.056 μs ┊ GC (mean ± σ): 5.67% ± 1.40%
▆█▁ ▅▂
▃███▅▂▂▂▂▂▂▂▂▂▁▂▂▂▂▁▂▂▃▃▂▂▄███▅▅▆█▅▄▄▅▄▃▃▄▃▃▂▂▂▂▃▃▃▂▂▃▂▂▂▂▂▂▂ ▃
388 ns Histogram: frequency by time 527 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.AlefeldPotraShi()))BenchmarkTools.Trial: 10000 samples with 192 evaluations per sample.
Range (min … max): 517.917 ns … 203.366 μs ┊ GC (min … max): 0.00% … 99.52%
Time (median): 582.333 ns ┊ GC (median): 0.00%
Time (mean ± σ): 598.661 ns ± 2.182 μs ┊ GC (mean ± σ): 4.72% ± 1.40%
▄█▆▃ ▂▂▁▁ ▅▇▇▆▅▄▂▁▂▃▅▄▄▃▁▁ ▁▃▃▂▂▁▁ ▁ ▂
█████▆▇▇▇▆▄▄▁▁▁▁▃▁▄▁███████████████████████▇▇▄▆████████▇█████ █
518 ns Histogram: log(frequency) by time 660 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
