stability_derivatives gives roll, pitch and yaw rate derivatives about solver.reference_point

[1-Bort-1]

TL;DR

stability_derivatives now also returns dp, dq and dr: the derivatives of CFx..CMz with respect to p̂ = pb/2V, q̂ = q c_ref/2V and r̂ = rb/2V, taken from linearize's omega columns with the body turning about solver.reference_point. Rate derivatives are only meaningful when the body turns about the same point the moments are taken about, which #353 makes possible.

What changed

• stability_derivatives stores solver.reference_point as the pivot on body_aero (set_va!(...; reference_point) from set_va! turns the body about a stored reference_point, and set_va!(body_aero, settings) applies yaw_rate #353), linearises in [va; omega] at body_aero.omega, and scales the omega columns by 2V/b, 2V/c_ref and 2V/b, with b the wing span and V wind_speed.
• coeffs_at_angles, which trim_angle brackets with, turns about solver.reference_point too. Before, a pitching body was trimmed about whatever pivot body_aero held: with omega = [0, 0.5, 0] and reference_point = [1, 0, 0], the returned trim at α = −1.12° had CMy 3.9e-6 turning about the origin and −0.083 turning about the reference point. It is now 1.3e-5, inside the bisection tolerance.

Choice of chord for q̂

The Python compute_rigid_body_stability_derivatives normalises q̂ with c = S/b while its moment coefficients use the max chord. Here q̂ uses body_aero.c_ref (max panel chord), the same length the moment coefficients are divided by, so one chord runs through every pitch term. For a tapered or curved kite dq therefore differs from Python's by the factor c_ref/(S/b).

Pivot matters

Same wing, α = 4°, β = 3°, reference point [0.25, 0.5, 0.1]: dCFz/dq (raw, per rad/s) is 0.1087 turning about the reference point against 0.1654 turning about the origin, and dCMx/dp is −0.510 against −0.456.

Stack

Base is #346. The branch also merges #353's head, because the pivot lives there; #353's commits show in this diff until it merges. Where the two meet, set_va!(body_aero, settings) takes #346's apparent_wind form with #353's yaw_rate omega, and #346's coeffs_at_angles (which now feeds both the trim bracketing and the tests) turns about solver.reference_point.

Depends-On: #353

Verification

• n/a to reproduce: new feature. The trim_angle pivot bug was reproduced first: CMy −0.083 at the returned trim
• test/solver/test_stability.jl red before, green after (juliaserver): new "rate derivatives match central differences of solve! about reference_point" errored on missing dp/dq/dr and failed on body_aero.reference_point == [0,0,0], now 6/6. Central differences are over omega = [0.1, −0.05, 0.08] ± 1e-4 at rtol 1e-4. New "pitching body: trim turning about the reference point" failed at 0.0835 < 1e-5·|slope|, now 8/8. The existing α/β testset passes 5/5
• On 1a38c0b (merged with Add stability_derivatives (angle of attack, sideslip) and trim_angle, built on linearize #346 b50a2b5 and set_va! turns the body about a stored reference_point, and set_va!(body_aero, settings) applies yaw_rate #353 6ecd77a), green in a fresh juliaserver session: test_stability.jl 22/22 (incl. Add stability_derivatives (angle of attack, sideslip) and trim_angle, built on linearize #346's NONLIN and SolveFailure trim testsets) · test_body_aerodynamics.jl 4888/4888 · test_forwarddiff.jl 10/10
• Docs build clean on 1a38c0b (makedocs without deploy)
• Local full suite on 1a38c0b: PASS (8 min, Julia 1.12.7, one matrix cell) · GitHub CI: Documentation PASS; the test matrix is skipped on this pull request because its base is Add stability_derivatives (angle of attack, sideslip) and trim_angle, built on linearize #346's branch, not main
• Risk: dp/dq/dr use wind_speed for V while the solver's q_ref uses a reference speed built from the panel inflow distribution, so at nonzero omega the p̂ scaling and the coefficients' dynamic pressure use slightly different speeds.

Scope

Against a plain merge of #346 and #353: +72 / −55 in src/stability.jl, test/solver/test_stability.jl, CHANGELOG.md, and src/body_aerodynamics.jl, where the last is only the set_va!(body_aero, settings) resolution. The rest of the diff against #346 is #353. Stack: 2/2 on #346.

Closes #345 · task VortexStepMethod.jl-345

[1-Bort-1]

Independent review (advisory)

Verdict: APPROVE WITH COMMENTS · 2 inline, 0 off the diff

Good

• The rate columns follow the plan in the card: linearize runs over [va; omega] with omega_idxs=4:6, and columns 4..6 are scaled by 2V/b, 2V/c_ref and 2V/b (src/stability.jl:18-30).
• The pivot reaches both AD paths: FiniteDiff calls set_va! on body_aero itself, which keeps its stored reference_point, and make_dual_shadow copies reference_point into the dual body (checked in linearize at src/solver.jl:1297-1320).
• The rate test checks against central differences of solve! about the same pivot, with a guard that no column is zero, so a missing or mis-scaled column fails it.
• The trim regression test solves again at the returned alpha about the reference point and checks CMy against the slope, which pins the pivot bug the card describes.
• set_va! now has one inflow formula, va - omega × (control_point - reference_point), and a reference_point setter; the old all(iszero, omega) branch and the per-wing index loop are gone.
• The card says why c_ref and not S/b is used for q̂, and names the difference from Python.
• The changed src and test files have no line over 92 characters (grep for .{93,}).

Not good

• src/stability.jl:17 — stability_derivatives permanently overwrites body_aero.reference_point with solver.reference_point, and only the test says so; the docstring says 'leaves body_aero at this inflow'. A caller who set their own pivot will see later set_va!/reinit! calls turn about a different point with no warning.
• test/solver/test_stability.jl:64 — coeffs_at_rate is a nested closure with logic that nearly repeats the file's top-level coeffs_at helper (set_va!, solve!, stack coefficients). Extending coeffs_at with omega/reference_point keywords gives one helper and follows §2 and §6.
• The card does not say what was searched for before adding the rate scaling and the new test helper (§2).
• The stability_derivatives docstring sentence 'dp, dq, dr with respect to the body rates about x, y and z as p̂ = …' is hard to parse; 'with respect to p̂ = pb/2V, q̂ = …, r̂ = …' would say it directly.
• V is wind_speed in the p̂/q̂/r̂ scaling, but the coefficients are normalised by the reference speed built from the panel inflow. The card flags this as a risk, but neither the docstring nor the CHANGELOG mentions it.
• pitch_moment_coeff has the same pivot-writing side effect as stability_derivatives, so everything trim_angle touches changes body_aero.reference_point too.
• The ForwardDiff testset builds a new BodyAerodynamics([wing]) on every iteration while body_aero already exists. This is cheap, but it shares wing with body_aero without saying so.
• In the ForwardDiff test, turns reuses the outer omega for its first entry and then names the loop variable omega_op, so two omegas are in scope.

_{claude, rubric CLEAN_CODE.md. A different lab from the implementer
on purpose: a reviewer sharing its blind spots would not flag its mistakes.}

[1-Bort-1]

MINOR: stability_derivatives permanently overwrites body_aero.reference_point with solver.reference_point, and only the test says so; the docstring says 'leaves body_aero at this inflow'. A caller who set their own pivot will see later set_va!/reinit! calls turn about a different point with no warning.

[1-Bort-1]

MINOR: coeffs_at_rate is a nested closure with logic that nearly repeats the file's top-level coeffs_at helper (set_va!, solve!, stack coefficients). Extending coeffs_at with omega/reference_point keywords gives one helper and follows §2 and §6.

[1-Bort-1]

Local full suite: PASS (7 min, Julia 1.12.7, one cell of the matrix)