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 264 evaluations.
Range (min … max): 296.966 ns … 31.867 μs ┊ GC (min … max): 0.00% … 98.67%
Time (median): 302.269 ns ┊ GC (median): 0.00%
Time (mean ± σ): 340.374 ns ± 621.954 ns ┊ GC (mean ± σ): 3.64% ± 1.97%
█▇
▃██▅▃▂▂▂▂▂▂▂▂▁▂▁▁▁▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▃▄▄▄▅▅▄▄▃▄▄▃▃▂▃▂▂▂▂▂▂▂▂▂ ▃
297 ns Histogram: frequency by time 385 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.Goat())BenchmarkTools.Trial: 10000 samples with 11 evaluations.
Range (min … max): 981.727 ns … 2.664 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.000 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.060 μs ± 96.976 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▇ █▅ ▁ ▄▆ ▅ ▄▃ ▂▁ ▁▁ ▂
██▁██▁█▁██▁▇▇▁▇▁▇▆▁▅▁▆▄▁▁▅▁▆▁██▁█▁██▁██▁█▁▇▆▁▁▁▇█▁██▁▇▁▇▇▁▇▆ █
982 ns Histogram: log(frequency) by time 1.32 μs <
Memory estimate: 128 bytes, allocs estimate: 1.using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Newton()))BenchmarkTools.Trial: 10000 samples with 209 evaluations.
Range (min … max): 361.244 ns … 43.344 μs ┊ GC (min … max): 0.00% … 98.73%
Time (median): 390.675 ns ┊ GC (median): 0.00%
Time (mean ± σ): 405.767 ns ± 740.518 ns ┊ GC (mean ± σ): 3.15% ± 1.71%
█▇
▂███▇▃▂▁▁▁▁▁▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▂▃▄▃▃▃▄▆▆▅▄▄▅▅▄▃▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁ ▂
361 ns Histogram: frequency by time 449 ns <
Memory estimate: 128 bytes, allocs estimate: 1.using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Thukral3B()))BenchmarkTools.Trial: 10000 samples with 200 evaluations.
Range (min … max): 410.500 ns … 46.184 μs ┊ GC (min … max): 0.00% … 98.70%
Time (median): 461.505 ns ┊ GC (median): 0.00%
Time (mean ± σ): 462.559 ns ± 783.482 ns ┊ GC (mean ± σ): 2.92% ± 1.71%
▅█
▁▁▂▃▇██▆▃▁▁▁▁▂▂▂▃▃▂▂▂▁▁▁▁▁▁▁▁▁▁▂▃▇▅▄▄▅▆▆▅▄▃▃▃▃▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁ ▂
410 ns Histogram: frequency by time 507 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.A42()))BenchmarkTools.Trial: 10000 samples with 20 evaluations.
Range (min … max): 954.950 ns … 2.130 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.055 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.032 μs ± 60.213 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█ ▇ ▁█ ▄▂
▂▅█▁█▆▁▃▂▁▂▂▂▁▂▂▁▂▂▁▂▂▂▁▂▃▁██▄██▇▁▅▃▁▃▂▂▂▂▂▁▂▃▁▃▃▂▃▃▂▁▂▂▁▂▂▂ ▃
955 ns Histogram: frequency by time 1.16 μs <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Bisection()))BenchmarkTools.Trial: 10000 samples with 9 evaluations.
Range (min … max): 2.055 μs … 4.622 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 2.133 μs ┊ GC (median): 0.00%
Time (mean ± σ): 2.182 μs ± 123.428 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃ █▄
▂▄█▁██▇▁▇▆▆▁▅▄▁▄▃▃▁▃▃▃▁▃▃▁▃▄▆▁▅▅▄▁▄▄▁▃▃▃▁▂▂▂▁▂▂▁▂▂▂▁▂▂▂▁▂▂▂ ▃
2.06 μs Histogram: frequency by time 2.53 μs <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.SuperHalley()))BenchmarkTools.Trial: 10000 samples with 206 evaluations.
Range (min … max): 371.354 ns … 45.776 μs ┊ GC (min … max): 0.00% … 98.72%
Time (median): 421.354 ns ┊ GC (median): 0.00%
Time (mean ± σ): 420.081 ns ± 781.735 ns ┊ GC (mean ± σ): 3.21% ± 1.71%
██
▃██▇▂▁▁▁▁▁▁▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▄▄▄▄▄▅▅▄▄▅▄▄▃▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁ ▂
371 ns Histogram: frequency by time 465 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Brent()))BenchmarkTools.Trial: 10000 samples with 116 evaluations.
Range (min … max): 760.336 ns … 81.160 μs ┊ GC (min … max): 0.00% … 98.86%
Time (median): 812.069 ns ┊ GC (median): 0.00%
Time (mean ± σ): 819.221 ns ± 1.130 μs ┊ GC (mean ± σ): 1.95% ± 1.40%
▁█▄ ▃ ▃█▅
▃███▆▂▁▁▂▂▂▂▃▂▂▁▁▁▁▃▆█▆▄▂▃▄▇███▆▄▄█▆▅▄▃▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
760 ns Histogram: frequency by time 876 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Order2()))BenchmarkTools.Trial: 10000 samples with 195 evaluations.
Range (min … max): 477.431 ns … 49.433 μs ┊ GC (min … max): 0.00% … 98.71%
Time (median): 533.328 ns ┊ GC (median): 0.00%
Time (mean ± σ): 531.058 ns ± 837.711 ns ┊ GC (mean ± σ): 2.72% ± 1.71%
██▁ ▃
▂████▃▂▁▁▁▁▂▃▃▃▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▃▆█▆▅▇▇▇▆▄▃▃▃▃▃▁▁▂▁▁▁▁▁▁▁▁▁ ▂
477 ns Histogram: frequency by time 574 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.AlefeldPotraShi()))BenchmarkTools.Trial: 10000 samples with 58 evaluations.
Range (min … max): 867.241 ns … 160.079 μs ┊ GC (min … max): 0.00% … 99.32%
Time (median): 925.862 ns ┊ GC (median): 0.00%
Time (mean ± σ): 940.303 ns ± 1.592 μs ┊ GC (mean ± σ): 1.69% ± 0.99%
█▁▆
▂▄▂▂▁▁▁▂▃▆▅▂▂▁▁▁▁▂▇████▃▃▂▂▃▄▂▃▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
867 ns Histogram: frequency by time 1.04 μs <
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)18238.74783586164058110436412157446724913271569356796169848820891801965891606768Comparison
