pub fn mean_motion(mu: Mu, a: Km) -> f64
Compute mean anomaly at a future time given initial mean anomaly, mean motion n, and elapsed time Δt.
M(t) = M₀ + n * Δt
Mean motion n = sqrt(μ/a³)