diff --git a/CHANGELOG.md b/CHANGELOG.md
index 14d79173..5b21776e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,8 +5,10 @@
 ### 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.
+  with respect to angle of attack, sideslip and the nondimensional roll, pitch and yaw
+  rates p̂ = pb/2V, q̂ = q c_ref/2V, r̂ = rb/2V, turning about `solver.reference_point`,
+  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, va)` 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
diff --git a/src/stability.jl b/src/stability.jl
index 2ba74529..7c0d8658 100644
--- a/src/stability.jl
+++ b/src/stability.jl
@@ -2,25 +2,33 @@
     stability_derivatives(solver, body_aero, alpha, beta, va; kwargs...)
 
 Aerodynamic coefficients `[CFx, CFy, CFz, CMx, CMy, CMz]` of `body_aero` at angle of attack
-`alpha` [rad], sideslip `beta` [rad] and apparent wind speed `va` [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.
+`alpha` [rad], sideslip `beta` [rad], apparent wind speed `va` [m/s] and rotation rate
+`body_aero.omega` [rad/s], and their derivatives: `dalpha` and `dbeta` [1/rad], and `dp`,
+`dq`, `dr` with respect to p̂ = pb/2V, q̂ = q c_ref/2V and r̂ = rb/2V, where p, q, r are the
+rates about body x, y, z, b the wing span, c_ref `body_aero.c_ref` and V is `va`, not the
+area-weighted panel inflow speed the coefficients are normalised by. Moments are about,
+and the body turns about, `solver.reference_point`. `kwargs` go to [`linearize`](@ref).
+Leaves `body_aero` at this inflow with `body_aero.reference_point = solver.reference_point`.
 
-Returns `(coeffs, dalpha, dbeta, converged)`.
+Returns `(coeffs, dalpha, dbeta, dp, dq, dr, converged)`.
 """
 function stability_derivatives(solver::Solver, body_aero::BodyAerodynamics, alpha, beta,
         va; kwargs...)
     va_vec = apparent_wind(alpha, beta, va)
-    jac, results, converged = linearize(solver, body_aero, va_vec;
-        theta_idxs=nothing, va_vec_idxs=1:3, aero_coeffs=true, kwargs...)
+    set_va!(body_aero, va_vec, body_aero.omega; reference_point=solver.reference_point)
+    jac, results, converged = linearize(solver, body_aero, [va_vec; body_aero.omega];
+        theta_idxs=nothing, va_vec_idxs=1:3, omega_idxs=4:6, aero_coeffs=true, kwargs...)
     dva_dalpha = ForwardDiff.derivative(
         angle -> apparent_wind(angle, beta, va), alpha)
     dva_dbeta = ForwardDiff.derivative(
         angle -> apparent_wind(alpha, angle, va), beta)
     coeff_jac = jac[1:6, :]
-    return (coeffs=results[1:6], dalpha=coeff_jac * dva_dalpha,
-        dbeta=coeff_jac * dva_dbeta, converged)
+    span = body_aero.wings[1].span
+    return (coeffs=results[1:6], dalpha=coeff_jac[:, 1:3] * dva_dalpha,
+        dbeta=coeff_jac[:, 1:3] * dva_dbeta,
+        dp=coeff_jac[:, 4] * 2va / span,
+        dq=coeff_jac[:, 5] * 2va / body_aero.c_ref,
+        dr=coeff_jac[:, 6] * 2va / span, converged)
 end
 
 """
@@ -32,7 +40,8 @@ sign between neighbouring entries of `alpha_range`, bisected to `alpha_tol` [rad
 sideslip `beta` [rad] and apparent wind speed `va` [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.
+[`SolveFailure`](@ref) if a solve misses the solver's tolerances. Leaves
+`body_aero.reference_point = solver.reference_point`.
 """
 function trim_angle(solver::Solver, body_aero::BodyAerodynamics, beta, va;
         alpha_range=deg2rad.(-5:2:15), alpha_tol=1e-5, backend=AutoForwardDiff())
@@ -55,11 +64,13 @@ end
 
 Aerodynamic coefficients `[CFx, CFy, CFz, CMx, CMy, CMz]` of `body_aero` solved at angle of
 attack `alpha` [rad], sideslip `beta` [rad] and apparent wind speed `va` [m/s], at the
-rotation rate `body_aero.omega`. Throws a [`SolveFailure`](@ref) if the solve misses the
+rotation rate `body_aero.omega`, which it turns about and stores as
+`body_aero.reference_point`. Throws a [`SolveFailure`](@ref) if the solve misses the
 solver's tolerances.
 """
 function coeffs_at_angles(solver, body_aero, alpha, beta, va)
-    set_va!(body_aero, apparent_wind(alpha, beta, va), body_aero.omega)
+    set_va!(body_aero, apparent_wind(alpha, beta, va), body_aero.omega;
+        reference_point=solver.reference_point)
     sol = solve!(solver, body_aero; throw_on_fail=true)
     return [sol.force_coeffs; sol.moment_coeffs]
 end
diff --git a/test/solver/test_stability.jl b/test/solver/test_stability.jl
index 55f82ec7..758b68a1 100644
--- a/test/solver/test_stability.jl
+++ b/test/solver/test_stability.jl
@@ -52,6 +52,30 @@ end
     @test derivatives.dbeta ≈ dbeta rtol = 1e-4 atol = 1e-6
 end
 
+@testset "rate derivatives match central differences of solve! about reference_point" begin
+    reference_point = [0.25, 0.5, 0.1]
+    body_aero, solver = trimmable_wing(0.05; reference_point, use_gamma_prev=false)
+    alpha, beta, va, step = deg2rad(4.0), deg2rad(3.0), 20.0, 1e-4
+    omega = [0.1, -0.05, 0.08]
+    body_aero.omega = omega
+
+    derivatives = stability_derivatives(solver, body_aero, alpha, beta, va)
+    @test derivatives.converged
+    @test body_aero.reference_point == reference_point
+
+    rate_scales = 2va ./ [body_aero.wings[1].span, body_aero.c_ref, body_aero.wings[1].span]
+    for (axis, derivative) in enumerate((derivatives.dp, derivatives.dq, derivatives.dr))
+        rate_step = step .* (1:3 .== axis)
+        body_aero.omega = omega + rate_step
+        coeffs_plus = coeffs_at_angles(solver, body_aero, alpha, beta, va)
+        body_aero.omega = omega - rate_step
+        coeffs_minus = coeffs_at_angles(solver, body_aero, alpha, beta, va)
+        central_difference = (coeffs_plus - coeffs_minus) / 2step * rate_scales[axis]
+        @test !iszero(central_difference)
+        @test derivative ≈ central_difference rtol = 1e-4 atol = 1e-6
+    end
+end
+
 @testset "trim_angle finds where CMy changes sign" begin
     beta, va = 0.0, 20.0
 
@@ -75,6 +99,18 @@ end
         @test trim.dCMy_dalpha > 0
     end
 
+    @testset "pitching body: trim turning about the reference point" begin
+        reference_point = [1.0, 0.0, 0.0]
+        body_aero, solver = trimmable_wing(-0.05; reference_point)
+        omega = [0.0, 0.5, 0.0]
+        set_va!(body_aero, apparent_wind(0.0, beta, va), omega)
+        trim = only(trim_angle(solver, body_aero, beta, va))
+        set_va!(body_aero, apparent_wind(trim.alpha, beta, va), omega;
+                reference_point)
+        cmy = solve!(solver, body_aero).moment_coeffs[2]
+        @test abs(cmy) < 1e-5 * abs(trim.dCMy_dalpha)
+    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)