fix(orbit): RA wrap, cos(dec) rate, asin clamp, robust Kepler (CORR-07, CORR-08)

output/builds/src/orbit_propagator.py

@@ -76,7 +76,8 @@ def v_mag(self) -> float:
     def to_ra_dec(self) -> Tuple[float, float]:
         """Convert ECI position to RA/Dec (degrees)."""
         r = self.r_mag
-        dec = math.asin(self.z / r) * RAD2DEG
+        ratio = max(-1.0, min(1.0, self.z / r)) if r > 0 else 0.0
+        dec = math.asin(ratio) * RAD2DEG
         ra = math.atan2(self.y, self.x) * RAD2DEG
         if ra < 0:
             ra += 360
@@ -190,13 +191,19 @@ def propagate_window(self, elem: OrbitalElements,
 
     def _solve_kepler(self, M: float, e: float, tol: float = 1e-12) -> float:
         """Solve Kepler's equation iteratively (Newton-Raphson)."""
-        E = M  # Initial guess
-        for _ in range(50):
-            dE = (E - e * math.sin(E) - M) / (1 - e * math.cos(E))
+        if not 0 <= e < 1:
+            raise ValueError("Eccentricity must be in the range [0, 1)")
+
+        # Work near the origin for numerical stability, then restore full turns.
+        wrapped_M = (M + math.pi) % (2 * math.pi) - math.pi
+        turns = M - wrapped_M
+        E = wrapped_M if e < 0.8 else math.copysign(math.pi, wrapped_M or 1.0)
+        for _ in range(100):
+            dE = (E - e * math.sin(E) - wrapped_M) / (1 - e * math.cos(E))
             E -= dE
             if abs(dE) < tol:
                 break
-        return E
+        return E + turns
 
     def estimate_orbit_regime(self, observations: list) -> Optional[OrbitalElements]:
         """
@@ -222,10 +229,10 @@ def estimate_orbit_regime(self, observations: list) -> Optional[OrbitalElements]
             return None
 
         # Approximate angular rate to classify orbit regime
-        rate_deg_sec = math.sqrt(
-            (obs[2]["ra_deg"] - obs[0]["ra_deg"])**2 +
-            (obs[2]["dec_deg"] - obs[0]["dec_deg"])**2
-        ) / dt
+        dra = ((obs[2]["ra_deg"] - obs[0]["ra_deg"] + 180) % 360 - 180)
+        dra *= math.cos(math.radians(obs[1]["dec_deg"]))
+        ddec = obs[2]["dec_deg"] - obs[0]["dec_deg"]
+        rate_deg_sec = math.hypot(dra, ddec) / dt
 
         # Higher rate = lower orbit (Kepler's third law)
         if rate_deg_sec > 0.5:

tests/test_orbit_propagator_fixes.py

@@ -0,0 +1,39 @@
+import math
+
+from output.builds.src import orbit_propagator
+
+
+def test_orbit_regime_rate_handles_ra_wrap_and_declination_scale():
+    observations = [
+        {"timestamp": 0.0, "ra_deg": 359.4, "dec_deg": 60.0},
+        {"timestamp": 1.0, "ra_deg": 0.0, "dec_deg": 60.0},
+        {"timestamp": 2.0, "ra_deg": 0.6, "dec_deg": 60.0},
+    ]
+
+    orbit = orbit_propagator.OrbitPropagator().estimate_orbit_regime(observations)
+
+    assert orbit is not None
+    assert math.isclose(orbit.altitude_km, 1000.0, abs_tol=1e-9)
+
+
+def test_state_vector_clamps_declination_roundoff():
+    class RoundedState(orbit_propagator.StateVector):
+        @property
+        def r_mag(self):
+            return abs(self.z) * (1.0 - 1e-15)
+
+    _, declination = RoundedState(z=1.0).to_ra_dec()
+
+    assert declination == 90.0
+
+
+def test_kepler_solver_converges_for_high_eccentricity():
+    mean_anomaly = 0.1
+    eccentricity = 0.95
+
+    eccentric_anomaly = orbit_propagator.OrbitPropagator()._solve_kepler(
+        mean_anomaly, eccentricity
+    )
+
+    residual = eccentric_anomaly - eccentricity * math.sin(eccentric_anomaly) - mean_anomaly
+    assert abs(residual) < 1e-12