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 641 evaluations per sample.
Range (min … max): 191.089 ns … 18.404 μs ┊ GC (min … max): 0.00% … 98.80%
Time (median): 208.102 ns ┊ GC (median): 0.00%
Time (mean ± σ): 235.105 ns ± 439.454 ns ┊ GC (mean ± σ): 5.80% ± 3.11%
▃█▇▂▁ ▄▄ ▃▄▅▆▆▅▅▄▃▂▁ ▁▂▃▃▂▂▂▁ ▂
█████▇▃▅▅▅▅▆███▇▇▆▆▄▁▃▃▃▅▄▄▄▃▁▁▁▁▁▁▁▁▁▁▄████████████████████▇ █
191 ns Histogram: log(frequency) by time 267 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.Goat())BenchmarkTools.Trial: 10000 samples with 160 evaluations per sample.
Range (min … max): 658.913 ns … 50.117 μs ┊ GC (min … max): 0.00% … 98.60%
Time (median): 708.194 ns ┊ GC (median): 0.00%
Time (mean ± σ): 709.038 ns ± 692.729 ns ┊ GC (mean ± σ): 1.38% ± 1.39%
█▅ ▂▁
▄██▆▃▂▂▂▂▂▂▂▂▁▁▂▁▂▂▁▂▁▂▅██▆▄▅▇▇▄▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▃▃▂▂▂▃▃▂▂▂▂ ▃
659 ns Histogram: frequency by time 781 ns <
Memory estimate: 128 bytes, allocs estimate: 1.using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Newton()))BenchmarkTools.Trial: 10000 samples with 383 evaluations per sample.
Range (min … max): 247.013 ns … 21.596 μs ┊ GC (min … max): 0.00% … 98.55%
Time (median): 296.871 ns ┊ GC (median): 0.00%
Time (mean ± σ): 288.909 ns ± 417.118 ns ┊ GC (mean ± σ): 2.88% ± 1.97%
█▇▁ ▂▄▁ ▂▄▄▆▆▆▅▄▃▁ ▁▂▂▂▂▁ ▂
█████▆▃▄▄▁▃▃███▇▆▆▅▅▃▁▄▄▄▃███████████▅▄▇▇███████▆▆▄▁▅▃▄▆▄▅▅▄▄ █
247 ns Histogram: log(frequency) by time 356 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 per sample.
Range (min … max): 275.413 ns … 27.052 μs ┊ GC (min … max): 0.00% … 98.80%
Time (median): 323.490 ns ┊ GC (median): 0.00%
Time (mean ± σ): 316.593 ns ± 375.856 ns ┊ GC (mean ± σ): 1.67% ± 1.40%
█▅ ▃▂ ▄▄▆▅▅▄▄▂ ▁▁▁▁ ▁
██▆▇▆▆▃▂▃▄▇██▇▆▅▆██████████▆▆███████▇▆▆▅▄▄▅▅▅▆▄▄▄▄▄▄▆▄▃▅▅▃▅▅▄ █
275 ns Histogram: log(frequency) by time 434 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.A42()))BenchmarkTools.Trial: 10000 samples with 113 evaluations per sample.
Range (min … max): 764.965 ns … 71.771 μs ┊ GC (min … max): 0.00% … 98.88%
Time (median): 823.124 ns ┊ GC (median): 0.00%
Time (mean ± σ): 822.985 ns ± 710.668 ns ┊ GC (mean ± σ): 0.86% ± 0.99%
▅█▇▂▁ ▄▇▆▄▂▂▆█▆▄▃▄▆▅▄▂▂▂▁ ▁▁▁ ▁▂▂▁▁ ▁▁ ▃
█████▅▅▅▆▇▆▄▃▃▁▃████████████████████▇▇▄▅▅▅▃▅▇████▇██████████▇ █
765 ns Histogram: log(frequency) by time 923 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Bisection()))BenchmarkTools.Trial: 10000 samples with 10 evaluations per sample.
Range (min … max): 1.301 μs … 3.585 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.313 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.379 μs ± 141.774 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▅ ▅▅▃▁ ▁▁ ▁
███▇▇▅▅▅▃▃▇████▇▆▅▅▅▇▇███▇▅▄▄▄▃▄▄▄▁▃▁▃▁▁▁▃▃▁▁▁▁▁▁▃▁▁▃▃▃▅▆▆▇ █
1.3 μs Histogram: log(frequency) by time 2.14 μs <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.SuperHalley()))BenchmarkTools.Trial: 10000 samples with 383 evaluations per sample.
Range (min … max): 247.433 ns … 21.506 μs ┊ GC (min … max): 0.00% … 98.68%
Time (median): 296.689 ns ┊ GC (median): 0.00%
Time (mean ± σ): 290.150 ns ± 464.982 ns ┊ GC (mean ± σ): 3.58% ± 2.20%
█▇▂ ▃▃▁ ▃▄▄▆▆▆▆▅▄▂▁▁ ▁▂▂▂▂▁ ▂
███▇██▆▃▃▃▃▃▁▃▁▅████▇▅▅▄▅▁▄▄▄▃▆▄▅█████████████▇▆▅▅▇█████████▇ █
247 ns Histogram: log(frequency) by time 332 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Brent()))BenchmarkTools.Trial: 10000 samples with 189 evaluations per sample.
Range (min … max): 515.884 ns … 43.925 μs ┊ GC (min … max): 0.00% … 98.72%
Time (median): 568.307 ns ┊ GC (median): 0.00%
Time (mean ± σ): 566.679 ns ± 608.950 ns ┊ GC (mean ± σ): 1.52% ± 1.40%
▃█ ▁▁ ▁
██▄▂▂▂▂▂▂▂▂▁▁▁▁▂▁▁▂▂▁▂▃▄▃▂▃██▆▄▅█▆▄▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▃▃▃▂▃▃▂▂▂▂ ▃
516 ns Histogram: frequency by time 630 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Order2()))BenchmarkTools.Trial: 10000 samples with 236 evaluations per sample.
Range (min … max): 315.377 ns … 51.315 μs ┊ GC (min … max): 0.00% … 99.21%
Time (median): 363.263 ns ┊ GC (median): 0.00%
Time (mean ± σ): 362.365 ns ± 704.316 ns ┊ GC (mean ± σ): 3.30% ± 1.72%
▇█▄ ▃▃ ▃▅▅▄▅▆▆▅▄▅▅▄▃▃▁ ▁▁▁▁▁▁ ▂
████▇▆▇▇▅▅▃▃▁▁▃▃▃▅▁▁▄████▇▇████████████████▇▅▆▅▅▅▇██████████▇ █
315 ns Histogram: log(frequency) by time 417 ns <
Memory estimate: 128 bytes, allocs estimate: 1.@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.AlefeldPotraShi()))BenchmarkTools.Trial: 10000 samples with 158 evaluations per sample.
Range (min … max): 666.430 ns … 52.780 μs ┊ GC (min … max): 0.00% … 98.62%
Time (median): 717.981 ns ┊ GC (median): 0.00%
Time (mean ± σ): 720.479 ns ± 724.196 ns ┊ GC (mean ± σ): 1.42% ± 1.39%
█▄ █ ▁
██▂▁▁▁▁▁▁▁▁▁▁▁▄█▅▄█▃▂▂▁▁▁▁▁▁▁▁▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
666 ns Histogram: frequency by time 872 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
