Fix two physics bugs the user caught: unbounded field and ignored current limits

1. Iron saturation now in the dynamics (was the "13 Tesla" bug). Flux
   linkage lambda(x,I) = L_air*I + L_iron*overlap(x)*g(I) with
   g(I)=I_sat*tanh(I/I_sat) saturating at the current where iron reaches
   B_sat. Both the circuit back-EMF (dlambda/dx) and the force (coenergy,
   dW'/dx) derive from the SAME lambda, so energy stays conserved (0.025%
   error) AND the field caps at B_sat instead of running to 13 T. Reduces
   to the old 0.5*I^2*dL/dx in the low-current limit. On the config that
   reported 60% efficiency / 11.6 T, it now gives 40% / 1.90 T.

2. Switch surge-current limit is now a hard feasibility constraint: a
   config whose peak discharge current exceeds the switch's pulse rating
   (continuous * surge factor per device kind) is infeasible -- otherwise
   the optimizer "wins" with configs that vaporize their own thyristor
   (e.g. 119 A through a 12 A BT151). Capacitor current stays a warning
   (electrolytics tolerate pulses; DB has continuous not pulse ratings).

Regression tests added for both (field cap near saturation, energy
conservation under strong discharge). 87 tests pass.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>

src/gausse/physics/circuit.py

@@ -1,14 +1,14 @@
 """ОДУ разряда конденсатора через ключ и катушку, связанное с механикой снаряда.
 
 Состояние: [Q, I, x, v] — заряд конденсатора, ток в катушке, положение и
-скорость снаряда. Индуктивность зависит от x (проход снаряда), поэтому
+скорость снаряда. Индуктивность железа насыщается по току, поэтому и член
-уравнение на dI/dt выведено из закона сохранения потокосцепления
+dI/dt, и наведённая ЭДС, и сила выводятся из ОДНОГО потокосцепления
-Ψ = L(x)*I:
+λ(x, I) (см. `inductance.CoilInductanceModel`), что сохраняет энергию:
 
     V_c = Q / C
     dQ/dt = -I
-    dI/dt = (V_c - I*R_total - I*(dL/dx)*v) / L(x)
+    dI/dt = (V_c - I*R_total - (∂λ/∂x)*v) / (∂λ/∂I)
-    dv/dt = F(x, I) / m
+    dv/dt = F(x, I) / m           (F из coenergy того же λ)
     dx/dt = v
 """
 
@@ -16,7 +16,6 @@ from dataclasses import dataclass
 
 import numpy as np
 
-from gausse.physics.force import force_on_slug_newtons
 from gausse.physics.inductance import CoilInductanceModel
 
 Q, I, X, V = range(4)
@@ -27,10 +26,6 @@ class StageCircuitParams:
     capacitance_f: float
     r_total_ohm: float
     mass_kg: float
-    mu_eff: float
-    total_turns: int
-    coil_length_m: float
-    b_sat_tesla: float
 
 
 def derivatives(
@@ -40,13 +35,13 @@ def derivatives(
     params: StageCircuitParams,
 ) -> list[float]:
     q, current, x, v = state
-    l_x = inductance_model.l_of_x(x)
+    dlambda_di = inductance_model.dlambda_di(x, current)
-    dl_dx = inductance_model.dl_dx(x)
+    dlambda_dx = inductance_model.dlambda_dx(x, current)
 
     v_c = q / params.capacitance_f
     d_q = -current
-    d_i = (v_c - current * params.r_total_ohm - current * dl_dx * v) / l_x
+    d_i = (v_c - current * params.r_total_ohm - dlambda_dx * v) / dlambda_di
-    force = force_on_slug_newtons(current_a=current, dl_dx=dl_dx)
+    force = inductance_model.force_newtons(x, current)
     d_v = force / params.mass_kg
     return [d_q, d_i, v, d_v]
 

src/gausse/physics/constants.py

@@ -4,3 +4,14 @@ MU_0 = 4e-7 * 3.141592653589793  # Гн/м, магнитная проницае
 
 RESISTIVITY_COPPER_OHM_M = 1.68e-8  # при ~20°C
 RESISTIVITY_ALUMINUM_OHM_M = 2.82e-8  # при ~20°C
+
+# Импульсный (surge) ток намного выше паспортного среднего/непрерывного,
+# который лежит в базе компонентов. Для одиночного коротко-импульсного
+# выстрела coilgun применяем эти множители к max_current_a детали и
+# отбраковываем конфигурацию, если пиковый ток разряда их превышает —
+# иначе оптимизатор выдаёт "победителей", которые сожгли бы свой ключ.
+# Значения — консервативные типовые (тиристор ITSM ~10× IT(AV), MOSFET
+# импульсный ~4× непрерывного, IGBT ~3×, электролит surge ~5× ripple).
+# Это оценки, а не паспортные surge-рейтинги конкретных деталей.
+SWITCH_SURGE_FACTOR = {"SCR": 10.0, "MOSFET": 4.0, "IGBT": 3.0}
+CAPACITOR_SURGE_FACTOR = 5.0

src/gausse/physics/inductance.py

@@ -12,6 +12,8 @@ from dataclasses import dataclass
 
 import numpy as np
 
+from gausse.physics.constants import MU_0
+
 INCH_M = 0.0254
 
 
@@ -76,12 +78,26 @@ def _sigmoid(z: np.ndarray) -> np.ndarray:
     return 0.5 * (1 + np.tanh(z / 2))
 
 
+def _ln_cosh(z: np.ndarray) -> np.ndarray:
+    """Численно устойчивый ln(cosh z) = |z| + ln((1+e^{-2|z|})/2)."""
+    az = np.abs(z)
+    return az + np.log1p(np.exp(-2 * az)) - np.log(2.0)
+
+
 @dataclass(frozen=True)
 class CoilInductanceModel:
-    """L(x) и аналитическая dL/dx для катушки с проходящим ферромагнитным снарядом.
+    """Потокосцепление λ(x, I) катушки с ферромагнитным снарядом и насыщением.
 
-    x — положение центра снаряда относительно центра катушки (м), катушка
+    x — положение центра снаряда относительно центра катушки (м).
-    считается центрированной в x=0 в своей локальной системе координат.
+
+    Вклад железа в потокосцепление насыщается по току:
+        λ(x, I) = L_air·I + L_iron·overlap(x)·g(I),
+        g(I) = I_sat·tanh(I/I_sat),  L_iron = L_air·(μ_eff−1)
+    где I_sat — ток, при котором поле в железе достигает B_sat. И сила, и
+    наведённая ЭДС в контуре выводятся из ОДНОГО и того же λ (метод
+    coenergy), поэтому энергия сохраняется, а поле не улетает за B_sat.
+    В пределе малого тока (g(I)≈I) всё сводится к линейной модели
+    F = 0.5·I²·dL/dx.
     """
 
     l_air_h: float
@@ -89,6 +105,20 @@ class CoilInductanceModel:
     slug_length_m: float
     mu_eff: float
     smoothing_width_m: float
+    total_turns: int = 0
+    b_sat_tesla: float = 1e9  # по умолчанию насыщение отключено (для геом. тестов)
+
+    @property
+    def l_iron_coeff(self) -> float:
+        return self.l_air_h * (self.mu_eff - 1.0)
+
+    @property
+    def saturation_current_a(self) -> float:
+        """Ток, при котором поле в железе достигает B_sat. inf = без насыщения."""
+        if self.total_turns <= 0 or self.mu_eff <= 1.0 or self.b_sat_tesla >= 1e8:
+            return float("inf")
+        turns_per_m = self.total_turns / self.coil_length_m
+        return self.b_sat_tesla / (MU_0 * self.mu_eff * turns_per_m)
 
     def _edges(self) -> tuple[float, float]:
         half_span = (self.coil_length_m + self.slug_length_m) / 2
@@ -106,8 +136,45 @@ class CoilInductanceModel:
         s2 = _sigmoid((e2 - x) / w)
         return (s1 * s2 / w) * (s2 - s1)
 
+    # --- насыщающиеся токовые функции: g, g', G (интеграл g) ---
+    def _g(self, i: np.ndarray) -> np.ndarray:
+        i_sat = self.saturation_current_a
+        if not np.isfinite(i_sat):
+            return np.asarray(i, dtype=float)
+        return i_sat * np.tanh(i / i_sat)
+
+    def _g_prime(self, i: np.ndarray) -> np.ndarray:
+        i_sat = self.saturation_current_a
+        if not np.isfinite(i_sat):
+            return np.ones_like(np.asarray(i, dtype=float))
+        return 1.0 - np.tanh(i / i_sat) ** 2  # sech^2, устойчиво
+
+    def _g_integral(self, i: np.ndarray) -> np.ndarray:
+        i_sat = self.saturation_current_a
+        if not np.isfinite(i_sat):
+            return 0.5 * np.asarray(i, dtype=float) ** 2
+        return i_sat**2 * _ln_cosh(i / i_sat)
+
+    # --- линейные (ненасыщенные) L(x), dL/dx — для геометрических тестов ---
     def l_of_x(self, x: np.ndarray) -> np.ndarray:
         return self.l_air_h * (1 + (self.mu_eff - 1) * self.overlap_fraction(x))
 
     def dl_dx(self, x: np.ndarray) -> np.ndarray:
-        return self.l_air_h * (self.mu_eff - 1) * self.d_overlap_dx(x)
+        return self.l_iron_coeff * self.d_overlap_dx(x)
+
+    # --- динамика с насыщением (используется контуром) ---
+    def dlambda_di(self, x: np.ndarray, i: np.ndarray) -> np.ndarray:
+        """Дифференциальная индуктивность ∂λ/∂I для члена dI/dt."""
+        return self.l_air_h + self.l_iron_coeff * self.overlap_fraction(x) * self._g_prime(i)
+
+    def dlambda_dx(self, x: np.ndarray, i: np.ndarray) -> np.ndarray:
+        """Пространственный член ∂λ/∂x для наведённой ЭДС."""
+        return self.l_iron_coeff * self.d_overlap_dx(x) * self._g(i)
+
+    def force_newtons(self, x: np.ndarray, i: np.ndarray) -> np.ndarray:
+        """Сила из coenergy: F = ∂/∂x ∫₀ᴵ λ dI' = L_iron·overlap'(x)·G(I)."""
+        return self.l_iron_coeff * self.d_overlap_dx(x) * self._g_integral(i)
+
+    def effective_field_current_a(self, i: float) -> float:
+        """'Насыщенный' эффективный ток для оценки поля (поле ~ g(I), кап B_sat)."""
+        return float(self._g(np.asarray(float(i))))

src/gausse/sim/stage.py

@@ -28,6 +28,7 @@ from gausse.physics.inductance import (
     effective_permeability,
     winding_geometry,
 )
+from gausse.physics.constants import SWITCH_SURGE_FACTOR
 from gausse.physics.force import saturation_scale, solenoid_field_estimate_tesla
 from gausse.physics.sensors import (
     inductive_signal_peak_v,
@@ -144,6 +145,8 @@ def run_stage(
         slug_length_m=projectile.length_m,
         mu_eff=mu_eff,
         smoothing_width_m=wire_od_m,
+        total_turns=geometry.total_turns,
+        b_sat_tesla=projectile.material.b_sat_tesla,
     )
 
     x_sensor_m = -stage.sensor_to_coil_distance_m
@@ -193,10 +196,6 @@ def run_stage(
         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(
@@ -231,11 +230,33 @@ def run_stage(
     # материала снаряда — численные КПД/скорость в этом случае, вероятно,
     # завышены относительно реального железа.
     peak_current_a = float(np.max(np.abs(discharge_sol.y[1])))
+    # поле оцениваем по НАСЫЩЕННОМУ эффективному току g(I) — теперь оно не
+    # улетает за B_sat, а упирается в него (см. inductance.py).
+    effective_field_current_a = inductance_model.effective_field_current_a(peak_current_a)
     peak_field_estimate_tesla = solenoid_field_estimate_tesla(
-        mu_eff, geometry.total_turns, geometry.coil_length_m, peak_current_a
+        mu_eff, geometry.total_turns, geometry.coil_length_m, effective_field_current_a
     )
-    saturation_warning = (
+    # предупреждение теперь означает реальный заход в насыщение (ток > I_sat)
-        saturation_scale(peak_field_estimate_tesla, projectile.material.b_sat_tesla) < 1.0
+    saturation_warning = peak_current_a > inductance_model.saturation_current_a
+
+    # Жёсткая отбраковка по импульсному току: полупроводниковый ключ реально
+    # выходит из строя за пределами surge-рейтинга, а электролит перегревается.
+    # Конфигурация, где пиковый ток разряда это превышает, — не пушка, а
+    # уничтоженная деталь на первом выстреле, поэтому feasible=False.
+    # Ключ — жёсткий предел: полупроводник реально выходит из строя за surge.
+    # Батарея конденсаторов — только предупреждение (см. build_detail):
+    # электролиты переносят высокий импульсный ток, а в базе — паспортный
+    # непрерывный/ripple рейтинг, поэтому жёстко браковать по нему нечестно.
+    switch_surge_a = stage.switch.max_current_a * SWITCH_SURGE_FACTOR.get(stage.switch.kind, 4.0)
+    if peak_current_a > switch_surge_a:
+        return StageResult(
+            feasible=False,
+            reason=(
+                f"пиковый ток {peak_current_a:.0f} А превышает импульсный предел ключа "
+                f"{stage.switch.part_number} ({switch_surge_a:.0f} А = {stage.switch.max_current_a:.0f} А × surge)"
+            ),
+            t_sensor_s=t_sensor_s,
+            t_fire_s=t_sensor_s + fire_delay_s,
         )
 
     return StageResult(

tests/test_circuit_rlc_analytical.py

@@ -28,15 +28,12 @@ def _make_frozen_model() -> CoilInductanceModel:
 
 
 def _integrate(r_ohm: float, t_span: tuple[float, float], t_eval: np.ndarray):
+    # mu_eff=1.0 -> вклад железа нулевой, λ = L_air·I, чистый RLC-контур
     model = _make_frozen_model()
     params = StageCircuitParams(
         capacitance_f=C_F,
         r_total_ohm=r_ohm,
         mass_kg=1.0,
-        mu_eff=1.0,
-        total_turns=100,
-        coil_length_m=0.05,
-        b_sat_tesla=1e9,  # эффективно отключает клэмп насыщения для этого теста
     )
     state0 = [Q0, 0.0, 0.0, 0.0]
     return solve_ivp(

tests/test_report.py

@@ -17,14 +17,15 @@ DB = ComponentDatabase.load()
 def _feasible_config() -> CoilgunConfig:
     wire = next(w for w in DB.wires if w.part_id == "cu-petv2-1.0mm")
     capacitor = next(c for c in DB.capacitors if c.part_number == "cap-470uf-400v")
-    switch = next(s for s in DB.switches if s.part_number == "IRFP250PBF")
+    # мощный ключ (180А, surge ~720А) — конфигурация с реальным пиковым током
+    switch = next(s for s in DB.switches if s.part_number == "IRFP4468PBF")
     sensor = next(s for s in DB.sensors if s.part_number == "TCST2103")
     steel = next(m for m in DB.projectile_materials if m.name == "Ст3 (конструкционная сталь)")
 
     projectile = ProjectileConfig(material=steel, diameter_m=0.008, length_m=0.02)
     stage = StageConfig(
         wire=wire, capacitor=capacitor, switch=switch, sensor=sensor, tube_od_m=0.01,
-        turns_per_layer=20, layers=4, sensor_to_coil_distance_m=0.02, charge_voltage_v=350.0,
+        turns_per_layer=20, layers=4, sensor_to_coil_distance_m=0.02, charge_voltage_v=200.0,
     )
     return CoilgunConfig(
         stages=[stage, stage], inter_stage_gaps_m=[0.05], projectile=projectile,

tests/test_stage_energy_conservation.py

@@ -30,8 +30,9 @@ def _wire(resistivity_ohm_m: float = 1.68e-8) -> WireSpec:
 
 
 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,
+        part_number="sw1", kind="MOSFET", max_current_a=500.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",
     )
@@ -87,29 +88,42 @@ def test_kinetic_energy_increases_for_approaching_slug():
     result = _run(_wire(), _switch())
     assert result.feasible, result.reason
     assert result.kinetic_energy_delta_j > 0
-    # рост кинетической энергии тут должен означать разгон ВПЕРЁД, а не
+    # рост кинетической энергии тут должен означать разгон ВПЕРЁД
-    # отброс назад (см. test_overpowered_discharge_can_fling_slug_backward)
     assert result.exit_v_mps > 0
 
 
-def test_overpowered_discharge_can_fling_slug_backward():
+def test_energy_conserved_under_strong_discharge():
-    """Честная фиксация реального режима отказа, а не только "успешных" случаев.
+    """Энергобаланс должен оставаться точным даже при сильном разряде большой
-
+    батареей (нелинейное насыщение железа не должно рвать сохранение энергии —
-    Слишком большая ёмкость (медленный разряд относительно скорости снаряда)
+    и сила, и ЭДС выведены из одного λ(x, I), см. inductance.py).
-    держит ток высоким уже после того, как снаряд проходит центр катушки —
-    зона x>0 тянет назад (dL/dx<0), и снаряд может выйти с отрицательной
-    скоростью. Энергобаланс при этом всё ещё точен (КЭ ~ v^2 не зависит от
-    знака) — это подтверждает, что находка отражает реальную физику плохо
-    настроенного расстояния/ёмкости, а не ошибку в уравнениях.
     """
-    overpowered_capacitor = CapacitorSpec(
+    big_capacitor = CapacitorSpec(
         part_number="c-big", capacitance_uf=1000.0, voltage_v=400.0, esr_ohm=0.05,
         max_current_a=300.0, price=300.0, source="test",
     )
-    result = _run(_wire(), _switch(), capacitor=overpowered_capacitor)
+    result = _run(_wire(), _switch(), capacitor=big_capacitor)
     assert result.feasible, result.reason
-    assert result.exit_v_mps < 0
     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_field_is_capped_near_saturation():
+    """Регрессия на баг '13 Тесла': поле не должно улетать за B_sat железа."""
+    strong_capacitor = CapacitorSpec(
+        part_number="c-strong", capacitance_uf=200.0, voltage_v=600.0, esr_ohm=0.05,
+        max_current_a=500.0, price=400.0, source="test",
+    )
+    high_v_switch = SwitchSpec(
+        part_number="sw-hi", kind="MOSFET", max_current_a=1000.0, max_voltage_v=800.0,
+        on_resistance_ohm=0.01, on_voltage_drop_v=None, turn_on_time_ns=50.0, price=50.0, source="test",
+    )
+    stage = StageConfig(
+        wire=_wire(), capacitor=strong_capacitor, switch=high_v_switch, sensor=OPTICAL_SENSOR,
+        tube_od_m=0.01, turns_per_layer=40, layers=5, sensor_to_coil_distance_m=0.02, charge_voltage_v=550.0,
+    )
+    result = run_stage(entry_x_m=-0.05, entry_v_mps=5.0, stage=stage, projectile=PROJECTILE)
+    assert result.feasible, result.reason
+    # поле упирается в B_sat снаряда (1.8 Тл в фикстуре), а не улетает в десятки Тл
+    assert result.peak_field_estimate_tesla <= STEEL.b_sat_tesla * 1.05