Improve reactor controls, reactions, and maintenance UI

AGENTS.md

@@ -23,6 +23,7 @@ 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.

artifacts/last_state.json

@@ -0,0 +1,104 @@
+{
+  "control": {
+    "setpoint_mw": 3000.0,
+    "rod_fraction": 0.3872181667930225
+  },
+  "plant": {
+    "core": {
+      "fuel_temperature": 633.2342060602825,
+      "neutron_flux": 31335131.43776027,
+      "reactivity_margin": 0.006241667151538859,
+      "power_output_mw": 2976.900110634399,
+      "burnup": 0.004945521266135161
+    },
+    "primary_loop": {
+      "temperature_in": 295.0,
+      "temperature_out": 338.7522062115579,
+      "pressure": 14.0,
+      "mass_flow_rate": 16200.0,
+      "steam_quality": 0.0
+    },
+    "secondary_loop": {
+      "temperature_in": 295.0,
+      "temperature_out": 360.05377003828596,
+      "pressure": 6.650537700382859,
+      "mass_flow_rate": 10880.0,
+      "steam_quality": 0.6505377003828593
+    },
+    "turbines": [
+      {
+        "steam_enthalpy": 3090.3226202297155,
+        "shaft_power_mw": 336.9056685758782,
+        "electrical_output_mw": 323.4294418328431,
+        "condenser_temperature": 305.0,
+        "load_demand_mw": 0.0,
+        "load_supplied_mw": 0.0
+      },
+      {
+        "steam_enthalpy": 3090.3226202297155,
+        "shaft_power_mw": 336.9056685758782,
+        "electrical_output_mw": 323.4294418328431,
+        "condenser_temperature": 305.0,
+        "load_demand_mw": 0.0,
+        "load_supplied_mw": 0.0
+      },
+      {
+        "steam_enthalpy": 3090.3226202297155,
+        "shaft_power_mw": 336.9056685758782,
+        "electrical_output_mw": 323.4294418328431,
+        "condenser_temperature": 305.0,
+        "load_demand_mw": 0.0,
+        "load_supplied_mw": 0.0
+      }
+    ],
+    "time_elapsed": 600.0
+  },
+  "metadata": {
+    "primary_pump_active": true,
+    "secondary_pump_active": true,
+    "turbine_active": true,
+    "turbine_units": [
+      true,
+      true,
+      true
+    ],
+    "shutdown": false,
+    "consumer": {
+      "online": false,
+      "demand_mw": 800.0,
+      "name": "Grid"
+    }
+  },
+  "health": {
+    "core": {
+      "name": "core",
+      "integrity": 0.9400000000000066,
+      "failed": false
+    },
+    "primary_pump": {
+      "name": "primary_pump",
+      "integrity": 0.5800000000000063,
+      "failed": false
+    },
+    "secondary_pump": {
+      "name": "secondary_pump",
+      "integrity": 0.1120000000000021,
+      "failed": false
+    },
+    "turbine_1": {
+      "name": "turbine_1",
+      "integrity": 0.610000000000003,
+      "failed": false
+    },
+    "turbine_2": {
+      "name": "turbine_2",
+      "integrity": 0.610000000000003,
+      "failed": false
+    },
+    "turbine_3": {
+      "name": "turbine_3",
+      "integrity": 0.610000000000003,
+      "failed": false
+    }
+  }
+}

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/constants.py

@@ -10,10 +10,16 @@ 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
+CONTROL_ROD_SPEED = 0.03  # fraction insertion per second
 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-16  # cm^2, tuned for simulation scale
+# 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/control.py

@@ -27,7 +27,8 @@ class ControlSystem:
         if self.manual_control:
             return self.rod_fraction
         error = (state.power_output_mw - self.setpoint_mw) / self.setpoint_mw
-        adjustment = -error * 0.3
+        # When power is low (negative error) withdraw rods; when high, insert them.
+        adjustment = error * 0.2
         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)

src/reactor_sim/dashboard.py

@@ -9,6 +9,7 @@ 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
@@ -50,8 +51,12 @@ class ReactorDashboard:
             DashboardKey("o", "Toggle secondary pump"),
             DashboardKey("t", "Toggle turbine"),
             DashboardKey("1/2/3", "Toggle turbine units 1-3"),
+            DashboardKey("y/u/i", "Maintain turbine 1/2/3"),
             DashboardKey("c", "Toggle consumer"),
             DashboardKey("r", "Reset & clear state"),
+            DashboardKey("m", "Maintain primary pump"),
+            DashboardKey("n", "Maintain secondary pump"),
+            DashboardKey("k", "Maintain core (requires shutdown)"),
             DashboardKey("+/-", "Withdraw/insert rods"),
             DashboardKey("[/]", "Adjust consumer demand −/+50 MW"),
             DashboardKey("s", "Setpoint −250 MW"),
@@ -142,6 +147,18 @@ class ReactorDashboard:
                 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"))
+            elif ch in (ord("n"), ord("N")):
+                self._queue_command(ReactorCommand.maintain("secondary_pump"))
+            elif ch in (ord("k"), ord("K")):
+                self._queue_command(ReactorCommand.maintain("core"))
+            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:
@@ -240,10 +257,12 @@ class ReactorDashboard:
                 ("Core Power", f"{state.core.power_output_mw:8.1f} MW"),
                 ("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}"),
             ],
         )
+        y = self._draw_section(win, y, "Key Poisons / Emitters", self._poison_lines(state))
         y = self._draw_section(
             win,
             y,
@@ -300,6 +319,7 @@ class ReactorDashboard:
             "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/k/y/u/i to request maintenance (stop equipment first).",
             "Press 'r' to reset/clear state if you want a cold start.",
             "Watch component health to avoid automatic trips.",
         ]
@@ -370,6 +390,24 @@ class ReactorDashboard:
     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) -> tuple[str, str]:
+            qty = inventory.get(symbol, 0.0)
+            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)"))
+        lines.append(fmt("Sm", "Sm (samarium)"))
+        lines.append(fmt("I", "I (iodine)"))
+        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 _draw_health_bar(self, win: "curses._CursesWindow", start_y: int) -> None:
         height, width = win.getmaxyx()
         if start_y >= height - 2:

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)
@@ -57,3 +58,52 @@ class FuelAssembly:
             power_mw,
         )
         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/neutronics.py

@@ -15,18 +15,19 @@ LOGGER = logging.getLogger(__name__)
 def temperature_feedback(temp: float) -> float:
     """Negative coefficient: higher temperature lowers reactivity."""
     reference = 900.0
-    coefficient = -2.5e-5
+    coefficient = -1.5e-5
     return coefficient * (temp - reference)
 
 
 def xenon_poisoning(flux: float) -> float:
-    return min(0.015, 2e-9 * flux)
+    return min(0.01, 5e-10 * flux)
 
 
 @dataclass
 class NeutronDynamics:
     beta_effective: float = 0.0065
     delayed_neutron_fraction: float = 0.0008
+    external_source_coupling: float = 1e-6
 
     def reactivity(self, state: CoreState, control_fraction: float) -> float:
         rho = (
@@ -37,14 +38,15 @@ class NeutronDynamics:
         )
         return rho
 
-    def flux_derivative(self, state: CoreState, rho: float) -> float:
+    def flux_derivative(self, state: CoreState, rho: float, external_source_rate: float = 0.0) -> 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
+        return ((rho - beta) / generation_time) * state.neutron_flux + 1e5 + source_term
 
-    def step(self, state: CoreState, control_fraction: float, dt: float) -> None:
+    def step(self, state: CoreState, control_fraction: float, dt: float, external_source_rate: float = 0.0) -> None:
         rho = self.reactivity(state, control_fraction)
-        d_flux = self.flux_derivative(state, rho)
+        d_flux = self.flux_derivative(state, rho, external_source_rate)
         state.neutron_flux = max(0.0, state.neutron_flux + d_flux * dt)
         state.reactivity_margin = rho
         LOGGER.debug(

src/reactor_sim/reactor.py

@@ -39,6 +39,7 @@ class Reactor:
     turbine_active: bool = True
     turbine_unit_active: list[bool] = field(default_factory=lambda: [True, True, True])
     shutdown: bool = False
+    poison_alerts: set[str] = field(default_factory=set)
 
     def __post_init__(self) -> None:
         if not self.turbines:
@@ -72,6 +73,8 @@ class Reactor:
             reactivity_margin=-0.02,
             power_output_mw=0.1,
             burnup=0.0,
+            fission_product_inventory={},
+            emitted_particles={},
         )
         primary = CoolantLoopState(
             temperature_in=ambient,
@@ -111,15 +114,28 @@ class Reactor:
             overrides = self._apply_command(command, state)
             rod_fraction = overrides.get("rod_fraction", rod_fraction)
 
-        self.neutronics.step(state.core, rod_fraction, dt)
+        decay_power, decay_neutron_source, decay_products, decay_particles = self.fuel.decay_reaction_effects(
+            state.core
+        )
+        self.neutronics.step(state.core, rod_fraction, dt, external_source_rate=decay_neutron_source)
 
         prompt_power, fission_rate, fission_event = self.fuel.prompt_energy_rate(
             state.core.neutron_flux, rod_fraction
         )
-        decay_power = decay_heat_fraction(state.core.burnup) * state.core.power_output_mw
+        decay_heat = decay_heat_fraction(state.core.burnup) * state.core.power_output_mw
-        total_power = prompt_power + decay_power
+        total_power = prompt_power + decay_heat + decay_power
         state.core.power_output_mw = total_power
         state.core.update_burnup(dt)
+        # Track fission products and emitted particles for diagnostics.
+        products: dict[str, float] = {}
+        for atom in fission_event.products:
+            products[atom.symbol] = products.get(atom.symbol, 0.0) + fission_rate * dt
+        for k, v in decay_products.items():
+            products[k] = products.get(k, 0.0) + v * dt
+        particles: dict[str, float] = {k: v * dt for k, v in decay_particles.items()}
+        state.core.add_products(products)
+        state.core.add_emitted_particles(particles)
+        self._check_poison_alerts(state)
 
         pump_demand = overrides.get("coolant_demand", self.control.coolant_demand(state.primary_loop))
         if self.primary_pump_active:
@@ -371,3 +387,11 @@ class Reactor:
                     )
                 )
         return plant
+
+    def _check_poison_alerts(self, state: PlantState) -> None:
+        inventory = state.core.fission_product_inventory or {}
+        for symbol, threshold in constants.KEY_POISON_THRESHOLDS.items():
+            amount = inventory.get(symbol, 0.0)
+            if amount >= threshold and symbol not in self.poison_alerts:
+                self.poison_alerts.add(symbol)
+                LOGGER.warning("Poison level high: %s inventory %.2e exceeds %.2e", symbol, amount, threshold)

src/reactor_sim/simulation.py

@@ -27,6 +27,29 @@ def _default_state_path() -> 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
@@ -113,6 +136,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/state.py

@@ -16,11 +16,21 @@ class CoreState:
     reactivity_margin: float  # delta rho
     power_output_mw: float  # MW thermal
     burnup: float  # fraction of fuel consumed
+    fission_product_inventory: dict[str, float] = field(default_factory=dict)
+    emitted_particles: dict[str, float] = field(default_factory=dict)
 
     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:
@@ -61,6 +71,8 @@ class PlantState:
             "primary_outlet_temp": self.primary_loop.temperature_out,
             "secondary_pressure": self.secondary_loop.pressure,
             "turbine_electric": self.total_electrical_output(),
+            "products": self.core.fission_product_inventory,
+            "particles": self.core.emitted_particles,
         }
 
     def total_electrical_output(self) -> float:
@@ -71,6 +83,9 @@ class PlantState:
 
     @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", {})
         turbines_blob = data.get("turbines")
         if turbines_blob is None:
             # Compatibility with previous single-turbine snapshots.
@@ -78,7 +93,7 @@ class PlantState:
             turbines_blob = [old_turbine] if old_turbine else []
         turbines = [TurbineState(**t) for t in turbines_blob]
         return cls(
-            core=CoreState(**data["core"]),
+            core=CoreState(**core_blob, fission_product_inventory=inventory, emitted_particles=particles),
             primary_loop=CoolantLoopState(**data["primary_loop"]),
             secondary_loop=CoolantLoopState(**data["secondary_loop"]),
             turbines=turbines,

src/reactor_sim/thermal.py

@@ -36,8 +36,12 @@ class ThermalSolver:
     def step_core(self, core: CoreState, primary: CoolantLoopState, power_mw: float, dt: float) -> None:
         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 any power not immediately convected away, and cools toward the primary outlet.
-        core.fuel_temperature = min(core.fuel_temperature, constants.MAX_CORE_TEMPERATURE)
+        heating = 0.005 * max(0.0, power_mw - temp_rise) * dt
+        cooling = 0.05 * max(0.0, core.fuel_temperature - primary.temperature_out) * dt
+        core.fuel_temperature += heating - cooling
+        # Keep fuel temperature bounded and never below the coolant outlet temperature.
+        core.fuel_temperature = min(max(primary.temperature_out, core.fuel_temperature), constants.MAX_CORE_TEMPERATURE)
         LOGGER.debug(
             "Primary loop: flow=%.0f kg/s temp_out=%.1fK core_temp=%.1fK",
             primary.mass_flow_rate,

src/reactor_sim/turbine.py

@@ -25,6 +25,7 @@ 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
 
     def step(
         self,
@@ -37,6 +38,9 @@ class Turbine:
         available_power = max(steam_power_mw, (enthalpy * mass_flow / 1_000.0) / 1_000.0)
         shaft_power_mw = available_power * self.mechanical_efficiency
         electrical = shaft_power_mw * self.generator_efficiency
+        if electrical > self.rated_output_mw:
+            electrical = self.rated_output_mw
+            shaft_power_mw = electrical / max(1e-6, self.generator_efficiency)
         condenser_temp = max(305.0, loop.temperature_in - 20.0)
         state.steam_enthalpy = enthalpy
         state.shaft_power_mw = shaft_power_mw

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_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