satkit 0.21.2

Satellite Toolkit
Documentation
# Keplerian Elements

This page describes the classical (Keplerian) orbital element set used by
[`satkit.kepler`](../api/kepler.md) (Rust: `satkit::kepler::Kepler`), the
conventions it follows, and what it does *not* do. For a worked notebook see
the [Keplerian Elements tutorial](../tutorials/Keplerian%20Elements.ipynb).

## The Element Set

A bound two-body orbit is described by six numbers. `satkit` stores them under
these names, in SI units and radians:

| Field   | Symbol     | Meaning                                              |
|---------|------------|------------------------------------------------------|
| `a`     | $a$        | semi-major axis, **meters**                          |
| `eccen` | $e$        | eccentricity, $0 \le e < 1$                          |
| `incl`  | $i$        | inclination, radians, $0 \le i \le \pi$              |
| `raan`  | $\Omega$   | right ascension of the ascending node, radians       |
| `w`     | $\omega$   | argument of perigee, radians                         |
| `nu`    | $\nu$      | true anomaly, radians                                |

The size and shape of the ellipse are $a$ and $e$; the orientation of the
orbital plane and of the ellipse within it are $i$, $\Omega$ and $\omega$; and
$\nu$ locates the satellite along the ellipse at the epoch of the elements.
The semiparameter (semi-latus rectum) $p = a(1 - e^2)$ is available as a
derived property, but the class is constructed from $a$, not $p$.

In Python the inclination property is spelled `inclination`; the constructor
argument and the Rust field are `incl`. Angles returned by `from_pv` are
reduced to $[0, 2\pi)$; angles you set are stored as given.

### Anomalies

Three angles can locate the satellite in its orbit
([Vallado 2013](references.md#vallado2013), §2.2):

- **True anomaly** $\nu$ — the angle at the focus (Earth's center) from
  perigee to the satellite. This is what is stored (`nu`).
- **Eccentric anomaly** $E$ — the angle at the *center* of the ellipse to the
  point on the auxiliary circle above the satellite. Related to $\nu$ by
  $\tan\frac{\nu}{2} = \sqrt{\frac{1+e}{1-e}}\tan\frac{E}{2}$.
- **Mean anomaly** $M$ — the angle that advances uniformly in time,
  $M = M_0 + n\,(t - t_0)$ with mean motion $n = \sqrt{\mu / a^3}$. Related to
  $E$ by Kepler's equation $M = E - e\sin E$.

Converting $M \to E$ requires solving Kepler's equation iteratively. `satkit`
uses Newton's method with the [Danby (1987)](references.md#danby1987) starting
value $E_0 = M + 0.85\,e\,\mathrm{sign}(\sin M)$ after reducing $M$ to
$[0, 2\pi)$, which converges in a handful of iterations for every $e < 1$
([Vallado 2013](references.md#vallado2013), Algorithm 2). The iteration is
capped, so a non-finite $M$ yields NaN rather than looping.

The class accepts any of the three anomalies when it is constructed, and
exposes all three as properties; `mean_anomaly` and `eccentric_anomaly` can
be assigned and are converted to `nu` on the spot.

## What the Elements Mean (and Don't)

**Osculating, not mean.** The elements are *osculating*: they describe the
two-body orbit that is tangent to the actual trajectory at the epoch of the
state. Under perturbations (oblateness, drag, third bodies, …) the osculating
elements vary continuously along the orbit — for a LEO satellite $a$ oscillates
by several kilometers within one revolution because of $J_2$ alone. They are
not the *mean* elements of an analytical theory. In particular, the elements
in a TLE are SGP4 mean elements and cannot be passed to `kepler` (or read
back from `from_pv`) without a significant, model-dependent error; use
[`satkit.TLE`](../api/tle.md) and the SGP4 propagator for those
(see [TLEs, SGP4 & OMMs](tle.md)).

**Frame.** `kepler` does no frame handling. `from_pv` interprets the position
and velocity you pass in whatever frame they are in, and `to_pv` returns the
state in that same frame. Elements are meaningful only in an inertial frame;
the rest of `satkit` assumes **GCRF**, so convert ITRF or TEME states with
[`frametransform`](../api/frametransform.md) first. Passing an Earth-fixed
state produces elements that are numerically valid but physically
meaningless.

**Central body.** The Earth's gravitational parameter
[`consts.MU_EARTH`](../api/consts.md) ($3.986004418 \times 10^{14}$
m³ s⁻²) is used everywhere: in the `from_pv` energy equation, in `to_pv`, and
in the mean motion, period and `propagate`. There is no way to use a
different $\mu$; for heliocentric or lunar orbits compute the elements
yourself.

**Closed orbits only.** `from_pv` returns an error (Python: `RuntimeError`)
for parabolic or hyperbolic states ($e \ge 1$) and for rectilinear states
(zero angular momentum, where the orbital plane is undefined). The
constructor itself does not validate $e$; supplying $e \ge 1$ produces
meaningless anomaly conversions.

**Singular cases.** $\Omega$ is undefined for an equatorial orbit and $\omega$
for a circular one. `from_pv` follows the conventions of
[Vallado (2013)](references.md#vallado2013), Algorithm 9: for a circular
inclined orbit `w` is 0 and `nu` holds the argument of latitude; for an
elliptical equatorial orbit `raan` is 0 and `w` holds the true longitude of
perigee; for a circular equatorial orbit both are 0 and `nu` holds the true
longitude. In each case `to_pv` reproduces the input state.

## Conversions

- **Elements → state** follows [Vallado (2013)]references.md#vallado2013,
  Algorithm 10: the state is formed in the perifocal (PQW) frame and rotated
  by $R_z(\Omega)\,R_x(i)\,R_z(\omega)$.
- **State → elements** follows Algorithm 9, with every angle extracted by
  `atan2` rather than `acos` so that near-zero inclinations (down to
  $10^{-9}$ rad) and anomalies near perigee/apogee are recovered to full
  precision; the round trip state → elements → state is accurate to better
  than $10^{-6}$ relative for $e \le 0.999$.
- **`propagate(dt)`** is pure two-body motion: only the mean anomaly advances,
  by $n\,\Delta t$. No perturbation is applied. For anything beyond a quick
  look, use the numerical propagator ([Force Model]forces.md).

## Examples

=== "Python"

    ```python
    import math
    import numpy as np
    import satkit as sk

    # Sun-synchronous-ish LEO, located by mean anomaly
    k = sk.kepler(
        a=7000.0e3,
        eccen=0.001,
        incl=math.radians(98.0),
        raan=math.radians(45.0),
        w=0.0,
        mean_anomaly=math.radians(30.0),
    )
    print(f"period = {k.period / 60:.2f} min, nu = {math.degrees(k.nu):.3f} deg")

    # Elements -> GCRF state -> elements
    r, v = k.to_pv()
    k2 = sk.kepler.from_pv(r, v)
    assert abs(k2.a - k.a) < 1e-3

    # Two-body propagation by a quarter period
    k3 = k.propagate(k.period / 4)
    print(f"mean anomaly after T/4 = {math.degrees(k3.mean_anomaly):.3f} deg")

    # An osculating snapshot of a numerically propagated state
    t0 = sk.time(2024, 1, 1)
    state = sk.satstate(t0, r, v)
    state1 = state.propagate(t0 + sk.duration.from_hours(1))
    k_osc = sk.kepler.from_pv(state1.pos, state1.vel)
    print(f"osculating a after 1 h: {k_osc.a / 1e3:.3f} km")
    ```

=== "Rust"

    ```rust
    use satkit::kepler::{Anomaly, Kepler};
    use satkit::Duration;

    let k = Kepler::new(
        7000.0e3,                 // a, m
        0.001,                    // eccen
        98.0_f64.to_radians(),    // incl
        45.0_f64.to_radians(),    // raan
        0.0,                      // w
        Anomaly::Mean(30.0_f64.to_radians()),
    );
    println!("period = {:.2} min, nu = {:.3} deg",
             k.period() / 60.0, k.nu.to_degrees());

    // Elements -> state -> elements
    let (r, v) = k.to_pv();
    let k2 = Kepler::from_pv(r, v)?;
    assert!((k2.a - k.a).abs() < 1e-3);

    // Two-body propagation by a quarter period
    let k3 = k.propagate(&Duration::from_seconds(k.period() / 4.0));
    println!("M after T/4 = {:.3} deg", k3.mean_anomaly().to_degrees());
    # Ok::<(), satkit::kepler::Error>(())
    ```

## References

- [Vallado, D. A. (2013)]references.md#vallado2013, *Fundamentals of Astrodynamics and Applications*, 4th ed., Microcosm Press. Algorithm 2 (Kepler's equation), Algorithm 9 (RV2COE), Algorithm 10 (COE2RV); §2.2–2.5.
- [Danby, J. M. A. (1987)]references.md#danby1987, "The solution of Kepler's equation, III," *Celestial Mechanics*, 40, 303–312. <https://doi.org/10.1007/BF01235847>