From 580cd2f1d0d98ab36490cc80992746347fc3127f Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Thu, 17 Sep 2026 01:15:18 +0200
Subject: [PATCH 1/3] linearize and make_dual_shadow take a body with several
 wings

theta_idxs and delta_idxs run over the unrefined sections of all wings in
order, the order calc_forces! fills moment_unrefined_dist in.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 CHANGELOG.md                                  |  5 ++
 src/body_aerodynamics.jl                      | 19 ++++++++
 src/solver.jl                                 | 37 +++++----------
 .../test_body_aerodynamics.jl                 | 46 +++++++++++++++++++
 4 files changed, 81 insertions(+), 26 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4c00c754..e79bcd2c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,11 @@
 
 ## Unreleased
 
+### Added
+
+- `linearize` takes a `BodyAerodynamics` with more than one wing; `theta_idxs` and
+  `delta_idxs` then run over the unrefined sections of all wings in order.
+
 ### Changed
 
 - Requires Julia 1.12 or 1.13; 1.10 and 1.11 keep resolving v5.1.1.
diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index d04a966d..b4e12ef2 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -252,6 +252,25 @@ function calculate_stall_angle_list!(stall_angles::AbstractVector,
     return nothing
 end
 
+"""
+    unrefined_deform!(body_aero::BodyAerodynamics, theta_angles, delta_angles)
+
+Deform each wing of `body_aero` by its entries of `theta_angles` and `delta_angles` [rad],
+which run over the unrefined sections of all wings in order; `nothing` leaves that angle
+unchanged. Call [`reinit!`](@ref) afterwards to update the panels.
+"""
+function unrefined_deform!(body_aero::BodyAerodynamics, theta_angles, delta_angles)
+    first_section = 1
+    for wing in body_aero.wings
+        wing_sections = first_section:first_section + wing.n_unrefined_sections - 1
+        unrefined_deform!(wing,
+            isnothing(theta_angles) ? nothing : view(theta_angles, wing_sections),
+            isnothing(delta_angles) ? nothing : view(delta_angles, wing_sections))
+        first_section += wing.n_unrefined_sections
+    end
+    return nothing
+end
+
 """
     reinit!(body_aero::BodyAerodynamics; init_aero, va, omega, refine_mesh, recompute_mapping, sort_sections)
 
diff --git a/src/solver.jl b/src/solver.jl
index c63d44cc..c3d6a429 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -1213,10 +1213,8 @@ buffers are freshly allocated as `TD`-typed.
 function make_dual_shadow(solver::Solver{P, U, Float64},
                           body_aero::BodyAerodynamics{P, W, Float64},
                           ::Type{TD}) where {P, U, W, TD}
-    length(body_aero.wings) == 1 || throw(ArgumentError(
-        "make_dual_shadow currently supports body_aero with one wing"))
-    wing_d = _wing_with_eltype(body_aero.wings[1], TD)
-    body_aero_d = BodyAerodynamics([wing_d];
+    wings_d = [_wing_with_eltype(wing, TD) for wing in body_aero.wings]
+    body_aero_d = BodyAerodynamics(wings_d;
         va = MVector{3, TD}(body_aero._va),
         omega = MVector{3, TD}(body_aero.omega),
     )
@@ -1251,8 +1249,8 @@ end
               backend=AutoForwardDiff(), kwargs...)
 
 Jacobian of aerodynamic outputs w.r.t. control and kinematic inputs at `y`. Each `*_idxs`
-selects which entries of `y` map to twist angles (one per unrefined section), trailing-edge
-deflections (one per unrefined section), apparent wind `(vx, vy, vz)`, and angular rate
+selects which entries of `y` map to twist angles and trailing-edge deflections (one per
+unrefined section, over all wings in order), apparent wind `(vx, vy, vz)`, and angular rate
 `(ωx, ωy, ωz)` respectively.
 
 `backend` accepts any `DifferentiationInterface` backend; `AutoForwardDiff()` (the default)
@@ -1263,7 +1261,7 @@ Returns `(jac, results, converged)` where `results` is `(F, M, moment_unrefined_
 or the corresponding coefficients when `aero_coeffs=true` — and `converged` is `false` (with a
 warning) if any internal solve missed the solver's tolerances.
 """
-function linearize(solver::Solver, body_aero::BodyAerodynamics, y::Vector{T};
+function linearize(solver::Solver{P, U}, body_aero::BodyAerodynamics, y::Vector{T};
         theta_idxs=1:4,
         delta_idxs=nothing,
         va_idxs=nothing,
@@ -1272,22 +1270,11 @@ function linearize(solver::Solver, body_aero::BodyAerodynamics, y::Vector{T};
         backend = AutoForwardDiff(),
         fd_absstep::Float64=1e-8,
         fd_relstep::Float64=1e-8,
-        kwargs...) where T
+        kwargs...) where {P, U, T}
 
-    !(length(body_aero.wings) == 1) && throw(ArgumentError("Linearization only works for a body_aero with one wing"))
-    wing = body_aero.wings[1]
-
-    # Validate that theta_idxs and delta_idxs match the number of unrefined sections
-    if !isnothing(theta_idxs) && wing.n_unrefined_sections > 0
-        length(theta_idxs) != wing.n_unrefined_sections && throw(ArgumentError(
-            "Length of theta_idxs ($(length(theta_idxs))) must match number of unrefined sections ($(wing.n_unrefined_sections))"))
-    end
-    if !isnothing(delta_idxs) && wing.n_unrefined_sections > 0
-        length(delta_idxs) != wing.n_unrefined_sections && throw(ArgumentError(
-            "Length of delta_idxs ($(length(delta_idxs))) must match number of unrefined sections ($(wing.n_unrefined_sections))"))
-    end
-    if wing.n_unrefined_sections == 0 && (!isnothing(theta_idxs) || !isnothing(delta_idxs))
-        throw(ArgumentError("Cannot use theta_idxs or delta_idxs when wing has no unrefined sections"))
+    for (name, idxs) in (("theta_idxs", theta_idxs), ("delta_idxs", delta_idxs))
+        isnothing(idxs) || length(idxs) == U || throw(ArgumentError(
+            "Length of $name ($(length(idxs))) must match number of unrefined sections ($U)"))
     end
 
     n_failed = Ref(0)
@@ -1298,22 +1285,20 @@ function linearize(solver::Solver, body_aero::BodyAerodynamics, y::Vector{T};
         if TI === Float64
             body_aero_c = body_aero
             solver_c = solver
-            wing_c = wing
         else
             shadow = shadow_ref[]
             if shadow === nothing || eltype(shadow[1]._va) !== TI
                 shadow_ref[] = make_dual_shadow(solver, body_aero, TI)
             end
             body_aero_c, solver_c = shadow_ref[]
-            wing_c = body_aero_c.wings[1]
         end
 
         @views theta_angles = isnothing(theta_idxs) ? nothing : y_in[theta_idxs]
         @views delta_angles = isnothing(delta_idxs) ? nothing : y_in[delta_idxs]
 
         if !isnothing(theta_angles) || !isnothing(delta_angles)
-            VortexStepMethod.unrefined_deform!(wing_c, theta_angles, delta_angles; smooth=false)
-            VortexStepMethod.reinit!(body_aero_c; init_aero=false)
+            unrefined_deform!(body_aero_c, theta_angles, delta_angles)
+            reinit!(body_aero_c; init_aero=false)
         end
 
         va = isnothing(va_idxs) ? MVector{3, TI}(body_aero_c._va) : y_in[va_idxs]
diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl
index 037a1b6f..cd88e5f1 100644
--- a/test/body_aerodynamics/test_body_aerodynamics.jl
+++ b/test/body_aerodynamics/test_body_aerodynamics.jl
@@ -512,6 +512,24 @@ function solve_wings(wings)
     return body_aero, solve!(Solver(body_aero), body_aero)
 end
 
+"""
+    linearize_wings(wings; kwargs...)
+
+`linearize` at zero twist and deflection of the body [`solve_wings`](@ref) builds from
+`wings`, over the twist and the deflection of every unrefined section, then the inflow
+and the angular rate.
+"""
+function linearize_wings(wings; kwargs...)
+    n_sections = sum(wing -> wing.n_unrefined_sections, wings)
+    body_aero = BodyAerodynamics(wings; va=[10.0, 0.0, 1.0])
+    solver = Solver(body_aero; use_gamma_prev=false, rtol=1e-10)
+    y0 = [zeros(2n_sections); body_aero.va; zeros(3)]
+    return VortexStepMethod.linearize(solver, body_aero, y0;
+        theta_idxs=1:n_sections, delta_idxs=n_sections+1:2n_sections,
+        va_idxs=2n_sections+1:2n_sections+3, omega_idxs=2n_sections+4:2n_sections+6,
+        kwargs...)
+end
+
 @testset "solve! on a two-wing body" begin
     n_panels = 6
     section_y = [2.0, 0.0, -2.0]
@@ -544,4 +562,32 @@ end
         @test gamma[n_panels] > 1.05single.gamma_distribution[n_panels]
         @test sol.force[3] > 2single.force[3]
     end
+
+    n_sections = 2length(section_y)
+    @testset "linearize: theta of each wing moves that wing's sections" begin
+        wings = [inviscid_wing(section_y; n_panels),
+                 inviscid_wing(section_y .- 1e4; n_panels)]
+        jac, _, converged = linearize_wings(wings)
+        wing_rows = (7:6+length(section_y), 7+length(section_y):6+n_sections)
+        wing_columns = (1:length(section_y), length(section_y)+1:n_sections)
+
+        @test converged
+        @test norm(jac[wing_rows[1], wing_columns[1]]) > 0
+        @test jac[wing_rows[2], wing_columns[2]] ≈ jac[wing_rows[1], wing_columns[1]] rtol=1e-4
+        @test norm(jac[wing_rows[1], wing_columns[2]]) < 1e-4norm(jac[wing_rows[1], wing_columns[1]])
+        @test norm(jac[wing_rows[2], wing_columns[1]]) < 1e-4norm(jac[wing_rows[1], wing_columns[1]])
+    end
+
+    @testset "linearize: AutoForwardDiff matches AutoFiniteDiff" begin
+        wings = [inviscid_wing(section_y; n_panels),
+                 inviscid_wing(section_y .- (span + 1.0); n_panels)]
+        jac_fwd, _, fwd_converged = linearize_wings(wings)
+        jac_fd, _, fd_converged = linearize_wings(wings; backend=nothing,
+            fd_absstep=1e-6, fd_relstep=1e-6)
+
+        @test fwd_converged
+        @test fd_converged
+        @test norm(jac_fwd[:, 1:n_sections]) > 0
+        @test jac_fwd ≈ jac_fd rtol=1e-4
+    end
 end

From 1d76212db1af0bae4ee355b8f1ef4533c7b635f2 Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Thu, 17 Sep 2026 01:18:39 +0200
Subject: [PATCH 2/3] Name linearize's locals and share the two-wing test setup

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/body_aerodynamics.jl                      | 10 ++--
 src/solver.jl                                 | 28 +++++----
 .../test_body_aerodynamics.jl                 | 59 ++++++++++---------
 3 files changed, 52 insertions(+), 45 deletions(-)

diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index b4e12ef2..f47934eb 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -260,13 +260,13 @@ which run over the unrefined sections of all wings in order; `nothing` leaves th
 unchanged. Call [`reinit!`](@ref) afterwards to update the panels.
 """
 function unrefined_deform!(body_aero::BodyAerodynamics, theta_angles, delta_angles)
-    first_section = 1
+    section_offset = 0
     for wing in body_aero.wings
-        wing_sections = first_section:first_section + wing.n_unrefined_sections - 1
+        section_idxs = section_offset .+ (1:wing.n_unrefined_sections)
         unrefined_deform!(wing,
-            isnothing(theta_angles) ? nothing : view(theta_angles, wing_sections),
-            isnothing(delta_angles) ? nothing : view(delta_angles, wing_sections))
-        first_section += wing.n_unrefined_sections
+            isnothing(theta_angles) ? nothing : view(theta_angles, section_idxs),
+            isnothing(delta_angles) ? nothing : view(delta_angles, section_idxs))
+        section_offset += wing.n_unrefined_sections
     end
     return nothing
 end
diff --git a/src/solver.jl b/src/solver.jl
index c3d6a429..80f67ef4 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -1177,8 +1177,10 @@ function _wing_with_eltype(wing::Wing{P, Float64}, ::Type{TD}) where {P, TD}
         wing.spanwise_distribution,
         PanelProperties{P, TD}(),
         MVector{3, TD}(wing.spanwise_direction),
-        Section{TD}[_section_with_eltype(s, TD) for s in wing.unrefined_sections],
-        Section{TD}[_section_with_eltype(s, TD) for s in wing.refined_sections],
+        Section{TD}[_section_with_eltype(section, TD)
+                    for section in wing.unrefined_sections],
+        Section{TD}[_section_with_eltype(section, TD)
+                    for section in wing.refined_sections],
         wing.remove_nan,
         wing.use_prior_polar,
         wing.billowing_percentage,
@@ -1186,7 +1188,8 @@ function _wing_with_eltype(wing::Wing{P, Float64}, ::Type{TD}) where {P, TD}
         copy(wing.refined_panel_mapping),
         copy(wing.refined_section_left_idx),
         Vector{TD}(wing.refined_section_weight),
-        Section{TD}[_section_with_eltype(s, TD) for s in wing.non_deformed_sections],
+        Section{TD}[_section_with_eltype(section, TD)
+                    for section in wing.non_deformed_sections],
         Vector{TD}(wing.theta_dist),
         Vector{TD}(wing.delta_dist),
         TD(wing.mass),
@@ -1261,7 +1264,7 @@ Returns `(jac, results, converged)` where `results` is `(F, M, moment_unrefined_
 or the corresponding coefficients when `aero_coeffs=true` — and `converged` is `false` (with a
 warning) if any internal solve missed the solver's tolerances.
 """
-function linearize(solver::Solver{P, U}, body_aero::BodyAerodynamics, y::Vector{T};
+function linearize(solver::Solver{<:Any, U}, body_aero::BodyAerodynamics, y::Vector{T};
         theta_idxs=1:4,
         delta_idxs=nothing,
         va_idxs=nothing,
@@ -1270,11 +1273,12 @@ function linearize(solver::Solver{P, U}, body_aero::BodyAerodynamics, y::Vector{
         backend = AutoForwardDiff(),
         fd_absstep::Float64=1e-8,
         fd_relstep::Float64=1e-8,
-        kwargs...) where {P, U, T}
+        kwargs...) where {U, T}
 
     for (name, idxs) in (("theta_idxs", theta_idxs), ("delta_idxs", delta_idxs))
         isnothing(idxs) || length(idxs) == U || throw(ArgumentError(
-            "Length of $name ($(length(idxs))) must match number of unrefined sections ($U)"))
+            "Length of $name ($(length(idxs))) must match number of unrefined sections " *
+            "($U)"))
     end
 
     n_failed = Ref(0)
@@ -1302,8 +1306,8 @@ function linearize(solver::Solver{P, U}, body_aero::BodyAerodynamics, y::Vector{
         end
 
         va = isnothing(va_idxs) ? MVector{3, TI}(body_aero_c._va) : y_in[va_idxs]
-        om = isnothing(omega_idxs) ? MVector{3, TI}(body_aero_c.omega) : y_in[omega_idxs]
-        set_va!(body_aero_c, va, om)
+        omega = isnothing(omega_idxs) ? MVector{3, TI}(body_aero_c.omega) : y_in[omega_idxs]
+        set_va!(body_aero_c, va, omega)
 
         solve!(solver_c, body_aero_c; kwargs...)
         solver_c.lr.converged || (n_failed[] += 1)
@@ -1319,13 +1323,13 @@ function linearize(solver::Solver{P, U}, body_aero::BodyAerodynamics, y::Vector{
         return nothing
     end
 
-    n_results = 3 + 3 + length(solver.sol.moment_unrefined_dist)
+    n_results = 3 + 3 + U
     jac = zeros(n_results, length(y))
     results = zeros(n_results)
-    be = backend === nothing ?
+    ad_backend = backend === nothing ?
         AutoFiniteDiff(absstep=fd_absstep, relstep=fd_relstep) : backend
-    prep = prepare_jacobian(calc_results!, results, be, y)
-    jacobian!(calc_results!, results, jac, prep, be, y)
+    prep = prepare_jacobian(calc_results!, results, ad_backend, y)
+    jacobian!(calc_results!, results, jac, prep, ad_backend, y)
     calc_results!(results, y)
     converged = n_failed[] == 0
     if !converged
diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl
index cd88e5f1..daae13df 100644
--- a/test/body_aerodynamics/test_body_aerodynamics.jl
+++ b/test/body_aerodynamics/test_body_aerodynamics.jl
@@ -513,15 +513,23 @@ function solve_wings(wings)
 end
 
 """
-    linearize_wings(wings; kwargs...)
+    wing_pair(section_y, n_panels, offset)
 
-`linearize` at zero twist and deflection of the body [`solve_wings`](@ref) builds from
-`wings`, over the twist and the deflection of every unrefined section, then the inflow
-and the angular rate.
+Two `inviscid_wing`s of `n_panels` panels, the second shifted by `-offset` [m] in y.
 """
-function linearize_wings(wings; kwargs...)
-    n_sections = sum(wing -> wing.n_unrefined_sections, wings)
-    body_aero = BodyAerodynamics(wings; va=[10.0, 0.0, 1.0])
+function wing_pair(section_y, n_panels, offset)
+    return [inviscid_wing(section_y; n_panels),
+            inviscid_wing(section_y .- offset; n_panels)]
+end
+
+"""
+    linearize_body(body_aero; kwargs...)
+
+`linearize` of `body_aero` at zero twist and deflection over the twist and the deflection of
+every unrefined section, then the inflow and the angular rate.
+"""
+function linearize_body(body_aero; kwargs...)
+    n_sections = sum(wing -> wing.n_unrefined_sections, body_aero.wings)
     solver = Solver(body_aero; use_gamma_prev=false, rtol=1e-10)
     y0 = [zeros(2n_sections); body_aero.va; zeros(3)]
     return VortexStepMethod.linearize(solver, body_aero, y0;
@@ -538,9 +546,7 @@ end
     @test single.solver_status == FEASIBLE
 
     @testset "wings far apart each act as the isolated wing" begin
-        wings = [inviscid_wing(section_y; n_panels),
-                 inviscid_wing(section_y .- 1e4; n_panels)]
-        body_aero, sol = solve_wings(wings)
+        body_aero, sol = solve_wings(wing_pair(section_y, n_panels, 1e4))
 
         @test length(body_aero.panels) == 2n_panels
         @test sol.solver_status == FEASIBLE
@@ -552,9 +558,7 @@ end
     end
 
     @testset "wings one chord apart induce on each other" begin
-        wings = [inviscid_wing(section_y; n_panels),
-                 inviscid_wing(section_y .- (span + 1.0); n_panels)]
-        _, sol = solve_wings(wings)
+        _, sol = solve_wings(wing_pair(section_y, n_panels, span + 1.0))
         gamma = sol.gamma_distribution
 
         @test sol.solver_status == FEASIBLE
@@ -563,31 +567,30 @@ end
         @test sol.force[3] > 2single.force[3]
     end
 
-    n_sections = 2length(section_y)
+    n_wing_sections = length(section_y)
     @testset "linearize: theta of each wing moves that wing's sections" begin
-        wings = [inviscid_wing(section_y; n_panels),
-                 inviscid_wing(section_y .- 1e4; n_panels)]
-        jac, _, converged = linearize_wings(wings)
-        wing_rows = (7:6+length(section_y), 7+length(section_y):6+n_sections)
-        wing_columns = (1:length(section_y), length(section_y)+1:n_sections)
+        body_aero, _ = solve_wings(wing_pair(section_y, n_panels, 1e4))
+        jac, _, converged = linearize_body(body_aero)
+        first_sections = 1:n_wing_sections
+        second_sections = n_wing_sections+1:2n_wing_sections
+        own_first = jac[6 .+ first_sections, first_sections]
 
         @test converged
-        @test norm(jac[wing_rows[1], wing_columns[1]]) > 0
-        @test jac[wing_rows[2], wing_columns[2]] ≈ jac[wing_rows[1], wing_columns[1]] rtol=1e-4
-        @test norm(jac[wing_rows[1], wing_columns[2]]) < 1e-4norm(jac[wing_rows[1], wing_columns[1]])
-        @test norm(jac[wing_rows[2], wing_columns[1]]) < 1e-4norm(jac[wing_rows[1], wing_columns[1]])
+        @test norm(own_first) > 0
+        @test jac[6 .+ second_sections, second_sections] ≈ own_first rtol=1e-4
+        @test norm(jac[6 .+ first_sections, second_sections]) < 1e-4norm(own_first)
+        @test norm(jac[6 .+ second_sections, first_sections]) < 1e-4norm(own_first)
     end
 
     @testset "linearize: AutoForwardDiff matches AutoFiniteDiff" begin
-        wings = [inviscid_wing(section_y; n_panels),
-                 inviscid_wing(section_y .- (span + 1.0); n_panels)]
-        jac_fwd, _, fwd_converged = linearize_wings(wings)
-        jac_fd, _, fd_converged = linearize_wings(wings; backend=nothing,
+        body_aero, _ = solve_wings(wing_pair(section_y, n_panels, span + 1.0))
+        jac_fwd, _, fwd_converged = linearize_body(body_aero)
+        jac_fd, _, fd_converged = linearize_body(body_aero; backend=nothing,
             fd_absstep=1e-6, fd_relstep=1e-6)
 
         @test fwd_converged
         @test fd_converged
-        @test norm(jac_fwd[:, 1:n_sections]) > 0
+        @test norm(jac_fwd[:, 1:2n_wing_sections]) > 0
         @test jac_fwd ≈ jac_fd rtol=1e-4
     end
 end

From 574b15b2b65d19f6b4feff28c2ebb376eb31daca Mon Sep 17 00:00:00 2001
From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com>
Date: Thu, 17 Sep 2026 09:46:16 +0200
Subject: [PATCH 3/3] Give each wing's unrefined-section range one source and
 test the angle split

unrefined_section_range(body_aero, wing_idx) replaces the running offsets in
calc_forces! and unrefined_deform!(body_aero, ...). A new testset checks that
each wing receives its own twist and deflection angles, which the INVISCID
Jacobian tests could not see, and the far-apart Jacobian bounds are tightened
to what they measure.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 docs/src/private_functions.md                 |  1 +
 src/body_aerodynamics.jl                      | 24 ++++++++++++----
 src/solver.jl                                 | 10 +++----
 .../test_body_aerodynamics.jl                 | 28 +++++++++++++++----
 4 files changed, 46 insertions(+), 17 deletions(-)

diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index 15f636ef..d6e97fef 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -83,6 +83,7 @@ calculate_filaments_for_plotting
 ### Mesh refinement and billowing
 ```@docs
 unrefined_deform!
+unrefined_section_range
 deform!
 compute_refined_panel_mapping!
 compute_refined_section_interpolation!
diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index 3ff1a9c2..e1fc69f3 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -252,6 +252,20 @@ function calculate_stall_angle_list!(stall_angles::AbstractVector,
     return nothing
 end
 
+"""
+    unrefined_section_range(body_aero::BodyAerodynamics, wing_idx)
+
+Indices of the unrefined sections of wing `wing_idx` in a distribution that runs over the
+unrefined sections of all wings in order, such as `moment_unrefined_dist`.
+"""
+function unrefined_section_range(body_aero::BodyAerodynamics, wing_idx)
+    offset = 0
+    for i in 1:wing_idx-1
+        offset += body_aero.wings[i].n_unrefined_sections
+    end
+    return offset .+ (1:body_aero.wings[wing_idx].n_unrefined_sections)
+end
+
 """
     unrefined_deform!(body_aero::BodyAerodynamics, theta_angles, delta_angles)
 
@@ -260,13 +274,11 @@ which run over the unrefined sections of all wings in order; `nothing` leaves th
 unchanged. Call [`reinit!`](@ref) afterwards to update the panels.
 """
 function unrefined_deform!(body_aero::BodyAerodynamics, theta_angles, delta_angles)
-    section_offset = 0
-    for wing in body_aero.wings
-        section_idxs = section_offset .+ (1:wing.n_unrefined_sections)
+    for (wing_idx, wing) in enumerate(body_aero.wings)
+        section_range = unrefined_section_range(body_aero, wing_idx)
         unrefined_deform!(wing,
-            isnothing(theta_angles) ? nothing : view(theta_angles, section_idxs),
-            isnothing(delta_angles) ? nothing : view(delta_angles, section_idxs))
-        section_offset += wing.n_unrefined_sections
+            isnothing(theta_angles) ? nothing : view(theta_angles, section_range),
+            isnothing(delta_angles) ? nothing : view(delta_angles, section_range))
     end
     return nothing
 end
diff --git a/src/solver.jl b/src/solver.jl
index a6f37048..6d177d28 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -481,13 +481,13 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
         fill!(unrefined_count_dist, 0)
 
         panel_idx = 1
-        unrefined_idx = 1
-        for wing in body_aero.wings
+        for (wing_idx, wing) in enumerate(body_aero.wings)
             if wing.n_unrefined_sections > 0
+                section_range = unrefined_section_range(body_aero, wing_idx)
                 for local_panel_idx in 1:wing.n_panels
                     panel = body_aero.panels[panel_idx]
                     original_section_idx = wing.refined_panel_mapping[local_panel_idx]
-                    target_unrefined_idx = unrefined_idx + original_section_idx - 1
+                    target_unrefined_idx = section_range[original_section_idx]
 
                     # Accumulate coefficients and moments
                     moment_unrefined_dist[target_unrefined_idx] += moment_dist[panel_idx]
@@ -511,8 +511,7 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
 
                 # Average coefficients and geometry. width and
                 # moment_coeff_unrefined_dist stay summed (extensive).
-                for i in 1:wing.n_unrefined_sections
-                    target_unrefined_idx = unrefined_idx + i - 1
+                for target_unrefined_idx in section_range
                     if unrefined_count_dist[target_unrefined_idx] > 0
                         count = unrefined_count_dist[target_unrefined_idx]
                         moment_unrefined_dist[target_unrefined_idx] /= count
@@ -529,7 +528,6 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
                         # sum of panel widths in the unrefined section
                     end
                 end
-                unrefined_idx += wing.n_unrefined_sections
             else
                 # Skip panels for wings with no unrefined sections
                 panel_idx += wing.n_panels
diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl
index daae13df..19cd914c 100644
--- a/test/body_aerodynamics/test_body_aerodynamics.jl
+++ b/test/body_aerodynamics/test_body_aerodynamics.jl
@@ -568,18 +568,36 @@ end
     end
 
     n_wing_sections = length(section_y)
+    first_sections = 1:n_wing_sections
+    second_sections = n_wing_sections+1:2n_wing_sections
+    @testset "unrefined_deform! hands each wing its own run of angles" begin
+        body_aero = BodyAerodynamics(wing_pair(section_y, n_panels, span + 1.0))
+        isolated = wing_pair(section_y, n_panels, span + 1.0)
+        theta = deg2rad.(1.0:2n_wing_sections)
+        delta = -2theta
+        VortexStepMethod.unrefined_deform!(body_aero, theta, delta)
+
+        for (wing_idx, section_range) in enumerate((first_sections, second_sections))
+            wing = isolated[wing_idx]
+            VortexStepMethod.unrefined_deform!(wing, theta[section_range],
+                delta[section_range])
+            @test VortexStepMethod.unrefined_section_range(body_aero, wing_idx) ==
+                  section_range
+            @test body_aero.wings[wing_idx].theta_dist == wing.theta_dist
+            @test body_aero.wings[wing_idx].delta_dist == wing.delta_dist
+        end
+    end
+
     @testset "linearize: theta of each wing moves that wing's sections" begin
         body_aero, _ = solve_wings(wing_pair(section_y, n_panels, 1e4))
         jac, _, converged = linearize_body(body_aero)
-        first_sections = 1:n_wing_sections
-        second_sections = n_wing_sections+1:2n_wing_sections
         own_first = jac[6 .+ first_sections, first_sections]
 
         @test converged
         @test norm(own_first) > 0
-        @test jac[6 .+ second_sections, second_sections] ≈ own_first rtol=1e-4
-        @test norm(jac[6 .+ first_sections, second_sections]) < 1e-4norm(own_first)
-        @test norm(jac[6 .+ second_sections, first_sections]) < 1e-4norm(own_first)
+        @test jac[6 .+ second_sections, second_sections] ≈ own_first rtol=1e-9
+        @test norm(jac[6 .+ first_sections, second_sections]) < 1e-7norm(own_first)
+        @test norm(jac[6 .+ second_sections, first_sections]) < 1e-7norm(own_first)
     end
 
     @testset "linearize: AutoForwardDiff matches AutoFiniteDiff" begin