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 626 evaluations.
Range (min … max): 193.877 ns … 17.662 μs ┊ GC (min … max): 0.00% … 98.65%
Time (median): 198.422 ns ┊ GC (median): 0.00%
Time (mean ± σ): 231.545 ns ± 402.007 ns ┊ GC (mean ± σ): 5.11% ± 2.95%
█▅▂▁ ▁ ▃▄▅▆▅▅▃▂▁ ▂
████▇▅▅▅▅▄▄██▇▆▆▆▄▃▄▁▁▁▄▁▃▃▄▁▃▁▁▃▁▃██████████▇▇██████▇▆▅▅▅▄▆▇ █
194 ns Histogram: log(frequency) by time 272 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.Goat())BenchmarkTools.Trial: 10000 samples with 163 evaluations.
Range (min … max): 659.577 ns … 45.950 μs ┊ GC (min … max): 0.00% … 98.52%
Time (median): 710.528 ns ┊ GC (median): 0.00%
Time (mean ± σ): 708.504 ns ± 640.549 ns ┊ GC (mean ± σ): 1.28% ± 1.39%
██▄ ▁ ▅█▅▇▅▅▃▂▁ ▁ ▁ ▂
█████▇▃▄▃▃▁███████████▇███████▇▅▅▅▅▅▅▁▃▁▅▃▄▁▄▃▁▅▃▁▁▁▄▃▃▃▅▃▅▆▆ █
660 ns Histogram: log(frequency) by time 916 ns <
Memory estimate: 128 bytes, allocs estimate: 1.using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Newton()))BenchmarkTools.Trial: 10000 samples with 394 evaluations.
Range (min … max): 245.000 ns … 19.502 μs ┊ GC (min … max): 0.00% … 98.41%
Time (median): 293.301 ns ┊ GC (median): 0.00%
Time (mean ± σ): 287.706 ns ± 461.935 ns ┊ GC (mean ± σ): 3.92% ± 2.41%
█▃▁ ▃▄▆▆▅▄▁ ▂
███▆▄▅▅▄██▇▇▆▅▄▄█████████████▇▇▆▅▅▅▅▅▄▁▅▃▃▃▃▃▃▃▃▄▁▁▁▃▁▁▃▁▁▄▅▅ █
245 ns Histogram: log(frequency) by time 414 ns <
Memory estimate: 128 bytes, allocs estimate: 1.using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Thukral3B()))BenchmarkTools.Trial: 10000 samples with 298 evaluations.
Range (min … max): 276.658 ns … 26.427 μs ┊ GC (min … max): 0.00% … 98.83%
Time (median): 323.945 ns ┊ GC (median): 0.00%
Time (mean ± σ): 318.779 ns ± 519.526 ns ┊ GC (mean ± σ): 3.26% ± 1.98%
▆█▅▃ ▄▄▂▄▆▅▄▅▅▄▄▃▂▂▁▁ ▂
██████▇█▇▅▃▃▃▄▄▁▃▄▁▄▁▃▇█▇▆▆▆▆▆▅▇███████████████████▇█▇▇▇▆▆▇▆▇ █
277 ns Histogram: log(frequency) by time 365 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.A42()))BenchmarkTools.Trial: 10000 samples with 106 evaluations.
Range (min … max): 775.792 ns … 75.566 μs ┊ GC (min … max): 0.00% … 98.89%
Time (median): 824.179 ns ┊ GC (median): 0.00%
Time (mean ± σ): 827.691 ns ± 748.277 ns ┊ GC (mean ± σ): 0.90% ± 0.99%
▄█▆▂▂ ▂▇▇▄▂▁▄▇▆▄▂▂▂▄▄▁ ▁▁ ▂
██████▅▅▂▆▆▅▄▃▃▂▂████████████████████▇▆▅▄▆█▆▇▅▅▆▆▆▅▆▇█▆▆▅▆▆▅▅ █
776 ns Histogram: log(frequency) by time 929 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.302 μs … 3.620 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.312 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.372 μs ± 124.219 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▇█▆▂ ▃▅▅▅▃▂▁ ▁▁ ▂
██████▇▆▆▆▆▆▇█▇▅▃▄▁▁▁▁▁▆████████▇▅▅▄▆▅▆▅▅▆▆▆▆▇▇▆▆████▇▆▄▄▃▃ █
1.3 μs Histogram: log(frequency) by time 1.68 μs <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.SuperHalley()))BenchmarkTools.Trial: 10000 samples with 379 evaluations.
Range (min … max): 246.715 ns … 21.068 μs ┊ GC (min … max): 0.00% … 98.69%
Time (median): 294.931 ns ┊ GC (median): 0.00%
Time (mean ± σ): 286.049 ns ± 412.018 ns ┊ GC (mean ± σ): 2.87% ± 1.97%
██▆▂▁▁ ▁ ▅▄▅▇▆▆▆▄▃▂▁ ▁ ▃
██████▇▆▄▅▅▅▄▃▃▁▅█▇▇▆▇▇▆▆▅▅▁▅█▇██████████████▇██▇█▇▇██▇▇▆▆▆▅▆ █
247 ns Histogram: log(frequency) by time 337 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Brent()))BenchmarkTools.Trial: 10000 samples with 194 evaluations.
Range (min … max): 495.768 ns … 41.319 μs ┊ GC (min … max): 0.00% … 98.70%
Time (median): 545.711 ns ┊ GC (median): 0.00%
Time (mean ± σ): 546.317 ns ± 699.474 ns ┊ GC (mean ± σ): 2.22% ± 1.71%
▅█▆▄▂ ▆▇▅▃▅▇▆▄▄▄▄▃▁▁▁▁ ▂
█████▇▆▇▇▆▅▆▃▁▁▁▁▅▃▁▁▃▃▁▃▇█████████████████▇▇█▇▆▇▇█▆▇██▇▆▇▇▇▇ █
496 ns Histogram: log(frequency) by time 606 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Order2()))BenchmarkTools.Trial: 10000 samples with 237 evaluations.
Range (min … max): 314.557 ns … 33.439 μs ┊ GC (min … max): 0.00% … 99.03%
Time (median): 360.502 ns ┊ GC (median): 0.00%
Time (mean ± σ): 356.455 ns ± 564.099 ns ┊ GC (mean ± σ): 2.74% ± 1.71%
█
▆█▅▃▂▂▂▂▂▂▂▂▁▂▂▂▂▂▁▁▂▂▂▁▂▂▂▂▃▅▃▃▄▆▄▃▄▅▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▂
315 ns Histogram: frequency by time 410 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.AlefeldPotraShi()))BenchmarkTools.Trial: 10000 samples with 153 evaluations.
Range (min … max): 679.438 ns … 52.432 μs ┊ GC (min … max): 0.00% … 98.65%
Time (median): 730.974 ns ┊ GC (median): 0.00%
Time (mean ± σ): 727.917 ns ± 723.331 ns ┊ GC (mean ± σ): 1.40% ± 1.40%
▆█▅▃▂ ▅▃▁▁▅▇▅▄▄▅▆▄▃▃▃▂▁ ▂
█████▆▅▆█▇▆▄▄▃▄▄▁▁▁███████████████████▇██▇▆▅▇██▆▇▆█▇▆▇█▇▆▆▇▇▇ █
679 ns Histogram: log(frequency) by time 805 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
