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 (minmax):  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 (minmax):  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 (minmax):  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 (minmax):  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 (minmax):  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 (minmax):  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 (minmax):  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 (minmax):  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 (minmax):  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 (minmax):  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