diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5ad0d3fc..b69c9a1b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,8 @@
 
 ### Added
 
+- `Solver(settings)` and `Solver(n_panels, n_unrefined_sections)` build a solver without
+  a `BodyAerodynamics`; keyword arguments override the settings.
 - `set_va!(body_aero, va_vec, omega; reference_point)` turns the body about
   `reference_point` [m] instead of the origin. The point is stored on
   `BodyAerodynamics`, starts at the origin, and is kept by later `set_va!`, `reinit!`
@@ -24,9 +26,15 @@
 - The Makie `plot!` methods for a `Panel` or a `BodyAerodynamics` return a
   `Vector{Makie.AbstractPlot}` instead of a `Vector{Any}`; for a `BodyAerodynamics`
   drawn as flat panels it is one flat list rather than a list per panel.
+- `Solver(body_aero; kwargs...)` and `Solver(body_aero, settings)` are deprecated and warn
+  on use; build the solver with `Solver(settings)` or
+  `Solver(n_panels, n_unrefined_sections)` instead.
 
 ### Fixed
 
+- `solve!` and `solve` throw a `DimensionMismatch` naming both sizes for a `body_aero` whose
+  panel or unrefined-section count differs from the solver's, where they failed on a
+  broadcast partway through or silently left section results at zero.
 - `set_va!(body_aero, settings)` applies `condition.yaw_rate` as a turn rate about the
   body z axis; it was read from the settings file and ignored.
 - The `VSMSolution` docstring gives `lift_dist`, `drag_dist` and `panel_moment_dist` in
diff --git a/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml b/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml
index e29eb6d5..1e6f80bb 100644
--- a/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml
+++ b/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml
@@ -44,7 +44,7 @@ condition:
 # Define wing geometry files and discretization parameters
 wings:
   - name: V3_Kite # Wing identifier for output labeling
-    geometry_file: data/TUDELFT_V3_KITE/aero_geometry.yaml
+    geometry_file: data/TUDELFT_V3_KITE/aero_geometry_coarse_discretisation.yaml
     n_panels: 54 # Total number of panels along wingspan
     spanwise_panel_distribution: SPLIT_PROVIDED # Panel spacing algorithm
     spanwise_direction: [0.0, 1.0, 0.0] # Unit vector defining wingspan direction
diff --git a/docs/src/examples.md b/docs/src/examples.md
index f9d2fa43..c214d2c7 100644
--- a/docs/src/examples.md
+++ b/docs/src/examples.md
@@ -84,8 +84,8 @@ julia> set_va!(body_aero, va_vec, [0, 0, 0.1])
 #### Step 5: Initialize solvers for both LLT and VSM methods
 
 ```julia
-julia> llt_solver = Solver(body_aero; aerodynamic_model_type=LLT)
-julia> vsm_solver = Solver(body_aero; aerodynamic_model_type=VSM)
+julia> llt_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=LLT)
+julia> vsm_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=VSM)
 ```
 
 #### Step 6: Solve using both methods
diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index eb94971d..ad67330d 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -12,6 +12,9 @@ settings_range
 
 ### Solver, forces and circulation
 ```@docs
+n_unrefined_sections
+solver_kwargs
+check_dimensions
 calculate_AIC_matrices!
 gamma_loop!
 build_spanwise_laplacian!
diff --git a/docs/src/settings.md b/docs/src/settings.md
index 90004454..ab17506b 100644
--- a/docs/src/settings.md
+++ b/docs/src/settings.md
@@ -11,7 +11,7 @@ settings = VSMSettings("my/vsm_settings.yaml"; data_prefix=false)  # as written
 
 wing = Wing(settings)
 body_aero = BodyAerodynamics([wing])
-solver = Solver(body_aero, settings)
+solver = Solver(settings)
 set_va!(body_aero, settings)
 ```
 
diff --git a/examples/V3_kite.jl b/examples/V3_kite.jl
index e3edb46d..443c88c1 100644
--- a/examples/V3_kite.jl
+++ b/examples/V3_kite.jl
@@ -42,7 +42,7 @@ wing = Wing(settings)
 refine!(wing)
 body_aero = BodyAerodynamics([wing])
 VortexStepMethod.reinit!(body_aero)
-solver = Solver(body_aero, settings)
+solver = Solver(settings)
 solver.reference_point .= REFERENCE_POINT
 
 if DEFORM
@@ -75,7 +75,7 @@ if NEURALFOIL
     refine!(wing_nf)
     body_nf = BodyAerodynamics([wing_nf])
     VortexStepMethod.reinit!(body_nf)
-    solver_nf = Solver(body_nf, settings_nf)
+    solver_nf = Solver(settings_nf)
     solver_nf.reference_point .= REFERENCE_POINT
 
     # Reading the generated directory instead of the OBJ shows the airfoils the polar
diff --git a/examples/V3_neuralfoil.jl b/examples/V3_neuralfoil.jl
index 25795625..245a2087 100644
--- a/examples/V3_neuralfoil.jl
+++ b/examples/V3_neuralfoil.jl
@@ -81,14 +81,18 @@ wing_cfd = Wing(settings_cfd)
 refine!(wing_cfd)
 body_cfd = BodyAerodynamics([wing_cfd])
 VortexStepMethod.reinit!(body_cfd)
-solver_cfd = Solver(body_cfd, settings_cfd)
+solver_cfd = Solver(settings_cfd)
 
 println("Creating wing with NeuralFoil polars...")
 wing_nf = Wing(nf_yaml; n_panels=50, spanwise_distribution=LINEAR)
 refine!(wing_nf)
 body_nf = BodyAerodynamics([wing_nf])
 VortexStepMethod.reinit!(body_nf)
-solver_nf = Solver(body_nf, settings_cfd)
+settings_nf = VSMSettings("TUDELFT_V3_KITE/vsm_settings.yaml")
+settings_nf.wings[1].geometry_file = nf_yaml
+settings_nf.solver_settings.relaxation_factor = RELAXATION
+settings_nf.solver_settings.artificial_damping = ARTIFICIAL_DAMPING
+solver_nf = Solver(settings_nf)
 
 # Compare CFD-polar and NeuralFoil-polar wings against published references
 # (Poland 2025 RANS CFD and wind tunnel). `plot_polars` sweeps each solver over the
diff --git a/examples/bench.jl b/examples/bench.jl
index be28e2d3..dc09f9a4 100644
--- a/examples/bench.jl
+++ b/examples/bench.jl
@@ -41,8 +41,8 @@ va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
 set_va!(body_aero, va_vec)
 
 # Step 4: Initialize solvers for both LLT and VSM methods
-llt_solver = Solver(body_aero; aerodynamic_model_type=LLT)
-vsm_solver = Solver(body_aero; aerodynamic_model_type=VSM)
+llt_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=LLT)
+vsm_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=VSM)
 
 # Step 5: Solve using both methods
 results_vsm = solve(vsm_solver, body_aero, nothing)
@@ -66,7 +66,7 @@ body_aero = BodyAerodynamics([wing])
 
 # Create solvers
 vsm_solver = Solver(
-    body_aero;
+    wing.n_panels, wing.n_unrefined_sections;
     aerodynamic_model_type=VSM,
     is_with_artificial_damping=false,
     solver_type=LOOP,
diff --git a/examples/billowing.jl b/examples/billowing.jl
index d6e956d3..e99d0d52 100644
--- a/examples/billowing.jl
+++ b/examples/billowing.jl
@@ -30,15 +30,12 @@ literature_paths = [
         "windtunnel_alpha_sweep_beta_00_0_Poland_2025_Rey_5e5.csv"),
 ]
 
-# Load solver settings (coarse: 54 panels, matches 10-section geometry)
-settings_data = VortexStepMethod.YAML.load_file(
-    joinpath(v3_dir, "vsm_settings_coarse.yaml"))
-condition_cfg = settings_data["condition"]
-solver_cfg = settings_data["solver_settings"]
-wing_cfg = settings_data["wings"][1]
-n_panels = wing_cfg["n_panels"]
+# Coarse settings: 54 panels on the 10-section geometry
+settings = VSMSettings(joinpath(v3_dir, "vsm_settings_coarse.yaml"); data_prefix=false)
+settings.wings[1].geometry_file = joinpath(project_dir, settings.wings[1].geometry_file)
+n_panels = settings.wings[1].n_panels
 
-BILLOWING_PCT = get(wing_cfg, "billowing_percentage", 0.0)
+BILLOWING_PCT = settings.wings[1].billowing_percentage
 
 labels = [
     "VSM flat",
@@ -49,9 +46,7 @@ labels = [
     "WindTunnel Re=5e5",
 ]
 
-# Load coarse geometry (10 structural rib sections)
-geom_data = VortexStepMethod.YAML.load_file(
-    joinpath(v3_dir, "aero_geometry_coarse_discretisation.yaml"))
+geom_data = VortexStepMethod.YAML.load_file(settings.wings[1].geometry_file)
 section_headers = geom_data["wing_sections"]["headers"]
 section_rows = geom_data["wing_sections"]["data"]
 
@@ -84,43 +79,13 @@ body_aero_bill = BodyAerodynamics([wing_bill])
 VortexStepMethod.reinit!(body_aero_bill)
 
 # --- Build solvers ---
-function make_solver(body_aero)
-    Solver(body_aero;
-        solver_type=(solver_cfg["solver_type"] == "NONLIN" ?
-            NONLIN : LOOP),
-        aerodynamic_model_type=getproperty(
-            VortexStepMethod,
-            Symbol(solver_cfg["aerodynamic_model_type"])),
-        density=solver_cfg["density"],
-        max_iterations=solver_cfg["max_iterations"],
-        rtol=solver_cfg["rtol"],
-        tol_reference_error=solver_cfg["tol_reference_error"],
-        relaxation_factor=solver_cfg["relaxation_factor"],
-        is_with_artificial_damping=solver_cfg["artificial_damping"],
-        artificial_damping=(
-            k2=solver_cfg["k2"], k4=solver_cfg["k4"]),
-        type_initial_gamma_distribution=getproperty(
-            VortexStepMethod,
-            Symbol(solver_cfg["type_initial_gamma_distribution"])),
-        use_gamma_prev=get(solver_cfg, "use_gamma_prev",
-                           get(solver_cfg, "use_gamme_prev", true)),
-        core_radius_fraction=solver_cfg["core_radius_fraction"],
-        mu=solver_cfg["mu"],
-        is_only_f_and_gamma_output=get(
-            solver_cfg, "calc_only_f_and_gamma", false),
-        correct_aoa=get(solver_cfg, "correct_aoa", false),
-        reference_point=get(solver_cfg, "reference_point",
-                            [0.422646, 0.0, 9.3667]),
-    )
-end
-
-solver_flat = make_solver(body_aero_flat)
-solver_bill = make_solver(body_aero_bill)
+solver_flat = Solver(settings; reference_point=[0.422646, 0.0, 9.3667])
+solver_bill = Solver(settings; reference_point=[0.422646, 0.0, 9.3667])
 
 # --- Set flight conditions ---
-wind_speed = condition_cfg["wind_speed"]
+wind_speed = settings.condition.wind_speed
 angle_of_attack_deg = 10.0
-sideslip_deg = condition_cfg["beta"]
+sideslip_deg = settings.condition.beta
 
 α0 = deg2rad(angle_of_attack_deg)
 β0 = deg2rad(sideslip_deg)
diff --git a/examples/linearize_check.jl b/examples/linearize_check.jl
index 63791e5b..b8bd47fc 100644
--- a/examples/linearize_check.jl
+++ b/examples/linearize_check.jl
@@ -23,7 +23,7 @@ yaml = obj_to_yaml(
 wing = Wing(yaml; n_panels=16)
 body_aero = BodyAerodynamics([wing])
 
-solver = Solver(body_aero;
+solver = Solver(wing.n_panels, wing.n_unrefined_sections;
     aerodynamic_model_type=VSM,
     is_with_artificial_damping=false,
     rtol=1e-7,
diff --git a/examples/obj_to_yaml_kite.jl b/examples/obj_to_yaml_kite.jl
index e756b2e2..3a4d9354 100644
--- a/examples/obj_to_yaml_kite.jl
+++ b/examples/obj_to_yaml_kite.jl
@@ -72,7 +72,8 @@ refine!(wing)
 body_aero = BodyAerodynamics([wing])
 VortexStepMethod.reinit!(body_aero)
 
-solver = Solver(body_aero; aerodynamic_model_type=VSM, rtol=1e-5, solver_type=LOOP)
+solver = Solver(wing.n_panels, wing.n_unrefined_sections;
+    aerodynamic_model_type=VSM, rtol=1e-5, solver_type=LOOP)
 set_va!(body_aero, [cos(deg2rad(8)) * va, 0.0, sin(deg2rad(8)) * va])
 results = VortexStepMethod.solve(solver, body_aero; log=true)
 
diff --git a/examples/pyramid_model.jl b/examples/pyramid_model.jl
index c6fe02c3..604ee950 100644
--- a/examples/pyramid_model.jl
+++ b/examples/pyramid_model.jl
@@ -14,7 +14,7 @@ vsm_settings = VSMSettings("pyramid_model/vsm_settings.yaml")
 wing = Wing(vsm_settings)
 refine!(wing)
 body_aero = BodyAerodynamics([wing])
-solver = Solver(body_aero, vsm_settings)
+solver = Solver(vsm_settings)
 
 # Set flight conditions from settings
 set_va!(body_aero, vsm_settings)
diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl
index dc09ad6f..2c9b790b 100644
--- a/examples/ram_air_kite.jl
+++ b/examples/ram_air_kite.jl
@@ -80,13 +80,15 @@ fig_audit = plot_slices_3d(joinpath("data", "ram_air_kite", "polars_xfoil");
                            delta=1.0, obj_path=obj_path)
 GLMakie.save("ram_air_slices_audit.png", fig_audit)
 body_xfoil = BodyAerodynamics([wing_xfoil])
-solver_xfoil = Solver(body_xfoil; aerodynamic_model_type=VSM, rtol=1e-5, solver_type=LOOP,
+solver_xfoil = Solver(wing_xfoil.n_panels, wing_xfoil.n_unrefined_sections;
+                      aerodynamic_model_type=VSM, rtol=1e-5, solver_type=LOOP,
                       relaxation_factor=RELAXATION, is_with_artificial_damping=ARTIFICIAL_DAMPING)
 
 println("Creating NeuralFoil wing...")
 wing_nf = matrix_wing(NF_SOLVER, "polars_neuralfoil")
 body_nf = BodyAerodynamics([wing_nf])
-solver_nf = Solver(body_nf; aerodynamic_model_type=VSM, rtol=1e-5, solver_type=LOOP,
+solver_nf = Solver(wing_nf.n_panels, wing_nf.n_unrefined_sections;
+                   aerodynamic_model_type=VSM, rtol=1e-5, solver_type=LOOP,
                    relaxation_factor=RELAXATION, is_with_artificial_damping=ARTIFICIAL_DAMPING)
 
 # Compare using plot_polars
diff --git a/examples/rectangular_wing.jl b/examples/rectangular_wing.jl
index 11aba7e3..15119a9a 100644
--- a/examples/rectangular_wing.jl
+++ b/examples/rectangular_wing.jl
@@ -45,8 +45,8 @@ va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
 set_va!(body_aero, va_vec, [0, 0, 0.1])
 
 # Step 4: Initialize solvers for both LLT and VSM methods
-llt_solver = Solver(body_aero; aerodynamic_model_type=LLT)
-vsm_solver = Solver(body_aero; aerodynamic_model_type=VSM)
+llt_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=LLT)
+vsm_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=VSM)
 
 # Step 5: Solve using both methods
 results_llt = solve(llt_solver, body_aero)
diff --git a/examples/stall_model.jl b/examples/stall_model.jl
index e504eff3..168fd3bb 100644
--- a/examples/stall_model.jl
+++ b/examples/stall_model.jl
@@ -46,11 +46,11 @@ refine!(CAD_wing)
 body_aero = BodyAerodynamics([CAD_wing])
 
 # Create solvers
-vsm_solver = Solver(body_aero;
+vsm_solver = Solver(CAD_wing.n_panels, CAD_wing.n_unrefined_sections;
     aerodynamic_model_type=VSM,
     is_with_artificial_damping=false
 )
-VSM_with_stall_correction = Solver(body_aero;
+VSM_with_stall_correction = Solver(CAD_wing.n_panels, CAD_wing.n_unrefined_sections;
     aerodynamic_model_type=VSM,
     is_with_artificial_damping=true
 )
diff --git a/mwes/mwe_warntype.jl b/mwes/mwe_warntype.jl
index 1e5f714c..664a0787 100644
--- a/mwes/mwe_warntype.jl
+++ b/mwes/mwe_warntype.jl
@@ -27,7 +27,7 @@ refine!(wing)
 body_aero = BodyAerodynamics([wing])
 va_vec = [cos(alpha), 0.0, sin(alpha)] .* 20.0
 set_va!(body_aero, va_vec)
-solver = Solver(body_aero)
+solver = Solver(wing.n_panels, wing.n_unrefined_sections)
 
 va_dist = ones(n_panels)
 va_unit_dist = ones(n_panels, 3)
diff --git a/src/solver.jl b/src/solver.jl
index 9f2e5ede..04bcc7f4 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -194,37 +194,82 @@ sol::VSMSolution = VSMSolution(): The result of calling [`solve!`](@ref)
     sol::VSMSolution{P, U, T} = VSMSolution{P, U, T}()
 end
 
-function Solver(body_aero::BodyAerodynamics{P, W, T}; reference_point=[0.0, 0.0, 0.0], kwargs...) where {P, W, T}
-    U = sum([wing.n_unrefined_sections for wing in body_aero.wings])
+"""
+    Solver(n_panels, n_unrefined_sections, T=Float64; reference_point=[0.0, 0.0, 0.0],
+           kwargs...)
+    Solver(settings::VSMSettings; kwargs...)
+
+Build a [`Solver`](@ref) for `n_panels` panels and `n_unrefined_sections` unrefined
+sections of element type `T`, with `kwargs` setting its fields. `settings` supplies the
+counts from the `n_panels` and `geometry_file` of its wings, and from its
+`solver_settings` the fields `kwargs` leaves unset. [`solve!`](@ref) throws a
+`DimensionMismatch` for a body of other counts.
+
+`Solver(body_aero; kwargs...)` and `Solver(body_aero, settings)` are deprecated.
+"""
+function Solver(n_panels::Integer, n_unrefined_sections::Integer, ::Type{T}=Float64;
+        reference_point=[0.0, 0.0, 0.0], kwargs...) where {T}
     reference_point_checked = check_reference_point(reference_point, T)
-    return Solver{P, U, T}(; reference_point=reference_point_checked, kwargs...)
+    return Solver{Int(n_panels), Int(n_unrefined_sections), T}(;
+        reference_point=reference_point_checked, kwargs...)
+end
+
+function Solver(settings::VSMSettings; kwargs...)
+    n_panels = sum(wing.n_panels for wing in settings.wings)
+    n_sections = sum(n_unrefined_sections, settings.wings)
+    return Solver(n_panels, n_sections; solver_kwargs(settings.solver_settings)...,
+        kwargs...)
 end
 
-function Solver(body_aero, settings::VSMSettings)
-    ss = settings.solver_settings
-    solver_type = ss.solver_type == "NONLIN" ? NONLIN : LOOP
-    reference_point = hasproperty(ss, :reference_point) ? ss.reference_point : [0.0, 0.0, 0.0]
-    Solver(body_aero;
-        solver_type,
-        aerodynamic_model_type=ss.aerodynamic_model_type,
-        density=ss.density,
-        max_iterations=ss.max_iterations,
-        rtol=ss.rtol,
-        tol_reference_error=ss.tol_reference_error,
-        relaxation_factor=ss.relaxation_factor,
-        is_with_artificial_damping=ss.artificial_damping,
-        artificial_damping=(k2=ss.k2, k4=ss.k4),
-        is_with_artificial_viscosity=ss.is_with_artificial_viscosity,
-        artificial_viscosity_factor=ss.artificial_viscosity_factor,
-        type_initial_gamma_distribution=ss.type_initial_gamma_distribution,
-        use_gamma_prev=ss.use_gamma_prev,
-        core_radius_fraction=ss.core_radius_fraction,
-        mu=ss.mu,
-        is_only_f_and_gamma_output=ss.calc_only_f_and_gamma,
-        correct_aoa=ss.correct_aoa,
-        flow_curvature=ss.flow_curvature,
-        is_with_viscous_drag_correction=ss.is_with_viscous_drag_correction,
-        reference_point=reference_point,
+function Solver(body_aero::BodyAerodynamics{P, W, T}; kwargs...) where {P, W, T}
+    Base.depwarn("`Solver(body_aero; kwargs...)` is deprecated, use " *
+        "`Solver(n_panels, n_unrefined_sections; kwargs...)` or " *
+        "`Solver(settings; kwargs...)`.", :Solver; force=true)
+    return Solver(P, n_unrefined_sections(body_aero), T; kwargs...)
+end
+
+function Solver(body_aero::BodyAerodynamics{P, W, T}, settings::VSMSettings
+        ) where {P, W, T}
+    Base.depwarn("`Solver(body_aero, settings)` is deprecated, use `Solver(settings)`.",
+        :Solver; force=true)
+    return Solver(P, n_unrefined_sections(body_aero), T;
+        solver_kwargs(settings.solver_settings)...)
+end
+
+"""
+    n_unrefined_sections(body_aero::BodyAerodynamics) -> Int
+
+Number of unrefined sections summed over the wings of `body_aero`.
+"""
+n_unrefined_sections(body_aero::BodyAerodynamics) =
+    sum(wing -> wing.n_unrefined_sections, body_aero.wings)
+
+"""
+    solver_kwargs(solver_settings::SolverSettings) -> NamedTuple
+
+The [`Solver`](@ref) fields that `solver_settings` sets, as keyword arguments.
+"""
+function solver_kwargs(solver_settings::SolverSettings)
+    return (
+        solver_type=solver_settings.solver_type == "NONLIN" ? NONLIN : LOOP,
+        aerodynamic_model_type=solver_settings.aerodynamic_model_type,
+        density=solver_settings.density,
+        max_iterations=solver_settings.max_iterations,
+        rtol=solver_settings.rtol,
+        tol_reference_error=solver_settings.tol_reference_error,
+        relaxation_factor=solver_settings.relaxation_factor,
+        is_with_artificial_damping=solver_settings.artificial_damping,
+        artificial_damping=(k2=solver_settings.k2, k4=solver_settings.k4),
+        is_with_artificial_viscosity=solver_settings.is_with_artificial_viscosity,
+        artificial_viscosity_factor=solver_settings.artificial_viscosity_factor,
+        type_initial_gamma_distribution=solver_settings.type_initial_gamma_distribution,
+        use_gamma_prev=solver_settings.use_gamma_prev,
+        core_radius_fraction=solver_settings.core_radius_fraction,
+        mu=solver_settings.mu,
+        is_only_f_and_gamma_output=solver_settings.calc_only_f_and_gamma,
+        correct_aoa=solver_settings.correct_aoa,
+        flow_curvature=solver_settings.flow_curvature,
+        is_with_viscous_drag_correction=solver_settings.is_with_viscous_drag_correction,
     )
 end
 
@@ -622,6 +667,21 @@ end
     end
 end
 
+"""
+    check_dimensions(solver::Solver, body_aero::BodyAerodynamics)
+
+Throw a `DimensionMismatch` unless `solver` was built for as many panels and unrefined
+sections as `body_aero` has.
+"""
+function check_dimensions(::Solver{P, U}, body_aero::BodyAerodynamics) where {P, U}
+    n_panels = length(body_aero.panels)
+    n_sections = n_unrefined_sections(body_aero)
+    n_panels == P && n_sections == U || throw(DimensionMismatch(
+        "Solver built for $P panels and $U unrefined sections is given a body_aero " *
+        "with $n_panels panels and $n_sections unrefined sections"))
+    return nothing
+end
+
 """
     solve_base!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=nothing;
                 log=false)
@@ -637,6 +697,7 @@ function solve_base!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma
                log=false) where {P, U, T}
     
     # check arguments
+    check_dimensions(solver, body_aero)
     isnothing(body_aero.panels[1].va) && throw(ArgumentError(
         "Inflow conditions are not set, use set_va!(body_aero, va_vec)"))
     
@@ -1218,7 +1279,7 @@ function make_dual_shadow(solver::Solver{P, U, Float64},
     body_aero_d = BodyAerodynamics(wings_d)
     set_va!(body_aero_d, MVector{3, TD}(body_aero._va), MVector{3, TD}(body_aero.omega);
             reference_point=body_aero.reference_point)
-    solver_d = Solver(body_aero_d;
+    solver_d = Solver(P, U, TD;
         solver_type = solver.solver_type,
         aerodynamic_model_type = solver.aerodynamic_model_type,
         density = TD(solver.density),
diff --git a/src/yaml_geometry.jl b/src/yaml_geometry.jl
index 11da48ac..456c3656 100644
--- a/src/yaml_geometry.jl
+++ b/src/yaml_geometry.jl
@@ -427,6 +427,18 @@ function Wing(settings::VSMSettings; sort_sections::Bool=true)
     end
 end
 
+"""
+    n_unrefined_sections(wing_settings::WingSettings) -> Int
+
+Number of unrefined sections the [`Wing`](@ref) built from `wing_settings` carries: the
+rows of its `geometry_file`'s `wing_sections`.
+"""
+function n_unrefined_sections(wing_settings::WingSettings)
+    isempty(wing_settings.geometry_file) && throw(ArgumentError(
+        "wing \"$(wing_settings.name)\" has no geometry_file to count its sections from"))
+    return length(YAML.load_file(wing_settings.geometry_file)["wing_sections"]["data"])
+end
+
 """
     ObjWing(obj_path[, dat_path]; n_panels, Re, alpha_range, delta_range,
             n_sections, spanwise_direction, aero_solver, remake, output_dir,
diff --git a/test/bench.jl b/test/bench.jl
index 887f25a0..197e3e33 100644
--- a/test/bench.jl
+++ b/test/bench.jl
@@ -72,8 +72,8 @@ using LinearAlgebra
     set_va!(body_aero, va_vec)
 
     # Initialize solvers for both LLT and VSM methods
-    solver = Solver(body_aero)
-    nonlin_solver = Solver(body_aero; solver_type=NONLIN)
+    solver = Solver(wing.n_panels, wing.n_unrefined_sections)
+    nonlin_solver = Solver(wing.n_panels, wing.n_unrefined_sections; solver_type=NONLIN)
 
     # Pre-allocate arrays
     gamma = rand(n_panels)
@@ -142,7 +142,7 @@ using LinearAlgebra
                 refine!(wing)
                 body_aero = BodyAerodynamics([wing])
                 
-                solver = Solver(body_aero;
+                solver = Solver(wing.n_panels, wing.n_unrefined_sections;
                     aerodynamic_model_type=model
                 )
                 solver.sol._va_dist .= va_vec_dist
diff --git a/test/bench_solve.jl b/test/bench_solve.jl
index 1ede63d8..a96832f3 100644
--- a/test/bench_solve.jl
+++ b/test/bench_solve.jl
@@ -38,7 +38,7 @@ va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
 set_va!(body_aero, va_vec)
 
 # Step 4: Initialize solvers for both LLT and VSM methods
-vsm_solver = Solver(body_aero; aerodynamic_model_type=VSM)
+vsm_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=VSM)
 
 # Step 5: Solve using both methods
 result = @benchmark  solve_base!($vsm_solver, $body_aero, nothing)  # 34 allocations
diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl
index 2fbe4b2c..7024b410 100644
--- a/test/body_aerodynamics/test_body_aerodynamics.jl
+++ b/test/body_aerodynamics/test_body_aerodynamics.jl
@@ -327,7 +327,7 @@ end
     set_va!(body_aero, va_vec)
 
     # Run analysis
-    loop_solver = Solver(body_aero;
+    loop_solver = Solver(wing.n_panels, wing.n_unrefined_sections;
         aerodynamic_model_type=model,
         core_radius_fraction=core_radius_fraction,
         solver_type=LOOP,
@@ -335,7 +335,7 @@ end
         atol=1e-8,
         rtol=1e-8
     )
-    nonlin_solver = Solver(body_aero;
+    nonlin_solver = Solver(wing.n_panels, wing.n_unrefined_sections;
         aerodynamic_model_type=model,
         core_radius_fraction=core_radius_fraction,
         solver_type=NONLIN,
@@ -548,7 +548,9 @@ The `BodyAerodynamics` built from `wings` in a 10 m/s inflow and its `solve!` so
 """
 function solve_wings(wings)
     body_aero = BodyAerodynamics(wings; va=[10.0, 0.0, 1.0])
-    return body_aero, solve!(Solver(body_aero), body_aero)
+    solver = Solver(sum(wing -> wing.n_panels, wings),
+                    sum(wing -> wing.n_unrefined_sections, wings))
+    return body_aero, solve!(solver, body_aero)
 end
 
 """
@@ -569,7 +571,8 @@ 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)
+    solver = Solver(sum(wing -> wing.n_panels, body_aero.wings), n_sections;
+                    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,
diff --git a/test/body_aerodynamics/test_results.jl b/test/body_aerodynamics/test_results.jl
index 2d3667d1..98fa96ad 100644
--- a/test/body_aerodynamics/test_results.jl
+++ b/test/body_aerodynamics/test_results.jl
@@ -32,7 +32,7 @@ end
 
     VortexStepMethod.unrefined_deform!(ram_wing, theta, delta; smooth=false)
     body_aero = BodyAerodynamics([ram_wing]; va=va_vec, omega)
-    solver = Solver(body_aero;
+    solver = Solver(ram_wing.n_panels, ram_wing.n_unrefined_sections;
         aerodynamic_model_type=VSM,
         is_with_artificial_damping=false,
         atol=1e-5,
diff --git a/test/plotting/test_plotting.jl b/test/plotting/test_plotting.jl
index 669bf42f..38354e8c 100644
--- a/test/plotting/test_plotting.jl
+++ b/test/plotting/test_plotting.jl
@@ -78,8 +78,11 @@ end
     end
 
     # Initialize the solvers
-    vsm_solver = Solver(body_aero; aerodynamic_model_type=VSM)
-    llt_solver = Solver(body_aero; aerodynamic_model_type=LLT)
+    wing = only(body_aero.wings)
+    vsm_solver = Solver(wing.n_panels, wing.n_unrefined_sections;
+        aerodynamic_model_type=VSM)
+    llt_solver = Solver(wing.n_panels, wing.n_unrefined_sections;
+        aerodynamic_model_type=LLT)
 
     # Solve the VSM and LLT
     results_vsm = solve(vsm_solver, body_aero)
diff --git a/test/solver/test_backend_comparison.jl b/test/solver/test_backend_comparison.jl
index 3cc47b63..849114ee 100644
--- a/test/solver/test_backend_comparison.jl
+++ b/test/solver/test_backend_comparison.jl
@@ -39,10 +39,13 @@ using Test
         return wing
     end
     function lift(matrices)
-        body = BodyAerodynamics([flat_wing(matrices)])
+        wing = flat_wing(matrices)
+        body = BodyAerodynamics([wing])
         aoa = deg2rad(4.0)
         set_va!(body, [cos(aoa), 0.0, sin(aoa)] * 15.0)
-        return collect(solve!(Solver(body; aerodynamic_model_type=VSM), body).force_coeffs)
+        solver = Solver(wing.n_panels, wing.n_unrefined_sections;
+            aerodynamic_model_type=VSM)
+        return collect(solve!(solver, body).force_coeffs)
     end
     fc_nf, fc_xf = lift(nf), lift(xf)
     @test fc_nf[3] > 0 && fc_xf[3] > 0          # positive lift at positive alpha
diff --git a/test/solver/test_flow_curvature.jl b/test/solver/test_flow_curvature.jl
index 7c895ca9..2639e344 100644
--- a/test/solver/test_flow_curvature.jl
+++ b/test/solver/test_flow_curvature.jl
@@ -102,8 +102,8 @@ end
         @test sign(dot(dv, panel.z_airf)) == sign(q_local)
     end
 
-    solver_off = Solver(body_aero; flow_curvature=false)
-    solver_on = Solver(body_aero; flow_curvature=true)
+    solver_off = Solver(wing.n_panels, wing.n_unrefined_sections; flow_curvature=false)
+    solver_on = Solver(wing.n_panels, wing.n_unrefined_sections; flow_curvature=true)
 
     function moment_at(solver, omega)
         set_va!(body_aero, [V, 0.0, 0.0], omega)
@@ -128,7 +128,7 @@ end
     end
 
     @testset "defaults to off" begin
-        @test Solver(body_aero).flow_curvature == false
+        @test Solver(wing.n_panels, wing.n_unrefined_sections).flow_curvature == false
         @test VortexStepMethod.SolverSettings().flow_curvature == false
     end
 
diff --git a/test/solver/test_forwarddiff.jl b/test/solver/test_forwarddiff.jl
index b23fa8b0..de4d6c8f 100644
--- a/test/solver/test_forwarddiff.jl
+++ b/test/solver/test_forwarddiff.jl
@@ -27,7 +27,7 @@ relative_error(jac, reference) = maximum(abs.(jac .- reference)) / maximum(abs,
         pivot_body = BodyAerodynamics([wing])
         set_va!(pivot_body, va_vec, omega_op; reference_point)
         y_op = [va_vec; omega_op]
-        solver = Solver(pivot_body;
+        solver = Solver(wing.n_panels, wing.n_unrefined_sections;
             use_gamma_prev=false,
             type_initial_gamma_distribution=ELLIPTIC)
 
@@ -49,7 +49,7 @@ relative_error(jac, reference) = maximum(abs.(jac .- reference)) / maximum(abs,
     end
 
     @testset "NONLIN+ForwardDiff is rejected" begin
-        solver_nl = Solver(body_aero; solver_type=NONLIN)
+        solver_nl = Solver(wing.n_panels, wing.n_unrefined_sections; solver_type=NONLIN)
         @test_throws ErrorException VortexStepMethod.linearize(
             solver_nl, body_aero, y0;
             theta_idxs=nothing, va_idxs=1:3, omega_idxs=4:6,
@@ -70,7 +70,7 @@ relative_error(jac, reference) = maximum(abs.(jac .- reference)) / maximum(abs,
             delta_range=deg2rad.(-3:3:3),
         )
         ram_body = BodyAerodynamics([ram_wing])
-        ram_solver = Solver(ram_body;
+        ram_solver = Solver(ram_wing.n_panels, ram_wing.n_unrefined_sections;
             aerodynamic_model_type=VSM,
             is_with_artificial_damping=false,
             rtol=1e-11,
diff --git a/test/solver/test_moment_units.jl b/test/solver/test_moment_units.jl
index baa67e7f..62bf2f9b 100644
--- a/test/solver/test_moment_units.jl
+++ b/test/solver/test_moment_units.jl
@@ -31,7 +31,10 @@ end
     @testset "solve!" begin
         small, large = map((1.0, k)) do scale
             body_aero = scaled_wing_aero(scale)
-            solve!(Solver(body_aero; reference_point=reference_point(scale)), body_aero)
+            wing = only(body_aero.wings)
+            solver = Solver(wing.n_panels, wing.n_unrefined_sections;
+                reference_point=reference_point(scale))
+            solve!(solver, body_aero)
         end
         @test small.solver_status == large.solver_status == FEASIBLE
         @test all(!iszero, small.panel_moment_dist)
@@ -48,7 +51,9 @@ end
     @testset "solve" begin
         small, large = map((1.0, k)) do scale
             body_aero = scaled_wing_aero(scale)
-            solve(Solver(body_aero), body_aero; reference_point=reference_point(scale))
+            wing = only(body_aero.wings)
+            solver = Solver(wing.n_panels, wing.n_unrefined_sections)
+            solve(solver, body_aero; reference_point=reference_point(scale))
         end
         for key in ("Mx", "My", "Mz", "M_distribution")
             @test large[key] ≈ k^3 .* small[key] rtol = 1e-6
diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl
index 9b119177..2dceee9d 100644
--- a/test/solver/test_solver.jl
+++ b/test/solver/test_solver.jl
@@ -18,11 +18,12 @@ end
             wing = Wing(settings)
             refine!(wing)
             body_aero = BodyAerodynamics([wing])
-            solver = Solver(body_aero, settings)
+            solver = Solver(settings)
 
             # Verify solver properties match settings
             @test solver.aerodynamic_model_type == VSM
             @test solver.density == 1.225
+            @test Solver(settings; density=1.0).density == 1.0
 
             # Test that the solver can solve
             va_vec = [10.0, 0.0, 0.0]
@@ -30,6 +31,37 @@ end
             sol = solve!(solver, body_aero)
             @test sol isa VSMSolution
 
+            @testset "body_aero constructors warn and match Solver(settings)" begin
+                with_settings = r"`Solver\(body_aero, settings\)` is deprecated"
+                solver_from_body = @test_logs((:warn, with_settings),
+                    Solver(body_aero, settings))
+                @test solver_from_body isa typeof(solver)
+                @test solver_from_body.density == solver.density
+                @test solve!(solver_from_body, body_aero).force ≈ sol.force
+
+                with_kwargs = r"`Solver\(body_aero; kwargs...\)` is deprecated"
+                solver_from_body = @test_logs((:warn, with_kwargs),
+                    Solver(body_aero; density=1.0))
+                @test solver_from_body isa typeof(solver)
+                @test solver_from_body.density == 1.0
+            end
+
+            @testset "Solver from panel and section counts" begin
+                n_sections = wing.n_unrefined_sections
+                solver_from_counts = Solver(wing.n_panels, n_sections)
+                @test solver_from_counts isa typeof(solver)
+                @test solve!(solver_from_counts, body_aero).force ≈ sol.force
+            end
+
+            @testset "solve refuses a body_aero sized for another solver" begin
+                n_sections = wing.n_unrefined_sections
+                for other in (Solver(wing.n_panels + 1, n_sections),
+                              Solver(wing.n_panels, n_sections + 1))
+                    @test_throws DimensionMismatch solve!(other, body_aero)
+                    @test_throws "Solver built for" solve!(other, body_aero)
+                    @test_throws DimensionMismatch solve(other, body_aero)
+                end
+            end
         finally
             # Cleanup
             rm(settings_file; force=true)
@@ -49,7 +81,7 @@ end
 
         body_aero = BodyAerodynamics([wing])
         solver = Solver(
-            body_aero;
+            wing.n_panels, wing.n_unrefined_sections;
             solver_type=NONLIN,
             aerodynamic_model_type=VSM,
             type_initial_gamma_distribution=ELLIPTIC,
@@ -82,10 +114,10 @@ end
         refine!(wing)
         body_aero = BodyAerodynamics([wing])
         va_vec = [10.0, 0.0, 5.0]   # 26.6 deg angle of attack, past stall
-        nonlin = Solver(body_aero; solver_type=NONLIN, aerodynamic_model_type=VSM,
-            type_initial_gamma_distribution=ELLIPTIC)
-        loop = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM,
-            type_initial_gamma_distribution=ELLIPTIC)
+        nonlin = Solver(wing.n_panels, wing.n_unrefined_sections; solver_type=NONLIN,
+            aerodynamic_model_type=VSM, type_initial_gamma_distribution=ELLIPTIC)
+        loop = Solver(wing.n_panels, wing.n_unrefined_sections; solver_type=LOOP,
+            aerodynamic_model_type=VSM, type_initial_gamma_distribution=ELLIPTIC)
 
         set_va!(body_aero, va_vec)
         sol_nonlin = solve!(nonlin, body_aero)
@@ -109,8 +141,9 @@ One unrelaxed fixed-point step `F(gamma)` of the LOOP iteration, so that
 """
 function unrelaxed_step(body_aero, gamma)
     # An infinite rtol accepts the single step, so solve_base! skips its retry.
-    probe = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM,
-        relaxation_factor=1.0, max_iterations=1, rtol=Inf)
+    wing = only(body_aero.wings)
+    probe = Solver(wing.n_panels, wing.n_unrefined_sections; solver_type=LOOP,
+        aerodynamic_model_type=VSM, relaxation_factor=1.0, max_iterations=1, rtol=Inf)
     VortexStepMethod.solve_base!(probe, body_aero, gamma)
     return copy(probe.lr.gamma_new)
 end
@@ -125,8 +158,8 @@ end
         wing = Wing(settings)
         refine!(wing)
         body_aero = BodyAerodynamics([wing])
-        solver = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM,
-            type_initial_gamma_distribution=ELLIPTIC)
+        solver = Solver(wing.n_panels, wing.n_unrefined_sections; solver_type=LOOP,
+            aerodynamic_model_type=VSM, type_initial_gamma_distribution=ELLIPTIC)
 
         # 0 deg, and 26.6 deg past stall
         for va_vec in ([10.0, 0.0, 0.0], [10.0, 0.0, 5.0])
@@ -154,7 +187,7 @@ calc_forces_allocs(solver, body_aero) =
         wing = Wing(settings)
         refine!(wing)
         body_aero = BodyAerodynamics([wing])
-        solver = Solver(body_aero, settings)
+        solver = Solver(settings)
         set_va!(body_aero, [10.0, 0.0, 0.0])
         solve!(solver, body_aero)
 
@@ -270,10 +303,10 @@ end
 
 @testset "solve! artificial viscosity: attached no-op, post-stall finite" begin
     body_aero = BodyAerodynamics([poststall_wing])
-    solver_off = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM,
-        is_with_artificial_viscosity=false)
-    solver_on = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM,
-        is_with_artificial_viscosity=true)
+    solver_off = Solver(poststall_wing.n_panels, poststall_wing.n_unrefined_sections;
+        solver_type=LOOP, aerodynamic_model_type=VSM, is_with_artificial_viscosity=false)
+    solver_on = Solver(poststall_wing.n_panels, poststall_wing.n_unrefined_sections;
+        solver_type=LOOP, aerodynamic_model_type=VSM, is_with_artificial_viscosity=true)
 
     # Attached flow: viscosity never fires, so results are bit-identical.
     set_va!(body_aero, [10.0, 0.0, 0.0])
@@ -290,8 +323,8 @@ end
 
 @testset "solve! reports a solve that missed the tolerances" begin
     body_aero = BodyAerodynamics([poststall_wing])
-    solver = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM,
-        max_iterations=1)
+    solver = Solver(poststall_wing.n_panels, poststall_wing.n_unrefined_sections;
+        solver_type=LOOP, aerodynamic_model_type=VSM, max_iterations=1)
     set_va!(body_aero, [10.0, 0.0, 0.0])
 
     sol = solve!(solver, body_aero)
@@ -302,7 +335,8 @@ end
     @test_throws "did not converge in 1 iterations" solve!(solver, body_aero;
         throw_on_fail=true)
 
-    converged = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM)
+    converged = Solver(poststall_wing.n_panels, poststall_wing.n_unrefined_sections;
+        solver_type=LOOP, aerodynamic_model_type=VSM)
     @test solve!(converged, body_aero; throw_on_fail=true) isa VSMSolution
 end
 
diff --git a/test/solver/test_unrefined_dist.jl b/test/solver/test_unrefined_dist.jl
index 08e955e7..6c676cc0 100644
--- a/test/solver/test_unrefined_dist.jl
+++ b/test/solver/test_unrefined_dist.jl
@@ -21,7 +21,7 @@ using Test
             wing = Wing(settings)
             refine!(wing)
             body_aero = BodyAerodynamics([wing])
-            solver = Solver(body_aero, settings)
+            solver = Solver(settings)
 
             # Set conditions and solve
             va_vec = [10.0, 0.0, 0.0]
@@ -96,7 +96,7 @@ using Test
             wing = Wing(settings)
             refine!(wing)
             body_aero = BodyAerodynamics([wing])
-            solver = Solver(body_aero, settings)
+            solver = Solver(settings)
 
             va_vec = [10.0, 0.0, 0.0]
             set_va!(body_aero, va_vec)
@@ -169,7 +169,7 @@ using Test
                 wing = Wing(settings)
                 refine!(wing)
                 body_aero = BodyAerodynamics([wing])
-                solver = Solver(body_aero, settings)
+                solver = Solver(settings)
 
                 va_vec = [10.0, 0.0, 0.0]
                 set_va!(body_aero, va_vec)
diff --git a/test/solver/test_viscous_drag_correction.jl b/test/solver/test_viscous_drag_correction.jl
index 8e8bf8e8..1358f00c 100644
--- a/test/solver/test_viscous_drag_correction.jl
+++ b/test/solver/test_viscous_drag_correction.jl
@@ -29,8 +29,8 @@ end
     refine!(wing)
     body_aero = BodyAerodynamics([wing])
 
-    solver_off = Solver(body_aero; use_gamma_prev=false)
-    solver_on = Solver(body_aero; use_gamma_prev=false,
+    solver_off = Solver(wing.n_panels, wing.n_unrefined_sections; use_gamma_prev=false)
+    solver_on = Solver(wing.n_panels, wing.n_unrefined_sections; use_gamma_prev=false,
                        is_with_viscous_drag_correction=true)
     va_vec_sideslip = V .* [cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta)]
     va_vec_straight = V .* [cos(alpha), 0.0, sin(alpha)]
@@ -96,7 +96,15 @@ end
     end
 
     @testset "defaults to off" begin
-        @test Solver(body_aero).is_with_viscous_drag_correction == false
+        @test solver_off.is_with_viscous_drag_correction == false
         @test VortexStepMethod.SolverSettings().is_with_viscous_drag_correction == false
     end
+
+    @testset "solver settings switch it on" begin
+        solver_settings = VortexStepMethod.SolverSettings(
+            is_with_viscous_drag_correction=true)
+        kwargs = VortexStepMethod.solver_kwargs(solver_settings)
+        @test Solver(wing.n_panels, wing.n_unrefined_sections;
+                     kwargs...).is_with_viscous_drag_correction
+    end
 end
diff --git a/test/verification/test_verification.jl b/test/verification/test_verification.jl
index c32d0e31..d90bdce8 100644
--- a/test/verification/test_verification.jl
+++ b/test/verification/test_verification.jl
@@ -56,8 +56,8 @@ settings of the Python verification cases.
 """
 function lift_drag_polar(wing, model, alphas; wind_speed, relaxation_factor)
     body_aero = BodyAerodynamics([wing])
-    solver = Solver(body_aero; aerodynamic_model_type=model, relaxation_factor,
-                    core_radius_fraction=1e-20)
+    solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=model,
+                    relaxation_factor, core_radius_fraction=1e-20)
     CL = zeros(length(alphas))
     CD = zeros(length(alphas))
     for (i, alpha) in enumerate(alphas)