Fix HIGH correctness bugs CORR-01..05

output/builds/src/conjunction_screener.py

@@ -121,15 +121,18 @@ def screen_pair(self,
         if miss_km > self.screen_distance_km:
             return None
 
-        # --- Combine covariances ---
-        # In RTN frame, covariances add (assuming independence)
-        C1 = cov1.to_matrix()
-        C2 = cov2.to_matrix()
-        C_combined = C1 + C2
+        # --- Combine covariances in a common frame ---
+        # Each RTN covariance is expressed in its object's local orbital frame.
+        # Rotate both to ECI before adding independent uncertainties.
+        M1 = self._rtn_to_eci_matrix(pos1, vel1)
+        M2 = self._rtn_to_eci_matrix(pos2, vel2)
+        C1_eci = M1 @ cov1.to_matrix() @ M1.T
+        C2_eci = M2 @ cov2.to_matrix() @ M2.T
+        C_combined_eci = C1_eci + C2_eci
 
         # --- Project into encounter plane ---
         C_2d, miss_2d = self._project_to_encounter_plane(
-            dr, dv, C_combined
+            dr, dv, C_combined_eci
         )
 
         # --- Compute collision probability ---
@@ -174,6 +177,21 @@ def screen_pair(self,
 
         return event
 
+    @staticmethod
+    def _rtn_to_eci_matrix(position: np.ndarray,
+                           velocity: np.ndarray) -> np.ndarray:
+        """Return the direction-cosine matrix that maps RTN vectors to ECI."""
+        r_norm = np.linalg.norm(position)
+        angular_momentum = np.cross(position, velocity)
+        h_norm = np.linalg.norm(angular_momentum)
+        if r_norm < 1e-12 or h_norm < 1e-12:
+            raise ValueError("Cannot construct RTN frame from degenerate state")
+
+        radial = position / r_norm
+        normal = angular_momentum / h_norm
+        transverse = np.cross(normal, radial)
+        return np.column_stack((radial, transverse, normal))
+
     def _project_to_encounter_plane(self,
                                       dr: np.ndarray,
                                       dv: np.ndarray,

output/builds/src/hardware/ascom_hal.py

@@ -519,7 +519,8 @@ def run_self_test():
     # 10. LEO tracking benchmark
     print("\n─── LEO Tracking Benchmark ───")
     altitudes = [200, 300, 420, 550, 800, 1200, 2000]
-    print(f"  {'Alt (km)':>10s}  {'Rate (°/s)':>10s}  {'RA (\"/s)':>10s}  {'Dec (\"/s)':>10s}  {'Trackable':>10s}")
+    ra_header, dec_header = 'RA ("/s)', 'Dec ("/s)'
+    print(f"  {'Alt (km)':>10s}  {'Rate (°/s)':>10s}  {ra_header:>10s}  {dec_header:>10s}  {'Trackable':>10s}")
     for alt in altitudes:
         r = TrackingRates.from_leo_pass(alt, 51.6, 60)
         ok = "✓ Yes" if r.total_rate_deg_per_sec <= 5.0 else "✗ Too fast"

output/builds/src/science/bayesian_iod.py

@@ -71,13 +71,22 @@ def compute_bounds(self, ra: float, dec: float,
 
         obs_r = np.linalg.norm(observer_pos)
 
-        # Minimum range: object at min altitude along LOS
-        rho_min = max(0.1, self.min_alt + R_EARTH - obs_r)
-        # Maximum range: object at max altitude along LOS
-        rho_max = self.max_alt + R_EARTH + obs_r
+        # Intersect the observation ray with the minimum/maximum altitude
+        # spheres: |observer_pos + rho * los|^2 = radius^2.
+        b = 2.0 * float(np.dot(observer_pos, los))
+        c_min = obs_r**2 - (R_EARTH + self.min_alt)**2
+        c_max = obs_r**2 - (R_EARTH + self.max_alt)**2
+        discriminant_min = b**2 - 4.0 * c_min
+        discriminant_max = b**2 - 4.0 * c_max
+        if discriminant_min < 0.0 or discriminant_max < 0.0:
+            raise ValueError("Line of sight does not intersect the admissible altitude shell")
+
+        rho_min = max(0.1, (-b + np.sqrt(discriminant_min)) / 2.0)
+        rho_max = (-b + np.sqrt(discriminant_max)) / 2.0
+        if rho_max <= rho_min:
+            raise ValueError("Line of sight has no positive admissible range")
 
         # Range-rate bounds from escape velocity constraint
-        r_max = obs_r + rho_max
         v_esc = np.sqrt(2 * MU_EARTH / (R_EARTH + self.min_alt))
         rhodot_min = -v_esc
         rhodot_max = v_esc

output/builds/src/science/non_gaussian_pc.py

@@ -81,6 +81,53 @@ class StateTransitionTensor:
     def __init__(self, order: int = 3):
         self.order = min(order, 4)
 
+    def _acceleration_eci(self, position: np.ndarray) -> np.ndarray:
+        """Return two-body plus J2 acceleration in the ECI frame."""
+        radius = np.linalg.norm(position)
+        if radius < 100.0:
+            raise ValueError("ECI position magnitude is too small")
+
+        x, y, z = position
+        radius_sq = radius * radius
+        z_ratio_sq = (z * z) / radius_sq
+        acceleration = -self.MU_EARTH * position / radius**3
+        j2_factor = (
+            1.5 * self.J2 * self.MU_EARTH * self.R_EARTH**2 / radius**5
+        )
+        acceleration += j2_factor * np.array([
+            x * (5.0 * z_ratio_sq - 1.0),
+            y * (5.0 * z_ratio_sq - 1.0),
+            z * (5.0 * z_ratio_sq - 3.0),
+        ])
+        return acceleration
+
+    def _variational_derivative(self, augmented_state: np.ndarray) -> np.ndarray:
+        """Evaluate ECI dynamics and d(Phi)/dt = A(t) Phi."""
+        state = augmented_state[:6]
+        position = state[:3]
+        velocity = state[3:]
+        phi = augmented_state[6:].reshape(6, 6)
+
+        radius = np.linalg.norm(position)
+        delta = max(1e-4, radius * 1e-7)
+        acceleration_jacobian = np.empty((3, 3))
+        for axis in range(3):
+            offset = np.zeros(3)
+            offset[axis] = delta
+            acceleration_jacobian[:, axis] = (
+                self._acceleration_eci(position + offset)
+                - self._acceleration_eci(position - offset)
+            ) / (2.0 * delta)
+
+        dynamics_jacobian = np.zeros((6, 6))
+        dynamics_jacobian[:3, 3:] = np.eye(3)
+        dynamics_jacobian[3:, :3] = acceleration_jacobian
+        return np.concatenate((
+            velocity,
+            self._acceleration_eci(position),
+            (dynamics_jacobian @ phi).ravel(),
+        ))
+
     def compute_stm(self, state: np.ndarray, dt: float) -> np.ndarray:
         """
         Compute 6x6 State Transition Matrix via two-body + J2 linearization.
@@ -89,44 +136,30 @@ def compute_stm(self, state: np.ndarray, dt: float) -> np.ndarray:
             state: [x, y, z, vx, vy, vz] in ECI (km, km/s)
             dt: propagation time (seconds)
         """
-        r = np.linalg.norm(state[:3])
-        if r < 100:
+        state = np.asarray(state, dtype=float)
+        if state.shape != (6,) or not np.all(np.isfinite(state)):
+            raise ValueError("state must be a finite six-element ECI vector")
+        if np.linalg.norm(state[:3]) < 100:
+            return np.eye(6)
+        if not math.isfinite(dt):
+            raise ValueError("dt must be finite")
+        if dt == 0.0:
             return np.eye(6)
 
-        # Two-body Keplerian STM (clohessy-wiltshire approximation for near-circular)
-        n = math.sqrt(self.MU_EARTH / r**3)
-        c, s = math.cos(n * dt), math.sin(n * dt)
-
-        # Linearized relative motion (CW equations)
-        phi = np.eye(6)
-        phi[0, 0] = 4 - 3*c
-        phi[0, 1] = 0
-        phi[0, 3] = s / n
-        phi[0, 4] = 2*(1 - c) / n
-        phi[1, 0] = 6*(s - n*dt)
-        phi[1, 1] = 1
-        phi[1, 3] = -2*(1 - c) / n
-        phi[1, 4] = (4*s - 3*n*dt) / n
-        phi[2, 2] = c
-        phi[2, 5] = s / n
-        phi[3, 0] = 3*n*s
-        phi[3, 3] = c
-        phi[3, 4] = 2*s
-        phi[4, 0] = -6*n*(1 - c)
-        phi[4, 3] = -2*s
-        phi[4, 4] = 4*c - 3
-        phi[5, 2] = -n*s
-        phi[5, 5] = c
-
-        # J2 secular correction to STM diagonal
-        a = r
-        inc_approx = math.acos(max(-1, min(1, state[2] / r)))
-        j2_rate = -1.5 * n * self.J2 * (self.R_EARTH / a)**2 * math.cos(inc_approx)
-        j2_correction = j2_rate * dt
-        phi[3, 0] += j2_correction * 0.1
-        phi[4, 1] += j2_correction * 0.05
-
-        return phi
+        # Integrate the absolute ECI state and its variational equations with
+        # RK4. A 60-second maximum step resolves typical Earth orbits while
+        # avoiding the invalid local-frame assumptions of Hill/CW dynamics.
+        augmented = np.concatenate((state, np.eye(6).ravel()))
+        step_count = max(1, math.ceil(abs(dt) / 60.0))
+        step = dt / step_count
+        for _ in range(step_count):
+            k1 = self._variational_derivative(augmented)
+            k2 = self._variational_derivative(augmented + 0.5 * step * k1)
+            k3 = self._variational_derivative(augmented + 0.5 * step * k2)
+            k4 = self._variational_derivative(augmented + step * k3)
+            augmented += step * (k1 + 2*k2 + 2*k3 + k4) / 6.0
+
+        return augmented[6:].reshape(6, 6)
 
     def compute_stt2(self, state: np.ndarray, dt: float) -> np.ndarray:
         """

output/builds/src/science/thermospheric_model.py

@@ -25,6 +25,8 @@
 import logging
 import numpy as np
 import math
+import os
+from pathlib import Path
 from typing import Tuple, List, Dict, Optional
 from dataclasses import dataclass, field
 
@@ -294,11 +296,39 @@ def forward(self, x):
                     h4 = self.dropout(h4)
                     return self.head(h4)
 
-            self.model = ResidualMLP(in_dim=8, hidden=hidden_dim)
-            self.model.eval()
-            self._has_torch = True
-            n_params = sum(p.numel() for p in self.model.parameters())
-            logger.info(f"DensityCorrectionNet initialized: {n_params:,} parameters")
+            configured_path = os.environ.get("DENSITY_MODEL_PATH")
+            weights_path = (
+                Path(configured_path).expanduser()
+                if configured_path
+                else Path(__file__).resolve().parents[4] / "models" / "density_mlp.pth"
+            )
+            if weights_path.is_file():
+                model = ResidualMLP(in_dim=8, hidden=hidden_dim)
+                try:
+                    state_dict = torch.load(
+                        weights_path, map_location="cpu", weights_only=True
+                    )
+                    model.load_state_dict(state_dict)
+                except Exception as exc:
+                    logger.warning(
+                        "Could not load density model weights from %s: %s; "
+                        "using empirical corrections",
+                        weights_path, exc,
+                    )
+                else:
+                    model.eval()
+                    self.model = model
+                    self._has_torch = True
+                    n_params = sum(p.numel() for p in model.parameters())
+                    logger.info(
+                        "DensityCorrectionNet loaded %s parameters from %s",
+                        f"{n_params:,}", weights_path,
+                    )
+            else:
+                logger.info(
+                    "Density model weights not found at %s — using empirical corrections only",
+                    weights_path,
+                )
         except ImportError:
             logger.info("PyTorch not available — using empirical corrections only")