fix: altitude-aware drag lifetime in cascade SIR, viewer connection error handling

- Add estimate_drag_lifetime_years() to atmosphere.py using exponential
  density model and representative debris ballistic coefficient
- Fix cascade_evolution_packets: use altitude-appropriate drag lifetime
  instead of hardcoded 25 yr, sensible debris seed (0.5% of sat count),
  R_EARTH_EQUATORIAL for altitude calculation
- Fix viewer_server do_POST: catch ConnectionAbortedError/BrokenPipeError
  when client disconnects during analysis generation
- Update README: new viewer screenshots, badge version/test count update

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: planet.jeroen

README.md

@@ -1,12 +1,12 @@
 # Humeris
 
-[![Version](https://img.shields.io/badge/version-1.28.1-blue.svg)](packages/core/pyproject.toml) [![Python](https://img.shields.io/badge/python-3.11_%7C_3.12_%7C_3.13-blue.svg)](packages/core/pyproject.toml) [![Tests](https://img.shields.io/badge/tests-3601_passing-brightgreen.svg)](tests/) [![License](https://img.shields.io/badge/license-MIT_(core)-green.svg)](LICENSE) [![License](https://img.shields.io/badge/license-Commercial_(pro)-red.svg)](COMMERCIAL-LICENSE.md) [![Architecture](https://img.shields.io/badge/architecture-hexagonal-purple.svg)](docs/architecture.md)
+[![Version](https://img.shields.io/badge/version-1.28.2-blue.svg)](packages/core/pyproject.toml) [![Python](https://img.shields.io/badge/python-3.11_%7C_3.12_%7C_3.13-blue.svg)](packages/core/pyproject.toml) [![Tests](https://img.shields.io/badge/tests-3653_passing-brightgreen.svg)](tests/) [![License](https://img.shields.io/badge/license-MIT_(core)-green.svg)](LICENSE) [![License](https://img.shields.io/badge/license-Commercial_(pro)-red.svg)](COMMERCIAL-LICENSE.md) [![Architecture](https://img.shields.io/badge/architecture-hexagonal-purple.svg)](docs/architecture.md)
 
 A Python library for satellite constellation analysis — from generating
 Walker shells to propagating orbits, screening conjunctions, and
 visualizing everything on a 3D globe. Pure Python, no compiled extensions.
 
-![Humeris interactive viewer](docs/interface.png)
+![Humeris constellation viewer — full orbital shell with ground tracks](docs/Pics/satellite_details_zoomed_out.png)
 
 ## What it does
 
@@ -22,6 +22,8 @@ visualizing everything on a 3D globe. Pure Python, no compiled extensions.
 - **Export everywhere** — CSV, GeoJSON, CZML, KML, Celestia, Blender, SpaceEngine, KSP, Universe Sandbox, Stellarium
 - **Interactive 3D viewer** — built-in CesiumJS viewer with 21 analysis layer types
 
+![Satellite detail view with orbital parameters and data table](docs/Pics/satelite_details.png)
+
 > **Disclaimer**: This library provides computational models for educational,
 > research, and engineering analysis purposes. It is not certified for
 > operational mission planning, safety-of-flight decisions, or regulatory
@@ -77,6 +79,10 @@ identity, and coordinate frame round-trips.
 The comparison artifacts — including cases where our results are imperfect
 or still evolving — are published alongside the code.
 
+![Kessler cascade debris cloud](docs/Pics/kessler_only.png)
+
+![Starlink conjunction screening — Kessler analysis](docs/Pics/starlink_kessler.png)
+
 Details: [Validation](docs/validation.md)
 
 ## Documentation

docs/Pics/kessler_only.png

@@ -47,7 +47,7 @@
 from humeris.domain.radiation import compute_l_shell
 from humeris.domain.eclipse import compute_beta_angle
 from humeris.domain.deorbit import assess_deorbit_compliance, DeorbitRegulation
-from humeris.domain.atmosphere import DragConfig
+from humeris.domain.atmosphere import DragConfig, estimate_drag_lifetime_years
 from humeris.domain.station_keeping import (
     compute_station_keeping_budget,
     StationKeepingConfig,
@@ -1864,15 +1864,19 @@ def cascade_evolution_packets(
     n_sats = len(states)
 
     # Shell parameters for SIR model
-    mean_alt_km = sum(s.semi_major_axis_m - OrbitalConstants.R_EARTH for s in states) / (n_sats * 1000.0)
-    shell_r = (OrbitalConstants.R_EARTH / 1000.0) + mean_alt_km
+    mean_alt_km = sum(s.semi_major_axis_m - OrbitalConstants.R_EARTH_EQUATORIAL for s in states) / (n_sats * 1000.0)
+    shell_r = (OrbitalConstants.R_EARTH_EQUATORIAL / 1000.0) + mean_alt_km
     shell_vol = 4.0 * math.pi * shell_r ** 2 * 50.0
+    drag_life = estimate_drag_lifetime_years(min(mean_alt_km, 900.0))
+    # Initial debris seed: 0.5% of satellite count (background debris)
+    debris_seed = max(10.0, n_sats * 0.005)
     duration_yr = max(0.1, total_seconds / (365.25 * 86400.0))
     step_yr = max(1e-6, duration_yr / max(num_steps, 1))
 
     sir = compute_cascade_sir(
-        shell_volume_km3=shell_vol, spatial_density_per_km3=n_sats / shell_vol,
+        shell_volume_km3=shell_vol, spatial_density_per_km3=debris_seed / shell_vol,
         mean_collision_velocity_ms=10000.0, satellite_count=n_sats,
+        drag_lifetime_years=drag_life,
         duration_years=duration_yr, step_years=step_yr,
     )
     sir_times, sir_s, sir_i = sir.time_series_years, sir.susceptible, sir.infected

docs/Pics/satelite_details.png

@@ -1682,9 +1682,14 @@ def _do_GET(self) -> None:
     def do_POST(self) -> None:
         try:
             self._do_POST()
+        except (ConnectionAbortedError, BrokenPipeError, ConnectionResetError):
+            logger.debug("Client disconnected during POST %s", self.path)
         except Exception:
             logger.exception("Unhandled error in POST %s", self.path)
-            self._error_response(500, "Internal server error")
+            try:
+                self._error_response(500, "Internal server error")
+            except (ConnectionAbortedError, BrokenPipeError, ConnectionResetError):
+                logger.debug("Client disconnected before error response")
 
     def _do_POST(self) -> None:
         base, param = self._route_path()

docs/Pics/satellite_details_zoomed_out.png

@@ -141,6 +141,65 @@ def atmospheric_density(
     return float(rho_base * np.exp(-(altitude_km - h_base) / scale_height))
 
 
+def estimate_drag_lifetime_years(
+    altitude_km: float,
+    model: AtmosphereModel = AtmosphereModel.VALLADO_4TH,
+) -> float:
+    """Estimate drag lifetime for debris at a given altitude.
+
+    Uses a representative debris ballistic coefficient (C_d=2.2, A=0.1 m²,
+    m=10 kg → B_c = 0.022 m²/kg) and the exponential atmosphere model to
+    estimate the time for an orbit to decay from the given altitude.
+
+    Approximation: lifetime ≈ H / (ρ * v * B_c * a) where H is the scale
+    height at the given altitude.
+
+    Args:
+        altitude_km: Altitude above Earth surface in km (100-2000).
+        model: Atmosphere model for density lookup.
+
+    Returns:
+        Estimated drag lifetime in years. Clamped to [0.01, 1000].
+    """
+    table = _MODEL_TABLES[model]
+    max_alt = table[-1][0]
+
+    # Above atmosphere table: essentially infinite lifetime
+    if altitude_km > max_alt:
+        return 1000.0
+
+    # Clamp to table minimum
+    alt = max(altitude_km, table[0][0])
+
+    # Representative debris: Cd=2.2, A=0.1 m², m=10 kg
+    b_c = 2.2 * 0.1 / 10.0  # 0.022 m²/kg
+
+    rho = atmospheric_density(alt, model=model)
+    a_m = OrbitalConstants.R_EARTH_EQUATORIAL + alt * 1000.0
+    v = float(np.sqrt(OrbitalConstants.MU_EARTH / a_m))
+
+    # Scale height at this altitude (interpolate from table)
+    lo, hi = 0, len(table) - 1
+    while lo < hi - 1:
+        mid = (lo + hi) // 2
+        if table[mid][0] <= alt:
+            lo = mid
+        else:
+            hi = mid
+    scale_height_km = table[lo][2]
+    scale_height_m = scale_height_km * 1000.0
+
+    # Lifetime ≈ H / (ρ * v * B_c * a)
+    da_dt = rho * v * b_c * a_m
+    if da_dt < 1e-30:
+        return 1000.0
+
+    lifetime_s = scale_height_m / da_dt
+    lifetime_yr = lifetime_s / (365.25 * 86400.0)
+
+    return max(0.01, min(lifetime_yr, 1000.0))
+
+
 def drag_acceleration(density: float, velocity: float, drag_config: DragConfig) -> float:
     """Drag acceleration magnitude.
 

docs/Pics/starlink_kessler.png

@@ -6,6 +6,7 @@
 import math
 
 from humeris.domain.cascade_analysis import CascadeSIR, compute_cascade_sir
+from humeris.domain.orbital_mechanics import OrbitalConstants
 
 
 class TestCascadeSIRSubcritical:
@@ -267,3 +268,116 @@ def test_time_series_length(self):
         assert len(result.susceptible) == expected_steps
         assert len(result.infected) == expected_steps
         assert len(result.recovered) == expected_steps
+
+
+# ── estimate_drag_lifetime_years ──────────────────────────────────
+
+class TestEstimateDragLifetime:
+    """Tests for altitude-dependent drag lifetime estimation."""
+
+    def test_200km_short_lifetime(self):
+        """At 200 km, drag lifetime is very short (< 2 years)."""
+        from humeris.domain.atmosphere import estimate_drag_lifetime_years
+
+        lifetime = estimate_drag_lifetime_years(200.0)
+        assert lifetime < 2.0
+
+    def test_550km_moderate_lifetime(self):
+        """At 550 km (Starlink altitude), drag lifetime is 5-20 years."""
+        from humeris.domain.atmosphere import estimate_drag_lifetime_years
+
+        lifetime = estimate_drag_lifetime_years(550.0)
+        assert 5.0 < lifetime < 20.0
+
+    def test_800km_long_lifetime(self):
+        """At 800 km, drag lifetime is > 30 years."""
+        from humeris.domain.atmosphere import estimate_drag_lifetime_years
+
+        lifetime = estimate_drag_lifetime_years(800.0)
+        assert lifetime > 30.0
+
+    def test_monotonic_increase(self):
+        """Drag lifetime increases monotonically with altitude."""
+        from humeris.domain.atmosphere import estimate_drag_lifetime_years
+
+        altitudes = [200, 300, 400, 500, 600, 700, 800, 900]
+        lifetimes = [estimate_drag_lifetime_years(a) for a in altitudes]
+        for i in range(len(lifetimes) - 1):
+            assert lifetimes[i] < lifetimes[i + 1], (
+                f"Lifetime decreased: {altitudes[i]} km={lifetimes[i]:.1f} yr "
+                f"> {altitudes[i+1]} km={lifetimes[i+1]:.1f} yr"
+            )
+
+    def test_meo_very_long(self):
+        """MEO altitudes (>1000 km) have very long drag lifetimes (>100 yr)."""
+        from humeris.domain.atmosphere import estimate_drag_lifetime_years
+
+        lifetime = estimate_drag_lifetime_years(1000.0)
+        assert lifetime > 100.0
+
+    def test_all_positive(self):
+        """Drag lifetime is always positive."""
+        from humeris.domain.atmosphere import estimate_drag_lifetime_years
+
+        for alt in [150, 200, 300, 500, 700, 900]:
+            assert estimate_drag_lifetime_years(alt) > 0.0
+
+
+# ── Cascade SIR R₀ with realistic parameters ─────────────────────
+
+class TestCascadeSIRRealisticR0:
+    """R₀ should be physically reasonable with altitude-appropriate drag lifetime."""
+
+    def test_starlink_r0_below_2(self):
+        """Starlink-like constellation at 550 km: R₀ should be < 2.0.
+
+        With ~6500 sats at 550 km, drag lifetime ~8-15 yr gives R₀ ~ 0.4-1.2.
+        The old hardcoded 25 yr gave R₀ ~ 1.7 — unrealistically high.
+        """
+        from humeris.domain.atmosphere import estimate_drag_lifetime_years
+
+        alt_km = 550.0
+        n_sats = 6500
+        shell_r = (OrbitalConstants.R_EARTH_EQUATORIAL / 1000.0) + alt_km
+        shell_vol = 4.0 * math.pi * shell_r ** 2 * 50.0
+        drag_life = estimate_drag_lifetime_years(alt_km)
+
+        result = compute_cascade_sir(
+            shell_volume_km3=shell_vol,
+            spatial_density_per_km3=0.005 * n_sats / shell_vol,
+            mean_collision_velocity_ms=10_000.0,
+            satellite_count=n_sats,
+            drag_lifetime_years=drag_life,
+            duration_years=1.0,
+        )
+        assert result.r_0 < 2.0, (
+            f"Starlink R₀={result.r_0:.2f} with drag_life={drag_life:.1f} yr "
+            f"— should be < 2.0"
+        )
+
+    def test_glonass_subcritical(self):
+        """GLONASS at ~19,100 km with 24 sats: R₀ should be << 1 (subcritical).
+
+        At MEO, drag is negligible, but with only 24 satellites the collision
+        rate is far too low for cascade.
+        """
+        from humeris.domain.atmosphere import estimate_drag_lifetime_years
+
+        alt_km = 19100.0
+        n_sats = 24
+        shell_r = (OrbitalConstants.R_EARTH_EQUATORIAL / 1000.0) + alt_km
+        shell_vol = 4.0 * math.pi * shell_r ** 2 * 200.0  # wider shell for MEO
+        # At MEO, clamp lifetime to table max
+        drag_life = min(estimate_drag_lifetime_years(min(alt_km, 900.0)), 500.0)
+
+        result = compute_cascade_sir(
+            shell_volume_km3=shell_vol,
+            spatial_density_per_km3=0.005 * n_sats / shell_vol,
+            mean_collision_velocity_ms=7_000.0,
+            satellite_count=n_sats,
+            drag_lifetime_years=drag_life,
+            duration_years=1.0,
+        )
+        assert result.r_0 < 0.1, (
+            f"GLONASS R₀={result.r_0:.4f} — should be << 1"
+        )