Take the panel normal from the quarter-chord step

CHANGELOG.md

@@ -10,6 +10,11 @@
 
 - The `VSMSolution` docstring gives `lift_dist`, `drag_dist` and `panel_moment_dist` in
   the per-unit-span units they hold, [N/m] and [Nm/m], instead of [N] and [Nm].
+- `panel_axes` takes the panel normal from the quarter-chord step, so the frame
+  closes as `z_airf = x_airf × y_airf` and `z_airf` is square to the bound
+  vortex. `alpha` is measured against that normal, so `cl`, `cd` and `cm` were
+  wrong on panels whose two sections have differently-directed chords — twist,
+  sweep or dihedral, not taper alone.
 
 ## VortexStepMethod v5.1.1 2026-09-12
 

src/panel_aerodynamics.jl

@@ -47,13 +47,11 @@ end
     panel_axes(le_1, te_1, le_2, te_2, chord_weight=0.5, orient=1)
 
 Airfoil frame and size of the panel between two sections, as
-`(; x_airf, y_airf, z_airf, chord, width)`.
-
-`chord_weight` ([`panel_chord_weight`](@ref)) enters as an offset from the
-midpoint rather than as `w·p₁ + (1-w)·p₂`: the two are equal, but the offset form
-leaves a constant term to fold, which a symbolic consumer builds several times
-faster. `orient` is `±1`, flipping `y_airf`/`z_airf` so the frame does not depend
-on section ordering.
+`(; x_airf, y_airf, z_airf, chord, width)`. `x_airf` is chordwise, `y_airf` runs
+along the quarter-chord line the bound vortex sits on — not square to `x_airf` on
+a swept panel — and `z_airf` is `x_airf × y_airf` normalised. `chord_weight`
+([`panel_chord_weight`](@ref)) is section 1's share of the chord-direction blend;
+`orient` is `±1` and flips `y_airf` and `z_airf` together.
 """
 @inline function panel_axes(le_1, te_1, le_2, te_2, chord_weight=0.5, orient=1)
     lean = chord_weight - 0.5
@@ -63,7 +61,7 @@ on section ordering.
     width = smooth_norm(span_vec)
     x_airf = chord_vec ./ smooth_norm(chord_vec)
     y_airf = orient .* (span_vec ./ width)
-    z_cross = cross(x_airf, le_1 .- le_2)
+    z_cross = cross(x_airf, span_vec)
     z_airf = orient .* (z_cross ./ smooth_norm(z_cross))
     return (; x_airf, y_airf, z_airf,
             chord=panel_chord(le_1, te_1, le_2, te_2), width)

test/panel/test_panel.jl

@@ -1,4 +1,5 @@
 using VortexStepMethod: Panel, Section, calculate_relative_alpha_and_relative_velocity, calculate_cl, calculate_cd_cm, reinit!, INVISCID, POLAR_VECTORS, MVec3
+using VortexStepMethod: panel_axes, panel_span_vector
 using Interpolations: linear_interpolation, Line
 using LinearAlgebra
 using Test
@@ -167,4 +168,40 @@ end
             @test isapprox(cm, expected_cm, rtol=1e-5)
         end
     end
-end
+end
+
+@testset "Panel frame on a swept, twisted panel" begin
+    twist = deg2rad(35)
+    le_1 = [0.0, 0.0, 0.0]
+    te_1 = [2.0, 0.0, 0.0]
+    le_2 = [1.2, 1.0, 0.0]
+    te_2 = le_2 .+ 1.4 .* [cos(twist), 0.0, -sin(twist)]
+    axes = panel_axes(le_1, te_1, le_2, te_2)
+    span_vec = panel_span_vector(le_1, te_1, le_2, te_2)
+
+    leading_edge_normal = cross(axes.x_airf, le_1 .- le_2)
+    leading_edge_normal ./= norm(leading_edge_normal)
+    separation = rad2deg(acos(clamp(abs(dot(leading_edge_normal, axes.z_airf)), -1, 1)))
+
+    @testset "the geometry discriminates the two definitions" begin
+        @test separation > 5
+    end
+
+    @testset "the normal is square to the bound vortex and the chord" begin
+        @test abs(dot(axes.z_airf, span_vec)) < 1e-12
+        @test abs(dot(axes.z_airf, axes.y_airf)) < 1e-12
+        @test abs(dot(axes.z_airf, axes.x_airf)) < 1e-12
+    end
+
+    @testset "the frame closes as z = x cross y" begin
+        closure = cross(axes.x_airf, axes.y_airf)
+        @test isapprox(closure ./ norm(closure), axes.z_airf; atol=1e-12)
+    end
+
+    @testset "orient flips y and z together" begin
+        flipped = panel_axes(le_1, te_1, le_2, te_2, 0.5, -1)
+        @test isapprox(flipped.y_airf, -axes.y_airf; atol=1e-12)
+        @test isapprox(flipped.z_airf, -axes.z_airf; atol=1e-12)
+        @test isapprox(flipped.x_airf, axes.x_airf; atol=1e-12)
+    end
+end

test/thesis_oriol_cayon.jl

@@ -194,11 +194,11 @@ function create_geometry_general(coordinates, Uinf, N, ring_geo, model)
         )
         push!(filaments, bound)
 
-        x_airf = cross(VSMpoint - LLpoint, section["p2"] - section["p1"])
+        z_airf = bound["x2"] - bound["x1"]
+        x_airf = cross(VSMpoint - LLpoint, z_airf)
         x_airf = x_airf / norm(x_airf)
         y_airf = VSMpoint - LLpoint
         y_airf = y_airf / norm(y_airf)
-        z_airf = bound["x2"] - bound["x1"]
         z_airf = z_airf / norm(z_airf)
         airf_coord = hcat(x_airf, y_airf, z_airf)