Add delayed kinetics and steam-drum balance

src/reactor_sim/neutronics.py

@@ -26,7 +26,7 @@ def xenon_poisoning(flux: float) -> float:
 @dataclass
 class NeutronDynamics:
     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
     shutdown_bias: float = -0.014
     iodine_yield: float = 1e-6  # inventory units per MW*s
@@ -55,12 +55,18 @@ class NeutronDynamics:
         return rho
 
     def flux_derivative(
-        self, state: CoreState, rho: float, external_source_rate: float = 0.0, baseline_source: float = 1e5
+        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
         source_term = self.external_source_coupling * external_source_rate
-        return ((rho - beta) / generation_time) * state.neutron_flux + baseline_source + source_term
+        prompt = ((rho - beta) / generation_time) * state.neutron_flux
+        return prompt + delayed_source + baseline_source + source_term
 
     def step(
         self,
@@ -77,9 +83,22 @@ class NeutronDynamics:
             rho = min(rho, -0.04)
         baseline = 0.0 if shutdown else 1e5
         source = 0.0 if shutdown else external_source_rate
-        d_flux = self.flux_derivative(state, rho, source, baseline_source=baseline)
+        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,
@@ -102,3 +121,38 @@ class NeutronDynamics:
 
     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