diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4ca2653a..db45738c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,8 @@
   doi:10.1088/1742-6596/2767/2/022068): each section gets a drag increment and a force
   along its span from the flow across it, in `solve!`, `solve` and `linearize`. Opt-in
   via `is_with_viscous_drag_correction` (default `false`) on the solver settings.
+- `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
 
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 edac27a7..d7123c47 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -252,6 +252,37 @@ 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)
+
+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)
+    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_range),
+            isnothing(delta_angles) ? nothing : view(delta_angles, section_range))
+    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 84cb5297..d032e697 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
@@ -1177,8 +1175,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 +1186,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),
@@ -1213,10 +1214,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),
     )
@@ -1252,8 +1251,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)
@@ -1264,7 +1263,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{<:Any, U}, body_aero::BodyAerodynamics, y::Vector{T};
         theta_idxs=1:4,
         delta_idxs=nothing,
         va_idxs=nothing,
@@ -1273,22 +1272,12 @@ 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 {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)
@@ -1299,27 +1288,25 @@ 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_vec = 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_vec, om)
+        omega = isnothing(omega_idxs) ? MVector{3, TI}(body_aero_c.omega) : y_in[omega_idxs]
+        set_va!(body_aero_c, va_vec, omega)
 
         solve!(solver_c, body_aero_c; kwargs...)
         solver_c.lr.converged || (n_failed[] += 1)
@@ -1335,13 +1322,13 @@ function linearize(solver::Solver, body_aero::BodyAerodynamics, y::Vector{T};
         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 28121e1d..c2c60e05 100644
--- a/test/body_aerodynamics/test_body_aerodynamics.jl
+++ b/test/body_aerodynamics/test_body_aerodynamics.jl
@@ -514,6 +514,32 @@ function solve_wings(wings)
     return body_aero, solve!(Solver(body_aero), body_aero)
 end
 
+"""
+    wing_pair(section_y, n_panels, offset)
+
+Two `inviscid_wing`s of `n_panels` panels, the second shifted by `-offset` [m] in y.
+"""
+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;
+        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]
@@ -522,9 +548,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
@@ -536,9 +560,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
@@ -546,4 +568,49 @@ end
         @test gamma[n_panels] > 1.05single.gamma_distribution[n_panels]
         @test sol.force[3] > 2single.force[3]
     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)
+        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-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
+        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:2n_wing_sections]) > 0
+        @test jac_fwd ≈ jac_fd rtol=1e-4
+    end
 end