From 39bfde64187824f16d17741cb34be7228609a108 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Fri, 21 Aug 2026 14:31:30 -0400 Subject: [PATCH] ForceFreeStates - BUGFIX! - Keep the sign of crit so fixed-boundary crossings are detected compute_smallest_eigenvalue returned findmin(abs, evals)[1], which is the smallest eigenvalue *magnitude* and therefore non-negative by construction. The caller detects conjugate points by testing crit_store[istep] * crit_store[istep-1] < 0, so that product could never be negative: no zero crossing was ever recorded and nzero was always zero, reporting every equilibrium as fixed-boundary stable. Return the signed eigenvalue at the smallest-magnitude index instead. This matches ode_output_get_crit in the Fortran, which sorts on -ABS(evalsi) to find the index and then returns evalsi at that index, keeping the sign. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_018YkwX4YMDVFaDi1NM9WcJD --- src/ForceFreeStates/FixedBoundaryStability.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/ForceFreeStates/FixedBoundaryStability.jl b/src/ForceFreeStates/FixedBoundaryStability.jl index 647cb45b2..c07ff5a65 100644 --- a/src/ForceFreeStates/FixedBoundaryStability.jl +++ b/src/ForceFreeStates/FixedBoundaryStability.jl @@ -92,7 +92,7 @@ end compute_smallest_eigenvalue(u) -> crit, nonherm Form the inverse plasma response matrix W⁻¹ using the solution matrix `u` and -returns its minimum eigenvalue by magnitude. Performs the same function as +returns its signed eigenvalue of smallest magnitude. Performs the same function as `ode_output_get_crit` in the Fortran code, except we explicitly form W⁻¹ here from U₁ * U₂⁻¹ using Julia's right division operator `/` instead of adj(adj(U₂)⁻¹ * adj(U₁)) as done in Fortran. We have also added a check to @@ -139,7 +139,8 @@ construction but may accumulate numerical noise during integration. # Enforce that W is Hermitian hermitianpart!(wp_inverse) # Overwrites W⁻¹ with (W⁻¹ + (W⁻¹)') / 2 - # Compute eigenvalues and return the smallest - crit = findmin(abs, eigvals!(Hermitian(wp_inverse)))[1] + # Return the eigenvalue of smallest magnitude, keeping its sign so the caller can detect crossings + evals = eigvals!(Hermitian(wp_inverse)) + crit = evals[findmin(abs, evals)[2]] return crit, nonherm end