Give kepler a position and a velocity and it tells you where the satellite will be later. It uses the two-body model: the Earth is a point mass and drag is ignored. That is fine over the short horizons a ground-station planner or an orbit viewer cares about.

There are two ways to propagate. The analytic one turns the state vector into orbital elements first. From there the mean anomaly grows linearly, M = M0 + n t, with n = sqrt(mu / a^3), where mu is the Earth's gravitational parameter and a the semi-major axis. To get back to a position you need the eccentric anomaly E, which comes from Kepler's equation M = E - e sin E (for e < 1). The solver uses Newton-Raphson, E <- E - (E - e sin E - M) / (1 - e cos E), starting at E = M + e sin M. For e < 0.97 three to five steps would do, but it always takes eight, so each position costs exactly eight sines and cosines. This is algorithm 2 in Vallado's Fundamentals of Astrodynamics.

The other way integrates the equations of motion with RK4. Orbits that are not bound (e >= 1) always go this way, even through the analytic entry point. bench/ times the two against each other.

Units are km and km/s in an Earth-centred frame. The Earth's gravitational parameter is the only one the library has; for another central body you supply your own.

Installation #

Install with opam:

$ opam install kepler

If opam cannot find the package, it is not in the public opam-repository yet. Add the overlay repository first:

$ opam repo add samoht https://tangled.org/gazagnaire.org/opam-overlay.git
$ opam update
$ opam install kepler

Usage #

let propagate pos vel =
  (* Analytic propagation with precomputed elements (fast path) *)
  let elems = Kepler.Analytic.precompute ~pos ~vel in
  let new_pos = Kepler.Analytic.at_precomputed elems ~dt:3600. in
  Printf.printf "Period: %.1f s\n" (Kepler.Analytic.period elems);
  (* Generate an orbit arc for visualization *)
  let arc =
    Kepler.Analytic.arc ~pos ~vel ~duration_s:5400. ~num_points:360
  in
  (new_pos, arc)

API Overview #

In Kepler.Analytic, call precompute once per satellite to turn its state vector into elements, then at_precomputed for each time you want. at skips the precomputation and takes a state vector directly (analytic for e < 1, RK4 otherwise). arc returns an array of positions over a time span centred on the epoch. period, eccentricity, semi_major_axis and inclination read the elements, and solve_kepler solves M = E - e sin E on its own.

Kepler.Propagate is the RK4 side: at steps about 30 seconds at a time, and arc uses the same integrator. mu is the Earth's gravitational parameter, 398600.4418 km^3/s^2.