Add boron reactivity trim and power measurement smoothing

src/reactor_sim/reactor.py

@@ -64,6 +64,8 @@ class Reactor:
         # Balance-of-plant controls
         self.feedwater_valve = 0.5
         self._last_steam_out_kg_s = 0.0
+        # Slow chemistry/boron trim control
+        self._boron_trim_active = True
         if not self.primary_pump_units or len(self.primary_pump_units) != 2:
             self.primary_pump_units = [True, True]
         if not self.secondary_pump_units or len(self.secondary_pump_units) != 2:
@@ -197,6 +199,9 @@ class Reactor:
             state.core
         )
         self.neutronics.update_poisons(state.core, dt)
+        # Apply soluble boron reactivity bias (slow trim).
+        boron_delta = state.boron_ppm - constants.CHEM_BORON_DEFAULT_PPM
+        self.neutronics.shutdown_bias = self.neutronics.base_shutdown_bias - boron_delta * constants.BORON_WORTH_PER_PPM
         self.neutronics.step(state.core, rod_fraction, dt, external_source_rate=decay_neutron_source, rod_banks=self.control.rod_banks)
 
         prompt_power, fission_rate, fission_event = self.fuel.prompt_energy_rate(
@@ -431,6 +436,7 @@ class Reactor:
         self._apply_secondary_boiloff(state, dt)
         self._update_secondary_level(state, dt)
         self._update_chemistry(state, dt)
+        self._apply_boron_trim(state, dt)
 
         steam_draw = self._step_turbine_bank(state, transferred, dt)
         if steam_draw > 0.0:
@@ -658,6 +664,20 @@ class Reactor:
         loop.level = min(1.2, max(0.0, loop.inventory_kg / nominal_mass))
         self._last_steam_out_kg_s = steam_out
 
+    def _apply_boron_trim(self, state: PlantState, dt: float) -> None:
+        """Slow soluble boron trim to hold power near setpoint; acts only near target."""
+        if not self._boron_trim_active or self.control.manual_control or self.shutdown:
+            return
+        if state.time_elapsed < 300.0:
+            return
+        if self.control.setpoint_mw <= 0.0:
+            return
+        error = (state.core.power_output_mw - self.control.setpoint_mw) / self.control.setpoint_mw
+        if abs(error) < 0.02:
+            return
+        delta = constants.BORON_TRIM_RATE_PPM_PER_S * error * dt
+        state.boron_ppm = min(constants.CHEM_MAX_PPM, max(0.0, state.boron_ppm + delta))
+
     def _update_chemistry(self, state: PlantState, dt: float) -> None:
         """Track dissolved species and fouling impacts on HX and condenser."""
         env = constants.ENVIRONMENT_TEMPERATURE