diff --git a/CHANGELOG.md b/CHANGELOG.md
index fd826c31..6cab7d10 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -62,6 +62,9 @@
 
 ### Fixed
 
+- `ELLIPTIC` initial circulation works on a body with more than one wing, where it threw
+  an `ArgumentError`: each wing gets an ellipse over its own span, along its own
+  `spanwise_direction` and centred on its own mid-span, also for a single wing off y = 0.
 - `get_lower_upper`, and with it the flap hinge in `deform_section`, takes the lower and
   upper surface heights where the contour crosses `x = crease_frac`. It took the nearest
   points below and above `y = 0`, which on a cambered section put the hinge near the
diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index 6cae4b09..46482608 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -36,6 +36,7 @@ calculate_relative_alpha_and_relative_velocity
 update_effective_angle_of_attack!
 calculate_stall_angle_list
 wing_span_flip
+spanwise_extent
 calculate_circulation_distribution_elliptical_wing
 _compute_reference_velocity_from_distribution
 smooth_circulation!
diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index ee513b15..0d628876 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -28,7 +28,6 @@ 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 kite body reference frame [m²]
 - `c_ref::Float64` = 1.0: Reference chord length (max panel chord) [m]
-- `y::MVector{P, Float64}` = MVector{P,Float64}(zeros(P))
 - `cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}}` = [LazyBufferCache() for _ in 1:15]
 """
 @with_kw mutable struct BodyAerodynamics{P, W<:AbstractWing, T, PN<:Panel{T}}
@@ -50,7 +49,6 @@ 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)
-    y::MVector{P, T} = zeros(MVector{P, T})
     cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:15]
 end
 
@@ -494,30 +492,26 @@ Returns: nothing
 end
 
 """
-    calculate_circulation_distribution_elliptical_wing(body_aero::BodyAerodynamics, gamma_0=1.0)
+    calculate_circulation_distribution_elliptical_wing(gamma_i, body_aero::BodyAerodynamics,
+                                                       gamma_0=1.0)
 
-Calculate circulation distribution for an elliptical wing.
-
-Returns: nothing
+Write into `gamma_i` an elliptic circulation of peak `gamma_0` over each wing, its control
+points measured along that wing's `spanwise_direction` from the wing's mid-span.
 """
-function calculate_circulation_distribution_elliptical_wing(gamma_i, body_aero::BodyAerodynamics, gamma_0=1.0)
-    length(body_aero.wings) == 1 || throw(ArgumentError("Multiple wings not yet implemented"))
-    
-    wing_span = body_aero.wings[1].span
-    @debug "Wing span: $wing_span"
-    
-    # Calculate y-coordinates of control points
-    y = body_aero.y
-    for (i, panel) in pairs(body_aero.panels) 
-        y[i] = panel.control_point[2] 
+function calculate_circulation_distribution_elliptical_wing(gamma_i,
+        body_aero::BodyAerodynamics, gamma_0=1.0)
+    panel_offset = 0
+    for wing in body_aero.wings
+        lo, hi = spanwise_extent(wing)
+        axis = normalize(wing.spanwise_direction)
+        for i in panel_offset .+ (1:wing.n_panels)
+            span_position = dot(body_aero.panels[i].control_point, axis) - (lo + hi) / 2
+            # Clamped: a control point can lie outside the span of the unrefined sections
+            gamma_i[i] = gamma_0 * sqrt(max(0.0, 1 - (2span_position / (hi - lo))^2))
+        end
+        panel_offset += wing.n_panels
     end
-    
-    # Calculate elliptical distribution (clamp to avoid sqrt of negative
-    # when control points lie outside the nominal span envelope)
-    gamma_i .= gamma_0 * sqrt.(max.(0.0, 1 .- (2 .* y ./ wing_span).^2))
-    
-    @debug "Calculated circulation distribution: $gamma_i"
-    nothing
+    return nothing
 end
 
 """
diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl
index 7329dcc2..c6b3e05d 100644
--- a/src/wing_geometry.jl
+++ b/src/wing_geometry.jl
@@ -1685,24 +1685,25 @@ function refine_mesh_with_billowing!(wing; reuse_aero_data::Bool=false)
 end
 
 """
-    calculate_span(wing::AbstractWing)
+    spanwise_extent(wing::AbstractWing)
 
-Calculate wing span along spanwise direction.
+Lowest and highest projection of the unrefined sections' LE and TE points on `wing`'s
+`spanwise_direction`, as `(lo, hi)` [m].
+"""
+function spanwise_extent(wing::AbstractWing)
+    axis = normalize(wing.spanwise_direction)
+    return extrema(dot(point, axis) for section in wing.unrefined_sections
+                   for point in (section.LE_point, section.TE_point))
+end
 
-Returns:
-    Float64: Wing span
+"""
+    calculate_span(wing::AbstractWing)
+
+Wing span along `spanwise_direction` [m].
 """
 function calculate_span(wing::AbstractWing)
-    # Normalize spanwise direction
-    vector_axis = wing.spanwise_direction ./ norm(wing.spanwise_direction)
-    
-    # Get all points
-    all_points = reduce(vcat, [[section.LE_point, section.TE_point] 
-                              for section in wing.unrefined_sections])
-    
-    # Project points and calculate span
-    projections = [dot(point, vector_axis) for point in all_points]
-    return maximum(projections) - minimum(projections)
+    lo, hi = spanwise_extent(wing)
+    return hi - lo
 end
 
 # Project point onto plane
diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl
index b36001ad..893853c4 100644
--- a/test/body_aerodynamics/test_body_aerodynamics.jl
+++ b/test/body_aerodynamics/test_body_aerodynamics.jl
@@ -616,6 +616,40 @@ end
         @test sol.force[3] > 2single.force[3]
     end
 
+    @testset "ELLIPTIC starts each wing on its own ellipse" begin
+        control_y = [panel.control_point[2] for panel in single_aero.panels]
+        ellipse = sqrt.(1 .- (2control_y ./ span) .^ 2)
+        fin = Wing(n_panels; spanwise_direction=[0.0, 0.0, 1.0])
+        for z in section_y .+ 5.0
+            add_section!(fin, [0.0, 0.0, z], [1.0, 0.0, z], INVISCID)
+        end
+        refine!(fin)
+        body_aero = BodyAerodynamics([wing_pair(section_y, n_panels, span + 1.0); fin])
+        gamma = zeros(3n_panels)
+        VortexStepMethod.calculate_circulation_distribution_elliptical_wing(gamma,
+            body_aero)
+
+        @test gamma ≈ repeat(ellipse, 3)
+
+        lone_wing_gamma = zeros(n_panels)
+        VortexStepMethod.calculate_circulation_distribution_elliptical_wing(
+            lone_wing_gamma, BodyAerodynamics([inviscid_wing(section_y .+ 5.0; n_panels)]))
+        @test lone_wing_gamma ≈ ellipse
+    end
+
+    @testset "ELLIPTIC and ZEROS converge to the same gamma" begin
+        wings = wing_pair(section_y, n_panels, span + 1.0)
+        _, from_zeros = solve_wings(wings)
+        gamma_zeros = copy(from_zeros.gamma_distribution)
+        body_aero = BodyAerodynamics(wings; va_vec=[10.0, 0.0, 1.0])
+        solver = Solver(2n_panels, 2length(section_y);
+                        type_initial_gamma_distribution=ELLIPTIC)
+        sol = solve!(solver, body_aero)
+
+        @test sol.solver_status == FEASIBLE
+        @test sol.gamma_distribution ≈ gamma_zeros rtol=1e-4
+    end
+
     n_wing_sections = length(section_y)
     first_sections = 1:n_wing_sections
     second_sections = n_wing_sections+1:2n_wing_sections