From 56e7dcc68314c662068a57e1faa0d052e450f0f1 Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Thu, 17 Sep 2026 00:31:33 +0200
Subject: [PATCH 1/3] Add stability_derivatives and trim_angle on linearize

Derivatives of the force and moment coefficients with respect to angle of
attack and sideslip, by the chain rule through linearize's va columns, and a
trim-angle search on CMy with the slope from those derivatives. The inflow
formula moves out of set_va!(body_aero, settings) into apparent_wind so both
share it.

Refs #330

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 CHANGELOG.md                  |  6 +++
 docs/src/functions.md         |  2 +
 docs/src/private_functions.md |  3 ++
 src/VortexStepMethod.jl       |  2 +
 src/body_aerodynamics.jl      | 38 ++++++----------
 src/stability.jl              | 81 +++++++++++++++++++++++++++++++++++
 test/runtests.jl              |  1 +
 test/solver/test_stability.jl | 80 ++++++++++++++++++++++++++++++++++
 8 files changed, 188 insertions(+), 25 deletions(-)
 create mode 100644 src/stability.jl
 create mode 100644 test/solver/test_stability.jl

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4c00c754..fbc80e7d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,12 @@
 
 ## Unreleased
 
+### Added
+
+- `stability_derivatives` gives the force and moment coefficients and their derivatives
+  with respect to angle of attack and sideslip, and `trim_angle` the angles of attack at
+  which `CMy` changes sign, with the slope that says whether each trim is stable.
+
 ### Changed
 
 - Requires Julia 1.12 or 1.13; 1.10 and 1.11 keep resolving v5.1.1.
diff --git a/docs/src/functions.md b/docs/src/functions.md
index fb805d29..76d8dc8a 100644
--- a/docs/src/functions.md
+++ b/docs/src/functions.md
@@ -102,6 +102,8 @@ solve!
 solve_base!
 reinit!(body_aero::BodyAerodynamics{P, W, T}) where {P, W, T}
 linearize
+stability_derivatives
+trim_angle
 calculate_results
 ```
 
diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index d7807c42..807d89b1 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -25,6 +25,9 @@ calculate_cd
 calculate_cm
 calculate_cd_cm
 set_pitch_rate_dist!
+apparent_wind
+pitch_moment_coeff
+bisect_trim
 calculate_relative_alpha_and_velocity
 calculate_relative_alpha_and_relative_velocity
 update_effective_angle_of_attack!
diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl
index 0521ddfc..ab30ae28 100644
--- a/src/VortexStepMethod.jl
+++ b/src/VortexStepMethod.jl
@@ -32,6 +32,7 @@ export slice_args, preview_args
 export ObjWing, Section, Wing, refine!, reinit!
 export BodyAerodynamics
 export Solver, VSMSolution, linearize, solve, solve!, solve_base!, calc_forces!
+export stability_derivatives, trim_angle
 export SolveFailure
 export calculate_results
 export add_section!, set_va!, section_pitch_rate
@@ -432,6 +433,7 @@ include("panel.jl")
 include("body_aerodynamics.jl")
 include("wake.jl")
 include("solver.jl")
+include("stability.jl")
 
 include("plotting_helpers.jl")
 
diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index d04a966d..6cd45002 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -1147,24 +1147,19 @@ function set_va!(body_aero::BodyAerodynamics, va_distribution::AbstractMatrix;
 end
 
 """
-    set_va!(body_aero::BodyAerodynamics, settings::VSMSettings)
-
-Set velocity array from VSM settings configuration.
+    apparent_wind(alpha, beta, wind_speed)
 
-This convenience method extracts flight conditions from VSMSettings and 
-constructs the velocity vector in the body reference frame based on:
-- Wind speed from settings.condition.wind_speed
-- Angle of attack from settings.condition.alpha (converted from degrees)
-- Sideslip angle from settings.condition.beta (converted from degrees)
+Apparent wind vector in the body frame [m/s] at angle of attack `alpha` [rad], sideslip
+`beta` [rad] and `wind_speed` [m/s].
+"""
+apparent_wind(alpha, beta, wind_speed) =
+    wind_speed .* [cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta)]
 
-The velocity vector is constructed as:
-- X_b (forward): wind_speed * cos(α) * cos(β)  
-- Y_b (right): wind_speed * sin(β)
-- Z_b (down): wind_speed * sin(α) * cos(β)
+"""
+    set_va!(body_aero::BodyAerodynamics, settings::VSMSettings)
 
-# Arguments
-- `body_aero::BodyAerodynamics`: The aerodynamic body to modify
-- `settings::VSMSettings`: Settings object containing flight conditions
+Set the uniform inflow of `body_aero` to the [`apparent_wind`](@ref) at the `alpha` and
+`beta` [°] and `wind_speed` [m/s] of `settings.condition`.
 
 # Example
 ```julia
@@ -1174,15 +1169,8 @@ set_va!(body_aero, settings)
 ```
 """
 function set_va!(body_aero::BodyAerodynamics, settings::VSMSettings)
-    α = deg2rad(settings.condition.alpha)
-    β = deg2rad(settings.condition.beta)
-    wind_speed = settings.condition.wind_speed
-    
-    va = wind_speed * [
-        cos(α)*cos(β),  # X_b (forward)
-        sin(β),         # Y_b (right)
-        sin(α)*cos(β)   # Z_b (down)
-    ]
-    
+    condition = settings.condition
+    va = apparent_wind(deg2rad(condition.alpha), deg2rad(condition.beta),
+        condition.wind_speed)
     set_va!(body_aero, va)
 end
diff --git a/src/stability.jl b/src/stability.jl
new file mode 100644
index 00000000..e384d9f4
--- /dev/null
+++ b/src/stability.jl
@@ -0,0 +1,81 @@
+"""
+    stability_derivatives(solver, body_aero, alpha, beta, wind_speed; kwargs...)
+
+Aerodynamic coefficients `[CFx, CFy, CFz, CMx, CMy, CMz]` of `body_aero` at angle of attack
+`alpha` [rad], sideslip `beta` [rad] and `wind_speed` [m/s], and their derivatives with
+respect to `alpha` and `beta` [1/rad], at the rotation rate `body_aero.omega` and with
+moments about `solver.reference_point`. `kwargs` go to [`linearize`](@ref), which leaves
+`body_aero` at this inflow.
+
+Returns `(coeffs, dalpha, dbeta, converged)`.
+"""
+function stability_derivatives(solver::Solver, body_aero::BodyAerodynamics, alpha, beta,
+        wind_speed; kwargs...)
+    va = apparent_wind(alpha, beta, wind_speed)
+    jac, results, converged = linearize(solver, body_aero, va;
+        theta_idxs=nothing, va_idxs=1:3, aero_coeffs=true, kwargs...)
+    dva_dalpha = ForwardDiff.derivative(a -> apparent_wind(a, beta, wind_speed), alpha)
+    dva_dbeta = ForwardDiff.derivative(b -> apparent_wind(alpha, b, wind_speed), beta)
+    coeff_jac = jac[1:6, :]
+    return (coeffs=results[1:6], dalpha=coeff_jac * dva_dalpha,
+        dbeta=coeff_jac * dva_dbeta, converged)
+end
+
+"""
+    trim_angle(solver, body_aero, beta, wind_speed; alpha_range=deg2rad.(-5:2:15),
+               alpha_tol=1e-5, kwargs...)
+
+Angles of attack [rad] at which `CMy` of `body_aero` about `solver.reference_point` changes
+sign between neighbouring entries of `alpha_range`, bisected to `alpha_tol` [rad], at
+sideslip `beta` [rad] and `wind_speed` [m/s]. Returns one `(alpha, dCMy_dalpha)` per trim,
+the slope [1/rad] from [`stability_derivatives`](@ref) with `kwargs`; a trim is statically
+stable where `dCMy_dalpha < 0`.
+"""
+function trim_angle(solver::Solver, body_aero::BodyAerodynamics, beta, wind_speed;
+        alpha_range=deg2rad.(-5:2:15), alpha_tol=1e-5, kwargs...)
+    trims = @NamedTuple{alpha::Float64, dCMy_dalpha::Float64}[]
+    nose_down = [pitch_moment_coeff(solver, body_aero, alpha, beta, wind_speed) < 0
+                 for alpha in alpha_range]
+    for i in 1:length(alpha_range)-1
+        nose_down[i] == nose_down[i+1] && continue
+        alpha = bisect_trim(solver, body_aero, alpha_range[i], alpha_range[i+1],
+            nose_down[i], beta, wind_speed, alpha_tol)
+        derivatives = stability_derivatives(solver, body_aero, alpha, beta, wind_speed;
+            kwargs...)
+        push!(trims, (alpha=alpha, dCMy_dalpha=derivatives.dalpha[5]))
+    end
+    return trims
+end
+
+"""
+    pitch_moment_coeff(solver, body_aero, alpha, beta, wind_speed)
+
+`CMy` of `body_aero` solved at angle of attack `alpha` [rad], sideslip `beta` [rad] and
+`wind_speed` [m/s], at the rotation rate `body_aero.omega`.
+"""
+function pitch_moment_coeff(solver, body_aero, alpha, beta, wind_speed)
+    set_va!(body_aero, apparent_wind(alpha, beta, wind_speed), body_aero.omega)
+    return solve!(solver, body_aero).moment_coeffs[2]
+end
+
+"""
+    bisect_trim(solver, body_aero, alpha_low, alpha_high, nose_down_low, beta, wind_speed,
+                alpha_tol)
+
+Bisect `[alpha_low, alpha_high]` [rad], across which the sign of `CMy` changes and
+`nose_down_low` is whether it is negative at `alpha_low`, to a width of `alpha_tol` [rad],
+and return its midpoint.
+"""
+function bisect_trim(solver, body_aero, alpha_low, alpha_high, nose_down_low, beta,
+        wind_speed, alpha_tol)
+    while alpha_high - alpha_low > alpha_tol
+        alpha = (alpha_low + alpha_high) / 2
+        nose_down = pitch_moment_coeff(solver, body_aero, alpha, beta, wind_speed) < 0
+        if nose_down == nose_down_low
+            alpha_low = alpha
+        else
+            alpha_high = alpha
+        end
+    end
+    return (alpha_low + alpha_high) / 2
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index aafe3fd8..14bf62da 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -57,6 +57,7 @@ function include_selected_tests()
     should_run_test("solver/test_forwarddiff.jl") && include("solver/test_forwarddiff.jl")
     should_run_test("solver/test_backend_comparison.jl") && include("solver/test_backend_comparison.jl")
     should_run_test("solver/test_unrefined_dist.jl") && include("solver/test_unrefined_dist.jl")
+    should_run_test("solver/test_stability.jl") && include("solver/test_stability.jl")
     should_run_test("verification/test_verification.jl") && include("verification/test_verification.jl")
     should_run_test("VortexStepMethod/test_VortexStepMethod.jl") && include("VortexStepMethod/test_VortexStepMethod.jl")
     should_run_test("wake/test_wake.jl") && include("wake/test_wake.jl")
diff --git a/test/solver/test_stability.jl b/test/solver/test_stability.jl
new file mode 100644
index 00000000..50a59590
--- /dev/null
+++ b/test/solver/test_stability.jl
@@ -0,0 +1,80 @@
+using VortexStepMethod
+using VortexStepMethod: apparent_wind
+using Test
+
+"""
+    trimmable_wing_aero(cm)
+
+A rectangular wing with lift slope 2π, `cl = 0.1` at zero angle of attack and a constant
+section `cm`.
+"""
+function trimmable_wing_aero(cm)
+    chord, span = 1.0, 6.0
+    alpha_range = deg2rad.(-10.0:5.0:20.0)
+    polar = (alpha_range, 2π .* alpha_range .+ 0.1, fill(0.02, length(alpha_range)),
+             fill(cm, length(alpha_range)))
+    wing = Wing(10)
+    add_section!(wing, [0.0, span / 2, 0.0], [chord, span / 2, 0.0], POLAR_VECTORS, polar)
+    add_section!(wing, [0.0, -span / 2, 0.0], [chord, -span / 2, 0.0], POLAR_VECTORS,
+                 polar)
+    refine!(wing)
+    return BodyAerodynamics([wing])
+end
+
+function coeffs_at(solver, body_aero, alpha, beta, wind_speed)
+    set_va!(body_aero, apparent_wind(alpha, beta, wind_speed))
+    sol = solve!(solver, body_aero)
+    return [sol.force_coeffs; sol.moment_coeffs]
+end
+
+@testset "stability_derivatives match central differences of solve!" begin
+    body_aero = trimmable_wing_aero(0.05)
+    solver = Solver(body_aero; reference_point=[0.25, 0.5, 0.1], use_gamma_prev=false)
+    alpha, beta, wind_speed, step = deg2rad(4.0), deg2rad(3.0), 20.0, 1e-4
+
+    derivatives = stability_derivatives(solver, body_aero, alpha, beta, wind_speed)
+    @test derivatives.converged
+    @test derivatives.coeffs ≈ coeffs_at(solver, body_aero, alpha, beta, wind_speed)
+
+    central_difference(coeffs_plus, coeffs_minus) = (coeffs_plus - coeffs_minus) / 2step
+    dalpha = central_difference(
+        coeffs_at(solver, body_aero, alpha + step, beta, wind_speed),
+        coeffs_at(solver, body_aero, alpha - step, beta, wind_speed))
+    dbeta = central_difference(
+        coeffs_at(solver, body_aero, alpha, beta + step, wind_speed),
+        coeffs_at(solver, body_aero, alpha, beta - step, wind_speed))
+    @test !iszero(dbeta)
+    @test derivatives.dalpha ≈ dalpha rtol = 1e-4 atol = 1e-6
+    @test derivatives.dbeta ≈ dbeta rtol = 1e-4 atol = 1e-6
+end
+
+@testset "trim_angle finds where CMy changes sign" begin
+    beta, wind_speed = 0.0, 20.0
+
+    @testset "moments about the leading edge: stable trim" begin
+        body_aero = trimmable_wing_aero(0.05)
+        solver = Solver(body_aero)
+        trims = trim_angle(solver, body_aero, beta, wind_speed)
+        @test length(trims) == 1
+        trim = only(trims)
+        @test abs(coeffs_at(solver, body_aero, trim.alpha, beta, wind_speed)[5]) < 1e-5
+        @test trim.dCMy_dalpha < 0
+        @test trim.dCMy_dalpha ≈
+              stability_derivatives(solver, body_aero, trim.alpha, beta, wind_speed).dalpha[5]
+    end
+
+    @testset "moments about the trailing edge: unstable trim" begin
+        body_aero = trimmable_wing_aero(-0.05)
+        solver = Solver(body_aero; reference_point=[1.0, 0.0, 0.0])
+        trim = only(trim_angle(solver, body_aero, beta, wind_speed))
+        @test abs(coeffs_at(solver, body_aero, trim.alpha, beta, wind_speed)[5]) < 1e-5
+        @test trim.dCMy_dalpha > 0
+    end
+
+    @testset "no sign change in alpha_range: no trim" begin
+        body_aero = trimmable_wing_aero(0.05)
+        solver = Solver(body_aero)
+        @test isempty(trim_angle(solver, body_aero, beta, wind_speed;
+                                 alpha_range=deg2rad.(4:2:12)))
+    end
+end

From 722752992d633300f767966af3f595fdb37ad1a7 Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Thu, 17 Sep 2026 00:33:11 +0200
Subject: [PATCH 2/3] Bisect the trim bracket with a generic sign-change helper

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 docs/src/private_functions.md |  2 +-
 src/stability.jl              | 40 +++++++++++++++++------------------
 test/solver/test_stability.jl |  4 ++--
 3 files changed, 23 insertions(+), 23 deletions(-)

diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index 807d89b1..812fdf01 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -27,7 +27,7 @@ calculate_cd_cm
 set_pitch_rate_dist!
 apparent_wind
 pitch_moment_coeff
-bisect_trim
+bisect_sign_change
 calculate_relative_alpha_and_velocity
 calculate_relative_alpha_and_relative_velocity
 update_effective_angle_of_attack!
diff --git a/src/stability.jl b/src/stability.jl
index e384d9f4..8046e935 100644
--- a/src/stability.jl
+++ b/src/stability.jl
@@ -14,8 +14,10 @@ function stability_derivatives(solver::Solver, body_aero::BodyAerodynamics, alph
     va = apparent_wind(alpha, beta, wind_speed)
     jac, results, converged = linearize(solver, body_aero, va;
         theta_idxs=nothing, va_idxs=1:3, aero_coeffs=true, kwargs...)
-    dva_dalpha = ForwardDiff.derivative(a -> apparent_wind(a, beta, wind_speed), alpha)
-    dva_dbeta = ForwardDiff.derivative(b -> apparent_wind(alpha, b, wind_speed), beta)
+    dva_dalpha = ForwardDiff.derivative(
+        angle -> apparent_wind(angle, beta, wind_speed), alpha)
+    dva_dbeta = ForwardDiff.derivative(
+        angle -> apparent_wind(alpha, angle, wind_speed), beta)
     coeff_jac = jac[1:6, :]
     return (coeffs=results[1:6], dalpha=coeff_jac * dva_dalpha,
         dbeta=coeff_jac * dva_dbeta, converged)
@@ -33,13 +35,14 @@ stable where `dCMy_dalpha < 0`.
 """
 function trim_angle(solver::Solver, body_aero::BodyAerodynamics, beta, wind_speed;
         alpha_range=deg2rad.(-5:2:15), alpha_tol=1e-5, kwargs...)
+    is_nose_down =
+        alpha -> pitch_moment_coeff(solver, body_aero, alpha, beta, wind_speed) < 0
+    nose_down = is_nose_down.(alpha_range)
     trims = @NamedTuple{alpha::Float64, dCMy_dalpha::Float64}[]
-    nose_down = [pitch_moment_coeff(solver, body_aero, alpha, beta, wind_speed) < 0
-                 for alpha in alpha_range]
     for i in 1:length(alpha_range)-1
         nose_down[i] == nose_down[i+1] && continue
-        alpha = bisect_trim(solver, body_aero, alpha_range[i], alpha_range[i+1],
-            nose_down[i], beta, wind_speed, alpha_tol)
+        alpha = bisect_sign_change(is_nose_down, alpha_range[i], alpha_range[i+1],
+            alpha_tol)
         derivatives = stability_derivatives(solver, body_aero, alpha, beta, wind_speed;
             kwargs...)
         push!(trims, (alpha=alpha, dCMy_dalpha=derivatives.dalpha[5]))
@@ -59,23 +62,20 @@ function pitch_moment_coeff(solver, body_aero, alpha, beta, wind_speed)
 end
 
 """
-    bisect_trim(solver, body_aero, alpha_low, alpha_high, nose_down_low, beta, wind_speed,
-                alpha_tol)
+    bisect_sign_change(predicate, low, high, tol)
 
-Bisect `[alpha_low, alpha_high]` [rad], across which the sign of `CMy` changes and
-`nose_down_low` is whether it is negative at `alpha_low`, to a width of `alpha_tol` [rad],
-and return its midpoint.
+Bisect `[low, high]`, across which the boolean `predicate` flips, to a width of `tol` and
+return the midpoint.
 """
-function bisect_trim(solver, body_aero, alpha_low, alpha_high, nose_down_low, beta,
-        wind_speed, alpha_tol)
-    while alpha_high - alpha_low > alpha_tol
-        alpha = (alpha_low + alpha_high) / 2
-        nose_down = pitch_moment_coeff(solver, body_aero, alpha, beta, wind_speed) < 0
-        if nose_down == nose_down_low
-            alpha_low = alpha
+function bisect_sign_change(predicate, low, high, tol)
+    predicate_low = predicate(low)
+    while high - low > tol
+        middle = (low + high) / 2
+        if predicate(middle) == predicate_low
+            low = middle
         else
-            alpha_high = alpha
+            high = middle
         end
     end
-    return (alpha_low + alpha_high) / 2
+    return (low + high) / 2
 end
diff --git a/test/solver/test_stability.jl b/test/solver/test_stability.jl
index 50a59590..07634e72 100644
--- a/test/solver/test_stability.jl
+++ b/test/solver/test_stability.jl
@@ -59,8 +59,8 @@ end
         trim = only(trims)
         @test abs(coeffs_at(solver, body_aero, trim.alpha, beta, wind_speed)[5]) < 1e-5
         @test trim.dCMy_dalpha < 0
-        @test trim.dCMy_dalpha ≈
-              stability_derivatives(solver, body_aero, trim.alpha, beta, wind_speed).dalpha[5]
+        derivatives = stability_derivatives(solver, body_aero, trim.alpha, beta, wind_speed)
+        @test trim.dCMy_dalpha ≈ derivatives.dalpha[5]
     end
 
     @testset "moments about the trailing edge: unstable trim" begin

From b50a2b5a4c4ae19547333a1903bd0e32c8698efc Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Thu, 17 Sep 2026 10:41:59 +0200
Subject: [PATCH 3/3] trim_angle throws on an unconverged solve and takes
 backend, not kwargs

The CMy sweep and bisection solve with throw_on_fail, as does the slope's
linearize, so a trim is never bisected on unconverged moments. trim_angle
takes an explicit backend instead of splatting kwargs into linearize, which
let a reference_point reach the slope but not the sweep. coeffs_at_angles
replaces pitch_moment_coeff and the test's copy of it.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 docs/src/private_functions.md |  3 ++-
 src/stability.jl              | 39 ++++++++++++++++++++++------------
 test/solver/test_stability.jl | 40 +++++++++++++++++++++++------------
 3 files changed, 53 insertions(+), 29 deletions(-)

diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index b45f7b15..c508a6f4 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -26,7 +26,8 @@ calculate_cm
 calculate_cd_cm
 set_pitch_rate_dist!
 apparent_wind
-pitch_moment_coeff
+coeffs_at_angles
+nose_down
 bisect_sign_change
 calculate_relative_alpha_and_velocity
 calculate_relative_alpha_and_relative_velocity
diff --git a/src/stability.jl b/src/stability.jl
index 28a3442e..9ca803cf 100644
--- a/src/stability.jl
+++ b/src/stability.jl
@@ -25,42 +25,53 @@ end
 
 """
     trim_angle(solver, body_aero, beta, wind_speed; alpha_range=deg2rad.(-5:2:15),
-               alpha_tol=1e-5, kwargs...)
+               alpha_tol=1e-5, backend=AutoForwardDiff())
 
 Angles of attack [rad] at which `CMy` of `body_aero` about `solver.reference_point` changes
 sign between neighbouring entries of `alpha_range`, bisected to `alpha_tol` [rad], at
 sideslip `beta` [rad] and `wind_speed` [m/s]. Returns one `(alpha, dCMy_dalpha)` per trim,
-the slope [1/rad] from [`stability_derivatives`](@ref) with `kwargs`; a trim is statically
-stable where `dCMy_dalpha < 0`.
+the slope [1/rad] from [`stability_derivatives`](@ref) with `backend`; a trim is statically
+stable where `dCMy_dalpha < 0`. Throws a [`SolveFailure`](@ref) if a solve misses the
+solver's tolerances.
 """
 function trim_angle(solver::Solver, body_aero::BodyAerodynamics, beta, wind_speed;
-        alpha_range=deg2rad.(-5:2:15), alpha_tol=1e-5, kwargs...)
-    is_nose_down =
-        alpha -> pitch_moment_coeff(solver, body_aero, alpha, beta, wind_speed) < 0
-    nose_down = is_nose_down.(alpha_range)
+        alpha_range=deg2rad.(-5:2:15), alpha_tol=1e-5, backend=AutoForwardDiff())
+    is_nose_down = alpha -> nose_down(solver, body_aero, alpha, beta, wind_speed)
+    nose_down_range = is_nose_down.(alpha_range)
     trims = @NamedTuple{alpha::Float64, dCMy_dalpha::Float64}[]
     for i in 1:length(alpha_range)-1
-        nose_down[i] == nose_down[i+1] && continue
+        nose_down_range[i] == nose_down_range[i+1] && continue
         alpha = bisect_sign_change(is_nose_down, alpha_range[i], alpha_range[i+1],
             alpha_tol)
         derivatives = stability_derivatives(solver, body_aero, alpha, beta, wind_speed;
-            kwargs...)
+            backend, throw_on_fail=true)
         push!(trims, (alpha=alpha, dCMy_dalpha=derivatives.dalpha[5]))
     end
     return trims
 end
 
 """
-    pitch_moment_coeff(solver, body_aero, alpha, beta, wind_speed)
+    coeffs_at_angles(solver, body_aero, alpha, beta, wind_speed)
 
-`CMy` of `body_aero` solved at angle of attack `alpha` [rad], sideslip `beta` [rad] and
-`wind_speed` [m/s], at the rotation rate `body_aero.omega`.
+Aerodynamic coefficients `[CFx, CFy, CFz, CMx, CMy, CMz]` of `body_aero` solved at angle of
+attack `alpha` [rad], sideslip `beta` [rad] and `wind_speed` [m/s], at the rotation rate
+`body_aero.omega`. Throws a [`SolveFailure`](@ref) if the solve misses the solver's
+tolerances.
 """
-function pitch_moment_coeff(solver, body_aero, alpha, beta, wind_speed)
+function coeffs_at_angles(solver, body_aero, alpha, beta, wind_speed)
     set_va!(body_aero, apparent_wind(alpha, beta, wind_speed), body_aero.omega)
-    return solve!(solver, body_aero).moment_coeffs[2]
+    sol = solve!(solver, body_aero; throw_on_fail=true)
+    return [sol.force_coeffs; sol.moment_coeffs]
 end
 
+"""
+    nose_down(solver, body_aero, alpha, beta, wind_speed)
+
+Whether `CMy` from [`coeffs_at_angles`](@ref) is negative.
+"""
+nose_down(solver, body_aero, alpha, beta, wind_speed) =
+    coeffs_at_angles(solver, body_aero, alpha, beta, wind_speed)[5] < 0
+
 """
     bisect_sign_change(predicate, low, high, tol)
 
diff --git a/test/solver/test_stability.jl b/test/solver/test_stability.jl
index 07634e72..3e26529f 100644
--- a/test/solver/test_stability.jl
+++ b/test/solver/test_stability.jl
@@ -1,5 +1,5 @@
 using VortexStepMethod
-using VortexStepMethod: apparent_wind
+using VortexStepMethod: coeffs_at_angles
 using Test
 
 """
@@ -21,12 +21,6 @@ function trimmable_wing_aero(cm)
     return BodyAerodynamics([wing])
 end
 
-function coeffs_at(solver, body_aero, alpha, beta, wind_speed)
-    set_va!(body_aero, apparent_wind(alpha, beta, wind_speed))
-    sol = solve!(solver, body_aero)
-    return [sol.force_coeffs; sol.moment_coeffs]
-end
-
 @testset "stability_derivatives match central differences of solve!" begin
     body_aero = trimmable_wing_aero(0.05)
     solver = Solver(body_aero; reference_point=[0.25, 0.5, 0.1], use_gamma_prev=false)
@@ -34,15 +28,16 @@ end
 
     derivatives = stability_derivatives(solver, body_aero, alpha, beta, wind_speed)
     @test derivatives.converged
-    @test derivatives.coeffs ≈ coeffs_at(solver, body_aero, alpha, beta, wind_speed)
+    @test derivatives.coeffs ≈
+          coeffs_at_angles(solver, body_aero, alpha, beta, wind_speed)
 
     central_difference(coeffs_plus, coeffs_minus) = (coeffs_plus - coeffs_minus) / 2step
     dalpha = central_difference(
-        coeffs_at(solver, body_aero, alpha + step, beta, wind_speed),
-        coeffs_at(solver, body_aero, alpha - step, beta, wind_speed))
+        coeffs_at_angles(solver, body_aero, alpha + step, beta, wind_speed),
+        coeffs_at_angles(solver, body_aero, alpha - step, beta, wind_speed))
     dbeta = central_difference(
-        coeffs_at(solver, body_aero, alpha, beta + step, wind_speed),
-        coeffs_at(solver, body_aero, alpha, beta - step, wind_speed))
+        coeffs_at_angles(solver, body_aero, alpha, beta + step, wind_speed),
+        coeffs_at_angles(solver, body_aero, alpha, beta - step, wind_speed))
     @test !iszero(dbeta)
     @test derivatives.dalpha ≈ dalpha rtol = 1e-4 atol = 1e-6
     @test derivatives.dbeta ≈ dbeta rtol = 1e-4 atol = 1e-6
@@ -57,7 +52,8 @@ end
         trims = trim_angle(solver, body_aero, beta, wind_speed)
         @test length(trims) == 1
         trim = only(trims)
-        @test abs(coeffs_at(solver, body_aero, trim.alpha, beta, wind_speed)[5]) < 1e-5
+        trim_coeffs = coeffs_at_angles(solver, body_aero, trim.alpha, beta, wind_speed)
+        @test abs(trim_coeffs[5]) < 1e-5
         @test trim.dCMy_dalpha < 0
         derivatives = stability_derivatives(solver, body_aero, trim.alpha, beta, wind_speed)
         @test trim.dCMy_dalpha ≈ derivatives.dalpha[5]
@@ -67,10 +63,26 @@ end
         body_aero = trimmable_wing_aero(-0.05)
         solver = Solver(body_aero; reference_point=[1.0, 0.0, 0.0])
         trim = only(trim_angle(solver, body_aero, beta, wind_speed))
-        @test abs(coeffs_at(solver, body_aero, trim.alpha, beta, wind_speed)[5]) < 1e-5
+        trim_coeffs = coeffs_at_angles(solver, body_aero, trim.alpha, beta, wind_speed)
+        @test abs(trim_coeffs[5]) < 1e-5
         @test trim.dCMy_dalpha > 0
     end
 
+    @testset "a NONLIN solver with backend=nothing finds the same trim" begin
+        body_aero = trimmable_wing_aero(0.05)
+        trim_loop = only(trim_angle(Solver(body_aero), body_aero, beta, wind_speed))
+        solver = Solver(body_aero; solver_type=NONLIN)
+        trim = only(trim_angle(solver, body_aero, beta, wind_speed; backend=nothing))
+        @test trim.alpha ≈ trim_loop.alpha atol = 1e-4
+        @test trim.dCMy_dalpha ≈ trim_loop.dCMy_dalpha rtol = 1e-4
+    end
+
+    @testset "a solve that misses the tolerances throws" begin
+        body_aero = trimmable_wing_aero(0.05)
+        solver = Solver(body_aero; max_iterations=2)
+        @test_throws SolveFailure trim_angle(solver, body_aero, beta, wind_speed)
+    end
+
     @testset "no sign change in alpha_range: no trim" begin
         body_aero = trimmable_wing_aero(0.05)
         solver = Solver(body_aero)