diff --git a/CHANGELOG.md b/CHANGELOG.md
index 511192e8..abb75765 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,12 +4,17 @@ Only recent releases are listed. Older entries are in this file's git history (`
 
 ## Unreleased
 
+### Added
+
+- Kepler: checked constructor `Kepler::try_new` / `Kepler::validate` (Python: the constructor and the `a`/`eccen`/`inclination`/`mu` setters raise `ValueError` for a non-finite value, `a <= 0`, `eccen` outside [0, 1), `incl` outside [0, π] or `mu <= 0` instead of producing NaN); per-instance gravitational parameter `mu` (`Kepler::with_mu`, `from_pv_with_mu`; Python `kepler(..., mu=)`, `kepler.mu`, `from_pv(..., mu=)`) so lunar and heliocentric elements have a correct period, `propagate` and `to_pv`; derived quantities `periapsis`, `apoapsis`, `specific_energy`, `angular_momentum`, `flight_path_angle`, `argument_of_latitude`, `true_longitude`; `SatState::from_kepler` / `satstate.from_kepler(time, kepler)`; a one-line `repr` for Python `kepler`; `Kepler` derives `PartialEq` and serde (`mu` defaults to Earth's when absent), `Anomaly` derives `Copy`/`PartialEq` ([#168](https://github.com/ssmichael1/satkit/pull/168))
+
 ### Fixed
 
 - A corrupt or truncated `tab5.2*.txt` in a data directory no longer panics the first frame transform: satkit warns and uses the compiled-in copy of the same IERS table (exact, not an approximation), and the parser now rejects text with no table header or fewer rows than declared — an HTML notice page saved under the table's name previously loaded as six empty series and silently dropped the nutation terms ([#166](https://github.com/ssmichael1/satkit/pull/166))
 
 ### Changed
 
+- **Breaking (Rust):** `Kepler.w` is renamed `Kepler.argp`, the struct gains a `mu` field (struct-literal construction must supply it — prefer `Kepler::new(...).with_mu(...)`), and `kepler::Error` is `#[non_exhaustive]` (new `InvalidElement` variant). Python: `argp` is the constructor parameter and property name; `w` still works as a constructor keyword and as a property (kept indefinitely) but emits `DeprecationWarning`. **Breaking (Python):** `kepler.from_pv` raises `ValueError` instead of `RuntimeError` for a hyperbolic/parabolic or rectilinear state, and every element setter (`a`, `eccen`, `inclination`, `raan`, `argp`, `nu`, `mu`) raises `ValueError` for an out-of-domain or non-finite value instead of accepting it ([#168](https://github.com/ssmichael1/satkit/pull/168))
 - `update_datafiles()` no longer downloads files that are compiled into the library: the IERS tables and gravity models are `default: false` in the manifest (still pinned and fetchable by name), and the unused `leap-seconds.list` (nothing ever read it — the runtime leap-second table is a compiled-in constant) is removed from the manifest entirely. The only static download left is the JPL ephemeris, alongside the daily EOP / space-weather / solar-cycle refreshes ([#163](https://github.com/ssmichael1/satkit/pull/163))
 
 ### Docs
diff --git a/docs/guide/kepler.md b/docs/guide/kepler.md
index dca3a8ae..bd5281e4 100644
--- a/docs/guide/kepler.md
+++ b/docs/guide/kepler.md
@@ -8,16 +8,18 @@ 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:
+these names, in SI units and radians, together with the gravitational
+parameter of the body they orbit:
 
 | Field   | Symbol     | Meaning                                              |
 |---------|------------|------------------------------------------------------|
-| `a`     | $a$        | semi-major axis, **meters**                          |
+| `a`     | $a$        | semi-major axis, **meters**, $a > 0$                 |
 | `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                         |
+| `argp`  | $\omega$   | argument of periapsis, radians                       |
 | `nu`    | $\nu$      | true anomaly, radians                                |
+| `mu`    | $\mu$      | gravitational parameter of the central body, m³ s⁻²; Earth's unless given |
 
 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
@@ -26,9 +28,45 @@ 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
+argument and the Rust field are `incl`. The argument of periapsis was called
+`w` before 0.22; in Python `w` still works as a constructor keyword and as a
+property (kept indefinitely), with a `DeprecationWarning`. Angles returned by `from_pv` are
 reduced to $[0, 2\pi)$; angles you set are stored as given.
 
+### Validation
+
+The Python constructor, `from_pv`, and every element setter reject an
+element outside its domain with `ValueError`: a non-finite value,
+$a \le 0$, $e \notin [0, 1)$, $i \notin [0, \pi]$ or $\mu \le 0$. The bounds
+are strict — $e = 1$ is not a closed orbit and is refused rather than
+producing NaN anomalies. The `mean_anomaly` / `eccentric_anomaly` setters
+likewise refuse a non-finite value, so an element set can never hold NaN
+through any setter. In Rust, `Kepler::try_new`
+performs the same checks (returning `kepler::Error::InvalidElement`, which
+names the offending element) and `Kepler::validate` re-checks an element set
+whose public fields were assigned directly; `Kepler::new` remains unchecked.
+
+### Derived quantities
+
+All are read-only properties in Python and methods in Rust:
+
+| Property                | Symbol / formula                                   | Units |
+|-------------------------|----------------------------------------------------|-------|
+| `semiparameter`         | $p = a(1 - e^2)$                                   | m     |
+| `periapsis`             | $r_p = a(1 - e)$                                   | m     |
+| `apoapsis`              | $r_a = a(1 + e)$                                   | m     |
+| `mean_motion`           | $n = \sqrt{\mu / a^3}$                             | rad/s |
+| `period`                | $T = 2\pi / n$                                     | s     |
+| `specific_energy`       | $\xi = -\mu / 2a$                                  | J/kg  |
+| `angular_momentum`      | $h = \sqrt{\mu p}$                                 | m²/s  |
+| `flight_path_angle`     | $\gamma = \operatorname{atan2}(e\sin\nu,\ 1 + e\cos\nu)$ — zero at periapsis and apoapsis, positive while climbing | rad |
+| `argument_of_latitude`  | $u = \omega + \nu$, reduced to $[0, 2\pi)$ — defined for circular orbits | rad |
+| `true_longitude`        | $\lambda = \Omega + \omega + \nu$, reduced to $[0, 2\pi)$ — defined for circular equatorial orbits | rad |
+
+`satstate.from_kepler(time, k)` (Rust `SatState::from_kepler`) builds a
+propagatable state from the two-body position and velocity of an element
+set, taken as GCRF at `time`.
+
 ### Anomalies
 
 Three angles can locate the satellite in its orbit
@@ -75,26 +113,31 @@ the rest of `satkit` assumes **GCRF**, so convert ITRF or TEME states with
 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`)
+**Central body.** Each element set carries its own $\mu$. By default it is
+the Earth's, [`consts.MU_EARTH`](../api/consts.md) ($3.986004418 \times
+10^{14}$ m³ s⁻²), and it is used everywhere: in the `from_pv` energy
+equation, in `to_pv`, and in the mean motion, period and `propagate`. For a
+lunar or heliocentric orbit pass `mu=` to the constructor or to `from_pv`
+(Rust: `Kepler::with_mu`, `Kepler::from_pv_with_mu`); the six geometric
+elements are unchanged, only the dynamics re-target the other body. Note
+that `from_pv` interprets a state with whatever $\mu$ it is given — a lunar
+state read with Earth's $\mu$ yields a valid-looking but wrong ellipse.
+
+**Closed orbits only.** `from_pv` returns an error (Python: `ValueError`)
 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.
+(zero angular momentum, where the orbital plane is undefined). The Python
+constructor and `Kepler::try_new` reject $e \ge 1$ up front (see
+[Validation](#validation)).
 
 **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
+inclined orbit `argp` is 0 and `nu` holds the argument of latitude; for an
+elliptical equatorial orbit `raan` is 0 and `argp` 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.
+longitude. In each case `to_pv` reproduces the input state. The
+`argument_of_latitude` and `true_longitude` properties give the well-defined
+combinations directly in every case.
 
 ## Conversions
 
@@ -125,7 +168,7 @@ longitude. In each case `to_pv` reproduces the input state.
         eccen=0.001,
         incl=math.radians(98.0),
         raan=math.radians(45.0),
-        w=0.0,
+        argp=0.0,
         mean_anomaly=math.radians(30.0),
     )
     print(f"period = {k.period / 60:.2f} min, nu = {math.degrees(k.nu):.3f} deg")
@@ -158,7 +201,7 @@ longitude. In each case `to_pv` reproduces the input state.
         0.001,                    // eccen
         98.0_f64.to_radians(),    // incl
         45.0_f64.to_radians(),    // raan
-        0.0,                      // w
+        0.0,                      // argp
         Anomaly::Mean(30.0_f64.to_radians()),
     );
     println!("period = {:.2} min, nu = {:.3} deg",
diff --git a/python/satkit/satkit.pyi b/python/satkit/satkit.pyi
index f95f5a4f..490f0e16 100644
--- a/python/satkit/satkit.pyi
+++ b/python/satkit/satkit.pyi
@@ -2321,8 +2321,14 @@ class kepler:
         - The class uses the semi-major axis (a), not the semiparameter
         - Elements are osculating and expressed in the frame of the input state
           (normally GCRF); the class does no frame handling
-        - The Earth gravitational parameter (MU_EARTH) is used throughout
-        - Only closed orbits are supported (0 <= eccen < 1)
+        - Each element set carries the gravitational parameter of its central
+          body, ``mu`` (m^3/s^2); Earth's (``satkit.consts.mu_earth``) unless
+          given, so lunar or heliocentric elements are supported by passing
+          ``mu=satkit.consts.mu_moon`` / ``mu_sun``
+        - Only closed orbits are supported (0 <= eccen < 1); the constructor,
+          ``from_pv`` and every element setter — the anomaly setters
+          included — raise ``ValueError`` for an element outside its domain
+          or a non-finite value, so an element set can never hold NaN
         - All angle units are radians
         - All length units are meters
         - All velocity units are meters / second
@@ -2336,25 +2342,38 @@ class kepler:
         eccen: float,
         incl: float,
         raan: float,
-        w: float,
+        argp: float | None = None,
         nu: float | None = None,
         *,
+        w: float | None = None,
         true_anomaly: float | None = None,
         eccentric_anomaly: float | None = None,
         mean_anomaly: float | None = None,
+        mu: float | None = None,
     ) -> None:
         """Create Keplerian element set object from input elements
 
         Args:
-            a: Semi-major axis, meters
+            a: Semi-major axis, meters (> 0)
             eccen: Eccentricity, unitless (0 <= eccen < 1)
-            incl: Inclination, radians
+            incl: Inclination, radians (0 <= incl <= pi)
             raan: Right ascension of ascending node, radians
-            w: Argument of perigee, radians
+            argp: Argument of periapsis, radians (5th positional argument)
             nu: True anomaly, radians (6th positional argument)
+            w: Argument of periapsis, radians — deprecated keyword alias of
+                ``argp`` (kept indefinitely); give one or the other, not both
             true_anomaly: True anomaly, radians (keyword alternative to nu)
             eccentric_anomaly: Eccentric anomaly, radians (keyword alternative to nu)
             mean_anomaly: Mean anomaly, radians (keyword alternative to nu)
+            mu: Gravitational parameter of the central body, m^3/s^2
+                (default ``satkit.consts.mu_earth``)
+
+        Raises:
+            ValueError: an element outside its domain (non-finite value,
+                ``a <= 0``, ``eccen`` outside [0, 1), ``incl`` outside
+                [0, pi], ``mu <= 0``); more or fewer than one anomaly given;
+                both ``argp`` and ``w`` given
+            TypeError: no argument of periapsis given
 
         Notes:
             Exactly one of ``nu``, ``true_anomaly``, ``eccentric_anomaly`` or
@@ -2372,12 +2391,17 @@ class kepler:
                 eccen=0.001,  # near-circular
                 incl=math.radians(51.6),
                 raan=math.radians(0),
-                w=math.radians(0),
+                argp=math.radians(0),
                 nu=math.radians(0),
             )
 
             # Same orbit, positional, located by mean anomaly instead
             k2 = satkit.kepler(6.781e6, 0.001, math.radians(51.6), 0, 0, mean_anomaly=1.0)
+
+            # A 100 km circular lunar orbit
+            k_moon = satkit.kepler(1837.4e3, 0.0, math.radians(90), 0, 0, 0,
+                                   mu=satkit.consts.mu_moon)
+            print(f"lunar period: {k_moon.period / 60:.1f} min")
             ```
         """
         ...
@@ -2438,7 +2462,11 @@ class kepler:
 
     @eccentric_anomaly.setter
     def eccentric_anomaly(self, value: float) -> None:
-        """Set the in-plane position by eccentric anomaly, radians"""
+        """Set the in-plane position by eccentric anomaly, radians
+
+        Converted to true anomaly (``nu``) on the spot. A non-finite value
+        raises ``ValueError`` and leaves the element set unchanged.
+        """
         ...
     @property
     def mean_anomaly(self) -> float:
@@ -2450,7 +2478,8 @@ class kepler:
         """Set the in-plane position by mean anomaly, radians
 
         Kepler's equation is solved for the eccentric anomaly and the result
-        stored as true anomaly (``nu``). A non-finite value yields NaN.
+        stored as true anomaly (``nu``). A non-finite value raises
+        ``ValueError`` and leaves the element set unchanged.
         """
         ...
     @property
@@ -2463,55 +2492,135 @@ class kepler:
         """Semiparameter (semi-latus rectum) p = a (1 - e^2), meters"""
         ...
 
+    @property
+    def periapsis(self) -> float:
+        """Radius of periapsis a (1 - e), meters"""
+        ...
+
+    @property
+    def apoapsis(self) -> float:
+        """Radius of apoapsis a (1 + e), meters"""
+        ...
+
+    @property
+    def specific_energy(self) -> float:
+        """Specific orbital energy -mu / (2 a), J/kg (m^2/s^2)"""
+        ...
+
+    @property
+    def angular_momentum(self) -> float:
+        """Magnitude of the specific angular momentum sqrt(mu p), m^2/s"""
+        ...
+
+    @property
+    def flight_path_angle(self) -> float:
+        """Flight-path angle, radians
+
+        atan2(e sin nu, 1 + e cos nu): the angle of the velocity above the
+        local horizontal — zero at periapsis and apoapsis, positive while
+        climbing.
+        """
+        ...
+
+    @property
+    def argument_of_latitude(self) -> float:
+        """Argument of latitude u = argp + nu, radians in [0, 2 pi)
+
+        Well defined for circular orbits, where ``argp`` and ``nu`` separately
+        are not.
+        """
+        ...
+
+    @property
+    def true_longitude(self) -> float:
+        """True longitude raan + argp + nu, radians in [0, 2 pi)
+
+        Well defined for circular equatorial orbits, where ``raan``, ``argp``
+        and ``nu`` separately are not.
+        """
+        ...
+
+    @property
+    def mu(self) -> float:
+        """Gravitational parameter of the central body, m^3/s^2
+
+        Setting it re-targets the dynamics (period, mean motion,
+        ``propagate``, ``to_pv``) at another body while keeping the six
+        geometric elements. Must be positive and finite (``ValueError``).
+        """
+        ...
+
+    @mu.setter
+    def mu(self, value: float) -> None: ...
     @property
     def a(self) -> float:
-        """Semi-major axis, meters"""
+        """Semi-major axis, meters (> 0; ``ValueError`` otherwise)"""
         ...
 
     @a.setter
     def a(self, value: float) -> None: ...
     @property
     def eccen(self) -> float:
-        """Eccentricity, unitless"""
+        """Eccentricity, unitless (0 <= eccen < 1; ``ValueError`` otherwise)"""
         ...
 
     @eccen.setter
     def eccen(self, value: float) -> None: ...
     @property
     def inclination(self) -> float:
-        """Inclination, radians"""
+        """Inclination, radians (0 <= incl <= pi; ``ValueError`` otherwise)"""
         ...
 
     @inclination.setter
     def inclination(self, value: float) -> None: ...
     @property
     def raan(self) -> float:
-        """Right ascension of ascending node, radians"""
+        """Right ascension of ascending node, radians (finite; ``ValueError`` otherwise)"""
         ...
 
     @raan.setter
     def raan(self, value: float) -> None: ...
     @property
     def nu(self) -> float:
-        """True anomaly, radians"""
+        """True anomaly, radians (finite; ``ValueError`` otherwise)"""
         ...
 
     @nu.setter
     def nu(self, value: float) -> None: ...
     @property
+    def argp(self) -> float:
+        """Argument of periapsis, radians (finite; ``ValueError`` otherwise)"""
+        ...
+
+    @argp.setter
+    def argp(self, value: float) -> None: ...
+    @property
     def w(self) -> float:
-        """Argument of perigee, radians"""
+        """Argument of periapsis, radians — deprecated alias of ``argp``
+
+        Kept indefinitely for compatibility; reading or assigning it emits
+        ``DeprecationWarning``. Validated like ``argp``.
+        """
         ...
 
     @w.setter
     def w(self, value: float) -> None: ...
     @staticmethod
-    def from_pv(pos: npt.NDArray[np.float64], vel: npt.NDArray[np.float64]) -> kepler:
+    def from_pv(
+        pos: npt.NDArray[np.float64],
+        vel: npt.NDArray[np.float64],
+        *,
+        mu: float | None = None,
+    ) -> kepler:
         """Create Keplerian element set from input position and velocity vectors
 
         Args:
-            pos: 3-element array representing position vector
-            vel: 3-element array representing velocity vector
+            pos: 3-element position vector, meters
+            vel: 3-element velocity vector, meters/second
+            mu: Gravitational parameter of the central body, m^3/s^2
+                (default ``satkit.consts.mu_earth``); the returned element
+                set carries it, so its period and ``to_pv`` refer to the same
+                body
 
         Returns:
             Keplerian element set object
@@ -2527,9 +2636,10 @@ class kepler:
             ```
 
         Raises:
-            RuntimeError: if the state is hyperbolic/parabolic (eccen >= 1) or
-                rectilinear (zero angular momentum), or if the inputs are not
-                3-element vectors.
+            ValueError: if the state is hyperbolic/parabolic (eccen >= 1) or
+                rectilinear (zero angular momentum), or if ``mu`` is not
+                positive and finite.
+            RuntimeError: if the inputs are not 3-element vectors.
         """
         ...
 
@@ -2911,6 +3021,30 @@ class satstate:
         """
         ...
 
+    @staticmethod
+    def from_kepler(time: _Time, kepler: kepler) -> satstate:
+        """Create a state from Keplerian elements
+
+        The two-body position (meters) and velocity (m/s) of the elements
+        (``kepler.to_pv()``) are taken as GCRF at ``time``; no covariance, no
+        maneuvers. The elements are treated as osculating GCRF elements, so
+        ``kepler.mu`` should be Earth's.
+
+        Args:
+            time (satkit.time): Epoch of the state
+            kepler (satkit.kepler): Osculating Keplerian elements
+
+        Returns:
+            satstate: state at ``time``
+
+        Example:
+            ```python
+            k = satkit.kepler(7000e3, 0.001, math.radians(98), 0, 0, 0)
+            state = satkit.satstate.from_kepler(satkit.time(2024, 1, 1), k)
+            ```
+        """
+        ...
+
     @property
     def pos(self) -> npt.NDArray[np.float64]:
         """Position in meters, GCRF frame (alias for pos_gcrf)"""
diff --git a/python/src/pykepler.rs b/python/src/pykepler.rs
index d1bf567d..f6acff98 100644
--- a/python/src/pykepler.rs
+++ b/python/src/pykepler.rs
@@ -17,17 +17,25 @@ use crate::pyutils::py_to_smatrix;
 /// ``eccentric_anomaly`` or ``mean_anomaly``.
 ///
 /// Args:
-///     a: semi-major axis, meters
-///     eccen: Eccentricity
-///     incl: Inclination, radians
+///     a: semi-major axis, meters (> 0)
+///     eccen: Eccentricity (0 <= eccen < 1)
+///     incl: Inclination, radians (0 <= incl <= pi)
 ///     raan: Right Ascension of the Ascending Node, radians
-///     w: Argument of Perigee, radians
+///     argp: Argument of Periapsis, radians (``w`` is accepted as an alias)
 ///     nu: True Anomaly, radians
 ///
 /// Keyword Args:
 ///     true_anomaly: True Anomaly, radians
 ///     eccentric_anomaly: Eccentric Anomaly, radians
 ///     mean_anomaly: Mean Anomaly, radians
+///     mu: Gravitational parameter of the central body, m^3/s^2
+///         (default: Earth, ``satkit.consts.mu_earth``)
+///
+/// Raises:
+///     ValueError: an element outside its domain (non-finite, ``a <= 0``,
+///         ``eccen`` outside [0, 1), ``incl`` outside [0, pi], ``mu <= 0``),
+///         or an ambiguous anomaly / ``argp``-``w`` specification. The
+///         element setters apply the same checks.
 ///
 /// Returns:
 ///     Kepler: Keplerian orbital elements
@@ -36,31 +44,83 @@ use crate::pyutils::py_to_smatrix;
 #[derive(Clone)]
 pub struct PyKepler(pub Kepler);
 
+/// Map a core validation failure to the Python exception the stubs promise.
+fn value_error(e: satkit::kepler::Error) -> PyErr {
+    pyo3::exceptions::PyValueError::new_err(e.to_string())
+}
+
 impl PyKepler {
     /// The same elements with the in-plane position replaced by `an`,
     /// converted to true anomaly by the core solver.
     fn with_anomaly(&self, an: Anomaly) -> Kepler {
         let k = &self.0;
-        Kepler::new(k.a, k.eccen, k.incl, k.raan, k.w, an)
+        Kepler::new(k.a, k.eccen, k.incl, k.raan, k.argp, an).with_mu(k.mu)
+    }
+
+    /// Assign `field` on a copy, validate, and commit only if the result is
+    /// a closed orbit; otherwise `ValueError` and the element set is unchanged.
+    fn set_validated(&mut self, field: impl FnOnce(&mut Kepler)) -> PyResult<()> {
+        let mut k = self.0;
+        field(&mut k);
+        k.validate().map_err(value_error)?;
+        self.0 = k;
+        Ok(())
+    }
+
+    /// `ValueError` for a non-finite anomaly, matching the message format of
+    /// the core validation errors.
+    fn check_anomaly(name: &str, val: f64) -> PyResult<()> {
+        if val.is_finite() {
+            Ok(())
+        } else {
+            Err(pyo3::exceptions::PyValueError::new_err(format!(
+                "invalid Keplerian element {name} = {val}: must be finite"
+            )))
+        }
+    }
+
+    fn warn_w_deprecated(py: Python) -> PyResult<()> {
+        let warning_type = py.get_type::<pyo3::exceptions::PyDeprecationWarning>();
+        PyErr::warn(
+            py,
+            warning_type.as_any(),
+            c"kepler.w is deprecated; use kepler.argp",
+            2,
+        )
     }
 }
 
 #[pymethods]
 impl PyKepler {
     #[new]
-    #[pyo3(signature = (a, eccen, incl, raan, w, nu=None, *, true_anomaly=None, eccentric_anomaly=None, mean_anomaly=None))]
+    #[pyo3(signature = (a, eccen, incl, raan, argp=None, nu=None, *, w=None, true_anomaly=None, eccentric_anomaly=None, mean_anomaly=None, mu=None))]
     #[allow(clippy::too_many_arguments)]
     fn new(
         a: f64,
         eccen: f64,
         incl: f64,
         raan: f64,
-        w: f64,
+        argp: Option<f64>,
         nu: Option<f64>,
+        w: Option<f64>,
         true_anomaly: Option<f64>,
         eccentric_anomaly: Option<f64>,
         mean_anomaly: Option<f64>,
+        mu: Option<f64>,
     ) -> PyResult<Self> {
+        let argp = match (argp, w) {
+            (Some(v), None) | (None, Some(v)) => v,
+            (None, None) => {
+                return Err(pyo3::exceptions::PyTypeError::new_err(
+                    "missing required argument: 'argp' (argument of periapsis, radians)",
+                ))
+            }
+            (Some(_), Some(_)) => {
+                return Err(pyo3::exceptions::PyValueError::new_err(
+                    "specify the argument of periapsis as argp or w, not both",
+                ))
+            }
+        };
         let an =
             match (nu, true_anomaly, eccentric_anomaly, mean_anomaly) {
                 (Some(v), None, None, None) | (None, Some(v), None, None) => Anomaly::True(v),
@@ -70,7 +130,12 @@ impl PyKepler {
                     "Specify exactly one of nu, true_anomaly, eccentric_anomaly, or mean_anomaly",
                 )),
             };
-        Ok(Self(Kepler::new(a, eccen, incl, raan, w, an)))
+        let mut k = Kepler::try_new(a, eccen, incl, raan, argp, an).map_err(value_error)?;
+        if let Some(mu) = mu {
+            k = k.with_mu(mu);
+            k.validate().map_err(value_error)?;
+        }
+        Ok(Self(k))
     }
 
     #[getter]
@@ -80,19 +145,19 @@ impl PyKepler {
     }
 
     #[setter(a)]
-    fn set_a(&mut self, val: f64) {
-        self.0.a = val;
+    fn set_a(&mut self, val: f64) -> PyResult<()> {
+        self.set_validated(|k| k.a = val)
     }
 
     #[getter]
-    /// Eccentricity
+    /// Eccentricity, unitless
     fn get_eccen(&self) -> f64 {
         self.0.eccen
     }
 
     #[setter(eccen)]
-    fn set_eccen(&mut self, val: f64) {
-        self.0.eccen = val;
+    fn set_eccen(&mut self, val: f64) -> PyResult<()> {
+        self.set_validated(|k| k.eccen = val)
     }
 
     #[getter]
@@ -102,8 +167,8 @@ impl PyKepler {
     }
 
     #[setter(inclination)]
-    fn set_inclination(&mut self, val: f64) {
-        self.0.incl = val;
+    fn set_inclination(&mut self, val: f64) -> PyResult<()> {
+        self.set_validated(|k| k.incl = val)
     }
 
     #[getter]
@@ -113,19 +178,35 @@ impl PyKepler {
     }
 
     #[setter(raan)]
-    fn set_raan(&mut self, val: f64) {
-        self.0.raan = val;
+    fn set_raan(&mut self, val: f64) -> PyResult<()> {
+        self.set_validated(|k| k.raan = val)
     }
 
     #[getter]
-    /// Argument of Perigee, radians
-    fn get_w(&self) -> f64 {
-        self.0.w
+    /// Argument of Periapsis, radians
+    fn get_argp(&self) -> f64 {
+        self.0.argp
+    }
+
+    #[setter(argp)]
+    fn set_argp(&mut self, val: f64) -> PyResult<()> {
+        self.set_validated(|k| k.argp = val)
+    }
+
+    #[getter]
+    /// Argument of Periapsis, radians
+    ///
+    /// Deprecated alias of ``argp`` (kept indefinitely); emits
+    /// ``DeprecationWarning``.
+    fn get_w(&self, py: Python) -> PyResult<f64> {
+        Self::warn_w_deprecated(py)?;
+        Ok(self.0.argp)
     }
 
     #[setter(w)]
-    fn set_w(&mut self, val: f64) {
-        self.0.w = val;
+    fn set_w(&mut self, py: Python, val: f64) -> PyResult<()> {
+        Self::warn_w_deprecated(py)?;
+        self.set_validated(|k| k.argp = val)
     }
 
     #[getter]
@@ -135,8 +216,19 @@ impl PyKepler {
     }
 
     #[setter(nu)]
-    fn set_nu(&mut self, val: f64) {
-        self.0.nu = val;
+    fn set_nu(&mut self, val: f64) -> PyResult<()> {
+        self.set_validated(|k| k.nu = val)
+    }
+
+    #[getter]
+    /// Gravitational parameter of the central body, m^3/s^2
+    fn get_mu(&self) -> f64 {
+        self.0.mu
+    }
+
+    #[setter(mu)]
+    fn set_mu(&mut self, val: f64) -> PyResult<()> {
+        self.set_validated(|k| k.mu = val)
     }
 
     /// Convert Keplerian elements to Cartesian
@@ -152,14 +244,33 @@ impl PyKepler {
     }
 
     /// Convert Cartesian elements to kepler
+    ///
+    /// Args:
+    ///     pos: 3-element position vector, meters
+    ///     vel: 3-element velocity vector, meters/second
+    ///
+    /// Keyword Args:
+    ///     mu: Gravitational parameter of the central body, m^3/s^2
+    ///         (default: Earth); the returned elements carry it
+    ///
+    /// Raises:
+    ///     ValueError: open (eccen >= 1) or rectilinear (zero angular
+    ///         momentum) state, or ``mu`` not positive and finite
+    ///     RuntimeError: inputs that are not 3-element vectors
     #[staticmethod]
-    fn from_pv(pos: &Bound<PyAny>, vel: &Bound<PyAny>) -> PyResult<Self> {
+    #[pyo3(signature = (pos, vel, *, mu=None))]
+    fn from_pv(pos: &Bound<PyAny>, vel: &Bound<PyAny>, mu: Option<f64>) -> PyResult<Self> {
         let pos = py_to_smatrix(pos)?;
         let vel = py_to_smatrix(vel)?;
-        match Kepler::from_pv(pos, vel) {
-            Ok(k) => Ok(Self(k)),
-            Err(e) => Err(pyo3::exceptions::PyRuntimeError::new_err(e.to_string())),
+        let mu = mu.unwrap_or(satkit::consts::MU_EARTH);
+        if !(mu.is_finite() && mu > 0.0) {
+            return Err(pyo3::exceptions::PyValueError::new_err(format!(
+                "invalid Keplerian element mu = {mu}: gravitational parameter must be positive"
+            )));
         }
+        Kepler::from_pv_with_mu(pos, vel, mu)
+            .map(Self)
+            .map_err(value_error)
     }
 
     /// Propagate the elements forward (or backward) in time
@@ -195,8 +306,10 @@ impl PyKepler {
     }
 
     #[setter(eccentric_anomaly)]
-    fn set_eccentric_anomaly(&mut self, val: f64) {
+    fn set_eccentric_anomaly(&mut self, val: f64) -> PyResult<()> {
+        Self::check_anomaly("eccentric_anomaly", val)?;
         self.0 = self.with_anomaly(Anomaly::Eccentric(val));
+        Ok(())
     }
 
     /// Return the mean motion of the satellite in radians/second
@@ -226,6 +339,52 @@ impl PyKepler {
         self.0.semiparameter()
     }
 
+    /// Radius of periapsis a (1 - e), meters
+    #[getter]
+    fn periapsis(&self) -> f64 {
+        self.0.periapsis()
+    }
+
+    /// Radius of apoapsis a (1 + e), meters
+    #[getter]
+    fn apoapsis(&self) -> f64 {
+        self.0.apoapsis()
+    }
+
+    /// Specific orbital energy -mu / 2a, J/kg
+    #[getter]
+    fn specific_energy(&self) -> f64 {
+        self.0.specific_energy()
+    }
+
+    /// Magnitude of the specific angular momentum sqrt(mu p), m^2/s
+    #[getter]
+    fn angular_momentum(&self) -> f64 {
+        self.0.angular_momentum()
+    }
+
+    /// Flight-path angle atan2(e sin nu, 1 + e cos nu), radians: the angle of
+    /// the velocity above the local horizontal, zero at periapsis and
+    /// apoapsis, positive while climbing
+    #[getter]
+    fn flight_path_angle(&self) -> f64 {
+        self.0.flight_path_angle()
+    }
+
+    /// Argument of latitude u = argp + nu, radians in [0, 2 pi); well defined
+    /// for circular orbits
+    #[getter]
+    fn argument_of_latitude(&self) -> f64 {
+        self.0.argument_of_latitude()
+    }
+
+    /// True longitude raan + argp + nu, radians in [0, 2 pi); well defined
+    /// for circular equatorial orbits
+    #[getter]
+    fn true_longitude(&self) -> f64 {
+        self.0.true_longitude()
+    }
+
     /// Return the mean anomaly of the satellite in radians
     ///
     /// Returns:
@@ -236,10 +395,13 @@ impl PyKepler {
     }
 
     #[setter(mean_anomaly)]
-    fn set_mean_anomaly(&mut self, val: f64) {
-        // Kepler's equation is solved by the core crate (range-reduced,
-        // Danby start, iteration-capped), so NaN or e >= 1 cannot hang here.
+    fn set_mean_anomaly(&mut self, val: f64) -> PyResult<()> {
+        // A non-finite M is refused up front so the element set can never
+        // hold a NaN anomaly; the core solver itself is iteration-capped and
+        // cannot hang on any input regardless.
+        Self::check_anomaly("mean_anomaly", val)?;
         self.0 = self.with_anomaly(Anomaly::Mean(val));
+        Ok(())
     }
 
     /// Return the true anomaly of the satellite in radians
@@ -256,18 +418,15 @@ impl PyKepler {
     }
 
     fn __repr__(&self) -> String {
-        self.__str__()
+        let k = &self.0;
+        format!(
+            "kepler(a={:.6e}, eccen={:.6e}, incl={:.6e}, raan={:.6e}, argp={:.6e}, nu={:.6e}, mu={:.6e})",
+            k.a, k.eccen, k.incl, k.raan, k.argp, k.nu, k.mu
+        )
     }
 
     fn __eq__(&self, other: &Self) -> bool {
-        let a = &self.0;
-        let b = &other.0;
-        a.a == b.a
-            && a.eccen == b.eccen
-            && a.incl == b.incl
-            && a.raan == b.raan
-            && a.w == b.w
-            && a.nu == b.nu
+        self.0 == other.0
     }
 
     fn __ne__(&self, other: &Self) -> bool {
@@ -282,20 +441,24 @@ impl PyKepler {
                 self.0.eccen,
                 self.0.incl,
                 self.0.raan,
-                self.0.w,
+                self.0.argp,
                 self.0.nu,
+                self.0.mu,
             ],
         )
     }
 
     fn __setstate__(&mut self, py: Python, state: Py<PyBytes>) -> PyResult<()> {
-        let [a, eccen, incl, raan, w, nu] = crate::pyutils::unpack_f64s(py, &state)?;
-        self.0.a = a;
-        self.0.eccen = eccen;
-        self.0.incl = incl;
-        self.0.raan = raan;
-        self.0.w = w;
-        self.0.nu = nu;
+        let [a, eccen, incl, raan, argp, nu, mu] = crate::pyutils::unpack_f64s(py, &state)?;
+        self.0 = Kepler {
+            a,
+            eccen,
+            incl,
+            raan,
+            argp,
+            nu,
+            mu,
+        };
         Ok(())
     }
 
diff --git a/python/src/pysatstate.rs b/python/src/pysatstate.rs
index bed6b855..7304e19b 100644
--- a/python/src/pysatstate.rs
+++ b/python/src/pysatstate.rs
@@ -80,6 +80,24 @@ impl PySatState {
         Ok(Self(state))
     }
 
+    /// Create a state from Keplerian elements
+    ///
+    /// The two-body position and velocity of the elements (``kepler.to_pv()``)
+    /// are taken as GCRF at ``time``; no covariance, no maneuvers. The
+    /// elements are treated as osculating GCRF elements, so ``kepler.mu``
+    /// should be Earth's.
+    ///
+    /// Args:
+    ///     time (satkit.time): epoch of the state
+    ///     kepler (satkit.kepler): osculating Keplerian elements
+    ///
+    /// Returns:
+    ///     satkit.satstate: state at ``time``
+    #[staticmethod]
+    fn from_kepler(time: &PyInstant, kepler: &crate::pykepler::PyKepler) -> Self {
+        Self(SatState::from_kepler(&time.0, &kepler.0))
+    }
+
     /// Set 1-sigma position uncertainty in a satellite-local or inertial frame.
     ///
     /// The uncertainty is interpreted as a diagonal 3x3 covariance in the
diff --git a/python/test/test_coordinates.py b/python/test/test_coordinates.py
index b2ca4ab5..6de5f574 100644
--- a/python/test/test_coordinates.py
+++ b/python/test/test_coordinates.py
@@ -21,7 +21,7 @@ def test_kepler_from_pv(self):
         assert kep.eccen == pytest.approx(0.83285, 1.0e-5)
         assert kep.inclination * rad2deg == pytest.approx(87.87, 1.0e-3)
         assert kep.raan * rad2deg == pytest.approx(227.89, 1.0e-3)
-        assert kep.w * rad2deg == pytest.approx(53.38, 1.0e-3)
+        assert kep.argp * rad2deg == pytest.approx(53.38, 1.0e-3)
         assert kep.nu * rad2deg == pytest.approx(92.335, 1.0e-3)
 
     def test_kepler_to_pv(self):
@@ -68,30 +68,40 @@ def test_kepler_kwargs_construction(self):
         with pytest.raises(TypeError):
             sk.kepler(7000e3, 0.1, 0.5, 1.0, 0.3, 0.7, bogus=1.0)
 
-    def test_mean_anomaly_setter_nan_does_not_hang(self):
-        """Regression: the setter used to spin forever on NaN with the GIL held."""
+    def test_anomaly_setters_reject_non_finite_and_high_e_converges(self):
+        """Regression for the setter that used to spin forever on NaN.
+
+        Non-finite input is now refused before reaching the solver (the
+        solver itself is iteration-capped; see the Rust test
+        ``test_mean2eccentric_nan_returns``), so no setter can leave a NaN
+        element behind. A high-eccentricity solve still runs, capped, in a
+        watchdog thread.
+        """
         import threading
 
+        k = sk.kepler(7000e3, 0.5, 0.5, 1.0, 0.3, 0.0)
+        k0 = sk.kepler(7000e3, 0.5, 0.5, 1.0, 0.3, 0.0)
+        for attr in ("mean_anomaly", "eccentric_anomaly"):
+            for bad in (float("nan"), float("inf"), float("-inf")):
+                with pytest.raises(ValueError, match=attr):
+                    setattr(k, attr, bad)
+        assert k == k0
+        with pytest.raises(ValueError):
+            k.eccen = 1.5
+        assert k == k0
+
         results = {}
 
         def worker():
-            k = sk.kepler(7000e3, 0.5, 0.5, 1.0, 0.3, 0.0)
-            k.mean_anomaly = float("nan")
-            results["nan"] = k.nu
-            k.mean_anomaly = float("inf")
-            results["inf"] = k.nu
-            # e >= 1 is outside the supported domain but must still return
-            k.eccen = 1.5
-            k.mean_anomaly = 1.0
-            results["hyper"] = k.nu
+            k.eccen = 0.999
+            k.mean_anomaly = 6.0
+            results["high_e"] = k.nu
 
         t = threading.Thread(target=worker, daemon=True)
         t.start()
         t.join(timeout=10.0)
         assert not t.is_alive(), "mean_anomaly setter hung"
-        assert m.isnan(results["nan"])
-        assert m.isnan(results["inf"])
-        assert "hyper" in results
+        assert m.isfinite(results["high_e"])
 
     def test_mean_anomaly_setter_roundtrip(self):
         """Set M, read it back: |dM| < 1e-12 at e = 0.9 over a full revolution."""
@@ -116,17 +126,149 @@ def test_mean_anomaly_setter_roundtrip(self):
     def test_kepler_pickle_roundtrip(self):
         import pickle
 
-        k = sk.kepler(7000e3, 0.1, 0.5, 1.0, 0.3, 0.7)
+        k = sk.kepler(7000e3, 0.1, 0.5, 1.0, 0.3, 0.7, mu=sk.consts.mu_moon)
         k2 = pickle.loads(pickle.dumps(k))
         assert k2 == k
-        assert (k2.a, k2.eccen, k2.inclination, k2.raan, k2.w, k2.nu) == (
+        assert (k2.a, k2.eccen, k2.inclination, k2.raan, k2.argp, k2.nu, k2.mu) == (
             k.a,
             k.eccen,
             k.inclination,
             k.raan,
-            k.w,
+            k.argp,
             k.nu,
+            k.mu,
         )
+        # mu is part of the state: a different central body is a different set
+        k_earth = sk.kepler(7000e3, 0.1, 0.5, 1.0, 0.3, 0.7)
+        assert k_earth != k
+        assert k_earth.mu == sk.consts.mu_earth
+
+    def test_kepler_validation_raises_value_error(self):
+        good = dict(a=7000e3, eccen=0.1, incl=0.5, raan=1.0, argp=0.3, nu=0.7)
+        sk.kepler(**good)
+        for bad in (
+            dict(a=0.0),
+            dict(a=-7000e3),
+            dict(a=float("nan")),
+            dict(eccen=1.0),
+            dict(eccen=-1e-3),
+            dict(eccen=float("inf")),
+            dict(incl=-1e-9),
+            dict(incl=m.pi + 1e-9),
+            dict(raan=float("nan")),
+            dict(nu=float("nan")),
+            dict(mu=0.0),
+            dict(mu=-1.0),
+        ):
+            with pytest.raises(ValueError):
+                sk.kepler(**{**good, **bad})
+        # Boundaries are inclusive where the domain is closed
+        sk.kepler(**{**good, "eccen": 0.0, "incl": 0.0})
+        sk.kepler(**{**good, "incl": m.pi})
+        # Setters validate too, and leave the element set unchanged on failure
+        k = sk.kepler(**good)
+        for attr, val in (
+            ("a", -1.0),
+            ("a", float("nan")),
+            ("eccen", 1.0),
+            ("inclination", 4.0),
+            ("mu", 0.0),
+            ("raan", float("nan")),
+            ("raan", float("inf")),
+            ("argp", float("nan")),
+            ("nu", float("-inf")),
+        ):
+            with pytest.raises(ValueError):
+                setattr(k, attr, val)
+        with pytest.warns(DeprecationWarning), pytest.raises(ValueError):
+            k.w = float("nan")
+        assert k == sk.kepler(**good)
+        # Finite angles of any size are accepted (stored as given)
+        k.raan = -10.0
+        k.argp = 100.0
+        k.nu = 7.0
+        assert (k.raan, k.argp, k.nu) == (-10.0, 100.0, 7.0)
+
+    def test_kepler_from_pv_open_or_rectilinear_is_value_error(self):
+        r = np.array([7000e3, 0.0, 0.0])
+        # Escape velocity and beyond: hyperbolic
+        v_esc = m.sqrt(2 * sk.consts.mu_earth / 7000e3)
+        with pytest.raises(ValueError, match="[Ee]ccentricity"):
+            sk.kepler.from_pv(r, np.array([0.0, 1.5 * v_esc, 0.0]))
+        # Parallel r and v: zero angular momentum
+        with pytest.raises(ValueError, match="angular momentum"):
+            sk.kepler.from_pv(r, np.array([1000.0, 0.0, 0.0]))
+        # Malformed input is still a RuntimeError (shape, not physics)
+        with pytest.raises(RuntimeError):
+            sk.kepler.from_pv(np.zeros(4), np.zeros(3))
+
+    def test_kepler_argp_and_deprecated_w(self):
+        k_argp = sk.kepler(7000e3, 0.1, 0.5, 1.0, argp=0.3, nu=0.7)
+        k_w = sk.kepler(7000e3, 0.1, 0.5, 1.0, w=0.3, nu=0.7)
+        k_pos = sk.kepler(7000e3, 0.1, 0.5, 1.0, 0.3, 0.7)
+        assert k_argp == k_w == k_pos
+        assert k_argp.argp == 0.3
+        with pytest.raises(ValueError):
+            sk.kepler(7000e3, 0.1, 0.5, 1.0, argp=0.3, w=0.3, nu=0.7)
+        with pytest.raises(TypeError):
+            sk.kepler(7000e3, 0.1, 0.5, 1.0, nu=0.7)
+        with pytest.warns(DeprecationWarning, match="argp"):
+            assert k_argp.w == 0.3
+        with pytest.warns(DeprecationWarning, match="argp"):
+            k_argp.w = 0.4
+        assert k_argp.argp == 0.4
+
+    def test_kepler_mu(self):
+        k = sk.kepler(2000e3, 0.05, 1.0, 0.2, 0.3, 0.4)
+        k_moon = sk.kepler(2000e3, 0.05, 1.0, 0.2, 0.3, 0.4, mu=sk.consts.mu_moon)
+        assert k.mu == sk.consts.mu_earth
+        assert k_moon.mu == sk.consts.mu_moon
+        assert k_moon.period / k.period == pytest.approx(
+            m.sqrt(sk.consts.mu_earth / sk.consts.mu_moon), rel=1e-12
+        )
+        assert k_moon.period == pytest.approx(8022.0, abs=5.0)
+        # from_pv with the same mu round-trips; with Earth's it does not
+        r, v = k_moon.to_pv()
+        back = sk.kepler.from_pv(r, v, mu=sk.consts.mu_moon)
+        assert back.mu == sk.consts.mu_moon
+        assert back.a == pytest.approx(k_moon.a, rel=1e-9)
+        assert abs(sk.kepler.from_pv(r, v).a - k_moon.a) / k_moon.a > 0.1
+        with pytest.raises(ValueError):
+            sk.kepler.from_pv(r, v, mu=0.0)
+        # propagate keeps mu, and the mu setter works
+        assert k_moon.propagate(10.0).mu == sk.consts.mu_moon
+        k.mu = sk.consts.mu_moon
+        assert k == k_moon
+
+    def test_kepler_derived_helpers(self):
+        a, e = 26600e3, 0.74
+        k = sk.kepler(a, e, 1.1, 5.0, 4.0, 0.0)
+        assert k.periapsis == pytest.approx(a * (1 - e))
+        assert k.apoapsis == pytest.approx(a * (1 + e))
+        assert k.specific_energy == pytest.approx(-sk.consts.mu_earth / (2 * a), rel=1e-12)
+        r, v = k.to_pv()
+        assert k.angular_momentum == pytest.approx(np.linalg.norm(np.cross(r, v)), rel=1e-12)
+        assert k.flight_path_angle == 0.0
+        k_out = sk.kepler(a, e, 1.1, 5.0, 4.0, 1.0)
+        r, v = k_out.to_pv()
+        gamma = m.asin(np.dot(r, v) / (np.linalg.norm(r) * np.linalg.norm(v)))
+        assert k_out.flight_path_angle == pytest.approx(gamma, abs=1e-12)
+        assert k_out.flight_path_angle > 0
+        assert k_out.argument_of_latitude == pytest.approx(5.0, abs=1e-12)
+        assert k_out.true_longitude == pytest.approx(10.0 - 2 * m.pi, abs=1e-12)
+
+    def test_kepler_repr_and_satstate_from_kepler(self):
+        k = sk.kepler(7000e3, 0.1, 0.5, 1.0, 0.3, 0.7)
+        text = repr(k)
+        assert text.startswith("kepler(a=7.000000e6") and "argp=" in text and "mu=" in text
+        assert "\n" not in text
+        assert str(k).startswith("Keplerian Elements:")
+        t0 = sk.time(2024, 1, 1)
+        s = sk.satstate.from_kepler(t0, k)
+        r, v = k.to_pv()
+        np.testing.assert_allclose(s.pos, r)
+        np.testing.assert_allclose(s.vel, v)
+        assert s.time == t0
 
     def test_kepler_propagate_accepts_int_and_duration(self):
         k = sk.kepler(7000e3, 0.1, 0.5, 1.0, 0.3, 0.7)
diff --git a/src/kepler.rs b/src/kepler.rs
index f7d8a5d3..52287881 100644
--- a/src/kepler.rs
+++ b/src/kepler.rs
@@ -5,6 +5,7 @@ use thiserror::Error;
 
 /// Errors that can occur while constructing or converting [`Kepler`] elements.
 #[derive(Debug, Error)]
+#[non_exhaustive]
 pub enum Error {
     /// Returned by [`Kepler::from_pv`] when the computed eccentricity is
     /// outside the valid range for an elliptical orbit.
@@ -16,6 +17,17 @@ pub enum Error {
     /// therefore inclination and RAAN — are undefined.
     #[error("Degenerate state: angular momentum is zero (rectilinear trajectory)")]
     Degenerate,
+
+    /// Returned by [`Kepler::try_new`] and [`Kepler::validate`] for an
+    /// element outside its domain: a non-finite value, `a <= 0`, `eccen`
+    /// outside `[0, 1)`, `incl` outside `[0, π]`, or `mu <= 0`.
+    #[error("invalid Keplerian element {name} = {value}: {reason}")]
+    #[non_exhaustive]
+    InvalidElement {
+        name: &'static str,
+        value: f64,
+        reason: &'static str,
+    },
 }
 
 /// Convenient type alias used throughout the `kepler` module.
@@ -45,32 +57,78 @@ pub type KeplerError = Error;
 ///   and perpendicular to the semimajor axis.  The eccentric anomaly is
 ///   a useful prerequisite to compute the mean anomaly
 ///
+#[derive(Debug, Clone, Copy, PartialEq)]
 pub enum Anomaly {
     Mean(f64),
     True(f64),
     Eccentric(f64),
 }
 
+impl Anomaly {
+    /// The angle carried by the variant, radians.
+    pub const fn value(self) -> f64 {
+        match self {
+            Self::Mean(v) | Self::True(v) | Self::Eccentric(v) => v,
+        }
+    }
+}
+
 // External library imports
 use crate::mathtypes::*;
 
 /// Keplerian Orbital Elements
 ///
-/// The 6 Keplerian orbital elements are:
-/// a: semi-major axis, meters
-/// eccen: Eccentricity
-/// incl: Inclination, radians
-/// RAAN: Right Ascension of the Ascending Node, radians
-/// w: Argument of Perigee, radians
-/// an: Anomaly of given type, radians
-#[derive(Debug, Clone, Copy)]
+/// The 6 Keplerian orbital elements, plus the gravitational parameter of the
+/// central body they refer to:
+///
+/// * `a`: semi-major axis, meters
+/// * `eccen`: eccentricity, `0 <= eccen < 1`
+/// * `incl`: inclination, radians, `0 <= incl <= π`
+/// * `raan`: right ascension of the ascending node, radians
+/// * `argp`: argument of periapsis, radians
+/// * `nu`: true anomaly, radians
+/// * `mu`: gravitational parameter of the central body, m³/s²
+///   ([`MU_EARTH`](crate::consts::MU_EARTH) unless set with [`Kepler::with_mu`])
+///
+/// The fields are public and may be assigned directly; nothing is validated
+/// on assignment. [`Kepler::try_new`] and [`Kepler::validate`] are the
+/// checked paths.
+#[derive(Debug, Clone, Copy, PartialEq, serde::Serialize, serde::Deserialize)]
 pub struct Kepler {
     pub a: f64,
     pub eccen: f64,
     pub incl: f64,
     pub raan: f64,
-    pub w: f64,
-    pub nu: f64, // True anomaly
+    pub argp: f64,
+    /// True anomaly
+    pub nu: f64,
+    /// Gravitational parameter, m³/s². Serialized element sets that predate
+    /// this field deserialize with Earth's value.
+    #[serde(default = "default_mu")]
+    pub mu: f64,
+}
+
+const fn default_mu() -> f64 {
+    crate::consts::MU_EARTH
+}
+
+/// `Ok(())` when `value` is finite and `ok` holds; the error names the element.
+fn check(name: &'static str, value: f64, ok: bool, reason: &'static str) -> Result<()> {
+    if !value.is_finite() {
+        return Err(Error::InvalidElement {
+            name,
+            value,
+            reason: "must be finite",
+        });
+    }
+    if !ok {
+        return Err(Error::InvalidElement {
+            name,
+            value,
+            reason,
+        });
+    }
+    Ok(())
 }
 
 // Convert mean to eccentric anomaly
@@ -138,18 +196,87 @@ impl Kepler {
     ///
     /// # Returns
     ///
-    /// * `Kepler` - A new Keplerian orbital element object
+    /// * `Kepler` - A new Keplerian orbital element object, with Earth's
+    ///   gravitational parameter (see [`Kepler::with_mu`])
+    ///
+    /// Nothing is validated: `eccen >= 1`, `a <= 0` or non-finite inputs
+    /// produce meaningless (NaN) anomaly conversions rather than an error.
+    /// Use [`Kepler::try_new`] for a checked constructor.
     pub fn new(a: f64, eccen: f64, i: f64, raan: f64, argp: f64, an: Anomaly) -> Self {
         Self {
             a,
             eccen,
             incl: i,
             raan,
-            w: argp,
+            argp,
             nu: to_trueanomaly(an, eccen),
+            mu: default_mu(),
         }
     }
 
+    /// Checked constructor: [`Kepler::new`] after validating every input.
+    ///
+    /// # Errors
+    ///
+    /// [`Error::InvalidElement`] when any input is non-finite, `a <= 0`,
+    /// `eccen` is outside `[0, 1)`, or `i` is outside `[0, π]`. The bounds
+    /// are strict: `eccen = 1` and `i = π + 1e-16` are rejected.
+    pub fn try_new(a: f64, eccen: f64, i: f64, raan: f64, argp: f64, an: Anomaly) -> Result<Self> {
+        check("raan", raan, true, "")?;
+        check("argp", argp, true, "")?;
+        check("anomaly", an.value(), true, "")?;
+        let k = Self::new(a, eccen, i, raan, argp, an);
+        k.validate()?;
+        Ok(k)
+    }
+
+    /// Check that the stored elements describe a closed orbit: every field
+    /// finite, `a > 0`, `0 <= eccen < 1`, `0 <= incl <= π`, `mu > 0`.
+    ///
+    /// This is the check [`Kepler::try_new`] applies; call it after
+    /// assigning fields directly.
+    pub fn validate(&self) -> Result<()> {
+        check(
+            "a",
+            self.a,
+            self.a > 0.0,
+            "semi-major axis must be positive",
+        )?;
+        check(
+            "eccen",
+            self.eccen,
+            (0.0..1.0).contains(&self.eccen),
+            "eccentricity must be in [0, 1) (closed orbits only)",
+        )?;
+        check(
+            "incl",
+            self.incl,
+            (0.0..=std::f64::consts::PI).contains(&self.incl),
+            "inclination must be in [0, π] radians",
+        )?;
+        check("raan", self.raan, true, "")?;
+        check("argp", self.argp, true, "")?;
+        check("nu", self.nu, true, "")?;
+        check(
+            "mu",
+            self.mu,
+            self.mu > 0.0,
+            "gravitational parameter must be positive",
+        )?;
+        Ok(())
+    }
+
+    /// The same elements referred to a central body with gravitational
+    /// parameter `mu` (m³/s²), e.g. [`MU_MOON`](crate::consts::MU_MOON).
+    ///
+    /// Only the dynamics change (mean motion, period, `propagate`, the
+    /// element ↔ state conversions); the six geometric elements are kept
+    /// as they are. Not validated; see [`Kepler::validate`].
+    pub const fn with_mu(mut self, mu: f64) -> Self {
+        self.mu = mu;
+        self
+    }
+
     /// Create a new Keplerian orbital element object with true anomaly
     ///
     /// # Arguments
@@ -220,14 +347,7 @@ impl Kepler {
         let n = self.mean_motion();
         let ma = n.mul_add(dt.as_seconds(), self.mean_anomaly());
         let nu = mean2true(ma, self.eccen);
-        Self {
-            a: self.a,
-            eccen: self.eccen,
-            incl: self.incl,
-            raan: self.raan,
-            w: self.w,
-            nu,
-        }
+        Self { nu, ..*self }
     }
 
     /// Return the eccentric anomaly of the satellite in radians
@@ -255,7 +375,7 @@ impl Kepler {
     ///
     /// * `f64` - Mean motion, radians/second
     pub fn mean_motion(&self) -> f64 {
-        (crate::consts::MU_EARTH / self.a.powi(3)).sqrt()
+        (self.mu / self.a.powi(3)).sqrt()
     }
 
     /// Return the period of the satellite in seconds
@@ -267,7 +387,52 @@ impl Kepler {
         2.0 * std::f64::consts::PI / self.mean_motion()
     }
 
-    /// Convert Cartesian coordinates to Keplerian orbital elements
+    /// Radius of periapsis `a (1 - e)`, meters
+    pub fn periapsis(&self) -> f64 {
+        self.a * (1.0 - self.eccen)
+    }
+
+    /// Radius of apoapsis `a (1 + e)`, meters
+    pub fn apoapsis(&self) -> f64 {
+        self.a * (1.0 + self.eccen)
+    }
+
+    /// Specific orbital energy `-μ / 2a`, J/kg (m²/s²)
+    pub fn specific_energy(&self) -> f64 {
+        -self.mu / (2.0 * self.a)
+    }
+
+    /// Magnitude of the specific angular momentum `√(μ p)`, m²/s
+    pub fn angular_momentum(&self) -> f64 {
+        (self.mu * self.semiparameter()).sqrt()
+    }
+
+    /// Flight-path angle `γ = atan2(e sin ν, 1 + e cos ν)`, radians: the
+    /// angle of the velocity above the local horizontal, zero at periapsis
+    /// and apoapsis, positive while climbing.
+    pub fn flight_path_angle(&self) -> f64 {
+        f64::atan2(
+            self.eccen * self.nu.sin(),
+            self.eccen.mul_add(self.nu.cos(), 1.0),
+        )
+    }
+
+    /// Argument of latitude `u = ω + ν`, radians, reduced to `[0, 2π)`.
+    /// Well defined for circular orbits, where ω and ν separately are not.
+    pub fn argument_of_latitude(&self) -> f64 {
+        (self.argp + self.nu).rem_euclid(std::f64::consts::TAU)
+    }
+
+    /// True longitude `λ = Ω + ω + ν`, radians, reduced to `[0, 2π)`.
+    /// Well defined for circular equatorial orbits, where Ω, ω and ν
+    /// separately are not.
+    pub fn true_longitude(&self) -> f64 {
+        (self.raan + self.argp + self.nu).rem_euclid(std::f64::consts::TAU)
+    }
+
+    /// Convert Cartesian coordinates to Keplerian orbital elements about
+    /// the Earth ([`MU_EARTH`](crate::consts::MU_EARTH)); see
+    /// [`Kepler::from_pv_with_mu`].
     ///
     /// # Arguments
     ///
@@ -279,8 +444,21 @@ impl Kepler {
     /// * `Kepler` - A new Keplerian orbital element object
     ///
     pub fn from_pv(r: Vector3, v: Vector3) -> Result<Self> {
+        Self::from_pv_with_mu(r, v, default_mu())
+    }
+
+    /// Convert Cartesian coordinates to Keplerian orbital elements about a
+    /// central body with gravitational parameter `mu` (m³/s²).
+    ///
+    /// The returned elements carry `mu`, so their period, `propagate` and
+    /// `to_pv` refer to the same body.
+    ///
+    /// # Errors
+    ///
+    /// [`Error::Degenerate`] for (near-)zero angular momentum,
+    /// [`Error::EccenOutOfBound`] for an open (parabolic/hyperbolic) state.
+    pub fn from_pv_with_mu(r: Vector3, v: Vector3, mu: f64) -> Result<Self> {
         use std::f64::consts::TAU;
-        let mu = crate::consts::MU_EARTH;
         let rmag = r.norm();
 
         let h = r.cross(&v);
@@ -345,7 +523,7 @@ impl Kepler {
             (raan, w, nu)
         };
 
-        Ok(Self::new(a, eccen, incl, raan, w, Anomaly::True(nu)))
+        Ok(Self::new(a, eccen, incl, raan, w, Anomaly::True(nu)).with_mu(mu))
     }
 
     /// Convert Keplerian orbital elements to Cartesian coordinates
@@ -359,9 +537,9 @@ impl Kepler {
         let r = p / self.eccen.mul_add(self.nu.cos(), 1.0);
         let r_pqw = numeris::vector![r * self.nu.cos(), r * self.nu.sin(), 0.0];
         let v_pqw = numeris::vector![-self.nu.sin(), self.eccen + self.nu.cos(), 0.0]
-            * (crate::consts::MU_EARTH / p).sqrt();
+            * (self.mu / p).sqrt();
         let q =
-            Quaternion::rotz(self.raan) * Quaternion::rotx(self.incl) * Quaternion::rotz(self.w);
+            Quaternion::rotz(self.raan) * Quaternion::rotx(self.incl) * Quaternion::rotz(self.argp);
         (q * r_pqw, q * v_pqw)
     }
 }
@@ -375,8 +553,8 @@ impl std::fmt::Display for Kepler {
         )?;
         write!(
             f,
-            "  Ω = {:.3} rad\n  ω = {:.3} rad\n  ν = {:.3} rad\n",
-            self.raan, self.w, self.nu
+            "  Ω = {:.3} rad\n  ω = {:.3} rad\n  ν = {:.3} rad\n  μ = {:.6e} m³/s²\n",
+            self.raan, self.argp, self.nu, self.mu
         )
     }
 }
@@ -425,7 +603,7 @@ mod tests {
         let k = Kepler::new(a, 0.0, 0.0, 0.0, 0.0, Anomaly::True(0.7));
         let (r, v) = k.to_pv();
         let k2 = Kepler::from_pv(r, v).unwrap();
-        assert!(k2.raan.is_finite() && k2.w.is_finite() && k2.nu.is_finite());
+        assert!(k2.raan.is_finite() && k2.argp.is_finite() && k2.nu.is_finite());
         let (r2, v2) = k2.to_pv();
         assert!((r - r2).norm() / r.norm() < 1.0e-9);
         assert!((v - v2).norm() / v.norm() < 1.0e-9);
@@ -439,7 +617,7 @@ mod tests {
         let k = Kepler::new(a, 0.1, 0.0, 0.0, 0.5, Anomaly::True(1.0));
         let (r, v) = k.to_pv();
         let k2 = Kepler::from_pv(r, v).unwrap();
-        assert!(k2.raan.is_finite() && k2.w.is_finite() && k2.nu.is_finite());
+        assert!(k2.raan.is_finite() && k2.argp.is_finite() && k2.nu.is_finite());
         let (r2, v2) = k2.to_pv();
         assert!((r - r2).norm() / r.norm() < 1.0e-9);
         assert!((v - v2).norm() / v.norm() < 1.0e-9);
@@ -592,7 +770,7 @@ mod tests {
                             && k2.eccen.is_finite()
                             && k2.incl.is_finite()
                             && k2.raan.is_finite()
-                            && k2.w.is_finite()
+                            && k2.argp.is_finite()
                             && k2.nu.is_finite(),
                         "non-finite element e={e} i={i} nu={nu}: {k2:?}"
                     );
@@ -640,6 +818,123 @@ mod tests {
         }
     }
 
+    #[test]
+    fn test_try_new_rejects_out_of_domain_elements() {
+        use std::f64::consts::PI;
+        let ok = Kepler::try_new(7000.0e3, 0.1, 0.5, 1.0, 0.3, Anomaly::True(0.7));
+        assert!(ok.is_ok());
+        let bad = [
+            (0.0, 0.1, 0.5, "a"),
+            (-7000.0e3, 0.1, 0.5, "a"),
+            (f64::NAN, 0.1, 0.5, "a"),
+            (7000.0e3, 1.0, 0.5, "eccen"),
+            (7000.0e3, -1.0e-3, 0.5, "eccen"),
+            (7000.0e3, f64::INFINITY, 0.5, "eccen"),
+            (7000.0e3, 0.1, -1.0e-9, "incl"),
+            (7000.0e3, 0.1, PI + 1.0e-9, "incl"),
+            (7000.0e3, 0.1, f64::NAN, "incl"),
+        ];
+        for (a, e, i, which) in bad {
+            match Kepler::try_new(a, e, i, 1.0, 0.3, Anomaly::True(0.7)) {
+                Err(Error::InvalidElement { name, .. }) => assert_eq!(name, which),
+                other => panic!("expected InvalidElement({which}), got {other:?}"),
+            }
+        }
+        // Non-finite angles, including the anomaly, are rejected too.
+        assert!(Kepler::try_new(7000.0e3, 0.1, 0.5, f64::NAN, 0.3, Anomaly::True(0.7)).is_err());
+        assert!(Kepler::try_new(7000.0e3, 0.1, 0.5, 1.0, 0.3, Anomaly::Mean(f64::NAN)).is_err());
+        // Boundaries: e = 0 and i ∈ {0, π} are valid.
+        assert!(Kepler::try_new(7000.0e3, 0.0, 0.0, 0.0, 0.0, Anomaly::True(0.0)).is_ok());
+        assert!(Kepler::try_new(7000.0e3, 0.0, PI, 0.0, 0.0, Anomaly::True(0.0)).is_ok());
+        // validate() catches a bad direct assignment, mu included.
+        let mut k = ok.unwrap();
+        k.mu = 0.0;
+        assert!(matches!(
+            k.validate(),
+            Err(Error::InvalidElement { name: "mu", .. })
+        ));
+    }
+
+    #[test]
+    fn test_mu_changes_dynamics_not_geometry() {
+        use crate::consts::{MU_EARTH, MU_MOON};
+        let k_earth = Kepler::new(2000.0e3, 0.05, 1.0, 0.2, 0.3, Anomaly::True(0.4));
+        assert_eq!(k_earth.mu, MU_EARTH);
+        let k_moon = k_earth.with_mu(MU_MOON);
+        // Same six elements …
+        assert_eq!(k_moon.a, k_earth.a);
+        assert_eq!(k_moon.nu, k_earth.nu);
+        // … different period: T ∝ 1/√μ.
+        let ratio = k_moon.period() / k_earth.period();
+        assert!((ratio - (MU_EARTH / MU_MOON).sqrt()).abs() < 1.0e-12);
+        // 2000 km about the Moon: ~ 2.1 h. Sanity check the magnitude.
+        assert!(
+            (k_moon.period() - 8022.0).abs() < 5.0,
+            "{}",
+            k_moon.period()
+        );
+        // The state round trip must use the same μ on both legs.
+        let (r, v) = k_moon.to_pv();
+        let back = Kepler::from_pv_with_mu(r, v, MU_MOON).unwrap();
+        assert_eq!(back.mu, MU_MOON);
+        assert!((back.a - k_moon.a).abs() / k_moon.a < 1.0e-9);
+        assert!((back.eccen - k_moon.eccen).abs() < 1.0e-9);
+        // Interpreting a lunar state with Earth's μ gives a different orbit.
+        let wrong = Kepler::from_pv(r, v).unwrap();
+        assert!((wrong.a - k_moon.a).abs() / k_moon.a > 0.1);
+        // propagate keeps mu.
+        assert_eq!(
+            k_moon.propagate(&crate::Duration::from_seconds(10.0)).mu,
+            MU_MOON
+        );
+    }
+
+    #[test]
+    fn test_derived_helpers() {
+        use std::f64::consts::{PI, TAU};
+        let a = 26_600.0e3;
+        let e = 0.74;
+        let k = Kepler::new(a, e, 1.1, 5.0, 4.0, Anomaly::True(0.0));
+        assert!((k.periapsis() - a * (1.0 - e)).abs() < 1.0e-6);
+        assert!((k.apoapsis() - a * (1.0 + e)).abs() < 1.0e-6);
+        assert!((k.periapsis() + k.apoapsis() - 2.0 * a).abs() < 1.0e-6);
+        // Vis-viva at periapsis agrees with the energy helper.
+        let (r, v) = k.to_pv();
+        let xi = v.norm_squared() / 2.0 - k.mu / r.norm();
+        assert!((xi - k.specific_energy()).abs() / xi.abs() < 1.0e-12);
+        // |r × v| agrees with the angular-momentum helper.
+        assert!((r.cross(&v).norm() - k.angular_momentum()).abs() / k.angular_momentum() < 1.0e-12);
+        // Flight-path angle: zero at periapsis and apoapsis, positive on the
+        // outbound leg, and equal to atan(e sinν / (1 + e cosν)) elsewhere.
+        assert_eq!(k.flight_path_angle(), 0.0);
+        let k_apo = Kepler::new(a, e, 1.1, 5.0, 4.0, Anomaly::True(PI));
+        assert!(k_apo.flight_path_angle().abs() < 1.0e-15);
+        let k_out = Kepler::new(a, e, 1.1, 5.0, 4.0, Anomaly::True(1.0));
+        assert!(k_out.flight_path_angle() > 0.0);
+        let (r, v) = k_out.to_pv();
+        let gamma = (r.dot(&v) / (r.norm() * v.norm())).asin();
+        assert!((gamma - k_out.flight_path_angle()).abs() < 1.0e-12);
+        // u and λ wrap into [0, 2π): ω + ν = 4 + 1 = 5; Ω + ω + ν = 10 → 10 - 2π.
+        assert!((k_out.argument_of_latitude() - 5.0).abs() < 1.0e-12);
+        assert!((k_out.true_longitude() - (10.0 - TAU)).abs() < 1.0e-12);
+        let k_neg = Kepler::new(a, e, 1.1, -1.0, -1.0, Anomaly::True(-1.0));
+        assert!((k_neg.argument_of_latitude() - (TAU - 2.0)).abs() < 1.0e-12);
+        assert!((k_neg.true_longitude() - (TAU - 3.0)).abs() < 1.0e-12);
+    }
+
+    #[test]
+    fn test_serde_roundtrip_and_missing_mu_defaults_to_earth() {
+        let k = Kepler::new(7000.0e3, 0.1, 0.5, 1.0, 0.3, Anomaly::True(0.7))
+            .with_mu(crate::consts::MU_MOON);
+        let json = serde_json::to_string(&k).unwrap();
+        let back: Kepler = serde_json::from_str(&json).unwrap();
+        assert_eq!(back, k);
+        // An element set serialized before `mu` existed still loads.
+        let legacy = r#"{"a":7000000.0,"eccen":0.1,"incl":0.5,"raan":1.0,"argp":0.3,"nu":0.7}"#;
+        let old: Kepler = serde_json::from_str(legacy).unwrap();
+        assert_eq!(old.mu, crate::consts::MU_EARTH);
+    }
+
     #[test]
     fn test_topv() {
         // Example 2-6 from Vallado
@@ -671,7 +966,7 @@ mod tests {
         assert!((k.eccen - 0.83285).abs() < 1e-3);
         assert!((k.incl - 87.87_f64.to_radians()).abs() < 1e-3);
         assert!((k.raan - 227.89_f64.to_radians()).abs() < 1e-3);
-        assert!((k.w - 53.38_f64.to_radians()).abs() < 1e-3);
+        assert!((k.argp - 53.38_f64.to_radians()).abs() < 1e-3);
         assert!((k.nu - 92.335_f64.to_radians()).abs() < 1e-3);
     }
 }
diff --git a/src/orbitprop/satstate.rs b/src/orbitprop/satstate.rs
index 58d90cb0..915b6fac 100644
--- a/src/orbitprop/satstate.rs
+++ b/src/orbitprop/satstate.rs
@@ -205,6 +205,17 @@ impl SatState {
         }
     }
 
+    /// Create a satellite state from Keplerian elements: the two-body
+    /// position and velocity of `kepler` ([`Kepler::to_pv`](crate::kepler::Kepler::to_pv))
+    /// taken as GCRF at `time`, with no covariance and no maneuvers.
+    ///
+    /// The elements are treated as osculating GCRF elements; `kepler.mu`
+    /// should therefore be Earth's.
+    pub fn from_kepler<T: TimeLike>(time: &T, kepler: &crate::kepler::Kepler) -> Self {
+        let (pos, vel) = kepler.to_pv();
+        Self::from_pv(time, &pos, &vel)
+    }
+
     /// Position vector in GCRF [meters]
     pub fn pos_gcrf(&self) -> Vector3 {
         self.pv.block::<3, 1>(0, 0)
diff --git a/src/tle/fitting.rs b/src/tle/fitting.rs
index d2dccca8..65d3d145 100644
--- a/src/tle/fitting.rs
+++ b/src/tle/fitting.rs
@@ -272,7 +272,7 @@ impl TLE {
             kepler.incl.to_degrees(),
             kepler.eccen,
             kepler.raan.to_degrees(),
-            kepler.w.to_degrees(),
+            kepler.argp.to_degrees(),
             kepler.mean_motion() * 60.0 * 60.0 * 24.0 / (2.0 * std::f64::consts::PI),
             kepler.mean_anomaly().to_degrees(),
             0.0,