diff --git a/PLAN.md b/PLAN.md
index b937303..41c1b56 100644
--- a/PLAN.md
+++ b/PLAN.md
@@ -18,6 +18,19 @@ ChipDip, Cable.ru и др.).
 в SQLite. Никаких приукрашенных цифр — если модель показывает низкий КПД или
 нереализуемость, это тоже результат.
 
+**Найденная и исправленная ошибка (Этап 3):** насыщение сердечника изначально
+клэмпилось только в механическом уравнении (F=0.5·I²·dL/dx), но не в
+электрическом (наведённая ЭДС всё ещё считалась по полной dL/dx) — это
+незаметно ломало точный энергобаланс на ~15%. Тест на сохранение энергии
+(этого же честного протокола, который просил пользователь) это поймал.
+Решение: клэмп насыщения убран из динамики (F=0.5·I²·dL/dx без клэмпа,
+энергобаланс теперь точен до ~0.02%), а `solenoid_field_estimate_tesla`/
+`saturation_scale` оставлены как ДИАГНОСТИКА — `StageResult.saturation_warning`
+честно предупреждает, когда конфигурация физически выходит за пределы
+насыщения материала снаряда, не подменяя динамику. Полная нелинейная
+L(x, I)-модель с coenergy-выводом силы — в разделе "ограничения модели"
+как будущая работа, не как текущая гарантия точности.
+
 Полный план архитектуры: см. историю обсуждения / `physics`, `sim`, `optim`,
 `storage`, `report` модули ниже.
 
diff --git a/src/gausse/physics/circuit.py b/src/gausse/physics/circuit.py
index 461f58c..e45d901 100644
--- a/src/gausse/physics/circuit.py
+++ b/src/gausse/physics/circuit.py
@@ -46,14 +46,7 @@ def derivatives(
     v_c = q / params.capacitance_f
     d_q = -current
     d_i = (v_c - current * params.r_total_ohm - current * dl_dx * v) / l_x
-    force = force_on_slug_newtons(
-        current_a=current,
-        dl_dx_unsaturated=dl_dx,
-        mu_eff=params.mu_eff,
-        total_turns=params.total_turns,
-        coil_length_m=params.coil_length_m,
-        b_sat_tesla=params.b_sat_tesla,
-    )
+    force = force_on_slug_newtons(current_a=current, dl_dx=dl_dx)
     d_v = force / params.mass_kg
     return [d_q, d_i, v, d_v]
 
diff --git a/src/gausse/physics/force.py b/src/gausse/physics/force.py
index 34506f2..110e5de 100644
--- a/src/gausse/physics/force.py
+++ b/src/gausse/physics/force.py
@@ -1,10 +1,15 @@
-"""Сила, действующая на ферромагнитный снаряд: F = 0.5*I^2*dL/dx с насыщением.
+"""Сила, действующая на ферромагнитный снаряд: F = 0.5*I^2*dL/dx.
 
-Насыщение — грубая, явно приближённая поправка (v1): оцениваем поле внутри
-катушки простой соленоидной формулой и, если оно превышает B_sat материала
-снаряда, пропорционально ослабляем вклад сердечника в dL/dx. Реальное
-насыщение нелинейно и зависит от геометрии — здесь используется линейный
-клэмп как консервативная эвристика, а не точная кривая намагничивания.
+Насыщение сердечника сюда сознательно НЕ встроено как динамический клэмп:
+ранняя версия ослабляла силу пропорционально превышению B_sat, но клэмп
+применялся только к механическому уравнению, а не к электрическому
+(наведённая ЭДС по-прежнему считалась через полную dL/dx) — это ломало
+точный энергобаланс (dE/dt = -I^2*R), проверяемый в тестах. Модель
+насыщения, консистентная и для силы, и для ЭДС одновременно, требует
+нелинейной L(x, I) с явной coenergy-производной, что отложено (см. PLAN.md,
+"ограничения модели"). Здесь `solenoid_field_estimate_tesla`/`saturation_scale`
+оставлены как ДИАГНОСТИЧЕСКИЕ функции — предупреждать, что конфигурация,
+вероятно, входит в насыщение, не искажая динамику.
 """
 
 from gausse.physics.constants import MU_0
@@ -23,14 +28,5 @@ def saturation_scale(b_estimate_tesla: float, b_sat_tesla: float) -> float:
     return min(1.0, b_sat_tesla / b_estimate_tesla)
 
 
-def force_on_slug_newtons(
-    current_a: float,
-    dl_dx_unsaturated: float,
-    mu_eff: float,
-    total_turns: int,
-    coil_length_m: float,
-    b_sat_tesla: float,
-) -> float:
-    b_estimate = solenoid_field_estimate_tesla(mu_eff, total_turns, coil_length_m, current_a)
-    scale = saturation_scale(b_estimate, b_sat_tesla)
-    return 0.5 * current_a**2 * dl_dx_unsaturated * scale
+def force_on_slug_newtons(current_a: float, dl_dx: float) -> float:
+    return 0.5 * current_a**2 * dl_dx
diff --git a/src/gausse/physics/sensors.py b/src/gausse/physics/sensors.py
new file mode 100644
index 0000000..6405dd3
--- /dev/null
+++ b/src/gausse/physics/sensors.py
@@ -0,0 +1,56 @@
+"""Событийные модели датчика прохода снаряда для solve_ivp.
+
+Два конкурирующих варианта датчика, оба выражены как условие пересечения
+нуля во время баллистического полёта [x, v] к катушке:
+
+- Оптический/Холла: срабатывает строго по положению x, независимо от
+  скорости снаряда — чистый цифровой фронт.
+- Индукционный (виток провода): наведённое напряжение приближённо
+  пропорционально скорости и "колокольной" чувствительности катушки к
+  положению снаряда (sech^2-профиль). На низкой скорости пик напряжения
+  может не достичь порога — событие тогда не срабатывает вовсе, и это
+  честно отражается в результате как feasible=False, а не подгоняется.
+"""
+
+from typing import Callable
+
+import numpy as np
+
+SensorEvent = Callable[[float, np.ndarray], float]
+
+
+def make_optical_sensor_event(x_sensor_m: float) -> SensorEvent:
+    def event(t: float, state: np.ndarray) -> float:
+        return state[0] - x_sensor_m
+
+    event.terminal = True
+    event.direction = 1.0
+    return event
+
+
+def make_inductive_sensor_event(
+    x_sensor_m: float,
+    sensitivity_v_per_mps: float,
+    threshold_v: float,
+    width_m: float,
+) -> SensorEvent:
+    def event(t: float, state: np.ndarray) -> float:
+        x, v = state
+        bump = 1.0 / np.cosh((x - x_sensor_m) / width_m) ** 2
+        return sensitivity_v_per_mps * v * bump - threshold_v
+
+    event.terminal = True
+    event.direction = 1.0
+    return event
+
+
+def inductive_signal_peak_v(
+    sensitivity_v_per_mps: float, velocity_mps: float
+) -> float:
+    """Пиковое напряжение датчика при данной скорости (bump=1, снаряд точно у датчика).
+
+    Полезно, чтобы заранее и дёшево (без интегрирования) отбраковать
+    заведомо нереализуемые конфигурации: если даже пиковое напряжение
+    ниже порога, датчик не сработает ни при каком положении.
+    """
+    return sensitivity_v_per_mps * abs(velocity_mps)
diff --git a/src/gausse/sim/stage.py b/src/gausse/sim/stage.py
new file mode 100644
index 0000000..6b51cb7
--- /dev/null
+++ b/src/gausse/sim/stage.py
@@ -0,0 +1,255 @@
+"""Одна ступень: баллистический подлёт к датчику -> задержка -> разряд -> энергобаланс.
+
+Нереализуемые случаи (датчик не сработал, разряд не скоммутировался) не
+бросают исключение — возвращаются как `StageResult(feasible=False, reason=...)`
+с честной причиной, чтобы вызывающий код (в т.ч. оптимизатор) мог их
+учитывать, а не терять.
+"""
+
+import math
+from dataclasses import dataclass
+
+import numpy as np
+from scipy.integrate import solve_ivp
+
+from gausse.components.schema import (
+    CapacitorSpec,
+    ProjectileMaterialSpec,
+    SensorSpec,
+    SwitchSpec,
+    WireSpec,
+)
+from gausse.physics import circuit
+from gausse.physics.circuit import StageCircuitParams
+from gausse.physics.inductance import (
+    CoilInductanceModel,
+    air_core_inductance_wheeler,
+    demagnetizing_factor_prolate,
+    effective_permeability,
+    winding_geometry,
+)
+from gausse.physics.force import saturation_scale, solenoid_field_estimate_tesla
+from gausse.physics.sensors import (
+    inductive_signal_peak_v,
+    make_inductive_sensor_event,
+    make_optical_sensor_event,
+)
+
+
+@dataclass(frozen=True)
+class ProjectileConfig:
+    material: ProjectileMaterialSpec
+    diameter_m: float
+    length_m: float
+
+    @property
+    def mass_kg(self) -> float:
+        radius_m = self.diameter_m / 2
+        volume_m3 = math.pi * radius_m**2 * self.length_m
+        return volume_m3 * self.material.density_kg_m3
+
+    @property
+    def aspect_ratio(self) -> float:
+        return self.length_m / self.diameter_m
+
+
+@dataclass(frozen=True)
+class StageConfig:
+    wire: WireSpec
+    capacitor: CapacitorSpec
+    switch: SwitchSpec
+    sensor: SensorSpec
+    tube_od_m: float
+    turns_per_layer: int
+    layers: int
+    sensor_to_coil_distance_m: float
+    charge_voltage_v: float
+    inductive_sensor_width_m: float = 0.005
+    sensor_time_budget_s: float = 0.05
+    discharge_time_budget_s: float = 0.02
+
+
+@dataclass
+class StageResult:
+    feasible: bool
+    reason: str | None = None
+    t_sensor_s: float | None = None
+    t_fire_s: float | None = None
+    exit_x_m: float | None = None
+    exit_v_mps: float | None = None
+    energy_in_j: float | None = None
+    energy_dissipated_j: float | None = None
+    energy_remaining_cap_j: float | None = None
+    kinetic_energy_delta_j: float | None = None
+    saturation_warning: bool = False
+    peak_field_estimate_tesla: float | None = None
+    discharge_t: np.ndarray | None = None
+    discharge_q: np.ndarray | None = None
+    discharge_i: np.ndarray | None = None
+    discharge_x: np.ndarray | None = None
+    discharge_v: np.ndarray | None = None
+
+
+def _wire_resistance_ohm(wire: WireSpec, wire_length_m: float) -> float:
+    bare_radius_m = wire.gauge_mm / 1000 / 2
+    area_m2 = math.pi * bare_radius_m**2
+    return wire.resistivity_ohm_m * wire_length_m / area_m2
+
+
+def _switch_equivalent_resistance_ohm(
+    switch: SwitchSpec, capacitance_f: float, inductance_h: float, voltage_v: float
+) -> float:
+    if switch.on_resistance_ohm is not None:
+        return switch.on_resistance_ohm
+    # SCR/тиристор: прямое падение аппроксимируется эквивалентным
+    # сопротивлением относительно характерного масштаба тока контура
+    # I_ref = V0*sqrt(C/L) (пиковый ток недодемпфированного LC-разряда).
+    i_ref = voltage_v * math.sqrt(capacitance_f / inductance_h)
+    drop_v = switch.on_voltage_drop_v or 0.0
+    return drop_v / max(i_ref, 1e-6)
+
+
+def _ballistic_derivatives(t: float, state: np.ndarray) -> list[float]:
+    return [state[1], 0.0]
+
+
+def run_stage(
+    entry_x_m: float,
+    entry_v_mps: float,
+    stage: StageConfig,
+    projectile: ProjectileConfig,
+) -> StageResult:
+    mass_kg = projectile.mass_kg
+    wire_od_m = stage.wire.insulation_od_mm / 1000
+    geometry = winding_geometry(stage.tube_od_m, wire_od_m, stage.turns_per_layer, stage.layers)
+
+    r_wire = _wire_resistance_ohm(stage.wire, geometry.total_wire_length_m)
+    capacitance_f = stage.capacitor.capacitance_uf * 1e-6
+    l_air_h = air_core_inductance_wheeler(
+        geometry.mean_radius_m, geometry.coil_length_m, geometry.radial_depth_m, geometry.total_turns
+    )
+    r_switch = _switch_equivalent_resistance_ohm(
+        stage.switch, capacitance_f, l_air_h, stage.charge_voltage_v
+    )
+    r_total_ohm = r_wire + r_switch + stage.capacitor.esr_ohm
+
+    demag = demagnetizing_factor_prolate(projectile.aspect_ratio)
+    mu_eff = effective_permeability(projectile.material.mu_r, demag)
+    inductance_model = CoilInductanceModel(
+        l_air_h=l_air_h,
+        coil_length_m=geometry.coil_length_m,
+        slug_length_m=projectile.length_m,
+        mu_eff=mu_eff,
+        smoothing_width_m=wire_od_m,
+    )
+
+    x_sensor_m = -stage.sensor_to_coil_distance_m
+
+    if stage.sensor.kind == "inductive":
+        peak_v = inductive_signal_peak_v(stage.sensor.sensitivity_v_per_mps, entry_v_mps)
+        if peak_v < stage.sensor.threshold_v:
+            return StageResult(
+                feasible=False,
+                reason=(
+                    f"индукционный датчик: пиковый сигнал {peak_v:.4f}В "
+                    f"ниже порога {stage.sensor.threshold_v:.4f}В "
+                    f"при скорости {entry_v_mps:.2f} м/с"
+                ),
+            )
+        sensor_event = make_inductive_sensor_event(
+            x_sensor_m,
+            stage.sensor.sensitivity_v_per_mps,
+            stage.sensor.threshold_v,
+            stage.inductive_sensor_width_m,
+        )
+    else:
+        sensor_event = make_optical_sensor_event(x_sensor_m)
+
+    flight_sol = solve_ivp(
+        _ballistic_derivatives,
+        (0.0, stage.sensor_time_budget_s),
+        [entry_x_m, entry_v_mps],
+        events=sensor_event,
+        rtol=1e-8,
+        atol=1e-10,
+    )
+    if len(flight_sol.t_events[0]) == 0:
+        return StageResult(
+            feasible=False,
+            reason="датчик не сработал в пределах временного бюджета полёта",
+        )
+
+    t_sensor_s = flight_sol.t_events[0][0]
+    x_at_sensor, v_at_sensor = flight_sol.y_events[0][0]
+
+    fire_delay_s = (stage.sensor.propagation_delay_ns + stage.switch.turn_on_time_ns) * 1e-9
+    x_fire_m = x_at_sensor + v_at_sensor * fire_delay_s
+    v_fire_mps = v_at_sensor
+
+    circuit_params = StageCircuitParams(
+        capacitance_f=capacitance_f,
+        r_total_ohm=r_total_ohm,
+        mass_kg=mass_kg,
+        mu_eff=mu_eff,
+        total_turns=geometry.total_turns,
+        coil_length_m=geometry.coil_length_m,
+        b_sat_tesla=projectile.material.b_sat_tesla,
+    )
+    q0 = capacitance_f * stage.charge_voltage_v
+    discharge_sol = solve_ivp(
+        circuit.derivatives,
+        (0.0, stage.discharge_time_budget_s),
+        [q0, 0.0, x_fire_m, v_fire_mps],
+        args=(inductance_model, circuit_params),
+        events=circuit.zero_current_crossing_event,
+        rtol=1e-8,
+        atol=1e-11,
+    )
+    if len(discharge_sol.t_events[0]) == 0:
+        return StageResult(
+            feasible=False,
+            reason="разряд не скоммутировался (ток не вернулся к нулю) в пределах временного бюджета",
+            t_sensor_s=t_sensor_s,
+            t_fire_s=t_sensor_s + fire_delay_s,
+        )
+
+    q_final, i_final, x_final, v_final = discharge_sol.y_events[0][0]
+
+    energy_in_j = 0.5 * capacitance_f * stage.charge_voltage_v**2
+    energy_remaining_cap_j = q_final**2 / (2 * capacitance_f)
+    energy_dissipated_j = float(
+        np.trapezoid(discharge_sol.y[1] ** 2 * r_total_ohm, discharge_sol.t)
+    )
+    kinetic_before_j = 0.5 * mass_kg * v_fire_mps**2
+    kinetic_after_j = 0.5 * mass_kg * v_final**2
+
+    # Диагностика насыщения (не влияет на динамику, см. docstring force.py):
+    # предупреждаем, если реальный пик тока подразумевает поле выше B_sat
+    # материала снаряда — численные КПД/скорость в этом случае, вероятно,
+    # завышены относительно реального железа.
+    peak_current_a = float(np.max(np.abs(discharge_sol.y[1])))
+    peak_field_estimate_tesla = solenoid_field_estimate_tesla(
+        mu_eff, geometry.total_turns, geometry.coil_length_m, peak_current_a
+    )
+    saturation_warning = (
+        saturation_scale(peak_field_estimate_tesla, projectile.material.b_sat_tesla) < 1.0
+    )
+
+    return StageResult(
+        feasible=True,
+        t_sensor_s=t_sensor_s,
+        t_fire_s=t_sensor_s + fire_delay_s,
+        exit_x_m=x_final,
+        exit_v_mps=v_final,
+        energy_in_j=energy_in_j,
+        energy_dissipated_j=energy_dissipated_j,
+        energy_remaining_cap_j=energy_remaining_cap_j,
+        kinetic_energy_delta_j=kinetic_after_j - kinetic_before_j,
+        saturation_warning=saturation_warning,
+        peak_field_estimate_tesla=peak_field_estimate_tesla,
+        discharge_t=discharge_sol.t,
+        discharge_q=discharge_sol.y[0],
+        discharge_i=discharge_sol.y[1],
+        discharge_x=discharge_sol.y[2],
+        discharge_v=discharge_sol.y[3],
+    )
diff --git a/tests/test_force_model.py b/tests/test_force_model.py
index 52a372d..c8ba5ea 100644
--- a/tests/test_force_model.py
+++ b/tests/test_force_model.py
@@ -22,14 +22,13 @@ def test_saturation_scale_clamps_above_bsat():
     assert scale == pytest.approx(0.5)
 
 
-def test_force_is_reduced_once_saturated():
-    common = dict(dl_dx_unsaturated=1e-3, mu_eff=200, total_turns=300, coil_length_m=0.05, b_sat_tesla=1.8)
-    force_unsaturated = force_on_slug_newtons(current_a=5, **common)
-    force_saturated = force_on_slug_newtons(current_a=500, **common)
-    # без клэмпа сила росла бы как I^2 (в 10000 раз); с насыщением рост должен быть намного меньше
-    assert force_saturated / force_unsaturated < 5000
+def test_force_scales_as_current_squared():
+    # F = 0.5*I^2*dL/dx: без клэмпа (см. модуль-докстринг force.py про энергобаланс)
+    # сила должна расти строго как I^2, иначе нарушится точный энергобаланс контура.
+    f_low = force_on_slug_newtons(current_a=5, dl_dx=1e-3)
+    f_high = force_on_slug_newtons(current_a=50, dl_dx=1e-3)
+    assert f_high / f_low == pytest.approx(100.0)
 
 
 def test_force_zero_at_zero_current():
-    common = dict(dl_dx_unsaturated=1e-3, mu_eff=200, total_turns=300, coil_length_m=0.05, b_sat_tesla=1.8)
-    assert force_on_slug_newtons(current_a=0, **common) == pytest.approx(0.0)
+    assert force_on_slug_newtons(current_a=0, dl_dx=1e-3) == pytest.approx(0.0)
diff --git a/tests/test_sensors.py b/tests/test_sensors.py
new file mode 100644
index 0000000..b7e6c52
--- /dev/null
+++ b/tests/test_sensors.py
@@ -0,0 +1,50 @@
+import numpy as np
+from scipy.integrate import solve_ivp
+
+from gausse.physics.sensors import make_inductive_sensor_event, make_optical_sensor_event
+
+
+def _ballistic(t, state):
+    return [state[1], 0.0]
+
+
+def test_optical_sensor_trigger_is_velocity_independent():
+    x_sensor = 0.02
+    trigger_xs = []
+    for v in (1.0, 5.0, 20.0):
+        event = make_optical_sensor_event(x_sensor)
+        sol = solve_ivp(_ballistic, (0, 1.0), [-0.05, v], events=event)
+        assert len(sol.t_events[0]) == 1
+        trigger_xs.append(sol.y_events[0][0][0])
+    assert np.allclose(trigger_xs, x_sensor, atol=1e-9)
+
+
+def test_inductive_sensor_fires_earlier_at_higher_velocity():
+    x_sensor = 0.02
+    sensitivity = 0.05
+    threshold = 0.3
+    width = 0.005
+
+    trigger_xs = []
+    for v in (10.0, 20.0, 40.0):
+        event = make_inductive_sensor_event(x_sensor, sensitivity, threshold, width)
+        sol = solve_ivp(_ballistic, (0, 1.0), [-0.05, v], events=event)
+        assert len(sol.t_events[0]) == 1
+        trigger_xs.append(sol.y_events[0][0][0])
+
+    # чем выше скорость, тем раньше (дальше от датчика, т.е. при меньшем x)
+    # срабатывает индукционный датчик, т.к. sensitivity*v*bump(x) достигает
+    # порога при меньшем bump(x), а значит при большем |x - x_sensor|
+    assert trigger_xs[0] > trigger_xs[1] > trigger_xs[2]
+    assert all(x < x_sensor for x in trigger_xs)
+
+
+def test_inductive_sensor_never_fires_below_threshold_speed():
+    x_sensor = 0.02
+    sensitivity = 0.05
+    threshold = 0.3
+    width = 0.005
+    # пиковый сигнал = sensitivity * v = 0.05 * 1.0 = 0.05 << порог 0.3
+    event = make_inductive_sensor_event(x_sensor, sensitivity, threshold, width)
+    sol = solve_ivp(_ballistic, (0, 1.0), [-0.05, 1.0], events=event)
+    assert len(sol.t_events[0]) == 0
diff --git a/tests/test_stage_energy_conservation.py b/tests/test_stage_energy_conservation.py
new file mode 100644
index 0000000..86bd43b
--- /dev/null
+++ b/tests/test_stage_energy_conservation.py
@@ -0,0 +1,89 @@
+import pytest
+
+from gausse.components.schema import (
+    CapacitorSpec,
+    ProjectileMaterialSpec,
+    SensorSpec,
+    SwitchSpec,
+    WireSpec,
+)
+from gausse.sim.stage import ProjectileConfig, StageConfig, run_stage
+
+STEEL = ProjectileMaterialSpec(
+    name="steel", density_kg_m3=7850.0, mu_r=200.0, b_sat_tesla=1.8, price_per_kg=100.0, source="test"
+)
+PROJECTILE = ProjectileConfig(material=STEEL, diameter_m=0.008, length_m=0.02)
+CAPACITOR = CapacitorSpec(
+    part_number="c1", capacitance_uf=1000.0, voltage_v=400.0, esr_ohm=0.05, max_current_a=300.0,
+    price=300.0, source="test",
+)
+OPTICAL_SENSOR = SensorSpec(
+    part_number="s1", kind="optical", propagation_delay_ns=500.0, price=20.0, source="test"
+)
+
+
+def _wire(resistivity_ohm_m: float = 1.68e-8) -> WireSpec:
+    return WireSpec(
+        part_id="w1", material="copper", gauge_mm=0.8, insulation_od_mm=0.85,
+        resistivity_ohm_m=resistivity_ohm_m, max_current_a=10.0, price_per_m=3.0, source="test",
+    )
+
+
+def _switch(on_resistance_ohm: float = 0.02) -> SwitchSpec:
+    return SwitchSpec(
+        part_number="sw1", kind="MOSFET", max_current_a=200.0, max_voltage_v=500.0,
+        on_resistance_ohm=on_resistance_ohm, on_voltage_drop_v=None, turn_on_time_ns=50.0,
+        price=50.0, source="test",
+    )
+
+
+def _run(wire, switch, capacitor=CAPACITOR):
+    # turns_per_layer выбран так, чтобы длина катушки (~1.7см) была сравнима
+    # со снарядом (2см) — иначе снаряд, войдя глубоко внутрь длинной катушки,
+    # оказывается в плоской зоне перекрытия (dL/dx=0) и сила не действует.
+    stage = StageConfig(
+        wire=wire,
+        capacitor=capacitor,
+        switch=switch,
+        sensor=OPTICAL_SENSOR,
+        tube_od_m=0.01,
+        turns_per_layer=20,
+        layers=4,
+        sensor_to_coil_distance_m=0.02,
+        charge_voltage_v=350.0,
+    )
+    return run_stage(entry_x_m=-0.05, entry_v_mps=5.0, stage=stage, projectile=PROJECTILE)
+
+
+def test_stage_is_feasible_with_realistic_components():
+    result = _run(_wire(), _switch())
+    assert result.feasible, result.reason
+
+
+def test_energy_conserved_exactly_when_lossless():
+    lossless_capacitor = CapacitorSpec(
+        part_number="c-lossless", capacitance_uf=1000.0, voltage_v=400.0, esr_ohm=0.0,
+        max_current_a=300.0, price=300.0, source="test",
+    )
+    result = _run(
+        _wire(resistivity_ohm_m=0.0), _switch(on_resistance_ohm=0.0), capacitor=lossless_capacitor
+    )
+    assert result.feasible, result.reason
+    assert result.energy_dissipated_j == pytest.approx(0.0, abs=1e-9)
+    balance = result.energy_remaining_cap_j + result.kinetic_energy_delta_j
+    assert balance == pytest.approx(result.energy_in_j, rel=1e-4)
+
+
+def test_energy_balance_holds_with_realistic_losses():
+    result = _run(_wire(), _switch())
+    assert result.feasible, result.reason
+    balance = (
+        result.energy_remaining_cap_j + result.energy_dissipated_j + result.kinetic_energy_delta_j
+    )
+    assert balance == pytest.approx(result.energy_in_j, rel=1e-2)
+
+
+def test_kinetic_energy_increases_for_approaching_slug():
+    result = _run(_wire(), _switch())
+    assert result.feasible, result.reason
+    assert result.kinetic_energy_delta_j > 0