Take the panel normal from the quarter-chord step

CHANGELOG.md

@@ -1,5 +1,24 @@
 # Changelog
 
+## Unreleased
+
+### Fixed
+
+- `panel_axes` builds the panel normal from the quarter-chord step rather than
+  the leading-edge step. `z_airf` was `cross(x_airf, le_1 - le_2)` while the
+  bound vortex, `width` and `y_airf` all come from `panel_span_vector`, so on a
+  panel whose two sections have differently-directed chords the normal was not
+  square to the vortex the loads are built from. `alpha` is measured against
+  that normal, so every panel coefficient inherited the error. The frame now
+  closes as `z_airf = x_airf × y_airf`, which is what the `Panel` field
+  docstring already described.
+
+  Taper alone never triggered it: `(le_1 - le_2) - span_vec` is a quarter of
+  `chord_vec_2 - chord_vec_1`, which is purely chordwise when the two chords are
+  parallel, and `cross(x_airf, ·)` annihilates it. Twist, sweep and dihedral do
+  trigger it — 0.84° of twist on the `:rectangular` reference wing, and
+  5.6–13.8° on the raked outboard bay of a leading-edge-inflatable kite.
+
 ## VortexStepMethod v4.3.1 2026-09-01
 
 ### Changed

src/panel_aerodynamics.jl

@@ -49,6 +49,11 @@ end
 Airfoil frame and size of the panel between two sections, as
 `(; x_airf, y_airf, z_airf, chord, width)`.
 
+`x_airf` is chordwise and `y_airf` runs along the quarter-chord line the bound
+vortex sits on, so `z_airf` is their cross product. The first two are not orthogonal
+on a swept panel; taking the normal from the quarter-chord step rather than the
+leading-edge step is what keeps it square to the vortex the loads are built from.
+
 `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
@@ -63,7 +68,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)