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 550 evaluations.
 Range (minmax):  209.227 ns 20.209 μs   GC (min … max): 0.00% … 98.93%
 Time  (median):     216.114 ns                GC (median):    0.00%
 Time  (mean ± σ):   246.739 ns ± 407.886 ns   GC (mean ± σ):  4.31% ±  2.60%

  ██▃ ▁          ▄▅▅▆▆▅▅▄▃▂▁▁          ▂
  ███▇█▇▄▅▃▄▄▁▃████▇▇▆▆▅▃▃▃▁▄▁▃▁▁▃▃▃▁▁▁▄██████████████▇██████ █
  209 ns        Histogram: log(frequency) by time        282 ns <

 Memory estimate: 128 bytes, allocs estimate: 1.
@benchmark orbitsolve(orb, t, PlanetOrbits.Goat())
BenchmarkTools.Trial: 10000 samples with 159 evaluations.
 Range (minmax):  676.107 ns 48.166 μs   GC (min … max): 0.00% … 98.47%
 Time  (median):     727.969 ns                GC (median):    0.00%
 Time  (mean ± σ):   726.302 ns ± 666.319 ns   GC (mean ± σ):  1.29% ±  1.39%

  ▆█▄         ▅▇▅▄▃▂       ▁                                  ▂
  █████▇▅▄▃▄▃███████████▆▇██▇█▇▇▆▅▆▆▄▄▄▄▃▄▄▁▁▄▄▁▃▃▄▃▃▁▃▃▄▄▅▃▆ █
  676 ns        Histogram: log(frequency) by time        924 ns <

 Memory estimate: 128 bytes, allocs estimate: 1.
using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Newton()))
BenchmarkTools.Trial: 10000 samples with 337 evaluations.
 Range (minmax):  262.448 ns 23.020 μs   GC (min … max): 0.00% … 98.60%
 Time  (median):     309.467 ns                GC (median):    0.00%
 Time  (mean ± σ):   304.754 ns ± 503.591 ns   GC (mean ± σ):  3.69% ±  2.21%

  █▅                                                             
  ██▃▂▂▂▂▂▂▂▂▂▂▁▁▂▂▁▂▂▂▂▂▂▂▂▂▂▂▃▃▄▆▄▅▇▅▄▄▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
  262 ns           Histogram: frequency by time          352 ns <

 Memory estimate: 128 bytes, allocs estimate: 1.
using Roots
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Thukral3B()))
BenchmarkTools.Trial: 10000 samples with 266 evaluations.
 Range (minmax):  293.861 ns 29.426 μs   GC (min … max): 0.00% … 98.97%
 Time  (median):     342.068 ns                GC (median):    0.00%
 Time  (mean ± σ):   341.493 ns ± 573.470 ns   GC (mean ± σ):  3.35% ±  1.98%

  █▆           ▅▇▆▆▄▃▂▁                                        ▂
  ████▇▄▄▃▃▃▇█▆█████████▇██▇▇▇▆▆▆▆▅▅▄▅▅▅▅▁▄▄▄▄▃▁▁▃▃▁▁▁▄▁▃▄▅▄▆▆ █
  294 ns        Histogram: log(frequency) by time        506 ns <

 Memory estimate: 128 bytes, allocs estimate: 1.
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.A42()))
BenchmarkTools.Trial: 10000 samples with 117 evaluations.
 Range (minmax):  757.137 ns 68.163 μs   GC (min … max): 0.00% … 98.78%
 Time  (median):     812.718 ns                GC (median):    0.00%
 Time  (mean ± σ):   813.044 ns ± 675.075 ns   GC (mean ± σ):  0.83% ±  0.99%

  █▇▂     ▅▇▂▇▄▆▆▃▃▃▁▁  ▁                                      ▃
  ███▆▆▇▅▁████████████▇▆██▇████▇▇█▇▆▇▇▅▅▅▅▅▅▃▄▁▃▁▃▁▄▃▁▆▅▆▆▆▆▅▆ █
  757 ns        Histogram: log(frequency) by time       1.04 μ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.295 μs  2.895 μs   GC (min … max): 0.00% … 0.00%
 Time  (median):     1.306 μs                GC (median):    0.00%
 Time  (mean ± σ):   1.364 μs ± 114.718 ns   GC (mean ± σ):  0.00% ± 0.00%

  ▇█               ▄▅▅▄▃▂▁                    ▁      ▂
  ███▇▆▄▅▆▆▄▁▄▆▃▃▁▁▁▁▁▁▁▁▇████████▇▇▅▅▃▅▄▄▁▅▅▁▃▅▄▄▄▇████▇▇▆ █
  1.3 μs       Histogram: log(frequency) by time      1.65 μs <

 Memory estimate: 128 bytes, allocs estimate: 1.
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.SuperHalley()))
BenchmarkTools.Trial: 10000 samples with 328 evaluations.
 Range (minmax):  266.320 ns 24.354 μs   GC (min … max): 0.00% … 98.67%
 Time  (median):     314.948 ns                GC (median):    0.00%
 Time  (mean ± σ):   311.908 ns ± 530.535 ns   GC (mean ± σ):  3.80% ±  2.21%

  █▄              ▁▆▆▅▄▃▂▁▁                                   ▂
  ██▆█▆▁▃▄▁▄▇█▆▆▅▅█████████████▇▇▇▆▅▅▅▅▅▃▄▁▅▃▁▄▄▄▁▄▁▁▃▁▃▄▅▅▅▆ █
  266 ns        Histogram: log(frequency) by time        432 ns <

 Memory estimate: 128 bytes, allocs estimate: 1.
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Brent()))
BenchmarkTools.Trial: 10000 samples with 8 evaluations.
 Range (minmax):  3.825 μs  7.783 μs   GC (min … max): 0.00% … 0.00%
 Time  (median):     3.876 μs                GC (median):    0.00%
 Time  (mean ± σ):   3.953 μs ± 206.535 ns   GC (mean ± σ):  0.00% ± 0.00%

  ▄▇█▅▂    ▁▄▅▅▄▂▁     ▁▁                                    ▂
  █████▇▅▄███████▆▆▆▆█████▆▅▅▃▄▄▄▃▃▁▃▄▄▃▃▄▅▄▅▅▃▄▄▄▁▃▃▆▇▇▇▇▇ █
  3.82 μs      Histogram: log(frequency) by time      4.97 μs <

 Memory estimate: 128 bytes, allocs estimate: 1.
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.Order2()))
BenchmarkTools.Trial: 10000 samples with 221 evaluations.
 Range (minmax):  333.561 ns 35.921 μs   GC (min … max): 0.00% … 99.02%
 Time  (median):     382.525 ns                GC (median):    0.00%
 Time  (mean ± σ):   377.373 ns ± 614.779 ns   GC (mean ± σ):  2.82% ±  1.71%

  ▁█                                                             
  ██▅▃▂▂▂▂▂▂▁▁▁▁▂▂▁▂▁▁▂▂▂▂▂▂▂▂▄▅▃▄▇▆▄▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
  334 ns           Histogram: frequency by time          433 ns <

 Memory estimate: 128 bytes, allocs estimate: 1.
@benchmark orbitsolve(orb, t, PlanetOrbits.RootsMethod(Roots.AlefeldPotraShi()))
BenchmarkTools.Trial: 10000 samples with 147 evaluations.
 Range (minmax):  696.068 ns 57.160 μs   GC (min … max): 0.00% … 98.74%
 Time  (median):     749.565 ns                GC (median):    0.00%
 Time  (mean ± σ):   754.526 ns ± 797.959 ns   GC (mean ± σ):  1.49% ±  1.40%

  █▆▁         ▅▆▂▆▅▆▄▂▂▂▁                                     ▂
  ███▆▅▇▆▅▅▄▃▃█████████████▇██▇▇█▇██▇▇▇▇▇▇▇▇▇▆▆▆▆▆▆▆▆▅▆▅▆▆▆▅▅ █
  696 ns        Histogram: log(frequency) by time        905 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)
25030.34323194365565831507041947847102607642627168092403020601399008600937133855

Comparison