diff --git a/output/builds/src/conjunction_screener.py b/output/builds/src/conjunction_screener.py index 562b839..06a2ea9 100644 --- a/output/builds/src/conjunction_screener.py +++ b/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, diff --git a/output/builds/src/hardware/ascom_hal.py b/output/builds/src/hardware/ascom_hal.py index 1f12ecb..31675db 100644 --- a/output/builds/src/hardware/ascom_hal.py +++ b/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" diff --git a/output/builds/src/science/bayesian_iod.py b/output/builds/src/science/bayesian_iod.py index 6a698d9..e566778 100644 --- a/output/builds/src/science/bayesian_iod.py +++ b/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 diff --git a/output/builds/src/science/non_gaussian_pc.py b/output/builds/src/science/non_gaussian_pc.py index 5c6a5b2..8fbe695 100644 --- a/output/builds/src/science/non_gaussian_pc.py +++ b/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: """ diff --git a/output/builds/src/science/thermospheric_model.py b/output/builds/src/science/thermospheric_model.py index 360be9d..b34a7f7 100644 --- a/output/builds/src/science/thermospheric_model.py +++ b/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")