Add cooldown telemetry, RHR cooldown path, and shutdown regression test

.gitignore

@@ -1,3 +1,4 @@
 .venv/
 **/__pycache__/
 *.egg-info/
+artifacts/

AGENTS.md

@@ -2,14 +2,18 @@
 
 ## Project Structure & Module Organization
 All source code lives under `src/reactor_sim`. Submodules map to plant systems: `fuel.py` and `neutronics.py` govern fission power, `thermal.py` and `coolant.py` cover heat transfer and pumps, `turbine.py` drives the steam cycle, and `consumer.py` represents external electrical loads. High-level coordination happens in `reactor.py` and `simulation.py`. The convenience runner `run_simulation.py` executes the default scenario; add notebooks or scenario scripts under `experiments/` (create as needed). Keep assets such as plots or exported telemetry inside `artifacts/`.
+Feature catalog lives in `FEATURES.md`; update it whenever behavior is added or changed. Long-term realism tasks are tracked in `TODO.md`; keep it in sync as work progresses.
 
 ## Build, Test, and Development Commands
-- `python -m reactor_sim.simulation` — run the default 10-minute transient and print JSON snapshots.
-- `python run_simulation.py` — same as above, handy for IDE launch configs.
-- `python -m pip install -e .[dev]` — install editable package with optional tooling when dev dependencies are defined.
+- `.venv/bin/python -m reactor_sim.simulation` — run the default 10-minute transient and print JSON snapshots using the checked-in virtualenv.
+- `.venv/bin/python run_simulation.py` — same as above, handy for IDE launch configs.
+- `.venv/bin/python -m pip install -e .[dev]` — install editable package with optional tooling when dev dependencies are defined.
+- When running pytest or ad-hoc scripts, prefer `.venv/bin/python -m pytest` so everyone uses the shared interpreter and dependency set that ships with the repo.
+- State snapshots now auto-save/load from `artifacts/last_state.json`; override the path via `FISSION_STATE_PATH` (or continue to use `FISSION_SAVE_STATE`/`FISSION_LOAD_STATE` for explicit overrides).
+- A git remote for `origin` is configured; push changes to `origin/main` once work is complete so dashboards stay in sync.
 
 ## Operations & Control Hooks
-Manual commands live in `reactor_sim.commands.ReactorCommand`. Pass a `command_provider` callable to `ReactorSimulation` to adjust rods, pumps, turbine, coolant demand, or the attached `ElectricalConsumer`. Use `ReactorCommand.scram_all()` for full shutdown, `ReactorCommand(consumer_online=True, consumer_demand=600)` to hook the grid, or toggle pumps (`primary_pump_on=False`) to simulate faults. Control helpers in `control.py` expose `set_rods`, `increment_rods`, and `scram`.
+Manual commands live in `reactor_sim.commands.ReactorCommand`. Pass a `command_provider` callable to `ReactorSimulation` to adjust rods, pumps, turbine, coolant demand, or the attached `ElectricalConsumer`. Use `ReactorCommand.scram_all()` for full shutdown, `ReactorCommand(consumer_online=True, consumer_demand=600)` to hook the grid, or toggle pumps (`primary_pump_on=False`) to simulate faults. Control helpers in `control.py` expose `set_rods`, `increment_rods`, and `scram`, and you can switch `set_manual_mode(True)` to pause the automatic rod controller. For hands-on runs, launch the curses dashboard (`FISSION_DASHBOARD=1 FISSION_REALTIME=1 python run_simulation.py`) and use the on-screen shortcuts (q quit/save, space SCRAM, p/o pumps, t turbine, 1/2/3 toggle individual turbines, y/u/i turbine maintenance, m/n pump maintenance, k core maintenance, c consumer, r reset/clear saved state, + insert rods / - withdraw rods in 0.05 steps, [/ ] consumer demand, s/d setpoint, `a` toggles auto/manual rods). Recommended startup: enable manual rods (`a`), withdraw to ~0.3 before ramping the turbine/consumer, then re-enable auto control when you want closed-loop operation.
 The plant now boots cold (ambient core temperature, idle pumps); scripts must sequence startup: enable pumps, gradually withdraw rods, connect the consumer after turbine spin-up, and use `ControlSystem.set_power_setpoint` to chase desired output. Set `FISSION_REALTIME=1` to run continuously with real-time pacing; optionally set `FISSION_SIM_DURATION=infinite` for indefinite runs and send SIGINT/Ctrl+C to stop. Use `FISSION_SIM_DURATION=600` (default) for bounded offline batches.
 
 ## Coding Style & Naming Conventions
@@ -20,9 +24,13 @@ Pytest is the preferred framework. Place tests under `tests/` mirroring the `rea
 
 ## Commit & Pull Request Guidelines
 Write commits in imperative mood (`Add turbine moisture separator model`). Squash small WIP commits before review. Pull requests must describe the scenario, attach log excerpts (`snapshots.json` diff) or plots for novel behavior, and link to any tracking issues. Include validation notes: commands run, expected trends (e.g., outlet temperature increase), and outstanding risks.
+- Agents should automatically create a commit immediately after code/config/text changes are finished; use a clear imperative subject and avoid piling multiple logical changes into one commit unless the task explicitly says otherwise.
 
 ## Safety & Configuration
 Sim parameters live in constructors; never hard-code environment-specific paths. When adding new physics, guard unstable calculations with clamps and highlight operating limits in comments. Sensitive experiments (e.g., beyond design basis) should default to disabled flags so scripted studies remain safe by default.
 
 ## Reliability & Persistence
-Component wear is tracked via `failures.py`; stress from overheating, pump starvation, or turbine imbalance will degrade integrity and eventually disable the affected subsystem with automatic SCRAM for core damage. Persist plant snapshots with `ControlSystem.save_state()`/`load_state()` (used by `Reactor.save_state`/`load_state`) so long-running studies can resume; `FISSION_LOAD_STATE`/`FISSION_SAVE_STATE` env vars wire this into the CLI.
+Component wear is tracked via `failures.py`; stress from overheating, pump starvation, or turbine imbalance will degrade integrity and eventually disable the affected subsystem with automatic SCRAM for core damage. Plant state now persists between runs to `artifacts/last_state.json` by default (override via `FISSION_STATE_PATH` or the explicit save/load env vars). In the dashboard, press `r` to clear the saved snapshot and reboot the reactor to a cold, green-field state whenever you need a clean slate.
+
+## Session Context
+See `CONTEXT_NOTES.md` for the latest behavioral changes, controls, and assumptions carried over between sessions.

CONTEXT_NOTES.md

@@ -0,0 +1,24 @@
+# Session Context Notes
+
+- Reactor model: two-loop but tuned to RBMK-like pressures (nominal ~7 MPa for both loops). Loop pressure clamps to saturation baseline when pumps are off; pumps ramp flow/pressure over spool time when stopping or starting.
+- Feedwater & level: secondary steam drum uses shrink/swell-aware level sensing to drive a feedwater valve; makeup flow scales with steam draw toward the target level instead of instant inventory clamps.
+- Chemistry/fouling: plant tracks dissolved O2/boron/sodium; impurities plus high temp/steam draw increase HX fouling (reduces UA) and add condenser fouling/back-pressure. Oxygen degasses with steam; impurity ingress accelerates when venting.
+- Reactivity bias: boron ppm now biases shutdown reactivity; a very slow trim nudges boron toward the power setpoint after ~300s of operation.
+- Turbines: produce zero output unless steam quality is present and effective steam flow is >10 kg/s. Turbine panel shows steam availability (enthalpy × quality × mass flow) and steam enthalpy instead of loop pressure; condenser pressure/temperature/fouling shown with nominal bounds.
+- Generators: two diesel units, rated 50 MW, spool time 10s. Auto mode default `False`; manual toggles b/v. Auto stops when no load. Relief valves toggles l (primary) / ; (secondary) and displayed per loop.
+- Pumps: per-unit controls g/h (primary), j/k (secondary). Flow/pressure ramp down over spool when pumps stop. Pump status thresholds use >0.1 kg/s to show STOPPING.
+- Maintenance hotkeys: p (core, requires shutdown), m/n (primary 1/2), ,/. (secondary 1/2), B/V (generator 1/2), y/u/i (turbine 1/2/3).
+- Dashboard: two-column layout, trends panel for fuel temp and core power (delta and rate). Power Stats show aux demand/supply, generator and turbine output. Turbine panel shows steam availability/enthalpy instead of loop pressure. Core/temp/power lines include nominal/max.
+- Thermal updates: primary/secondary inlet temps now back-computed; when secondary flow is near zero, loops cool toward ambient over time. Relief venting removes mass/enthalpy with a multi-second ramp toward ~1 MPa and quenches superheat toward target-pressure saturation. Pumps cap target pressure when reliefs are open. Condenser modeled with vacuum pump drawdown, cooling-sink temperature, and fouling-driven back-pressure penalty.
+- Meltdown threshold: 2873 K. Auto rod control clears shutdown when set to auto and adjusts rods. Control rod worth/tuning currently unchanged.
+- Tests: `pytest` passing after all changes. Key regression additions include generator manual mode, turbine no-steam output, auto rod control, and passive cool-down.
+- Coolant demand fixed: demand increases when primary outlet is above target (sign was flipped before), so hot loops ramp flow instead of backing off.
+- Pump spin-down: pressure/flow ramp down over pump spool time with STOPPING->OFF threshold at 0.1 kg/s; prevents instant drop when last pump stops.
+- Steam pressure display shows 0 unless steam quality ≥0.05 and flow ≥100 kg/s to avoid showing pump head as steam pressure.
+- Passive cool-down: when secondary flow ~0, loops cool toward ambient; primary inlet/outlet back-propagated from transferred heat and ambient.
+- Relief valves: l (primary) and ; (secondary) vent with mass/enthalpy loss, ramp pressure toward ~1 MPa over several seconds, cap pump targets while open; status displayed per loop.
+- Generator behavior: starting/running only produce power when load is present; auto off by default; manual toggles b/v; auto stops with no load; base aux drops to 0 when idle/cold.
+- Pressure tying: loop pressure floors to saturation(temp) when pumps off; pump targets aim for ~7 MPa nominal RBMK-like setpoints.
+- Turbines/governor: require meaningful steam flow/quality; otherwise zero output. Steam supply pressure in turbine panel reads 0 when no steam. Throttle now biases toward load demand with a governor term, and overspeed/overload >105% rated electrical trips a unit and spools it down.
+- Rod control now supports three banks with weighted worth; xenon/iodine tracked with decay and burn-out; new UA·ΔT_lm steam-generator model and pump head/flow curves.
+- Dashboard shows heat-exchanger ΔT/efficiency and protections; pumps and HX changes documented in FEATURES.md / TODO.md.

FEATURES.md

@@ -0,0 +1,10 @@
+## C.O.R.E. feature set
+
+- **Core physics**: point-kinetics with per-bank delayed neutron precursors, temperature feedback, fuel burnup penalty, xenon/iodine buildup with decay and burn-out, and rod-bank worth curves.
+- **Rod control**: three rod banks with weighted worth; auto controller chases 3 GW setpoint with safety backoff and filtered power feedback; manual mode with staged bank motion and SCRAM; state persists across runs. Soluble boron bias contributes slow negative reactivity and trims toward the setpoint.
+- **Coolant & hydraulics**: primary/secondary pumps with head/flow curves, power draw scaling, wear tracking; pressure floors tied to saturation; auxiliary power model with generator auto-start.
+- **Heat transfer**: steam-generator UA·ΔT_lm model with a pinch cap to keep the primary outlet hotter than the secondary, coolant heating uses total fission power with fuel heating decoupled from exchanger draw, and the secondary thermal solver includes passive cool-down plus steam-drum mass/energy balance with latent heat and a shrink/swell-aware feedwater valve controller; dissolved oxygen/sodium drive HX fouling that reduces effective UA.
+- **Pressurizer & inventory**: primary pressurizer trims toward 7 MPa with level tracking, loop inventories/levels steer flow availability, secondary steam boil-off draws down level with auto makeup, and pumps reduce flow/status to `CAV` when NPSH is insufficient.
+- **Steam cycle**: three turbines with spool dynamics, throttle mapping and a simple governor with overspeed/overload trip, condenser vacuum/back-pressure with fouling and cooling sink temperature, chemistry-driven fouling/back-pressure penalties, load dispatch to consumer, steam quality gating for output, generator states with batteries/spool, and steam enthalpy-driven availability readout on the dashboard.
+- **Protections & failures**: health monitor degrading components under stress, automatic SCRAM on core or heat-sink loss plus DNB/subcool and secondary level/pressure trips, relief valves per loop with venting/mass loss and pump pressure caps, maintenance actions to restore integrity.
+- **Persistence & ops**: snapshots auto-save/load to `artifacts/last_state.json`; dashboard with live metrics, protections/warnings, heat-exchanger telemetry, component health, and control shortcuts.

TODO.md

@@ -0,0 +1,24 @@
+## Future realism upgrades
+
+- [x] Steam generator UA·ΔT_lm heat exchange and pump head/flow curves; keep validating temps under nominal load.
+- [x] Rod banks with worth curves, xenon/samarium buildup, and delayed-group kinetics per bank.
+- [x] Pressurizer behavior, primary/secondary inventory and level effects, and pump NPSH/cavitation checks.
+- [x] Model feedwater/steam-drum mass-energy balance, turbine throttle/efficiency maps, and condenser back-pressure.
+- [x] Introduce CHF/DNB margin, clad/fuel split temps, and SCRAM matrix for subcooling loss or SG level/pressure trips.
+- [x] Flesh out condenser behavior: vacuum pump limits, cooling water temperature coupling, and dynamic back-pressure with fouling.
+- [x] Dashboard polish: compact turbine/generator rows, color critical warnings (SCRAM/heat-sink), and reduce repeated log noise.
+- [ ] Dashboard multi-page view (F1/F2): retain numeric view on F1; future F2 schematic should mirror real PWR layout with ASCII art, flow/relief status, and minimal animations; add help/status hints and size checks; keep perf sane.
+- [x] Core thermal realism: add a simple radial fuel model (pellet/rim/clad temps) with burnup-driven conductivity drop and gap conductance; upgrade CHF/DNB correlation (e.g., Groeneveld/Chen) parameterized by pressure, mass flux, and quality, then calibrate margins.
+- [ ] Transient protection ladder: add dP/dt and dT/dt trips for SG overfill/depressurization and rod-run-in alarms; implement graded warn/arm/trip stages surfaced on the dashboard.
+- [x] Chemistry & fouling: track dissolved oxygen/boron and corrosion/fouling that degrade HX efficiency and condenser vacuum; let feedwater temperature/chemistry affect steam purity/back-pressure.
+- [x] Balance-of-plant dynamics: steam-drum level controller with shrink/swell, feedwater valve model, turbine throttle governor/overspeed trip, and improved load-follow tied to grid demand ramps.
+- [ ] Neutronics/feedback smoothing: add detector/measurement lag, fuel→clad→coolant transport delays, shared boron trim for fine regulation, and retune rod gains/rate limits to reduce power hunting while keeping safety margins intact.
+- [ ] Component wear & maintenance: make wear depend on duty cycle and off-nominal conditions (cavitation, high ΔP, high temp); add preventive maintenance scheduling and show next-due in the dashboard.
+- [ ] Scenarios & tooling: presets for cold start, load-follow, and fault injection (pump fail, relief stuck) with seedable randomness; snapshot diff tooling to compare saved states.
+- [ ] Incremental realism plan:
+  - [x] Add stored enthalpy for primary/secondary loops and a steam-drum mass/energy balance (sensible + latent) while keeping existing pump logic and tests passing. Target representative PWR conditions: primary 15–16 MPa, 290–320 °C inlet/320–330 °C outlet, secondary saturation ~6–7 MPa with boil at ~490–510 K.
+  - [x] Adjust HX/pressure handling to use stored energy (saturation clamp and pressure rise) and validate steam formation with both pumps at ~3 GW. Use realistic tube-side material assumptions (Inconel 690/SS cladding) and clamp steam quality to phase-equilibrium enthalpy.
+  - [x] Update turbine power mapping to consume steam enthalpy/quality and align protection trips with real steam presence; drive inlet steam around 6–7 MPa, quality/enthalpy-based flow to ~550–600 MW(e) per machine class if steam is available.
+  - [x] Add integration test: cold start → gens/pumps 2/2 → ramp to ~3 GW → confirm steam quality threshold at the secondary drum → enable all turbines and require electrical output. Include a step that tolerates one secondary pump off for a period to prove redundancy still yields steam.
+- [x] Dashboard follow-ups: replace turbine “Steam P” with a more useful steam availability signal (enthalpy × steam flow).
+- [x] Relief modeling: vent both loops gradually to ~1 MPa when reliefs are open, removing steam enthalpy/mass and capping pump targets to prevent instant repressurization.

pyproject.toml

@@ -10,6 +10,10 @@ dependencies = []
 
 [project.optional-dependencies]
 dev = ["pytest>=7.0"]
+dashboard = [
+    "rich>=13.7.0",
+    "textual>=0.50.0",
+]
 
 [build-system]
 requires = ["setuptools>=61"]

src/reactor_sim/__init__.py

@@ -5,6 +5,7 @@ from __future__ import annotations
 from .atomic import Atom, AtomicPhysics, FissionEvent
 from .commands import ReactorCommand
 from .consumer import ElectricalConsumer
+from .dashboard import ReactorDashboard
 from .failures import HealthMonitor
 from .logging_utils import configure_logging
 from .state import CoreState, CoolantLoopState, TurbineState, PlantState
@@ -24,5 +25,6 @@ __all__ = [
     "ReactorCommand",
     "ElectricalConsumer",
     "HealthMonitor",
+    "ReactorDashboard",
     "configure_logging",
 ]

src/reactor_sim/atomic.py

@@ -164,6 +164,14 @@ class FissionEvent:
     energy_mev: float
 
 
+@dataclass(frozen=True)
+class ReactionEvent:
+    parent: Atom
+    products: tuple[Atom, ...]
+    emitted_particles: tuple[str, ...]
+    energy_mev: float
+
+
 def make_atom(protons: int, neutrons: int) -> Atom:
     return Atom(symbol=element_symbol(protons), protons=protons, neutrons=neutrons)
 
@@ -235,6 +243,64 @@ class AtomicPhysics:
         )
         return event
 
+    def alpha_decay(self, atom: Atom) -> ReactionEvent:
+        if atom.protons <= 2 or atom.neutrons <= 2:
+            return ReactionEvent(parent=atom, products=(atom,), emitted_particles=(), energy_mev=0.0)
+        daughter = make_atom(atom.protons - 2, atom.neutrons - 2)
+        alpha = make_atom(2, 2)
+        energy = 5.0  # Typical alpha decay energy
+        return ReactionEvent(parent=atom, products=(daughter, alpha), emitted_particles=("alpha",), energy_mev=energy)
+
+    def beta_minus_decay(self, atom: Atom) -> ReactionEvent:
+        if atom.neutrons <= 0:
+            return ReactionEvent(parent=atom, products=(atom,), emitted_particles=(), energy_mev=0.0)
+        daughter = make_atom(atom.protons + 1, atom.neutrons - 1)
+        energy = 1.0
+        return ReactionEvent(
+            parent=atom,
+            products=(daughter,),
+            emitted_particles=("beta-", "anti-nu_e"),
+            energy_mev=energy,
+        )
+
+    def beta_plus_decay(self, atom: Atom) -> ReactionEvent:
+        if atom.protons <= 1:
+            return ReactionEvent(parent=atom, products=(atom,), emitted_particles=(), energy_mev=0.0)
+        daughter = make_atom(atom.protons - 1, atom.neutrons + 1)
+        energy = 1.0
+        return ReactionEvent(
+            parent=atom,
+            products=(daughter,),
+            emitted_particles=("beta+", "nu_e"),
+            energy_mev=energy,
+        )
+
+    def electron_capture(self, atom: Atom) -> ReactionEvent:
+        if atom.protons <= 1:
+            return ReactionEvent(parent=atom, products=(atom,), emitted_particles=(), energy_mev=0.0)
+        daughter = make_atom(atom.protons - 1, atom.neutrons + 1)
+        energy = 0.5
+        return ReactionEvent(
+            parent=atom,
+            products=(daughter,),
+            emitted_particles=("nu_e", "gamma"),
+            energy_mev=energy,
+        )
+
+    def gamma_emission(self, atom: Atom, energy_mev: float = 0.2) -> ReactionEvent:
+        return ReactionEvent(parent=atom, products=(atom,), emitted_particles=("gamma",), energy_mev=energy_mev)
+
+    def spontaneous_fission(self, atom: Atom) -> ReactionEvent:
+        # Reuse the generic splitter to keep mass/charge balanced.
+        split = self._generic_split(atom)
+        neutrons = ("n",) * max(0, split.emitted_neutrons)
+        return ReactionEvent(
+            parent=atom,
+            products=split.products,
+            emitted_particles=neutrons,
+            energy_mev=split.energy_mev,
+        )
+
     def _generic_split(self, atom: Atom) -> FissionEvent:
         heavy_mass = max(1, math.floor(atom.mass_number * 0.55))
         light_mass = atom.mass_number - heavy_mass

src/reactor_sim/commands.py

@@ -19,8 +19,22 @@ class ReactorCommand:
     power_setpoint: float | None = None
     consumer_online: bool | None = None
     consumer_demand: float | None = None
+    rod_manual: bool | None = None
+    turbine_units: dict[int, bool] | None = None
+    primary_pumps: dict[int, bool] | None = None
+    secondary_pumps: dict[int, bool] | None = None
+    generator_units: dict[int, bool] | None = None
+    generator_auto: bool | None = None
+    primary_relief: bool | None = None
+    secondary_relief: bool | None = None
+    maintenance_components: tuple[str, ...] = tuple()
 
     @classmethod
     def scram_all(cls) -> "ReactorCommand":
         """Convenience constructor for an emergency shutdown."""
         return cls(scram=True, turbine_on=False, primary_pump_on=True, secondary_pump_on=True)
+
+    @classmethod
+    def maintain(cls, *components: str) -> "ReactorCommand":
+        """Request maintenance for the provided component names."""
+        return cls(maintenance_components=tuple(components))

src/reactor_sim/constants.py

@@ -7,13 +7,91 @@ MEGAWATT = 1_000_000.0
 NEUTRON_LIFETIME = 0.1  # seconds, prompt neutron lifetime surrogate
 FUEL_ENERGY_DENSITY = 200.0 * MEGAWATT  # J/kg released as heat
 COOLANT_HEAT_CAPACITY = 4_200.0  # J/(kg*K) for water/steam
-COOLANT_DENSITY = 700.0  # kg/m^3 averaged between phases
-MAX_CORE_TEMPERATURE = 1_800.0  # K
-MAX_PRESSURE = 15.0  # MPa typical PWR primary loop limit
-CONTROL_ROD_SPEED = 0.05  # fraction insertion per second
+COOLANT_DENSITY = 720.0  # kg/m^3 averaged between phases
+STEAM_LATENT_HEAT = 2_200_000.0  # J/kg approximate latent heat of vaporization
+CORE_MELTDOWN_TEMPERATURE = 2_873.0  # K (approx 2600C) threshold for irreversible meltdown
+MAX_CORE_TEMPERATURE = CORE_MELTDOWN_TEMPERATURE  # Allow simulation to approach meltdown temperature
+MAX_PRESSURE = 16.0  # MPa PWR primary loop limit
+CLAD_MAX_TEMPERATURE = 1_200.0  # K clad softening / DNB concern
+CHF_MASS_FLUX_REF = 1_500.0  # kg/m2-s reference mass flux surrogate
+CHF_PRESSURE_REF_MPA = 7.0  # MPa reference pressure for CHF surrogate
+CONTROL_ROD_SPEED = 0.03  # fraction insertion per second
+CONTROL_ROD_WORTH = 0.042  # delta rho contribution when fully withdrawn
+CONTROL_ROD_BANK_WEIGHTS = (0.4, 0.35, 0.25)
+ROD_MANUAL_STEP = 0.025
 STEAM_TURBINE_EFFICIENCY = 0.34
 GENERATOR_EFFICIENCY = 0.96
 ENVIRONMENT_TEMPERATURE = 295.0  # K
 AMU_TO_KG = 1.660_539_066_60e-27
 MEV_TO_J = 1.602_176_634e-13
-ELECTRON_FISSION_CROSS_SECTION = 5e-23  # cm^2, tuned for simulation scale
+ELECTRON_FISSION_CROSS_SECTION = 5e-16  # cm^2, tuned for simulation scale
+PUMP_SPOOL_TIME = 5.0  # seconds to reach commanded flow
+PRIMARY_PUMP_SHUTOFF_HEAD_MPA = 17.0  # approximate shutoff head for primary pumps
+SECONDARY_PUMP_SHUTOFF_HEAD_MPA = 8.0
+TURBINE_SPOOL_TIME = 12.0  # seconds to reach steady output
+
+# Turbine/condenser parameters
+TURBINE_THROTTLE_MIN = 0.1
+TURBINE_THROTTLE_MAX = 1.0
+TURBINE_THROTTLE_EFFICIENCY_DROP = 0.15  # efficiency loss when at minimum throttle
+CONDENSER_BASE_PRESSURE_MPA = 0.01
+CONDENSER_MAX_PRESSURE_MPA = 0.3
+CONDENSER_BACKPRESSURE_PENALTY = 0.35  # fractional power loss at max back-pressure
+CONDENSER_VACUUM_PUMP_RATE = 0.05  # MPa per second drawdown toward base when below max load
+CONDENSER_COOLING_WATER_TEMP_K = 295.0  # cooling sink temperature
+CONDENSER_FOULING_RATE = 0.00002  # incremental penalty per second of hot operation
+CONDENSER_FOULING_MAX_PENALTY = 0.2  # max additional backpressure penalty from fouling
+CONDENSER_CHEM_FOULING_RATE = 0.0005  # per-second fouling increment scaled by impurity ppm
+CONDENSER_CHEM_BACKPRESSURE_FACTOR = 0.0002  # MPa increase per ppm impurities toward condenser pressure
+GENERATOR_SPOOL_TIME = 10.0  # seconds to reach full output
+# Auxiliary power assumptions
+PUMP_POWER_MW = 12.0  # MW draw per pump unit
+BASE_AUX_LOAD_MW = 5.0  # control, instrumentation, misc.
+NORMAL_CORE_POWER_MW = 3_000.0
+TEST_MAX_POWER_MW = 4_000.0
+PRIMARY_OUTLET_TARGET_K = 580.0
+SECONDARY_OUTLET_TARGET_K = 520.0
+PRIMARY_NOMINAL_PRESSURE = 15.5  # MPa PWR primary pressure
+SECONDARY_NOMINAL_PRESSURE = 6.5  # MPa steam drum/steam line pressure surrogate
+STEAM_GENERATOR_UA_MW_PER_K = 150.0  # overall UA for steam generator (MW/K)
+# Loop volume / inventory assumptions
+PRIMARY_LOOP_VOLUME_M3 = 350.0
+SECONDARY_LOOP_VOLUME_M3 = 320.0
+PRIMARY_PRESSURIZER_SETPOINT_MPA = 15.5
+PRIMARY_PRESSURIZER_DEADBAND_MPA = 0.2
+PRIMARY_PRESSURIZER_HEAT_RATE_MPA_PER_S = 0.08
+PRIMARY_PRESSURIZER_SPRAY_RATE_MPA_PER_S = 0.12
+PRIMARY_PRESSURIZER_LEVEL_DRAW_PER_S = 0.002
+PRIMARY_PRESSURIZER_LEVEL_FILL_PER_S = 0.001
+LOOP_INVENTORY_CORRECTION_RATE = 0.08  # fraction of nominal mass restored per second toward target
+PRIMARY_INVENTORY_TARGET = 0.95
+SECONDARY_INVENTORY_TARGET = 0.9
+SECONDARY_STEAM_LOSS_FRACTION = 0.02  # fraction of steam mass that leaves the loop each second
+NPSH_REQUIRED_MPA = 0.25
+LOW_LEVEL_FLOW_FLOOR = 0.05
+# Chemistry & fouling
+CHEM_MAX_PPM = 5_000.0
+CHEM_OXYGEN_DEFAULT_PPM = 50.0  # deoxygenated feedwater target (ppb -> ppm surrogate)
+CHEM_BORON_DEFAULT_PPM = 500.0
+CHEM_SODIUM_DEFAULT_PPM = 5.0
+HX_FOULING_RATE = 1e-5  # fouling increment per second scaled by impurities/temp
+HX_FOULING_HEAL_RATE = 5e-6  # cleaning/settling when cool/low steam
+HX_FOULING_MAX_PENALTY = 0.25  # fractional UA loss cap
+BORON_WORTH_PER_PPM = 8e-6  # delta rho per ppm relative to baseline boron
+BORON_TRIM_RATE_PPM_PER_S = 0.04  # slow boron trim toward setpoint when near target
+# Mild thermal/measurement lags
+FUEL_TO_CLAD_TIME_CONSTANT = 0.3  # seconds (mild lag)
+CLAD_TO_COOLANT_TIME_CONSTANT = 0.2  # seconds (mild lag)
+POWER_MEASUREMENT_TIME_CONSTANT = 10.0  # seconds
+# Passive cooldown
+PASSIVE_COOL_RATE_PRIMARY = 0.05  # K/s toward ambient when low/no transfer
+PASSIVE_COOL_RATE_SECONDARY = 0.08  # K/s toward ambient when no steam/heat sink
+RHR_ACTIVE = True
+RHR_CUTOFF_POWER_MW = 5.0
+RHR_COOL_RATE = 0.2  # K/s forced cooldown when power is near zero
+# Threshold inventories (event counts) for flagging common poisons in diagnostics.
+KEY_POISON_THRESHOLDS = {
+    "Xe": 1e20,  # xenon
+    "I": 1e20,  # iodine precursor to xenon
+    "Sm": 5e19,  # samarium
+}

src/reactor_sim/consumer.py

@@ -14,6 +14,7 @@ class ElectricalConsumer:
     demand_mw: float
     online: bool = False
     power_received_mw: float = 0.0
+    _under_supply_logged: bool = False
 
     def request_power(self) -> float:
         return self.demand_mw if self.online else 0.0
@@ -30,11 +31,16 @@ class ElectricalConsumer:
 
     def update_power_received(self, supplied_mw: float) -> None:
         self.power_received_mw = supplied_mw
-        if supplied_mw < self.request_power():
+        if supplied_mw + 1e-6 < self.request_power():
+            if not self._under_supply_logged:
                 LOGGER.warning(
                     "%s under-supplied: %.1f/%.1f MW",
                     self.name,
                     supplied_mw,
                     self.request_power(),
                 )
-
+                self._under_supply_logged = True
+        else:
+            if self._under_supply_logged:
+                LOGGER.info("%s demand satisfied", self.name)
+            self._under_supply_logged = False

src/reactor_sim/control.py

@@ -2,7 +2,7 @@
 
 from __future__ import annotations
 
-from dataclasses import dataclass
+from dataclasses import dataclass, field
 import json
 import logging
 from pathlib import Path
@@ -21,27 +21,104 @@ def clamp(value: float, lo: float, hi: float) -> float:
 class ControlSystem:
     setpoint_mw: float = 3_000.0
     rod_fraction: float = 0.5
+    manual_control: bool = False
+    rod_banks: list[float] = field(default_factory=lambda: [0.5, 0.5, 0.5])
+    rod_target: float = 0.5
+    _filtered_power_mw: float = 0.0
+    _integral_error: float = 0.0
+    _ramp_start_mw: float = 0.0
+    _ramp_progress_mw: float = 0.0
+    _last_manual: bool = False
 
     def update_rods(self, state: CoreState, dt: float) -> float:
-        error = (state.power_output_mw - self.setpoint_mw) / self.setpoint_mw
-        adjustment = -error * 0.3
-        adjustment = clamp(adjustment, -constants.CONTROL_ROD_SPEED * dt, constants.CONTROL_ROD_SPEED * dt)
-        previous = self.rod_fraction
-        self.rod_fraction = clamp(self.rod_fraction + adjustment, 0.0, 0.95)
-        LOGGER.debug("Control rods %.3f -> %.3f (error=%.3f)", previous, self.rod_fraction, error)
+        if not self.rod_banks or len(self.rod_banks) != len(constants.CONTROL_ROD_BANK_WEIGHTS):
+            self.rod_banks = [self.rod_fraction] * len(constants.CONTROL_ROD_BANK_WEIGHTS)
+        # Keep manual tweaks in sync with the target.
+        self.rod_target = clamp(self.rod_target, 0.0, 0.95)
+        if self.manual_control:
+            if abs(self.rod_fraction - self.effective_insertion()) > 1e-6:
+                self.rod_target = clamp(self.rod_fraction, 0.0, 0.95)
+            self._advance_banks(self.rod_target, dt)
+            self._last_manual = True
+            return self.rod_fraction
+
+        # On transition from manual -> auto, start a gentle setpoint ramp from current power.
+        if self._last_manual:
+            self._ramp_start_mw = max(0.0, state.power_output_mw)
+            self._ramp_progress_mw = 0.0
+        self._last_manual = False
+
+        # Setpoint ramp: close the gap at a limited MW/s.
+        ramp_rate = 5.0  # MW per second toward setpoint
+        effective_setpoint = self.setpoint_mw
+        if self._ramp_start_mw > 0.0 and self.setpoint_mw > self._ramp_start_mw:
+            self._ramp_progress_mw = min(
+                self.setpoint_mw - self._ramp_start_mw,
+                self._ramp_progress_mw + ramp_rate * dt,
+            )
+            effective_setpoint = self._ramp_start_mw + self._ramp_progress_mw
+        else:
+            effective_setpoint = self.setpoint_mw
+
+        raw_power = state.power_output_mw
+        # Begin filtering once we're in the vicinity of the setpoint to avoid chasing noise.
+        if self._filtered_power_mw <= 0.0:
+            self._filtered_power_mw = raw_power
+        if raw_power > 0.7 * self.setpoint_mw:
+            tau = constants.POWER_MEASUREMENT_TIME_CONSTANT
+            alpha = clamp(dt / max(1e-6, tau), 0.0, 1.0)
+            self._filtered_power_mw += alpha * (raw_power - self._filtered_power_mw)
+            measured_power = self._filtered_power_mw
+        else:
+            measured_power = raw_power
+
+        error = (measured_power - effective_setpoint) / max(1e-6, effective_setpoint)
+        near = measured_power > 0.9 * effective_setpoint
+        # Deadband near setpoint to prevent dithering; only apply when we're close to target.
+        if near and abs(error) < 0.01:
+            self._advance_banks(self.rod_target, dt)
+            return self.rod_fraction
+        # Integrate a bit to remove steady-state error.
+        self._integral_error = clamp(self._integral_error + error * dt, -0.05, 0.05)
+        p_gain = 0.35 if not near else 0.2
+        i_gain = 0.015 if not near else 0.01
+        adjustment = p_gain * error + i_gain * self._integral_error
+        speed = constants.CONTROL_ROD_SPEED * dt
+        # Hard high-band clamp: above 5% over setpoint, freeze withdrawals; above 10%, insert faster.
+        if measured_power > 1.1 * effective_setpoint:
+            speed *= 0.4
+            adjustment = -speed  # force insertion at capped rate
+        elif measured_power > 1.05 * effective_setpoint:
+            adjustment = min(adjustment, 0.0)
+            speed *= 0.6
+        # Soft clamp above setpoint: if we're >5% high, enforce insertion at capped rate.
+        if error > 0.1:
+            adjustment = min(adjustment, 0.0)
+            speed *= 0.5
+        elif error > 0.05:
+            adjustment = min(adjustment, adjustment)
+            speed *= 0.7
+        if near:
+            speed *= 0.5  # slow down when near setpoint to avoid overshoot
+        adjustment = clamp(adjustment, -speed, speed)
+        self.rod_target = clamp(self.rod_target + adjustment, 0.0, 0.95)
+        self._advance_banks(self.rod_target, dt)
+        LOGGER.debug("Control rod target=%.3f (error=%.3f)", self.rod_target, error)
         return self.rod_fraction
 
     def set_rods(self, fraction: float) -> float:
-        previous = self.rod_fraction
-        self.rod_fraction = clamp(fraction, 0.0, 0.95)
-        LOGGER.info("Manual rod set %.3f -> %.3f", previous, self.rod_fraction)
-        return self.rod_fraction
+        self.rod_target = self._quantize_manual(fraction)
+        self._advance_banks(self.rod_target, 0.0)
+        LOGGER.info("Manual rod target set to %.3f", self.rod_target)
+        return self.rod_target
 
     def increment_rods(self, delta: float) -> float:
         return self.set_rods(self.rod_fraction + delta)
 
     def scram(self) -> float:
-        self.rod_fraction = 0.95
+        self.rod_target = 0.95
+        self.rod_banks = [0.95 for _ in self.rod_banks]
+        self._sync_fraction()
         LOGGER.warning("SCRAM: rods fully inserted")
         return self.rod_fraction
 
@@ -50,13 +127,91 @@ class ControlSystem:
         self.setpoint_mw = clamp(megawatts, 100.0, 4_000.0)
         LOGGER.info("Power setpoint %.0f -> %.0f MW", previous, self.setpoint_mw)
 
-    def coolant_demand(self, primary: CoolantLoopState) -> float:
-        desired_temp = 580.0
-        error = (primary.temperature_out - desired_temp) / 100.0
-        demand = clamp(0.8 - error, 0.0, 1.0)
-        LOGGER.debug("Coolant demand %.2f for outlet %.1fK", demand, primary.temperature_out)
+    def set_manual_mode(self, manual: bool) -> None:
+        if self.manual_control != manual:
+            self.manual_control = manual
+        LOGGER.info("Rod control %s", "manual" if manual else "automatic")
+
+    def safety_backoff(self, subcooling_margin: float | None, dnb_margin: float | None, dt: float) -> None:
+        """Insert rods proactively when thermal margins are thin."""
+        if self.manual_control:
+            return
+        severity = 0.0
+        if subcooling_margin is not None:
+            severity = max(severity, max(0.0, 3.0 - subcooling_margin) / 3.0)
+        if dnb_margin is not None:
+            severity = max(severity, max(0.0, 0.7 - dnb_margin) / 0.7)
+        if severity <= 0.0:
+            return
+        backoff = (0.001 + 0.01 * severity) * dt
+        self.rod_target = clamp(self.rod_target + backoff, 0.0, 0.95)
+        self._advance_banks(self.rod_target, dt)
+        LOGGER.debug("Safety backoff applied: target=%.3f severity=%.2f", self.rod_target, severity)
+
+    def coolant_demand(
+        self,
+        primary: CoolantLoopState,
+        core_power_mw: float | None = None,
+        electrical_output_mw: float | None = None,
+    ) -> float:
+        desired_temp = constants.PRIMARY_OUTLET_TARGET_K
+        # Increase demand when outlet is hotter than desired, reduce when cooler.
+        temp_error = (primary.temperature_out - desired_temp) / 100.0
+        demand = 0.8 + temp_error
+        # Keep a light power-proportional floor so both pumps stay spinning without flooding the loop.
+        power_floor = 0.0
+        if core_power_mw is not None:
+            power_fraction = clamp(core_power_mw / constants.NORMAL_CORE_POWER_MW, 0.0, 1.5)
+            power_floor = 0.2 + 0.25 * power_fraction
+        demand = max(demand, power_floor)
+        # At power, keep primary pumps near full speed to preserve pressure/subcooling.
+        if core_power_mw is not None and core_power_mw > 500.0:
+            demand = max(demand, 0.8)
+        elif core_power_mw is not None and core_power_mw > 100.0:
+            demand = max(demand, 0.6)
+        demand = clamp(demand, 0.0, 1.0)
+        LOGGER.debug(
+            "Coolant demand %.2f (temp_error=%.2f, power_floor=%.2f) for outlet %.1fK power %.1f MW elec %.1f MW",
+            demand,
+            temp_error,
+            power_floor,
+            primary.temperature_out,
+            core_power_mw or 0.0,
+            electrical_output_mw or 0.0,
+        )
         return demand
 
+    def effective_insertion(self) -> float:
+        if not self.rod_banks:
+            return self.rod_fraction
+        weights = constants.CONTROL_ROD_BANK_WEIGHTS
+        total = sum(weights)
+        effective = sum(w * b for w, b in zip(weights, self.rod_banks)) / total
+        return clamp(effective, 0.0, 0.95)
+
+    def _advance_banks(self, target: float, dt: float) -> None:
+        speed = constants.CONTROL_ROD_SPEED * dt
+        new_banks: list[float] = []
+        for idx, pos in enumerate(self.rod_banks):
+            direction = 1 if target > pos else -1
+            step = direction * speed
+            updated = clamp(pos + step, 0.0, 0.95)
+            # Avoid overshoot
+            if (direction > 0 and updated > target) or (direction < 0 and updated < target):
+                updated = target
+            new_banks.append(updated)
+        self.rod_banks = new_banks
+        self._sync_fraction()
+
+    def _sync_fraction(self) -> None:
+        self.rod_fraction = self.effective_insertion()
+
+    def _quantize_manual(self, fraction: float) -> float:
+        step = constants.ROD_MANUAL_STEP
+        quantized = round(fraction / step) * step
+        return clamp(quantized, 0.0, 0.95)
+
+
     def save_state(
         self,
         filepath: str,
@@ -68,6 +223,9 @@ class ControlSystem:
             "control": {
                 "setpoint_mw": self.setpoint_mw,
                 "rod_fraction": self.rod_fraction,
+                "manual_control": self.manual_control,
+                "rod_banks": self.rod_banks,
+                "rod_target": self.rod_target,
             },
             "plant": plant_state.to_dict(),
             "metadata": metadata or {},
@@ -85,6 +243,10 @@ class ControlSystem:
         control = data.get("control", {})
         self.setpoint_mw = control.get("setpoint_mw", self.setpoint_mw)
         self.rod_fraction = control.get("rod_fraction", self.rod_fraction)
+        self.manual_control = control.get("manual_control", self.manual_control)
+        self.rod_banks = control.get("rod_banks", self.rod_banks) or self.rod_banks
+        self.rod_target = control.get("rod_target", self.rod_fraction)
+        self._sync_fraction()
         plant = PlantState.from_dict(data["plant"])
         LOGGER.info("Loaded plant state from %s", path)
         return plant, data.get("metadata", {}), data.get("health")

src/reactor_sim/coolant.py

@@ -5,6 +5,7 @@ from __future__ import annotations
 from dataclasses import dataclass
 import logging
 
+from . import constants
 from .state import CoolantLoopState
 
 LOGGER = logging.getLogger(__name__)
@@ -14,11 +15,22 @@ LOGGER = logging.getLogger(__name__)
 class Pump:
     nominal_flow: float
     efficiency: float = 0.9
+    shutoff_head_mpa: float = 6.0
+    spool_time: float = constants.PUMP_SPOOL_TIME
 
     def flow_rate(self, demand: float) -> float:
         demand = max(0.0, min(1.0, demand))
         return self.nominal_flow * (0.2 + 0.8 * demand) * self.efficiency
 
+    def performance(self, demand: float) -> tuple[float, float]:
+        """Return (flow_kg_s, head_mpa) at the given demand using a simple pump curve."""
+        demand = max(0.0, min(1.0, demand))
+        flow = self.flow_rate(demand)
+        flow_frac = flow / max(1e-3, self.nominal_flow)
+        # Keep a healthy head near nominal flow; fall off gently beyond the rated point.
+        head = self.shutoff_head_mpa * max(0.2, 1.0 - 0.4 * max(0.0, flow_frac))
+        return flow, head
+
     def step(self, loop: CoolantLoopState, demand: float) -> None:
         loop.mass_flow_rate = self.flow_rate(demand)
         loop.pressure = 12.0 * demand + 2.0

src/reactor_sim/dashboard.py

@@ -0,0 +1,941 @@
+"""Interactive curses dashboard for real-time reactor operation."""
+
+from __future__ import annotations
+
+import curses
+import logging
+from collections import deque
+from dataclasses import dataclass
+from pathlib import Path
+from typing import Optional
+
+from . import constants
+from .commands import ReactorCommand
+from .reactor import Reactor
+from .simulation import ReactorSimulation
+from .state import PlantState, PumpState
+
+LOGGER = logging.getLogger(__name__)
+
+
+def _build_numpad_mapping() -> dict[int, float]:
+    # Use keypad matrix constants when available; skip missing ones to avoid import errors on some terminals.
+    mapping: dict[int, float] = {}
+    table = {
+        "KEY_C1": 0.1,  # numpad 1
+        "KEY_C2": 0.2,  # numpad 2
+        "KEY_C3": 0.3,  # numpad 3
+        "KEY_B1": 0.4,  # numpad 4
+        "KEY_B2": 0.5,  # numpad 5
+        "KEY_B3": 0.6,  # numpad 6
+        "KEY_A1": 0.7,  # numpad 7
+        "KEY_A2": 0.8,  # numpad 8
+        "KEY_A3": 0.9,  # numpad 9
+        # Common keypad aliases when NumLock is on
+        "KEY_END": 0.1,
+        "KEY_DOWN": 0.2,
+        "KEY_NPAGE": 0.3,
+        "KEY_LEFT": 0.4,
+        "KEY_B2": 0.5,  # center stays 0.5
+        "KEY_RIGHT": 0.6,
+        "KEY_HOME": 0.7,
+        "KEY_UP": 0.8,
+        "KEY_PPAGE": 0.9,
+    }
+    for name, value in table.items():
+        code = getattr(curses, name, None)
+        if code is not None:
+            mapping[code] = value
+    return mapping
+
+
+_NUMPAD_ROD_KEYS = _build_numpad_mapping()
+
+
+@dataclass
+class DashboardKey:
+    key: str
+    description: str
+
+
+class ReactorDashboard:
+    """Minimal TUI dashboard allowing user control in real-time."""
+
+    def __init__(
+        self,
+        reactor: Reactor,
+        start_state: Optional[PlantState],
+        timestep: float = 1.0,
+        save_path: Optional[str] = None,
+    ) -> None:
+        self.reactor = reactor
+        self.start_state = start_state
+        self.timestep = timestep
+        self.save_path = save_path
+        self.pending_command: Optional[ReactorCommand] = None
+        self.sim: Optional[ReactorSimulation] = None
+        self.quit_requested = False
+        self.reset_requested = False
+        self.page = 1  # 1=metrics, 2=schematic (placeholder)
+        self._last_state: Optional[PlantState] = None
+        self._trend_history: deque[tuple[float, float, float]] = deque(maxlen=120)
+        self.log_buffer: deque[str] = deque(maxlen=8)
+        self._log_handler: Optional[logging.Handler] = None
+        self._previous_handlers: list[logging.Handler] = []
+        self._logger = logging.getLogger("reactor_sim")
+        self.help_sections: list[tuple[str, list[DashboardKey]]] = [
+            (
+                "Reactor / Safety",
+                [
+                    DashboardKey("q", "Quit & save"),
+                    DashboardKey("space", "SCRAM"),
+                    DashboardKey("r", "Reset & clear state"),
+                    DashboardKey("a", "Toggle auto rod control"),
+                    DashboardKey("F1/F2", "Metrics / schematic views"),
+                    DashboardKey("+/-", "Withdraw/insert rods"),
+                    DashboardKey("1-9 / Numpad", "Set rods to 0.1 … 0.9 (manual)"),
+                    DashboardKey("[/]", "Adjust consumer demand −/+50 MW"),
+                    DashboardKey("s/d", "Setpoint −/+250 MW"),
+                    DashboardKey("p", "Maintain core (shutdown required)"),
+                ],
+            ),
+            (
+                "Pumps",
+                [
+                    DashboardKey("g/h", "Toggle primary pump 1/2"),
+                    DashboardKey("j/k", "Toggle secondary pump 1/2"),
+                    DashboardKey("m/n", "Maintain primary pumps 1/2"),
+                    DashboardKey(",/.", "Maintain secondary pumps 1/2"),
+                ],
+            ),
+            (
+                "Generators",
+                [
+                    DashboardKey("b/v", "Toggle generator 1/2"),
+                    DashboardKey("x", "Toggle generator auto"),
+                    DashboardKey("l/;", "Toggle relief primary/secondary"),
+                    DashboardKey("B/V", "Maintain generator 1/2"),
+                ],
+            ),
+            (
+                "Turbines / Grid",
+                [
+                    DashboardKey("t", "Toggle turbine bank"),
+                    DashboardKey("Shift+1/2/3", "Toggle turbine units 1-3"),
+                    DashboardKey("y/u/i", "Maintain turbine 1/2/3"),
+                    DashboardKey("c", "Toggle consumer"),
+                ],
+            ),
+        ]
+
+    def run(self) -> None:
+        curses.wrapper(self._main)
+
+    def _main(self, stdscr: "curses._CursesWindow") -> None:  # type: ignore[name-defined]
+        curses.curs_set(0)
+        curses.start_color()
+        curses.use_default_colors()
+        curses.init_pair(1, curses.COLOR_CYAN, -1)
+        curses.init_pair(2, curses.COLOR_YELLOW, -1)
+        curses.init_pair(3, curses.COLOR_GREEN, -1)
+        curses.init_pair(4, curses.COLOR_RED, -1)
+        stdscr.nodelay(True)
+        stdscr.keypad(True)
+        self._install_log_capture()
+        try:
+            while True:
+                self.sim = ReactorSimulation(
+                    self.reactor,
+                    timestep=self.timestep,
+                    duration=None,
+                    realtime=True,
+                    command_provider=self._next_command,
+                )
+                self.sim.start_state = self.start_state
+                try:
+                    for state in self.sim.run():
+                        self._last_state = state
+                        self._draw(stdscr, state)
+                        self._handle_input(stdscr)
+                        if self.quit_requested or self.reset_requested:
+                            # Persist the latest state if we are exiting early.
+                            if self.sim:
+                                self.sim.last_state = state
+                            self.sim.stop()
+                            break
+                finally:
+                    if self.save_path and self.sim and self.sim.last_state and not self.reset_requested:
+                        self.reactor.save_state(self.save_path, self.sim.last_state)
+                if self.quit_requested:
+                    break
+                if self.reset_requested:
+                    self._clear_saved_state()
+                    self._reset_to_greenfield()
+                    continue
+                break
+        finally:
+            self._restore_logging()
+
+    def _handle_input(self, stdscr: "curses._CursesWindow") -> None:
+        while True:
+            ch = stdscr.getch()
+            if ch == -1:
+                break
+            keyname = None
+            try:
+                keyname = curses.keyname(ch)
+            except curses.error:
+                keyname = None
+            if ch in (ord("q"), ord("Q")):
+                self.quit_requested = True
+                return
+            if ch == ord(" "):
+                self._queue_command(ReactorCommand.scram_all())
+            elif ch in (ord("o"), ord("O")):
+                continue
+            elif ch in (ord("g"), ord("G")):
+                self._toggle_primary_pump_unit(0)
+            elif ch in (ord("h"), ord("H")):
+                self._toggle_primary_pump_unit(1)
+            elif ch in (ord("j"), ord("J")):
+                self._toggle_secondary_pump_unit(0)
+            elif ch in (ord("k"), ord("K")):
+                self._toggle_secondary_pump_unit(1)
+            elif ch in (ord("p"), ord("P")):
+                self._queue_command(ReactorCommand.maintain("core"))
+            elif ch in (ord("b"), ord("B")):
+                self._toggle_generator_unit(0)
+            elif ch in (ord("v"), ord("V")):
+                self._toggle_generator_unit(1)
+            elif ch == ord("l"):
+                self._queue_command(ReactorCommand(primary_relief=not self.reactor.primary_relief_open))
+            elif ch == ord(";"):
+                self._queue_command(ReactorCommand(secondary_relief=not self.reactor.secondary_relief_open))
+            elif ch in (ord("x"), ord("X")):
+                self._queue_command(ReactorCommand(generator_auto=not self.reactor.generator_auto))
+            elif ch in (ord("t"), ord("T")):
+                self._queue_command(ReactorCommand(turbine_on=not self.reactor.turbine_active))
+            elif keyname and keyname.decode(errors="ignore") in ("!", "@", "#", '"'):
+                name = keyname.decode(errors="ignore")
+                turbine_hotkeys = {"!": 0, "@": 1, "#": 2, '"': 1}
+                self._toggle_turbine_unit(turbine_hotkeys[name])
+            elif ch in (ord("!"), ord("@"), ord("#"), ord('"')):
+                turbine_hotkeys = {ord("!"): 0, ord("@"): 1, ord("#"): 2, ord('"'): 1}
+                self._toggle_turbine_unit(turbine_hotkeys[ch])
+            elif keyname and keyname.startswith(b"KP_") and keyname[-1:] in b"123456789":
+                target = (keyname[-1] - ord("0")) / 10.0  # type: ignore[arg-type]
+                self._queue_command(ReactorCommand(rod_position=target, rod_manual=True))
+            elif ord("1") <= ch <= ord("9"):
+                target = (ch - ord("0")) / 10.0
+                self._queue_command(ReactorCommand(rod_position=target, rod_manual=True))
+            elif ch in _NUMPAD_ROD_KEYS:
+                self._queue_command(ReactorCommand(rod_position=_NUMPAD_ROD_KEYS[ch], rod_manual=True))
+            elif curses.KEY_F1 <= ch <= curses.KEY_F9:
+                target = (ch - curses.KEY_F1 + 1) / 10.0
+                self._queue_command(ReactorCommand(rod_position=target, rod_manual=True))
+            elif ch in (ord("+"), ord("=")):
+                # Insert rods (increase fraction)
+                self._queue_command(ReactorCommand(rod_position=self._clamped_rod(constants.ROD_MANUAL_STEP)))
+            elif ch == ord("-"):
+                # Withdraw rods (decrease fraction)
+                self._queue_command(ReactorCommand(rod_position=self._clamped_rod(-constants.ROD_MANUAL_STEP)))
+            elif ch == ord("["):
+                demand = self._current_demand() - 50.0
+                self._queue_command(ReactorCommand(consumer_demand=max(0.0, demand)))
+            elif ch == ord("]"):
+                demand = self._current_demand() + 50.0
+                self._queue_command(ReactorCommand(consumer_demand=demand))
+            elif ch in (ord("c"), ord("C")):
+                online = not (self.reactor.consumer.online if self.reactor.consumer else False)
+                self._queue_command(ReactorCommand(consumer_online=online))
+            elif ch in (ord("r"), ord("R")):
+                self._request_reset()
+            elif ch in (ord("s"), ord("S")):
+                self._queue_command(ReactorCommand(power_setpoint=self.reactor.control.setpoint_mw - 250.0))
+            elif ch in (ord("d"), ord("D")):
+                self._queue_command(ReactorCommand(power_setpoint=self.reactor.control.setpoint_mw + 250.0))
+            elif ch in (ord("a"), ord("A")):
+                self._queue_command(ReactorCommand(rod_manual=not self.reactor.control.manual_control))
+            elif ch in (ord("m"), ord("M")):
+                self._queue_command(ReactorCommand.maintain("primary_pump_1"))
+            elif ch in (ord("n"), ord("N")):
+                self._queue_command(ReactorCommand.maintain("primary_pump_2"))
+            elif ch == ord(","):
+                self._queue_command(ReactorCommand.maintain("secondary_pump_1"))
+            elif ch == ord("."):
+                self._queue_command(ReactorCommand.maintain("secondary_pump_2"))
+            elif ch in (ord("B"),):
+                self._queue_command(ReactorCommand.maintain("generator_1"))
+            elif ch in (ord("V"),):
+                self._queue_command(ReactorCommand.maintain("generator_2"))
+            elif ch in (ord("y"), ord("Y")):
+                self._queue_command(ReactorCommand.maintain("turbine_1"))
+            elif ch in (ord("u"), ord("U")):
+                self._queue_command(ReactorCommand.maintain("turbine_2"))
+            elif ch in (ord("i"), ord("I")):
+                self._queue_command(ReactorCommand.maintain("turbine_3"))
+
+    def _queue_command(self, command: ReactorCommand) -> None:
+        if self.pending_command is None:
+            self.pending_command = command
+        else:
+            for field in command.__dataclass_fields__:  # type: ignore[attr-defined]
+                value = getattr(command, field)
+                if value is None or value is False:
+                    continue
+                if field == "rod_step" and getattr(self.pending_command, field) is not None:
+                    setattr(
+                        self.pending_command,
+                        field,
+                        getattr(self.pending_command, field) + value,
+                    )
+                else:
+                    setattr(self.pending_command, field, value)
+        if self.pending_command.rod_position is not None and self.pending_command.rod_manual is None:
+            self.pending_command.rod_manual = True
+
+    def _toggle_turbine_unit(self, index: int) -> None:
+        if index < 0 or index >= len(self.reactor.turbine_unit_active):
+            return
+        current = self.reactor.turbine_unit_active[index]
+        self._queue_command(ReactorCommand(turbine_units={index + 1: not current}))
+
+    def _toggle_primary_pump_unit(self, index: int) -> None:
+        if index < 0 or index >= len(self.reactor.primary_pump_units):
+            return
+        current = self.reactor.primary_pump_units[index]
+        self._queue_command(ReactorCommand(primary_pumps={index + 1: not current}))
+
+    def _toggle_secondary_pump_unit(self, index: int) -> None:
+        if index < 0 or index >= len(self.reactor.secondary_pump_units):
+            return
+        current = self.reactor.secondary_pump_units[index]
+        self._queue_command(ReactorCommand(secondary_pumps={index + 1: not current}))
+
+    def _toggle_generator_unit(self, index: int) -> None:
+        current = False
+        if self._last_state and index < len(self._last_state.generators):
+            gen = self._last_state.generators[index]
+            current = gen.running or gen.starting
+        self._queue_command(ReactorCommand(generator_units={index + 1: not current}))
+
+    def _request_reset(self) -> None:
+        self.reset_requested = True
+        if self.sim:
+            self.sim.stop()
+
+    def _clear_saved_state(self) -> None:
+        if not self.save_path:
+            return
+        path = Path(self.save_path)
+        try:
+            path.unlink()
+            LOGGER.info("Cleared saved state at %s", path)
+        except FileNotFoundError:
+            LOGGER.info("No saved state to clear at %s", path)
+
+    def _reset_to_greenfield(self) -> None:
+        LOGGER.info("Resetting reactor to initial state")
+        self.reactor = Reactor.default()
+        self.start_state = None
+        self.pending_command = None
+        self._last_state = None
+        self.reset_requested = False
+        self.log_buffer.clear()
+
+
+    def _next_command(self, _: float, __: PlantState) -> Optional[ReactorCommand]:
+        cmd = self.pending_command
+        self.pending_command = None
+        return cmd
+
+    def _draw(self, stdscr: "curses._CursesWindow", state: PlantState) -> None:
+        stdscr.erase()
+        height, width = stdscr.getmaxyx()
+        min_status = 6
+        if height < min_status + 12 or width < 70:
+            stdscr.addstr(
+                0,
+                0,
+                "Terminal too small; try >=70x16 or reduce font size.".ljust(width),
+                curses.color_pair(4),
+            )
+            stdscr.refresh()
+            return
+
+        log_rows = max(2, min(len(self.log_buffer) + 2, 8))
+        status_height = min(height - 1, max(min_status, min_status + log_rows))
+        data_height = max(1, height - status_height)
+        gap = 2
+        right_width = max(28, width // 3)
+        left_width = width - right_width - gap
+        if left_width < 50:
+            left_width = min(50, width - (18 + gap))
+            right_width = width - left_width - gap
+
+        data_height = max(1, data_height)
+        left_width = max(1, left_width)
+        right_width = max(1, right_width)
+
+        data_win = stdscr.derwin(data_height, left_width, 0, 0)
+        help_win = stdscr.derwin(data_height, right_width, 0, left_width + gap)
+        status_win = stdscr.derwin(status_height, width, data_height, 0)
+
+        self._draw_data_panel(data_win, state)
+        self._draw_help_panel(help_win)
+        self._draw_status_panel(status_win, state)
+        stdscr.refresh()
+
+    def _draw_data_panel(self, win: "curses._CursesWindow", state: PlantState) -> None:
+        win.erase()
+        win.box()
+        win.addstr(0, 2, " Plant Overview ", curses.color_pair(1) | curses.A_BOLD)
+        self._update_trends(state)
+        height, width = win.getmaxyx()
+        inner_height = height - 2
+        inner_width = width - 2
+        left_width = max(28, inner_width // 2)
+        right_width = inner_width - left_width
+        left_win = win.derwin(inner_height, left_width, 1, 1)
+        right_win = win.derwin(inner_height, right_width, 1, 1 + left_width)
+        for row in range(1, height - 1):
+            win.addch(row, 1 + left_width, curses.ACS_VLINE)
+            if left_width + 1 < width - 1:
+                win.addch(row, 1 + left_width + 1, curses.ACS_VLINE)
+
+        left_win.erase()
+        right_win.erase()
+
+        left_y = 0
+        left_y = self._draw_section(
+            left_win,
+            left_y,
+            "Core",
+            [
+                (
+                    "Fuel Temp",
+                    f"{state.core.fuel_temperature:6.1f} K (Max {constants.CORE_MELTDOWN_TEMPERATURE:4.0f})",
+                ),
+                (
+                    "Core Power",
+                    f"{state.core.power_output_mw:6.1f} MW (Nom {constants.NORMAL_CORE_POWER_MW:4.0f}/Max {constants.TEST_MAX_POWER_MW:4.0f})",
+                ),
+                ("Neutron Flux", f"{state.core.neutron_flux:10.2e}"),
+                ("Rods", f"{self.reactor.control.rod_fraction:.3f}"),
+                ("Rod Mode", "AUTO" if not self.reactor.control.manual_control else "MANUAL"),
+                ("Setpoint", f"{self.reactor.control.setpoint_mw:7.0f} MW"),
+                ("Reactivity", f"{state.core.reactivity_margin:+.4f}"),
+                ("Boron", f"{state.boron_ppm:7.1f} ppm"),
+                ("P(meas)", f"{self._measured_power(state):6.1f} MW"),
+            ],
+        )
+        left_y = self._draw_section(
+            left_win,
+            left_y,
+            "Trends",
+            self._trend_lines(state),
+        )
+        left_y = self._draw_section(left_win, left_y, "Key Poisons / Emitters", self._poison_lines(state))
+        left_y = self._draw_section(
+            left_win,
+            left_y,
+            "Primary Loop",
+            [
+                ("Pump1", self._pump_status(state.primary_pumps, 0)),
+                ("Pump2", self._pump_status(state.primary_pumps, 1)),
+                (
+                    "Flow",
+                    f"{state.primary_loop.mass_flow_rate:7.0f}/{self.reactor.primary_pump.nominal_flow * len(self.reactor.primary_pump_units):.0f} kg/s",
+                ),
+                ("Level", f"{state.primary_loop.level*100:6.1f}%"),
+                ("Inlet Temp", f"{state.primary_loop.temperature_in:7.1f} K"),
+                ("Outlet Temp", f"{state.primary_loop.temperature_out:7.1f} K (Target {constants.PRIMARY_OUTLET_TARGET_K:4.0f})"),
+                ("Pressure", f"{state.primary_loop.pressure:5.2f}/{constants.MAX_PRESSURE:4.1f} MPa"),
+                ("Pressurizer", f"{self.reactor.pressurizer_level*100:6.1f}% @ {constants.PRIMARY_PRESSURIZER_SETPOINT_MPA:4.1f} MPa"),
+                ("Relief", "OPEN" if self.reactor.primary_relief_open else "CLOSED"),
+            ],
+        )
+        self._draw_section(
+            left_win,
+            left_y,
+            "Secondary Loop",
+            [
+                ("Pump1", self._pump_status(state.secondary_pumps, 0)),
+                ("Pump2", self._pump_status(state.secondary_pumps, 1)),
+                (
+                    "Flow",
+                    f"{state.secondary_loop.mass_flow_rate:7.0f}/{self.reactor.secondary_pump.nominal_flow * len(self.reactor.secondary_pump_units):.0f} kg/s",
+                ),
+                ("Level", f"{state.secondary_loop.level*100:6.1f}% (Target {constants.SECONDARY_INVENTORY_TARGET*100:4.0f}%)"),
+                (
+                    "Feedwater",
+                    f"valve {self.reactor.feedwater_valve*100:5.1f}% steam {state.secondary_loop.mass_flow_rate * max(0.0, state.secondary_loop.steam_quality):6.0f} kg/s",
+                ),
+                ("Inlet Temp", f"{state.secondary_loop.temperature_in:7.1f} K (Target {constants.SECONDARY_OUTLET_TARGET_K:4.0f})"),
+                ("Outlet Temp", f"{state.secondary_loop.temperature_out:7.1f} K (Target {constants.SECONDARY_OUTLET_TARGET_K:4.0f})"),
+                ("Pressure", f"{state.secondary_loop.pressure:5.2f}/{constants.MAX_PRESSURE:4.1f} MPa"),
+                ("Steam Quality", f"{state.secondary_loop.steam_quality:5.2f}/1.00"),
+                ("Relief", "OPEN" if self.reactor.secondary_relief_open else "CLOSED"),
+            ],
+        )
+
+        right_y = 0
+        consumer_status = "n/a"
+        consumer_demand = 0.0
+        if self.reactor.consumer:
+            consumer_status = "ONLINE" if self.reactor.consumer.online else "OFF"
+            consumer_demand = self.reactor.consumer.demand_mw
+        right_y = self._draw_section(
+            right_win,
+            right_y,
+            "Turbine / Grid",
+            [
+                ("Turbines", " ".join(self._turbine_status_lines())),
+                ("Rated Elec", f"{len(self.reactor.turbines)*self.reactor.turbines[0].rated_output_mw:7.1f} MW"),
+                (
+                    "Steam",
+                    f"h={state.turbines[0].steam_enthalpy:5.0f} kJ/kg avail {self._steam_available_power(state):6.1f} MW "
+                    f"flow {state.secondary_loop.mass_flow_rate * max(0.0, state.secondary_loop.steam_quality):6.0f} kg/s"
+                    if state.turbines
+                    else "n/a",
+                ),
+                (
+                    "Units Elec",
+                    " ".join([f"{t.electrical_output_mw:6.1f}MW" for t in state.turbines]) if state.turbines else "n/a",
+                ),
+                (
+                    "Governor",
+                    (
+                        f"thr {self.reactor.turbines[0].throttle:4.2f}→{self._desired_throttle(state.turbines[0]):4.2f} "
+                        f"ΔP {(state.turbines[0].load_demand_mw - state.turbines[0].electrical_output_mw):6.1f} MW"
+                    )
+                    if state.turbines
+                    else "n/a",
+                ),
+                (
+                    "Condenser",
+                    (
+                        f"P={state.turbines[0].condenser_pressure:4.2f}/{constants.CONDENSER_MAX_PRESSURE_MPA:4.2f} MPa "
+                        f"T={state.turbines[0].condenser_temperature:6.1f}K Foul={state.turbines[0].fouling_penalty*100:4.1f}%"
+                    )
+                    if state.turbines
+                    else "n/a",
+                ),
+                ("Electrical", f"{state.total_electrical_output():7.1f} MW | Load {self._total_load_supplied(state):6.1f}/{self._total_load_demand(state):6.1f} MW"),
+                ("Consumer", f"{consumer_status} demand {consumer_demand:6.1f} MW"),
+            ],
+        )
+        right_y = self._draw_section(right_win, right_y, "Generators", self._generator_lines(state))
+        right_y = self._draw_section(right_win, right_y, "Power Stats", self._power_lines(state))
+        right_y = self._draw_section(right_win, right_y, "Heat Exchanger", self._heat_exchanger_lines(state))
+        right_y = self._draw_section(right_win, right_y, "Protections / Warnings", self._protection_lines(state))
+        right_y = self._draw_section(right_win, right_y, "Maintenance", self._maintenance_lines())
+        self._draw_health_bars(right_win, right_y)
+
+    def _draw_help_panel(self, win: "curses._CursesWindow") -> None:
+        def _add_safe(row: int, col: int, text: str, attr: int = 0) -> bool:
+            max_y, max_x = win.getmaxyx()
+            if row >= max_y - 1 or col >= max_x - 1:
+                return False
+            clipped = text[: max(0, max_x - col - 1)]
+            try:
+                win.addstr(row, col, clipped, attr)
+            except curses.error:
+                return False
+            return True
+
+        win.erase()
+        win.box()
+        _add_safe(0, 2, " Controls ", curses.color_pair(1) | curses.A_BOLD)
+        y = 2
+        for title, entries in self.help_sections:
+            if not _add_safe(y, 2, title, curses.color_pair(1) | curses.A_BOLD):
+                return
+            y += 1
+            for entry in entries:
+                if not _add_safe(y, 4, f"{entry.key:<8} {entry.description}"):
+                    return
+                y += 1
+            y += 1
+        if not _add_safe(y, 2, "Tips:", curses.color_pair(2) | curses.A_BOLD):
+            return
+        tips = [
+            "Start pumps before withdrawing rods.",
+            "Bring turbine and consumer online after thermal stabilization.",
+            "Toggle turbine units (1/2/3) for staggered maintenance.",
+            "Use m/n/,/. for pump maintenance; B/V for generators.",
+            "Press 'r' to reset/clear state if you want a cold start.",
+            "Watch component health, DNB margin, and subcooling to avoid automatic trips.",
+        ]
+        for idx, tip in enumerate(tips, start=y + 2):
+            if not _add_safe(idx, 4, f"- {tip}"):
+                break
+
+    def _draw_status_panel(self, win: "curses._CursesWindow", state: PlantState) -> None:
+        win.erase()
+        win.hline(0, 0, curses.ACS_HLINE, win.getmaxyx()[1])
+        turbine_text = " ".join(self._turbine_status_lines())
+        msg = (
+            f"Time {state.time_elapsed:7.1f}s | Rods {self.reactor.control.rod_fraction:.3f} | "
+            f"Primary {'ON' if self.reactor.primary_pump_active else 'OFF'} | "
+            f"Secondary {'ON' if self.reactor.secondary_pump_active else 'OFF'} | "
+            f"Turbines {turbine_text} | Page {'Metrics' if self.page == 1 else 'Schematic'}"
+        )
+        win.addstr(1, 1, msg, curses.color_pair(3))
+        if self.reactor.health_monitor.failure_log:
+            win.addstr(
+                3,
+                1,
+                f"Failures: {', '.join(self.reactor.health_monitor.failure_log)}",
+                curses.color_pair(4) | curses.A_BOLD,
+            )
+            log_y = 4
+        else:
+            log_y = 3
+        for record in list(self.log_buffer):
+            if log_y >= win.getmaxyx()[0] - 1:
+                break
+            win.addstr(log_y, 1, record[: win.getmaxyx()[1] - 2], curses.color_pair(2))
+            log_y += 1
+        win.addstr(log_y, 1, "Press 'q' to exit and persist the current snapshot.", curses.color_pair(2))
+
+    def _draw_section(
+        self,
+        win: "curses._CursesWindow",
+        start_y: int,
+        title: str,
+        lines: list[tuple[str, str] | tuple[str, str, int] | str],
+    ) -> int:
+        height, width = win.getmaxyx()
+        inner_width = width - 4
+        if start_y >= height - 2:
+            return height - 2
+        win.addstr(start_y, 2, title[:inner_width], curses.A_BOLD | curses.color_pair(1))
+        row = start_y + 1
+        for line in lines:
+            if row >= height - 1:
+                break
+            attr = 0
+            if isinstance(line, tuple):
+                if len(line) == 3:
+                    label, value, attr = line
+                else:
+                    label, value = line
+                text = f"{label:<18}: {value}"
+            else:
+                text = line
+            win.addstr(row, 4, text[:inner_width], attr)
+            row += 1
+        return row + 1
+
+    def _flow_arrow(self, flow: float) -> str:
+        if flow > 15000:
+            return "====>"
+        if flow > 5000:
+            return "===>"
+        if flow > 500:
+            return "->"
+        return "--"
+
+    def _pump_glyph(self, pump_state: PumpState | None) -> str:
+        if pump_state is None:
+            return "·"
+        status = getattr(pump_state, "status", "OFF")
+        if status == "RUN":
+            return "▶"
+        if status == "CAV":
+            return "!"
+        if status == "STARTING":
+            return ">"
+        if status == "STOPPING":
+            return "-"
+        return "·"
+
+
+    def _turbine_status_lines(self) -> list[str]:
+        if not self.reactor.turbine_unit_active:
+            return ["n/a"]
+        lines: list[str] = []
+        for idx, active in enumerate(self.reactor.turbine_unit_active):
+            label = f"{idx + 1}:"
+            status = "ON" if active else "OFF"
+            if idx < len(getattr(self._last_state, "turbines", [])):
+                t_state = self._last_state.turbines[idx]
+                status = getattr(t_state, "status", status)
+            lines.append(f"{label}{status}")
+        return lines
+
+    def _total_load_supplied(self, state: PlantState) -> float:
+        return sum(t.load_supplied_mw for t in state.turbines)
+
+    def _total_load_demand(self, state: PlantState) -> float:
+        return sum(t.load_demand_mw for t in state.turbines)
+
+    def _poison_lines(self, state: PlantState) -> list[tuple[str, str]]:
+        inventory = state.core.fission_product_inventory or {}
+        particles = state.core.emitted_particles or {}
+        lines: list[tuple[str, str]] = []
+        def fmt(symbol: str, label: str, qty: float) -> tuple[str, str]:
+            threshold = constants.KEY_POISON_THRESHOLDS.get(symbol)
+            flag = " !" if threshold is not None and qty >= threshold else ""
+            return (f"{label}{flag}", f"{qty:9.2e}")
+
+        lines.append(fmt("Xe", "Xe (xenon)", getattr(state.core, "xenon_inventory", 0.0)))
+        lines.append(fmt("Sm", "Sm (samarium)", inventory.get("Sm", 0.0)))
+        lines.append(fmt("I", "I (iodine)", getattr(state.core, "iodine_inventory", 0.0)))
+        try:
+            xe_penalty = -self.reactor.neutronics.xenon_penalty(state.core)
+            lines.append(("Xe Δρ", f"{xe_penalty:+.4f}"))
+        except Exception:
+            pass
+        lines.append(("Neutrons (src)", f"{particles.get('n', 0.0):9.2e}"))
+        lines.append(("Gammas", f"{particles.get('gamma', 0.0):9.2e}"))
+        lines.append(("Alphas", f"{particles.get('alpha', 0.0):9.2e}"))
+        return lines
+
+    def _health_lines(self) -> list[tuple[str, str]]:
+        comps = self.reactor.health_monitor.components
+        lines: list[tuple[str, str]] = []
+        for name, comp in comps.items():
+            pct = f"{comp.integrity*100:5.1f}%"
+            state = "FAILED" if comp.failed else pct
+            lines.append((name, state))
+        return lines
+
+    def _maintenance_lines(self) -> list[tuple[str, str]]:
+        if not self.reactor.maintenance_active:
+            return [("Active", "None")]
+        return [(comp, "IN PROGRESS") for comp in sorted(self.reactor.maintenance_active)]
+
+    def _generator_lines(self, state: PlantState) -> list[tuple[str, str]]:
+        if not state.generators:
+            return [("Status", "n/a")]
+        lines: list[tuple[str, str]] = []
+        control = "AUTO" if self.reactor.generator_auto else "MANUAL"
+        lines.append(("Control", control))
+        for idx, gen in enumerate(state.generators):
+            status = "RUN" if gen.running else "START" if gen.starting else "OFF"
+            spool = f" spool {gen.spool_remaining:4.1f}s" if gen.starting else ""
+            lines.append((f"Gen{idx + 1}", f"{status} {gen.power_output_mw:6.1f}/{self.reactor.generators[idx].rated_output_mw:4.0f} MW{spool}"))
+            lines.append((f"  Battery", f"{gen.battery_charge*100:5.1f}% out {gen.battery_output_mw:4.1f} MW"))
+        return lines
+
+    def _power_lines(self, state: PlantState) -> list[tuple[str, str]]:
+        draws = getattr(state, "aux_draws", {}) or {}
+        primary_nom = constants.PUMP_POWER_MW * len(self.reactor.primary_pump_units)
+        secondary_nom = constants.PUMP_POWER_MW * len(self.reactor.secondary_pump_units)
+        lines = [
+            ("Base Aux", f"{draws.get('base', 0.0):6.1f}/{constants.BASE_AUX_LOAD_MW:4.1f} MW"),
+            ("Primary Aux", f"{draws.get('primary_pumps', 0.0):6.1f}/{primary_nom:4.1f} MW"),
+            ("Secondary Aux", f"{draws.get('secondary_pumps', 0.0):6.1f}/{secondary_nom:4.1f} MW"),
+            ("Aux Demand", f"{draws.get('total_demand', 0.0):6.1f} MW"),
+            ("Aux Supplied", f"{draws.get('supplied', 0.0):6.1f} MW"),
+            ("Gen Output", f"{draws.get('generator_output', 0.0):6.1f} MW"),
+            ("Turbine Elec", f"{draws.get('turbine_output', 0.0):6.1f} MW"),
+        ]
+        return lines
+
+    def _heat_exchanger_lines(self, state: PlantState) -> list[tuple[str, str]]:
+        delta_t = getattr(state, "primary_to_secondary_delta_t", 0.0)
+        eff = getattr(state, "heat_exchanger_efficiency", 0.0)
+        hx_fouling = getattr(state, "hx_fouling", 0.0)
+        return [
+            ("ΔT (pri-sec)", f"{delta_t:6.1f} K"),
+            ("HX Eff", f"{eff*100:6.1f}%"),
+            ("Chem/Foul", f"O2 {getattr(state, 'dissolved_oxygen_ppm', 0.0):5.1f} ppm Na {getattr(state, 'sodium_ppm', 0.0):5.1f} ppm Foul {hx_fouling*100:5.1f}%"),
+        ]
+
+    def _protection_lines(self, state: PlantState) -> list[tuple[str, str]]:
+        lines: list[tuple[str, str] | tuple[str, str, int]] = []
+        lines.append(("SCRAM", "ACTIVE" if self.reactor.shutdown else "CLEAR", curses.color_pair(4) if self.reactor.shutdown else 0))
+        if self.reactor.meltdown:
+            lines.append(("Meltdown", "IN PROGRESS", curses.color_pair(4) | curses.A_BOLD))
+        cooldown_status = "Normal"
+        if state.core.power_output_mw <= constants.RHR_CUTOFF_POWER_MW and (self.reactor.primary_pump_active or self.reactor.secondary_pump_active):
+            cooldown_status = "RHR/Passive"
+        lines.append(("Cooldown", cooldown_status))
+        sec_flow_low = state.secondary_loop.mass_flow_rate <= 1.0 or not self.reactor.secondary_pump_active
+        heat_sink_risk = sec_flow_low and state.core.power_output_mw > 50.0
+        if heat_sink_risk:
+            heat_text = "TRIPPED low secondary flow >50 MW"
+            heat_attr = curses.color_pair(4) | curses.A_BOLD
+        elif sec_flow_low:
+            heat_text = "ARMED (secondary off/low flow)"
+            heat_attr = curses.color_pair(2) | curses.A_BOLD
+        else:
+            heat_text = "OK"
+            heat_attr = curses.color_pair(3)
+        lines.append(("Heat sink", heat_text, heat_attr))
+
+        draws = getattr(state, "aux_draws", {}) or {}
+        demand = draws.get("total_demand", 0.0)
+        supplied = draws.get("supplied", 0.0)
+        if demand > 0.1 and supplied + 1e-6 < demand:
+            aux_text = f"DEFICIT {supplied:5.1f}/{demand:5.1f} MW"
+            aux_attr = curses.color_pair(2) | curses.A_BOLD
+        elif demand > 0.1:
+            aux_text = f"OK {supplied:5.1f}/{demand:5.1f} MW"
+            aux_attr = curses.color_pair(3)
+        else:
+            aux_text = "Idle"
+            aux_attr = 0
+        lines.append(("Aux power", aux_text, aux_attr))
+
+        reliefs = []
+        if self.reactor.primary_relief_open:
+            reliefs.append("Primary")
+        if self.reactor.secondary_relief_open:
+            reliefs.append("Secondary")
+        relief_attr = curses.color_pair(2) | curses.A_BOLD if reliefs else 0
+        lines.append(("Relief valves", ", ".join(reliefs) if reliefs else "Closed", relief_attr))
+        lines.append(("DNB margin", f"{state.core.dnb_margin:5.2f}" if state.core.dnb_margin is not None else "n/a"))
+        lines.append(("Subcooling", f"{state.core.subcooling_margin:5.1f} K" if state.core.subcooling_margin is not None else "n/a"))
+        lines.append(
+            (
+                "SG level",
+                f"{state.secondary_loop.level*100:5.1f}%",
+            )
+        )
+        lines.append(("SG pressure", f"{state.secondary_loop.pressure:5.2f}/{constants.MAX_PRESSURE:4.1f} MPa"))
+        lines.append(
+            (
+                "SCRAM trips",
+                "DNB<0.5 | Subcool<2K | SG lvl<5/>98% | SG P>15.2MPa",
+            )
+        )
+        return lines
+
+    def _steam_available_power(self, state: PlantState) -> float:
+        mass_flow = state.secondary_loop.mass_flow_rate * max(0.0, state.secondary_loop.steam_quality)
+        if mass_flow <= 1.0:
+            return 0.0
+        if state.turbines:
+            enthalpy_kjkg = max(0.0, state.turbines[0].steam_enthalpy)
+        else:
+            enthalpy_kjkg = (constants.STEAM_LATENT_HEAT / 1_000.0)
+        return (enthalpy_kjkg * mass_flow) / 1_000.0
+
+    def _measured_power(self, state: PlantState) -> float:
+        ctl = getattr(self, "reactor", None)
+        if ctl and getattr(ctl, "control", None):
+            filtered = getattr(ctl.control, "_filtered_power_mw", 0.0)
+            if filtered > 0.0:
+                return filtered
+        return state.core.power_output_mw
+
+    def _desired_throttle(self, turbine_state) -> float:
+        if not self.reactor.turbines:
+            return 0.0
+        turbine = self.reactor.turbines[0]
+        demand = turbine_state.load_demand_mw
+        return 0.4 if demand <= 0 else min(1.0, 0.4 + demand / max(1e-6, turbine.rated_output_mw))
+
+    def _update_trends(self, state: PlantState) -> None:
+        self._trend_history.append((state.time_elapsed, state.core.fuel_temperature, state.core.power_output_mw))
+
+    def _trend_lines(self, state: PlantState) -> list[tuple[str, str]]:
+        if len(self._trend_history) < 2:
+            return [("Fuel Temp Δ", "n/a"), ("Core Power Δ", "n/a"), ("P(meas)", "n/a")]
+        start_t, start_temp, start_power = self._trend_history[0]
+        end_t, end_temp, end_power = self._trend_history[-1]
+        duration = max(1.0, end_t - start_t)
+        temp_delta = end_temp - start_temp
+        power_delta = end_power - start_power
+        temp_rate = temp_delta / duration
+        power_rate = power_delta / duration
+        measured = getattr(self.reactor.control, "_filtered_power_mw", 0.0) if hasattr(self.reactor, "control") else 0.0
+        return [
+            ("Fuel Temp Δ", f"{end_temp:7.1f} K (Δ{temp_delta:+6.1f} / {duration:4.0f}s, {temp_rate:+5.2f}/s)"),
+            ("Core Power Δ", f"{end_power:7.1f} MW (Δ{power_delta:+6.1f} / {duration:4.0f}s, {power_rate:+5.2f}/s)"),
+            ("P(meas)", f"{measured:7.1f} MW" if measured > 0 else "n/a"),
+        ]
+
+    def _draw_health_bars(self, win: "curses._CursesWindow", start_y: int) -> int:
+        height, width = win.getmaxyx()
+        inner_width = width - 4
+        if start_y >= height - 2:
+            return height - 2
+        win.addstr(start_y, 2, "Component Health", curses.A_BOLD | curses.color_pair(1))
+        bar_width = max(8, min(inner_width - 18, 40))
+        label_width = 16
+        row = start_y + 1
+        for name, comp in self.reactor.health_monitor.components.items():
+            if row >= height - 1:
+                break
+            label = f"{name:<{label_width}}"
+            target = 0.0 if comp.failed else comp.integrity
+            filled = int(bar_width * max(0.0, min(1.0, target)))
+            bar = "#" * filled + "-" * (bar_width - filled)
+            color = 3 if comp.integrity > 0.5 else 2 if comp.integrity > 0.2 else 4
+            win.addstr(row, 4, f"{label}:")
+            bar_start = 4 + label_width + 1
+            win.addstr(row, bar_start, bar[:bar_width], curses.color_pair(color))
+            percent_text = "FAILED" if comp.failed else f"{comp.integrity*100:5.1f}%"
+            percent_x = min(width - len(percent_text) - 2, bar_start + bar_width + 2)
+            win.addstr(row, percent_x, percent_text, curses.color_pair(color))
+            row += 1
+        return row + 1
+
+    def _pump_status(self, pumps: list, index: int) -> str:
+        if index >= len(pumps):
+            return "n/a"
+        state = pumps[index]
+        status = getattr(state, "status", "ON" if state.active else "OFF")
+        return f"{status:<8} {state.flow_rate:6.0f} kg/s"
+
+    def _current_demand(self) -> float:
+        if self.reactor.consumer:
+            return self.reactor.consumer.demand_mw
+        return 0.0
+
+    def _clamped_rod(self, delta: float) -> float:
+        new_fraction = self.reactor.control.rod_fraction + delta
+        step = constants.ROD_MANUAL_STEP
+        quantized = round(new_fraction / step) * step
+        return max(0.0, min(0.95, quantized))
+
+    def _install_log_capture(self) -> None:
+        if self._log_handler:
+            return
+        self._previous_handlers = list(self._logger.handlers)
+        for handler in self._previous_handlers:
+            self._logger.removeHandler(handler)
+        handler = _DashboardLogHandler(self.log_buffer)
+        handler.setFormatter(logging.Formatter("%(levelname)s | %(message)s"))
+        self._logger.addHandler(handler)
+        self._logger.propagate = False
+        self._log_handler = handler
+
+    def _restore_logging(self) -> None:
+        if not self._log_handler:
+            return
+        self._logger.removeHandler(self._log_handler)
+        for handler in self._previous_handlers:
+            self._logger.addHandler(handler)
+        self._log_handler = None
+        self._previous_handlers = []
+
+
+class _DashboardLogHandler(logging.Handler):
+    def __init__(self, buffer: deque[str]) -> None:
+        super().__init__()
+        self.buffer = buffer
+        self._last_msg: str | None = None
+        self._repeat_count: int = 0
+
+    def emit(self, record: logging.LogRecord) -> None:
+        msg = self.format(record)
+        if msg == self._last_msg:
+            self._repeat_count += 1
+            if self.buffer and self.buffer[-1].startswith(self._last_msg):
+                try:
+                    self.buffer[-1] = f"{self._last_msg} (x{self._repeat_count + 1})"
+                except Exception:
+                    self.buffer.append(f"{msg} (x{self._repeat_count + 1})")
+            else:
+                self.buffer.append(f"{msg} (x{self._repeat_count + 1})")
+            return
+        else:
+            self._last_msg = msg
+            self._repeat_count = 0
+        self.buffer.append(msg)

src/reactor_sim/failures.py

@@ -30,6 +30,15 @@ class ComponentHealth:
             self.failed = True
             LOGGER.error("Component %s has failed", self.name)
 
+    def restore(self, amount: float) -> None:
+        if amount <= 0.0:
+            return
+        previous = self.integrity
+        self.integrity = min(1.0, self.integrity + amount)
+        if self.integrity > 0.05:
+            self.failed = False
+        LOGGER.info("Maintenance on %s: %.3f -> %.3f", self.name, previous, self.integrity)
+
     def snapshot(self) -> dict:
         return asdict(self)
 
@@ -41,13 +50,20 @@ class ComponentHealth:
 class HealthMonitor:
     """Tracks component wear and signals failures."""
 
-    def __init__(self) -> None:
+    def __init__(self, disable_degradation: bool = False) -> None:
+        self.disable_degradation = disable_degradation
         self.components: Dict[str, ComponentHealth] = {
             "core": ComponentHealth("core"),
-            "primary_pump": ComponentHealth("primary_pump"),
-            "secondary_pump": ComponentHealth("secondary_pump"),
-            "turbine": ComponentHealth("turbine"),
+            "primary_pump_1": ComponentHealth("primary_pump_1"),
+            "primary_pump_2": ComponentHealth("primary_pump_2"),
+            "secondary_pump_1": ComponentHealth("secondary_pump_1"),
+            "secondary_pump_2": ComponentHealth("secondary_pump_2"),
+            "generator_1": ComponentHealth("generator_1"),
+            "generator_2": ComponentHealth("generator_2"),
         }
+        for idx in range(3):
+            name = f"turbine_{idx + 1}"
+            self.components[name] = ComponentHealth(name)
         self.failure_log: list[str] = []
 
     def component(self, name: str) -> ComponentHealth:
@@ -56,43 +72,76 @@ class HealthMonitor:
     def evaluate(
         self,
         state: PlantState,
-        primary_active: bool,
-        secondary_active: bool,
-        turbine_active: bool,
+        primary_units: Iterable[bool],
+        secondary_units: Iterable[bool],
+        turbine_active: Iterable[bool],
+        generator_states: Iterable,
         dt: float,
     ) -> List[str]:
+        if self.disable_degradation:
+            return []
         events: list[str] = []
+        turbine_flags = list(turbine_active)
         core = self.component("core")
         core_temp = state.core.fuel_temperature
-        temp_stress = max(0.0, (core_temp - 900.0) / (constants.MAX_CORE_TEMPERATURE - 900.0))
+        temp_stress = max(0.0, (core_temp - 900.0) / max(1e-6, (constants.MAX_CORE_TEMPERATURE - 900.0)))
         base_degrade = 0.0001 * dt
         core.degrade(base_degrade + temp_stress * 0.01 * dt)
 
-        if primary_active:
-            primary_flow = state.primary_loop.mass_flow_rate
-            flow_ratio = 0.0 if primary_flow <= 0 else min(1.0, primary_flow / 18_000.0)
-            self.component("primary_pump").degrade((0.0002 + (1 - flow_ratio) * 0.005) * dt)
+        prim_units = list(primary_units)
+        sec_units = list(secondary_units)
+        prim_states = state.primary_pumps or []
+        sec_states = state.secondary_pumps or []
+        primary_temp = getattr(state.primary_loop, "temperature_out", 295.0)
+        secondary_temp = getattr(state.secondary_loop, "temperature_out", 295.0)
+        primary_heat_factor = max(0.0, (primary_temp - 600.0) / 400.0)  # harsher wear only when very hot
+        secondary_heat_factor = max(0.0, (secondary_temp - 520.0) / 300.0)
+        for idx, active in enumerate(prim_units):
+            comp = self.component(f"primary_pump_{idx + 1}")
+            if idx < len(prim_states) and active:
+                flow = prim_states[idx].flow_rate
+                flow_ratio = 0.0 if flow <= 0 else min(1.0, flow / 9_000.0)
+                low_flow_penalty = (1 - flow_ratio) * 0.001
+                rate = (0.0001 + low_flow_penalty) * (1 + primary_heat_factor)
+                comp.degrade(rate * dt)
             else:
-            self.component("primary_pump").degrade(0.0005 * dt)
+                comp.degrade(0.0)
+        for idx, active in enumerate(sec_units):
+            comp = self.component(f"secondary_pump_{idx + 1}")
+            if idx < len(sec_states) and active:
+                flow = sec_states[idx].flow_rate
+                flow_ratio = 0.0 if flow <= 0 else min(1.0, flow / 8_000.0)
+                low_flow_penalty = (1 - flow_ratio) * 0.0008
+                rate = (0.0001 + low_flow_penalty) * (1 + secondary_heat_factor)
+                comp.degrade(rate * dt)
+            else:
+                comp.degrade(0.0)
 
-        if secondary_active:
-            secondary_flow = state.secondary_loop.mass_flow_rate
-            flow_ratio = 0.0 if secondary_flow <= 0 else min(1.0, secondary_flow / 16_000.0)
-            self.component("secondary_pump").degrade((0.0002 + (1 - flow_ratio) * 0.004) * dt)
+        for idx, gen_state in enumerate(generator_states):
+            comp = self.component(f"generator_{idx + 1}")
+            running = getattr(gen_state, "running", False) or getattr(gen_state, "starting", False)
+            if running:
+                comp.degrade(0.00015 * dt)
             else:
-            self.component("secondary_pump").degrade(0.0005 * dt)
+                comp.degrade(0.0)
 
-        if turbine_active:
-            electrical = state.turbine.electrical_output_mw
-            load_ratio = (
-                0.0
-                if state.turbine.load_demand_mw <= 0
-                else min(1.0, electrical / max(1e-6, state.turbine.load_demand_mw))
-            )
-            stress = 0.0002 + abs(1 - load_ratio) * 0.003
-            self.component("turbine").degrade(stress * dt)
+        turbines = state.turbines if hasattr(state, "turbines") else []
+        for idx, active in enumerate(turbine_flags):
+            name = f"turbine_{idx + 1}"
+            component = self.component(name)
+            if active and idx < len(turbines):
+                turbine_state = turbines[idx]
+                demand = turbine_state.load_demand_mw
+                supplied = turbine_state.load_supplied_mw
+                if demand <= 0:
+                    load_ratio = 0.0
                 else:
-            self.component("turbine").degrade(0.0001 * dt)
+                    load_ratio = min(1.0, supplied / max(1e-6, demand))
+                mismatch = abs(1 - load_ratio)
+                stress = 0.00005 + mismatch * 0.0006
+            else:
+                stress = 0.00002
+            component.degrade(stress * dt)
 
         for name, component in self.components.items():
             if component.failed and name not in self.failure_log:
@@ -105,6 +154,19 @@ class HealthMonitor:
 
     def load_snapshot(self, data: dict) -> None:
         for name, comp_data in data.items():
-            if name in self.components:
-                self.components[name] = ComponentHealth.from_snapshot(comp_data)
+            mapped = "turbine_1" if name == "turbine" else name
+            if mapped in self.components:
+                self.components[mapped] = ComponentHealth.from_snapshot(comp_data)
+            elif mapped == "primary_pump":
+                self.components["primary_pump_1"] = ComponentHealth.from_snapshot(comp_data)
+                self.components["primary_pump_2"] = ComponentHealth.from_snapshot(comp_data)
+            elif mapped == "secondary_pump":
+                self.components["secondary_pump_1"] = ComponentHealth.from_snapshot(comp_data)
 
+    def maintain(self, component: str, amount: float = 0.05) -> bool:
+        comp = self.components.get(component)
+        if not comp:
+            LOGGER.warning("Maintenance requested for unknown component %s", component)
+            return False
+        comp.restore(amount)
+        return True

src/reactor_sim/fuel.py

@@ -29,6 +29,7 @@ class FuelAssembly:
     mass_kg: float
     fissile_atom: Atom = field(default_factory=lambda: make_atom(92, 143))
     atomic_physics: AtomicPhysics = field(default_factory=AtomicPhysics)
+    decay_activity_base: float = 1e16  # nominal decay events per second at burnup 0
 
     def available_energy_j(self, state: CoreState) -> float:
         fraction_remaining = max(0.0, 1.0 - state.burnup)
@@ -42,7 +43,10 @@ class FuelAssembly:
         event = self.simulate_electron_hit()
         effective_flux = max(0.0, flux * max(0.0, 1.0 - control_fraction))
         atoms = self.mass_kg / self.fissile_atom.atomic_mass_kg
-        event_rate = effective_flux * constants.ELECTRON_FISSION_CROSS_SECTION * atoms * self.enrichment
+        event_rate = max(
+            0.0,
+            effective_flux * constants.ELECTRON_FISSION_CROSS_SECTION * atoms * self.enrichment,
+        )
         power_watts = event_rate * event.energy_mev * constants.MEV_TO_J
         power_mw = power_watts / constants.MEGAWATT
         LOGGER.debug(
@@ -53,4 +57,53 @@ class FuelAssembly:
             event.products[1].mass_number,
             power_mw,
         )
-        return max(0.0, power_mw), event
+        return max(0.0, power_mw), event_rate, event
+
+    def decay_reaction_effects(
+        self, state: CoreState
+    ) -> tuple[float, float, dict[str, float], dict[str, float]]:
+        """Return (power MW, neutron source rate, product rates, particle rates)."""
+        # Scale activity with burnup (fission products) and a small flux contribution.
+        activity = self.decay_activity_base * (0.05 + state.burnup) * (1.0 + 0.00001 * state.neutron_flux)
+        activity = max(0.0, min(activity, 1e20))
+
+        reactions = [
+            self.atomic_physics.alpha_decay(self.fissile_atom),
+            self.atomic_physics.beta_minus_decay(self.fissile_atom),
+            self.atomic_physics.beta_plus_decay(self.fissile_atom),
+            self.atomic_physics.electron_capture(self.fissile_atom),
+            self.atomic_physics.gamma_emission(self.fissile_atom),
+            self.atomic_physics.spontaneous_fission(self.fissile_atom),
+        ]
+        weights = [1.0, 1.2, 0.5, 0.4, 1.6, 0.05]
+        total_weight = sum(weights)
+
+        power_watts = 0.0
+        neutron_source = 0.0
+        product_rates: dict[str, float] = {}
+        particle_rates: dict[str, float] = {}
+
+        for event, weight in zip(reactions, weights):
+            rate = activity * (weight / total_weight)
+            power_watts += rate * event.energy_mev * constants.MEV_TO_J
+            for atom in event.products:
+                product_rates[atom.symbol] = product_rates.get(atom.symbol, 0.0) + rate
+            for particle in event.emitted_particles:
+                particle_rates[particle] = particle_rates.get(particle, 0.0) + rate
+                if particle == "n":
+                    neutron_source += rate
+
+        power_mw = power_watts / constants.MEGAWATT
+        capped_power = min(power_mw, 0.1 * max(state.power_output_mw, 1.0))
+        LOGGER.debug(
+            "Decay reactions contribute %.2f MW (raw %.2f MW) with neutron source %.2e 1/s",
+            capped_power,
+            power_mw,
+            neutron_source,
+        )
+        return max(0.0, capped_power), neutron_source, product_rates, particle_rates
+
+    def decay_reaction_power(self, state: CoreState) -> float:
+        """Compatibility shim: return only power contribution."""
+        power_mw, _, _, _ = self.decay_reaction_effects(state)
+        return power_mw

src/reactor_sim/generator.py

@@ -0,0 +1,81 @@
+"""Auxiliary diesel generator model with spool dynamics."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+import logging
+
+from . import constants
+
+LOGGER = logging.getLogger(__name__)
+
+
+@dataclass
+class GeneratorState:
+    running: bool
+    starting: bool
+    spool_remaining: float
+    power_output_mw: float
+    battery_charge: float
+    status: str = "OFF"
+    battery_output_mw: float = 0.0
+
+
+@dataclass
+class DieselGenerator:
+    rated_output_mw: float = 50.0
+    spool_time: float = constants.GENERATOR_SPOOL_TIME
+
+    def start(self, state: GeneratorState) -> None:
+        if state.running or state.starting:
+            return
+        if state.battery_charge <= 0.05:
+            LOGGER.warning("Generator start failed: insufficient battery")
+            return
+        state.starting = True
+        state.spool_remaining = self.spool_time
+        state.status = "STARTING"
+        LOGGER.info("Generator starting (spool %.0fs)", self.spool_time)
+
+    def stop(self, state: GeneratorState) -> None:
+        if not (state.running or state.starting):
+            return
+        state.running = False
+        state.starting = False
+        state.spool_remaining = 0.0
+        state.power_output_mw = 0.0
+        state.status = "OFF"
+        LOGGER.info("Generator stopped")
+
+    def step(self, state: GeneratorState, load_demand_mw: float, dt: float) -> float:
+        """Advance generator dynamics and return delivered power."""
+        state.battery_output_mw = 0.0
+        if state.starting:
+            state.spool_remaining = max(0.0, state.spool_remaining - dt)
+            progress = 1.0 - state.spool_remaining / max(self.spool_time, 1e-6)
+            available = self.rated_output_mw * progress
+            delivered = min(available, max(0.0, load_demand_mw))
+            state.power_output_mw = delivered
+            state.battery_output_mw = min(0.5, delivered) if delivered > 0 else 0.0
+            if state.spool_remaining <= 0.0:
+                state.starting = False
+                state.running = True
+                state.status = "RUN"
+                LOGGER.info("Generator online at %.1f MW", self.rated_output_mw)
+        elif state.running:
+            available = self.rated_output_mw
+            delivered = min(available, max(0.0, load_demand_mw))
+            state.power_output_mw = delivered
+            state.status = "RUN" if state.power_output_mw > 0 else "IDLE"
+        else:
+            state.power_output_mw = 0.0
+            state.status = "OFF"
+
+        if state.running:
+            state.battery_charge = min(1.0, state.battery_charge + 0.02 * dt)
+        elif state.starting:
+            state.battery_charge = max(0.0, state.battery_charge - 0.003 * dt)
+        else:
+            state.battery_charge = max(0.0, state.battery_charge - 0.00005 * dt)
+
+        return state.power_output_mw

src/reactor_sim/neutronics.py

@@ -7,7 +7,7 @@ import logging
 
 from . import constants
 from .fuel import fuel_reactivity_penalty
-from .state import CoreState
+from .state import CoreState, clamp
 
 LOGGER = logging.getLogger(__name__)
 
@@ -15,41 +15,145 @@ LOGGER = logging.getLogger(__name__)
 def temperature_feedback(temp: float) -> float:
     """Negative coefficient: higher temperature lowers reactivity."""
     reference = 900.0
-    coefficient = -5e-5
+    coefficient = -1.7e-5
     return coefficient * (temp - reference)
 
 
 def xenon_poisoning(flux: float) -> float:
-    return min(0.05, 1e-8 * flux)
+    return min(0.01, 5e-10 * flux)
 
 
 @dataclass
 class NeutronDynamics:
+    base_shutdown_bias: float = -0.014
+    shutdown_bias: float = -0.014
     beta_effective: float = 0.0065
-    delayed_neutron_fraction: float = 0.0008
+    delayed_decay_const: float = 0.08  # 1/s effective precursor decay
+    external_source_coupling: float = 1e-6
+    iodine_yield: float = 1e-6  # inventory units per MW*s
+    iodine_decay_const: float = 1.0 / 66000.0  # ~18h
+    xenon_decay_const: float = 1.0 / 33000.0  # ~9h
+    xenon_burnout_coeff: float = 1e-13  # per n/cm2
+    xenon_reactivity_coeff: float = 0.05
 
-    def reactivity(self, state: CoreState, control_fraction: float) -> float:
+    def reactivity(self, state: CoreState, control_fraction: float, rod_banks: list[float] | None = None) -> float:
+        if rod_banks:
+            weights = constants.CONTROL_ROD_BANK_WEIGHTS
+            worth = 0.0
+            total = sum(weights)
+            for w, pos in zip(weights, rod_banks):
+                worth += w * (1.0 - clamp(pos, 0.0, 0.95) / 0.95)
+            rod_term = constants.CONTROL_ROD_WORTH * worth / total
+        else:
+            rod_term = constants.CONTROL_ROD_WORTH * (1.0 - control_fraction)
         rho = (
-            0.02 * (1.0 - control_fraction)
+            self.shutdown_bias +
+            rod_term
             + temperature_feedback(state.fuel_temperature)
             - fuel_reactivity_penalty(state.burnup)
-            - xenon_poisoning(state.neutron_flux)
+            - self.xenon_penalty(state)
         )
         return rho
 
-    def flux_derivative(self, state: CoreState, rho: float) -> float:
+    def flux_derivative(
+        self,
+        state: CoreState,
+        rho: float,
+        delayed_source: float,
+        external_source_rate: float = 0.0,
+        baseline_source: float = 1e5,
+    ) -> float:
         generation_time = constants.NEUTRON_LIFETIME
         beta = self.beta_effective
-        return ((rho - beta) / generation_time) * state.neutron_flux + 1e5
+        source_term = self.external_source_coupling * external_source_rate
+        prompt = ((rho - beta) / generation_time) * state.neutron_flux
+        return prompt + delayed_source + baseline_source + source_term
 
-    def step(self, state: CoreState, control_fraction: float, dt: float) -> None:
-        rho = self.reactivity(state, control_fraction)
-        d_flux = self.flux_derivative(state, rho)
+    def step(
+        self,
+        state: CoreState,
+        control_fraction: float,
+        dt: float,
+        external_source_rate: float = 0.0,
+        rod_banks: list[float] | None = None,
+    ) -> None:
+        rho = self.reactivity(state, control_fraction, rod_banks)
+        rho = min(rho, 0.02)
+        shutdown = control_fraction >= 0.95
+        if shutdown:
+            rho = min(rho, -0.04)
+        baseline = 0.0 if shutdown else 1e5
+        source = 0.0 if shutdown else external_source_rate
+        rod_positions = rod_banks if rod_banks else [control_fraction] * len(constants.CONTROL_ROD_BANK_WEIGHTS)
+        self._ensure_precursors(state, len(rod_positions))
+        bank_factors = self._bank_factors(rod_positions)
+        bank_betas = self._bank_betas(len(bank_factors))
+        delayed_source = self._delayed_source(state, bank_factors)
+
+        d_flux = self.flux_derivative(
+            state,
+            rho,
+            delayed_source,
+            external_source_rate=source,
+            baseline_source=baseline,
+        )
         state.neutron_flux = max(0.0, state.neutron_flux + d_flux * dt)
         state.reactivity_margin = rho
+        self._update_precursors(state, bank_factors, bank_betas, dt)
         LOGGER.debug(
             "Neutronics: rho=%.5f, flux=%.2e n/cm2/s, d_flux=%.2e",
             rho,
             state.neutron_flux,
             d_flux,
         )
+
+    def update_poisons(self, state: CoreState, dt: float) -> None:
+        prod_I = max(0.0, state.power_output_mw) * self.iodine_yield
+        decay_I = state.iodine_inventory * self.iodine_decay_const
+        state.iodine_inventory = max(0.0, state.iodine_inventory + (prod_I - decay_I) * dt)
+        prod_Xe = decay_I
+        burn_Xe = state.neutron_flux * self.xenon_burnout_coeff
+        decay_Xe = state.xenon_inventory * self.xenon_decay_const
+        state.xenon_inventory = max(0.0, state.xenon_inventory + (prod_Xe - decay_Xe - burn_Xe) * dt)
+
+    def xenon_penalty(self, state: CoreState) -> float:
+        """Return delta-rho penalty from xenon inventory (positive magnitude)."""
+        return self._xenon_penalty(state)
+
+    def _xenon_penalty(self, state: CoreState) -> float:
+        return min(0.05, state.xenon_inventory * self.xenon_reactivity_coeff)
+
+    def _bank_betas(self, bank_count: int) -> list[float]:
+        weights = list(constants.CONTROL_ROD_BANK_WEIGHTS)
+        if bank_count != len(weights):
+            weights = [1.0 for _ in range(bank_count)]
+        total = sum(weights) if weights else 1.0
+        return [self.beta_effective * (w / total) for w in weights]
+
+    def _bank_factors(self, positions: list[float]) -> list[float]:
+        factors: list[float] = []
+        for pos in positions:
+            insertion = clamp(pos, 0.0, 0.95)
+            factors.append(max(0.0, 1.0 - insertion / 0.95))
+        return factors
+
+    def _ensure_precursors(self, state: CoreState, bank_count: int) -> None:
+        if not state.delayed_precursors or len(state.delayed_precursors) != bank_count:
+            state.delayed_precursors = [0.0 for _ in range(bank_count)]
+
+    def _delayed_source(self, state: CoreState, bank_factors: list[float]) -> float:
+        decay = self.delayed_decay_const
+        return sum(decay * precursor * factor for precursor, factor in zip(state.delayed_precursors, bank_factors))
+
+    def _update_precursors(
+        self, state: CoreState, bank_factors: list[float], bank_betas: list[float], dt: float
+    ) -> None:
+        generation_time = constants.NEUTRON_LIFETIME
+        decay = self.delayed_decay_const
+        new_pools: list[float] = []
+        for precursor, factor, beta in zip(state.delayed_precursors, bank_factors, bank_betas):
+            production = (beta / generation_time) * state.neutron_flux * factor
+            loss = decay * precursor
+            updated = max(0.0, precursor + (production - loss) * dt)
+            new_pools.append(updated)
+        state.delayed_precursors = new_pools

src/reactor_sim/rich_dashboard.py

@@ -0,0 +1,133 @@
+"""Alternate dashboard using rich Live rendering (non-interactive)."""
+
+from __future__ import annotations
+
+import logging
+from typing import Optional
+
+from rich.console import Group
+from rich.live import Live
+from rich.panel import Panel
+from rich.table import Table
+from rich.layout import Layout
+from rich.text import Text
+
+from .reactor import Reactor
+from .simulation import ReactorSimulation
+from .state import PlantState
+from . import constants
+
+LOGGER = logging.getLogger(__name__)
+
+
+def _table(title: str) -> Table:
+    tbl = Table.grid(expand=True, padding=(0, 1))
+    tbl.title = title
+    tbl.title_style = "bold cyan"
+    return tbl
+
+
+class RichDashboard:
+    """Read-only dashboard that refreshes plant metrics using rich."""
+
+    def __init__(
+        self,
+        reactor: Reactor,
+        start_state: Optional[PlantState],
+        timestep: float = 1.0,
+        save_path: Optional[str] = None,
+    ) -> None:
+        self.reactor = reactor
+        self.start_state = start_state
+        self.timestep = timestep
+        self.save_path = save_path
+        self.sim: Optional[ReactorSimulation] = None
+
+    def _render(self, state: PlantState) -> Layout:
+        layout = Layout()
+        layout.split_column(Layout(name="upper"), Layout(name="lower"))
+        layout["upper"].split_row(Layout(name="core"), Layout(name="loops"))
+        layout["lower"].split_row(Layout(name="turbine"), Layout(name="misc"))
+
+        core = _table("Core")
+        core.add_row(f"Power {state.core.power_output_mw:6.1f} MW", f"Fuel {state.core.fuel_temperature:6.1f} K")
+        core.add_row(f"Rods {self.reactor.control.rod_fraction:.3f}", f"Mode {'AUTO' if not self.reactor.control.manual_control else 'MAN'}")
+        core.add_row(f"Setpoint {self.reactor.control.setpoint_mw:5.0f} MW", f"Reactivity {state.core.reactivity_margin:+.4f}")
+        layout["upper"]["core"].update(Panel(core, padding=0))
+
+        loops = _table("Loops")
+        loops.add_row(
+            f"P pri {state.primary_loop.pressure:4.1f}/{constants.MAX_PRESSURE:4.1f} MPa",
+            f"Tin {state.primary_loop.temperature_in:5.1f}K",
+            f"Tout {state.primary_loop.temperature_out:5.1f}K",
+            f"Flow {state.primary_loop.mass_flow_rate:6.0f} kg/s",
+        )
+        loops.add_row(
+            f"P sec {state.secondary_loop.pressure:4.1f}/{constants.MAX_PRESSURE:4.1f} MPa",
+            f"Tin {state.secondary_loop.temperature_in:5.1f}K",
+            f"Tout {state.secondary_loop.temperature_out:5.1f}K",
+            f"q {state.secondary_loop.steam_quality:4.2f}",
+        )
+        loops.add_row(
+            f"HX ΔT {state.primary_to_secondary_delta_t:4.0f}K",
+            f"Eff {state.heat_exchanger_efficiency*100:5.1f}%",
+            f"Relief pri {'OPEN' if self.reactor.primary_relief_open else 'CLOSED'}",
+            f"Relief sec {'OPEN' if self.reactor.secondary_relief_open else 'CLOSED'}",
+        )
+        layout["upper"]["loops"].update(Panel(loops, padding=0))
+
+        turbine = _table("Turbines / Grid")
+        avail = getattr(self, "_steam_available_power", lambda s: 0.0)(state)  # type: ignore[arg-type]
+        steam_h = state.turbines[0].steam_enthalpy if state.turbines else 0.0
+        turbine.add_row(
+            f"Steam avail {avail:5.1f} MW",
+            f"h {steam_h:5.0f} kJ/kg",
+            f"Cond P {state.turbines[0].condenser_pressure:4.2f} MPa" if state.turbines else "Cond P n/a",
+        )
+        turbine.add_row(
+            f"Unit1 {state.turbines[0].electrical_output_mw:5.1f} MW" if state.turbines else "Unit1 n/a",
+            f"Unit2 {state.turbines[1].electrical_output_mw:5.1f} MW" if len(state.turbines) > 1 else "Unit2 n/a",
+            f"Unit3 {state.turbines[2].electrical_output_mw:5.1f} MW" if len(state.turbines) > 2 else "Unit3 n/a",
+        )
+        turbine.add_row(
+            f"Electrical {state.total_electrical_output():5.1f} MW",
+            f"Demand {self._total_load_demand(state):5.1f} MW",
+            f"Supplied {self._total_load_supplied(state):5.1f} MW",
+        )
+        layout["lower"]["turbine"].update(Panel(turbine, padding=0))
+
+        misc_group = []
+        misc_group.append(Text(f"Time {state.time_elapsed:6.1f}s", style="cyan"))
+        misc_group.append(Text(f"Primary pumps: {[p.status for p in state.primary_pumps] if state.primary_pumps else []}"))
+        misc_group.append(Text(f"Secondary pumps: {[p.status for p in state.secondary_pumps] if state.secondary_pumps else []}"))
+        misc_group.append(Text(f"Pressurizer level {self.reactor.pressurizer_level*100:5.1f}% @ {constants.PRIMARY_PRESSURIZER_SETPOINT_MPA:4.1f} MPa"))
+        misc_group.append(Text(f"Reliefs: pri={'OPEN' if self.reactor.primary_relief_open else 'CLOSED'} sec={'OPEN' if self.reactor.secondary_relief_open else 'CLOSED'}"))
+        if self.reactor.health_monitor.failure_log:
+            misc_group.append(Text(f"Failures: {', '.join(self.reactor.health_monitor.failure_log)}", style="bold red"))
+        layout["lower"]["misc"].update(Panel(Group(*misc_group), padding=0))
+        return layout
+
+    def _total_load_supplied(self, state: PlantState) -> float:
+        return sum(t.load_supplied_mw for t in state.turbines)
+
+    def _total_load_demand(self, state: PlantState) -> float:
+        return sum(t.load_demand_mw for t in state.turbines)
+
+    def run(self) -> None:
+        self.sim = ReactorSimulation(
+            self.reactor,
+            timestep=self.timestep,
+            duration=None,
+            realtime=True,
+        )
+        self.sim.start_state = self.start_state
+        try:
+            with Live(console=None, refresh_per_second=4) as live:
+                for state in self.sim.run():
+                    live.update(self._render(state))
+        except KeyboardInterrupt:
+            if self.sim:
+                self.sim.stop()
+        finally:
+            if self.save_path and self.sim and self.sim.last_state:
+                self.reactor.save_state(self.save_path, self.sim.last_state)

src/reactor_sim/simulation.py

@@ -7,6 +7,7 @@ import json
 import logging
 import os
 from dataclasses import dataclass, field
+from pathlib import Path
 from threading import Event
 import time
 from typing import Callable, Iterable, Optional
@@ -21,6 +22,34 @@ LOGGER = logging.getLogger(__name__)
 CommandProvider = Callable[[float, PlantState], Optional[ReactorCommand]]
 
 
+def _default_state_path() -> Path:
+    custom = os.getenv("FISSION_STATE_PATH")
+    return Path(custom) if custom else Path("artifacts/last_state.json")
+
+
+def _poison_log_path() -> Path:
+    custom = os.getenv("FISSION_POISON_LOG")
+    return Path(custom) if custom else Path("artifacts/poisons.json")
+
+
+def _write_poison_log(path: Path, state: PlantState | None, snapshots: list[dict[str, float]] | None = None) -> None:
+    if state is None and snapshots:
+        latest = snapshots[-1]
+        inventory = latest.get("products", {})
+        particles = latest.get("particles", {})
+        time_elapsed = latest.get("time_elapsed", 0.0)
+    elif state is not None:
+        inventory = state.core.fission_product_inventory
+        particles = state.core.emitted_particles
+        time_elapsed = state.time_elapsed
+    else:
+        return
+    path.parent.mkdir(parents=True, exist_ok=True)
+    payload = {"time_elapsed": time_elapsed, "products": inventory, "particles": particles}
+    path.write_text(json.dumps(payload, indent=2))
+    LOGGER.info("Wrote poison snapshot to %s", path)
+
+
 @dataclass
 class ReactorSimulation:
     reactor: Reactor
@@ -66,19 +95,68 @@ def main() -> None:
     log_file = os.getenv("FISSION_LOG_FILE")
     configure_logging(log_level, log_file)
     realtime = os.getenv("FISSION_REALTIME", "0") == "1"
+    alternate_dashboard = os.getenv("ALTERNATE_DASHBOARD", "0") == "1"
+    snapshot_at_env = os.getenv("FISSION_SNAPSHOT_AT")
+    snapshot_at = float(snapshot_at_env) if snapshot_at_env else None
+    snapshot_path = os.getenv("FISSION_SNAPSHOT_PATH", "artifacts/snapshot.txt")
     duration_env = os.getenv("FISSION_SIM_DURATION")
     if duration_env:
         duration = None if duration_env.lower() in {"none", "infinite"} else float(duration_env)
     else:
         duration = None if realtime else 600.0
     reactor = Reactor.default()
-    sim = ReactorSimulation(reactor, timestep=5.0, duration=duration, realtime=realtime)
+    dashboard_mode = os.getenv("FISSION_DASHBOARD", "0") == "1"
+    timestep_env = os.getenv("FISSION_TIMESTEP")
+    if timestep_env:
+        timestep = float(timestep_env)
+    else:
+        timestep = 1.0 if dashboard_mode or realtime else 5.0
+
+    sim = ReactorSimulation(reactor, timestep=timestep, duration=duration, realtime=realtime)
+    state_path = _default_state_path()
     load_path = os.getenv("FISSION_LOAD_STATE")
-    save_path = os.getenv("FISSION_SAVE_STATE")
+    if not load_path and state_path.exists():
+        load_path = str(state_path)
+    save_path = os.getenv("FISSION_SAVE_STATE") or str(state_path)
     if load_path:
         sim.start_state = reactor.load_state(load_path)
+    if dashboard_mode and snapshot_at is None:
+        if alternate_dashboard:
             try:
-        if realtime:
+                from .textual_dashboard import run_textual_dashboard
+            except ImportError as exc:  # pragma: no cover - optional dependency
+                LOGGER.error("Textual dashboard requested but dependency missing: %s", exc)
+                return
+            run_textual_dashboard(
+                reactor,
+                start_state=sim.start_state,
+                timestep=sim.timestep,
+                save_path=save_path,
+            )
+            return
+        else:
+            from .dashboard import ReactorDashboard
+
+            dashboard = ReactorDashboard(
+                reactor,
+                start_state=sim.start_state,
+                timestep=sim.timestep,
+                save_path=save_path,
+            )
+            dashboard.run()
+            return
+    try:
+        if snapshot_at is not None:
+            sim.duration = snapshot_at
+            LOGGER.info("Running headless to t=%.1fs for snapshot...", snapshot_at)
+            for _ in sim.run():
+                pass
+            if sim.last_state:
+                from .snapshot import write_snapshot
+
+                write_snapshot(snapshot_path, reactor, sim.last_state)
+                LOGGER.info("Snapshot written to %s", snapshot_path)
+        elif realtime:
             LOGGER.info("Running in real-time mode (Ctrl+C to stop)...")
             for _ in sim.run():
                 pass
@@ -86,6 +164,7 @@ def main() -> None:
             snapshots = sim.log()
             LOGGER.info("Captured %d snapshots", len(snapshots))
             print(json.dumps(snapshots[-5:], indent=2))
+            _write_poison_log(_poison_log_path(), sim.last_state, snapshots)
     except KeyboardInterrupt:
         sim.stop()
         LOGGER.warning("Simulation interrupted by user")

src/reactor_sim/snapshot.py

@@ -0,0 +1,53 @@
+"""Snapshot formatting helpers shared across dashboards."""
+
+from __future__ import annotations
+
+from pathlib import Path
+from typing import Iterable
+
+from . import constants
+from .reactor import Reactor
+from .state import PlantState
+
+
+def snapshot_lines(reactor: Reactor, state: PlantState) -> list[str]:
+    core = state.core
+    prim = state.primary_loop
+    sec = state.secondary_loop
+    t0 = state.turbines[0] if state.turbines else None
+    lines: list[str] = [
+        f"Time {state.time_elapsed:6.1f}s",
+        f"Core: {core.power_output_mw:6.1f}MW fuel {core.fuel_temperature:6.1f}K rods {reactor.control.rod_fraction:.3f} ({'AUTO' if not reactor.control.manual_control else 'MAN'})",
+        f"Primary: P={prim.pressure:4.1f}/{constants.MAX_PRESSURE:4.1f}MPa Tin={prim.temperature_in:6.1f}K Tout={prim.temperature_out:6.1f}K Flow={prim.mass_flow_rate:6.0f}kg/s",
+        f"Secondary: P={sec.pressure:4.1f}/{constants.MAX_PRESSURE:4.1f}MPa Tin={sec.temperature_in:6.1f}K Tout={sec.temperature_out:6.1f}K q={sec.steam_quality:4.2f} Flow={sec.mass_flow_rate:6.0f}kg/s",
+        f"HX ΔT={state.primary_to_secondary_delta_t:4.0f}K Eff={state.heat_exchanger_efficiency*100:5.1f}%",
+    ]
+    if t0:
+        lines.append(
+            f"Turbines: h={t0.steam_enthalpy:5.0f}kJ/kg avail={_steam_available_power(state):5.1f}MW "
+            f"CondP={t0.condenser_pressure:4.2f}/{constants.CONDENSER_MAX_PRESSURE_MPA:4.2f}MPa "
+            f"CondT={t0.condenser_temperature:6.1f}K"
+        )
+        lines.append("Outputs: " + " ".join([f"T{idx+1}:{t.electrical_output_mw:5.1f}MW" for idx, t in enumerate(state.turbines)]))
+    failures = ", ".join(reactor.health_monitor.failure_log) if reactor.health_monitor.failure_log else "None"
+    lines.append(
+        f"Status: pumps pri {[p.status for p in state.primary_pumps]} sec {[p.status for p in state.secondary_pumps]} "
+        f"relief pri={'OPEN' if reactor.primary_relief_open else 'CLOSED'} sec={'OPEN' if reactor.secondary_relief_open else 'CLOSED'} "
+        f"failures={failures}"
+    )
+    return lines
+
+
+def write_snapshot(path: Path | str, reactor: Reactor, state: PlantState) -> None:
+    p = Path(path)
+    p.parent.mkdir(parents=True, exist_ok=True)
+    p.write_text("\n".join(snapshot_lines(reactor, state)))
+
+
+def _steam_available_power(state: PlantState) -> float:
+    mass_flow = state.secondary_loop.mass_flow_rate * max(0.0, state.secondary_loop.steam_quality)
+    if mass_flow <= 1.0:
+        return 0.0
+    enthalpy = state.turbines[0].steam_enthalpy if state.turbines else (constants.STEAM_LATENT_HEAT / 1_000.0)
+    return (enthalpy * mass_flow) / 1_000.0
+

src/reactor_sim/state.py

@@ -4,6 +4,10 @@ from __future__ import annotations
 
 from dataclasses import dataclass, field, asdict
 
+from . import constants
+
+from .generator import GeneratorState
+
 
 def clamp(value: float, min_value: float, max_value: float) -> float:
     return max(min_value, min(max_value, value))
@@ -16,11 +20,39 @@ class CoreState:
     reactivity_margin: float  # delta rho
     power_output_mw: float  # MW thermal
     burnup: float  # fraction of fuel consumed
+    clad_temperature: float | None = None  # Kelvin
+    pellet_center_temperature: float | None = None  # Kelvin, peak centerline surrogate
+    gap_conductance: float = 1.0  # surrogate factor 0-1
+    dnb_margin: float | None = None  # ratio to critical heat flux surrogate
+    subcooling_margin: float | None = None  # K until boiling
+    xenon_inventory: float = 0.0
+    iodine_inventory: float = 0.0
+    delayed_precursors: list[float] = field(default_factory=list)
+    fission_product_inventory: dict[str, float] = field(default_factory=dict)
+    emitted_particles: dict[str, float] = field(default_factory=dict)
+
+    def __post_init__(self) -> None:
+        if self.clad_temperature is None:
+            self.clad_temperature = self.fuel_temperature
+        if self.pellet_center_temperature is None:
+            self.pellet_center_temperature = self.fuel_temperature
+        if self.dnb_margin is None:
+            self.dnb_margin = 1.0
+        if self.subcooling_margin is None:
+            self.subcooling_margin = 0.0
 
     def update_burnup(self, dt: float) -> None:
         produced_energy_mwh = self.power_output_mw * (dt / 3600.0)
         self.burnup = clamp(self.burnup + produced_energy_mwh * 1e-5, 0.0, 0.99)
 
+    def add_products(self, products: dict[str, float]) -> None:
+        for element, amount in products.items():
+            self.fission_product_inventory[element] = self.fission_product_inventory.get(element, 0.0) + amount
+
+    def add_emitted_particles(self, particles: dict[str, float]) -> None:
+        for name, amount in particles.items():
+            self.emitted_particles[name] = self.emitted_particles.get(name, 0.0) + amount
+
 
 @dataclass
 class CoolantLoopState:
@@ -29,6 +61,9 @@ class CoolantLoopState:
     pressure: float  # MPa
     mass_flow_rate: float  # kg/s
     steam_quality: float  # fraction of vapor
+    inventory_kg: float = 0.0  # bulk mass of coolant
+    level: float = 1.0  # fraction full relative to nominal volume
+    energy_j: float = 0.0  # stored thermal/latent energy for the loop
 
     def average_temperature(self) -> float:
         return 0.5 * (self.temperature_in + self.temperature_out)
@@ -40,8 +75,19 @@ class TurbineState:
     shaft_power_mw: float
     electrical_output_mw: float
     condenser_temperature: float
+    condenser_pressure: float = constants.CONDENSER_BASE_PRESSURE_MPA
+    fouling_penalty: float = 0.0
     load_demand_mw: float = 0.0
     load_supplied_mw: float = 0.0
+    status: str = "OFF"
+
+
+@dataclass
+class PumpState:
+    active: bool
+    flow_rate: float
+    pressure: float
+    status: str = "OFF"
 
 
 @dataclass
@@ -49,7 +95,17 @@ class PlantState:
     core: CoreState
     primary_loop: CoolantLoopState
     secondary_loop: CoolantLoopState
-    turbine: TurbineState
+    turbines: list[TurbineState]
+    primary_pumps: list[PumpState] = field(default_factory=list)
+    secondary_pumps: list[PumpState] = field(default_factory=list)
+    generators: list[GeneratorState] = field(default_factory=list)
+    aux_draws: dict[str, float] = field(default_factory=dict)
+    heat_exchanger_efficiency: float = 0.0
+    primary_to_secondary_delta_t: float = 0.0
+    dissolved_oxygen_ppm: float = constants.CHEM_OXYGEN_DEFAULT_PPM
+    boron_ppm: float = constants.CHEM_BORON_DEFAULT_PPM
+    sodium_ppm: float = constants.CHEM_SODIUM_DEFAULT_PPM
+    hx_fouling: float = 0.0
     time_elapsed: float = field(default=0.0)
 
     def snapshot(self) -> dict[str, float]:
@@ -60,18 +116,73 @@ class PlantState:
             "neutron_flux": self.core.neutron_flux,
             "primary_outlet_temp": self.primary_loop.temperature_out,
             "secondary_pressure": self.secondary_loop.pressure,
-            "turbine_electric": self.turbine.electrical_output_mw,
+            "turbine_electric": self.total_electrical_output(),
+            "products": self.core.fission_product_inventory,
+            "particles": self.core.emitted_particles,
+            "primary_pumps": [pump.active for pump in self.primary_pumps],
+            "secondary_pumps": [pump.active for pump in self.secondary_pumps],
+            "generators": [gen.running or gen.starting for gen in self.generators],
         }
 
+    def total_electrical_output(self) -> float:
+        return sum(t.electrical_output_mw for t in self.turbines)
+
     def to_dict(self) -> dict:
         return asdict(self)
 
     @classmethod
     def from_dict(cls, data: dict) -> "PlantState":
+        core_blob = dict(data["core"])
+        inventory = core_blob.pop("fission_product_inventory", {})
+        particles = core_blob.pop("emitted_particles", {})
+        # Backwards/forwards compatibility for optional core fields.
+        core_blob.pop("dnb_margin", None)
+        core_blob.pop("subcooling_margin", None)
+        core_blob.setdefault("clad_temperature", core_blob.get("fuel_temperature", 295.0))
+        core_blob.setdefault("pellet_center_temperature", core_blob.get("fuel_temperature", 295.0))
+        core_blob.setdefault("gap_conductance", 1.0)
+        turbines_blob = data.get("turbines")
+        if turbines_blob is None:
+            # Compatibility with previous single-turbine snapshots.
+            old_turbine = data.get("turbine")
+            turbines_blob = [old_turbine] if old_turbine else []
+        turbines = [TurbineState(**t) for t in turbines_blob]
+        prim_pumps_blob = data.get("primary_pumps", [])
+        sec_pumps_blob = data.get("secondary_pumps", [])
+        generators_blob = data.get("generators", [])
+        generators = [GeneratorState(**g) for g in generators_blob]
+        aux_draws = data.get("aux_draws", {})
+        hx_eff = data.get("heat_exchanger_efficiency", 0.0)
+        delta_t = data.get("primary_to_secondary_delta_t", 0.0)
+        dissolved_oxygen = data.get("dissolved_oxygen_ppm", constants.CHEM_OXYGEN_DEFAULT_PPM)
+        boron_ppm = data.get("boron_ppm", constants.CHEM_BORON_DEFAULT_PPM)
+        sodium_ppm = data.get("sodium_ppm", constants.CHEM_SODIUM_DEFAULT_PPM)
+        hx_fouling = data.get("hx_fouling", 0.0)
         return cls(
-            core=CoreState(**data["core"]),
-            primary_loop=CoolantLoopState(**data["primary_loop"]),
-            secondary_loop=CoolantLoopState(**data["secondary_loop"]),
-            turbine=TurbineState(**data["turbine"]),
+            core=CoreState(**core_blob, fission_product_inventory=inventory, emitted_particles=particles),
+            primary_loop=CoolantLoopState(**_with_energy(data["primary_loop"])),
+            secondary_loop=CoolantLoopState(**_with_energy(data["secondary_loop"])),
+            turbines=turbines,
+            primary_pumps=[PumpState(**p) for p in prim_pumps_blob],
+            secondary_pumps=[PumpState(**p) for p in sec_pumps_blob],
+            generators=generators,
+            aux_draws=aux_draws,
+            heat_exchanger_efficiency=hx_eff,
+            primary_to_secondary_delta_t=delta_t,
+            dissolved_oxygen_ppm=dissolved_oxygen,
+            boron_ppm=boron_ppm,
+            sodium_ppm=sodium_ppm,
+            hx_fouling=hx_fouling,
             time_elapsed=data.get("time_elapsed", 0.0),
         )
+
+
+def _with_energy(loop_blob: dict) -> dict:
+    """Backwards compatibility: derive energy if missing."""
+    if "energy_j" in loop_blob:
+        return loop_blob
+    energy = 0.5 * (loop_blob.get("temperature_in", 295.0) + loop_blob.get("temperature_out", 295.0))
+    energy *= loop_blob.get("inventory_kg", 0.0) * constants.COOLANT_HEAT_CAPACITY
+    out = dict(loop_blob)
+    out["energy_j"] = energy
+    return out

src/reactor_sim/textual_dashboard.py

@@ -0,0 +1,603 @@
+"""Interactive Textual-based dashboard mirroring the curses layout."""
+
+from __future__ import annotations
+
+import logging
+import os
+from collections import deque
+from pathlib import Path
+from typing import Optional
+
+from textual.app import App, ComposeResult
+from textual.containers import Grid, Horizontal, HorizontalScroll, Vertical, VerticalScroll
+from textual.widgets import Button, Footer, Header, Static
+from textual.timer import Timer
+
+from . import constants
+from .reactor import Reactor
+from .simulation import ReactorSimulation
+from .state import PlantState
+from .commands import ReactorCommand
+from .snapshot import snapshot_lines
+
+LOGGER = logging.getLogger(__name__)
+
+
+def _bar(label: str, value: float, width: int = 24) -> str:
+    filled = int(max(0.0, min(1.0, value)) * width)
+    return f"{label:<14} [{'#'*filled}{'-'*(width-filled)}] {value*100:5.1f}%"
+
+
+class TextualDashboard(App):
+    """Textual dashboard with controls and sections similar to the curses view."""
+
+    CSS = """
+    Screen { layout: vertical; }
+    #controls { padding: 0 0; height: auto; }
+    #controls .row { height: auto; }
+    #controls Button { min-width: 8; padding: 0 1; }
+    .columns { height: 1fr; }
+    .col { width: 1fr; height: 1fr; }
+    .panel { padding: 0 1; border: round $primary; }
+    """
+
+    BINDINGS = [
+        ("q", "quit", "Quit"),
+        ("space", "scram", "SCRAM"),
+        ("g", "toggle_primary_p1", "P1"),
+        ("h", "toggle_primary_p2", "P2"),
+        ("j", "toggle_secondary_p1", "S1"),
+        ("k", "toggle_secondary_p2", "S2"),
+        ("b", "toggle_generator1", "Gen1"),
+        ("v", "toggle_generator2", "Gen2"),
+        ("t", "toggle_turbine_bank", "Turbines"),
+        ("1", "toggle_turbine_unit('1')", "T1"),
+        ("2", "toggle_turbine_unit('2')", "T2"),
+        ("3", "toggle_turbine_unit('3')", "T3"),
+        ("+", "rods_insert", "+ rods"),
+        ("-", "rods_withdraw", "- rods"),
+        ("[", "demand_down", "Demand -"),
+        ("]", "demand_up", "Demand +"),
+        ("s", "setpoint_down", "SP -"),
+        ("d", "setpoint_up", "SP +"),
+        ("l", "toggle_primary_relief", "Relief pri"),
+        (";", "toggle_secondary_relief", "Relief sec"),
+        ("c", "toggle_consumer", "Consumer"),
+        ("a", "toggle_auto_rods", "Auto rods"),
+        ("f12", "snapshot", "Snapshot"),
+        # Maintenance is mouse-only here; keyboard remains in curses.
+    ]
+
+    timestep: float = 1.0
+
+    def __init__(
+        self,
+        reactor: Reactor,
+        start_state: Optional[PlantState],
+        timestep: float = 1.0,
+        save_path: Optional[str] = None,
+    ) -> None:
+        super().__init__()
+        self.reactor = reactor
+        self.state = start_state or self.reactor.initial_state()
+        self.timestep = timestep
+        self.save_path = save_path
+        self._pending: deque[ReactorCommand] = deque()
+        self._timer: Optional[Timer] = None
+        self._trend_history: deque[tuple[float, float, float]] = deque(maxlen=120)
+        self.log_buffer: deque[str] = deque(maxlen=8)
+        self._log_handler: Optional[logging.Handler] = None
+        self._previous_handlers: list[logging.Handler] = []
+        snap_at_env = os.getenv("FISSION_SNAPSHOT_AT")
+        self.snapshot_at = float(snap_at_env) if snap_at_env else None
+        self.snapshot_done = False
+        self.snapshot_path = Path(os.getenv("FISSION_SNAPSHOT_PATH", "artifacts/textual_snapshot.txt"))
+
+        # Panels
+        self.core_panel = Static(classes="panel")
+        self.trend_panel = Static(classes="panel")
+        self.poison_panel = Static(classes="panel")
+        self.primary_panel = Static(classes="panel")
+        self.secondary_panel = Static(classes="panel")
+        self.turbine_panel = Static(classes="panel")
+        self.generator_panel = Static(classes="panel")
+        self.power_panel = Static(classes="panel")
+        self.hx_panel = Static(classes="panel")
+        self.protection_panel = Static(classes="panel")
+        self.maintenance_panel = Static(classes="panel")
+        self.health_panel = Static(classes="panel")
+        self.help_panel = Static(classes="panel")
+        self.log_panel = Static(classes="panel")
+        self.status_panel = Static(classes="panel")
+        self.controls_panel: Grid | None = None
+
+    def compose(self) -> ComposeResult:
+        yield Header()
+        controls = self._build_controls()
+        yield controls
+        with Horizontal(classes="columns"):
+            with VerticalScroll(classes="col"):
+                yield self.core_panel
+                yield self.trend_panel
+                yield self.poison_panel
+                yield self.primary_panel
+                yield self.secondary_panel
+                yield self.power_panel
+                yield self.generator_panel
+                yield self.status_panel
+            with VerticalScroll(classes="col"):
+                yield self.turbine_panel
+                yield self.hx_panel
+                yield self.protection_panel
+                yield self.maintenance_panel
+                yield self.health_panel
+                yield self.log_panel
+                yield self.help_panel
+        yield Footer()
+
+    def on_mount(self) -> None:
+        self._install_log_capture()
+        self._timer = self.set_interval(self.timestep, self._tick, pause=False)
+        self._render_panels()
+
+    def action_quit(self) -> None:
+        if self._timer:
+            self._timer.pause()
+        if self.save_path and self.state:
+            self.reactor.save_state(self.save_path, self.state)
+        self._restore_logging()
+        self.exit()
+
+    # Command helpers
+    def _enqueue(self, cmd: ReactorCommand) -> None:
+        self._pending.append(cmd)
+
+    def action_scram(self) -> None:
+        self._enqueue(ReactorCommand.scram_all())
+
+    def action_toggle_primary_p1(self) -> None:
+        self._enqueue(ReactorCommand(primary_pumps={1: not self.reactor.primary_pump_units[0]}))
+
+    def action_toggle_primary_p2(self) -> None:
+        self._enqueue(ReactorCommand(primary_pumps={2: not self.reactor.primary_pump_units[1]}))
+
+    def action_toggle_secondary_p1(self) -> None:
+        self._enqueue(ReactorCommand(secondary_pumps={1: not self.reactor.secondary_pump_units[0]}))
+
+    def action_toggle_secondary_p2(self) -> None:
+        self._enqueue(ReactorCommand(secondary_pumps={2: not self.reactor.secondary_pump_units[1]}))
+
+    def action_toggle_generator1(self) -> None:
+        self._enqueue(ReactorCommand(generator_units={1: True}))
+
+    def action_toggle_generator2(self) -> None:
+        self._enqueue(ReactorCommand(generator_units={2: True}))
+
+    def action_toggle_turbine_bank(self) -> None:
+        self._enqueue(ReactorCommand(turbine_on=not self.reactor.turbine_active))
+
+    def action_toggle_turbine_unit(self, unit: str) -> None:
+        idx = int(unit)
+        current = self.reactor.turbine_unit_active[idx - 1] if idx - 1 < len(self.reactor.turbine_unit_active) else False
+        self._enqueue(ReactorCommand(turbine_units={idx: not current}))
+
+    def action_rods_insert(self) -> None:
+        self._enqueue(ReactorCommand(rod_position=min(0.95, self.reactor.control.rod_fraction + constants.ROD_MANUAL_STEP)))
+
+    def action_rods_withdraw(self) -> None:
+        self._enqueue(ReactorCommand(rod_position=max(0.0, self.reactor.control.rod_fraction - constants.ROD_MANUAL_STEP)))
+
+    def action_demand_down(self) -> None:
+        if self.reactor.consumer:
+            self._enqueue(ReactorCommand(consumer_demand=max(0.0, self.reactor.consumer.demand_mw - 50.0)))
+
+    def action_demand_up(self) -> None:
+        if self.reactor.consumer:
+            self._enqueue(ReactorCommand(consumer_demand=self.reactor.consumer.demand_mw + 50.0))
+
+    def action_setpoint_down(self) -> None:
+        self._enqueue(ReactorCommand(power_setpoint=self.reactor.control.setpoint_mw - 250.0))
+
+    def action_setpoint_up(self) -> None:
+        self._enqueue(ReactorCommand(power_setpoint=self.reactor.control.setpoint_mw + 250.0))
+
+    def action_toggle_primary_relief(self) -> None:
+        self._enqueue(ReactorCommand(primary_relief=not self.reactor.primary_relief_open))
+
+    def action_toggle_secondary_relief(self) -> None:
+        self._enqueue(ReactorCommand(secondary_relief=not self.reactor.secondary_relief_open))
+
+    def action_toggle_consumer(self) -> None:
+        if self.reactor.consumer:
+            self._enqueue(ReactorCommand(consumer_online=not self.reactor.consumer.online))
+
+    def action_toggle_auto_rods(self) -> None:
+        self._enqueue(ReactorCommand(rod_manual=not self.reactor.control.manual_control))
+
+    # Maintenance (mouse-driven)
+    def action_maintain_core(self) -> None:
+        self._enqueue(ReactorCommand.maintain("core"))
+
+    def action_maintain_primary_p1(self) -> None:
+        self._enqueue(ReactorCommand.maintain("primary_pump_1"))
+
+    def action_maintain_primary_p2(self) -> None:
+        self._enqueue(ReactorCommand.maintain("primary_pump_2"))
+
+    def action_maintain_secondary_p1(self) -> None:
+        self._enqueue(ReactorCommand.maintain("secondary_pump_1"))
+
+    def action_maintain_secondary_p2(self) -> None:
+        self._enqueue(ReactorCommand.maintain("secondary_pump_2"))
+
+    def action_maintain_turbine_1(self) -> None:
+        self._enqueue(ReactorCommand.maintain("turbine_1"))
+
+    def action_maintain_turbine_2(self) -> None:
+        self._enqueue(ReactorCommand.maintain("turbine_2"))
+
+    def action_maintain_turbine_3(self) -> None:
+        self._enqueue(ReactorCommand.maintain("turbine_3"))
+
+    def action_maintain_generator_1(self) -> None:
+        self._enqueue(ReactorCommand.maintain("generator_1"))
+
+    def action_maintain_generator_2(self) -> None:
+        self._enqueue(ReactorCommand.maintain("generator_2"))
+
+    def action_snapshot(self) -> None:
+        self._save_snapshot(auto=False)
+
+    def _merge_commands(self) -> Optional[ReactorCommand]:
+        if not self._pending:
+            return None
+        cmd = ReactorCommand()
+        while self._pending:
+            nxt = self._pending.popleft()
+            for field in nxt.__dataclass_fields__:  # type: ignore[attr-defined]
+                val = getattr(nxt, field)
+                if val is None or val is False:
+                    continue
+                setattr(cmd, field, val)
+        return cmd
+
+    def _tick(self) -> None:
+        cmd = self._merge_commands()
+        self.reactor.step(self.state, self.timestep, cmd)
+        self._render_panels()
+        if self.snapshot_at is not None and not self.snapshot_done and self.state.time_elapsed >= self.snapshot_at:
+            self._save_snapshot(auto=True)
+
+    def _build_controls(self) -> Vertical:
+        row1 = Horizontal(
+            Button("SCRAM", id="scram"),
+            Button("P1", id="p1"),
+            Button("P2", id="p2"),
+            Button("S1", id="s1"),
+            Button("S2", id="s2"),
+            Button("RelP", id="relief_pri"),
+            Button("RelS", id="relief_sec"),
+            Button("Turb", id="turbines"),
+            Button("T1", id="t1"),
+            Button("T2", id="t2"),
+            Button("T3", id="t3"),
+            Button("G1", id="g1"),
+            Button("G2", id="g2"),
+            Button("Grid", id="consumer"),
+            Button("AutoR", id="auto_rods"),
+            Button("+Rod", id="rods_plus"),
+            Button("-Rod", id="rods_minus"),
+            classes="row",
+        )
+        row2 = Horizontal(
+            Button("D+50", id="demand_up"),
+            Button("D-50", id="demand_down"),
+            Button("SP+250", id="sp_up"),
+            Button("SP-250", id="sp_down"),
+            Button("Snap", id="snapshot"),
+            classes="row",
+        )
+        row3 = Horizontal(
+            Button("MCore", id="m_core"),
+            Button("MP1", id="m_p1"),
+            Button("MP2", id="m_p2"),
+            Button("MS1", id="m_s1"),
+            Button("MS2", id="m_s2"),
+            Button("MT1", id="m_t1"),
+            Button("MT2", id="m_t2"),
+            Button("MT3", id="m_t3"),
+            Button("MG1", id="m_g1"),
+            Button("MG2", id="m_g2"),
+            classes="row",
+        )
+        return Vertical(
+            HorizontalScroll(row1, classes="row"),
+            HorizontalScroll(row2, classes="row"),
+            HorizontalScroll(row3, classes="row"),
+            id="controls",
+        )
+
+    def _render_panels(self) -> None:
+        self._trend_history.append((self.state.time_elapsed, self.state.core.fuel_temperature, self.state.core.power_output_mw))
+
+        self.core_panel.update(
+            "[bold cyan]Core[/bold cyan]\n"
+            f"Power {self.state.core.power_output_mw:6.1f} MW (Nom {constants.NORMAL_CORE_POWER_MW:4.0f}/Max {constants.TEST_MAX_POWER_MW:4.0f})\n"
+            f"Fuel {self.state.core.fuel_temperature:6.1f} K  Rods {self.reactor.control.rod_fraction:.3f} ({'AUTO' if not self.reactor.control.manual_control else 'MAN'})\n"
+            f"Setpoint {self.reactor.control.setpoint_mw:5.0f} MW  Reactivity {self.state.core.reactivity_margin:+.4f}\n"
+            f"DNB {self.state.core.dnb_margin:4.2f} Subcool {self.state.core.subcooling_margin:4.1f}K"
+        )
+        self.trend_panel.update(self._trend_text())
+        self.poison_panel.update(self._poison_text())
+        self.primary_panel.update(
+            "[bold cyan]Primary Loop[/bold cyan]\n"
+            f"Flow {self.state.primary_loop.mass_flow_rate:7.0f}/{self.reactor.primary_pump.nominal_flow * len(self.reactor.primary_pump_units):.0f} kg/s\n"
+            f"Level {self.state.primary_loop.level*100:6.1f}%\n"
+            f"Tin {self.state.primary_loop.temperature_in:7.1f} K  Tout {self.state.primary_loop.temperature_out:7.1f} K (Target {constants.PRIMARY_OUTLET_TARGET_K:4.0f})\n"
+            f"P {self.state.primary_loop.pressure:5.2f}/{constants.MAX_PRESSURE:4.1f} MPa  Pressurizer {self.reactor.pressurizer_level*100:6.1f}% @ {constants.PRIMARY_PRESSURIZER_SETPOINT_MPA:4.1f} MPa\n"
+            f"Relief {'OPEN' if self.reactor.primary_relief_open else 'CLOSED'} Pumps {[p.status for p in self.state.primary_pumps]}"
+        )
+        self.secondary_panel.update(
+            "[bold cyan]Secondary Loop[/bold cyan]\n"
+            f"Flow {self.state.secondary_loop.mass_flow_rate:7.0f}/{self.reactor.secondary_pump.nominal_flow * len(self.reactor.secondary_pump_units):.0f} kg/s\n"
+            f"Level {self.state.secondary_loop.level*100:6.1f}%\n"
+            f"Tin {self.state.secondary_loop.temperature_in:7.1f} K  Tout {self.state.secondary_loop.temperature_out:7.1f} K (Target {constants.SECONDARY_OUTLET_TARGET_K:4.0f})\n"
+            f"P {self.state.secondary_loop.pressure:5.2f}/{constants.MAX_PRESSURE:4.1f} MPa  q {self.state.secondary_loop.steam_quality:5.2f}/1.00\n"
+            f"Relief {'OPEN' if self.reactor.secondary_relief_open else 'CLOSED'} Pumps {[p.status for p in self.state.secondary_pumps]}"
+        )
+        self.turbine_panel.update(self._turbine_text())
+        self.generator_panel.update(self._generator_text())
+        self.power_panel.update(self._power_text())
+        self.hx_panel.update(
+            "[bold cyan]Heat Exchanger[/bold cyan]\n"
+            f"ΔT (pri-sec) {self.state.primary_to_secondary_delta_t:4.0f} K\n"
+            f"Efficiency {self.state.heat_exchanger_efficiency*100:5.1f}%"
+        )
+        self.protection_panel.update(self._protection_text())
+        self.maintenance_panel.update(self._maintenance_text())
+        self.health_panel.update(self._health_text())
+        self.log_panel.update(self._log_text())
+        self.help_panel.update(self._help_text())
+        failures = ", ".join(self.reactor.health_monitor.failure_log) if self.reactor.health_monitor.failure_log else "None"
+        self.status_panel.update(
+            "[bold cyan]Status[/bold cyan]\n"
+            f"Time {self.state.time_elapsed:6.1f}s\n"
+            f"Consumer {'ON' if (self.reactor.consumer and self.reactor.consumer.online) else 'OFF'} Demand {self.reactor.consumer.demand_mw if self.reactor.consumer else 0.0:5.1f} MW\n"
+            f"Failures: {failures}"
+        )
+
+    def _trend_text(self) -> str:
+        if len(self._trend_history) < 2:
+            return "[bold cyan]Trends[/bold cyan]\nFuel Temp Δ n/a\nCore Power Δ n/a"
+        start_t, start_temp, start_power = self._trend_history[0]
+        end_t, end_temp, end_power = self._trend_history[-1]
+        duration = max(1.0, end_t - start_t)
+        temp_delta = end_temp - start_temp
+        power_delta = end_power - start_power
+        temp_rate = temp_delta / duration
+        power_rate = power_delta / duration
+        return (
+            "[bold cyan]Trends[/bold cyan]\n"
+            f"Fuel Temp Δ {end_temp:7.1f} K (Δ{temp_delta:+6.1f} / {duration:4.0f}s, {temp_rate:+5.2f}/s)\n"
+            f"Core Power Δ {end_power:7.1f} MW (Δ{power_delta:+6.1f} / {duration:4.0f}s, {power_rate:+5.2f}/s)"
+        )
+
+    def _poison_text(self) -> str:
+        inventory = self.state.core.fission_product_inventory or {}
+        particles = self.state.core.emitted_particles or {}
+        xe = getattr(self.state.core, "xenon_inventory", 0.0)
+        sm = inventory.get("Sm", 0.0)
+        iodine = inventory.get("I", 0.0)
+        xe_drho = getattr(self.state.core, "reactivity_margin", 0.0)
+        return (
+            "[bold cyan]Key Poisons / Emitters[/bold cyan]\n"
+            f"Xe (xenon): {xe:9.2e}\n"
+            f"Sm (samarium): {sm:9.2e}\n"
+            f"I (iodine): {iodine:9.2e}\n"
+            f"Xe Δρ: {xe_drho:+.4f}\n"
+            f"Neutrons (src): {particles.get('neutrons', 0.0):.2e}"
+        )
+
+    def _turbine_text(self) -> str:
+        steam_avail = self._steam_available_power(self.state)
+        enthalpy = self.state.turbines[0].steam_enthalpy if self.state.turbines else 0.0
+        cond = ""
+        if self.state.turbines:
+            cond = (
+                f"Cond P {self.state.turbines[0].condenser_pressure:4.2f}/{constants.CONDENSER_MAX_PRESSURE_MPA:4.2f} MPa "
+                f"T {self.state.turbines[0].condenser_temperature:6.1f}K"
+            )
+        lines = [
+            "[bold cyan]Turbine / Grid[/bold cyan]",
+            f"Turbines {' '.join(self._turbine_status_lines()) if self._turbine_status_lines() else 'n/a'}",
+            f"Rated Elec {len(self.reactor.turbines)*self.reactor.turbines[0].rated_output_mw:7.1f} MW" if self.reactor.turbines else "Rated Elec n/a",
+            f"Steam Avail {steam_avail:5.1f} MW h={enthalpy:5.0f} kJ/kg {cond}",
+        ]
+        if self.state.turbines:
+            lines.append("Unit Elec " + " ".join([f"{t.electrical_output_mw:6.1f}MW" for t in self.state.turbines]))
+        lines.append(f"Throttle {self.reactor.turbines[0].throttle:5.2f}" if self.reactor.turbines else "Throttle n/a")
+        lines.append(f"Electrical {self.state.total_electrical_output():7.1f} MW Load {self._total_load_supplied(self.state):7.1f}/{self._total_load_demand(self.state):7.1f} MW")
+        if self.reactor.consumer:
+            lines.append(f"Consumer {'ONLINE' if self.reactor.consumer.online else 'OFF'} Demand {self.reactor.consumer.demand_mw:7.1f} MW")
+        return "\n".join(lines)
+
+    def _generator_text(self) -> str:
+        lines = ["[bold cyan]Generators[/bold cyan]"]
+        for idx, g in enumerate(self.state.generators):
+            status = "RUN" if g.running else "OFF"
+            if g.starting:
+                status = "START"
+            lines.append(f"Gen{idx+1}: {status:5} {g.power_output_mw:5.1f} MW batt {g.battery_charge*100:5.1f}%")
+        return "\n".join(lines)
+
+    def _power_text(self) -> str:
+        draws = getattr(self.state, "aux_draws", {}) or {}
+        base = draws.get("base", 0.0)
+        prim = draws.get("primary_pumps", 0.0)
+        sec = draws.get("secondary_pumps", 0.0)
+        demand = draws.get("total_demand", 0.0)
+        supplied = draws.get("supplied", 0.0)
+        gen_out = draws.get("generator_output", 0.0)
+        turb_out = draws.get("turbine_output", 0.0)
+        return (
+            "[bold cyan]Power Stats[/bold cyan]\n"
+            f"Base Aux {base:5.1f} MW  Prim Aux {prim:5.1f} MW  Sec Aux {sec:5.1f} MW\n"
+            f"Aux demand {demand:5.1f} MW supplied {supplied:5.1f} MW\n"
+            f"Gen out {gen_out:5.1f} MW Turbine out {turb_out:5.1f} MW"
+        )
+
+    def _protection_text(self) -> str:
+        lines = ["[bold cyan]Protections / Warnings[/bold cyan]"]
+        lines.append(f"SCRAM {'ACTIVE' if self.reactor.shutdown else 'CLEAR'}")
+        if self.reactor.meltdown:
+            lines.append("Meltdown IN PROGRESS")
+        sec_flow_low = self.state.secondary_loop.mass_flow_rate <= 1.0 or not self.reactor.secondary_pump_active
+        heat_sink_risk = sec_flow_low and self.state.core.power_output_mw > 50.0
+        if heat_sink_risk:
+            heat_text = "TRIPPED low secondary flow >50 MW"
+        elif sec_flow_low:
+            heat_text = "ARMED (secondary off/low flow)"
+        else:
+            heat_text = "OK"
+        lines.append(f"Heat sink {heat_text}")
+        lines.append(f"DNB margin {self.state.core.dnb_margin:4.2f}")
+        lines.append(f"Subcooling {self.state.core.subcooling_margin:5.1f} K")
+        lines.append(f"Reliefs pri={'OPEN' if self.reactor.primary_relief_open else 'CLOSED'} sec={'OPEN' if self.reactor.secondary_relief_open else 'CLOSED'}")
+        return "\n".join(lines)
+
+    def _maintenance_text(self) -> str:
+        active = list(self.reactor.maintenance_active)
+        return "[bold cyan]Maintenance[/bold cyan]\nActive: " + (", ".join(active) if active else "None")
+
+    def _health_text(self) -> str:
+        lines = ["[bold cyan]Component Health[/bold cyan]"]
+        for name, comp in self.reactor.health_monitor.components.items():
+            lines.append(_bar(name, comp.integrity))
+        return "\n".join(lines)
+
+    def _help_text(self) -> str:
+        tips = [
+            "Start pumps before withdrawing rods.",
+            "Bring turbine/consumer online after thermal stabilization.",
+            "Toggle turbine units 1/2/3 individually.",
+            "Use m/n/,/. in curses; mapped to j/k etc here.",
+            "F12 saves a snapshot; set FISSION_SNAPSHOT_AT for auto.",
+        ]
+        return "[bold cyan]Controls & Tips[/bold cyan]\n" + "\n".join(f"- {t}" for t in tips)
+
+    def _log_text(self) -> str:
+        lines = ["[bold cyan]Logs[/bold cyan]"]
+        if not self.log_buffer:
+            lines.append("No recent logs.")
+        else:
+            lines.extend(list(self.log_buffer))
+        return "\n".join(lines)
+
+    def _turbine_status_lines(self) -> list[str]:
+        if not self.reactor.turbine_unit_active:
+            return []
+        lines: list[str] = []
+        for idx, active in enumerate(self.reactor.turbine_unit_active):
+            label = f"{idx + 1}:"
+            status = "ON" if active else "OFF"
+            if idx < len(getattr(self.state, "turbines", [])):
+                t_state = self.state.turbines[idx]
+                status = getattr(t_state, "status", status)
+            lines.append(f"{label}{status}")
+        return lines
+
+    def _total_load_supplied(self, state: PlantState) -> float:
+        return sum(t.load_supplied_mw for t in state.turbines)
+
+    def _total_load_demand(self, state: PlantState) -> float:
+        return sum(t.load_demand_mw for t in state.turbines)
+
+    def _steam_available_power(self, state: PlantState) -> float:
+        mass_flow = state.secondary_loop.mass_flow_rate * max(0.0, state.secondary_loop.steam_quality)
+        if mass_flow <= 1.0:
+            return 0.0
+        enthalpy = state.turbines[0].steam_enthalpy if state.turbines else (constants.STEAM_LATENT_HEAT / 1_000.0)
+        return (enthalpy * mass_flow) / 1_000.0
+
+    def _snapshot_lines(self) -> list[str]:
+        return snapshot_lines(self.reactor, self.state)
+
+    def _save_snapshot(self, auto: bool = False) -> None:
+        try:
+            self.snapshot_path.parent.mkdir(parents=True, exist_ok=True)
+            self.snapshot_path.write_text("\n".join(self._snapshot_lines()))
+            self.snapshot_done = True
+            LOGGER.info("Saved dashboard snapshot to %s%s", self.snapshot_path, " (auto)" if auto else "")
+        except Exception as exc:  # pragma: no cover
+            LOGGER.error("Failed to save snapshot: %s", exc)
+
+    def _install_log_capture(self) -> None:
+        # Silence existing root handlers to avoid spewing logs over the UI.
+        root = logging.getLogger()
+        root.handlers = []
+        # Capture reactor_sim logs into the on-screen log buffer.
+        logger = logging.getLogger("reactor_sim")
+        logger.propagate = False
+        self._previous_handlers = list(logger.handlers)
+        handler = logging.StreamHandler()
+        handler.setLevel(logging.INFO)
+
+        def emit(record: logging.LogRecord) -> None:
+            msg = handler.format(record)
+            self.log_buffer.append(msg)
+
+        handler.emit = emit  # type: ignore[assignment]
+        logger.handlers = [handler]
+        logger.setLevel(logging.INFO)
+        self._log_handler = handler
+
+    def _restore_logging(self) -> None:
+        logger = logging.getLogger("reactor_sim")
+        if self._previous_handlers:
+            logger.handlers = self._previous_handlers
+        if self._log_handler and self._log_handler in logger.handlers:
+            logger.removeHandler(self._log_handler)
+
+    def on_button_pressed(self, event: Button.Pressed) -> None:  # type: ignore[override]
+        mapping = {
+            "scram": self.action_scram,
+            "p1": self.action_toggle_primary_p1,
+            "p2": self.action_toggle_primary_p2,
+            "s1": self.action_toggle_secondary_p1,
+            "s2": self.action_toggle_secondary_p2,
+            "relief_pri": self.action_toggle_primary_relief,
+            "relief_sec": self.action_toggle_secondary_relief,
+            "turbines": self.action_toggle_turbine_bank,
+            "t1": lambda: self.action_toggle_turbine_unit("1"),
+            "t2": lambda: self.action_toggle_turbine_unit("2"),
+            "t3": lambda: self.action_toggle_turbine_unit("3"),
+            "g1": self.action_toggle_generator1,
+            "g2": self.action_toggle_generator2,
+            "consumer": self.action_toggle_consumer,
+            "auto_rods": self.action_toggle_auto_rods,
+            "rods_plus": self.action_rods_insert,
+            "rods_minus": self.action_rods_withdraw,
+            "demand_up": self.action_demand_up,
+            "demand_down": self.action_demand_down,
+            "sp_up": self.action_setpoint_up,
+            "sp_down": self.action_setpoint_down,
+            "snapshot": lambda: self._save_snapshot(auto=False),
+            "m_core": self.action_maintain_core,
+            "m_p1": self.action_maintain_primary_p1,
+            "m_p2": self.action_maintain_primary_p2,
+            "m_s1": self.action_maintain_secondary_p1,
+            "m_s2": self.action_maintain_secondary_p2,
+            "m_t1": self.action_maintain_turbine_1,
+            "m_t2": self.action_maintain_turbine_2,
+            "m_t3": self.action_maintain_turbine_3,
+            "m_g1": self.action_maintain_generator_1,
+            "m_g2": self.action_maintain_generator_2,
+        }
+        handler = mapping.get(event.button.id or "")
+        if handler:
+            handler()
+
+
+def run_textual_dashboard(reactor: Reactor, start_state: Optional[PlantState], timestep: float, save_path: Optional[str]) -> None:
+    app = TextualDashboard(reactor, start_state, timestep=timestep, save_path=save_path)
+    app.run()

src/reactor_sim/thermal.py

@@ -5,18 +5,45 @@ from __future__ import annotations
 from dataclasses import dataclass
 import logging
 
+import math
+
 from . import constants
 from .state import CoolantLoopState, CoreState
 
 LOGGER = logging.getLogger(__name__)
 
 
-def heat_transfer(primary: CoolantLoopState, secondary: CoolantLoopState, core_power_mw: float) -> float:
+def heat_transfer(
+    primary: CoolantLoopState, secondary: CoolantLoopState, core_power_mw: float, fouling_factor: float = 0.0
+) -> float:
     """Return MW transferred to the secondary loop."""
-    delta_t = max(0.0, primary.temperature_out - secondary.temperature_in)
-    conductance = 0.05  # steam generator effectiveness
-    transferred = min(core_power_mw, conductance * delta_t)
-    LOGGER.debug("Heat transfer %.2f MW with ΔT=%.1fK", transferred, delta_t)
+    if primary.mass_flow_rate <= 0.0 or secondary.mass_flow_rate <= 0.0:
+        return 0.0
+    delta_t1 = max(1e-3, primary.temperature_out - secondary.temperature_in)
+    delta_t2 = max(1e-3, primary.temperature_in - secondary.temperature_out)
+    if delta_t1 <= 0.0 or delta_t2 <= 0.0:
+        return 0.0
+    if abs(delta_t1 - delta_t2) < 1e-6:
+        lmtd = delta_t1
+    else:
+        lmtd = (delta_t1 - delta_t2) / math.log(delta_t1 / delta_t2)
+    fouling = max(0.0, min(constants.HX_FOULING_MAX_PENALTY, fouling_factor))
+    ua = constants.STEAM_GENERATOR_UA_MW_PER_K * (1.0 - fouling)
+    ua_limited = ua * lmtd
+
+    # Prevent the heat exchanger from over-transferring and inverting the outlet temperatures.
+    primary_capacity = primary.mass_flow_rate * constants.COOLANT_HEAT_CAPACITY
+    secondary_capacity = secondary.mass_flow_rate * constants.COOLANT_HEAT_CAPACITY
+    approach = 2.0  # K minimum approach between loop outlets
+    pinch_limited = 0.0
+    if primary_capacity > 0.0 and secondary_capacity > 0.0:
+        temp_gap = primary.temperature_out - secondary.temperature_in - approach
+        if temp_gap > 0.0:
+            pinch_watts = temp_gap / ((1.0 / primary_capacity) + (1.0 / secondary_capacity))
+            pinch_limited = max(0.0, pinch_watts / constants.MEGAWATT)
+
+    transferred = max(0.0, min(core_power_mw, ua_limited, pinch_limited))
+    LOGGER.debug("Heat transfer %.2f MW with LMTD=%.1fK (ΔT1=%.1f ΔT2=%.1f)", transferred, lmtd, delta_t1, delta_t2)
     return transferred
 
 
@@ -26,15 +53,108 @@ def temperature_rise(power_mw: float, mass_flow: float) -> float:
     return (power_mw * constants.MEGAWATT) / (mass_flow * constants.COOLANT_HEAT_CAPACITY)
 
 
+def saturation_pressure(temp_k: float) -> float:
+    """Return approximate saturation pressure (MPa) for water at temp_k."""
+    temp_c = max(0.0, temp_k - 273.15)
+    # Antoine equation parameters for water in 1-374C range, pressure in mmHg.
+    a, b, c = 8.14019, 1810.94, 244.485
+    psat_mmHg = 10 ** (a - (b / (c + temp_c)))
+    psat_mpa = psat_mmHg * 133.322e-6
+    return min(constants.MAX_PRESSURE, max(0.01, psat_mpa))
+
+
+def saturation_temperature(pressure_mpa: float) -> float:
+    """Approximate saturation temperature (K) for water at the given pressure."""
+    target = max(0.01, min(constants.MAX_PRESSURE, pressure_mpa))
+    low, high = 273.15, 900.0
+    for _ in range(40):
+        mid = 0.5 * (low + high)
+        if saturation_pressure(mid) < target:
+            low = mid
+        else:
+            high = mid
+    return high
+
+
 @dataclass
 class ThermalSolver:
     primary_volume_m3: float = 300.0
 
-    def step_core(self, core: CoreState, primary: CoolantLoopState, power_mw: float, dt: float) -> None:
+    def _resolve_secondary_state(self, secondary: CoolantLoopState) -> None:
+        """Project stored energy onto temperature, quality, and pressure."""
+        cp = constants.COOLANT_HEAT_CAPACITY
+        mass = max(1e-6, secondary.inventory_kg)
+        secondary.energy_j = max(0.0, secondary.energy_j)
+        sat_temp = saturation_temperature(max(0.05, secondary.pressure))
+        liquid_energy = mass * cp * sat_temp
+        available = secondary.energy_j
+
+        if available <= liquid_energy:
+            temp = available / (mass * cp)
+            secondary.temperature_out = max(constants.ENVIRONMENT_TEMPERATURE, temp)
+            secondary.steam_quality = 0.0
+        else:
+            latent_energy = min(available - liquid_energy, mass * constants.STEAM_LATENT_HEAT)
+            quality = latent_energy / (mass * constants.STEAM_LATENT_HEAT)
+            superheat_energy = max(0.0, available - liquid_energy - latent_energy)
+            superheat_temp = superheat_energy / (mass * cp) if quality >= 1.0 else 0.0
+            secondary.temperature_out = sat_temp + superheat_temp
+            secondary.steam_quality = max(0.0, min(1.0, quality))
+            secondary.energy_j = liquid_energy + latent_energy + superheat_energy
+
+        secondary.pressure = min(
+            constants.MAX_PRESSURE, max(secondary.pressure, saturation_pressure(secondary.temperature_out))
+        )
+
+    def step_core(
+        self,
+        core: CoreState,
+        primary: CoolantLoopState,
+        power_mw: float,
+        dt: float,
+        residual_power_mw: float | None = None,
+    ) -> None:
+        def _lag(prev: float, new: float, tau: float) -> float:
+            if tau <= 0.0:
+                return new
+            alpha = min(1.0, max(0.0, dt / max(1e-6, tau)))
+            return prev + alpha * (new - prev)
+
+        if residual_power_mw is None:
+            residual_power_mw = power_mw
+        prev_fuel = core.fuel_temperature
+        prev_clad = core.clad_temperature or primary.temperature_out
         temp_rise = temperature_rise(power_mw, primary.mass_flow_rate)
         primary.temperature_out = primary.temperature_in + temp_rise
-        core.fuel_temperature += 0.1 * (power_mw - temp_rise) * dt
+        # Fuel heats from total fission power (even when most is convected) plus any residual left in the coolant.
+        heating = (0.003 * max(0.0, power_mw) + 0.01 * max(0.0, residual_power_mw)) * dt
+        cooling = 0.025 * max(0.0, core.fuel_temperature - primary.temperature_out) * dt
+        core.fuel_temperature += heating - cooling
+        # Simple radial split: fuel -> clad conduction with burnup-driven conductivity drop and gap penalty.
+        gap_penalty = 1.0 + 2.0 * core.burnup
+        conductivity = 0.03 / gap_penalty
+        conduction = conductivity * max(0.0, core.fuel_temperature - (core.clad_temperature or primary.temperature_out)) * dt
+        clad = core.clad_temperature or primary.temperature_out
+        clad_cooling = 0.06 * max(0.0, clad - primary.temperature_out) * dt
+        clad = max(primary.temperature_out, clad + conduction - clad_cooling)
+        core.fuel_temperature = max(primary.temperature_out, core.fuel_temperature - conduction)
+        # Peak pellet centerline surrogate based on power and burnup conductivity loss.
+        peak_delta = 20.0 * (core.power_output_mw / max(1.0, constants.NORMAL_CORE_POWER_MW)) * (1.0 + core.burnup)
+        core.pellet_center_temperature = min(constants.MAX_CORE_TEMPERATURE, core.fuel_temperature + peak_delta)
+        # Keep temperatures bounded and never below coolant outlet.
         core.fuel_temperature = min(core.fuel_temperature, constants.MAX_CORE_TEMPERATURE)
+        core.clad_temperature = min(clad, constants.MAX_CORE_TEMPERATURE)
+        # Apply mild lags so heat moves from fuel to clad to coolant over a short time constant.
+        core.fuel_temperature = _lag(prev_fuel, core.fuel_temperature, constants.FUEL_TO_CLAD_TIME_CONSTANT)
+        core.clad_temperature = _lag(prev_clad, core.clad_temperature or prev_clad, constants.CLAD_TO_COOLANT_TIME_CONSTANT)
+        core.subcooling_margin = max(0.0, saturation_temperature(primary.pressure) - primary.temperature_out)
+        chf = self._critical_heat_flux(primary)
+        heat_flux = (power_mw * constants.MEGAWATT) / max(1.0, self._core_surface_area())
+        core.dnb_margin = max(0.0, chf / max(1e-6, heat_flux))
+        avg_temp = 0.5 * (primary.temperature_in + primary.temperature_out)
+        primary.energy_j = max(
+            0.0, primary.inventory_kg * constants.COOLANT_HEAT_CAPACITY * avg_temp
+        )
         LOGGER.debug(
             "Primary loop: flow=%.0f kg/s temp_out=%.1fK core_temp=%.1fK",
             primary.mass_flow_rate,
@@ -42,14 +162,50 @@ class ThermalSolver:
             core.fuel_temperature,
         )
 
-    def step_secondary(self, secondary: CoolantLoopState, transferred_mw: float) -> None:
-        delta_t = temperature_rise(transferred_mw, secondary.mass_flow_rate)
-        secondary.temperature_out = secondary.temperature_in + delta_t
-        secondary.steam_quality = min(1.0, max(0.0, delta_t / 100.0))
-        secondary.pressure = min(constants.MAX_PRESSURE, 6.0 + delta_t * 0.01)
+    def step_secondary(self, secondary: CoolantLoopState, transferred_mw: float, dt: float = 1.0) -> None:
+        """Update secondary loop using a stored-energy steam-drum balance."""
+        cp = constants.COOLANT_HEAT_CAPACITY
+        mass = max(1e-6, secondary.inventory_kg)
+        if secondary.energy_j <= 0.0:
+            secondary.energy_j = mass * cp * secondary.average_temperature()
+
+        # Add transferred heat; if no heat, bleed toward ambient.
+        if transferred_mw > 0.0:
+            secondary.energy_j += transferred_mw * constants.MEGAWATT * dt
+        else:
+            bleed = 0.01 * (secondary.temperature_out - constants.ENVIRONMENT_TEMPERATURE)
+            secondary.energy_j = max(
+                mass * cp * constants.ENVIRONMENT_TEMPERATURE, secondary.energy_j - max(0.0, bleed) * mass * cp * dt
+            )
+
+        self._resolve_secondary_state(secondary)
         LOGGER.debug(
-            "Secondary loop: transferred=%.1fMW temp_out=%.1fK quality=%.2f",
+            "Secondary loop: transferred=%.1fMW temp_out=%.1fK quality=%.2f energy=%.1eJ",
             transferred_mw,
             secondary.temperature_out,
             secondary.steam_quality,
+            secondary.energy_j,
         )
+
+    def remove_steam_energy(self, secondary: CoolantLoopState, steam_power_mw: float, dt: float) -> None:
+        """Remove steam enthalpy consumed by turbines and rebalance the drum."""
+        if steam_power_mw <= 0.0:
+            return
+        secondary.energy_j = max(0.0, secondary.energy_j - steam_power_mw * constants.MEGAWATT * dt)
+        self._resolve_secondary_state(secondary)
+
+    def _critical_heat_flux(self, primary: CoolantLoopState) -> float:
+        """CHF surrogate with pressure, mass flux, and subcooling influence (Chen-like)."""
+        area = self._core_surface_area()
+        mass_flux = max(1.0, primary.mass_flow_rate / max(1.0, area))  # kg/s per m2 surrogate
+        pressure_factor = 0.6 + 0.4 * min(1.0, primary.pressure / max(0.1, constants.MAX_PRESSURE))
+        subcool = max(0.0, saturation_temperature(primary.pressure) - primary.temperature_out)
+        subcool_factor = max(0.2, min(1.0, subcool / 30.0))
+        flux_factor = max(0.4, min(3.0, (mass_flux / 500.0) ** 0.5))
+        base_chf = 1.2e7  # W/m2 surrogate baseline
+        return base_chf * flux_factor * pressure_factor * subcool_factor
+
+    def _core_surface_area(self) -> float:
+        # Simple surrogate: area scaling with volume^(2/3)
+        volume = self.primary_volume_m3
+        return max(1.0, (volume ** (2.0 / 3.0)) * 10.0)

src/reactor_sim/turbine.py

@@ -6,10 +6,8 @@ from dataclasses import dataclass
 import logging
 
 from . import constants
-from typing import Optional
-
+from .thermal import saturation_temperature, saturation_pressure
 from .state import CoolantLoopState, TurbineState
-from .consumer import ElectricalConsumer
 
 LOGGER = logging.getLogger(__name__)
 
@@ -28,42 +26,83 @@ class SteamGenerator:
 class Turbine:
     generator_efficiency: float = constants.GENERATOR_EFFICIENCY
     mechanical_efficiency: float = constants.STEAM_TURBINE_EFFICIENCY
+    rated_output_mw: float = 400.0  # cap per unit electrical output
+    spool_time: float = constants.TURBINE_SPOOL_TIME
+    throttle: float = 1.0  # 0-1 valve position
 
     def step(
         self,
         loop: CoolantLoopState,
         state: TurbineState,
-        consumer: Optional[ElectricalConsumer] = None,
+        steam_power_mw: float = 0.0,
+        dt: float = 1.0,
     ) -> None:
-        enthalpy = 2_700.0 + loop.steam_quality * 600.0
-        mass_flow = loop.mass_flow_rate * 0.6
-        shaft_power_mw = (enthalpy * mass_flow / 1_000.0) * self.mechanical_efficiency / 1_000.0
+        effective_mass_flow = loop.mass_flow_rate * max(0.0, loop.steam_quality)
+        if steam_power_mw <= 0.0 and (loop.steam_quality <= 0.01 or effective_mass_flow <= 10.0):
+            # No steam available; turbine should idle.
+            state.shaft_power_mw = 0.0
+            state.electrical_output_mw = 0.0
+            state.load_demand_mw = 0.0
+            state.load_supplied_mw = 0.0
+            state.steam_enthalpy = 0.0
+            state.condenser_temperature = max(constants.CONDENSER_COOLING_WATER_TEMP_K, loop.temperature_in - 20.0)
+            state.condenser_pressure = max(constants.CONDENSER_BASE_PRESSURE_MPA, state.condenser_pressure - 0.01 * dt)
+            state.fouling_penalty = max(0.0, state.fouling_penalty - 0.0001 * dt)
+            return
+
+        throttle = min(constants.TURBINE_THROTTLE_MAX, max(constants.TURBINE_THROTTLE_MIN, self.throttle))
+        throttle_eff = 1.0 - constants.TURBINE_THROTTLE_EFFICIENCY_DROP * (constants.TURBINE_THROTTLE_MAX - throttle)
+
+        sat_temp = saturation_temperature(max(0.05, loop.pressure))
+        superheat = max(0.0, loop.temperature_out - sat_temp)
+        enthalpy = (constants.STEAM_LATENT_HEAT / 1_000.0) + (superheat * constants.COOLANT_HEAT_CAPACITY / 1_000.0)
+        mass_flow = effective_mass_flow * 0.6 * throttle
+        computed_power = (enthalpy * mass_flow) / 1_000.0  # MW from enthalpy flow
+        available_power = steam_power_mw if steam_power_mw > 0 else computed_power
+        available_power = min(available_power, computed_power)
+        backpressure_loss = 1.0 - _backpressure_penalty(state)
+        shaft_power_mw = available_power * self.mechanical_efficiency * throttle_eff * backpressure_loss
         electrical = shaft_power_mw * self.generator_efficiency
-        if consumer:
-            load_demand = consumer.request_power()
-            supplied = min(electrical, load_demand)
-            consumer.update_power_received(supplied)
-            LOGGER.debug(
-                "Consumer %s demand %.1f -> supplied %.1f MW",
-                consumer.name,
-                load_demand,
-                supplied,
+        if electrical > self.rated_output_mw:
+            electrical = self.rated_output_mw
+            shaft_power_mw = electrical / max(1e-6, self.generator_efficiency)
+        condenser_temp = max(constants.CONDENSER_COOLING_WATER_TEMP_K, loop.temperature_in - 20.0)
+        # Vacuum pump tends toward base pressure; fouling raises it slowly when hot.
+        target_pressure = constants.CONDENSER_BASE_PRESSURE_MPA
+        if condenser_temp > constants.CONDENSER_COOLING_WATER_TEMP_K + 20.0:
+            state.fouling_penalty = min(
+                constants.CONDENSER_FOULING_MAX_PENALTY,
+                state.fouling_penalty + constants.CONDENSER_FOULING_RATE * dt,
+            )
+        state.condenser_pressure = max(
+            target_pressure,
+            min(constants.CONDENSER_MAX_PRESSURE_MPA, state.condenser_pressure - constants.CONDENSER_VACUUM_PUMP_RATE * dt),
         )
-        else:
-            load_demand = 0.0
-            supplied = 0.0
-        condenser_temp = max(305.0, loop.temperature_in - 20.0)
         state.steam_enthalpy = enthalpy
-        state.shaft_power_mw = shaft_power_mw
-        state.electrical_output_mw = electrical
+        state.shaft_power_mw = _ramp(state.shaft_power_mw, shaft_power_mw, dt, self.spool_time)
+        state.electrical_output_mw = _ramp(state.electrical_output_mw, electrical, dt, self.spool_time)
         state.condenser_temperature = condenser_temp
-        state.load_demand_mw = load_demand
-        state.load_supplied_mw = supplied
         LOGGER.debug(
-            "Turbine output: shaft=%.1fMW electrical=%.1fMW condenser=%.1fK load %.1f/%.1f",
+            "Turbine output: shaft=%.1fMW electrical=%.1fMW condenser=%.1fK",
             shaft_power_mw,
             electrical,
             condenser_temp,
-            supplied,
-            load_demand,
         )
+
+
+def _ramp(current: float, target: float, dt: float, time_constant: float) -> float:
+    if time_constant <= 0.0:
+        return target
+    alpha = min(1.0, max(0.0, dt / max(1e-6, time_constant)))
+    return current + (target - current) * alpha
+
+
+def _backpressure_penalty(state: TurbineState) -> float:
+    base = constants.CONDENSER_BASE_PRESSURE_MPA
+    max_p = constants.CONDENSER_MAX_PRESSURE_MPA
+    pressure = max(base, min(max_p, state.condenser_pressure))
+    if pressure <= base:
+        return min(constants.CONDENSER_BACKPRESSURE_PENALTY, state.fouling_penalty)
+    frac = (pressure - base) / max(1e-6, max_p - base)
+    penalty = frac * constants.CONDENSER_BACKPRESSURE_PENALTY
+    return min(constants.CONDENSER_BACKPRESSURE_PENALTY + state.fouling_penalty, penalty + state.fouling_penalty)

tests/test_atomic.py

@@ -0,0 +1,54 @@
+from reactor_sim.atomic import AtomicPhysics, make_atom
+
+
+def test_alpha_decay_reduces_mass_and_charge():
+    physics = AtomicPhysics(seed=1)
+    parent = make_atom(92, 143)
+    event = physics.alpha_decay(parent)
+    daughter, alpha = event.products
+    assert daughter.protons == 90
+    assert daughter.neutrons == 141
+    assert alpha.protons == 2 and alpha.neutrons == 2
+    assert "alpha" in event.emitted_particles
+    assert event.energy_mev > 0
+
+
+def test_beta_minus_conserves_mass_number():
+    physics = AtomicPhysics(seed=2)
+    parent = make_atom(26, 30)
+    event = physics.beta_minus_decay(parent)
+    (daughter,) = event.products
+    assert daughter.mass_number == parent.mass_number
+    assert daughter.protons == parent.protons + 1
+    assert "beta-" in event.emitted_particles
+
+
+def test_beta_plus_and_electron_capture_lower_proton_count():
+    physics = AtomicPhysics(seed=3)
+    parent = make_atom(15, 16)
+    beta_plus = physics.beta_plus_decay(parent)
+    capture = physics.electron_capture(parent)
+    assert beta_plus.products[0].protons == parent.protons - 1
+    assert capture.products[0].protons == parent.protons - 1
+    assert "beta+" in beta_plus.emitted_particles
+    assert "gamma" in capture.emitted_particles
+
+
+def test_gamma_emission_keeps_nuclide_same():
+    physics = AtomicPhysics(seed=4)
+    parent = make_atom(8, 8)
+    event = physics.gamma_emission(parent, energy_mev=0.3)
+    (product,) = event.products
+    assert product == parent
+    assert "gamma" in event.emitted_particles
+    assert event.energy_mev == 0.3
+
+
+def test_spontaneous_fission_creates_two_fragments():
+    physics = AtomicPhysics(seed=5)
+    parent = make_atom(92, 143)
+    event = physics.spontaneous_fission(parent)
+    assert len(event.products) == 2
+    total_mass = sum(atom.mass_number for atom in event.products)
+    assert total_mass <= parent.mass_number
+    assert event.energy_mev > 0

tests/test_aux_power.py

@@ -0,0 +1,18 @@
+from reactor_sim.reactor import Reactor
+from reactor_sim.commands import ReactorCommand
+
+
+def test_pump_stays_off_without_aux_power():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.shutdown = False
+    reactor.primary_pump_active = True
+    reactor.primary_pump_units = [True, False]
+    reactor.secondary_pump_active = False
+    reactor.generator_auto = False
+    reactor.turbine_unit_active = [False, False, False]
+
+    reactor.step(state, dt=1.0, command=ReactorCommand(primary_pumps={1: True}))
+
+    assert state.primary_loop.mass_flow_rate == 0.0
+    assert state.primary_pumps[0].status in {"OFF", "STOPPING"}

tests/test_control.py

@@ -0,0 +1,40 @@
+import pytest
+
+from reactor_sim.control import ControlSystem
+from reactor_sim import constants
+from reactor_sim.state import CoreState
+
+
+def _core_state() -> CoreState:
+    return CoreState(
+        fuel_temperature=300.0,
+        neutron_flux=1e5,
+        reactivity_margin=0.0,
+        power_output_mw=0.0,
+        burnup=0.0,
+    )
+
+
+def test_manual_rods_quantized_to_step():
+    control = ControlSystem()
+    control.manual_control = True
+    core = _core_state()
+
+    control.set_rods(0.333)
+    assert control.rod_target == 0.325
+    control.update_rods(core, dt=100.0)
+    assert control.rod_fraction == pytest.approx(0.325, rel=1e-6)
+
+    control.increment_rods(0.014)
+    assert control.rod_target == pytest.approx(0.35)
+    control.update_rods(core, dt=100.0)
+    assert control.rod_fraction == pytest.approx(0.35, rel=1e-6)
+
+    # Clamp upper bound
+    control.set_rods(1.0)
+    control.update_rods(core, dt=100.0)
+    assert control.rod_fraction == 0.95
+
+
+def test_dashboard_step_constant_exposed():
+    assert constants.ROD_MANUAL_STEP == 0.025

tests/test_fuel.py

@@ -0,0 +1,40 @@
+from reactor_sim.atomic import AtomicPhysics
+from reactor_sim.fuel import FuelAssembly
+from reactor_sim.state import CoreState
+
+
+def _core_state(
+    temp: float = 600.0,
+    flux: float = 1e6,
+    burnup: float = 0.05,
+    power: float = 100.0,
+) -> CoreState:
+    return CoreState(
+        fuel_temperature=temp,
+        neutron_flux=flux,
+        reactivity_margin=0.0,
+        power_output_mw=power,
+        burnup=burnup,
+    )
+
+
+def test_decay_reaction_power_positive_and_capped():
+    physics = AtomicPhysics(seed=7)
+    fuel = FuelAssembly(enrichment=0.045, mass_kg=80_000.0, atomic_physics=physics)
+    state = _core_state()
+    power = fuel.decay_reaction_power(state)
+    assert power > 0
+    # Should stay under 10% of current power_output_mw.
+    assert power <= state.power_output_mw * 0.1 + 1e-6
+
+
+def test_decay_reaction_power_scales_with_burnup_and_flux():
+    physics = AtomicPhysics(seed=9)
+    fuel = FuelAssembly(enrichment=0.045, mass_kg=80_000.0, atomic_physics=physics)
+    low = _core_state(burnup=0.01, flux=1e5, power=50.0)
+    high = _core_state(burnup=0.5, flux=5e6, power=200.0)
+    low_power = fuel.decay_reaction_power(low)
+    high_power = fuel.decay_reaction_power(high)
+    assert high_power > low_power
+    # Still capped to a reasonable fraction of the prevailing power.
+    assert high_power <= high.power_output_mw * 0.1 + 1e-6

tests/test_neutronics.py

@@ -0,0 +1,57 @@
+import pytest
+
+from reactor_sim.neutronics import NeutronDynamics
+from reactor_sim.state import CoreState
+
+
+def _core_state(
+    temperature: float = 950.0,
+    flux: float = 2e7,
+    burnup: float = 0.01,
+    power: float = 500.0,
+) -> CoreState:
+    return CoreState(
+        fuel_temperature=temperature,
+        neutron_flux=flux,
+        reactivity_margin=0.0,
+        power_output_mw=power,
+        burnup=burnup,
+    )
+
+
+def test_reactivity_increases_with_rod_withdrawal():
+    dynamics = NeutronDynamics()
+    state = _core_state()
+    rho_full_out = dynamics.reactivity(state, control_fraction=0.0)
+    rho_half = dynamics.reactivity(state, control_fraction=0.5)
+    assert rho_full_out > rho_half
+
+
+def test_poisons_accumulate_under_power():
+    dynamics = NeutronDynamics()
+    state = _core_state(power=800.0, flux=1e6)
+    dynamics.update_poisons(state, dt=100.0)
+    dynamics.update_poisons(state, dt=100.0)
+    assert state.iodine_inventory > 0.0
+    assert state.xenon_inventory > 0.0
+
+
+def test_xenon_penalty_caps():
+    dynamics = NeutronDynamics()
+    state = _core_state()
+    state.xenon_inventory = 50.0
+    assert dynamics.xenon_penalty(state) == 0.05
+    state.xenon_inventory = 0.2
+    assert dynamics.xenon_penalty(state) == pytest.approx(0.01)
+
+
+def test_delayed_precursors_follow_rod_banks():
+    dynamics = NeutronDynamics()
+    state = _core_state(power=600.0, flux=5e6)
+    dynamics.step(state, control_fraction=0.2, dt=1.0, rod_banks=[0.2, 0.2, 0.2])
+    initial_sum = sum(state.delayed_precursors)
+    assert len(state.delayed_precursors) == 3
+    assert initial_sum > 0.0
+
+    dynamics.step(state, control_fraction=0.95, dt=2.0, rod_banks=[0.95, 0.95, 0.95])
+    assert sum(state.delayed_precursors) < initial_sum

tests/test_pressurizer.py

@@ -0,0 +1,36 @@
+from reactor_sim.reactor import Reactor
+from reactor_sim.commands import ReactorCommand
+
+
+def test_pressurizer_raises_pressure_with_level():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    state.primary_loop.pressure = 0.5
+    reactor.pressurizer_level = 0.8
+    reactor.primary_pump_active = False
+    reactor.secondary_pump_active = False
+    reactor.shutdown = False
+
+    reactor.step(state, dt=1.0, command=ReactorCommand())
+
+    assert state.primary_loop.pressure > 0.5
+    assert reactor.pressurizer_level < 0.8
+
+
+def test_low_npsh_limits_primary_flow():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.shutdown = False
+    reactor.control.manual_control = True
+    reactor.control.rod_fraction = 0.0
+    reactor.primary_pump_active = True
+    reactor.primary_pump_units = [True, True]
+    reactor.secondary_pump_active = False
+    state.primary_loop.pressure = 0.05  # near-vacuum to force cavitation
+    state.primary_loop.temperature_in = 400.0
+    state.primary_loop.temperature_out = 600.0
+
+    reactor.step(state, dt=1.0, command=ReactorCommand(generator_units={1: True}))
+
+    assert state.primary_pumps[0].status == "CAV"
+    assert state.primary_loop.mass_flow_rate < 100.0

tests/test_simulation.py

@@ -4,6 +4,8 @@ from pathlib import Path
 import pytest
 
 from reactor_sim import constants
+from reactor_sim.commands import ReactorCommand
+from reactor_sim.failures import HealthMonitor
 from reactor_sim.reactor import Reactor
 from reactor_sim.simulation import ReactorSimulation
 
@@ -13,11 +15,12 @@ def test_reactor_initial_state_is_cold():
     state = reactor.initial_state()
     assert state.core.fuel_temperature == constants.ENVIRONMENT_TEMPERATURE
     assert state.primary_loop.mass_flow_rate == 0.0
-    assert state.turbine.electrical_output_mw == 0.0
+    assert state.total_electrical_output() == 0.0
 
 
 def test_state_save_and_load_roundtrip(tmp_path: Path):
     reactor = Reactor.default()
+    reactor.control.manual_control = True
     sim = ReactorSimulation(reactor, timestep=5.0, duration=15.0)
     sim.log()
     save_path = tmp_path / "plant_state.json"
@@ -30,13 +33,346 @@ def test_state_save_and_load_roundtrip(tmp_path: Path):
         sim.last_state.core.fuel_temperature
     )
     assert restored_reactor.control.rod_fraction == reactor.control.rod_fraction
+    assert restored_reactor.control.manual_control == reactor.control.manual_control
+    assert len(restored_state.primary_pumps) == 2
+    assert len(restored_state.secondary_pumps) == 2
 
 
 def test_health_monitor_flags_core_failure():
     reactor = Reactor.default()
     state = reactor.initial_state()
     state.core.fuel_temperature = constants.MAX_CORE_TEMPERATURE
-    failures = reactor.health_monitor.evaluate(state, True, True, True, dt=200.0)
+    failures = reactor.health_monitor.evaluate(
+        state, [True, True], [True, True], [True, True, True], state.generators, dt=200.0
+    )
     assert "core" in failures
     reactor._handle_failure("core")
     assert reactor.shutdown is True
+
+
+def test_maintenance_recovers_component_health():
+    monitor = HealthMonitor()
+    pump = monitor.component("secondary_pump_1")
+    pump.integrity = 0.3
+    pump.fail()
+    restored = monitor.maintain("secondary_pump_1", amount=0.5)
+    assert restored is True
+    assert pump.integrity == pytest.approx(0.8)
+    assert pump.failed is False
+
+
+def test_secondary_pump_loss_triggers_scram_and_no_steam():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    # Make sure some power is present to trigger heat-sink logic.
+    state.core.power_output_mw = 500.0
+    reactor.secondary_pump_active = False
+    reactor.control.manual_control = True
+    reactor.control.rod_fraction = 0.1
+    reactor.step(state, dt=1.0)
+    assert reactor.shutdown is True
+    assert all(t.electrical_output_mw == 0.0 for t in state.turbines)
+
+
+def test_cold_shutdown_stays_subcritical():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.control.manual_control = True
+    reactor.control.rod_fraction = 0.95
+    reactor.primary_pump_active = False
+    reactor.secondary_pump_active = False
+    initial_power = state.core.power_output_mw
+    for _ in range(10):
+        reactor.step(state, dt=1.0)
+    assert state.core.power_output_mw <= initial_power + 0.5
+    assert reactor.shutdown is True
+
+
+def test_toggle_maintenance_progresses_until_restored():
+    reactor = Reactor.default()
+    reactor.primary_pump_units = [False, False]
+    reactor.primary_pump_active = False
+    pump = reactor.health_monitor.component("primary_pump_1")
+    pump.integrity = 0.2
+
+    def provider(t: float, _state):
+        if t == 0:
+            return ReactorCommand.maintain("primary_pump_1")
+        return None
+
+    sim = ReactorSimulation(reactor, timestep=1.0, duration=50.0, command_provider=provider)
+    sim.log()
+    assert pump.integrity >= 0.99
+    assert "primary_pump_1" not in reactor.maintenance_active
+
+
+def test_primary_pump_unit_toggle_updates_active_flag():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.primary_pump_active = True
+    reactor.primary_pump_units = [True, True]
+
+    reactor.step(state, dt=1.0, command=ReactorCommand(primary_pumps={1: False}))
+    assert reactor.primary_pump_units == [False, True]
+    assert reactor.primary_pump_active is True
+
+    reactor.step(state, dt=1.0, command=ReactorCommand(primary_pumps={2: False}))
+    assert reactor.primary_pump_units == [False, False]
+    assert reactor.primary_pump_active is False
+
+
+def test_secondary_pump_unit_toggle_can_restart_pump():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.secondary_pump_active = False
+    reactor.secondary_pump_units = [False, False]
+
+    cmd = ReactorCommand(secondary_pumps={1: True}, generator_units={1: True})
+    for _ in range(5):
+        reactor.step(state, dt=1.0, command=cmd)
+        cmd = ReactorCommand(coolant_demand=0.75)
+
+    assert reactor.secondary_pump_units == [True, False]
+    assert reactor.secondary_pump_active is True
+    assert state.secondary_loop.mass_flow_rate > 0.0
+
+
+def test_primary_pumps_spool_up_over_seconds():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.secondary_pump_units = [False, False]
+    # Enable both pumps and command full flow; spool should take multiple steps.
+    target_flow = reactor.primary_pump.flow_rate(1.0) * len(reactor.primary_pump_units)
+    cmd = ReactorCommand(primary_pumps={1: True, 2: True}, generator_units={1: True}, coolant_demand=1.0)
+    reactor.step(state, dt=1.0, command=cmd)
+    reactor.step(state, dt=1.0, command=ReactorCommand(coolant_demand=1.0))
+    first_flow = state.primary_loop.mass_flow_rate
+    assert 0.0 < first_flow < target_flow
+
+    for _ in range(10):
+        reactor.step(state, dt=1.0, command=ReactorCommand(coolant_demand=1.0))
+
+    assert state.primary_loop.mass_flow_rate == pytest.approx(target_flow, rel=0.15)
+
+
+def test_full_rod_withdrawal_reaches_gigawatt_power():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.shutdown = False
+    reactor.control.manual_control = True
+    reactor.control.rod_fraction = 0.0
+    reactor.primary_pump_active = True
+    reactor.secondary_pump_active = True
+    reactor.primary_pump_units = [True, True]
+    reactor.secondary_pump_units = [True, True]
+    reactor.generator_auto = True
+    reactor.step(state, dt=1.0, command=ReactorCommand(generator_units={1: True}))
+
+    early_power = 0.0
+    for step in range(60):
+        reactor.step(state, dt=1.0)
+        if step == 10:
+            early_power = state.core.power_output_mw
+    assert state.core.power_output_mw > max(2_000.0, early_power * 2)
+    assert state.core.fuel_temperature > 600.0
+
+
+def test_partially_inserted_rods_hold_near_three_gw():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.shutdown = False
+    reactor.control.set_manual_mode(False)
+    reactor.primary_pump_active = True
+    reactor.secondary_pump_active = True
+    reactor.primary_pump_units = [True, True]
+    reactor.secondary_pump_units = [True, True]
+    reactor.generator_auto = True
+    reactor.step(state, dt=1.0, command=ReactorCommand(generator_units={1: True}))
+
+    for _ in range(180):
+        reactor.step(state, dt=1.0)
+
+    assert 2_500.0 < state.core.power_output_mw < 3_500.0
+    assert 0.05 < reactor.control.rod_fraction < 0.9
+
+
+def test_generator_spools_and_powers_pumps():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.shutdown = False
+    reactor.control.manual_control = True
+    reactor.control.rod_fraction = 0.95  # keep power low; focus on aux power
+    reactor.turbine_unit_active = [False, False, False]
+    reactor.secondary_pump_units = [False, False]
+
+    for step in range(12):
+        cmd = ReactorCommand(generator_units={1: True}, primary_pumps={1: True}) if step == 0 else None
+        reactor.step(state, dt=1.0, command=cmd)
+
+    assert state.generators and state.generators[0].running is True
+    assert state.generators[0].power_output_mw > 0.0
+    assert state.primary_loop.mass_flow_rate > 0.0
+
+
+def test_generator_manual_mode_allows_single_unit_and_stop():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.shutdown = False
+    reactor.control.manual_control = True
+    reactor.control.rod_fraction = 0.95
+    reactor.generator_auto = False
+    reactor.primary_pump_units = [True, False]
+    reactor.secondary_pump_units = [False, False]
+
+    reactor.step(state, dt=1.0, command=ReactorCommand(generator_units={1: True}, primary_pumps={1: True}))
+    assert state.generators[0].starting or state.generators[0].running
+
+    for _ in range(15):
+        reactor.step(state, dt=1.0)
+
+    assert state.generators[0].running is True
+    reactor.step(state, dt=1.0, command=ReactorCommand(generator_units={1: False}))
+    for _ in range(5):
+        reactor.step(state, dt=1.0)
+    assert state.generators[0].running is False
+    assert state.generators[1].running is False
+
+
+def test_meltdown_triggers_shutdown():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.shutdown = False
+    reactor.control.manual_control = True
+    reactor.control.rod_fraction = 0.0
+    reactor.primary_pump_active = True
+    reactor.secondary_pump_active = True
+    state.core.fuel_temperature = constants.CORE_MELTDOWN_TEMPERATURE + 50.0
+
+    reactor.step(state, dt=1.0)
+
+    assert reactor.shutdown is True
+    assert reactor.meltdown is True
+
+
+def test_auto_control_resets_shutdown_and_moves_rods():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.shutdown = True
+    reactor.control.manual_control = True
+    reactor.control.rod_fraction = 0.95
+
+    reactor.step(state, dt=1.0, command=ReactorCommand(rod_manual=False))
+
+    assert reactor.shutdown is False
+    assert reactor.control.manual_control is False
+    assert reactor.control.rod_fraction < 0.95
+
+
+def test_chemistry_builds_fouling_and_backpressure():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    # Push impurities high to accelerate fouling dynamics.
+    state.dissolved_oxygen_ppm = 200.0
+    state.sodium_ppm = 100.0
+    state.secondary_loop.mass_flow_rate = 20_000.0
+    state.secondary_loop.steam_quality = 0.3
+    state.secondary_loop.temperature_out = 600.0
+    state.secondary_loop.temperature_in = 560.0
+    base_hx = state.hx_fouling
+    base_foul = state.turbines[0].fouling_penalty if state.turbines else 0.0
+    base_pressure = state.turbines[0].condenser_pressure if state.turbines else constants.CONDENSER_BASE_PRESSURE_MPA
+
+    reactor._update_chemistry(state, dt=20.0)
+
+    assert state.hx_fouling > base_hx
+    if state.turbines:
+        assert state.turbines[0].fouling_penalty > base_foul
+        assert state.turbines[0].condenser_pressure >= base_pressure
+
+
+def test_full_power_reaches_steam_and_turbine_output():
+    """Integration: ramp to full power with staged rod control and verify sustained steam/electric output."""
+    reactor = Reactor.default()
+    reactor.health_monitor.disable_degradation = True
+    reactor.allow_external_aux = True
+    reactor.relaxed_npsh = True
+    reactor.control.set_power_setpoint(3_000.0)
+    state = reactor.initial_state()
+    reactor.step(
+        state,
+        dt=1.0,
+        command=ReactorCommand(
+            generator_units={1: True, 2: True},
+            primary_pumps={1: True, 2: True},
+            secondary_pumps={1: True, 2: True},
+            rod_manual=True,
+            rod_position=0.55,
+        ),
+    )
+    checkpoints = {300, 600, 900, 1800, 2700, 3600}
+    results = {}
+    turbines_started = False
+    for i in range(3600):
+        cmd = None
+        if state.core.power_output_mw >= 2_500.0 and reactor.control.manual_control:
+            cmd = ReactorCommand(rod_manual=False)
+        if (
+            not turbines_started
+            and state.secondary_loop.steam_quality > 0.02
+            and state.secondary_loop.pressure > 1.0
+        ):
+            cmd = ReactorCommand(turbine_on=True, turbine_units={1: True, 2: True, 3: True})
+            turbines_started = True
+        if i == 600 and not turbines_started:
+            cmd = ReactorCommand(turbine_on=True, turbine_units={1: True, 2: True, 3: True})
+            turbines_started = True
+        reactor.step(state, dt=1.0, command=cmd)
+        if state.time_elapsed in checkpoints:
+            results[state.time_elapsed] = {
+                "quality": state.secondary_loop.steam_quality,
+                "electric": state.total_electrical_output(),
+                "core_temp": state.core.fuel_temperature,
+            }
+
+    # At or after 10 minutes of operation, ensure we have meaningful steam and electrical output.
+    assert results[600]["quality"] > 0.05
+    assert results[600]["electric"] > 50.0
+    assert results[3600]["quality"] > 0.1
+    assert results[3600]["electric"] > 150.0
+    # No runaway core temperatures.
+    assert results[3600]["core_temp"] < constants.CORE_MELTDOWN_TEMPERATURE
+
+
+def test_cooldown_reaches_ambient_without_runaway():
+    """Shutdown with pumps running should cool loops toward ambient, no runaway."""
+    reactor = Reactor.default()
+    reactor.health_monitor.disable_degradation = True
+    reactor.allow_external_aux = True
+    reactor.relaxed_npsh = True
+    state = reactor.initial_state()
+    # Start hot.
+    reactor.step(
+        state,
+        dt=1.0,
+        command=ReactorCommand(
+            generator_units={1: True, 2: True},
+            primary_pumps={1: True, 2: True},
+            secondary_pumps={1: True, 2: True},
+            rod_manual=True,
+            rod_position=0.55,
+        ),
+    )
+    turbines_started = False
+    for i in range(1800):
+        cmd = None
+        if not turbines_started and state.secondary_loop.steam_quality > 0.02 and state.secondary_loop.pressure > 1.0:
+            cmd = ReactorCommand(turbine_on=True, turbine_units={1: True, 2: True, 3: True})
+            turbines_started = True
+        if i == 900:
+            cmd = ReactorCommand(rod_position=0.95, turbine_on=False, turbine_units={1: False, 2: False, 3: False})
+        reactor.step(state, dt=1.0, command=cmd)
+
+    assert not reactor.meltdown
+    assert state.core.power_output_mw < 1.0
+    assert state.primary_loop.temperature_out < 320.0
+    assert state.secondary_loop.temperature_out < 320.0

tests/test_thermal.py

@@ -0,0 +1,50 @@
+import pytest
+
+from reactor_sim import constants
+from reactor_sim.state import CoolantLoopState
+from reactor_sim.thermal import ThermalSolver, heat_transfer, saturation_temperature
+
+
+def _secondary_loop(temp_in: float = 350.0, pressure: float = 0.5, flow: float = 200.0) -> CoolantLoopState:
+    nominal_mass = constants.SECONDARY_LOOP_VOLUME_M3 * constants.COOLANT_DENSITY
+    return CoolantLoopState(
+        temperature_in=temp_in,
+        temperature_out=temp_in,
+        pressure=pressure,
+        mass_flow_rate=flow,
+        steam_quality=0.0,
+        inventory_kg=nominal_mass,
+        level=constants.SECONDARY_INVENTORY_TARGET,
+    )
+
+
+def test_secondary_heats_to_saturation_before_boiling():
+    solver = ThermalSolver()
+    loop = _secondary_loop(temp_in=330.0, flow=180.0, pressure=0.5)
+    sat_temp = saturation_temperature(loop.pressure)
+    solver.step_secondary(loop, transferred_mw=50.0, dt=1.0)
+    assert loop.steam_quality == pytest.approx(0.0)
+    assert loop.temperature_out < sat_temp
+
+
+def test_secondary_generates_steam_when_energy_exceeds_sensible_heat():
+    solver = ThermalSolver()
+    loop = _secondary_loop(temp_in=330.0, flow=180.0, pressure=0.5)
+    loop.inventory_kg *= 0.1  # reduce mass to let boil-up happen quickly
+    sat_temp = saturation_temperature(loop.pressure)
+    solver.step_secondary(loop, transferred_mw=120.0, dt=100.0)
+    assert loop.temperature_out == pytest.approx(sat_temp, rel=0.05)
+    assert loop.steam_quality > 0.0
+    assert loop.steam_quality < 1.0
+
+
+def test_heat_transfer_reduced_by_fouling():
+    primary = CoolantLoopState(
+        temperature_in=360.0, temperature_out=380.0, pressure=15.0, mass_flow_rate=50_000.0, steam_quality=0.0
+    )
+    secondary = CoolantLoopState(
+        temperature_in=320.0, temperature_out=330.0, pressure=6.5, mass_flow_rate=50_000.0, steam_quality=0.1
+    )
+    clean = heat_transfer(primary, secondary, core_power_mw=7_000.0, fouling_factor=0.0)
+    fouled = heat_transfer(primary, secondary, core_power_mw=7_000.0, fouling_factor=0.25)
+    assert fouled < clean

tests/test_turbine.py

@@ -0,0 +1,57 @@
+import pytest
+
+from reactor_sim.state import CoolantLoopState, TurbineState
+from reactor_sim.turbine import Turbine
+
+
+def test_turbine_spools_toward_target_output():
+    turbine = Turbine()
+    turbine.throttle = 1.0
+    loop = CoolantLoopState(
+        temperature_in=600.0,
+        temperature_out=650.0,
+        pressure=0.02,
+        mass_flow_rate=20_000.0,
+        steam_quality=0.9,
+    )
+    state = TurbineState(
+        steam_enthalpy=0.0,
+        shaft_power_mw=0.0,
+        electrical_output_mw=0.0,
+        condenser_temperature=300.0,
+    )
+    target_electric = min(
+        turbine.rated_output_mw,
+        300.0 * turbine.mechanical_efficiency * turbine.generator_efficiency,
+    )
+
+    dt = 1.0
+    turbine.step(loop, state, steam_power_mw=300.0, dt=dt)
+    assert 0.0 < state.electrical_output_mw < target_electric
+
+    for _ in range(60):
+        turbine.step(loop, state, steam_power_mw=300.0, dt=dt)
+
+    assert state.electrical_output_mw == pytest.approx(target_electric, rel=0.05)
+
+
+def test_turbine_produces_no_power_without_steam():
+    turbine = Turbine()
+    loop = CoolantLoopState(
+        temperature_in=295.0,
+        temperature_out=295.0,
+        pressure=6.0,
+        mass_flow_rate=20_000.0,
+        steam_quality=0.0,
+    )
+    state = TurbineState(
+        steam_enthalpy=0.0,
+        shaft_power_mw=0.0,
+        electrical_output_mw=0.0,
+        condenser_temperature=300.0,
+    )
+
+    turbine.step(loop, state, steam_power_mw=0.0, dt=1.0)
+
+    assert state.electrical_output_mw == 0.0
+    assert state.shaft_power_mw == 0.0