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.

Figures

The TU Delft V3 kite from examples/V3_kite.jl, at 20 m/s, with moments about the origin. trim_angle sweeps the dots (2° apart), bisects the bracket where CMy changes sign, and returns one trim at 7.40° with slope −2.71 /rad, so the trim is stable. The dashed line is that slope. The ragged curve past about 12° is solve! in stall, not the trim search.

The same kite at β = 5°. Lines are stability_derivatives and dots are central differences of solve! with a step of 1e-4 rad. The α derivatives agree to 0.006 /rad everywhere. The β derivatives agree except at α = 15°, where every coefficient's finite difference jumps. The cause is solve!, not the derivative: at that α, stepping β by 1e-4 rad moves CFz between 0.9889 and 0.9982 and back, so the solve is landing on two stall branches. The ForwardDiff derivative there is −0.40 /rad, in line with its neighbours.

The two test wings in test_stability.jl

The figures come from a one-off script that is not committed. It calls only the exported API and coeffs_at_angles.

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, and exported. The examples, the README and the docs pages build their inflow with it instead of writing the formula out; they all ran at zero sideslip, where the vector is bit-identical, so no example's output changes. linearize_check.jl, bench.jl and stall_model.jl had left cos β off z, which only mattered had anyone set a sideslip. The old docstring labelled the 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 rather than steering the bisection with a number the solver does not stand behind.
• 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.

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, open as stability_derivatives gives roll, pitch and yaw rate derivatives about solver.reference_point #356 on top of this branch. They need the rotation taken about the reference point, which set_va! turns the body about a stored reference_point, and set_va!(body_aero, settings) applies yaw_rate #353 has since landed on main.
• linearize and make_dual_shadow accept a BodyAerodynamics with more... #344 let linearize and make_dual_shadow take a body with more than one wing, and landed as linearize and make_dual_shadow take a BodyAerodynamics with several wings #355. stability_derivatives and trim_angle take one too, since neither adds a restriction of its own; the multi-wing coverage is main's, in test_body_aerodynamics.jl.

Found on the way (not changed)

• CI does not stop at the first failure. .github/workflows/CI.yml sets fail-fast: false. Changing that is its own cleanup: PR.
• solve! jumps between branches in stall. On the V3 kite at α = 15°, β = 5°, steps of 1e-4 rad in β move CFz by about 0.009 and back (Figures). A trim_angle bracket there would bisect onto the jump rather than onto a smooth zero. The fix belongs in the solver's stall handling, not here.
• 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, and still 15/15 on the tree merged with main at 9d0f44e (Julia 1.13, 21 s).
  • apparent_wind points along body x, tilts up with alpha, right with beta (c8264df): 2 errored while apparent_wind was unexported, 3/3 after; the file is 18/18. examples/bench.jl runs to the end, and the docs build with apparent_wind on the public functions page.
  • The a solve that misses the tolerances throws testset was red on the code before b50a2b5: it returned a trim instead of throwing.
  • The 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: 4903/4903, and test/solver/test_solver.jl green, both on the merged tree. The first is where the set_va!(body_aero, settings) merge resolution is checked — set_va! with VSMSettings applies the yaw rate about body z is 6/6 here, and 1/6 with only this branch's side of that conflict, body_aero.omega coming out [0, 0, 0] against the expected [0, 0, 0.5236]. The second is the file Build a Solver from VSMSettings or from panel and section counts, deprecate the body_aero constructors, and check the sizes in solve! #340 rewrote, including Solver Constructor Tests 19/19.
• Merged with main, up to date at 9d0f44e. Two merges this round: Build a Solver from VSMSettings or from panel and section counts, deprecate the body_aero constructors, and check the sizes in solve! #340 deprecated Solver(body_aero; kwargs...), so the test helper now returns the body and a solver built from wing.n_panels and wing.n_unrefined_sections, as the rest of the suite does; no source file needed a change. The only conflict was the Added list in CHANGELOG.md, where both sides' entries are kept. Earlier merges also conflicted in set_va!(body_aero, settings), which now builds its inflow from apparent_wind and turns the body at settings.condition.yaw_rate from set_va! turns the body about a stored reference_point, and set_va!(body_aero, settings) applies yaw_rate #353.
• Docs build: clean on the merged tree, no missing docstrings; the only warnings are Documenter's page-size threshold and its deploy auto-detection, neither from this branch. No REUSE in this repo.
• Local CI mirror (Julia 1.13, fail-fast) on c8264df: PASS in 11 min, Aqua included. GitHub CI on c8264df: all seven checks green (35596033044), Julia 1.12 windows included. The test_forwarddiff.jl's POLAR_MATRICES check flakes at ~4%, and the oper #287 POLAR_MATRICES flake turned it red once on 0da60e3; dispatched rather than re-run, that same commit was green on all five jobs.
• Figures: rendered on c8264df in the worktree's session, from the examples environment.
• 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

+229 / −53 across 17 files:

• src/stability.jl (92 lines, about half of them docstrings)
• test/solver/test_stability.jl (97 lines)
• src/body_aerodynamics.jl: apparent_wind moved out of set_va!(settings), and its docstring cut from 29 lines to 4, net −15
• one line each in runtests.jl, the export list and the include list
• six examples, the README and two docs pages calling apparent_wind instead of the formula (net −12)
• docs and CHANGELOG entries

Closes #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.}

[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 (9 min, Julia 1.13.0, one cell of the matrix)

[1-Bart-1]

Examples still write the inflow formula by hand

Fix this in this PR

[1-Bort-1]

Fixed in c8264df: apparent_wind is exported, and the six examples, the README and the docs pages call it instead of writing out the formula.

[1-Bart-1]

I would like to see some figs in the PR body to make it easier to review.

[1-Bort-1]

Added three figures to the #346 description: the trim search on the V3 kite, the V3 derivatives against central differences, and the two test wings (collapsed).