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 646 evaluations.
Range (min … max): 191.238 ns … 17.318 μs ┊ GC (min … max): 0.00% … 98.66%
Time (median): 199.124 ns ┊ GC (median): 0.00%
Time (mean ± σ): 231.055 ns ± 414.814 ns ┊ GC (mean ± σ): 5.58% ± 3.11%
█▇▂▁ ▁▁ ▃▄▆▆▆▅▅▄▃▁ ▂
█████▇▆▆▆▃▅███▇▅▅▃▄▄▅▆▁▄▄▁▃▃▃▃▃▁▁▁▃▅████████████████▇█▇▅▅▅▇▆▇ █
191 ns Histogram: log(frequency) by time 271 ns <
Memory estimate: 128 bytes, allocs estimate: 1.
@benchmark orbitsolve(orb, t, PlanetOrbits.Goat())
BenchmarkTools.Trial: 10000 samples with 163 evaluations.
Range (min … max): 650.043 ns … 48.111 μs ┊ GC (min … max): 0.00% … 98.60%
Time (median): 699.337 ns ┊ GC (median): 0.00%
Time (mean ± σ): 703.324 ns ± 672.750 ns ┊ GC (mean ± σ): 1.35% ± 1.39%
█▆ ▁ ▇▆▆▅▃▂ ▁ ▂
████▅▄▃▇████████▇█████▆▆▆▅▆▅▅▆▅▅▅▄▄▄▁▃▃▄▃▄▃▁▁▁▄▁▃▄▄▄▃▆▆▆▇▇▇▆▆ █
650 ns Histogram: log(frequency) by time 1.01 μs <
Memory estimate: 128 bytes, allocs estimate: 1.
using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Newton()))
BenchmarkTools.Trial: 10000 samples with 379 evaluations.
Range (min … max): 248.934 ns … 20.466 μs ┊ GC (min … max): 0.00% … 98.54%
Time (median): 296.040 ns ┊ GC (median): 0.00%
Time (mean ± σ): 286.804 ns ± 400.750 ns ┊ GC (mean ± σ): 2.79% ± 1.97%
▇█▄▂ ▁ ▄▃▃▇▆▅▆▅▄▃▂▁ ▂
████▇▇█▇▅▄▄▄▄▄▁▄▃▅▇█▇▆▅▅▄▄▄▅▁▃▄▇▇██████████████▇▆▆▆▆▇▇▆▇█▇▇▇█ █
249 ns Histogram: log(frequency) by time 332 ns <
Memory estimate: 128 bytes, allocs estimate: 1.
using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Thukral3B()))
BenchmarkTools.Trial: 10000 samples with 301 evaluations.
Range (min … max): 277.993 ns … 25.636 μs ┊ GC (min … max): 0.00% … 98.81%
Time (median): 324.492 ns ┊ GC (median): 0.00%
Time (mean ± σ): 319.800 ns ± 504.902 ns ┊ GC (mean ± σ): 3.15% ± 1.97%
▆█▄▂ ▄▄▂▅▆▄▄▆▅▄▄▃▂▁▁ ▂
▅████▇▆▇▇▆▆▆▄▄▄▄▄▄▁█▇▅▅▇█▆▅▄▆▅▇█████████████████▇▅▇▇▅▇▅▇▆▇▆▆▇ █
278 ns Histogram: log(frequency) by time 367 ns <
Memory estimate: 128 bytes, allocs estimate: 1.
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.A42()))
BenchmarkTools.Trial: 10000 samples with 108 evaluations.
Range (min … max): 773.667 ns … 71.302 μs ┊ GC (min … max): 0.00% … 98.53%
Time (median): 820.880 ns ┊ GC (median): 0.00%
Time (mean ± σ): 828.744 ns ± 706.150 ns ┊ GC (mean ± σ): 0.85% ± 0.99%
▂█ ▆▂ ▇
██▃▂▂▂▂▂▂▂██▃▆█▄▃▆▄▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
774 ns Histogram: frequency by time 1.03 μs <
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.301 μs … 3.989 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.309 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.370 μs ± 126.844 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▇▄▁ ▃▅▅▄▄▂▁ ▁ ▂
█████▇▇▇▅▅▆▅█▆▅▃▃▁▁▃▁▁▁████████▇▆▅▅▄▆▄▅▃▃▄▅▃▅▆▄▇████▇▇▆▆▅▃▄ █
1.3 μs Histogram: log(frequency) by time 1.69 μs <
Memory estimate: 128 bytes, allocs estimate: 1.
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.SuperHalley()))
BenchmarkTools.Trial: 10000 samples with 394 evaluations.
Range (min … max): 246.881 ns … 20.244 μs ┊ GC (min … max): 0.00% … 98.53%
Time (median): 296.236 ns ┊ GC (median): 0.00%
Time (mean ± σ): 290.099 ns ± 484.201 ns ┊ GC (mean ± σ): 4.08% ± 2.41%
▅█▅▂ ▃▂▂▆▅▄▆▅▄▄▄▃▂▁ ▂
████▇▇█▇▄▄▄▃▅▄▄▃▅▄██▇▆▄▅▄▄▄▄▄▃▁▄▅███████████████▇▇█▆▆▆▇▇▇▇▇▇▇ █
247 ns Histogram: log(frequency) by time 330 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): 493.861 ns … 41.254 μs ┊ GC (min … max): 0.00% … 98.69%
Time (median): 541.990 ns ┊ GC (median): 0.00%
Time (mean ± σ): 539.265 ns ± 574.419 ns ┊ GC (mean ± σ): 1.50% ± 1.40%
█ ▁
▇█▅▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▂▁▂▂██▄▃▆█▆▄▃▄▄▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
494 ns Histogram: 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 238 evaluations.
Range (min … max): 315.294 ns … 33.721 μs ┊ GC (min … max): 0.00% … 98.93%
Time (median): 360.588 ns ┊ GC (median): 0.00%
Time (mean ± σ): 357.660 ns ± 572.170 ns ┊ GC (mean ± σ): 2.77% ± 1.71%
█
▆█▄▃▂▂▂▂▂▂▂▂▂▂▂▁▂▁▂▁▁▂▂▂▂▂▃▄▃▃▅▆▄▃▅▅▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▂
315 ns Histogram: frequency by time 416 ns <
Memory estimate: 128 bytes, allocs estimate: 1.
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.AlefeldPotraShi()))
BenchmarkTools.Trial: 10000 samples with 154 evaluations.
Range (min … max): 678.474 ns … 52.005 μs ┊ GC (min … max): 0.00% … 98.65%
Time (median): 728.630 ns ┊ GC (median): 0.00%
Time (mean ± σ): 722.631 ns ± 513.997 ns ┊ GC (mean ± σ): 0.71% ± 0.99%
█▅ ▃▂
██▄▃▂▂▂▂▂▂▂▂▂▂▁▃▄▂▂██▄▃█▇▄▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
678 ns Histogram: frequency by time 834 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.95447361446997022436542900109713321175234563222583018047585031821959181256
Comparison