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 made possible.

What changed

• stability_derivatives stores solver.reference_point as the pivot on body_aero (set_va!(...; reference_point)), 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 the apparent wind speed va.
• stability_derivatives, trim_angle and coeffs_at_angles leave body_aero.reference_point = solver.reference_point, and their docstrings say so: the omega columns are only rate derivatives when the body turns about the point the moments are taken about, so a separate pivot on body_aero cannot be kept through the call.
• 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).

Searched before writing

rg for span, c_ref, 2V, omega_idxs, rate, pb/2 and nondimensional in src/ and test/: nothing scaled rates already, and linearize's omega_idxs columns were the only rate Jacobian. The rate test reuses coeffs_at_angles (stepping body_aero.omega) rather than a helper of its own.

Found on the way

• test/solver/test_forwarddiff.jl's pivot testset (from set_va! turns the body about a stored reference_point, and set_va!(body_aero, settings) applies yaw_rate #353, now on main, outside this diff) builds a new BodyAerodynamics([wing]) per case sharing wing with body_aero, and has omega and loop variable omega_op in scope together. Test-only tidying; I would fold it into a cleanup: PR rather than widen this one.

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

#346 and #353 have both merged, so this now targets main and the diff is only src/stability.jl, its test and the changelog. Merging main after #368 renamed the speed argument wind_speed → va and the linearize keyword va_idxs → va_vec_idxs: src/stability.jl and CHANGELOG conflicted, resolved by keeping the rate derivatives and the pivot and taking #368's names; the new testsets were renamed the same way.

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 inside 11/11 for the trim testset. The existing α/β testset passes 5/5
• After merging main at 8754b6d (Add stability_derivatives (angle of attack, sideslip) and trim_angle, built on linearize #346, BREAKING: name the public apparent wind va, va_vec and va_dist, with no deprecation aliases #368): test_stability, test_forwarddiff, test_body_aerodynamics 4940/4940 pass on 1a30832 · docs build clean (docs/make.jl; only the existing functions.md / private_functions.md size warnings)
• After the review round 953b2d1: test_stability 27/27, docs build clean · Local full suite started on 1a30832: PASS (7 min, Julia 1.13.0, one matrix cell; 953b2d1 changes only src/stability.jl docstrings and test_stability.jl) · GitHub CI on 953b2d1: PASS (Documentation, Julia 1.12 on ubuntu/macOS/windows, Julia 1.13 on ubuntu, setup test, codecov)
• Risk: dp/dq/dr use va for V while the coefficients are normalised by the area-weighted RMS panel inflow speed, so at nonzero omega the p̂ scaling and the dynamic pressure use slightly different speeds. The docstring now states it.

Scope

+63 / −14 against main, across src/stability.jl, test/solver/test_stability.jl and CHANGELOG.md. Stack: 2/2 after #346, which merged.

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]

Kept the overwrite: the rates are only the rate derivatives when the body turns about the moment point. The docstrings of stability_derivatives, trim_angle and coeffs_at_angles now say they leave body_aero.reference_point = solver.reference_point, in 953b2d1.

[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]

Removed coeffs_at_rate: the test now steps body_aero.omega and calls coeffs_at_angles, the helper that replaced coeffs_at, in 953b2d1.

[1-Bort-1]

Local full suite: PASS (6 min, Julia 1.13.0, one cell of the matrix)

[1-Bart-1]

Resolve conflicts and address review comments

[1-Bort-1]

Conflicts with main (#346, #368 renames) resolved in 1a30832. Review 5229755456 addressed in 953b2d1: docstrings state the stored pivot and that V is va rather than the panel inflow speed, the rate test reuses coeffs_at_angles, dp/dq/dr wording rewritten. The two ForwardDiff-test points are about test_forwarddiff.jl, which came in with #353 and is not in this diff; noted in the description.