From 7403315dcd603dd2d6ea107a378194e50b84b720 Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Thu, 24 Sep 2026 10:36:50 +0200
Subject: [PATCH 1/3] Measure alpha in the plane square to the bound vortex

Every alpha the solver forms (panel_inflow, the circulation loop, the
geometric and aerodynamic-centre angles) and the inflow direction that turns
the force now go through inflow_angle, which dots the velocity with the chord
projected onto the plane square to y_airf: Gaunaa, Li & Pirrung (2026), CP1.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/body_aerodynamics.jl   | 59 ++++++-----------------------
 src/panel_aerodynamics.jl  | 17 +++++++--
 src/solver.jl              | 77 +++++++-------------------------------
 test/panel/test_panel.jl   | 15 +++++++-
 test/solver/test_solver.jl | 33 ++++++++++++++++
 5 files changed, 85 insertions(+), 116 deletions(-)

diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index bd390826..82ec509d 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -523,20 +523,16 @@ end
 
 """
     update_effective_angle_of_attack!(alpha_corrected, body_aero::BodyAerodynamics, gamma,
-                                      core_radius_fraction, z_airf_dist, x_airf_dist,
-                                      va_vec_dist, va_dist, va_unit_dist)
+                                      core_radius_fraction, va_vec_dist, va_dist,
+                                      va_unit_dist)
 
-Update angle of attack at aerodynamic center for VSM method.
-
-Returns:
-    nothing
+Write each panel's [`inflow_angle`](@ref) at its aerodynamic centre into
+`alpha_corrected`, the induced velocity taken from `gamma` on the LLT influence matrix.
 """
 function update_effective_angle_of_attack!(alpha_corrected,
-    body_aero::BodyAerodynamics, 
+    body_aero::BodyAerodynamics,
     gamma,
     core_radius_fraction,
-    z_airf_dist,
-    x_airf_dist,
     va_vec_dist,
     va_dist,
     va_unit_dist)
@@ -551,31 +547,11 @@ function update_effective_angle_of_attack!(alpha_corrected,
         mul!(view(induced_velocity, :, k), view(body_aero.AIC_aero_center, :, :, k), gamma)
     end
 
-    # In-place relative velocity calculation
-    relative_velocity = body_aero.cache[2][va_vec_dist]
-    relative_velocity .= va_vec_dist .+ induced_velocity
-
-    # Preallocate and compute dot products manually
-    n = size(relative_velocity, 1)
-    v_normal     = body_aero.cache[3][relative_velocity]
-    v_tangential = body_aero.cache[4][relative_velocity]
-    
-    @inbounds for i in 1:n
-        vn = 0.0
-        vt = 0.0
-        for j in 1:3
-            vn += z_airf_dist[i, j] * relative_velocity[i, j]
-            vt += x_airf_dist[i, j] * relative_velocity[i, j]
-        end
-        v_normal[i] = vn
-        v_tangential[i] = vt
-    end
-
-    # Direct angle calculation without temporary arrays
-    @inbounds for i in 1:n
-        alpha_corrected[i] = atan(v_normal[i], v_tangential[i])
+    @inbounds for (i, panel) in enumerate(body_aero.panels)
+        v_rel = matrix_row(va_vec_dist, i) .+ matrix_row(induced_velocity, i)
+        alpha_corrected[i] = inflow_angle(v_rel, SVector{3}(panel.y_airf),
+                                          SVector{3}(panel.z_airf))
     end
-
     nothing
 end
 
@@ -861,19 +837,8 @@ function calculate_results(
                 body_aero.pitch_rate_dist[i], chord_dist[i], v_rel_dist[i])
         end
         panel_width_dist[i] = panel.width
-        va = va_dist[i]
-        x_norm = norm3(panel.x_airf)
-        z_norm = norm3(panel.z_airf)
-        if va == 0.0 || x_norm == 0.0 || z_norm == 0.0
-            alpha_geometric[i] = NaN
-        else
-            inv_va = 1.0 / va
-            v_tangential = -dot3(panel.x_airf, panel.va_vec) *
-                           inv_va / x_norm
-            v_normal = -dot3(panel.z_airf, panel.va_vec) *
-                       inv_va / z_norm
-            alpha_geometric[i] = atan(-v_normal, -v_tangential)
-        end
+        alpha_geometric[i] = iszero(va_dist[i]) ? NaN : inflow_angle(panel.va_vec,
+            SVector{3}(panel.y_airf), SVector{3}(panel.z_airf))
     end
 
     # Calculate alpha corrections based on model type
@@ -883,8 +848,6 @@ function calculate_results(
             body_aero,
             gamma_new,
             core_radius_fraction,
-            z_airf_dist,
-            x_airf_dist,
             va_vec_dist,
             va_dist,
             va_unit_dist
diff --git a/src/panel_aerodynamics.jl b/src/panel_aerodynamics.jl
index 46cfba89..edf26e44 100644
--- a/src/panel_aerodynamics.jl
+++ b/src/panel_aerodynamics.jl
@@ -95,6 +95,15 @@ already sampled at three-quarter chord.
     ifelse(v_rel > 0, -0.25π * pitch_rate * chord / (2 * smooth_norm(v_rel)),
            zero(pitch_rate * chord))
 
+"""
+    inflow_angle(velocity, y_airf, z_airf)
+
+Angle of attack of `velocity` in the plane square to `y_airf`: from the chord's
+projection on that plane, `y_airf × z_airf`, towards `z_airf`.
+"""
+@inline inflow_angle(velocity, y_airf, z_airf) =
+    atan(dot(velocity, z_airf), dot(velocity, cross(y_airf, z_airf)))
+
 """
     effective_alpha(alpha, deficiency)
 
@@ -115,9 +124,9 @@ leaves it zero. `deficiency` feeds [`effective_alpha`](@ref).
 """
 @inline function panel_inflow(axes, va_vec_1, va_vec_2, v_ind, dva_vec_1=nothing,
                               dva_vec_2=nothing, deficiency=0)
-    (; x_airf, y_airf, z_airf, chord) = axes
+    (; y_airf, z_airf, chord) = axes
     v_eff = 0.5 .* (va_vec_1 .+ va_vec_2) .+ v_ind
-    alpha = atan(dot(v_eff, z_airf), dot(v_eff, x_airf))
+    alpha = inflow_angle(v_eff, y_airf, z_airf)
     v_span = cross(v_eff, y_airf)
     pitch_rate = isnothing(dva_vec_1) ? zero(chord) :
         section_pitch_rate(0.5 .* (dva_vec_1 .+ dva_vec_2), z_airf, chord)
@@ -143,8 +152,8 @@ Lift and drag unit vectors, `(; dir_lift, dir_drag)`. `alpha_dir` is the angle t
 turns the force, `spanwise` the wing's spanwise direction.
 """
 @inline function panel_force_directions(axes, alpha_dir, spanwise)
-    (; x_airf, y_airf, z_airf) = axes
-    dir_inflow = cos(alpha_dir) .* x_airf .+ sin(alpha_dir) .* z_airf
+    (; y_airf, z_airf) = axes
+    dir_inflow = cos(alpha_dir) .* cross(y_airf, z_airf) .+ sin(alpha_dir) .* z_airf
     lift_cross = cross(dir_inflow, y_airf)
     dir_lift = lift_cross ./ smooth_norm(lift_cross)
     drag_cross = cross(spanwise, dir_lift)
diff --git a/src/solver.jl b/src/solver.jl
index d6f07f34..ac48e75c 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -178,7 +178,7 @@ sol::VSMSolution = VSMSolution(): The result of calling [`solve!`](@ref)
     # Intermediate results
     lr::LoopResult{P, T} = LoopResult{P, T}()
     br::BaseResult{P, T} = BaseResult{P, T}()
-    cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:10]
+    cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:8]
     cache_base::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}}  = [LazyBufferCache()]
     cache_lin::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:4]
 
@@ -373,35 +373,9 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
         end
         width_dist[i] = panel.width
 
-        # Geometric AoA using panel-local axes and prescribed
-        # freestream — scalar ops to avoid allocations
-        begin
-            va1 = solver.sol.va_vec_dist[i,1]
-            va2 = solver.sol.va_vec_dist[i,2]
-            va3 = solver.sol.va_vec_dist[i,3]
-            va = sqrt(va1^2 + va2^2 + va3^2)
-            x1 = solver.sol._x_airf_dist[i,1]
-            x2 = solver.sol._x_airf_dist[i,2]
-            x3 = solver.sol._x_airf_dist[i,3]
-            x_norm = sqrt(x1^2 + x2^2 + x3^2)
-            z1 = solver.sol._z_airf_dist[i,1]
-            z2 = solver.sol._z_airf_dist[i,2]
-            z3 = solver.sol._z_airf_dist[i,3]
-            z_norm = sqrt(z1^2 + z2^2 + z3^2)
-            if va == 0 || x_norm == 0 || z_norm == 0
-                alpha_geometric_dist[i] = NaN
-            else
-                inv_va = -1.0 / va
-                vu1 = va1 * inv_va
-                vu2 = va2 * inv_va
-                vu3 = va3 * inv_va
-                v_tangential = (x1*vu1+x2*vu2+x3*vu3) / x_norm
-                v_normal = (z1*vu1+z2*vu2+z3*vu3) / z_norm
-                alpha_geometric_dist[i] = atan(
-                    -v_normal, -v_tangential)
-            end
-        end
-
+        va_vec = matrix_row(solver.sol.va_vec_dist, i)
+        alpha_geometric_dist[i] = iszero(va_vec) ? NaN : inflow_angle(va_vec,
+            matrix_row(solver.sol._y_airf_dist, i), matrix_row(solver.sol._z_airf_dist, i))
     end
 
     lift = solver.sol.lift_dist
@@ -415,8 +389,6 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
             body_aero,
             gamma_new,
             solver.core_radius_fraction,
-            solver.sol._z_airf_dist,
-            solver.sol._x_airf_dist,
             solver.sol.va_vec_dist,
             solver.br.va_dist,
             solver.br.va_unit_dist
@@ -756,6 +728,9 @@ end
 
 @inline smooth_sqrt(x) = sqrt(x + 1e-30)
 
+"""Row `i` of an `n×3` matrix as a 3-vector."""
+@inline matrix_row(matrix, i) = SVector(matrix[i, 1], matrix[i, 2], matrix[i, 3])
+
 @inline function update_gamma_candidate!(
     gamma_out,
     gamma_in,
@@ -775,9 +750,6 @@ end
     relative_velocity_crossz,
     v_acrossz_dist,
     z_airf_dist,
-    x_airf_dist,
-    v_normal_dist,
-    v_tangential_dist,
     va_magw_dist,
     cl_dist,
     chord_dist,
@@ -806,20 +778,11 @@ end
     end
 
     @inbounds for i in 1:n_panels
-        v_normal_dist[i] =
-            z_airf_dist[i,1]*relative_velocity_dist[i,1] +
-            z_airf_dist[i,2]*relative_velocity_dist[i,2] +
-            z_airf_dist[i,3]*relative_velocity_dist[i,3]
-        v_tangential_dist[i] =
-            x_airf_dist[i,1]*relative_velocity_dist[i,1] +
-            x_airf_dist[i,2]*relative_velocity_dist[i,2] +
-            x_airf_dist[i,3]*relative_velocity_dist[i,3]
-        solver.lr.v_span_dist[i] =
-            y_airf_dist[i,1]*relative_velocity_dist[i,1] +
-            y_airf_dist[i,2]*relative_velocity_dist[i,2] +
-            y_airf_dist[i,3]*relative_velocity_dist[i,3]
+        v_rel = matrix_row(relative_velocity_dist, i)
+        y_airf = matrix_row(y_airf_dist, i)
+        solver.lr.alpha_dist[i] = inflow_angle(v_rel, y_airf, matrix_row(z_airf_dist, i))
+        solver.lr.v_span_dist[i] = dot(y_airf, v_rel)
     end
-    solver.lr.alpha_dist .= atan.(v_normal_dist, v_tangential_dist)
 
     @inbounds for i in 1:n_panels
         solver.lr.v_rel_dist[i] = smooth_sqrt(
@@ -945,7 +908,6 @@ function gamma_loop!(
 ) where {P, U, T}
     va_vec_dist = solver.sol.va_vec_dist
     chord_dist = solver.sol._chord_dist
-    x_airf_dist = solver.sol._x_airf_dist
     y_airf_dist = solver.sol._y_airf_dist
     z_airf_dist = solver.sol._z_airf_dist
     solver.lr.converged   = false
@@ -961,8 +923,6 @@ function gamma_loop!(
     relative_velocity_crossz = solver.cache[6][va_vec_dist]
     v_acrossz_dist           = solver.cache[7][va_vec_dist]
     cl_dist                  = solver.cache[8][solver.lr.gamma_new]
-    v_normal_dist            = solver.cache[9][solver.lr.gamma_new]
-    v_tangential_dist        = solver.cache[10][solver.lr.gamma_new]
 
     AIC_x = @view body_aero.AIC[:, :, 1]
     AIC_y = @view body_aero.AIC[:, :, 2]
@@ -991,9 +951,7 @@ function gamma_loop!(
             velocity_view_x, velocity_view_y, velocity_view_z,
             va_vec_dist, induced_velocity_all, relative_velocity_dist,
             y_airf_dist, relative_velocity_crossz, v_acrossz_dist,
-            z_airf_dist, x_airf_dist,
-            v_normal_dist, v_tangential_dist,
-            va_magw_dist, cl_dist, chord_dist,
+            z_airf_dist, va_magw_dist, cl_dist, chord_dist,
         )
         @inbounds for i in 1:n_panels
             residual[i] -= gamma_iter[i]
@@ -1013,9 +971,7 @@ function gamma_loop!(
                     velocity_view_x, velocity_view_y, velocity_view_z,
                     va_vec_dist, induced_velocity_all, relative_velocity_dist,
                     y_airf_dist, relative_velocity_crossz, v_acrossz_dist,
-                    z_airf_dist, x_airf_dist,
-                    v_normal_dist, v_tangential_dist,
-                    va_magw_dist, cl_dist, chord_dist,
+                    z_airf_dist, va_magw_dist, cl_dist, chord_dist,
                 )
                 inv_step = 1.0 / step
                 for i in 1:n_panels
@@ -1042,9 +998,7 @@ function gamma_loop!(
                     velocity_view_x, velocity_view_y, velocity_view_z,
                     va_vec_dist, induced_velocity_all, relative_velocity_dist,
                     y_airf_dist, relative_velocity_crossz, v_acrossz_dist,
-                    z_airf_dist, x_airf_dist,
-                    v_normal_dist, v_tangential_dist,
-                    va_magw_dist, cl_dist, chord_dist,
+                    z_airf_dist, va_magw_dist, cl_dist, chord_dist,
                 )
                 @inbounds for i in 1:n_panels
                     residual_perturbed[i] -= gamma_perturbed[i]
@@ -1108,9 +1062,6 @@ function gamma_loop!(
                 relative_velocity_crossz,
                 v_acrossz_dist,
                 z_airf_dist,
-                x_airf_dist,
-                v_normal_dist,
-                v_tangential_dist,
                 va_magw_dist,
                 cl_dist,
                 chord_dist,
diff --git a/test/panel/test_panel.jl b/test/panel/test_panel.jl
index 63fc8231..dca6d354 100644
--- a/test/panel/test_panel.jl
+++ b/test/panel/test_panel.jl
@@ -1,5 +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 VortexStepMethod: panel_axes, panel_span_vector, panel_inflow
 using Interpolations: linear_interpolation, Line
 using LinearAlgebra
 using Test
@@ -215,4 +215,17 @@ end
         @test isapprox(flipped.z_airf, -axes.z_airf; atol=1e-12)
         @test isapprox(flipped.x_airf, axes.x_airf; atol=1e-12)
     end
+
+    @testset "alpha is measured in the plane square to the bound vortex" begin
+        chord_in_plane = normalize(axes.x_airf .- dot(axes.x_airf, axes.y_airf) .* axes.y_airf)
+        for v_eff in ([1.0, 0.0, 0.0], [1.0, 0.3, 0.1]), orient in (1, -1)
+            oriented = panel_axes(le_1, te_1, le_2, te_2, 0.5, orient)
+            alpha = panel_inflow(oriented, v_eff, v_eff, zeros(3)).alpha
+            @test alpha ≈ atan(dot(v_eff, oriented.z_airf), dot(v_eff, chord_in_plane))
+            along_span = v_eff .+ 2.0 .* axes.y_airf
+            @test panel_inflow(oriented, along_span, along_span, zeros(3)).alpha ≈ alpha
+        end
+        @test rad2deg(panel_inflow(axes, [1.0, 0, 0], [1.0, 0, 0], zeros(3)).alpha) ≈
+              -20.097 atol = 1e-3
+    end
 end
diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl
index 089b3353..adcb1898 100644
--- a/test/solver/test_solver.jl
+++ b/test/solver/test_solver.jl
@@ -328,3 +328,36 @@ end
     @test VortexStepMethod.finite_full(ForwardDiff.Dual(1.0, 2.0))
     @test !VortexStepMethod.finite_full(ForwardDiff.Dual(1.0, Inf))
 end
+
+@testset "alpha is the crossflow angle on a swept, tapered wing" begin
+    polar = (deg2rad.([-20.0, 0.0, 20.0]), [-1.8, 0.2, 2.2], fill(0.02, 3), fill(-0.05, 3))
+    wing = Wing(8; spanwise_distribution=LINEAR)
+    add_section!(wing, [1.5, 4.0, 0.3], [2.2, 4.0, 0.1], POLAR_VECTORS, polar)
+    add_section!(wing, [0.0, 0.0, 0.0], [2.0, 0.0, 0.0], POLAR_VECTORS, polar)
+    add_section!(wing, [1.5, -4.0, 0.3], [2.2, -4.0, 0.1], POLAR_VECTORS, polar)
+    refine!(wing)
+    body_aero = BodyAerodynamics([wing])
+    va = [15.0, 1.0, 1.5]
+    set_va!(body_aero, va)
+
+    function crossflow_alpha(panel, velocity)
+        chord_in_plane = panel.x_airf .- dot(panel.x_airf, panel.y_airf) .* panel.y_airf
+        return atan(dot(velocity, panel.z_airf), dot(velocity, normalize(chord_in_plane)))
+    end
+    induced(AIC, gamma, i) = [dot(AIC[i, :, k], gamma) for k in 1:3]
+
+    for correct_aoa in (false, true)
+        solver = Solver(length(body_aero.panels), 3; correct_aoa)
+        sol = solve!(solver, body_aero)
+        @test sol.solver_status == FEASIBLE
+        gamma = solver.lr.gamma_new
+        for (i, panel) in enumerate(body_aero.panels)
+            @test solver.lr.alpha_dist[i] ≈
+                  crossflow_alpha(panel, va .+ induced(body_aero.AIC, gamma, i)) atol = 1e-6
+            @test sol.alpha_geometric_dist[i] ≈ crossflow_alpha(panel, va) atol = 1e-10
+            correct_aoa || continue
+            @test sol.alpha_dist[i] ≈ crossflow_alpha(panel,
+                va .+ induced(body_aero.AIC_aero_center, gamma, i)) atol = 1e-8
+        end
+    end
+end

From 8aa6c079f5e4a3fb4db281da2cf536e99549278d Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Thu, 24 Sep 2026 10:46:02 +0200
Subject: [PATCH 2/3] Document inflow_angle and matrix_row; changelog the
 crossflow alpha

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                  | 4 ++++
 docs/src/private_functions.md | 2 ++
 test/panel/test_panel.jl      | 3 ++-
 3 files changed, 8 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index dd9c7964..cf3905c0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,10 @@
   the right tip, z = x × y up. The separate `KB` definition, whose `Z = Y × X` pointed
   down, is gone; the solver's frame is unchanged. CL, CD and CS are documented as
   wind-axis coefficients, `cfx`, `cfy` and `cfz` as the body-axis ones.
+- The angle of attack is measured in the plane square to the quarter-chord line, from
+  the chord's projection on that plane (Gaunaa, Li & Pirrung 2026, CP1), everywhere it
+  is formed and where it turns the lift direction. Swept panels read their polars at a
+  different angle than before; unswept panels are unchanged.
 
 ## VortexStepMethod v6.0.0 2026-09-23
 
diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index f7482533..c5b1c873 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -41,6 +41,7 @@ calculate_circulation_distribution_elliptical_wing
 _compute_reference_velocity_from_distribution
 smooth_distribution!
 make_dual_shadow
+matrix_row
 ```
 
 ### Panel aerodynamics
@@ -54,6 +55,7 @@ smooth_norm
 panel_span_vector
 panel_chord_weight
 panel_axes
+inflow_angle
 effective_alpha
 panel_inflow
 dynamic_pressure
diff --git a/test/panel/test_panel.jl b/test/panel/test_panel.jl
index dca6d354..113b75f4 100644
--- a/test/panel/test_panel.jl
+++ b/test/panel/test_panel.jl
@@ -217,7 +217,8 @@ end
     end
 
     @testset "alpha is measured in the plane square to the bound vortex" begin
-        chord_in_plane = normalize(axes.x_airf .- dot(axes.x_airf, axes.y_airf) .* axes.y_airf)
+        chord_in_plane = normalize(axes.x_airf .-
+                                   dot(axes.x_airf, axes.y_airf) .* axes.y_airf)
         for v_eff in ([1.0, 0.0, 0.0], [1.0, 0.3, 0.1]), orient in (1, -1)
             oriented = panel_axes(le_1, te_1, le_2, te_2, 0.5, orient)
             alpha = panel_inflow(oriented, v_eff, v_eff, zeros(3)).alpha

From 746e79bb8120ec10edfd006fac974997ff41073e Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Thu, 24 Sep 2026 11:52:05 +0200
Subject: [PATCH 3/3] Trim the results cache, share the crossflow oracle, test
 the lift direction

BodyAerodynamics.cache drops the three slots nothing reads since alpha
went through inflow_angle, matrix_row moves next to inflow_angle, and the
unused duplicate calculate_relative_alpha_and_velocity goes. The panel
test now checks that dir_lift is square to the crossflow velocity, and
both alpha tests take crossflow_alpha from test_data_utils.jl.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 docs/src/private_functions.md |  1 -
 src/body_aerodynamics.jl      | 30 +++++++++++++++---------------
 src/panel.jl                  | 10 ----------
 src/panel_aerodynamics.jl     |  3 +++
 src/solver.jl                 |  3 ---
 test/panel/test_panel.jl      | 16 +++++++++++-----
 test/solver/test_solver.jl    | 17 +++++++----------
 test/test_data_utils.jl       | 11 +++++++++++
 8 files changed, 47 insertions(+), 44 deletions(-)

diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index c5b1c873..a498dff0 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -31,7 +31,6 @@ set_pitch_rate_dist!
 coeffs_at_angles
 nose_down
 bisect_sign_change
-calculate_relative_alpha_and_velocity
 calculate_relative_alpha_and_relative_velocity
 update_effective_angle_of_attack!
 calculate_stall_angle_list
diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index 82ec509d..e37caa70 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -28,7 +28,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru
                         influence coefficients, used only for the corrected angle of attack
 - `projected_area::Float64` = 1.0: The area projected onto the xy-plane of the KA body frame [m²]
 - `c_ref::Float64` = 1.0: Reference chord length (max panel chord) [m]
-- `cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}}` = [LazyBufferCache() for _ in 1:15]
+- `cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}}` = [LazyBufferCache() for _ in 1:12]
 """
 @with_kw mutable struct BodyAerodynamics{P, W<:AbstractWing, T, PN<:Panel{T}}
     panels::Vector{PN}
@@ -49,7 +49,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru
     AIC_aero_center::Array{T, 3} = zeros(T, P, P, 3)
     projected_area::T = one(T)
     c_ref::T = one(T)
-    cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:15]
+    cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:12]
 end
 
 """
@@ -808,21 +808,21 @@ function calculate_results(
 )
 
     n_panels = length(panels)
-    if length(body_aero.cache) < 15
-        append!(body_aero.cache, [LazyBufferCache() for _ in 1:(15 - length(body_aero.cache))])
+    if length(body_aero.cache) < 12
+        append!(body_aero.cache, [LazyBufferCache() for _ in 1:(12 - length(body_aero.cache))])
     end
 
-    cl_dist = body_aero.cache[5][alpha_dist]
-    cd_dist = body_aero.cache[6][alpha_dist]
-    cm_dist = body_aero.cache[7][alpha_dist]
-    panel_width_dist = body_aero.cache[8][alpha_dist]
-    alpha_corrected = body_aero.cache[9][alpha_dist]
-    cl_prescribed_va = body_aero.cache[10][alpha_dist]
-    cd_prescribed_va = body_aero.cache[11][alpha_dist]
-    cs_prescribed_va = body_aero.cache[12][alpha_dist]
-    f_body_3D = body_aero.cache[13][alpha_dist, (3, length(alpha_dist))]
-    m_body_3D = body_aero.cache[14][alpha_dist, (3, length(alpha_dist))]
-    alpha_geometric = body_aero.cache[15][alpha_dist]
+    cl_dist = body_aero.cache[2][alpha_dist]
+    cd_dist = body_aero.cache[3][alpha_dist]
+    cm_dist = body_aero.cache[4][alpha_dist]
+    panel_width_dist = body_aero.cache[5][alpha_dist]
+    alpha_corrected = body_aero.cache[6][alpha_dist]
+    cl_prescribed_va = body_aero.cache[7][alpha_dist]
+    cd_prescribed_va = body_aero.cache[8][alpha_dist]
+    cs_prescribed_va = body_aero.cache[9][alpha_dist]
+    f_body_3D = body_aero.cache[10][alpha_dist, (3, length(alpha_dist))]
+    m_body_3D = body_aero.cache[11][alpha_dist, (3, length(alpha_dist))]
+    alpha_geometric = body_aero.cache[12][alpha_dist]
 
     fill!(f_body_3D, 0.0)
     fill!(m_body_3D, 0.0)
diff --git a/src/panel.jl b/src/panel.jl
index 3a211857..531ee1be 100644
--- a/src/panel.jl
+++ b/src/panel.jl
@@ -479,16 +479,6 @@ function calculate_relative_alpha_and_relative_velocity(
     return flow.alpha, flow.v_eff
 end
 
-"""
-    calculate_relative_alpha_and_velocity(panel::Panel, induced_velocity)
-
-Calculate relative angle of attack and relative velocity of the panel.
-"""
-function calculate_relative_alpha_and_velocity(panel::Panel, induced_velocity)
-    flow = panel_inflow(panel_axes(panel), panel.va_vec, panel.va_vec, induced_velocity)
-    return flow.alpha, flow.v_eff
-end
-
 """
     calculate_cl(panel::Panel, alpha)
     calculate_cl(panel::Panel, alpha, delta)
diff --git a/src/panel_aerodynamics.jl b/src/panel_aerodynamics.jl
index edf26e44..8aabc7e8 100644
--- a/src/panel_aerodynamics.jl
+++ b/src/panel_aerodynamics.jl
@@ -104,6 +104,9 @@ projection on that plane, `y_airf × z_airf`, towards `z_airf`.
 @inline inflow_angle(velocity, y_airf, z_airf) =
     atan(dot(velocity, z_airf), dot(velocity, cross(y_airf, z_airf)))
 
+"""Row `i` of an `n×3` matrix as a 3-vector."""
+@inline matrix_row(matrix, i) = SVector(matrix[i, 1], matrix[i, 2], matrix[i, 3])
+
 """
     effective_alpha(alpha, deficiency)
 
diff --git a/src/solver.jl b/src/solver.jl
index ac48e75c..508de43f 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -728,9 +728,6 @@ end
 
 @inline smooth_sqrt(x) = sqrt(x + 1e-30)
 
-"""Row `i` of an `n×3` matrix as a 3-vector."""
-@inline matrix_row(matrix, i) = SVector(matrix[i, 1], matrix[i, 2], matrix[i, 3])
-
 @inline function update_gamma_candidate!(
     gamma_out,
     gamma_in,
diff --git a/test/panel/test_panel.jl b/test/panel/test_panel.jl
index 113b75f4..05c8b42a 100644
--- a/test/panel/test_panel.jl
+++ b/test/panel/test_panel.jl
@@ -1,8 +1,11 @@
 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, panel_inflow
+using VortexStepMethod: panel_axes, panel_span_vector, panel_inflow, panel_force_directions
 using Interpolations: linear_interpolation, Line
 using LinearAlgebra
 using Test
+if !@isdefined(test_data_path)
+    include("../test_data_utils.jl")
+end
 
 function create_panel(section1::Section, section2::Section)
     # Calculate panel geometry
@@ -217,14 +220,17 @@ end
     end
 
     @testset "alpha is measured in the plane square to the bound vortex" begin
-        chord_in_plane = normalize(axes.x_airf .-
-                                   dot(axes.x_airf, axes.y_airf) .* axes.y_airf)
         for v_eff in ([1.0, 0.0, 0.0], [1.0, 0.3, 0.1]), orient in (1, -1)
             oriented = panel_axes(le_1, te_1, le_2, te_2, 0.5, orient)
+            (; x_airf, y_airf, z_airf) = oriented
             alpha = panel_inflow(oriented, v_eff, v_eff, zeros(3)).alpha
-            @test alpha ≈ atan(dot(v_eff, oriented.z_airf), dot(v_eff, chord_in_plane))
-            along_span = v_eff .+ 2.0 .* axes.y_airf
+            @test alpha ≈ crossflow_alpha(v_eff, x_airf, y_airf, z_airf)
+            along_span = v_eff .+ 2.0 .* y_airf
             @test panel_inflow(oriented, along_span, along_span, zeros(3)).alpha ≈ alpha
+            dir_lift = panel_force_directions(oriented, alpha, y_airf).dir_lift
+            v_crossflow = v_eff .- dot(v_eff, y_airf) .* y_airf
+            @test abs(dot(dir_lift, v_crossflow)) < 1e-12
+            @test abs(dot(dir_lift, y_airf)) < 1e-12
         end
         @test rad2deg(panel_inflow(axes, [1.0, 0, 0], [1.0, 0, 0], zeros(3)).alpha) ≈
               -20.097 atol = 1e-3
diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl
index adcb1898..25ebf514 100644
--- a/test/solver/test_solver.jl
+++ b/test/solver/test_solver.jl
@@ -340,24 +340,21 @@ end
     va = [15.0, 1.0, 1.5]
     set_va!(body_aero, va)
 
-    function crossflow_alpha(panel, velocity)
-        chord_in_plane = panel.x_airf .- dot(panel.x_airf, panel.y_airf) .* panel.y_airf
-        return atan(dot(velocity, panel.z_airf), dot(velocity, normalize(chord_in_plane)))
-    end
-    induced(AIC, gamma, i) = [dot(AIC[i, :, k], gamma) for k in 1:3]
-
     for correct_aoa in (false, true)
         solver = Solver(length(body_aero.panels), 3; correct_aoa)
         sol = solve!(solver, body_aero)
         @test sol.solver_status == FEASIBLE
         gamma = solver.lr.gamma_new
         for (i, panel) in enumerate(body_aero.panels)
+            frame = (panel.x_airf, panel.y_airf, panel.z_airf)
+            induced = body_aero.AIC[i, :, :]' * gamma
             @test solver.lr.alpha_dist[i] ≈
-                  crossflow_alpha(panel, va .+ induced(body_aero.AIC, gamma, i)) atol = 1e-6
-            @test sol.alpha_geometric_dist[i] ≈ crossflow_alpha(panel, va) atol = 1e-10
+                  crossflow_alpha(va .+ induced, frame...) atol = 1e-6
+            @test sol.alpha_geometric_dist[i] ≈ crossflow_alpha(va, frame...) atol = 1e-10
             correct_aoa || continue
-            @test sol.alpha_dist[i] ≈ crossflow_alpha(panel,
-                va .+ induced(body_aero.AIC_aero_center, gamma, i)) atol = 1e-8
+            induced_aero_center = body_aero.AIC_aero_center[i, :, :]' * gamma
+            @test sol.alpha_dist[i] ≈
+                  crossflow_alpha(va .+ induced_aero_center, frame...) atol = 1e-8
         end
     end
 end
diff --git a/test/test_data_utils.jl b/test/test_data_utils.jl
index 00cdf7f6..4ef182b8 100644
--- a/test/test_data_utils.jl
+++ b/test/test_data_utils.jl
@@ -213,3 +213,14 @@ function get_complete_settings_file(module_name)
         return create_temp_wing_settings(module_name, basename(get_standard_wing_file(module_name)))
     end
 end
+
+"""
+    crossflow_alpha(velocity, x_airf, y_airf, z_airf)
+
+Angle of `velocity` in the plane square to `y_airf`, from the chord `x_airf` projected
+onto that plane towards `z_airf`: Gaunaa et al. 2026 Eqs. (1) and (3), formed directly.
+"""
+function crossflow_alpha(velocity, x_airf, y_airf, z_airf)
+    chord_in_plane = normalize(x_airf .- dot(x_airf, y_airf) .* y_airf)
+    return atan(dot(velocity, z_airf), dot(velocity, chord_in_plane))
+end