diff --git a/CHANGELOG.md b/CHANGELOG.md
index f65c25f9..31558fd1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,11 @@
 
 ### 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.
+- `apparent_wind(alpha, beta, wind_speed)` gives the body-frame inflow vector at an angle
+  of attack and sideslip, as `set_va!(body_aero, settings)` sets it.
 - `Solver(settings)` and `Solver(n_panels, n_unrefined_sections)` build a solver without
   a `BodyAerodynamics`; keyword arguments override the settings.
 - `set_va!(body_aero, va_vec, omega; reference_point)` turns the body about
diff --git a/README.md b/README.md
index 42bf46c5..265bbea6 100644
--- a/README.md
+++ b/README.md
@@ -131,7 +131,7 @@ refine!(wing)
 body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
-va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+va_vec = apparent_wind(alpha, 0.0, va)
 set_va!(body_aero, va_vec)
 ```
 It is possible to import the wing geometry using an `.obj` file as shown in the example `ram_air_kite.jl`. During the import the polars are calculated automatically using XFoil. This approach is valid for rigid wings and ram-air kites, but not for leading edge inflatable kites.
diff --git a/docs/src/examples.md b/docs/src/examples.md
index c214d2c7..14794e24 100644
--- a/docs/src/examples.md
+++ b/docs/src/examples.md
@@ -77,7 +77,7 @@ multiple wings.
 ###### Set inflow conditions
 
 ```julia
-julia> va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+julia> va_vec = apparent_wind(alpha, 0.0, va)
 julia> set_va!(body_aero, va_vec, [0, 0, 0.1])
 ```
 
diff --git a/docs/src/functions.md b/docs/src/functions.md
index fb805d29..8f934a66 100644
--- a/docs/src/functions.md
+++ b/docs/src/functions.md
@@ -96,12 +96,15 @@ CurrentModule = VortexStepMethod
 ## Setting the inflow conditions and solving
 ```@docs
 set_va!
+apparent_wind
 section_pitch_rate
 solve
 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/index.md b/docs/src/index.md
index fc1cfb2f..bf5a4d6d 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -137,7 +137,7 @@ refine!(wing)
 body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
-va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+va_vec = apparent_wind(alpha, 0.0, va)
 set_va!(body_aero, va_vec)
 ```
 
diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index 43016830..6b60ab7e 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -28,6 +28,9 @@ calculate_cd
 calculate_cm
 calculate_cd_cm
 set_pitch_rate_dist!
+coeffs_at_angles
+nose_down
+bisect_sign_change
 calculate_relative_alpha_and_velocity
 calculate_relative_alpha_and_relative_velocity
 update_effective_angle_of_attack!
diff --git a/examples/bench.jl b/examples/bench.jl
index dc09f9a4..08c4f116 100644
--- a/examples/bench.jl
+++ b/examples/bench.jl
@@ -37,7 +37,7 @@ refine!(wing)
 body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
-va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+va_vec = apparent_wind(alpha, 0.0, va)
 set_va!(body_aero, va_vec)
 
 # Step 4: Initialize solvers for both LLT and VSM methods
@@ -76,13 +76,8 @@ vsm_solver = Solver(
 va = 15.0
 aoa = 15.0
 side_slip = 0.0
-yaw_rate = 0.0
 aoa_rad = deg2rad(aoa)
-va_vec = [
-    cos(aoa_rad) * cos(side_slip),
-    sin(side_slip),
-    sin(aoa_rad)
-] * va
+va_vec = apparent_wind(aoa_rad, side_slip, va)
 set_va!(body_aero, va_vec)
 
 # Solving
diff --git a/examples/billowing.jl b/examples/billowing.jl
index e99d0d52..d78b8831 100644
--- a/examples/billowing.jl
+++ b/examples/billowing.jl
@@ -89,7 +89,7 @@ sideslip_deg = settings.condition.beta
 
 α0 = deg2rad(angle_of_attack_deg)
 β0 = deg2rad(sideslip_deg)
-va_vec = wind_speed .* [cos(α0) * cos(β0), sin(β0), sin(α0) * cos(β0)]
+va_vec = apparent_wind(α0, β0, wind_speed)
 set_va!(body_aero_flat, va_vec)
 set_va!(body_aero_bill, va_vec)
 
diff --git a/examples/linearize_check.jl b/examples/linearize_check.jl
index b8bd47fc..4117cf30 100644
--- a/examples/linearize_check.jl
+++ b/examples/linearize_check.jl
@@ -35,11 +35,7 @@ va         = 15.0
 aoa_deg    = 10.0
 aoa_rad    = deg2rad(aoa_deg)
 side_slip  = 0.0
-va_vec_b_0 = [
-    cos(aoa_rad) * cos(side_slip),
-    sin(side_slip),
-    sin(aoa_rad),
-] * va
+va_vec_b_0 = apparent_wind(aoa_rad, side_slip, va)
 omega_b_0  = zeros(3)
 theta_0    = zeros(n_unrefined)
 
diff --git a/examples/obj_to_yaml_kite.jl b/examples/obj_to_yaml_kite.jl
index 3a4d9354..3856533b 100644
--- a/examples/obj_to_yaml_kite.jl
+++ b/examples/obj_to_yaml_kite.jl
@@ -74,7 +74,7 @@ VortexStepMethod.reinit!(body_aero)
 
 solver = Solver(wing.n_panels, wing.n_unrefined_sections;
     aerodynamic_model_type=VSM, rtol=1e-5, solver_type=LOOP)
-set_va!(body_aero, [cos(deg2rad(8)) * va, 0.0, sin(deg2rad(8)) * va])
+set_va!(body_aero, apparent_wind(deg2rad(8), 0.0, va))
 results = VortexStepMethod.solve(solver, body_aero; log=true)
 
 if PLOT
diff --git a/examples/rectangular_wing.jl b/examples/rectangular_wing.jl
index 15119a9a..badb3a89 100644
--- a/examples/rectangular_wing.jl
+++ b/examples/rectangular_wing.jl
@@ -41,7 +41,7 @@ refine!(wing)
 body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
-va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+va_vec = apparent_wind(alpha, 0.0, va)
 set_va!(body_aero, va_vec, [0, 0, 0.1])
 
 # Step 4: Initialize solvers for both LLT and VSM methods
diff --git a/examples/stall_model.jl b/examples/stall_model.jl
index 168fd3bb..6b6acb3a 100644
--- a/examples/stall_model.jl
+++ b/examples/stall_model.jl
@@ -61,11 +61,7 @@ aoa = 17.0
 side_slip = 0.0
 yaw_rate = 0.0
 aoa_rad = deg2rad(aoa)
-va_vec = [
-    cos(aoa_rad) * cos(side_slip),
-    sin(side_slip),
-    sin(aoa_rad)
-] * va
+va_vec = apparent_wind(aoa_rad, side_slip, va)
 set_va!(body_aero, va_vec)
 
 # Plotting geometry
diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl
index f28a73f5..4c89927b 100644
--- a/src/VortexStepMethod.jl
+++ b/src/VortexStepMethod.jl
@@ -32,9 +32,10 @@ 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
+export add_section!, apparent_wind, set_va!, section_pitch_rate
 export calculate_projected_area, calculate_span
 export MVec3
 
@@ -437,6 +438,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 be9f9156..22e97679 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -1158,26 +1158,20 @@ function set_va!(body_aero::BodyAerodynamics, va_vec_dist::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)
-- Yaw rate from settings.condition.yaw_rate (converted from °/s), applied as `omega`
-  about Z_b and turning the body about `body_aero.reference_point`
+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`, turning the body about
+`body_aero.reference_point` at its `yaw_rate` [°/s] about Z_b.
 
 # Example
 ```julia
@@ -1187,16 +1181,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_vec = wind_speed * [
-        cos(α)*cos(β),  # X_b (forward)
-        sin(β),         # Y_b (right)
-        sin(α)*cos(β)   # Z_b (down)
-    ]
-    omega = [0.0, 0.0, deg2rad(settings.condition.yaw_rate)]
-
-    set_va!(body_aero, va_vec, omega)
+    condition = settings.condition
+    va_vec = apparent_wind(deg2rad(condition.alpha), deg2rad(condition.beta),
+        condition.wind_speed)
+    set_va!(body_aero, va_vec, [0.0, 0.0, deg2rad(condition.yaw_rate)])
 end
diff --git a/src/stability.jl b/src/stability.jl
new file mode 100644
index 00000000..9ca803cf
--- /dev/null
+++ b/src/stability.jl
@@ -0,0 +1,92 @@
+"""
+    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_vec = apparent_wind(alpha, beta, wind_speed)
+    jac, results, converged = linearize(solver, body_aero, va_vec;
+        theta_idxs=nothing, va_idxs=1:3, aero_coeffs=true, kwargs...)
+    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)
+end
+
+"""
+    trim_angle(solver, body_aero, beta, wind_speed; alpha_range=deg2rad.(-5:2:15),
+               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 `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, 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_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;
+            backend, throw_on_fail=true)
+        push!(trims, (alpha=alpha, dCMy_dalpha=derivatives.dalpha[5]))
+    end
+    return trims
+end
+
+"""
+    coeffs_at_angles(solver, body_aero, alpha, beta, wind_speed)
+
+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 coeffs_at_angles(solver, body_aero, alpha, beta, wind_speed)
+    set_va!(body_aero, apparent_wind(alpha, beta, wind_speed), body_aero.omega)
+    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)
+
+Bisect `[low, high]`, across which the boolean `predicate` flips, to a width of `tol` and
+return the midpoint.
+"""
+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
+            high = middle
+        end
+    end
+    return (low + high) / 2
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index b77976d4..3077e22b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -58,6 +58,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..7145322d
--- /dev/null
+++ b/test/solver/test_stability.jl
@@ -0,0 +1,97 @@
+using VortexStepMethod
+using VortexStepMethod: coeffs_at_angles
+using Test
+
+"""
+    trimmable_wing(cm; kwargs...)
+
+A rectangular wing with lift slope 2π, `cl = 0.1` at zero angle of attack and a constant
+section `cm`, as a `BodyAerodynamics` and a `Solver` sized for it. Keyword arguments go to
+the solver.
+"""
+function trimmable_wing(cm; kwargs...)
+    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]),
+           Solver(wing.n_panels, wing.n_unrefined_sections; kwargs...)
+end
+
+@testset "apparent_wind points along body x, tilts up with alpha, right with beta" begin
+    @test apparent_wind(0.0, 0.0, 10.0) == [10.0, 0.0, 0.0]
+    va_vec = apparent_wind(deg2rad(30), deg2rad(20), 10.0)
+    @test va_vec ≈ 10.0 .* [cosd(30) * cosd(20), sind(20), sind(30) * cosd(20)]
+    @test hypot(va_vec...) ≈ 10.0
+end
+
+@testset "stability_derivatives match central differences of solve!" begin
+    body_aero, solver = trimmable_wing(0.05; 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_angles(solver, body_aero, alpha, beta, wind_speed)
+
+    central_difference(coeffs_plus, coeffs_minus) = (coeffs_plus - coeffs_minus) / 2step
+    dalpha = central_difference(
+        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_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
+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, solver = trimmable_wing(0.05)
+        trims = trim_angle(solver, body_aero, beta, wind_speed)
+        @test length(trims) == 1
+        trim = only(trims)
+        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]
+    end
+
+    @testset "moments about the trailing edge: unstable trim" begin
+        body_aero, solver = trimmable_wing(-0.05; reference_point=[1.0, 0.0, 0.0])
+        trim = only(trim_angle(solver, body_aero, beta, wind_speed))
+        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
+        loop_body, loop_solver = trimmable_wing(0.05)
+        body_aero, solver = trimmable_wing(0.05; solver_type=NONLIN)
+        trim_loop = only(trim_angle(loop_solver, loop_body, beta, wind_speed))
+        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, solver = trimmable_wing(0.05; 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, solver = trimmable_wing(0.05)
+        @test isempty(trim_angle(solver, body_aero, beta, wind_speed;
+                                 alpha_range=deg2rad.(4:2:12)))
+    end
+end