From fee98850edad266a8da8ab259ac6094e935c4d30 Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Wed, 23 Sep 2026 11:15:44 +0200
Subject: [PATCH 1/8] Default correct_aoa to true and exclude only the
 section's own bound vortex at the aero centre

Bart decided on #318 to default correct_aoa to true. The aero-centre
matrix used for the force directions now keeps every other panel's bound
filament (Gaunaa, Li & Pirrung 2026, Eq. 2) instead of zeroing all of them.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 .../318-correct-aoa-defaults-to-false-so-the-sol.md        | 3 +++
 data/TUDELFT_V3_KITE/vsm_settings.yaml                     | 2 +-
 data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml              | 4 +---
 src/body_aerodynamics.jl                                   | 7 ++++---
 src/settings.jl                                            | 6 +++---
 src/solver.jl                                              | 2 +-
 6 files changed, 13 insertions(+), 11 deletions(-)
 create mode 100644 changelog.d/318-correct-aoa-defaults-to-false-so-the-sol.md

diff --git a/changelog.d/318-correct-aoa-defaults-to-false-so-the-sol.md b/changelog.d/318-correct-aoa-defaults-to-false-so-the-sol.md
new file mode 100644
index 00000000..e5a13dd6
--- /dev/null
+++ b/changelog.d/318-correct-aoa-defaults-to-false-so-the-sol.md
@@ -0,0 +1,3 @@
+### Changed
+
+- `correct_aoa` defaults to `true`: the polars are still read at the 3/4-chord control point, and each section force is turned by the flow at its 1/4-chord aerodynamic centre (Gaunaa, Li & Pirrung, TORQUE 2026). Every VSM result moves, the induced drag most. That flow now leaves out only the section's own bound vortex, not every panel's, which also moves `LLT` results on swept or curved wings.
diff --git a/data/TUDELFT_V3_KITE/vsm_settings.yaml b/data/TUDELFT_V3_KITE/vsm_settings.yaml
index 5477aed9..5b599db0 100644
--- a/data/TUDELFT_V3_KITE/vsm_settings.yaml
+++ b/data/TUDELFT_V3_KITE/vsm_settings.yaml
@@ -53,7 +53,7 @@ solver_settings:
     # --- Core Aerodynamic Model ---
     aerodynamic_model_type: VSM       # VSM=3D vortex method, LLT=lifting line theory
     solver_type: LOOP                 # LOOP=fixed-point iteration, NEWTON=Newton-Raphson
-    #correct_aoa: false               # Defaults to false, meaning the computed AoA at 3/4c is used
+    #correct_aoa: true                # true=force directions from the flow at 1/4 chord, false=at 3/4 chord
     
     # --- Physical Properties ---
     density: 1.225                   # [kg/m³] Air density (ISA sea level)
diff --git a/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml b/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml
index d3fb98e8..3d689609 100644
--- a/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml
+++ b/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml
@@ -57,7 +57,7 @@ solver_settings:
   # --- Core Aerodynamic Model ---
   aerodynamic_model_type: VSM # VSM=3D vortex method, LLT=lifting line theory
   solver_type: LOOP # LOOP=fixed-point iteration, NONLIN=nonlinear solver
-  #correct_aoa: false               # Defaults to false, meaning the computed AoA at 3/4c is used
+  correct_aoa: true # true=force directions from the flow at 1/4 chord, false=at 3/4 chord
 
   # --- Physical Properties ---
   density: 1.225 # [kg/m³] Air density (ISA sea level)
@@ -78,5 +78,3 @@ solver_settings:
 
   # --- Output Control ---
   calc_only_f_and_gamma: false # true=compute only forces & circulation
-
-  correct_aoa: false
diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index 67a650c9..c203d5e2 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -429,7 +429,9 @@ end
     calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::Model, core_radius_fraction,
                             va_dist, va_unit_dist, target=body_aero.AIC)
 
-Calculate Aerodynamic Influence Coefficient matrices.
+Calculate the velocity each panel's unit-strength vortex ring induces at every control
+point (`VSM`, less the panel's own 2D bound vortex) or aerodynamic centre (`LLT`, less
+the panel's own bound filament).
 
 See also: [`BodyAerodynamics`](@ref), [`Model`](@ref)
 
@@ -442,7 +444,6 @@ Returns: nothing
                               target::AbstractArray{T, 3}=body_aero.AIC) where {P, W, T}
     # Determine evaluation point based on model
     evaluation_point = model == VSM ? :control_point : :aero_center
-    evaluation_point_on_bound = model == LLT
 
     # Allocate work vectors for this function (separate from those used by child functions)
     velocity_induced = zeros(MVector{3, T})
@@ -479,7 +480,7 @@ Returns: nothing
                 tempvel,
                 filaments,
                 ep,
-                evaluation_point_on_bound,
+                model == LLT && icp == jring,
                 va,
                 va_unit,
                 one(T),
diff --git a/src/settings.jl b/src/settings.jl
index 4062d56d..8e34492a 100644
--- a/src/settings.jl
+++ b/src/settings.jl
@@ -174,8 +174,8 @@ Solver configuration, used within [`VSMSettings`](@ref).
 - `calc_only_f_and_gamma`: Leave the analysis fields of [`VSMSolution`](@ref) (`lift`,
     `cl`, their distributions, span and centers of pressure) at their last values
     (default `false`)
-- `correct_aoa`: Perform angle of attack correction (VSM model only)
-    (default `false`)
+- `correct_aoa`: Turn each section force by the flow at its aerodynamic centre
+    rather than at its control point, VSM model only (default `true`)
 - `flow_curvature`: Add the thin-airfoil pitch-rate moment increment to each
     section (default `false`)
 - `is_with_viscous_drag_correction`: Add the spanwise-flow viscous drag and side
@@ -199,7 +199,7 @@ Solver configuration, used within [`VSMSettings`](@ref).
     core_radius_fraction::Float64 = 0.05
     mu::Float64 = 1.81e-5                   # dynamic viscosity [N·s/m²]
     calc_only_f_and_gamma::Bool=false       # skip the analysis fields of VSMSolution
-    correct_aoa::Bool=false                 # perform aoa correction
+    correct_aoa::Bool=true                  # force directions from the aero-centre flow
     flow_curvature::Bool=false              # thin-airfoil pitch-rate moment increment
     is_with_viscous_drag_correction::Bool=false # spanwise-flow viscous force (Gaunaa 2024)
     is_with_attached_trailed_force::Bool=false # chordwise bound-vortex force (Gaunaa 2026)
diff --git a/src/solver.jl b/src/solver.jl
index 2cbe9a77..dd747aa8 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -211,7 +211,7 @@ sol::VSMSolution = VSMSolution(): The result of calling [`solve!`](@ref)
     core_radius_fraction::T = T(0.05)
     mu::T = T(1.81e-5)
     is_only_f_and_gamma_output::Bool = false
-    correct_aoa::Bool = false
+    correct_aoa::Bool = true
     flow_curvature::Bool = false
     is_with_viscous_drag_correction::Bool = false
     is_with_attached_trailed_force::Bool = false

From 40e63d9446dfb97fb4c86e8ecc0ad99ee9b0a76b Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Wed, 30 Sep 2026 10:47:19 +0200
Subject: [PATCH 2/8] Test the elliptic AR 4 Oswald efficiency and the
 aero-centre matrix's bound filaments

Also list correct_aoa in the Solver docstring.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/solver.jl                                 |  2 +
 .../test_body_aerodynamics.jl                 | 45 ++++++++++++++++++-
 test/verification/test_verification.jl        | 11 +++++
 3 files changed, 57 insertions(+), 1 deletion(-)

diff --git a/src/solver.jl b/src/solver.jl
index dd747aa8..4fcdcba0 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -172,6 +172,8 @@ Main solver structure for the Vortex Step Method. See also: [`solve!`](@ref)
 - mu::Float64 = 1.81e-5: Dynamic viscosity [N·s/m²]
 - `is_only_f_and_gamma_output`::Bool = false: Whether `solve!` skips the fields of
   [`VSMSolution`](@ref) that only analysis reads, see [`SolverSettings`](@ref)
+- `correct_aoa`::Bool = true: Turn each section force by the flow at its aerodynamic
+    centre rather than at its control point, VSM model only
 - `flow_curvature`::Bool = false: Add the thin-airfoil pitch-rate moment
     increment `-(π/4) q̂` to each section, see: [`flow_curvature_cm`](@ref)
 - `is_with_viscous_drag_correction`::Bool = false: Add the spanwise-flow viscous drag
diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl
index a47650a6..e1ec5579 100644
--- a/test/body_aerodynamics/test_body_aerodynamics.jl
+++ b/test/body_aerodynamics/test_body_aerodynamics.jl
@@ -1,5 +1,7 @@
 using VortexStepMethod
-using VortexStepMethod: calculate_cl, calculate_cd_cm, calculate_projected_area, calculate_AIC_matrices!
+using VortexStepMethod: calculate_cl, calculate_cd_cm, calculate_projected_area, calculate_AIC_matrices!,
+    velocity_3D_bound_vortex!, velocity_3D_trailing_vortex!,
+    velocity_3D_trailing_vortex_semiinfinite!
 using LinearAlgebra
 using Test
 using Logging
@@ -133,6 +135,47 @@ end
     end
 end
 
+@testset "LLT matrix leaves out only each panel's own bound filament" begin
+    radius = 4.0
+    wing = Wing(3)
+    for theta in range(-π / 4, π / 4, length=4)
+        section_point = [0.0, radius * sin(theta), radius * cos(theta)]
+        add_section!(wing, section_point, section_point .+ [1.0, 0.0, 0.0], INVISCID)
+    end
+    refine!(wing)
+    body_aero = BodyAerodynamics([wing])
+    va_vec = 10.0 .* [cosd(5), 0.0, sind(5)]
+    set_va!(body_aero, va_vec)
+    n_panels = length(body_aero.panels)
+    core_radius_fraction = 0.05
+    calculate_AIC_matrices!(body_aero, LLT, core_radius_fraction,
+        fill(norm(va_vec), n_panels), repeat(normalize(va_vec)', n_panels))
+
+    work_vectors = body_aero.work_vectors
+    velocity = zeros(3)
+    for (i, panel_i) in enumerate(body_aero.panels), (j, panel_j) in enumerate(body_aero.panels)
+        filaments = panel_j.filaments
+        expected = zeros(3)
+        for filament in filaments[2:3]
+            velocity_3D_trailing_vortex!(velocity, filament, panel_i.aero_center, 1.0,
+                norm(va_vec), work_vectors)
+            expected .+= velocity
+        end
+        for filament in filaments[4:5]
+            velocity_3D_trailing_vortex_semiinfinite!(velocity, filament, normalize(va_vec),
+                panel_i.aero_center, 1.0, norm(va_vec), work_vectors)
+            expected .+= velocity
+        end
+        if i != j
+            velocity_3D_bound_vortex!(velocity, filaments[1], panel_i.aero_center, 1.0,
+                core_radius_fraction, work_vectors)
+            @test norm(velocity) > 1e-3
+            expected .+= velocity
+        end
+        @test body_aero.AIC[i, j, :] ≈ expected atol = 1e-12
+    end
+end
+
 
 @testset "Wing Geometry Creation" begin
     @testset "Origin Translation" begin
diff --git a/test/verification/test_verification.jl b/test/verification/test_verification.jl
index cdf39467..7d51a3d0 100644
--- a/test/verification/test_verification.jl
+++ b/test/verification/test_verification.jl
@@ -105,6 +105,17 @@ max_error(actual, expected) = maximum(abs.(actual .- expected))
         @test max_error(CD_llt, CD_vsm) < 1e-2
     end
 
+    @testset "elliptic AR 4 wing reaches Oswald efficiency 1 (Gaunaa et al. 2026)" begin
+        aspect_ratio = 4.0
+        span = aspect_ratio * π / 4
+        coordinates = generate_coordinates_el_wing(1.0, span, 41, "cos")
+        wing = wing_from_coordinates(coordinates, INVISCID)
+        CL, CD = lift_drag_polar(wing, VSM, [10.0]; va=10.0, relaxation_factor=0.05)
+        oswald_efficiency = only(CL)^2 / (π * aspect_ratio * only(CD))
+
+        @test oswald_efficiency ≈ 1 atol = 0.01
+    end
+
     @testset "curved Clark Y wing matches RANS" begin
         alphas = [9.0]
         coordinates = generate_coordinates_curved_wing(2.18, 6.969, π / 4, 4.673, 60,

From 92aed6b012f7ba9c42361fbd5526c37b429dac64 Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Wed, 30 Sep 2026 10:49:48 +0200
Subject: [PATCH 3/8] Say what flow update_effective_angle_of_attack! measures

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/body_aerodynamics.jl                         | 8 +++-----
 test/body_aerodynamics/test_body_aerodynamics.jl | 7 ++++---
 2 files changed, 7 insertions(+), 8 deletions(-)

diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index c203d5e2..e72688a5 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -442,7 +442,6 @@ Returns: nothing
                               va_dist::AbstractVector{T},
                               va_unit_dist::AbstractMatrix{T},
                               target::AbstractArray{T, 3}=body_aero.AIC) where {P, W, T}
-    # Determine evaluation point based on model
     evaluation_point = model == VSM ? :control_point : :aero_center
 
     # Allocate work vectors for this function (separate from those used by child functions)
@@ -527,10 +526,9 @@ end
                                       core_radius_fraction, z_airf_dist, x_airf_dist,
                                       va_vec_dist, va_dist, va_unit_dist)
 
-Update angle of attack at aerodynamic center for VSM method.
-
-Returns:
-    nothing
+Write into `alpha_corrected` [rad] the angle of attack at each panel's aerodynamic centre,
+from the inflow plus the velocity `gamma` induces there without the panel's own bound
+filament.
 """
 function update_effective_angle_of_attack!(alpha_corrected,
     body_aero::BodyAerodynamics, 
diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl
index e1ec5579..4e87a56c 100644
--- a/test/body_aerodynamics/test_body_aerodynamics.jl
+++ b/test/body_aerodynamics/test_body_aerodynamics.jl
@@ -1,6 +1,6 @@
 using VortexStepMethod
-using VortexStepMethod: calculate_cl, calculate_cd_cm, calculate_projected_area, calculate_AIC_matrices!,
-    velocity_3D_bound_vortex!, velocity_3D_trailing_vortex!,
+using VortexStepMethod: calculate_cl, calculate_cd_cm, calculate_projected_area,
+    calculate_AIC_matrices!, velocity_3D_bound_vortex!, velocity_3D_trailing_vortex!,
     velocity_3D_trailing_vortex_semiinfinite!
 using LinearAlgebra
 using Test
@@ -153,7 +153,8 @@ end
 
     work_vectors = body_aero.work_vectors
     velocity = zeros(3)
-    for (i, panel_i) in enumerate(body_aero.panels), (j, panel_j) in enumerate(body_aero.panels)
+    panels = body_aero.panels
+    for (i, panel_i) in enumerate(panels), (j, panel_j) in enumerate(panels)
         filaments = panel_j.filaments
         expected = zeros(3)
         for filament in filaments[2:3]

From c5b8229025c2564e1574e556a153e9c02e652144 Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Wed, 30 Sep 2026 10:55:09 +0200
Subject: [PATCH 4/8] Keep calculate_AIC_matrices! allocation-free and a vortex
 segment finite on its own axis

The correct_aoa default now calls calculate_AIC_matrices! from calc_forces!
on every solve, so its temporaries come from the body's buffers. With a
near-zero core, a point on a filament's line fell through the colinearity
test on rounding noise and came back NaN.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/body_aerodynamics.jl                        | 15 ++++++++-------
 src/filament.jl                                 |  9 +++++----
 src/solver.jl                                   |  2 +-
 .../body_aerodynamics/test_body_aerodynamics.jl | 17 +++++++++++++++++
 4 files changed, 31 insertions(+), 12 deletions(-)

diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index e72688a5..780abf71 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -444,18 +444,19 @@ Returns: nothing
                               target::AbstractArray{T, 3}=body_aero.AIC) where {P, W, T}
     evaluation_point = model == VSM ? :control_point : :aero_center
 
-    # Allocate work vectors for this function (separate from those used by child functions)
     velocity_induced = zeros(MVector{3, T})
-    tempvel = zeros(MVector{3, T})
-    va_unit = zeros(MVector{3, T})
     U_2D = zeros(MVector{3, T})
+    # Slots 1-7 are the filament kernels' scratch; these two escape into them.
+    tempvel = body_aero.work_vectors[8]
+    va_unit = body_aero.work_vectors[9]
 
     # Python parity: one shared area-weighted wake vector for all panels.
-    panel_areas = [panel.chord * panel.width for panel in body_aero.panels]
-    va_vec_dist = zeros(T, length(body_aero.panels), 3)
-    @inbounds for i in 1:length(body_aero.panels), k in 1:3
-        va_vec_dist[i, k] = va_unit_dist[i, k] * va_dist[i]
+    panel_areas = body_aero.cache[5][va_dist]
+    for (i, panel) in enumerate(body_aero.panels)
+        panel_areas[i] = panel.chord * panel.width
     end
+    va_vec_dist = body_aero.cache[6][va_unit_dist]
+    va_vec_dist .= va_unit_dist .* va_dist
     wake_velocity = _compute_reference_velocity_from_distribution(
         va_vec_dist,
         length(body_aero.panels),
diff --git a/src/filament.jl b/src/filament.jl
index 023addec..68c6e8be 100644
--- a/src/filament.jl
+++ b/src/filament.jl
@@ -94,7 +94,7 @@ end
 
 Calculate the Biot–Savart velocity induced by a straight vortex segment at `XVP`.
 Inside the core radius `epsilon` the velocity is evaluated on the core boundary and
-scaled linearly with the distance to the axis.
+scaled linearly with the distance to the axis; on the axis, to rounding, it is zero.
 """
 @inline function velocity_3D_vortex_segment!(
     vel,
@@ -112,7 +112,10 @@ scaled linearly with the distance to the axis.
 
     cross3!(r1Xr0, r1, r0)
     nr1Xr0 = norm3(r1Xr0)
-    if nr1Xr0 / nr0 > epsilon
+    axis_distance = nr1Xr0 / nr0
+    if axis_distance < 1e-12 * max(epsilon, norm3(r1))
+        vel .= 0.0
+    elseif axis_distance > epsilon
         cross3!(r1Xr2, r1, r2)
         nr1 = norm3(r1)
         nr2 = norm3(r2)
@@ -124,8 +127,6 @@ scaled linearly with the distance to the axis.
         @inbounds for k in 1:3
             vel[k] = coeff * r1Xr2[k]
         end
-    elseif nr1Xr0 / nr0 < 1e-12 * epsilon
-        vel .= 0.0
     else
         nr0sq = nr0 * nr0
         d_r1_r0 = dot3(r1, r0)
diff --git a/src/solver.jl b/src/solver.jl
index 4fcdcba0..fd44ad57 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -453,7 +453,7 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
     panel_moment_dist = solver.sol.panel_moment_dist
 
     # Calculate alpha corrections based on model type
-    if solver.correct_aoa && aerodynamic_model_type == VSM      # 64 bytes
+    if solver.correct_aoa && aerodynamic_model_type == VSM
         update_effective_angle_of_attack!(
             alpha_corrected,
             body_aero,
diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl
index 4e87a56c..b21ca2a0 100644
--- a/test/body_aerodynamics/test_body_aerodynamics.jl
+++ b/test/body_aerodynamics/test_body_aerodynamics.jl
@@ -177,6 +177,23 @@ end
     end
 end
 
+@testset "LLT matrix is finite on a swept wing without vortex core" begin
+    wing = Wing(8)
+    for y in range(-2.0, 2.0, length=9)
+        leading_edge = [0.2abs(y), y, -0.5abs(y)]
+        add_section!(wing, leading_edge, leading_edge .+ [1.0, 0.0, 0.0], INVISCID)
+    end
+    refine!(wing)
+    body_aero = BodyAerodynamics([wing])
+    va_vec = 10.0 .* [cosd(5), 0.0, sind(5)]
+    set_va!(body_aero, va_vec)
+    n_panels = length(body_aero.panels)
+    calculate_AIC_matrices!(body_aero, LLT, 1e-20, fill(norm(va_vec), n_panels),
+        repeat(normalize(va_vec)', n_panels))
+
+    @test all(isfinite, body_aero.AIC)
+end
+
 
 @testset "Wing Geometry Creation" begin
     @testset "Origin Translation" begin

From 98775b116de2f123176cab3572be818ecc6eae18 Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Wed, 30 Sep 2026 11:00:20 +0200
Subject: [PATCH 5/8] Evaluate the vortex core in closed form so it stays
 differentiable on the axis

With the aero-centre matrix keeping neighbouring bound filaments, an aero
centre on a refined straight segment lies on its neighbours' axes. The core
branch normalised the radial vector there and the old guard returned a hard
zero, so ForwardDiff lost the slope FiniteDiff sees. The closed form is the
same velocity (2e-11 relative) and linear through the axis.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/filament.jl                      | 37 ++++++++--------------------
 test/filament/test_bound_filament.jl | 17 +++++++++++++
 2 files changed, 27 insertions(+), 27 deletions(-)

diff --git a/src/filament.jl b/src/filament.jl
index 68c6e8be..04ebc6af 100644
--- a/src/filament.jl
+++ b/src/filament.jl
@@ -94,7 +94,7 @@ end
 
 Calculate the Biot–Savart velocity induced by a straight vortex segment at `XVP`.
 Inside the core radius `epsilon` the velocity is evaluated on the core boundary and
-scaled linearly with the distance to the axis; on the axis, to rounding, it is zero.
+scaled linearly with the distance to the axis. Without a core, it is zero on the axis.
 """
 @inline function velocity_3D_vortex_segment!(
     vel,
@@ -104,16 +104,16 @@ scaled linearly with the distance to the axis; on the axis, to rounding, it is z
     epsilon,
     work_vectors
 )
-    r1, r2, r1Xr2, r1Xr0, r1r2norm, r1_proj, r2_proj = work_vectors
+    r1, r2, r1Xr2, r1Xr0, r1r2norm = work_vectors
     r0 = filament.r0
     nr0 = filament.length
     r1 .= XVP .- filament.x1
     r2 .= XVP .- filament.x2
 
     cross3!(r1Xr0, r1, r0)
-    nr1Xr0 = norm3(r1Xr0)
-    axis_distance = nr1Xr0 / nr0
-    if axis_distance < 1e-12 * max(epsilon, norm3(r1))
+    axis_distance = norm3(r1Xr0) / nr0
+    rounding = 1e-12 * norm3(r1)
+    if epsilon <= rounding && axis_distance <= rounding
         vel .= 0.0
     elseif axis_distance > epsilon
         cross3!(r1Xr2, r1, r2)
@@ -128,32 +128,15 @@ scaled linearly with the distance to the axis; on the axis, to rounding, it is z
             vel[k] = coeff * r1Xr2[k]
         end
     else
+        # The core-boundary velocity scaled by axis_distance / epsilon, in closed form.
         nr0sq = nr0 * nr0
         d_r1_r0 = dot3(r1, r0)
         d_r2_r0 = dot3(r2, r0)
-        r_rad = r1Xr0
+        d_sum = d_r1_r0 / sqrt(d_r1_r0^2 / nr0sq + epsilon^2) -
+                d_r2_r0 / sqrt(d_r2_r0^2 / nr0sq + epsilon^2)
+        coeff = -(gamma / (4π)) * d_sum / (epsilon^2 * nr0sq)
         @inbounds for k in 1:3
-            r_rad[k] = r1[k] - d_r1_r0 * r0[k] / nr0sq
-        end
-        nr_rad = norm3(r_rad)
-        @inbounds for k in 1:3
-            r1_proj[k] = d_r1_r0 * r0[k] / nr0sq +
-                         epsilon * r_rad[k] / nr_rad
-            r2_proj[k] = d_r2_r0 * r0[k] / nr0sq +
-                         epsilon * r_rad[k] / nr_rad
-        end
-        cross3!(r1Xr2, r1_proj, r2_proj)
-        nr1_proj = norm3(r1_proj)
-        nr2_proj = norm3(r2_proj)
-        d_sum = 0.0
-        @inbounds for k in 1:3
-            d_sum += r0[k] * (r1_proj[k]/nr1_proj -
-                              r2_proj[k]/nr2_proj)
-        end
-        scale = nr1Xr0 / (nr0 * epsilon)
-        coeff = scale * (gamma / (4π)) / (norm3(r1Xr2)^2) * d_sum
-        @inbounds for k in 1:3
-            vel[k] = coeff * r1Xr2[k]
+            vel[k] = coeff * r1Xr0[k]
         end
     end
     nothing
diff --git a/test/filament/test_bound_filament.jl b/test/filament/test_bound_filament.jl
index 397eaf99..694ed86e 100644
--- a/test/filament/test_bound_filament.jl
+++ b/test/filament/test_bound_filament.jl
@@ -1,5 +1,6 @@
 using VortexStepMethod: BoundFilament, velocity_3D_bound_vortex!,
     velocity_3D_trailing_vortex!, reinit!, ALPHA0, NU
+using ForwardDiff
 using LinearAlgebra
 using Test
 
@@ -260,4 +261,20 @@ end
             @test v[3] > 0
         end
     end
+
+    @testset "ForwardDiff sees the core's slope on the extended axis" begin
+        function off_axis_velocity(offset)
+            T = typeof(offset)
+            filament = BoundFilament{T}()
+            reinit!(filament, zeros(T, 3), T[1, 0, 0])
+            velocity = zeros(T, 3)
+            velocity_3D_bound_vortex!(velocity, filament, [1.5, offset, 0.0], gamma,
+                core_radius_fraction, ntuple(_ -> zeros(T, 3), 10))
+            return velocity[3]
+        end
+        slope = ForwardDiff.derivative(off_axis_velocity, 0.0)
+
+        @test slope != 0
+        @test slope ≈ off_axis_velocity(1e-4) / 1e-4 rtol = 1e-9
+    end
 end

From 800c244abcb138949ff1ac7a68d04e2587ba0cf7 Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Wed, 30 Sep 2026 11:02:53 +0200
Subject: [PATCH 6/8] Name the angle the polars are read at, and bisect trims
 below the CMy bound

cl_dist is read at alpha_uncorrected, which equals alpha_dist only without
correct_aoa. trim_angle's default alpha_tol bounds |CMy| by 3.3e-5 on the
trimmable wing, so the 1e-5 check passed by where the bisection landed.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 test/solver/test_polar_vectors.jl | 2 +-
 test/solver/test_stability.jl     | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/solver/test_polar_vectors.jl b/test/solver/test_polar_vectors.jl
index 48439163..2544223a 100644
--- a/test/solver/test_polar_vectors.jl
+++ b/test/solver/test_polar_vectors.jl
@@ -38,7 +38,7 @@ end
         sol = solve_three_panel_wing(POLAR_VECTORS,
             linear_polar(4π, 0.1, -0.02), linear_polar(2π, 0.05, -0.01))
         @test sol.solver_status == FEASIBLE
-        @test sol.cl_dist ≈ [4π, 3π, 2π] .* sol.alpha_dist
+        @test sol.cl_dist ≈ [4π, 3π, 2π] .* sol.alpha_uncorrected
         @test sol.cd_dist ≈ [0.1, 0.075, 0.05]
         @test sol.cm_dist ≈ [-0.02, -0.015, -0.01]
         @testset "the stronger +y side rolls the wing positive about x" begin
diff --git a/test/solver/test_stability.jl b/test/solver/test_stability.jl
index 758b68a1..180dcec8 100644
--- a/test/solver/test_stability.jl
+++ b/test/solver/test_stability.jl
@@ -81,7 +81,7 @@ end
 
     @testset "moments about the leading edge: stable trim" begin
         body_aero, solver = trimmable_wing(0.05)
-        trims = trim_angle(solver, body_aero, beta, va)
+        trims = trim_angle(solver, body_aero, beta, va; alpha_tol=1e-7)
         @test length(trims) == 1
         trim = only(trims)
         trim_coeffs = coeffs_at_angles(solver, body_aero, trim.alpha, beta, va)
@@ -93,7 +93,7 @@ end
 
     @testset "moments about the trailing edge: unstable trim" begin
         body_aero, solver = trimmable_wing(-0.05; reference_point=[1.0, 0.0, 0.0])
-        trim = only(trim_angle(solver, body_aero, beta, va))
+        trim = only(trim_angle(solver, body_aero, beta, va; alpha_tol=1e-7))
         trim_coeffs = coeffs_at_angles(solver, body_aero, trim.alpha, beta, va)
         @test abs(trim_coeffs[5]) < 1e-5
         @test trim.dCMy_dalpha > 0

From 7c904f0271b2161f0e2d257a68cb330a7423b01b Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Wed, 30 Sep 2026 11:06:46 +0200
Subject: [PATCH 7/8] Tidy the changed functions and tests after review

Drop comments that restate code, name the kernels' work-vector slots
correctly, share the uniform-inflow AIC call across the matrix tests and
hoist the ForwardDiff probe out of its testset.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 ...orrect-aoa-defaults-to-false-so-the-sol.md |  4 ++
 data/TUDELFT_V3_KITE/vsm_settings.yaml        |  2 +-
 data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml |  2 +-
 src/body_aerodynamics.jl                      | 14 ++----
 src/filament.jl                               |  1 -
 src/solver.jl                                 |  3 +-
 .../test_body_aerodynamics.jl                 | 44 +++++++------------
 test/filament/test_bound_filament.jl          | 30 ++++++++-----
 8 files changed, 45 insertions(+), 55 deletions(-)

diff --git a/changelog.d/318-correct-aoa-defaults-to-false-so-the-sol.md b/changelog.d/318-correct-aoa-defaults-to-false-so-the-sol.md
index e5a13dd6..19470cb6 100644
--- a/changelog.d/318-correct-aoa-defaults-to-false-so-the-sol.md
+++ b/changelog.d/318-correct-aoa-defaults-to-false-so-the-sol.md
@@ -1,3 +1,7 @@
 ### Changed
 
 - `correct_aoa` defaults to `true`: the polars are still read at the 3/4-chord control point, and each section force is turned by the flow at its 1/4-chord aerodynamic centre (Gaunaa, Li & Pirrung, TORQUE 2026). Every VSM result moves, the induced drag most. That flow now leaves out only the section's own bound vortex, not every panel's, which also moves `LLT` results on swept or curved wings.
+
+### Fixed
+
+- A vortex segment evaluated on its own axis returns zero rather than NaN when `core_radius_fraction` is near zero, and ForwardDiff sees the velocity's slope there rather than zero.
diff --git a/data/TUDELFT_V3_KITE/vsm_settings.yaml b/data/TUDELFT_V3_KITE/vsm_settings.yaml
index 5b599db0..9326b2bf 100644
--- a/data/TUDELFT_V3_KITE/vsm_settings.yaml
+++ b/data/TUDELFT_V3_KITE/vsm_settings.yaml
@@ -53,7 +53,7 @@ solver_settings:
     # --- Core Aerodynamic Model ---
     aerodynamic_model_type: VSM       # VSM=3D vortex method, LLT=lifting line theory
     solver_type: LOOP                 # LOOP=fixed-point iteration, NEWTON=Newton-Raphson
-    #correct_aoa: true                # true=force directions from the flow at 1/4 chord, false=at 3/4 chord
+    #correct_aoa: true                # true=forces turned by the 1/4-chord flow
     
     # --- Physical Properties ---
     density: 1.225                   # [kg/m³] Air density (ISA sea level)
diff --git a/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml b/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml
index 3d689609..9c44a7a9 100644
--- a/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml
+++ b/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml
@@ -57,7 +57,7 @@ solver_settings:
   # --- Core Aerodynamic Model ---
   aerodynamic_model_type: VSM # VSM=3D vortex method, LLT=lifting line theory
   solver_type: LOOP # LOOP=fixed-point iteration, NONLIN=nonlinear solver
-  correct_aoa: true # true=force directions from the flow at 1/4 chord, false=at 3/4 chord
+  correct_aoa: true # true=forces turned by the 1/4-chord flow
 
   # --- Physical Properties ---
   density: 1.225 # [kg/m³] Air density (ISA sea level)
diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index 780abf71..bb3b57e3 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -442,11 +442,9 @@ Returns: nothing
                               va_dist::AbstractVector{T},
                               va_unit_dist::AbstractMatrix{T},
                               target::AbstractArray{T, 3}=body_aero.AIC) where {P, W, T}
-    evaluation_point = model == VSM ? :control_point : :aero_center
-
     velocity_induced = zeros(MVector{3, T})
     U_2D = zeros(MVector{3, T})
-    # Slots 1-7 are the filament kernels' scratch; these two escape into them.
+    # Slots 1-5 are the filament kernels' scratch; these two escape into them.
     tempvel = body_aero.work_vectors[8]
     va_unit = body_aero.work_vectors[9]
 
@@ -473,8 +471,7 @@ Returns: nothing
         filaments = panel_jring.filaments
         for icp in eachindex(body_aero.panels)
             panel_icp = body_aero.panels[icp]
-            ep = evaluation_point == :control_point ? panel_icp.control_point :
-                 panel_icp.aero_center
+            ep = model == VSM ? panel_icp.control_point : panel_icp.aero_center
             calculate_velocity_induced_single_ring_semiinfinite!(
                 velocity_induced,
                 tempvel,
@@ -488,7 +485,6 @@ Returns: nothing
                 body_aero.work_vectors
             )
                       
-            # Subtract 2D induced velocity for VSM
             if icp == jring && model == VSM
                 calculate_velocity_induced_bound_2D!(U_2D, panel_jring, ep, body_aero.work_vectors)
                 velocity_induced .-= U_2D
@@ -541,8 +537,7 @@ function update_effective_angle_of_attack!(alpha_corrected,
     va_dist,
     va_unit_dist)
 
-    # Its own buffer: `AIC` holds the control-point matrix the circulation was solved
-    # against, so overwriting it here would leave post-solve readers on the LLT one.
+    # Not `AIC`: that holds the control-point matrix post-solve readers expect.
     calculate_AIC_matrices!(body_aero, LLT, core_radius_fraction, va_dist,
                             va_unit_dist, body_aero.AIC_aero_center)
 
@@ -551,11 +546,9 @@ 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]
@@ -571,7 +564,6 @@ function update_effective_angle_of_attack!(alpha_corrected,
         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])
     end
diff --git a/src/filament.jl b/src/filament.jl
index 04ebc6af..aeb22c16 100644
--- a/src/filament.jl
+++ b/src/filament.jl
@@ -128,7 +128,6 @@ scaled linearly with the distance to the axis. Without a core, it is zero on the
             vel[k] = coeff * r1Xr2[k]
         end
     else
-        # The core-boundary velocity scaled by axis_distance / epsilon, in closed form.
         nr0sq = nr0 * nr0
         d_r1_r0 = dot3(r1, r0)
         d_r2_r0 = dot3(r2, r0)
diff --git a/src/solver.jl b/src/solver.jl
index fd44ad57..fcd39a4e 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -407,7 +407,7 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
     aerodynamic_model_type = solver.aerodynamic_model_type
 
     # Calculate coefficients for each panel
-    for (i, panel) in enumerate(panels)                                               # zero bytes
+    for (i, panel) in enumerate(panels)
         cl_dist[i] = calculate_cl(panel, alpha_dist[i])
         cd_dist[i], cm_dist[i] = calculate_cd_cm(panel, alpha_dist[i])
         if solver.flow_curvature
@@ -452,7 +452,6 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
     drag = solver.sol.drag_dist
     panel_moment_dist = solver.sol.panel_moment_dist
 
-    # Calculate alpha corrections based on model type
     if solver.correct_aoa && aerodynamic_model_type == VSM
         update_effective_angle_of_attack!(
             alpha_corrected,
diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl
index b21ca2a0..95121071 100644
--- a/test/body_aerodynamics/test_body_aerodynamics.jl
+++ b/test/body_aerodynamics/test_body_aerodynamics.jl
@@ -26,6 +26,18 @@ function inviscid_wing(section_y; n_panels=length(section_y) - 1)
     return wing
 end
 
+"""
+    uniform_AIC!(body_aero, model, core_radius_fraction, va_vec)
+
+Write into `body_aero.AIC` the influence matrices of `model` under the uniform inflow
+`va_vec` [m/s].
+"""
+function uniform_AIC!(body_aero, model, core_radius_fraction, va_vec)
+    n_panels = length(body_aero.panels)
+    calculate_AIC_matrices!(body_aero, model, core_radius_fraction,
+        fill(norm(va_vec), n_panels), repeat(normalize(va_vec)', n_panels))
+end
+
 @testset "Induction Matrix Creation" begin
     # Setup
     n_panels = 3
@@ -79,17 +91,7 @@ end
             LLT
         )
 
-        # Calculate new matrices
-        va_dist = fill(norm(va_vec), length(body_aero.panels))
-        va_unit_dist = repeat(reshape(va_vec ./ norm(va_vec), 1, 3),
-                              length(body_aero.panels))
-        calculate_AIC_matrices!(
-            body_aero,
-            LLT,
-            core_radius_fraction,
-            va_dist,
-            va_unit_dist
-        )
+        uniform_AIC!(body_aero, LLT, core_radius_fraction, va_vec)
         AIC_x, AIC_y, AIC_z = @views body_aero.AIC[:, :, 1], body_aero.AIC[:, :, 2], body_aero.AIC[:, :, 3]
 
         # Compare matrices
@@ -115,17 +117,7 @@ end
             VSM
         )
 
-        # Calculate new matrices
-        va_dist = fill(norm(va_vec), length(body_aero.panels))
-        va_unit_dist = repeat(reshape(va_vec ./ norm(va_vec), 1, 3),
-                              length(body_aero.panels))
-        calculate_AIC_matrices!(
-            body_aero,
-            VSM,
-            core_radius_fraction,
-            va_dist,
-            va_unit_dist
-        )
+        uniform_AIC!(body_aero, VSM, core_radius_fraction, va_vec)
         AIC_x, AIC_y, AIC_z = body_aero.AIC[:, :, 1], body_aero.AIC[:, :, 2], body_aero.AIC[:, :, 3]
 
         # Compare matrices with higher precision for VSM
@@ -146,10 +138,8 @@ end
     body_aero = BodyAerodynamics([wing])
     va_vec = 10.0 .* [cosd(5), 0.0, sind(5)]
     set_va!(body_aero, va_vec)
-    n_panels = length(body_aero.panels)
     core_radius_fraction = 0.05
-    calculate_AIC_matrices!(body_aero, LLT, core_radius_fraction,
-        fill(norm(va_vec), n_panels), repeat(normalize(va_vec)', n_panels))
+    uniform_AIC!(body_aero, LLT, core_radius_fraction, va_vec)
 
     work_vectors = body_aero.work_vectors
     velocity = zeros(3)
@@ -187,9 +177,7 @@ end
     body_aero = BodyAerodynamics([wing])
     va_vec = 10.0 .* [cosd(5), 0.0, sind(5)]
     set_va!(body_aero, va_vec)
-    n_panels = length(body_aero.panels)
-    calculate_AIC_matrices!(body_aero, LLT, 1e-20, fill(norm(va_vec), n_panels),
-        repeat(normalize(va_vec)', n_panels))
+    uniform_AIC!(body_aero, LLT, 1e-20, va_vec)
 
     @test all(isfinite, body_aero.AIC)
 end
diff --git a/test/filament/test_bound_filament.jl b/test/filament/test_bound_filament.jl
index 694ed86e..e786d247 100644
--- a/test/filament/test_bound_filament.jl
+++ b/test/filament/test_bound_filament.jl
@@ -27,6 +27,22 @@ function analytical_solution(control_point, gamma)
            dot(r0, r1/norm(r1) - r2/norm(r2))
 end
 
+"""
+    off_axis_velocity(offset, gamma, core_radius_fraction)
+
+z velocity [m/s] induced by the unit filament from the origin along x at the point `offset`
+[m] off its axis, half a length beyond its end.
+"""
+function off_axis_velocity(offset, gamma, core_radius_fraction)
+    T = typeof(offset)
+    filament = BoundFilament{T}()
+    reinit!(filament, zeros(T, 3), T[1, 0, 0])
+    velocity = zeros(T, 3)
+    velocity_3D_bound_vortex!(velocity, filament, [1.5, offset, 0.0], gamma,
+        core_radius_fraction, ntuple(_ -> zeros(T, 3), 10))
+    return velocity[3]
+end
+
 @testset "BoundFilament Tests" begin
     gamma = 1.0
     core_radius_fraction = 0.01
@@ -263,18 +279,10 @@ end
     end
 
     @testset "ForwardDiff sees the core's slope on the extended axis" begin
-        function off_axis_velocity(offset)
-            T = typeof(offset)
-            filament = BoundFilament{T}()
-            reinit!(filament, zeros(T, 3), T[1, 0, 0])
-            velocity = zeros(T, 3)
-            velocity_3D_bound_vortex!(velocity, filament, [1.5, offset, 0.0], gamma,
-                core_radius_fraction, ntuple(_ -> zeros(T, 3), 10))
-            return velocity[3]
-        end
-        slope = ForwardDiff.derivative(off_axis_velocity, 0.0)
+        velocity_at(offset) = off_axis_velocity(offset, gamma, core_radius_fraction)
+        slope = ForwardDiff.derivative(velocity_at, 0.0)
 
         @test slope != 0
-        @test slope ≈ off_axis_velocity(1e-4) / 1e-4 rtol = 1e-9
+        @test slope ≈ velocity_at(1e-4) / 1e-4 rtol = 1e-9
     end
 end

From a6362600b189a4c675de65a003394bece304d07a Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Wed, 30 Sep 2026 11:56:49 +0200
Subject: [PATCH 8/8] Keep the AIC scratch in work_vectors and pin
 correct_aoa=false and the zero-core fallback

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/body_aerodynamics.jl                         |  6 +++---
 src/filament.jl                                  |  3 ++-
 test/body_aerodynamics/test_body_aerodynamics.jl |  3 +++
 test/verification/test_verification.jl           | 14 ++++++++------
 4 files changed, 16 insertions(+), 10 deletions(-)

diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index bb3b57e3..ef504700 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -442,9 +442,9 @@ Returns: nothing
                               va_dist::AbstractVector{T},
                               va_unit_dist::AbstractMatrix{T},
                               target::AbstractArray{T, 3}=body_aero.AIC) where {P, W, T}
-    velocity_induced = zeros(MVector{3, T})
-    U_2D = zeros(MVector{3, T})
-    # Slots 1-5 are the filament kernels' scratch; these two escape into them.
+    # Slots 1-5 are the filament kernels' scratch.
+    velocity_induced = body_aero.work_vectors[6]
+    U_2D = body_aero.work_vectors[7]
     tempvel = body_aero.work_vectors[8]
     va_unit = body_aero.work_vectors[9]
 
diff --git a/src/filament.jl b/src/filament.jl
index aeb22c16..d79f612d 100644
--- a/src/filament.jl
+++ b/src/filament.jl
@@ -94,7 +94,8 @@ end
 
 Calculate the Biot–Savart velocity induced by a straight vortex segment at `XVP`.
 Inside the core radius `epsilon` the velocity is evaluated on the core boundary and
-scaled linearly with the distance to the axis. Without a core, it is zero on the axis.
+scaled linearly with the distance to the axis. Without a core, it is zero within
+1e-12 of the axis, relative to the distance from `x1`.
 """
 @inline function velocity_3D_vortex_segment!(
     vel,
diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl
index 95121071..2cd8d749 100644
--- a/test/body_aerodynamics/test_body_aerodynamics.jl
+++ b/test/body_aerodynamics/test_body_aerodynamics.jl
@@ -177,9 +177,12 @@ end
     body_aero = BodyAerodynamics([wing])
     va_vec = 10.0 .* [cosd(5), 0.0, sind(5)]
     set_va!(body_aero, va_vec)
+    uniform_AIC!(body_aero, LLT, 1e-10, va_vec)
+    AIC_small_core = copy(body_aero.AIC)
     uniform_AIC!(body_aero, LLT, 1e-20, va_vec)
 
     @test all(isfinite, body_aero.AIC)
+    @test body_aero.AIC ≈ AIC_small_core
 end
 
 
diff --git a/test/verification/test_verification.jl b/test/verification/test_verification.jl
index 7d51a3d0..7ca63f41 100644
--- a/test/verification/test_verification.jl
+++ b/test/verification/test_verification.jl
@@ -49,15 +49,15 @@ function wing_from_coordinates(coordinates, aero_model, aero_data=nothing;
 end
 
 """
-    lift_drag_polar(wing, model, alphas; va, relaxation_factor)
+    lift_drag_polar(wing, model, alphas; va, relaxation_factor, correct_aoa=true)
 
 Wing `CL` and `CD` at each angle of attack in `alphas` [deg], solved with the
 settings of the Python verification cases.
 """
-function lift_drag_polar(wing, model, alphas; va, relaxation_factor)
+function lift_drag_polar(wing, model, alphas; va, relaxation_factor, correct_aoa=true)
     body_aero = BodyAerodynamics([wing])
     solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=model,
-                    relaxation_factor, core_radius_fraction=1e-20)
+                    relaxation_factor, core_radius_fraction=1e-20, correct_aoa)
     CL = zeros(length(alphas))
     CD = zeros(length(alphas))
     for (i, alpha) in enumerate(alphas)
@@ -110,10 +110,12 @@ max_error(actual, expected) = maximum(abs.(actual .- expected))
         span = aspect_ratio * π / 4
         coordinates = generate_coordinates_el_wing(1.0, span, 41, "cos")
         wing = wing_from_coordinates(coordinates, INVISCID)
-        CL, CD = lift_drag_polar(wing, VSM, [10.0]; va=10.0, relaxation_factor=0.05)
-        oswald_efficiency = only(CL)^2 / (π * aspect_ratio * only(CD))
+        oswald_efficiency(CL, CD) = only(CL)^2 / (π * aspect_ratio * only(CD))
+        polar(correct_aoa) = lift_drag_polar(wing, VSM, [10.0]; va=10.0,
+                                             relaxation_factor=0.05, correct_aoa)
 
-        @test oswald_efficiency ≈ 1 atol = 0.01
+        @test oswald_efficiency(polar(true)...) ≈ 1 atol = 0.01
+        @test oswald_efficiency(polar(false)...) < 0.8
     end
 
     @testset "curved Clark Y wing matches RANS" begin