diff --git a/src/sepflows/absorption/co2_capture.py b/src/sepflows/absorption/co2_capture.py
index 4fcd30e..1c5e2d0 100644
--- a/src/sepflows/absorption/co2_capture.py
+++ b/src/sepflows/absorption/co2_capture.py
@@ -182,10 +182,7 @@ def solve(self) -> CO2CaptureResult:
             * (self._c_amine * self._mw_mea + (1 - self._c_amine) * self._mw_water)
             / 1000.0
         )  # kJ/(mol·K)
-        lrhx_duty_kw = (
-            self._eps_hx * amine_flow * cp_kj_mol_k * delta_t_hx / 3600.0
-        ) * 1000.0  # kJ/h → kW  (× 1000/3600)
-        # Actually: kJ/mol·K × mol/h × K ÷ 3.6 → kW
+        # Q̇ = ε × ṁ_amine × Cp × ΔT.  Units: kJ/(mol·K) × mol/h × K / 3.6 → kW.
         lrhx_duty_kw = self._eps_hx * amine_flow * cp_kj_mol_k * delta_t_hx / 3.6
 
         # ── Absorber height estimate ──────────────────────────────────────────
diff --git a/src/sepflows/asu/cryogenic.py b/src/sepflows/asu/cryogenic.py
index 9b1010e..1bee9d8 100644
--- a/src/sepflows/asu/cryogenic.py
+++ b/src/sepflows/asu/cryogenic.py
@@ -288,13 +288,20 @@ def _boiling_temp_at_pressure(component: str, pressure_bar: float) -> float:
         from sepflows.utils.thermodynamics import antoine_pressure
 
         tb_normal = NORMAL_BOILING_POINTS.get(component, 90.0)  # K at 1 atm
+        if pressure_bar <= 0.0:
+            raise ValueError(f"pressure_bar must be positive, got {pressure_bar}")
         # Iterate: find T such that P_ant(T) = pressure_bar × 1e5
         target_pa = pressure_bar * 1e5
         t = tb_normal
         for _ in range(50):
             p = antoine_pressure(component, t)
             ratio = target_pa / max(p, 1.0)
-            t *= ratio**0.15
             if abs(ratio - 1.0) < 1e-6:
                 break
+            # Clamp the update factor so a pathological Antoine extrapolation
+            # can't drive t towards zero (antoine_pressure raises on t<=0) or
+            # off to thousands of kelvin in a single step.
+            step = ratio**0.15
+            step = min(max(step, 0.5), 2.0)
+            t = max(t * step, 1.0)
         return t
diff --git a/src/sepflows/distillation/rigorous.py b/src/sepflows/distillation/rigorous.py
index 7d20801..adece81 100644
--- a/src/sepflows/distillation/rigorous.py
+++ b/src/sepflows/distillation/rigorous.py
@@ -380,8 +380,13 @@ def _thomas_solve(
         # Forward sweep
         c_prime = np.zeros(n)
         d_prime = np.zeros(n)
-        c_prime[0] = c_diag[0] / b_diag[0]
-        d_prime[0] = d_rhs[0] / b_diag[0]
+        # Guard against a zero pivot on the first row.  This can occur in
+        # pathological initialisations (e.g. v_flows[0] = 0 during a cold
+        # start or a degenerate K-value), and without a guard the solver
+        # silently produces inf/NaN that corrupts the whole column profile.
+        b0 = b_diag[0] if abs(b_diag[0]) > 1e-30 else -1e-30
+        c_prime[0] = c_diag[0] / b0
+        d_prime[0] = d_rhs[0] / b0
         for j in range(1, n):
             denom = b_diag[j] - a_diag[j] * c_prime[j - 1]
             if abs(denom) < 1e-30:
diff --git a/src/sepflows/distillation/shortcut.py b/src/sepflows/distillation/shortcut.py
index 0a2a650..a3ca92c 100644
--- a/src/sepflows/distillation/shortcut.py
+++ b/src/sepflows/distillation/shortcut.py
@@ -249,9 +249,19 @@ def _gilliland(self, n_min: float, r_min: float, r: float) -> float:
         where X = (R − Rmin)/(R + 1) and Y = (N − Nmin)/(N + 1).
         """
         x_g = (r - r_min) / (r + 1.0)
+        # The Molokanov fit contains x_g**0.5 in the denominator and is
+        # only defined for x_g > 0 (i.e. r strictly greater than r_min).
+        # Callers already bump r to r_min + 0.01 but r_min can be negative
+        # for near-ideal splits, pushing x_g onto or below zero.  Clamp to
+        # a small positive value so the correlation stays finite and the
+        # resulting N stays a real, monotone function of the inputs.
+        x_g = max(x_g, 1e-6)
         # Molokanov equation
         y_g = 1.0 - math.exp((1.0 + 54.4 * x_g) / (11.0 + 117.2 * x_g) * (x_g - 1.0) / (x_g**0.5))
-        n = (y_g + n_min) / (1.0 - y_g)
+        denom = 1.0 - y_g
+        if denom <= 0.0:
+            denom = 1e-9
+        n = (y_g + n_min) / denom
         return n
 
     def _kirkbride_feed_stage(
diff --git a/src/sepflows/utils/thermodynamics.py b/src/sepflows/utils/thermodynamics.py
index 91d6a08..0e13060 100644
--- a/src/sepflows/utils/thermodynamics.py
+++ b/src/sepflows/utils/thermodynamics.py
@@ -244,6 +244,18 @@ def _residual(t: float) -> float:
             t_hi += 20.0
             f_hi = _residual(t_hi)
 
+    if f_lo * f_hi > 0:
+        # Bracket expansion failed to produce a sign change.  Bisecting
+        # anyway would silently return whichever endpoint happened to
+        # win the comparison, masking the real problem (e.g. operating
+        # pressure above the critical pressure of every component, or
+        # an Antoine extrapolation that can't reach P_op in [150, 700] K).
+        raise RuntimeError(
+            "Bubble-point bracket failed to bracket a root in "
+            f"[{t_lo:.1f}, {t_hi:.1f}] K at P={pressure_pa:.0f} Pa. "
+            "Check feed composition and pressure, or provide a better t_init_k."
+        )
+
     # Bisection
     for iteration in range(max_iter):
         t_mid = (t_lo + t_hi) / 2.0
@@ -304,6 +316,13 @@ def _residual(t: float) -> float:
             t_hi += 20.0
             f_hi = _residual(t_hi)
 
+    if f_lo * f_hi > 0:
+        raise RuntimeError(
+            "Dew-point bracket failed to bracket a root in "
+            f"[{t_lo:.1f}, {t_hi:.1f}] K at P={pressure_pa:.0f} Pa. "
+            "Check vapour composition and pressure, or provide a better t_init_k."
+        )
+
     for iteration in range(max_iter):
         t_mid = (t_lo + t_hi) / 2.0
         f_mid = _residual(t_mid)