diff --git a/docs/src/conventions.md b/docs/src/conventions.md
index f0a97b2777..154afed892 100644
--- a/docs/src/conventions.md
+++ b/docs/src/conventions.md
@@ -379,7 +379,7 @@ Here:
 
 - ``\omega_b`` and ``\omega_D`` are the bounce and magnetic-precession frequency **coefficients at
   thermal energy** (``x = 1``): ``\omega_b \propto v_\mathrm{th}`` and
-  ``\omega_D \propto v_\mathrm{th}^2`` (`wbhat`, `wdhat` in `Torque.jl`). Their velocity-space
+  ``\omega_D \propto v_\mathrm{th}^2`` (`wbbar`, `wdbar` times `bhat`, `dhat` in `BounceAveraging.jl`). Their velocity-space
   dependence enters explicitly as ``\sqrt{x}`` (bounce ``\propto v``) and ``x`` (precession
   ``\propto E``).
 - ``\ell_\mathrm{eff}`` is the effective bounce harmonic (`PitchIntegration.jl`):
diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml
index 2ad60ab72d..fa41a0dd50 100644
--- a/examples/DIIID-like_SLAYER_example/gpec.toml
+++ b/examples/DIIID-like_SLAYER_example/gpec.toml
@@ -18,7 +18,7 @@ psihigh = 0.9995               # Upper limit of normalized flux coordinate
 mpsi = 0                       # Number of radial grid points (0 = auto-compute from psi_accuracy)
 psi_accuracy = 0.001           # Target absolute error in q for auto-mpsi
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 
 # Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning]
diff --git a/examples/DIIID-like_error_field_example/gpec.toml b/examples/DIIID-like_error_field_example/gpec.toml
index c7738d769b..070ffa8172 100644
--- a/examples/DIIID-like_error_field_example/gpec.toml
+++ b/examples/DIIID-like_error_field_example/gpec.toml
@@ -22,7 +22,7 @@ psihigh = 0.995                # Upper limit of normalized poloidal flux (captur
 mpsi = 0                       # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 impurity_charge = 6                     # Impurity charge (carbon)
 impurity_mass = 12                      # Impurity mass [proton masses]
diff --git a/examples/DIIID-like_gal_resistive_example/gpec.toml b/examples/DIIID-like_gal_resistive_example/gpec.toml
index 12d3c8e646..19fc8cfa2f 100644
--- a/examples/DIIID-like_gal_resistive_example/gpec.toml
+++ b/examples/DIIID-like_gal_resistive_example/gpec.toml
@@ -17,7 +17,7 @@ psihigh = 0.993                # Upper limit of normalized flux coordinate (0.99
 mpsi = 128                     # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 
 # Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning]
diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
index 73c500c5ce..d9d48d1c98 100644
--- a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
+++ b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
@@ -18,7 +18,7 @@ psihigh = 0.993                # Upper limit of normalized flux coordinate (0.99
 mpsi = 128                     # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 
 # Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning]
diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml
index 56705328d9..f7feb30261 100644
--- a/examples/DIIID-like_ideal_example/gpec.toml
+++ b/examples/DIIID-like_ideal_example/gpec.toml
@@ -17,7 +17,7 @@ psihigh = 0.995                # Upper limit of normalized poloidal flux (captur
 mpsi = 0                       # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 impurity_charge = 6                     # Impurity charge (carbon)
 impurity_mass = 12                      # Impurity mass [proton masses]
diff --git a/examples/DIIID-like_ideal_multi_n_example/gpec.toml b/examples/DIIID-like_ideal_multi_n_example/gpec.toml
index 8db24007c8..c6bfa903ee 100644
--- a/examples/DIIID-like_ideal_multi_n_example/gpec.toml
+++ b/examples/DIIID-like_ideal_multi_n_example/gpec.toml
@@ -17,7 +17,7 @@ psilow = 1e-4                  # Lower limit of normalized poloidal flux
 psihigh = 0.9                  # Upper limit of normalized poloidal flux
 mpsi = 128                     # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 
 # Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning]
diff --git a/examples/DIIID-like_riccati_deltaprime_example/gpec.toml b/examples/DIIID-like_riccati_deltaprime_example/gpec.toml
index 77f4960279..31457af0d4 100644
--- a/examples/DIIID-like_riccati_deltaprime_example/gpec.toml
+++ b/examples/DIIID-like_riccati_deltaprime_example/gpec.toml
@@ -17,7 +17,7 @@ psihigh = 0.995                # Upper limit of normalized poloidal flux (captur
 mpsi = 0                       # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 
 # Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning]
diff --git a/examples/DIIID-like_truncation_example/gpec.toml b/examples/DIIID-like_truncation_example/gpec.toml
index 90c47c52e5..eb6539369c 100644
--- a/examples/DIIID-like_truncation_example/gpec.toml
+++ b/examples/DIIID-like_truncation_example/gpec.toml
@@ -20,7 +20,7 @@ psihigh = 0.9999               # Upper limit of normalized poloidal flux (past t
 mpsi = 0                       # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 
 [ForceFreeStates]
diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml
index 7cc873c46e..f4fabe17a9 100644
--- a/examples/Solovev_ideal_example/gpec.toml
+++ b/examples/Solovev_ideal_example/gpec.toml
@@ -13,7 +13,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7                # Relative tolerance of the ODE solves that trace each flux surface
 
 [Wall]
diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml
index 8eeab5d6f2..a7b68e8d62 100644
--- a/examples/Solovev_ideal_example_3D/gpec.toml
+++ b/examples/Solovev_ideal_example_3D/gpec.toml
@@ -12,7 +12,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7                # Relative tolerance of the ODE solves that trace each flux surface
 
 # Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning]
diff --git a/examples/Solovev_ideal_example_multi_n/gpec.toml b/examples/Solovev_ideal_example_multi_n/gpec.toml
index 0509be3445..c8e51b8c78 100644
--- a/examples/Solovev_ideal_example_multi_n/gpec.toml
+++ b/examples/Solovev_ideal_example_multi_n/gpec.toml
@@ -13,7 +13,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7                # Relative tolerance of the ODE solves that trace each flux surface
 
 [Wall]
diff --git a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml
index af79fd46ff..0c272a9a0a 100644
--- a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml
+++ b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml
@@ -13,7 +13,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7                # Relative tolerance of the ODE solves that trace each flux surface
 
 [Wall]
diff --git a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml
index e051a232ea..f200df8c2a 100644
--- a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml
+++ b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml
@@ -13,7 +13,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7                # Relative tolerance of the ODE solves that trace each flux surface
 
 [Wall]
diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml
index bd48a5331c..24b0c408f1 100644
--- a/examples/a10_kinetic_example/gpec.toml
+++ b/examples/a10_kinetic_example/gpec.toml
@@ -19,7 +19,7 @@ psilow = 1e-3                  # Lower limit of normalized poloidal flux
 psihigh = 0.99                 # Upper limit of normalized poloidal flux
 mpsi = 16                      # Number of radial grid points (low value for fast iteration)
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7       # Relative tolerance of the ODE solves that trace each flux surface
 impurity_charge = 6            # Impurity charge
 impurity_mass = 12             # Impurity mass
diff --git a/src/Deck.jl b/src/Deck.jl
index 361cd183ad..ea0f4b17ce 100644
--- a/src/Deck.jl
+++ b/src/Deck.jl
@@ -49,6 +49,8 @@ const _MOVED_KEYS = Dict{String,Dict{String,String}}(
         "impurity_charge" => "moved to [Equilibrium]",
         "impurity_mass" => "moved to [Equilibrium]",
         "nn" => "removed; KineticForces takes the toroidal mode number from the run",
+        "rlar_flag" => "removed; the large-aspect-ratio methods are not ported, use fgar_flag or tgar_flag",
+        "clar_flag" => "removed; the large-aspect-ratio methods are not ported, use fgar_flag or tgar_flag",
         "kinetic_file" => _KINETIC_FILE_HOME),
     "SLAYER" => Dict{String,String}(
         "enabled" => "a [Tearing] section with inner_model = \"SLAYER\" runs the tearing analysis",
diff --git a/src/Equilibrium/Axisymmetric/Direct.jl b/src/Equilibrium/Axisymmetric/Direct.jl
index 9eb785ffcf..e0c2b93e7a 100644
--- a/src/Equilibrium/Axisymmetric/Direct.jl
+++ b/src/Equilibrium/Axisymmetric/Direct.jl
@@ -609,17 +609,13 @@ Direct equilibrium driver: locates the axis, traces each flux surface with `fiel
     if q0 <= 0.0
         @warn "q0 extrapolation to axis gives q0 = $(@sprintf("%.3f", q0)) ≤ 0 — likely a spline artifact from psilow being too large; check psilow or use newq0 to override."
     end
-    # newq0 = -1 flips the extrapolated q0; resolved locally so the config keeps the user's request.
     newq0 = equil_params.newq0
-    if newq0 == -1
-        newq0 = -q0
-    end
     if newq0 != 0.0
         @info "Revising q-profile for newq0 = $(@sprintf("%.3f", newq0))"
         f0 = profiles.F_spline.y[1] - profiles.F_deriv(psi_nodes[1]; hint=Ref(1)) * psi_nodes[1]
         f0fac = f0^2 * ((newq0 / q0)^2 - 1.0)
         for i in 1:(mpsi+1)
-            ffac = sqrt(1.0 + f0fac / profiles.F_spline.y[i]^2) * sign(newq0)
+            ffac = sqrt(1.0 + f0fac / profiles.F_spline.y[i]^2)
             sq_nodes[i, 1] *= ffac
             sq_nodes[i, 4] *= ffac
             rzphi_nodes[i, :, 3] .*= ffac
diff --git a/src/Equilibrium/Axisymmetric/Inverse.jl b/src/Equilibrium/Axisymmetric/Inverse.jl
index e8cec7ddb5..440bee1b04 100644
--- a/src/Equilibrium/Axisymmetric/Inverse.jl
+++ b/src/Equilibrium/Axisymmetric/Inverse.jl
@@ -283,16 +283,13 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No
     f1_sq_lo = sq_deriv(sq_xs[1])
     f1_sq_hi = sq_deriv(sq_xs[end])
     q0 = f_sq[1, 4] - f1_sq_lo[4] * sq_xs[1]
-    if newq0 == -1
-        newq0 = -q0
-    end
 
     if newq0 != 0
-        f0 = f_sq[1, 2] - f1_sq_lo[2] * sq_xs[1]
+        f0 = f_sq[1, 1] - f1_sq_lo[1] * sq_xs[1]
         f0fac = f0^2 * ((newq0 / q0)^2 - 1)
         q0 = newq0
         for ipsi in 0:mpsi
-            ffac = sqrt(1 + f0fac / f_sq[ipsi+1, 1]^2) * sign(newq0)
+            ffac = sqrt(1 + f0fac / f_sq[ipsi+1, 1]^2)
             sq_fs[ipsi+1, 1] *= ffac
             sq_fs[ipsi+1, 4] *= ffac
             rzphi_fs[ipsi+1, :, 3] *= ffac
diff --git a/src/Equilibrium/Config.jl b/src/Equilibrium/Config.jl
index ac091a75dd..646c1da2ea 100644
--- a/src/Equilibrium/Config.jl
+++ b/src/Equilibrium/Config.jl
@@ -53,7 +53,7 @@ specified in the input.
     two-pass auto grid (knot count scales as τ^(-1/3))
   - `mtheta::Int` - Number of poloidal grid points
   - `newq0::Float64` - Target on-axis safety factor q(0); the q and F profiles are rescaled to
-    meet it (0 = use input value, -1 = use the axis extrapolation with its sign flipped)
+    meet it (0 = use input value); q is positive, so negative values are rejected
   - `flux_surface_rtol::Float64` - Relative tolerance of the ODE solves that trace each flux surface; their absolute tolerance
     is `equil_abstol(flux_surface_rtol)`, i.e. `flux_surface_rtol` capped at `EQUIL_ABSTOL_MAX`
   - `imas_cocos::Int` - Expected COCOS convention of an IMAS `dd.equilibrium` input (11 = IMAS standard, 2 = GPEC internal)
@@ -174,6 +174,7 @@ specified in the input.
             error("Cannot recognize jac_type = $(jac_type)")
         end
         (isfinite(flux_surface_rtol) && flux_surface_rtol > 0) || error("flux_surface_rtol = $flux_surface_rtol must be a positive, finite relative tolerance")
+        newq0 >= 0 || error("newq0 = $newq0 is not supported: q must be positive. Use newq0 = 0 (keep the input q) or newq0 > 0.")
         if psihigh > 1.0
             @warn "psihigh = $psihigh exceeds 1.0 (separatrix); clamping to 1.0"
         end
diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl
index f7fb5ef392..dd910df689 100644
--- a/src/Equilibrium/KineticProfiles.jl
+++ b/src/Equilibrium/KineticProfiles.jl
@@ -20,20 +20,40 @@ const _KEV_J = 1.602e-16         # 1 keV in J (temperature normalization in the
 const _NU_PREFAC = 3.5e17        # Krook collision-frequency prefactor
 const _MP = 1.672_614e-27        # proton mass [kg]
 const _ME = 9.109_1e-31          # electron mass [kg]
+const _T_FLOOR_KEV = 1.0e-6      # collision-model temperature floor (1 meV), as Fortran PENTRC
 
-# NRL-formulary Coulomb logarithm (natural log); n_e in m⁻³, T_e in J.
-_coulomb_log(ne, Te) = 17.3 - 0.5 * log(ne / 1.0e20) + 1.5 * log(Te / _KEV_J)
+# Temperature in keV for the collision model, floored so a cold edge node keeps a finite, positive ν.
+_collision_tkev(T) = max(T / _KEV_J, _T_FLOOR_KEV)
+
+# NRL-formulary Coulomb logarithm (natural log) at the floored temperature, floored at 1; n_e in m⁻³, T_e in J.
+_coulomb_log(ne, Te) = max(17.3 - 0.5 * log(ne / 1.0e20) + 1.5 * log(_collision_tkev(Te)), 1.0)
+
+# Warn once per kinetic source, as Fortran PENTRC, when a cold node engages the temperature or Coulomb-logarithm floor.
+function _warn_collision_floors(source, ne, Te, Ti)
+    any(T -> T / _KEV_J < _T_FLOOR_KEV, Te) | any(T -> T / _KEV_J < _T_FLOOR_KEV, Ti) &&
+        @warn "$source: flooring T at 1 meV in the collision frequencies" maxlog = 1 _id = Symbol("T_floor ", source)
+    any(i -> 17.3 - 0.5 * log(ne[i] / 1.0e20) + 1.5 * log(_collision_tkev(Te[i])) < 1, eachindex(ne)) &&
+        @warn "$source: flooring the Coulomb logarithm at 1 in the collision frequencies" maxlog = 1 _id = Symbol("lnLambda_floor ", source)
+    return nothing
+end
 
 # Names kinetic profiles in error messages by their provenance, when they carry one.
 _describe(data) = isempty(data.provenance) ? "kinetic profiles" : "kinetic profiles ($(data.provenance))"
 
-# Hard stop on an unphysical negative density (bad input or a cubic-resample overshoot).
-function _assert_nonneg_density(source, label, arr, psi_reg)
+# Hard stop on an unphysical negative density or temperature (bad input or an extrapolation past the table).
+function _assert_nonneg_profile(source, label, arr, psi_reg)
     i = findfirst(<(0), arr)
     i === nothing ||
         error("$source: negative $label = $(arr[i]) at ψ_n = $(psi_reg[i]) — unphysical")
 end
 
+# Hard stop on a non-positive electron or main-ion density, which the collision model divides by.
+function _assert_positive_density(source, label, arr, psi_reg)
+    i = findfirst(<=(0), arr)
+    i === nothing ||
+        error("$source: $label = $(arr[i]) at ψ_n = $(psi_reg[i]) — densities must be positive")
+end
+
 """
     ResolvedNTVSpecies{P}
 
@@ -235,13 +255,16 @@ function resolve_ntv_species(data::KineticProfileData, ion_species::AbstractVect
 
     nkin = _NKIN
     psi_reg = collect(0:nkin) ./ nkin
-    ne = _cubic_resample(psi_in, ne_in, psi_reg)
-    ni_total = _cubic_resample(psi_in, ni_in, psi_reg)
-    Ti = _cubic_resample(psi_in, Ti_in, psi_reg) .* _EV_J
-    Te = _cubic_resample(psi_in, Te_in, psi_reg) .* _EV_J
-    omegaE = _cubic_resample(psi_in, omE_in, psi_reg)
-    _assert_nonneg_density(source, "n_e", ne, psi_reg)
-    _assert_nonneg_density(source, "n_i (total)", ni_total, psi_reg)
+    ne = _pchip_resample(psi_in, ne_in, psi_reg)
+    ni_total = _pchip_resample(psi_in, ni_in, psi_reg)
+    Ti = _pchip_resample(psi_in, Ti_in, psi_reg) .* _EV_J
+    Te = _pchip_resample(psi_in, Te_in, psi_reg) .* _EV_J
+    omegaE = _pchip_resample(psi_in, omE_in, psi_reg)
+    _assert_positive_density(source, "n_e", ne, psi_reg)
+    _assert_positive_density(source, "n_i (total)", ni_total, psi_reg)
+    _assert_nonneg_profile(source, "T_i", Ti, psi_reg)
+    _assert_nonneg_profile(source, "T_e", Te, psi_reg)
+    _warn_collision_floors(source, ne, Te, Ti)
 
     # Fractions are shares of the total n_i: never above 1, and exactly 1 when every species is fraction-specified.
     fracs = [s.fraction for s in ion_species if !isnan(s.fraction)]
@@ -275,11 +298,11 @@ function resolve_ntv_species(data::KineticProfileData, ion_species::AbstractVect
                     "ion_species[$si]: density profile `$(s.density)` not found in the kinetic profiles " *
                     "(available: $(data.species_densities === nothing ? "none — ASCII files carry only n_i" : join(keys(data.species_densities), ", ")))"
                 )
-            ns[si] = _cubic_resample(psi_in, data.species_densities[s.density], psi_reg)
+            ns[si] = _pchip_resample(psi_in, data.species_densities[s.density], psi_reg)
         else
             ns[si] = s.fraction .* ni_total
         end
-        _assert_nonneg_density(source, "ion_species[$si] density", ns[si], psi_reg)
+        _assert_nonneg_profile(source, "ion_species[$si] density", ns[si], psi_reg)
     end
 
     # Shared composition + collisionality (natural-log Coulomb log), per grid point.
@@ -296,18 +319,18 @@ function resolve_ntv_species(data::KineticProfileData, ion_species::AbstractVect
         zpitch[i] = zp
         n_main[i] = nm
         n_imp[i] = nimp
-        loglam[i] = (ne[i] > 0 && Te[i] > 0) ? _coulomb_log(ne[i], Te[i]) : 0.0
-        nue[i] = Te[i] > 0 ? (zp / _NU_PREFAC) * ne[i] * loglam[i] / (sqrt(_ME / _MP) * (Te[i] / _KEV_J)^1.5) : 0.0
+        loglam[i] = _coulomb_log(ne[i], Te[i])
+        nue[i] = (zp / _NU_PREFAC) * ne[i] * loglam[i] / (sqrt(_ME / _MP) * _collision_tkev(Te[i])^1.5)
     end
     # A negative impurity density means the main ions over-neutralize (Σ z_s n_s > n_e) — unphysical.
-    _assert_nonneg_density(source, "impurity density (n_e < Σ z_s n_s)", n_imp, psi_reg)
+    _assert_nonneg_profile(source, "impurity density (n_e < Σ z_s n_s)", n_imp, psi_reg)
     # Approaching the zpitch pole at Zeff = zimp signals a dilute main-ion mix taken up by impurity.
     maximum(zeff) > 0.9 * zimp &&
         @warn "max(Zeff) = $(maximum(zeff)) is within 10% of zimp = $zimp — zpitch closure near its pole; check ion_species densities/fractions"
 
     # Per-species Krook ν (Logan & Park 2013 Eq. 6) with a test-particle z_s²; `zpitch` is a main-ion closure,
     # applied approximately to the impurity/electron test species.
-    _nu(zsp, msp) = [Ti[i] > 0 ? (zpitch[i] / _NU_PREFAC) * zsp^2 * n_main[i] * loglam[i] / (sqrt(Float64(msp)) * (Ti[i] / _KEV_J)^1.5) : 0.0 for i in 1:npts]
+    _nu(zsp, msp) = [(zpitch[i] / _NU_PREFAC) * zsp^2 * n_main[i] * loglam[i] / (sqrt(Float64(msp)) * _collision_tkev(Ti[i])^1.5) for i in 1:npts]
     _view(dens, nu) = KineticProfileSplines(psi_reg, dens, ne, Ti, Te, omegaE, loglam, nu, nue, zeff)
 
     # Full NTV species set: main ions, the neutrality-closing impurity, and (optionally) electrons.
@@ -341,6 +364,16 @@ Whether `data` carries the profiles the NTV view needs (`n_e`, `T_i`, `T_e` and
 """
 has_ntv_fields(data::KineticProfileData) = all(!isnothing, (data.n_e, data.T_i, data.T_e, data.omega_E))
 
+"""
+    ntv_profiles_usable(data::KineticProfileData) -> Bool
+
+Whether `data` has the NTV fields with positive densities and non-negative temperatures at its table points,
+the inputs `kinetic_profile_splines` accepts.
+"""
+ntv_profiles_usable(data::KineticProfileData) =
+    has_ntv_fields(data) && all(>(0), data.n_e) && (data.n_i === nothing || all(>(0), data.n_i)) &&
+    all(>=(0), data.T_i) && all(>=(0), data.T_e)
+
 """
     kinetic_profile_splines(data::KineticProfileData;
                             zi::Int=1, zimp::Int=6, mi::Int=2, mimp::Int=12,
@@ -351,7 +384,7 @@ has_ntv_fields(data::KineticProfileData) = all(!isnothing, (data.n_e, data.T_i,
 
 The NTV view of raw kinetic profiles (e.g. `equil.kinetic`): interpolate onto a regular
 101-point ψ grid, optionally apply profile scaling knobs, derive collisional / Z_eff
-diagnostics, and return a `KineticProfileSplines` with independent named cubic splines.
+diagnostics, and return a `KineticProfileSplines` of named splines (pchip for lnΛ and ν, cubic otherwise).
 
 Needs `n_e`, `T_i`, `T_e` and `omega_E`; `n_i` defaults to `n_e`. See
 [`read_kinetic_file`](@ref) for the file formats.
@@ -383,6 +416,8 @@ function kinetic_profile_splines(data::KineticProfileData;
     ExB_rotation_factor::Float64=1.0, toroidal_rotation_factor::Float64=1.0,
     chi1::Union{Nothing,Float64}=nothing)
 
+    density_factor > 0 && temperature_factor > 0 ||
+        throw(ArgumentError("density_factor = $density_factor and temperature_factor = $temperature_factor must be positive"))
     # n_i defaults to n_e (quasineutrality); omega_tor and the diffusivities are unused here.
     _need(field, name) = field === nothing ? error("$(_describe(data)) missing required '$name'") : field
     psi_input = data.psi
@@ -395,12 +430,16 @@ function kinetic_profile_splines(data::KineticProfileData;
     nkin = _NKIN
     psi_reg = collect(0:nkin) ./ nkin
 
-    # Cubic-spline resample as in Fortran PENTRC; linear interpolation misplaces ω_E zero crossings and so the resonance.
-    ni = _cubic_resample(psi_input, ni_input, psi_reg)
-    ne = _cubic_resample(psi_input, ne_input, psi_reg)
-    Ti = _cubic_resample(psi_input, Ti_input_eV, psi_reg) .* _EV_J
-    Te = _cubic_resample(psi_input, Te_input_eV, psi_reg) .* _EV_J
-    omegaE = _cubic_resample(psi_input, omegaE_input, psi_reg)
+    # Monotone resample as in Fortran PENTRC; linear interpolation misplaces ω_E zero crossings and so the resonance.
+    ni = _pchip_resample(psi_input, ni_input, psi_reg)
+    ne = _pchip_resample(psi_input, ne_input, psi_reg)
+    Ti = _pchip_resample(psi_input, Ti_input_eV, psi_reg) .* _EV_J
+    Te = _pchip_resample(psi_input, Te_input_eV, psi_reg) .* _EV_J
+    omegaE = _pchip_resample(psi_input, omegaE_input, psi_reg)
+    _assert_positive_density(_describe(data), "n_e", ne, psi_reg)
+    _assert_positive_density(_describe(data), "n_i", ni, psi_reg)
+    _assert_nonneg_profile(_describe(data), "T_i", Ti, psi_reg)
+    _assert_nonneg_profile(_describe(data), "T_e", Te, psi_reg)
 
     needs_rotation_reform = density_factor != 1.0 || temperature_factor != 1.0 || toroidal_rotation_factor != 1.0
     any_scaling = needs_rotation_reform || ExB_rotation_factor != 1.0
@@ -442,6 +481,8 @@ function kinetic_profile_splines(data::KineticProfileData;
     if ExB_rotation_factor != 1.0
         omegaE .*= ExB_rotation_factor
     end
+    # After scaling, so a small temperature_factor that engages the floor is reported.
+    _warn_collision_floors(_describe(data), ne, Te, Ti)
 
     loglam = zeros(Float64, nkin + 1)
     nui = zeros(Float64, nkin + 1)
@@ -454,17 +495,15 @@ function kinetic_profile_splines(data::KineticProfileData;
         T_i = Ti[i]
         T_e = Te[i]
 
-        z = n_e > 0 ? zimp - (n_i / n_e) * zi * (zimp - zi) : Float64(zimp)
+        z = zimp - (n_i / n_e) * zi * (zimp - zi)
         zpitch = 1.0 + (1.0 + mimp) / (2.0 * mimp) * zimp * (z - 1.0) / (zimp - z)
 
         ll = _coulomb_log(n_e, T_e)
         loglam[i] = ll
 
         # Krook ν (Logan & Park 2013 Eq. 6) with a test-particle z²; the single-species limit of resolve_ntv_species.
-        nui[i] = T_i > 0 ?
-                 (zpitch / _NU_PREFAC) * zi^2 * n_i * ll / (sqrt(1.0 * mi) * (T_i / _KEV_J)^1.5) : 0.0
-        nue[i] = T_e > 0 ?
-                 (zpitch / _NU_PREFAC) * n_e * ll / (sqrt(_ME / _MP) * (T_e / _KEV_J)^1.5) : 0.0
+        nui[i] = (zpitch / _NU_PREFAC) * zi^2 * n_i * ll / (sqrt(1.0 * mi) * _collision_tkev(T_i)^1.5)
+        nue[i] = (zpitch / _NU_PREFAC) * n_e * ll / (sqrt(_ME / _MP) * _collision_tkev(T_e)^1.5)
         zeff[i] = z
     end
 
@@ -496,11 +535,11 @@ function _read_kinetic_table(kinetic_file::AbstractString)
 end
 
 """
-Cubic spline of `y(x)` on the irregular input grid, evaluated at `x_new` with smooth cubic
-extrapolation, as Fortran PENTRC `spline_fit(..., "extrap")`.
+Monotone (pchip) fit of `y(x)` on the irregular input grid, evaluated at `x_new` with
+extrapolation, as Fortran PENTRC `spline_fit_pchip`: no overshoot between table points.
 """
-function _cubic_resample(x::AbstractVector, y::AbstractVector, x_new::AbstractVector)
-    spl = cubic_interp(collect(Float64, x), collect(Float64, y); extrap=ExtendExtrap())
+function _pchip_resample(x::AbstractVector, y::AbstractVector, x_new::AbstractVector)
+    spl = pchip_interp(collect(Float64, x), collect(Float64, y); extrap=ExtendExtrap())
     return [spl(xv) for xv in x_new]
 end
 
diff --git a/src/Equilibrium/PlasmaEquilibrium.jl b/src/Equilibrium/PlasmaEquilibrium.jl
index fed3b79963..dd1737bf66 100644
--- a/src/Equilibrium/PlasmaEquilibrium.jl
+++ b/src/Equilibrium/PlasmaEquilibrium.jl
@@ -222,9 +222,10 @@ KineticProfileData(data::KineticProfileData; kwargs...) =
 """
     KineticProfileSplines
 
-Named 1D cubic spline interpolants for kinetic profiles (densities, temperatures, ExB rotation,
+Named 1D spline interpolants for kinetic profiles (densities, temperatures, ExB rotation,
 collisional diagnostics) on the regular kinetic ψ grid, with derivative views of the densities and
-temperatures for the diamagnetic frequencies of the NTV kernel.
+temperatures for the diamagnetic frequencies of the NTV kernel. The Coulomb logarithm and collision
+frequencies are monotone (pchip) interpolants, so a steep cold edge cannot ring below zero.
 
 # Fields
 
@@ -237,7 +238,7 @@ temperatures for the diamagnetic frequencies of the NTV kernel.
   - `zeff_spline`: Effective charge Z_eff
   - `ni_deriv`, `ne_deriv`, `Ti_deriv`, `Te_deriv`: Derivative views for ψ-derivative access
 """
-struct KineticProfileSplines{S,D}
+struct KineticProfileSplines{S,D,C}
     xs::Vector{Float64}
     npts::Int
     npts_minus_1::Int
@@ -246,9 +247,9 @@ struct KineticProfileSplines{S,D}
     Ti_spline::S
     Te_spline::S
     omegaE_spline::S
-    loglam_spline::S
-    nui_spline::S
-    nue_spline::S
+    loglam_spline::C
+    nui_spline::C
+    nue_spline::C
     zeff_spline::S
     ni_deriv::D
     ne_deriv::D
@@ -261,7 +262,8 @@ end
                           extrap=ExtendExtrap())
 
 Build the kinetic profile splines from arrays of values defined on the shared
-ψ grid `xs`. Uses CubicFit boundary conditions to match `ProfileSplines`.
+ψ grid `xs`. Uses CubicFit boundary conditions to match `ProfileSplines`, and pchip for the
+Coulomb logarithm and collision frequencies as Fortran PENTRC.
 Temperatures must already be in Joules.
 """
 function KineticProfileSplines(xs::Vector{Float64},
@@ -278,9 +280,9 @@ function KineticProfileSplines(xs::Vector{Float64},
     Ti_spline = cubic_interp(xs, Ti; extrap=extrap)
     Te_spline = cubic_interp(xs, Te; extrap=extrap)
     omegaE_spline = cubic_interp(xs, omegaE; extrap=extrap)
-    loglam_spline = cubic_interp(xs, loglam; extrap=extrap)
-    nui_spline = cubic_interp(xs, nui; extrap=extrap)
-    nue_spline = cubic_interp(xs, nue; extrap=extrap)
+    loglam_spline = pchip_interp(xs, loglam; extrap=extrap)
+    nui_spline = pchip_interp(xs, nui; extrap=extrap)
+    nue_spline = pchip_interp(xs, nue; extrap=extrap)
     zeff_spline = cubic_interp(xs, zeff; extrap=extrap)
 
     ni_deriv = deriv1(ni_spline)
@@ -288,7 +290,7 @@ function KineticProfileSplines(xs::Vector{Float64},
     Ti_deriv = deriv1(Ti_spline)
     Te_deriv = deriv1(Te_spline)
 
-    KineticProfileSplines{typeof(ni_spline),typeof(ni_deriv)}(
+    KineticProfileSplines{typeof(ni_spline),typeof(ni_deriv),typeof(nui_spline)}(
         xs, npts, npts - 1,
         ni_spline, ne_spline, Ti_spline, Te_spline,
         omegaE_spline, loglam_spline, nui_spline, nue_spline, zeff_spline,
diff --git a/src/Equilibrium/Sources/Chease.jl b/src/Equilibrium/Sources/Chease.jl
index 734981e641..085532955f 100644
--- a/src/Equilibrium/Sources/Chease.jl
+++ b/src/Equilibrium/Sources/Chease.jl
@@ -1,3 +1,7 @@
+# CHEASE supplies q directly and GPEC requires q > 0, so a non-positive q stops here, as Fortran's chease reader does.
+_chease_q(zq, config) =
+    minimum(zq) > 0 ? zq : error("CHEASE input $(config.eq_filename) has q = $(minimum(zq)) ≤ 0; GPEC requires q > 0")
+
 """
     read_chease_binary(config) -> InverseRunInput
 
@@ -40,7 +44,7 @@ function read_chease_binary(config::EquilibriumConfig)
         fs = zeros(npsi1, 4)
         fs[:, 1] .= ztmf
         fs[:, 2] .= zcppr
-        fs[:, 3] .= zq
+        fs[:, 3] .= _chease_q(zq, config)
 
         # Pressure from the cumulative integral of the P′ column.
         itp_pressure = cubic_interp(xs, fs[:, 2])
@@ -205,7 +209,7 @@ function read_chease_ascii(config::EquilibriumConfig)
     # ztmf equals zq .* zfb (both are F).
     fs[:, 1] .= ztmf
     fs[:, 2] .= zcppr # normalized Pressure
-    fs[:, 3] .= zq # q profile
+    fs[:, 3] .= _chease_q(zq, config) # q profile
     # Pressure from the cumulative integral of the P′ column.
     itp_pressure = cubic_interp(xs, fs[:, 2])
     fsi_pressure = FastInterpolations.cumulative_integrate(itp_pressure)
diff --git a/src/KineticForces/BounceAveraging.jl b/src/KineticForces/BounceAveraging.jl
index b3a54c3aac..c8e6975295 100644
--- a/src/KineticForces/BounceAveraging.jl
+++ b/src/KineticForces/BounceAveraging.jl
@@ -834,8 +834,8 @@ bounce point, so a monotone interval between stationary points of B contains one
 `bo/λ` lies strictly between its endpoint B values; the crossing is then solved from the
 cell's cubic, with `hints` carrying each interval's cell over from the previous λ.
 
-Returns roots sorted descending, the order the deepest-well and marginally-trapped logic
-downstream assumes; fewer than two roots signals a degenerate λ and sends the caller to
+Returns roots sorted ascending, as Fortran `spline_roots`, so the deepest-well search pairs neighbouring
+roots; fewer than two roots signals a degenerate λ and sends the caller to
 the fallback. `hints` is sized `length(bf.theta) + 1`, the extra slot being the far half
 of the interval that wraps through the seam.
 """
@@ -870,7 +870,7 @@ function _bounce_points_at_lambda!(
         push!(bpts, mod(root, 1.0))
     end
 
-    sort!(bpts; rev=true)
+    sort!(bpts)
     return bpts
 end
 
@@ -885,14 +885,14 @@ function _find_bounce_points_and_grid(
     ntheta::Int, bf::SurfaceBField, bpts_buf::Vector{Float64}, hints::Vector{Int}
 )
     if sigma == 0  # trapped
-        # Bounce points: roots of v_par(θ) in (0,1), sorted descending as the well logic below assumes.
+        # Bounce points: roots of v_par(θ) in (0,1), sorted ascending so neighbours bound one well or barrier.
         bpts = _bounce_points_at_lambda!(bpts_buf, hints, bf, lmda, bo)
 
         if length(bpts) < 2
             # Degenerate λ (v_par tangent to zero at an extremum, or < 2 stationary points of B): rare,
             # so fall back to an adaptive whole-interval root scan.
             vpar_fn = θ -> _vpar_from_spline(B_vpar, lmda, bo, θ)
-            bpts = sort!(Roots.find_zeros(vpar_fn, 0.0, 1.0); rev=true)
+            bpts = sort!(Roots.find_zeros(vpar_fn, 0.0, 1.0))
         end
 
         nbpts = length(bpts)
@@ -1023,10 +1023,8 @@ function _bounce_samples!(
                 continue
             else
                 # After midpoint: hold the previous sample to the end of the grid.
-                # The wd slot is deliberately held from the wb integrand (g_wb),
-                # reproducing the Fortran behavior for parity.
                 fill!(view(g_wb, i:ntheta), g_wb[i-1])
-                fill!(view(g_wd, i:ntheta), g_wb[i-1])
+                fill!(view(g_wd, i:ntheta), g_wd[i-1])
                 fill!(view(jvtheta, i:ntheta), jvtheta[i-1])
                 break
             end
diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl
index 710aae592e..7ffea77ce2 100644
--- a/src/KineticForces/CalculatedKineticMatrices.jl
+++ b/src/KineticForces/CalculatedKineticMatrices.jl
@@ -171,14 +171,15 @@ end
 
 # The force-balance solve's two-pass grid inputs: the kinetic profiles as the KineticForces model views them (its
 # control when the source is one, else the default control) and the kinetic-resonance surfaces to pin. The model
-# needs the NTV view, so a missing profile errors; otherwise profiles read for another consumer (e.g. SLAYER without
-# omega_E) are left out and the grid follows the equilibrium alone.
+# needs the NTV view, so an unusable profile errors; otherwise profiles read for another consumer (e.g. SLAYER without
+# omega_E, or with a zero edge density) are left out and the grid follows the equilibrium alone.
 function ForceFreeStates.kinetic_grid_inputs(source, ctrl, equil, intr; rationals)
     kf_ctrl = source isa KineticForcesControl ? source : KineticForcesControl()
     data = equil.kinetic
-    ntv = source isa KineticForcesControl || (data !== nothing && Equilibrium.has_ntv_fields(data))
+    ntv = source isa KineticForcesControl || (data !== nothing && Equilibrium.ntv_profiles_usable(data))
     data !== nothing && !ntv &&
-        @info "Two-pass psi grid: the kinetic profiles lack an NTV field (n_e, T_i, T_e, omega_E), so the grid follows the equilibrium alone"
+        @info "Two-pass psi grid: the kinetic profiles lack an NTV field (n_e, T_i, T_e, omega_E) or have a non-positive density " *
+              "or negative temperature, so the grid follows the equilibrium alone"
     kinetic_profiles = ntv ? ntv_profiles(kf_ctrl, equil) : nothing
     species = ntv ? ntv_species(kf_ctrl, equil) : nothing
     pinned = resonance_grid_nodes(ctrl, kf_ctrl, kinetic_profiles, species, equil, intr; rationals)
diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl
index 0c806e64c0..dd087dc941 100644
--- a/src/KineticForces/Compute.jl
+++ b/src/KineticForces/Compute.jl
@@ -34,10 +34,12 @@ end
 
 Warn when the ψ torque quadrature terminated without satisfying the requested tolerances
 (hit `maxevals_psi`), or when a nonzero user-set `atol_psi` dominated termination — the
-silent-garbage scenario for weak applied fields, since NTV scales as δB².
+silent-garbage scenario for weak applied fields, since NTV scales as δB². Error when the
+torque or its error estimate is not finite, which the tolerance comparisons cannot catch.
 """
 function check_psi_quadrature_convergence(total::ComplexF64, quad_err::Float64,
     ctrl::KineticForcesControl, method::String)
+    isfinite(total) && isfinite(quad_err) || error("ψ torque quadrature ($method) returned a non-finite torque $total (error estimate $quad_err).")
     if quad_err > max(ctrl.atol_psi, ctrl.rtol_psi * abs(total))
         @warn "ψ torque quadrature ($method) did not converge within maxevals_psi=$(ctrl.maxevals_psi): " *
               "error estimate $quad_err vs |T|=$(abs(total)) N·m. Raise maxevals_psi or loosen rtol_psi."
diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl
index 7a91500f09..33cadb66be 100644
--- a/src/KineticForces/KineticForcesStructs.jl
+++ b/src/KineticForces/KineticForcesStructs.jl
@@ -7,8 +7,7 @@ Single source of truth for the NTV calculation methods. Each entry is a NamedTup
   - `name`  — short method identifier used as the HDF5 group key and in `intr.method`
   - `flag`  — the `KineticForcesControl` field symbol that enables the method
   - `kind`  — dispatch routing tag consumed by `method_kind` / `Torque.jl`
-    (`:gar` for the GAR/matrix family, `:fcgl`/`:rlar`/`:clar` for the
-    three special-cased methods)
+    (`:gar` for the GAR/matrix family, `:fcgl` for the fluid CGL method)
   - `doc`   — one-line description printed in verbose output
 
 To add a method: append an entry here and add the matching `*_flag` field
@@ -18,8 +17,6 @@ const METHOD_REGISTRY = (
     (name="fgar", flag=:fgar_flag, kind=:gar, doc="Full general-aspect-ratio calculation"),
     (name="tgar", flag=:tgar_flag, kind=:gar, doc="Trapped particle general-aspect-ratio calculation"),
     (name="pgar", flag=:pgar_flag, kind=:gar, doc="Passing particle general-aspect-ratio calculation"),
-    (name="rlar", flag=:rlar_flag, kind=:rlar, doc="Trapped particle large-aspect-ratio calculation"),
-    (name="clar", flag=:clar_flag, kind=:clar, doc="Trapped particle cylindrical large-aspect-ratio calculation"),
     (name="fcgl", flag=:fcgl_flag, kind=:fcgl, doc="Fluid Chew–Goldberger–Low calculation"),
     (name="fwmm", flag=:fwmm_flag, kind=:gar, doc="Full energy calculation using MXM Euler–Lagrange matrix"),
     (name="twmm", flag=:twmm_flag, kind=:gar, doc="Trapped energy calculation using MXM Euler–Lagrange matrix"),
@@ -37,7 +34,7 @@ const METHOD_REGISTRY = (
 """
     method_kind(name) -> Symbol
 
-Return the dispatch routing tag (`:gar`, `:fcgl`, `:rlar`, `:clar`) for the NTV
+Return the dispatch routing tag (`:gar` or `:fcgl`) for the NTV
 method `name`, looked up from [`METHOD_REGISTRY`]. Errors on an unknown method.
 """
 function method_kind(name::AbstractString)
@@ -80,8 +77,6 @@ Immutable: vary a field by building a new control.
     fgar_flag::Bool = true          # Full general-aspect-ratio
     tgar_flag::Bool = false         # Trapped particle GAR
     pgar_flag::Bool = false         # Passing particle GAR
-    rlar_flag::Bool = false         # Trapped particle large-aspect-ratio
-    clar_flag::Bool = false         # Trapped cylindrical LAR
     fcgl_flag::Bool = false         # Fluid Chew-Goldberger-Low
     fwmm_flag::Bool = false         # Full energy via MXM E-L matrix
     twmm_flag::Bool = false         # Trapped energy via MXM E-L matrix
diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl
index ff3b38aaf1..82e6102ac5 100644
--- a/src/KineticForces/Torque.jl
+++ b/src/KineticForces/Torque.jl
@@ -16,7 +16,7 @@ Toroidal torque from nonambipolar transport at one ψ and bounce harmonic `l`, w
   - `wdfac::Float64`: Drift factor
   - `divxfac::Float64`: Divergence factor
   - `electron::Bool`: Calculate quantities for electrons (zi,mi ignored)
-  - `method::String`: Integration method (RLAR, CLAR, *GAR, *TMM, *WMM, *KMM)
+  - `method::String`: Integration method (FCGL, *GAR, *TMM, *WMM, *KMM)
     where * = F,T,P for full,trapped,passing
   - `equil`: PlasmaEquilibrium with 2D interpolants and named profile/geometry splines
   - `intr::KineticForcesInternal`: Internal state with mode indexing and perturbation splines
@@ -58,10 +58,7 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int,
            _torque_surface(psi, zi, mi, electron, method, equil, intr, kinetic_profiles;
         need_matrices=(!isnothing(op_wmats)), n, wdfac, divxfac) :
            surface
-    (; chrg, mass, tspl, bmax, bmin, q, n_s, T_s, welec, wdian, wdiat, wtran, wgyro, nuk, epsr, dVdpsi, dbob_m_f, divx_m_f) = surf
-    wbhat = (π / 4) * sqrt(epsr / 2) * wtran
-    wdhat = q^3 * wtran^2 / (4 * epsr * wgyro) * wdfac
-    nueff = nuk / (2 * epsr)
+    (; chrg, mass, tspl, bmax, bmin, q, n_s, T_s, welec, wdian, wdiat, nuk, epsr, dbob_m_f, divx_m_f) = surf
 
     # Method selection — route on the registry dispatch tag (errors on unknown method)
     kind = method_kind(method)
@@ -69,16 +66,6 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int,
         tpsi_var[] = calculate_fcgl(psi, n, l, tspl, dbob_m_f, divx_m_f, divxfac,
             n_s, T_s, equil, intr)
 
-    elseif kind == :rlar
-        tpsi_var[] = calculate_rlar(psi, n, l, q, epsr, wdian, wdiat, welec,
-            wdhat, wbhat, nueff, dVdpsi, n_s, T_s,
-            dbob_m_f, intr.bo, bmin)
-
-    elseif kind == :clar
-        tpsi_var[] = calculate_clar(psi, n, l, q, epsr, wdian, wdiat, welec,
-            nuk, intr.bo, bmax, bmin, n_s, T_s, mass, chrg,
-            tspl, dbob_m_f, divx_m_f, divxfac, wdfac)
-
     else  # :gar — fgar/tgar/pgar + all *mm matrix methods
         tpsi_var[] = calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec,
             nuk, intr.bo, bmax, bmin, n_s, T_s, mass, chrg,
@@ -167,63 +154,6 @@ function calculate_fcgl(psi, n, l, tspl, dbob_m_f, divx_m_f, divxfac::Float64,
 end
 
 
-"""
-    calculate_rlar(psi, n, l, q, epsr, wdian, wdiat, welec, wdhat, wbhat,
-                   nueff, dVdpsi, n_s, T_s, dbob_m_f, bo, bmin)::ComplexF64
-
-RLAR (Reduced Large Aspect Ratio) torque: energy-space integration with a pitch-averaged
-bounce response. Valid for large aspect ratio (ε ≪ 1); the κ bounce average is a placeholder estimate.
-
-Reference: [Logan et al., Phys. Plasmas, 2013]
-"""
-function calculate_rlar(psi, n, l, q, epsr, wdian, wdiat, welec,
-    wdhat, wbhat, nueff, dVdpsi::Float64,
-    n_s::Float64, T_s::Float64,
-    dbob_m_f, bo, bmin=0.5)::ComplexF64
-
-    # Setup parameters for energy integration
-    lnq = Float64(l)  # Resonant mode for trapped particles
-
-    # Maximum pitch angle parameter (use safe estimate)
-    lmdamax = min(1.0 / (1 - epsr), bo / bmin)
-
-    # Energy-space quadrature ∫ K(x) dx of the resonance operator K
-    xint = integrate_energy(wdian, wdiat, welec, wdhat, wbhat, nueff,
-        l, lnq, n, psi, lmdamax, "rlar")
-
-    # Kappa/bounce averaging (effect of field perturbations)
-    # Placeholder: simplified estimate
-    kappaint_val = sqrt(mean(abs.(dbob_m_f) .^ 2))
-
-    # dV/dpsi gradient term from Clebsch coordinate Jacobian
-    psi_factor = dVdpsi
-
-    # Normalization: √(ε/(2π³)) * n² * n_s * T_s
-    norm = sqrt(epsr / (2.0 * π^3)) * Float64(n)^2 * n_s * T_s
-
-    # dT/dψ = -0.5 · dV/dψ · κ_int · x_int · norm
-    result = psi_factor * kappaint_val * 0.5 * (-xint) * norm
-
-    return result
-end
-
-
-"""
-    calculate_clar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo,
-                   bmax, bmin, n_s, T_s, mass, chrg, tspl, dbob_m_f, divx_m_f,
-                   divxfac, wdfac)::ComplexF64
-
-CLAR (Circular Large Aspect Ratio) torque. Not implemented: warns once and returns zero.
-"""
-function calculate_clar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo,
-    bmax, bmin, n_s::Float64, T_s::Float64, mass, chrg, tspl,
-    dbob_m_f, divx_m_f, divxfac, wdfac)::ComplexF64
-
-    @warn "CLAR method not yet fully implemented, returning zero" maxlog = 1
-    return ComplexF64(0.0, 0.0)
-end
-
-
 """
     calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo, bmax,
                   bmin, n_s, T_s, mass, chrg, tspl, dbob_m_f, divx_m_f,
@@ -521,8 +451,6 @@ function _setup_surface_state(
 
     wdian = -twopi * T_s * dn_s_dpsi / (chrg * intr.chi1 * n_s)
     wdiat = -twopi * dT_s_dpsi / (chrg * intr.chi1)
-    wtran = sqrt(2 * T_s / mass) / (q * intr.ro)
-    wgyro = chrg * intr.bo / mass
     nuk = nu_s
 
     rsquared_bmin = equil.rzphi_rsquared((psi, theta_bmin))
@@ -532,7 +460,6 @@ function _setup_surface_state(
     avg_r = equil.geometry.avg_r_spline(psi)
     avg_R = equil.geometry.avg_R_spline(psi)
     epsr = avg_r / avg_R
-    dVdpsi = equil.profiles.dVdpsi_spline(psi)
 
     # Cells and stationary points of B(θ), reused by every λ of every harmonic.
     bf = _surface_b_field(B_vpar)
@@ -541,8 +468,8 @@ function _setup_surface_state(
         chrg, mass,
         tspl, B_vpar, bf, bmax, bmin, theta_bmax,
         q, n_s, T_s, welec,
-        wdian, wdiat, wtran, wgyro, nuk,
-        epsr, dVdpsi
+        wdian, wdiat, nuk,
+        epsr
     )
 end
 
diff --git a/src/KineticForces/Utils.jl b/src/KineticForces/Utils.jl
index 0f17f7925d..60c3aa03e8 100644
--- a/src/KineticForces/Utils.jl
+++ b/src/KineticForces/Utils.jl
@@ -56,7 +56,7 @@ density peaks (Logan & Park, Phys. Plasmas 20, 122507 (2013), §IV–V). `xeval`
 peak of the Maxwellian-weighted drive `x^2.5·e^−x`. The ℓ = 0 node is the ω_d-shifted ExB
 (superbanana-plateau) resonance.
 
-The frequencies are the cheap large-aspect-ratio closed forms of the `rlar` method,
+The frequencies are the cheap large-aspect-ratio closed forms
 `ω_b = (π/4)·√(ε/2)·wtran` and `ω_d = q·T_s/(2·ε·R₀²·Z·e·B₀)·wdfac`, with ε = ⟨r⟩/⟨R⟩ clamped
 away from the axis; no bounce averaging, which is enough for panel placement.
 """
diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl
index 7ff86b067a..8378c52cd5 100644
--- a/test/runtests_equil.jl
+++ b/test/runtests_equil.jl
@@ -116,7 +116,7 @@ using TOML
 
     @testset "newq0 q-profile revision" begin
         # newq0 is a target q(0), not an index: non-integer requests must survive the config,
-        # and neither the -1 sentinel nor an explicit target may write back to the frozen config.
+        # and q is positive, so a negative target is rejected.
         newq0_config(newq0) = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(;
             eq_filename=joinpath(data_dir, "EQDSK_COCOS_02"),
             eq_type="efit",
@@ -128,15 +128,7 @@ using TOML
             newq0=newq0
         )
 
-        base_config = newq0_config(0)
-        eq_base = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(base_config)
-
-        # -1 means "flip the sign of the axis extrapolation": f0fac vanishes, so ffac is
-        # exactly -1 and the revised q-profile is the negated baseline.
-        flip_config = newq0_config(-1)
-        eq_flip = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(flip_config)
-        @test flip_config.newq0 == -1.0
-        @test eq_flip.profiles.q_spline.y == -eq_base.profiles.q_spline.y
+        @test_throws ErrorException newq0_config(-1)
 
         # A non-integer target used to throw InexactError at construction (newq0 was ::Int)
         target_config = newq0_config(1.05)
@@ -148,6 +140,17 @@ using TOML
         xs = eq_target.profiles.xs
         q_axis = eq_target.profiles.q_spline.y[1] - eq_target.profiles.q_deriv(xs[1]; hint=Ref(1)) * xs[1]
         @test isapprox(q_axis, 1.05; rtol=0.02)
+
+        # The inverse path rescales F to meet the target too.
+        Eq = GeneralizedPerturbedEquilibrium.Equilibrium
+        eq_inv = Eq.PlasmaEquilibrium(
+            Eq.EquilibriumConfig(; eq_type="chease_binary", eq_filename=joinpath(data_dir, "INP1_binary"),
+                jac_type="boozer", mpsi=16, psilow=0.01, psihigh=0.99, r0_norm=6.8, b0_norm=7.4, newq0=1.5)
+        )
+        xs = eq_inv.profiles.xs
+        @test isapprox(eq_inv.profiles.q_spline.y[1] - eq_inv.profiles.q_deriv(xs[1]; hint=Ref(1)) * xs[1], 1.5; rtol=0.02)
+        # CHEASE supplies q directly, so a non-positive q stops the read.
+        @test_throws r"GPEC requires q > 0" Eq._chease_q([1.2, -0.1], (; eq_filename="INP1"))
     end
 
     @testset "flux_surface_rtol must be positive and finite" begin
diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl
index 1b03704e1f..2be7e099d7 100644
--- a/test/runtests_kinetic.jl
+++ b/test/runtests_kinetic.jl
@@ -273,7 +273,7 @@ using TOML
                 for btarget in range(blo, bhi; length=400)[2:end-1]
                     lmda = bo / btarget
                     roots = copy(KF._bounce_points_at_lambda!(buf, hints, bf2, lmda, bo))
-                    ref = sort!(KF.Roots.find_zeros(θ -> KF._vpar_from_spline(Bv2, lmda, bo, θ), 0.0, 1.0); rev=true)
+                    ref = sort!(KF.Roots.find_zeros(θ -> KF._vpar_from_spline(Bv2, lmda, bo, θ), 0.0, 1.0))
                     @test length(roots) == length(ref)
                     # The crossing inside the shared cell lies strictly between the two
                     # stationary points that bound it.
@@ -295,13 +295,13 @@ using TOML
             hints = ones(Int, length(bf.theta) + 1)
             for lmda in range(lmdatpb * (1 + 1e-6), lmdamax * (1 - 1e-6); length=200)
                 new_roots = copy(KF._bounce_points_at_lambda!(buf, hints, bf, lmda, bo))
-                ref_roots = sort!(KF.Roots.find_zeros(θ -> KF._vpar_from_spline(B_vpar, lmda, bo, θ), 0.0, 1.0); rev=true)
+                ref_roots = sort!(KF.Roots.find_zeros(θ -> KF._vpar_from_spline(B_vpar, lmda, bo, θ), 0.0, 1.0))
                 @test length(new_roots) == length(ref_roots)
                 if length(new_roots) == length(ref_roots)
                     @test all(abs.(new_roots .- ref_roots) .< 1e-9)
                 end
-                # Roots are returned in the descending order downstream assumes.
-                @test issorted(new_roots; rev=true)
+                # Roots are returned ascending, as Fortran spline_roots, so the well search pairs neighbours.
+                @test issorted(new_roots)
                 # v_par vanishes at each root.
                 for θ in new_roots
                     # Solved directly from the cell's cubic, so v_par vanishes to
@@ -311,6 +311,44 @@ using TOML
             end
         end
 
+        @testset "deepest well of a four-root surface lies between neighbouring roots" begin
+            # Two B minima give four bounce points at λ between the minima's B values; every well must be a
+            # neighbouring pair, so the chosen window holds no other root and v_par > 0 at its midpoint.
+            xs4 = collect(range(0.0, 1.0; length=257))
+            B4 = @. (1 + 0.1cos(4pi * xs4 + 1.1)) / (1 + 0.05cos(2pi * xs4))
+            B4[end] = B4[1]
+            Bv4 = cubic_interp(xs4, B4; bc=PeriodicBC())
+            bf4 = KF._surface_b_field(Bv4)
+            buf = Float64[]
+            hints = ones(Int, length(bf4.theta) + 1)
+            nfour = 0
+            for lmda in range(1 / maximum(B4), 1 / minimum(B4); length=200)[2:end-1]
+                roots = copy(KF._bounce_points_at_lambda!(buf, hints, bf4, lmda, 1.0))
+                length(roots) == 4 || continue
+                nfour += 1
+                t1, t2 = KF._find_deepest_well(roots, Bv4, lmda, 1.0)
+                @test !any(r -> t1 < r < t2 || t1 < r + 1 < t2, roots)
+                @test KF._vpar_from_spline(Bv4, lmda, 1.0, mod(0.5 * (t1 + t2), 1.0)) > 0
+            end
+            @test nfour > 0
+        end
+
+        @testset "past the midpoint the drift integrand holds its own last value" begin
+            # v_par turns negative at θ = 0.75: the bounce and drift integrands each hold their own previous sample.
+            ntheta = 32
+            pts = collect(range(0.0, 1.0; length=ntheta))
+            wts = fill(1 / (ntheta - 1), ntheta)
+            tspl = (out, θ) -> (out[1] = 1.0; out[2] = 0.5; out[4] = 1.0; out[5] = 0.3; out)
+            Bstep = θ -> θ < 0.75 ? 1.0 : 3.0
+            scr = KF.BounceScratch(ntheta, 1)
+            KF._bounce_samples!(zeros(ntheta), zeros(ComplexF64, ntheta), nothing, nothing, pts, wts, 0.5, 1, 1.0, 1.0, tspl, Bstep, 1.0,
+                [1], [0.0im], [0.0im], 1.0, false, 1, scr)
+            k = findfirst(>=(0.75), pts)
+            @test scr.g_wd[end] == scr.g_wd[k-1] != 0
+            @test scr.g_wb[end] == scr.g_wb[k-1]
+            @test scr.g_wd[end] != scr.g_wb[end]
+        end
+
         @testset "resolves the pair straddling an extremum" begin
             # Approaching a boundary, the two bounce points collapse onto the extremum
             # of B from either side. Assert that directly rather than against another
@@ -712,6 +750,9 @@ using TOML
         # Nonzero atol_psi dominating a small torque: the silent-garbage scenario
         ctrl_atol = KF.KineticForcesControl(; atol_psi=1e-2)
         @test_logs (:warn, r"atol_psi") KF.check_psi_quadrature_convergence(1e-3 + 0.0im, 1e-3, ctrl_atol, "fgar")
+        # A NaN torque fails every tolerance comparison, so it must error instead of passing silently
+        @test_throws ErrorException KF.check_psi_quadrature_convergence(complex(NaN), NaN, ctrl, "fgar")
+        @test_throws ErrorException KF.check_psi_quadrature_convergence(total, NaN, ctrl, "fgar")
     end
 
     @testset "METHOD_REGISTRY" begin
@@ -722,7 +763,7 @@ using TOML
         for entry in KF.METHOD_REGISTRY
             @test entry.flag in fieldnames(KF.KineticForcesControl)
             @test endswith(string(entry.flag), "_flag")
-            @test entry.kind in (:gar, :fcgl, :rlar, :clar)
+            @test entry.kind in (:gar, :fcgl)
             @test KF.method_kind(entry.name) == entry.kind
         end
     end
diff --git a/test/runtests_multiion.jl b/test/runtests_multiion.jl
index ed2e78d9ee..7a1e62c8c6 100644
--- a/test/runtests_multiion.jl
+++ b/test/runtests_multiion.jl
@@ -178,6 +178,42 @@ using HDF5
         rm(h5)
     end
 
+    @testset "cold edge keeps the collision model finite and monotone" begin
+        psi = collect(0.0:0.05:1.0)
+        cold(a) = @. a * (1 - psi)^2
+        data = EQ.KineticProfileData(; psi=psi, n_i=fill(1.0e19, length(psi)), n_e=fill(1.1e19, length(psi)),
+            T_i=cold(2000.0), T_e=cold(2500.0), omega_E=fill(1.0e4, length(psi)))
+        fine = 0.0:0.001:1.0
+        for kp in (EQ.kinetic_profile_splines(data; zi=1, zimp=6, mi=2, mimp=12),
+            only(s.profiles for s in EQ.resolve_ntv_species(data, [KF.IonSpecies(; z=1, m=2, fraction=1.0)]; zimp=6, mimp=12) if s.z == 1))
+            @test all(>=(1.0), kp.loglam_spline.(fine))
+            for nu in (kp.nui_spline, kp.nue_spline)
+                @test all(v -> isfinite(v) && v > 0, nu.(fine))
+                # The ν rise toward the cold edge stays monotone between grid nodes
+                @test issorted(nu.(0.5:0.001:1.0))
+            end
+        end
+        @test_throws ErrorException EQ.kinetic_profile_splines(EQ.KineticProfileData(data; T_e=cold(2500.0) .- 1.0))
+        @test_throws ErrorException EQ.resolve_ntv_species(EQ.KineticProfileData(data; T_i=cold(2000.0) .- 1.0), [KF.IonSpecies(; z=1, m=2, fraction=1.0)])
+        # The collision model divides by the densities, so a zero density stops the run as in Fortran PENTRC.
+        @test_throws r"densities must be positive" EQ.kinetic_profile_splines(EQ.KineticProfileData(data; n_e=[fill(1.1e19, length(psi) - 1); 0.0]))
+        # ...but only for the NTV view: the two-pass grid leaves such a file out instead of stopping a deck that does not use NTV.
+        @test EQ.ntv_profiles_usable(data)
+        @test !EQ.ntv_profiles_usable(EQ.KineticProfileData(data; n_e=[fill(1.1e19, length(psi) - 1); 0.0]))
+        # Scaling factors must keep the profiles physical.
+        @test_throws ArgumentError EQ.kinetic_profile_splines(data; temperature_factor=-1.0, chi1=1.0)
+        @test_throws ArgumentError EQ.kinetic_profile_splines(data; density_factor=0.0, chi1=1.0)
+
+        # A pedestal table: the cubic resample overshoots the top, the monotone one stays within the table's range.
+        ped_psi = [0.0, 0.3, 0.6, 0.85, 0.9, 0.95, 1.0]
+        ped_T = [2000.0, 1900.0, 1700.0, 1500.0, 400.0, 100.0, 50.0]
+        reg = collect(0:100) ./ 100
+        @test maximum(EQ.cubic_interp(ped_psi, ped_T; extrap=EQ.ExtendExtrap()).(reg)) > maximum(ped_T) ||
+              minimum(EQ.cubic_interp(ped_psi, ped_T; extrap=EQ.ExtendExtrap()).(reg)) < minimum(ped_T)
+        resampled = EQ._pchip_resample(ped_psi, ped_T, reg)
+        @test minimum(ped_T) <= minimum(resampled) && maximum(resampled) <= maximum(ped_T)
+    end
+
     @testset "combine_species_states" begin
         # Combined profile = Σ species dT/dψ interpolated onto the union grid; total = Σ total.
         mk(psi, dt) = (
diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl
index 78f35b5d6a..975a2b57ed 100644
--- a/test/runtests_solve_api.jl
+++ b/test/runtests_solve_api.jl
@@ -563,6 +563,13 @@ using HDF5
             @test_throws "script-only" parsed(Dict{String,Any}("kinetic" => "TestKineticMatrices"))
             @test_throws "not a kinetic source" parsed(Dict{String,Any}("kinetic" => "fixed"))
             @test_throws "no longer a key" parsed(Dict{String,Any}("kinetic_source" => "fixed"))
+            # The removed large-aspect-ratio torque methods stop the parse whatever their value.
+            for key in ("rlar_flag", "clar_flag"), value in (true, false)
+                mktempdir() do dir
+                    deck_dir = _stage_deck(dir; extra_sections=Dict{String,Any}("KineticForces" => Dict{String,Any}(key => value)))
+                    @test_throws "large-aspect-ratio methods are not ported" GPEC.parse_deck(TOML.parsefile(joinpath(deck_dir, "gpec.toml")), deck_dir)
+                end
+            end
         end
 
         @testset "kinetic torque" begin
diff --git a/test/test_data/regression_solovev_ideal_example/gpec.toml b/test/test_data/regression_solovev_ideal_example/gpec.toml
index 63eaaf40bc..b2663d180e 100644
--- a/test/test_data/regression_solovev_ideal_example/gpec.toml
+++ b/test/test_data/regression_solovev_ideal_example/gpec.toml
@@ -9,7 +9,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 16                               # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7                # Relative tolerance of the ODE solves that trace each flux surface
 
 [Wall]
diff --git a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml
index 7bab16e4ef..0571479748 100644
--- a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml
+++ b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml
@@ -9,7 +9,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 16                               # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7                # Relative tolerance of the ODE solves that trace each flux surface
 
 [Wall]
diff --git a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml
index 8769042f13..9285035b38 100644
--- a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml
+++ b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml
@@ -10,7 +10,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 16                               # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7                # Relative tolerance of the ODE solves that trace each flux surface
 
 [Wall]
diff --git a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml
index ac207f8e42..2e8954e36a 100644
--- a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml
+++ b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml
@@ -11,7 +11,7 @@ psilow = 1e-4                           # Lower limit of normalized flux coordin
 psihigh = 0.9995                        # Upper limit of normalized flux coordinate
 mpsi = 16                               # Number of radial grid points
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7                # Relative tolerance of the ODE solves that trace each flux surface
 
 [Wall]