Add stability_derivatives (angle of attack, sideslip) and trim_angle, built on linearize

[1-Bort-1]

TL;DR

New stability_derivatives(solver, body_aero, alpha, beta, wind_speed) returns the six force and moment coefficients and their derivatives with respect to angle of attack and sideslip. New trim_angle(solver, body_aero, beta, wind_speed) finds where CMy changes sign and returns the slope dCMy/dalpha at each trim, so you can see whether it is stable. The Python package has both and this package had neither. Both are built on linearize, as #330 proposes, rather than on a second finite-difference path.

How it works

• Derivatives. linearize(...; va_idxs=1:3, aero_coeffs=true) gives the coefficient Jacobian with respect to va. That is multiplied by dva/dalpha and dva/dbeta, which ForwardDiff takes from apparent_wind. Moments are about solver.reference_point, the rotation rate is whatever body_aero.omega holds, and extra keywords go to linearize (e.g. backend=nothing for a NONLIN solver, since ForwardDiff only runs on LOOP).
• Inflow formula. apparent_wind(alpha, beta, wind_speed) is the formula set_va!(body_aero, settings) already used, V·[cosα cosβ, sinβ, sinα cosβ], moved out so both callers share it. The old docstring labelled those components forward/right/down; the body frame is x back, z up, so the docstring now just names the formula.
• Trim search. It sweeps alpha_range (default −5° to 15° in 2° steps, in radians), bisects every bracket where CMy < 0 flips to alpha_tol (default 1e-5 rad), and takes each slope from stability_derivatives. It returns every trim found, and an empty vector if there is none.
• Trim failures. Every solve behind a trim runs with throw_on_fail: the sweep, the bisection and the slope's linearize. A solve that misses the solver's tolerances throws SolveFailure instead of steering the bisection. On the test wing with max_iterations=2, the version before the review silently returned a trim at 1.00° with slope −1.56; converged, the answer is 1.666° and −1.11.
• Trim keywords. trim_angle takes backend explicitly rather than passing keywords through to linearize. With pass-through, a reference_point= reached the slope but not the CMy sweep, so the trim and its stability verdict could be about two different points. Both review comments are addressed in b50a2b5.

Where this differs from the Python package

• Units. Angles are in radians, like the rest of this Julia API. Python takes degrees.
• Sideslip. The inflow keeps cos β on z, as set_va!(settings) has it. Python main uses sin α without it, so the numbers differ from Python at nonzero sideslip.
• Trim results. Python raises when it finds more than one stable trim, and falls back to "closest to zero" when it finds none. Here every sign change is returned, and none means an empty vector.
• Aircraft frame and Malz table. Not ported. The aircraft-frame mapping is six sign flips a caller can apply. The Malz table is broken on Python main (Add rigid-body stability derivatives and a trim-angle search, as in the Python package #330 lists why) and nothing here asked for it.

Split

#330 holds more than one idea, so this is part 1:

• stability_derivatives gives roll, pitch and yaw rate derivatives... #345 adds the p/q/r rate derivatives. They need the rotation taken about the reference point (set_va! rotates the body about the origin only, and condition.yaw_rate is read but never applied #329). Until then linearize's omega columns rotate about the origin while moments are taken about solver.reference_point, so porting the rates now would mean a second copy of ω × (r − r0).
• linearize and make_dual_shadow accept a BodyAerodynamics with more... #344 lets linearize and make_dual_shadow take more than one wing. Until it lands, stability_derivatives throws on a multi-wing body, and so does trim_angle once it has found a trim.

Found on the way (not changed)

• Examples still write the inflow formula by hand. examples/billowing.jl has the same formula. examples/linearize_check.jl leaves cos β off z. apparent_wind is private, so I left the examples alone.
• CI does not stop at the first failure. .github/workflows/CI.yml sets fail-fast: false. Changing that is its own cleanup: PR.
• Before searching: I searched src/, test/, examples/ and docs/ for stability, derivative, trim, malz and bisect. Nothing to reuse; the one bisection, in obj_slice.jl, carries its own per-step payload.

Verification

• Reproduction: n/a, new feature.
• test/solver/test_stability.jl, red before, green after.
  • Before the implementation, include gave 0 passed, 0 failed, 1 errored because stability_derivatives, trim_angle and apparent_wind were undefined.
  • After: 15/15 in juliaserver.
  • The new a solve that misses the tolerances throws testset was red on the code before b50a2b5: it returned a trim instead of throwing.
  • The new a NONLIN solver with backend=nothing finds the same trim testset covers the finite-difference path. It matches the LOOP + ForwardDiff trim to 1e-4.
  • On a rectangular POLAR_VECTORS wing, the derivatives agree with central differences of solve! (step 1e-4 rad) to about 1e-8 relative, e.g. dCFz/dα 4.414345384 against 4.414345369. Tested to rtol 1e-4.
  • With section cm = 0.05 and moments about the leading edge, it finds one trim at 1.666° with slope −1.11 /rad. With cm = −0.05 and moments about the trailing edge, one trim at −0.050° with slope +3.34 /rad.
• test/body_aerodynamics/test_body_aerodynamics.jl: green, including set_va! with VSMSettings (5/5), which covers the moved formula.
• Merged with main, up to date at 83867a3.
  • Add the Gaunaa et al. 2024 spanwise-flow viscous drag correction as an opt-in solver setting #352 (viscous drag correction): the only conflict was the two Added entries in CHANGELOG.md; both are kept. test_viscous_drag_correction.jl passed 90/90 on that merge.
  • Name the apparent wind va, va_vec and va_dist everywhere except the public API #349 (va rename): the conflict was the body of set_va!(body_aero, settings). It now calls apparent_wind, and the vector local is va_vec, as in stability_derivatives. test_body_aerodynamics.jl is green on the merged tree.
• Docs build: clean at b50a2b5, no missing docstrings. No REUSE in this repo.
• Local CI mirror on b50a2b5 (Julia 1.12, fail-fast): PASS in 9 min, Aqua included. GitHub CI on b50a2b5: PASS (all six jobs and codecov/patch).
• jetls is not installed on the box, so it was not run.
• Risk: dCFy/dβ is exactly zero on a flat wing, so the side-force derivative is only checked against the finite difference where it is zero. The other five sideslip derivatives are nonzero and checked.

Scope

+209 / −25 across 8 files:

• src/stability.jl (92 lines, about half of them docstrings)
• test/solver/test_stability.jl (92 lines)
• src/body_aerodynamics.jl: apparent_wind moved out of set_va!(settings), and its docstring cut from 25 lines to 3, net −12
• one line each in runtests.jl, the export list and the include list
• docs and CHANGELOG entries

Stack: 1/3 → #345 (rate derivatives, blocked on #329), #344 (multi-wing linearize).

Refs #330 · task VortexStepMethod.jl-330

[1-Bort-1]

Independent review (advisory)

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

Good

• The derivatives are built on linearize, as the card promises: va_idxs=1:3 with theta_idxs=nothing, and omega read from body_aero.omega at src/solver.jl:1319-1321, so there is no second finite-difference path.
• The formula in apparent_wind matches the removed body of set_va!(settings) term for term, and that method now calls it, so the inflow has one source.
• dva/dalpha and dva/dbeta come from ForwardDiff over apparent_wind, so the formula's derivative is not written out a second time.
• The docstring fix fits the geometry: the test wing runs from the leading edge at x=0 to the trailing edge at x=1, and docs/src/settings.md:54 has z = up, so the old forward/down labels were wrong.
• The tests name what they protect, check both slope signs and a range with no trim, and compare against central differences of solve! rather than against linearize itself.
• Public symbols are in functions.md, private helpers in private_functions.md, and CHANGELOG and export and include lines are added, matching the plan in the card.

Not good

• src/stability.jl:46 — kwargs reach linearize and so solve! for the slope, but not the CMy sweep in pitch_moment_coeff. A caller passing reference_point= (a solve! keyword) gets trims found about solver.reference_point but slopes about another point, so the stable/unstable verdict can be silently wrong.
• src/stability.jl:61 — pitch_moment_coeff never checks solver.lr.converged, and line 48 throws away derivatives.converged. A trim bisected on unconverged CMy values comes back looking just like a good one.
• coeffs_at in the test is pitch_moment_coeff returning all six coefficients; one private src helper returning all six, with trim_angle taking [5], would serve both.
• The logic of the trim predicate sits in the closure is_nose_down (§6); a named helper, or pitch_moment_coeff(...) < 0 passed straight in, reads the same without it.
• The trim tests use the default Solver, where use_gamma_prev is on, so each bisection step starts from the last step's gamma; the derivative test turns it off, and the trim tests do not say why they leave it on.
• The card says to pass backend=nothing for a NONLIN solver, but no test covers that path through stability_derivatives or trim_angle.
• trims is fixed to Float64, so a Float32 or dual alpha_range is converted without warning; fine for now, but worth a line in the docstring if it is intended.

_{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: kwargs reach linearize and so solve! for the slope, but not the CMy sweep in pitch_moment_coeff. A caller passing reference_point= (a solve! keyword) gets trims found about solver.reference_point but slopes about another point, so the stable/unstable verdict can be silently wrong.

[1-Bort-1]

Fixed in b50a2b5: trim_angle takes an explicit backend instead of splatting kwargs, so the sweep, the bisection and the slope all take moments about solver.reference_point.

[1-Bort-1]

MINOR: pitch_moment_coeff never checks solver.lr.converged, and line 48 throws away derivatives.converged. A trim bisected on unconverged CMy values comes back looking just like a good one.

[1-Bort-1]

Fixed in b50a2b5: the sweep, bisection and slope solves all run with throw_on_fail, so an unconverged solve throws SolveFailure instead of returning a trim (new testset, red on the old code: max_iterations=2 returned a trim at 1.00° with slope -1.56).

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.

📢 Thoughts on this report? Let us know!

[1-Bort-1]

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