diff --git a/CHANGELOG.md b/CHANGELOG.md
index bccd9932..ec23fb6b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,42 @@ SPDX-FileCopyrightText: 2025 Uwe Fechner, Bart van de Lint
 SPDX-License-Identifier: MPL-2.0
 -->
 
+# v0.7.0 DD-02-2026
+
+## Changed
+- BREAKING: Julia version requirement raised from 1.10 to 1.11, 1.12.
+- `reinit!()` uses a unified code path for all wing types, calling
+  `match_aero_sections_to_structure!` and
+  `compute_spatial_group_mapping!` during VSM rebuild.
+- `test_bench.jl` refactored from ad-hoc benchmarks into a proper
+  `@testset` suite with `setup_bench_sam()` helper.
+- Added `[workspace]` configuration in `Project.toml` for docs, examples,
+  scripts, and test sub-projects.
+- Manifest files renamed to `.default` suffix and gitignored.
+
+## Added
+- Asymmetric aero/structural section counts: aerodynamic and structural
+  meshes can now have different numbers of sections. When counts differ,
+  `match_aero_sections_to_structure!()` rebuilds unrefined
+  sections from structural LE/TE positions while `use_prior_polar=true`
+  preserves existing refined panel polars. Opt-in via
+  `use_prior_polar=true` on the VortexStepMethod wing.
+- `identify_wing_segments()` — identifies LE/TE pairs from groups
+  (preferred) or via a consecutive-pair heuristic.
+- `compute_spatial_group_mapping!()` — maps groups to VSM sections by
+  spatial proximity, supporting n_groups != n_aero_sections.
+- REFINE wings can now have groups (used for LE/TE pair identification).
+- QUATERNION wings can now have `wing_segments` for structural geometry
+  locking.
+- YAML loader fallback LE/TE detection in
+  `update_aero_yaml_from_struc_yaml!()` when no groups are defined
+  (consecutive-pair heuristic with x-coordinate check).
+- `test_match_aero_sections.jl` — tests geometry matching and polar
+  interpolation for both REFINE and QUATERNION wings, including
+  mismatched section counts.
+- Helper scripts: `bin/install` (environment setup, Julia version detection)
+  and `bin/run_julia` (launcher with system image support).
+
 # v0.6.1 23-02-2026
 
 ## Fixed
diff --git a/Project.toml b/Project.toml
index 2d2108f0..145a222d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "SymbolicAWEModels"
 uuid = "9c9a347c-5289-41db-a9b9-25ccc76c3360"
-version = "0.6.1"
+version = "0.7.0"
 authors = ["Bart van de Lint <bart@vandelint.net> and contributors"]
 
 [deps]
diff --git a/README.md b/README.md
index 2d9a8c0d..0cd9876f 100644
--- a/README.md
+++ b/README.md
@@ -216,6 +216,7 @@ julia --project=test test/test_point.jl
 | `test_principal_body_frame` | Principal vs body frame separation |
 | `test_heading_calculation` | Kite heading from tether geometry |
 | `test_section_alignment` | VSM section ↔ structural point mapping |
+| `test_match_aero_sections` | Asymmetric aero/structural section matching, polar interpolation |
 | `test_profile_law` | Atmospheric wind profile verification |
 | `test_bench` | Performance regression tracking |
 
diff --git a/data/2plate_kite/aero_geometry.yaml b/data/2plate_kite/aero_geometry.yaml
index 7634fb5a..ced517f8 100644
--- a/data/2plate_kite/aero_geometry.yaml
+++ b/data/2plate_kite/aero_geometry.yaml
@@ -7,13 +7,13 @@ wing_sections:
   headers: [airfoil_id, LE_x, LE_y, LE_z, TE_x, TE_y, TE_z]
   data:
     # Right section: points 2 (LE) and 3 (TE)
-    - [1, -0.5, 1.0, 2.0, 0.5, 1.0, 2.2]
+    - [1, -0.5, 1.0, 2.0, 0.5, 1.0, 2.3]
 
     # Center section: points 4 (LE) and 5 (TE)
-    - [1, -0.5, 0.0, 2.5, 0.5, 0.0, 2.7]
+    - [1, -0.5, 0.0, 2.5, 0.5, 0.0, 2.8]
 
     # Left section: points 6 (LE) and 7 (TE)
-    - [1, -0.5, -1.0, 2.0, 0.5, -1.0, 2.2]
+    - [1, -0.5, -1.0, 2.0, 0.5, -1.0, 2.3]
 
 wing_airfoils:
   alpha_range: [-180, 180, 1]  # deg
diff --git a/docs/Project.toml b/docs/Project.toml
index 29d24313..810cce87 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -11,3 +11,7 @@ VortexStepMethod = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9"
 
 [compat]
 Documenter = "1.11"
+
+[sources]
+SymbolicAWEModels = {path = ".."}
+
diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index 82fabb07..d0b65537 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -133,6 +133,8 @@ SymbolicAWEModels.resolve_ref_spec
 SymbolicAWEModels.validate_sys_struct
 SymbolicAWEModels.build_name_dict
 SymbolicAWEModels.identify_wing_segments
+SymbolicAWEModels.match_aero_sections_to_structure!
+SymbolicAWEModels.compute_spatial_group_mapping!
 SymbolicAWEModels.init_untransformed_components!
 SymbolicAWEModels.adjust_vsm_panels_to_origin!
 SymbolicAWEModels.apply_aero_z_offset!
diff --git a/docs/src/tutorial_julia.md b/docs/src/tutorial_julia.md
index 7bdf6a01..f60eaca2 100644
--- a/docs/src/tutorial_julia.md
+++ b/docs/src/tutorial_julia.md
@@ -2,6 +2,9 @@
 EditURL = "literate/tutorial_julia.jl"
 ```
 
+Copyright (c) 2025 Bart van de Lint, Jelle Poland
+SPDX-License-Identifier: MPL-2.0
+
 ```@meta
 CurrentModule = SymbolicAWEModels
 ```
diff --git a/docs/src/vsm_coupling.md b/docs/src/vsm_coupling.md
index 70222154..d732bead 100644
--- a/docs/src/vsm_coupling.md
+++ b/docs/src/vsm_coupling.md
@@ -1,28 +1,48 @@
 # VSM coupling
 
-This document explains how SymbolicAWEModels couples with the Vortex Step Method (VSM) for aerodynamic force computation. The coupling differs between the two wing types: [`QUATERNION` and `REFINE`](@ref WingType).
+This document explains how SymbolicAWEModels couples with the
+Vortex Step Method (VSM) for aerodynamic force computation. The
+coupling is configured by two orthogonal choices:
+[`WingType`](@ref) (structural representation) and
+[`AeroMode`](@ref) (force computation strategy).
 
 ## Overview
 
-Both wing types use **unrefined sections** as the fundamental structural element that maps to VSM geometry:
+Both wing types use **unrefined sections** as the fundamental
+element that maps to VSM geometry:
 
-- **Unrefined section**: A structural element defined by two points (leading edge and trailing edge) along the wing span
-- **Refined panel**: VSM can subdivide unrefined sections into multiple panels for higher aerodynamic fidelity
-- **`refined_panel_mapping`**: Maps each refined panel back to its parent unrefined section
+- **Unrefined section**: A structural element defined by two
+  points (leading edge and trailing edge) along the wing span
+- **Refined panel**: VSM can subdivide unrefined sections into
+  multiple panels for higher aerodynamic fidelity
+- **`refined_panel_mapping`**: Maps each refined panel back to
+  its parent unrefined section
 
-The VSM solver computes aerodynamic coefficients (cl, cd, cm, alpha) at both refined and unrefined levels. SymbolicAWEModels uses the **unrefined-level coefficients** to map forces back to the structural model.
+The VSM solver computes aerodynamic coefficients (cl, cd, cm,
+alpha) at both refined and unrefined levels. SymbolicAWEModels
+uses the **unrefined-level coefficients** to map forces back to
+the structural model.
 
 ## Wing types
 
-### REFINE wing type
+[`WingType`](@ref) controls the **structural representation** of
+the wing — how it deforms and how forces are distributed to
+structural degrees of freedom.
 
-The `REFINE` wing type creates the most direct coupling between structure and aerodynamics.
+### REFINE
+
+The `REFINE` wing type creates the most direct coupling between
+structure and aerodynamics.
 
 #### Structural model
 
-- Wing structure consists of [`WING`](@ref DynamicsType)-type points organized in leading edge (LE) and trailing edge (TE) pairs
-- Each consecutive pair (point i, point i+1) forms a structural segment (strut)
-- Points can move independently - the wing can deform structurally
+- Wing structure consists of [`WING`](@ref DynamicsType)-type
+  points organized in leading edge (LE) and trailing edge (TE)
+  pairs
+- Each consecutive pair (point i, point i+1) forms a structural
+  segment (strut)
+- Points can move independently — the wing can deform
+  structurally
 - Number of structural segments = (number of WING points) / 2
 
 #### VSM mapping
@@ -35,20 +55,16 @@ Structural points:  [LE₁, TE₁]  [LE₂, TE₂]  [LE₃, TE₃]
 Unrefined sections:   Sec₁         Sec₂         Sec₃
 ```
 
-Each VSM section is defined by its LE and TE point positions, taken directly from the structural point positions. VSM can subdivide these into refined panels for higher fidelity; `refined_panel_mapping` maps each refined panel back to its parent section.
-
-#### Force distribution
-
-Forces are computed and distributed in several steps:
-
-1. **VSM solve**: Computes aerodynamic coefficients and per-panel forces/moments in body frame
-2. **Panel → section mapping**: Each refined panel maps to its parent unrefined section via `refined_panel_mapping`
-3. **Moment-preserving LE/TE split** (`compute_aerostruc_loads`): Each panel's force and moment are split into LE and TE contributions that preserve the total moment about a reference point
-4. **Accumulation at structural points**: LE/TE forces are accumulated at the corresponding structural points via the `point_to_vsm_point` mapping
+Each VSM section is defined by its LE and TE point positions,
+taken directly from the structural point positions. VSM can
+subdivide these into refined panels for higher fidelity;
+`refined_panel_mapping` maps each refined panel back to its
+parent section.
 
 #### Geometry update
 
-Each timestep, structural point positions update VSM section geometry:
+Each timestep, structural point positions update VSM section
+geometry:
 
 ```julia
 # For each structural point mapped to a VSM section point:
@@ -56,22 +72,43 @@ pos_b = R_b_to_w' * (point.pos_w - wing.origin)
 section.LE_point = pos_b  # or section.TE_point
 ```
 
-This bidirectional coupling allows structural deformation to affect aerodynamics.
+This bidirectional coupling allows structural deformation to
+affect aerodynamics.
+
+#### Force distribution
+
+Per-panel forces are distributed to structural points:
+
+1. **Panel → section mapping**: Each refined panel maps to its
+   parent unrefined section via `refined_panel_mapping`
+2. **Moment-preserving LE/TE split** (`compute_aerostruc_loads`):
+   Each panel's force and moment are split into LE and TE
+   contributions that preserve the total moment about a reference
+   point
+3. **Accumulation at structural points**: LE/TE forces are
+   accumulated at the corresponding structural points via the
+   `point_to_vsm_point` mapping
 
-### QUATERNION wing type
+### QUATERNION
 
-The `QUATERNION` wing type uses a rigid body representation with optional deformable groups.
+The `QUATERNION` wing type uses a rigid body representation with
+optional deformable groups.
 
 #### Structural model
 
 - Wing treated as a rigid body with quaternion-based orientation
-- No per-point wing structure - aerodynamic forces applied to wing center of mass
-- Optional **groups** represent deformable sections with twist degrees of freedom
-- Groups control segment twist angles but don't affect primary aerodynamic geometry
+- No per-point wing structure — aerodynamic forces applied to
+  wing center of mass
+- Optional **groups** represent deformable sections with twist
+  degrees of freedom
+- Groups control segment twist angles. With
+  `use_prior_polar=true`, group LE/TE positions also define the
+  aerodynamic section geometry
 
 #### VSM mapping
 
-VSM still uses unrefined sections, but they don't correspond to individual structural points:
+VSM still uses unrefined sections, but they don't correspond to
+individual structural points:
 
 ```
 Unrefined sections: [Sec₁,  Sec₂,  Sec₃,  Sec₄,  Sec₅]
@@ -80,42 +117,129 @@ Groups:              [─ Group₁ ─]    [─ Group₂ ─]
                      twist DOF θ₁    twist DOF θ₂
 ```
 
-Multiple unrefined sections can be combined into a single group for twist control. The mapping is configured via:
+Multiple unrefined sections can be combined into a single group
+for twist control. The mapping is configured via:
 
 ```julia
 group.unrefined_section_idxs = [start_idx:end_idx]
 ```
 
-#### Force computation
-
-VSM linearization returns integrated aerodynamic coefficients per wing:
-
-1. **VSM linearize**: Computes baseline coefficients and Jacobian around the current operating point
-   - Output state `vsm_x = [C_F(3), C_M(3), section_moments(n_unrefined)]`
-   - Input state `vsm_y = [va_b(3), twist_angles(n_unrefined), ω_b(3)]`
-2. **Symbolic linearization**: The ODE uses `F = q∞ · A · (C₀ + J · Δstate)` where `Δstate = state - state₀`
-   - `C₀[1:3]` → total force coefficient, `C₀[4:6]` → total moment coefficient
-   - `C₀[7:end]` → per-section twist moment coefficients, summed per group
-3. **Rigid body dynamics**: Integrated force/moment drive quaternion dynamics
+#### Force distribution
 
-Groups affect twist deformation applied to VSM sections before the solve. Each group's aerodynamic moment is the sum of its unrefined section moments, driving the twist DOF.
+Integrated force and moment coefficients are applied to the rigid
+body, driving quaternion dynamics. Each group's aerodynamic moment
+is the sum of its unrefined section moments, driving the twist
+DOF.
+
+## Aero modes
+
+[`AeroMode`](@ref) controls **how aerodynamic forces enter the
+ODE system** — orthogonal to the wing type choice.
+
+### AERO_DIRECT
+
+The VSM solver runs a full nonlinear solve. The resulting forces
+are stored in the wing struct and read by registered symbolic
+functions during ODE evaluation:
+
+1. `VortexStepMethod.linearize()` computes baseline coefficients
+2. Physical forces are computed: `F = q∞ · A · C₀`
+3. Forces are stored in `wing.aero_force_b` and
+   `wing.aero_moment_b` (QUATERNION) or per-point via
+   `distribute_panel_forces_to_points!` (REFINE)
+4. Between VSM updates (controlled by `vsm_interval`), forces
+   are held constant
+
+For QUATERNION wings, the symbolic equations read the stored
+forces via `get_aero_force_override` / `get_aero_moment_override`.
+For REFINE wings, per-point forces are read via
+`get_point_aero_force`.
+
+### AERO_LINEARIZED
+
+A first-order Taylor expansion using the Jacobian from VSM
+linearization. This enables the ODE solver to see smooth
+force variations between VSM updates:
+
+1. `VortexStepMethod.linearize()` computes the Jacobian and
+   baseline coefficients around the current operating point
+   - Output state
+     `vsm_x = [C_F(3), C_M(3), section_moments(n_unrefined)]`
+   - Input state
+     `vsm_y = [va_b(3), twist_angles(n_unrefined), ω_b(3)]`
+2. The symbolic ODE uses
+   `F = q∞ · A · (C₀ + J · Δstate)` where
+   `Δstate = state - state₀`
+   - `C₀[1:3]` → total force coefficient,
+     `C₀[4:6]` → total moment coefficient
+   - `C₀[7:end]` → per-section twist moment coefficients,
+     summed per group
+3. Between VSM updates, the Jacobian extrapolates forces as
+   the state evolves
+
+### AERO_NONE
+
+Returns zero forces. Useful for debugging rigid body dynamics
+without aerodynamic coupling.
+
+### Compatibility
+
+| Wing type | Default aero mode | Supported modes |
+|-----------|-------------------|-----------------|
+| **QUATERNION** | `AERO_LINEARIZED` | `AERO_LINEARIZED`, `AERO_DIRECT`, `AERO_NONE` |
+| **REFINE** | `AERO_DIRECT` | `AERO_DIRECT`, `AERO_NONE` |
+
+`REFINE` + `AERO_LINEARIZED` is not yet implemented (raises an
+error at runtime).
+
+## Aligning aero sections to structure
+
+When the number of aerodynamic sections differs from the number
+of structural LE/TE pairs, `match_aero_sections_to_structure!`
+rebuilds the unrefined sections so their geometry matches the
+structure. This applies to both wing types and requires
+`use_prior_polar=true` on the VortexStepMethod wing.
+
+The steps are:
+
+1. **Find structural LE/TE pairs**: `identify_wing_segments`
+   extracts pairs from groups (preferred) or uses a
+   consecutive-pair heuristic
+2. **Rebuild unrefined sections**: For each structural pair, a
+   new `Section` is created with LE/TE positions from the
+   structural points (in body frame). Its airfoil data
+   (`aero_model`, `aero_data`) is copied from the nearest
+   original unrefined section by span index
+3. **Re-refine**: `refine!` updates refined panel geometry from
+   the rebuilt unrefined sections. Because `use_prior_polar=true`
+   and `n_panels` is unchanged, existing refined panel polars are
+   preserved — only positions are re-interpolated
+4. **Resize linearization state**: For non-REFINE wings, `vsm_y`,
+   `vsm_x`, and `vsm_jac` are resized to match the new section
+   count
 
 ## Refined panel mapping
 
-Both wing types use `refined_panel_mapping` to handle VSM mesh refinement:
+Both wing types use `refined_panel_mapping` to handle VSM mesh
+refinement:
 
 ### Purpose
 
-VSM can subdivide unrefined sections into multiple refined panels for higher aerodynamic fidelity. The mapping tracks which parent unrefined section each refined panel belongs to.
+VSM can subdivide unrefined sections into multiple refined panels
+for higher aerodynamic fidelity. The mapping tracks which parent
+unrefined section each refined panel belongs to.
 
 ### Computation
 
-After VSM refinement, `compute_refined_panel_mapping!` finds the closest unrefined section for each refined panel by comparing center positions:
+After VSM refinement, `compute_refined_panel_mapping!` finds the
+closest unrefined section for each refined panel by comparing
+center positions:
 
 ```julia
 for each refined_panel in wing.refined_sections
     center = compute_center(refined_panel)
-    closest_section = argmin(distance(center, unrefined_section_centers))
+    closest_section = argmin(
+        distance(center, unrefined_section_centers))
     refined_panel_mapping[refined_panel_idx] = closest_section
 end
 ```
@@ -124,27 +248,36 @@ end
 
 The mapping enables:
 
-1. **Group twist angles**: Applying the correct twist angle from groups to refined panels via their parent section
-2. **Force distribution (REFINE)**: Accumulating refined panel forces at the structural points of their parent section
-3. **Linearization (QUATERNION)**: Propagating state perturbations through the correct sections
+1. **Group twist angles**: Applying the correct twist angle from
+   groups to refined panels via their parent section
+2. **Force distribution (REFINE)**: Accumulating refined panel
+   forces at the structural points of their parent section
+3. **Linearization (QUATERNION + AERO_LINEARIZED)**: Propagating
+   state perturbations through the correct sections
 
-## Key differences summary
+## Wing type summary
 
 | Aspect | REFINE | QUATERNION |
 |--------|--------|------------|
-| **Structural representation** | Individual WING points | Rigid body with quaternion |
-| **Unrefined section count** | = number of structural LE/TE pairs | Independent of structure |
-| **Force distribution** | Per-point via moment-preserving LE/TE split | Integrated force/moment coefficients |
-| **Deformation coupling** | Direct: point motion → VSM geometry | Indirect: group twists → VSM sections |
-| **Computational cost** | Higher (full VSM solve per step) | Lower (linearized) |
-| **Fidelity** | Higher (aeroelastic coupling) | Lower (rigid body) |
+| **Structural repr.** | Individual WING points | Rigid body + quaternion |
+| **Section count** | = structural LE/TE pairs | Independent; optionally rebuilt via `use_prior_polar` |
+| **Force distribution** | Per-point moment-preserving LE/TE split | Integrated force/moment on body |
+| **Deformation** | Direct: point motion → VSM geometry | Indirect: group twists → sections |
+| **Default aero mode** | `AERO_DIRECT` | `AERO_LINEARIZED` |
 
 ## Implementation files
 
-- `src/vsm_refine.jl`: REFINE wing force distribution and geometry updates
-- `src/system_structure/types.jl`: Component type definitions (Point, Segment, etc.)
-- `src/system_structure/wing.jl`: Wing and VSMWing type definitions, group-to-section mapping
-- `src/generate_system/vsm_eqs.jl`: Symbolic VSM equation generation
-- `src/generate_system/wing_eqs.jl`: Wing dynamics equation generation
-- `src/linearize.jl`: VSM linearization and Jacobian updates
-- VortexStepMethod.jl `src/wing_geometry.jl`: `refined_panel_mapping` computation
+- `src/vsm_refine.jl`: Aero-to-structure alignment (all wing
+  types), REFINE force distribution, and geometry updates
+- `src/system_structure/types.jl`: Component type definitions
+  including `WingType` and `AeroMode` enums
+- `src/system_structure/wing.jl`: Wing and VSMWing type
+  definitions, group-to-section mapping
+- `src/generate_system/vsm_eqs.jl`: Symbolic VSM equation
+  generation (all wing type × aero mode combinations)
+- `src/generate_system/wing_eqs.jl`: Wing dynamics equation
+  generation
+- `src/linearize.jl`: VSM update dispatch — linearization
+  (QUATERNION) and nonlinear solve (REFINE)
+- VortexStepMethod.jl `src/wing_geometry.jl`:
+  `refined_panel_mapping` computation
diff --git a/examples/coupled_2plate_kite.jl b/examples/coupled_2plate_kite.jl
index af6d37ff..1664854d 100644
--- a/examples/coupled_2plate_kite.jl
+++ b/examples/coupled_2plate_kite.jl
@@ -20,7 +20,7 @@ pkg_root = dirname(@__DIR__)
 set_data_path(joinpath(pkg_root, "data", MODEL_NAME))
 
 struc_yaml = joinpath(get_data_path(),
-                      "quat_struc_geometry.yaml")
+                      "refine_struc_geometry.yaml")
 aero_yaml = joinpath(get_data_path(), "aero_geometry.yaml")
 update_aero_yaml_from_struc_yaml!(struc_yaml, aero_yaml)
 
diff --git a/src/system_structure/system_structure_core.jl b/src/system_structure/system_structure_core.jl
index ab7a38c9..25285e69 100644
--- a/src/system_structure/system_structure_core.jl
+++ b/src/system_structure/system_structure_core.jl
@@ -398,6 +398,82 @@ function calc_inertia_y_rotation(I_tensor)
     return I_diag, Ry
 end
 
+"""
+    compute_spatial_group_mapping!(the_wing, groups, points)
+
+Map groups to unrefined sections using spatial proximity.
+Each group is assigned to the closest unrefined section
+based on distance between centers.
+"""
+function compute_spatial_group_mapping!(
+    the_wing::VSMWing,
+    groups::AbstractVector{Group},
+    points::AbstractVector{Point}
+)
+    the_vsm_wing = the_wing.vsm_wing
+    n_unrefined = the_vsm_wing.n_unrefined_sections
+    n_groups = length(the_wing.base.group_idxs)
+
+    # Compute group centers in body frame
+    group_centers = Vector{MVec3}(undef, n_groups)
+    for (local_idx, group_idx) in
+            enumerate(the_wing.base.group_idxs)
+        group = groups[group_idx]
+        center = zeros(3)
+        for pt_idx in group.point_idxs
+            center += the_wing.base.R_b_to_c' *
+                (points[pt_idx].pos_cad -
+                 the_wing.base.pos_cad)
+        end
+        group_centers[local_idx] =
+            center / length(group.point_idxs)
+    end
+
+    # Compute unrefined section centers
+    unrefined_centers = Vector{MVec3}(
+        undef, n_unrefined)
+    for i in 1:n_unrefined
+        le_point =
+            the_vsm_wing.unrefined_sections[i].LE_point
+        te_point =
+            the_vsm_wing.unrefined_sections[i].TE_point
+        unrefined_centers[i] =
+            (le_point + te_point) / 2
+    end
+
+    # Map each group to closest unrefined section
+    for (local_idx, group_idx) in
+            enumerate(the_wing.base.group_idxs)
+        group = groups[group_idx]
+        min_dist = Inf
+        closest_idx = 1
+        for unrefined_idx in 1:n_unrefined
+            dist = norm(
+                group_centers[local_idx] -
+                unrefined_centers[unrefined_idx])
+            if dist < min_dist
+                min_dist = dist
+                closest_idx = unrefined_idx
+            end
+        end
+        group.unrefined_section_idxs =
+            Int64[closest_idx]
+    end
+
+    # Validate: check all sections are covered
+    assigned = Set{Int64}()
+    for group_idx in the_wing.base.group_idxs
+        union!(assigned,
+            groups[group_idx].unrefined_section_idxs)
+    end
+    if length(assigned) != n_unrefined
+        unassigned = setdiff(1:n_unrefined, assigned)
+        @warn "Wing $(the_wing.base.idx): " *
+            "$(length(unassigned)) unrefined sections " *
+            "not assigned to any group: $unassigned"
+    end
+end
+
 # ==================== CONSTRUCTOR ==================== #
 
 """
@@ -459,14 +535,6 @@ function SystemStructure(name, set;
      tether_names_dict, winch_names_dict, wing_names_dict, transform_names_dict) =
         assign_indices_and_resolve!(points, groups, segments, pulleys, tethers, winches, wings, transforms)
 
-    # For REFINE wings, clear group_idxs - groups are not used with REFINE
-    # (REFINE applies panel forces directly to structural points)
-    for wing in wings
-        if wing.wing_type == REFINE && !isempty(wing.group_idxs)
-            empty!(wing.group_idxs)
-        end
-    end
-
     # If no wings defined, convert WING points to STATIC
     if isempty(wings)
         wing_point_idxs = findall(p -> p.type == WING, points)
@@ -682,12 +750,14 @@ function SystemStructure(name, set;
                 # Use integer as name for auto-created groups
                 group_name = group_idx
 
-                # Only TE point, moment_frac=0.0 (rotate around LE)
-                new_group = Group(group_name, [te_idx], DYNAMIC, 0.0)
+                # Both LE and TE points (matches YAML convention)
+                new_group = Group(group_name,
+                    [le_idx, te_idx], DYNAMIC, 0.0)
 
-                # Assign idx and resolve point_refs since these are dynamically created
+                # Assign idx and resolve point_refs since
+                # these are dynamically created
                 new_group.idx = group_idx
-                new_group.point_idxs = [te_idx]
+                new_group.point_idxs = [le_idx, te_idx]
 
                 push!(groups, new_group)
                 push!(new_group_idxs, Int64(group_idx))
@@ -709,59 +779,24 @@ function SystemStructure(name, set;
         end
     end
 
-    """
-        compute_spatial_group_mapping!(the_wing, groups, points)
-
-    Map groups to unrefined sections using spatial proximity.
-    Each group is assigned to the closest unrefined section based on distance between centers.
-    """
-    function compute_spatial_group_mapping!(the_wing::VSMWing, groups::AbstractVector{Group}, points::AbstractVector{Point})
-        the_vsm_wing = the_wing.vsm_wing
-        n_unrefined = the_vsm_wing.n_unrefined_sections
-        n_groups = length(the_wing.base.group_idxs)
-
-        # Compute group centers in body frame (average of all attach points)
-        group_centers = Vector{MVec3}(undef, n_groups)
-        for (local_idx, group_idx) in enumerate(the_wing.base.group_idxs)
-            group = groups[group_idx]
-            center = zeros(3)
-            for pt_idx in group.point_idxs
-                center += the_wing.base.R_b_to_c' * (points[pt_idx].pos_cad - the_wing.base.pos_cad)
-            end
-            group_centers[local_idx] = center / length(group.point_idxs)
-        end
-
-        # Compute unrefined section centers
-        unrefined_centers = Vector{MVec3}(undef, n_unrefined)
-        for i in 1:n_unrefined
-            le_point = the_vsm_wing.unrefined_sections[i].LE_point
-            te_point = the_vsm_wing.unrefined_sections[i].TE_point
-            unrefined_centers[i] = (le_point + te_point) / 2
-        end
-
-        # Map each group to closest unrefined section
-        for (local_idx, group_idx) in enumerate(the_wing.base.group_idxs)
-            group = groups[group_idx]
-            min_dist = Inf
-            closest_idx = 1
-            for unrefined_idx in 1:n_unrefined
-                dist = norm(group_centers[local_idx] - unrefined_centers[unrefined_idx])
-                if dist < min_dist
-                    min_dist = dist
-                    closest_idx = unrefined_idx
-                end
-            end
-            group.unrefined_section_idxs = Int64[closest_idx]
-        end
+    # Match aero sections to structural LE/TE for ALL
+    # VSMWing types (runs after auto-group creation so
+    # identify_wing_segments can use groups).
+    for wing in wings
+        isa(wing, VSMWing) || continue
+        wing.aero_mode == AERO_NONE && continue
+        match_aero_sections_to_structure!(
+            wing, points; groups=groups)
+    end
 
-        # Validate: check all sections are covered
-        assigned = Set{Int64}()
-        for group_idx in the_wing.base.group_idxs
-            union!(assigned, groups[group_idx].unrefined_section_idxs)
-        end
-        if length(assigned) != n_unrefined
-            unassigned = setdiff(1:n_unrefined, assigned)
-            @warn "Wing $(the_wing.base.idx): $(length(unassigned)) unrefined sections not assigned to any group: $unassigned"
+    # Clear REFINE wing.group_idxs — groups were used
+    # for LE/TE identification but REFINE doesn't use
+    # them for aerodynamics.  Groups stay in sys_struct
+    # (useful for structural info / future linearization).
+    for wing in wings
+        if wing.wing_type == REFINE &&
+           !isempty(wing.group_idxs)
+            empty!(wing.group_idxs)
         end
     end
 
@@ -858,7 +893,9 @@ function SystemStructure(name, set;
             # Identify wing segments (LE/TE pairs)
             if isnothing(wing.wing_segments)
                 wing.wing_segments =
-                    identify_wing_segments(wing_points)
+                    identify_wing_segments(
+                        wing_points; groups=groups,
+                        wing_group_idxs=wing.group_idxs)
             end
 
             # Set default reference points if not provided
diff --git a/src/system_structure/utilities.jl b/src/system_structure/utilities.jl
index 27ad3992..8592db95 100644
--- a/src/system_structure/utilities.jl
+++ b/src/system_structure/utilities.jl
@@ -369,35 +369,39 @@ function reinit!(sys_struct::SystemStructure, set::Settings;
             # Transform sections: CAD → body frame
             # (must match SystemStructure constructor)
             vsm_wing = wing.vsm_wing
-            if wing.wing_type == QUATERNION
-                vsm_wing.T_cad_body .= wing.pos_cad
-                adjust_vsm_panels_to_origin!(
-                    vsm_wing, wing.pos_cad)
-                rotate_vsm_sections!(
-                    vsm_wing, wing.R_b_to_c')
-                vsm_wing.R_cad_body .= wing.R_b_to_c
+            vsm_wing.T_cad_body .= wing.pos_cad
+            adjust_vsm_panels_to_origin!(
+                vsm_wing, wing.pos_cad)
+            rotate_vsm_sections!(
+                vsm_wing, wing.R_b_to_c')
+            vsm_wing.R_cad_body .= wing.R_b_to_c
+            if wing.wing_type != REFINE
                 apply_aero_z_offset!(
                     vsm_wing, wing.aero_z_offset)
-                VortexStepMethod.reinit!(
-                    wing.vsm_aero)
-            elseif wing.wing_type == REFINE
-                vsm_wing.T_cad_body .= wing.pos_cad
-                adjust_vsm_panels_to_origin!(
-                    vsm_wing, wing.pos_cad)
-                rotate_vsm_sections!(
-                    vsm_wing, wing.R_b_to_c')
-                vsm_wing.R_cad_body .= wing.R_b_to_c
-                VortexStepMethod.reinit!(
-                    wing.vsm_aero)
-                if !isnothing(wing.point_to_vsm_point)
-                    wing_point_idxs = collect(
-                        keys(wing.point_to_vsm_point))
-                    wing_points = [points[idx]
-                        for idx in wing_point_idxs]
-                    wing.point_to_vsm_point =
-                        build_point_to_vsm_point_mapping(
-                            wing_points, vsm_wing)
-                end
+            end
+            VortexStepMethod.reinit!(wing.vsm_aero)
+
+            # Match aero sections to structure (all types)
+            match_aero_sections_to_structure!(
+                wing, points; groups=groups)
+
+            # Recompute group→section mapping
+            if wing.wing_type == QUATERNION &&
+               !isempty(wing.group_idxs)
+                compute_spatial_group_mapping!(
+                    wing, groups, points)
+            end
+
+            # REFINE-only: rebuild point mapping
+            if wing.wing_type == REFINE &&
+               !isnothing(wing.point_to_vsm_point)
+                wing_point_idxs = collect(
+                    keys(wing.point_to_vsm_point))
+                wing_pts = [points[idx]
+                    for idx in wing_point_idxs]
+                wing.point_to_vsm_point =
+                    build_point_to_vsm_point_mapping(
+                        wing_pts, vsm_wing)
             end
         end
     end
diff --git a/src/system_structure/wing.jl b/src/system_structure/wing.jl
index 6fdb4908..49f51a6a 100644
--- a/src/system_structure/wing.jl
+++ b/src/system_structure/wing.jl
@@ -316,7 +316,8 @@ Creates vsm_wing, vsm_aero, and vsm_solver internally.
 - `y_damping::SimFloat=150.0`: Lateral damping coefficient.
 - `wing_type::WingType=QUATERNION`: Aerodynamic model type.
 - `point_to_vsm_point`: 1:1 structural point to VSM point mapping (REFINE only).
-- `wing_segments`: LE/TE pairs (REFINE only).
+- `wing_segments`: LE/TE pairs (populated for all VSM wing types by
+  `match_aero_sections_to_structure!`).
 - `z_ref_points`: Chord direction reference points (REFINE only, names or indices).
 - `y_ref_points`: Span direction reference points (REFINE only, names or indices).
 - `origin`: Reference to origin point (REFINE only, name or index).
@@ -346,8 +347,6 @@ function VSMWing(name, set::Settings,
 
     # Validation
     if wing_type == REFINE
-        @assert length(groups) == 0
-            "REFINE wings cannot have groups"
         @assert !isnothing(origin)
             "REFINE wings require origin to define KCU position"
         if !isnothing(pos_cad)
@@ -359,8 +358,6 @@ function VSMWing(name, set::Settings,
     else
         @assert isnothing(point_to_vsm_point)
             "QUATERNION wings: no point_to_vsm_point"
-        @assert isnothing(wing_segments)
-            "QUATERNION wings: no wing_segments"
         # origin, z_ref_points, y_ref_points are now
         # accepted for QUATERNION wings (body frame
         # defined by structural ref points)
diff --git a/src/vsm_refine.jl b/src/vsm_refine.jl
index c85df302..c4f59ff3 100644
--- a/src/vsm_refine.jl
+++ b/src/vsm_refine.jl
@@ -2,8 +2,11 @@
 # SPDX-License-Identifier: MPL-2.0
 
 """
-Helper functions for REFINE wing type that applies VSM panel forces directly to
-structural points.
+Helper functions for VSM wing types.
+
+REFINE-specific functions (panel force distribution, structural geometry
+updates) are at the bottom of this file.  The shared
+`match_aero_sections_to_structure!` works for all VSMWing types.
 """
 
 # Baseline chord-based aero scaling for REFINE wings.
@@ -11,37 +14,81 @@ structural points.
 const AERO_SCALE_CHORD = 0.0
 
 """
-    identify_wing_segments(wing_points::AbstractVector{Point})
+    identify_wing_segments(wing_points; groups=nothing, wing_group_idxs=nothing)
 
 Identify wing segments (LE/TE pairs) from WING-type points.
 
-Assumes points are organized in pairs along the span, with even-numbered points
-being leading edge and odd-numbered being trailing edge (or vice versa).
+When `groups` and `wing_group_idxs` are provided, uses group `point_idxs`
+to determine LE (`point_idxs[1]`) and TE (`point_idxs[end]`) for each
+section. Falls back to a consecutive-pair heuristic (sorted by point index)
+when groups are unavailable.
+
+In both paths an x-coordinate check swaps LE/TE if needed (LE has
+smaller `pos_cad[1]`).
 
 # Arguments
-- `wing_points::AbstractVector{Point}`: All WING-type points for a wing (sorted by index)
+- `wing_points::AbstractVector{Point}`: WING-type points for a wing.
+
+# Keyword Arguments
+- `groups::Union{Nothing, AbstractVector{Group}}`: All groups in the
+  system (indexed by `wing_group_idxs`).
+- `wing_group_idxs::Union{Nothing, AbstractVector{<:Integer}}`:
+  Indices into `groups` belonging to this wing.
 
 # Returns
-- `Vector{Tuple{Int64, Int64}}`: Vector of (le_point_idx, te_point_idx) pairs defining segments
+- `Vector{Tuple{Int64, Int64}}`: (le_point_idx, te_point_idx) pairs.
 """
-function identify_wing_segments(wing_points::AbstractVector{Point})
-    # Sort points by index to ensure consistent ordering
+function identify_wing_segments(
+    wing_points::AbstractVector{Point};
+    groups::Union{Nothing, AbstractVector{Group}}=nothing,
+    wing_group_idxs::Union{Nothing,
+        AbstractVector{<:Integer}}=nothing
+)
+    use_groups = !isnothing(groups) &&
+        !isnothing(wing_group_idxs) &&
+        !isempty(wing_group_idxs)
+
+    if use_groups
+        segments = Tuple{Int64, Int64}[]
+        for g_idx in wing_group_idxs
+            group = groups[g_idx]
+            length(group.point_idxs) >= 2 || error(
+                "Group $(group.name): need at least " *
+                "2 point_idxs (LE/TE), got " *
+                "$(length(group.point_idxs))")
+            le_idx = group.point_idxs[1]
+            te_idx = group.point_idxs[end]
+            le_point = wing_points[findfirst(
+                p -> p.idx == le_idx, wing_points)]
+            te_point = wing_points[findfirst(
+                p -> p.idx == te_idx, wing_points)]
+            # Safety: swap if LE actually has larger x
+            if le_point.pos_cad[1] < te_point.pos_cad[1]
+                push!(segments, (le_idx, te_idx))
+            else
+                push!(segments, (te_idx, le_idx))
+            end
+        end
+        return segments
+    end
+
+    # Fallback: consecutive-pair heuristic
     sorted_points = sort(wing_points, by=p->p.idx)
 
     n_points = length(sorted_points)
-    @assert n_points % 2 == 0 "Wing must have even number of points (LE/TE pairs)"
+    @assert n_points % 2 == 0 (
+        "Wing must have even number of points " *
+        "(LE/TE pairs)")
 
     n_segments = n_points ÷ 2
     segments = Tuple{Int64, Int64}[]
 
-    # Group consecutive pairs: (point[1], point[2]), (point[3], point[4]), ...
     for i in 1:n_segments
         le_idx = 2*i - 1
         te_idx = 2*i
         le_point = sorted_points[le_idx]
         te_point = sorted_points[te_idx]
 
-        # Determine which is LE and which is TE by x-coordinate (LE has smaller x)
         if le_point.pos_cad[1] < te_point.pos_cad[1]
             push!(segments, (le_point.idx, te_point.idx))
         else
@@ -52,6 +99,164 @@ function identify_wing_segments(wing_points::AbstractVector{Point})
     return segments
 end
 
+"""
+    match_aero_sections_to_structure!(wing, points; groups)
+
+Rebuild unrefined sections to match structural LE/TE positions,
+preserving refined panel polars via `use_prior_polar`.
+
+Works for **all** VSMWing types (QUATERNION and REFINE).  When the
+structural and aerodynamic section counts match the rebuild is a 1:1
+copy (`source_idx == i`) that ensures positions exactly match
+structural points.  When they differ, `use_prior_polar` and non-empty
+`refined_sections` are required.
+
+For non-REFINE wings whose section count changed, the linearization
+vectors (`vsm_y`, `vsm_x`, `vsm_jac`) are resized to match the new
+`n_unrefined_sections`.
+
+# Keyword Arguments
+- `groups::AbstractVector{Group}`: Groups in the system (used for
+  group-based LE/TE identification via [`identify_wing_segments`](@ref)).
+"""
+function match_aero_sections_to_structure!(
+    wing::VSMWing,
+    points::AbstractVector{Point};
+    groups::AbstractVector{Group}=Group[]
+)
+    wing_points = [
+        p for p in points if
+        p.type == WING && p.wing_idx == wing.idx
+    ]
+
+    wing_group_idxs = isempty(groups) ? nothing :
+        wing.group_idxs
+    grps = isempty(groups) ? nothing : groups
+    has_groups = !isnothing(grps) &&
+        !isnothing(wing_group_idxs) &&
+        !isempty(wing_group_idxs)
+
+    if has_groups
+        n_struct_sections = length(wing_group_idxs)
+        # REFINE: each group is a 2-point strut (LE/TE)
+        if wing.wing_type == REFINE
+            for g_idx in wing_group_idxs
+                g = grps[g_idx]
+                length(g.point_idxs) == 2 || error(
+                    "REFINE wing $(wing.idx): group " *
+                    "$(g.name) must have exactly 2 " *
+                    "points (LE/TE pair), got " *
+                    "$(length(g.point_idxs))")
+            end
+        end
+    else
+        n_points = length(wing_points)
+        n_points % 2 == 0 || error(
+            "Wing $(wing.idx): no groups and odd " *
+            "number of WING points " *
+            "($(n_points)). Define groups to " *
+            "specify LE/TE pairs.")
+        n_struct_sections = n_points ÷ 2
+    end
+
+    n_aero_sections =
+        length(wing.vsm_wing.unrefined_sections)
+    counts_differ = n_struct_sections != n_aero_sections
+
+    if counts_differ
+        wing.vsm_wing.use_prior_polar || error(
+            "Wing $(wing.idx): structural sections " *
+            "($(n_struct_sections)) do not match " *
+            "aerodynamic sections " *
+            "($(n_aero_sections)). Set " *
+            "use_prior_polar=true to rebuild " *
+            "unrefined sections from structural " *
+            "geometry."
+        )
+
+        isempty(wing.vsm_wing.refined_sections) &&
+            error(
+                "Wing $(wing.idx): cannot rebuild " *
+                "unrefined sections because no " *
+                "refined sections exist to " *
+                "preserve polars from."
+            )
+    end
+
+    wing_segments = identify_wing_segments(
+        wing_points; groups=grps,
+        wing_group_idxs=wing_group_idxs)
+    wing.wing_segments = wing_segments
+    length(wing_segments) == n_struct_sections || error(
+        "Wing $(wing.idx): failed to identify " *
+        "structural LE/TE pairs."
+    )
+
+    original_sections = wing.vsm_wing.unrefined_sections
+    n_original = length(original_sections)
+    n_original > 0 || error(
+        "Wing $(wing.idx): aerodynamic geometry " *
+        "has zero unrefined sections."
+    )
+    R_b_to_c = wing.R_b_to_c
+    origin_cad = wing.pos_cad
+    new_sections = Vector{VortexStepMethod.Section}(
+        undef, n_struct_sections)
+
+    for (i, (le_idx, te_idx)) in enumerate(wing_segments)
+        source_idx = if counts_differ
+            n_struct_sections == 1 ? 1 :
+                round(Int,
+                    1 + (i - 1) * (n_original - 1) /
+                    (n_struct_sections - 1))
+        else
+            i  # 1:1 copy when counts match
+        end
+        source_section = original_sections[source_idx]
+
+        le_body = R_b_to_c' *
+            (points[le_idx].pos_cad - origin_cad)
+        te_body = R_b_to_c' *
+            (points[te_idx].pos_cad - origin_cad)
+
+        section = VortexStepMethod.Section()
+        if isnothing(source_section.aero_data)
+            VortexStepMethod.reinit!(
+                section, le_body, te_body,
+                source_section.aero_model
+            )
+        else
+            VortexStepMethod.reinit!(
+                section, le_body, te_body,
+                source_section.aero_model,
+                source_section.aero_data
+            )
+        end
+        new_sections[i] = section
+    end
+
+    wing.vsm_wing.unrefined_sections = new_sections
+    wing.vsm_wing.n_unrefined_sections =
+        Int16(n_struct_sections)
+
+    refine!(wing.vsm_wing;
+        recompute_mapping=true, sort_sections=false)
+    VortexStepMethod.reinit!(wing.vsm_aero)
+
+    # Resize linearization vectors for non-REFINE wings
+    # when section count changed.
+    if counts_differ && wing.wing_type != REFINE
+        n_un = Int(n_struct_sections)
+        ny = 3 + n_un + 3
+        nx = 3 + 3 + n_un
+        wing.vsm_y = zeros(SimFloat, ny)
+        wing.vsm_x = zeros(SimFloat, nx)
+        wing.vsm_jac = zeros(SimFloat, nx, ny)
+    end
+
+    return nothing
+end
+
 """
     build_point_to_vsm_point_mapping(wing_points::AbstractVector{Point}, vsm_wing::VortexStepMethod.AbstractWing)
 
diff --git a/src/yaml_loader.jl b/src/yaml_loader.jl
index f3e49ea8..cd76e0e1 100644
--- a/src/yaml_loader.jl
+++ b/src/yaml_loader.jl
@@ -1066,6 +1066,7 @@ function update_aero_yaml_from_struc_yaml!(source_struc_yaml::AbstractString,
     # Parse struc YAML to extract WING point positions and group LE/TE pairs
     struc_lines = readlines(struc_full_path)
     wing_pos_dict = Dict{String, Vector{Float64}}()  # name => [x, y, z]
+    wing_names_ordered = String[]                     # preserve YAML order
     group_le_te = Vector{Tuple{String, String}}()     # (le_name, te_name)
 
     current_section = :none
@@ -1097,10 +1098,12 @@ function update_aero_yaml_from_struc_yaml!(source_struc_yaml::AbstractString,
             # Format: - [name, [x, y, z], TYPE, ...]
             m = match(r"^\s*-\s*\[(\w+)\s*,\s*\[([-+]?\d+\.?\d*(?:[eE][-+]?\d+)?)\s*,\s*([-+]?\d+\.?\d*(?:[eE][-+]?\d+)?)\s*,\s*([-+]?\d+\.?\d*(?:[eE][-+]?\d+)?)\s*\]\s*,\s*(\w+)", line)
             if m !== nothing && m.captures[5] == "WING"
-                wing_pos_dict[m.captures[1]] = [
+                name = m.captures[1]
+                wing_pos_dict[name] = [
                     parse(Float64, m.captures[2]),
                     parse(Float64, m.captures[3]),
                     parse(Float64, m.captures[4])]
+                push!(wing_names_ordered, name)
             end
         elseif current_section == :groups
             # Format: - [name, [pt1, pt2, ...], ...]
@@ -1114,10 +1117,26 @@ function update_aero_yaml_from_struc_yaml!(source_struc_yaml::AbstractString,
         end
     end
 
-    # Validate
+    # Fallback: consecutive-pair heuristic when no groups
     if isempty(group_le_te)
-        error("No groups with point_idxs found in $struc_full_path")
+        n_wp = length(wing_names_ordered)
+        n_wp % 2 == 0 || error(
+            "No groups in $struc_full_path and odd " *
+            "number of WING points ($n_wp). " *
+            "Define groups to specify LE/TE pairs.")
+        for i in 1:2:n_wp
+            a = wing_names_ordered[i]
+            b = wing_names_ordered[i + 1]
+            # LE has smaller x
+            if wing_pos_dict[a][1] <
+                    wing_pos_dict[b][1]
+                push!(group_le_te, (a, b))
+            else
+                push!(group_le_te, (b, a))
+            end
+        end
     end
+
     for (le, te) in group_le_te
         haskey(wing_pos_dict, le) || error(
             "Group LE point '$le' not found in WING points")
diff --git a/test/test_bench.jl b/test/test_bench.jl
index 8fc528b9..cfd82c11 100644
--- a/test/test_bench.jl
+++ b/test/test_bench.jl
@@ -1,11 +1,13 @@
 # SPDX-FileCopyrightText: 2025 Bart van de Lint
 # SPDX-License-Identifier: MPL-2.0
 
-# Benchmark ODE RHS for 2plate_kite model.
-# Tests: raw RHS, registered functions directly, allocation profile.
+# Benchmark tests for ODE RHS of 2plate_kite model.
+# Verifies allocation counts for registered functions and
+# ensures no allocations originate from package source code.
 #
 # Usage: julia --project=test test/test_bench.jl
 
+using Test
 using SymbolicAWEModels
 using SymbolicAWEModels: VortexStepMethod, SystemStructure, KVec3
 using KiteUtils
@@ -14,301 +16,158 @@ using Statistics
 using Printf
 using Profile
 
-function setup_sam()
+function setup_bench_sam()
     pkg_root = dirname(@__DIR__)
     src_data = joinpath(pkg_root, "data", "2plate_kite")
-
     tmpdir = mktempdir()
-    data_path = joinpath(tmpdir, "2plate_kite")
-    cp(src_data, data_path; force=true)
-
-    set_data_path(data_path)
+    dp = joinpath(tmpdir, "2plate_kite")
+    cp(src_data, dp; force=true)
+    set_data_path(dp)
     set = Settings("system.yaml")
-
     vsm_set = VortexStepMethod.VSMSettings(
-        joinpath(data_path, "vsm_settings.yaml");
+        joinpath(dp, "vsm_settings.yaml");
         data_prefix=false)
-
-    struc_yaml = joinpath(data_path,
-                          "quat_struc_geometry.yaml")
+    struc_yaml = joinpath(dp, "quat_struc_geometry.yaml")
     sys = load_sys_struct_from_yaml(struc_yaml;
         system_name="bench", set, vsm_set)
     sys.winches[:main_winch].brake = true
-
     sam = SymbolicAWEModel(set, sys)
-    init!(sam; remake=false, prn=true)
+    init!(sam; remake=false, prn=false)
     return sam
 end
 
-# --- 1. RHS Benchmark ---
-function bench_rhs(sam)
-    f = sam.integrator.f
-    u = copy(sam.integrator.u)
-    p = sam.integrator.p
-    t = sam.integrator.t
-    du = similar(u)
-    f(du, u, p, t)
-    return @benchmark $f($du, $u, $p, $t) samples=10
+function get_pkg_src_files()
+    pkg_root = dirname(@__DIR__)
+    src_files = Set{String}()
+    for (_, _, files) in walkdir(joinpath(pkg_root, "src"))
+        for f in files
+            endswith(f, ".jl") && push!(src_files, f)
+        end
+    end
+    return src_files
 end
 
-# --- 2. Direct registered function calls ---
-function bench_registered_functions(sys::SystemStructure)
+@testset verbose = true "Benchmarks" begin
+    sam = setup_bench_sam()
+    sys = sam.sys_struct
     idx = Int64(1)
-    # Warmup
-    SymbolicAWEModels.get_pos_w(sys, idx)
-    SymbolicAWEModels.get_Q_b_to_w(sys, idx)
-    SymbolicAWEModels.get_l0(sys, idx)
-    SymbolicAWEModels.get_extra_mass(sys, idx)
-
-    header = @sprintf("  %-20s %6s\n", "Function", "Allocs")
-    sep = "  " * "-"^30
-
-    # Raw field access
-    println("\n  Raw field access (no @register):")
-    println(sep)
-    print(header)
-    println(sep)
-
-    a = @allocations sys.segments[idx].l0
-    @printf("  %-20s %6d\n", "seg.l0", a)
-
-    a = @allocations sys.points[idx].extra_mass
-    @printf("  %-20s %6d\n", "pt.extra_mass", a)
-
-    a = @allocations sys.points[idx].pos_w
-    @printf("  %-20s %6d\n", "pt.pos_w", a)
-
-    a = @allocations sys.wings[idx].Q_b_to_w
-    @printf("  %-20s %6d\n", "wing.Q_b_to_w", a)
-
-    a = @allocations sys.wings[idx].com_w
-    @printf("  %-20s %6d\n", "wing.com_w", a)
-
-    a = @allocations sys.wings[idx].R_b_to_p
-    @printf("  %-20s %6d\n", "wing.R_b_to_p", a)
-
-    a = @allocations sys.wings[idx].inertia_principal
-    @printf("  %-20s %6d\n", "wing.inertia_princ", a)
-
-    # @register_symbolic calls
-    println("\n  @register_symbolic calls:")
-    println(sep)
-    print(header)
-    println(sep)
-
-    a = @allocations SymbolicAWEModels.get_l0(sys, idx)
-    @printf("  %-20s %6d\n", "get_l0", a)
-
-    a = @allocations SymbolicAWEModels.get_extra_mass(sys, idx)
-    @printf("  %-20s %6d\n", "get_extra_mass", a)
-
-    a = @allocations SymbolicAWEModels.get_pos_w(sys, idx)
-    @printf("  %-20s %6d\n", "get_pos_w", a)
-
-    a = @allocations SymbolicAWEModels.get_Q_b_to_w(sys, idx)
-    @printf("  %-20s %6d\n", "get_Q_b_to_w", a)
-
-    a = @allocations SymbolicAWEModels.get_com_w(sys, idx)
-    @printf("  %-20s %6d\n", "get_com_w", a)
-
-    a = @allocations SymbolicAWEModels.get_R_b_to_p(sys, idx)
-    @printf("  %-20s %6d\n", "get_R_b_to_p", a)
-
-    a = @allocations SymbolicAWEModels.get_inertia_principal(
-        sys, idx)
-    @printf("  %-20s %6d\n", "get_inertia_princ", a)
-
-    # VSM accessors
-    println("\n  VSM accessors:")
-    println(sep)
-    print(header)
-    println(sep)
-
-    a = @allocations SymbolicAWEModels.get_vsm_y(sys, idx, 1)
-    @printf("  %-20s %6d\n", "get_vsm_y", a)
 
-    a = @allocations SymbolicAWEModels.get_vsm_x(sys, idx, 1)
-    @printf("  %-20s %6d\n", "get_vsm_x", a)
-
-    a = @allocations SymbolicAWEModels.get_vsm_jac(
-        sys, idx, 1, 1)
-    @printf("  %-20s %6d\n", "get_vsm_jac", a)
-
-    a = @allocations SymbolicAWEModels.get_aero_force_override(
-        sys, idx, 1)
-    @printf("  %-20s %6d\n", "get_aero_f_override", a)
-
-    a = @allocations SymbolicAWEModels.get_aero_moment_override(
-        sys, idx, 1)
-    @printf("  %-20s %6d\n", "get_aero_m_override", a)
-end
-
-# --- 3. Allocation profile ---
-function profile_rhs_allocs(sam)
-    f = sam.integrator.f
-    u = copy(sam.integrator.u)
-    p = sam.integrator.p
-    t = sam.integrator.t
-    du = similar(u)
-    f(du, u, p, t)
-
-    Profile.Allocs.clear()
-    Profile.Allocs.@profile sample_rate=1.0 f(du, u, p, t)
-    results = Profile.Allocs.fetch()
-
-    type_counts = Dict{Any,Int}()
-    type_bytes = Dict{Any,Int}()
-
-    # Helper: is this a meaningful Julia source frame?
-    function is_julia_src(frame)
-        frame.line <= 0 && return false
-        file = string(frame.file)
-        # Skip C runtime, GC, profiler, sysimage
-        contains(file, "gc-") && return false
-        contains(file, "Profile") && return false
-        contains(file, "datatype.c") && return false
-        endswith(file, ".c") && return false
-        endswith(file, ".h") && return false
-        file == ":-1" && return false
-        return true
+    @testset "ODE RHS" begin
+        f = sam.integrator.f
+        u = copy(sam.integrator.u)
+        p = sam.integrator.p
+        t = sam.integrator.t
+        du = similar(u)
+        f(du, u, p, t)
+
+        rhs = @benchmark $f($du, $u, $p, $t) samples = 10
+        med_ms = median(rhs.times) / 1e6
+        @printf("  median: %8.3f ms | allocs: %d (%d B)\n",
+                med_ms, rhs.allocs, rhs.memory)
+        @test rhs.allocs <= 200
     end
 
-    # Aggregate by source location (first Julia src frame)
-    loc_counts = Dict{String,Int}()
-    loc_bytes = Dict{String,Int}()
-
-    for a in results.allocs
-        T = a.type
-        type_counts[T] = get(type_counts, T, 0) + 1
-        type_bytes[T] = get(type_bytes, T, 0) + a.size
-
-        loc = "unknown"
-        for frame in a.stacktrace
-            is_julia_src(frame) || continue
-            file = string(frame.file)
-            # Skip base Julia math (sys.so compiled)
-            endswith(file, ".so") && continue
-            basename(file) in (
-                "int.jl", "float.jl", "promotion.jl",
-                "number.jl", "boot.jl") && continue
-            loc = "$(basename(file)):$(frame.line)"
-            break
-        end
-        loc_counts[loc] = get(loc_counts, loc, 0) + 1
-        loc_bytes[loc] = get(loc_bytes, loc, 0) + a.size
+    # --- warmup all accessors once ---
+    SymbolicAWEModels.get_l0(sys, idx)
+    SymbolicAWEModels.get_extra_mass(sys, idx)
+    SymbolicAWEModels.get_pos_w(sys, idx)
+    SymbolicAWEModels.get_Q_b_to_w(sys, idx)
+    SymbolicAWEModels.get_com_w(sys, idx)
+    SymbolicAWEModels.get_R_b_to_p(sys, idx)
+    SymbolicAWEModels.get_inertia_principal(sys, idx)
+    SymbolicAWEModels.get_vsm_y(sys, idx, 1)
+    SymbolicAWEModels.get_vsm_x(sys, idx, 1)
+    SymbolicAWEModels.get_vsm_jac(sys, idx, 1, 1)
+    SymbolicAWEModels.get_aero_force_override(sys, idx, 1)
+    SymbolicAWEModels.get_aero_moment_override(sys, idx, 1)
+
+    @testset "@register_symbolic" begin
+        a = @allocations SymbolicAWEModels.get_l0(sys, idx)
+        @test a == 0
+        a = @allocations SymbolicAWEModels.get_extra_mass(
+            sys, idx)
+        @test a == 0
+        a = @allocations SymbolicAWEModels.get_pos_w(
+            sys, idx)
+        @test a <= 1
+        a = @allocations SymbolicAWEModels.get_Q_b_to_w(
+            sys, idx)
+        @test a == 0
+        a = @allocations SymbolicAWEModels.get_com_w(
+            sys, idx)
+        @test a == 0
+        a = @allocations SymbolicAWEModels.get_R_b_to_p(
+            sys, idx)
+        @test a == 0
+        a = @allocations SymbolicAWEModels.get_inertia_principal(
+            sys, idx)
+        @test a == 0
     end
 
-    # Aggregate by call stack (top 5 Julia src frames,
-    # skipping base math)
-    stack_counts = Dict{String,Int}()
-    stack_bytes = Dict{String,Int}()
-    for a in results.allocs
-        frames = String[]
-        for frame in a.stacktrace
-            is_julia_src(frame) || continue
-            fname = basename(string(frame.file))
-            func = string(frame.func)
-            # Truncate long generated function names
-            if length(func) > 40
-                func = func[1:37] * "..."
-            end
-            push!(frames, "$func ($fname:$(frame.line))")
-            length(frames) >= 5 && break
-        end
-        key = isempty(frames) ? "unknown" :
-            join(frames, "\n       <- ")
-        stack_counts[key] = get(stack_counts, key, 0) + 1
-        stack_bytes[key] = get(stack_bytes, key, 0) + a.size
+    @testset "VSM accessors" begin
+        a = @allocations SymbolicAWEModels.get_vsm_y(
+            sys, idx, 1)
+        @test a == 0
+        a = @allocations SymbolicAWEModels.get_vsm_x(
+            sys, idx, 1)
+        @test a == 0
+        a = @allocations SymbolicAWEModels.get_vsm_jac(
+            sys, idx, 1, 1)
+        @test a <= 2
+        a = @allocations SymbolicAWEModels.get_aero_force_override(
+            sys, idx, 1)
+        @test a <= 4
+        a = @allocations SymbolicAWEModels.get_aero_moment_override(
+            sys, idx, 1)
+        @test a == 0
     end
 
-    return (; type_counts, type_bytes,
-              loc_counts, loc_bytes,
-              stack_counts, stack_bytes)
-end
-
-# --- Reporting ---
-function print_rhs_bench(rhs)
-    println("\n", "="^60)
-    println("  ODE RHS Benchmark")
-    println("="^60)
-    @printf("  median: %8.3f ms\n",
-            median(rhs.times) / 1e6)
-    @printf("  allocs: %8d  (%d bytes)\n",
-            rhs.allocs, rhs.memory)
-    println("="^60)
-end
-
-function print_alloc_profile(prof)
-    (; type_counts, type_bytes,
-       loc_counts, loc_bytes,
-       stack_counts, stack_bytes) = prof
-    sorted = sort(collect(type_counts); by=last, rev=true)
-    total = sum(values(type_counts))
-
-    println("\n", "="^60)
-    println("  Allocation Profile by Type ($total total)")
-    println("="^60)
-    @printf("  %-35s %7s %10s\n", "Type", "Count", "Bytes")
-    println("  ", "-"^56)
-    for (T, count) in sorted[1:min(15, end)]
-        name = string(T)
-        if length(name) > 35
-            name = name[1:32] * "..."
+    @testset "No package allocations in RHS" begin
+        f = sam.integrator.f
+        u = copy(sam.integrator.u)
+        p = sam.integrator.p
+        t = sam.integrator.t
+        du = similar(u)
+        f(du, u, p, t)
+
+        Profile.Allocs.clear()
+        Profile.Allocs.@profile sample_rate = 1.0 begin
+            f(du, u, p, t)
         end
-        bytes = get(type_bytes, T, 0)
-        @printf("  %-35s %7d %10d\n", name, count, bytes)
-    end
-    println("="^60)
-
-    # --- By source location ---
-    loc_sorted = sort(collect(loc_counts);
-                      by=last, rev=true)
-    println("\n", "="^70)
-    println("  Top Allocating Source Locations")
-    println("="^70)
-    @printf("  %-50s %7s %10s\n",
-            "Location", "Count", "Bytes")
-    println("  ", "-"^67)
-    for (loc, count) in loc_sorted[1:min(20, end)]
-        name = loc
-        if length(name) > 50
-            name = "..." * name[end-46:end]
+        results = Profile.Allocs.fetch()
+        src_files = get_pkg_src_files()
+
+        skip_bases = ("int.jl", "float.jl", "promotion.jl",
+                      "number.jl", "boot.jl")
+
+        pkg_locs = String[]
+        for a in results.allocs
+            for frame in a.stacktrace
+                frame.line <= 0 && continue
+                file = string(frame.file)
+                any(p -> contains(file, p),
+                    ("gc-", "Profile", "datatype.c")) &&
+                    continue
+                (endswith(file, ".c") ||
+                    endswith(file, ".h")) && continue
+                file == ":-1" && continue
+                fname = basename(file)
+                endswith(fname, ".so") && continue
+                fname in skip_bases && continue
+                if fname in src_files
+                    loc = "$fname:$(frame.line)"
+                    loc ∉ pkg_locs && push!(pkg_locs, loc)
+                end
+                break  # only check first meaningful frame
+            end
         end
-        bytes = get(loc_bytes, loc, 0)
-        @printf("  %-50s %7d %10d\n", name, count, bytes)
-    end
-    println("="^70)
 
-    # --- By call stack ---
-    stack_sorted = sort(collect(stack_counts);
-                        by=last, rev=true)
-    println("\n", "="^70)
-    println("  Top Allocating Call Stacks (top 3 frames)")
-    println("="^70)
-    for (i, (stack, count)) in enumerate(
-            stack_sorted[1:min(15, end)])
-        bytes = get(stack_bytes, stack, 0)
-        @printf("  #%-2d  %5d allocs (%d bytes)\n",
-                i, count, bytes)
-        for frame in split(stack, " <- ")
-            println("       ", frame)
+        if !isempty(pkg_locs)
+            println("  Allocating locations in package src:")
+            for loc in pkg_locs
+                println("    ", loc)
+            end
         end
-        println()
+        @test isempty(pkg_locs)
     end
-    println("="^70)
 end
-
-# --- Run ---
-println("\n>> Setting up 2plate_kite model...")
-sam = setup_sam()
-
-println("\n>> Benchmarking ODE RHS...")
-rhs = bench_rhs(sam)
-print_rhs_bench(rhs)
-
-println("\n>> Benchmarking registered functions directly...")
-bench_registered_functions(sam.sys_struct)
-
-println("\n>> Profiling RHS allocations...")
-prof = profile_rhs_allocs(sam)
-print_alloc_profile(prof)
diff --git a/test/test_match_aero_sections.jl b/test/test_match_aero_sections.jl
new file mode 100644
index 00000000..41780ea7
--- /dev/null
+++ b/test/test_match_aero_sections.jl
@@ -0,0 +1,310 @@
+# Copyright (c) 2025 Bart van de Lint
+# SPDX-License-Identifier: MPL-2.0
+
+# test_match_aero_sections.jl
+# Tests for match_aero_sections_to_structure!:
+# verifies geometry alignment and that use_prior_polar
+# preserves refined panel polars across section count
+# changes (no re-interpolation).
+
+using Test
+using SymbolicAWEModels
+using SymbolicAWEModels: KVec3, VortexStepMethod, WING,
+    REFINE, QUATERNION, SimFloat,
+    match_aero_sections_to_structure!
+using KiteUtils
+using LinearAlgebra
+
+# ── helpers ──────────────────────────────────────────────
+pkg_root = dirname(@__DIR__)
+src_data = joinpath(pkg_root, "data", "2plate_kite")
+tmpdir   = mktempdir()
+data_path = joinpath(tmpdir, "2plate_kite")
+cp(src_data, data_path; force=true)
+set_data_path(data_path)
+
+struc_yaml = joinpath(data_path,
+    "quat_struc_geometry.yaml")
+refine_yaml = joinpath(data_path,
+    "refine_struc_geometry.yaml")
+aero_yaml = joinpath(data_path, "aero_geometry.yaml")
+
+SymbolicAWEModels.update_aero_yaml_from_struc_yaml!(
+    struc_yaml, aero_yaml)
+
+set = Settings("system.yaml")
+set.g_earth = 0.0
+vsm_set_path = joinpath(data_path, "vsm_settings.yaml")
+vsm_set = VortexStepMethod.VSMSettings(
+    vsm_set_path; data_prefix=false)
+
+# ─────────────────────────────────────────────────────────
+@testset "match_aero_sections — REFINE" begin
+
+    @testset "geometry: LE/TE match structural points" begin
+        sys = SymbolicAWEModels.load_sys_struct_from_yaml(
+            refine_yaml;
+            system_name="refine_geom", set, vsm_set)
+        wing = sys.wings[1]
+        points = sys.points
+        vsm_w = wing.vsm_wing
+
+        n_struct = count(
+            p -> p.type == WING &&
+                 p.wing_idx == wing.idx,
+            points) ÷ 2
+
+        # Add extra section to create mismatch
+        extra = deepcopy(vsm_w.unrefined_sections[1])
+        push!(vsm_w.unrefined_sections, extra)
+        vsm_w.n_unrefined_sections =
+            Int16(length(vsm_w.unrefined_sections))
+        @test vsm_w.n_unrefined_sections != n_struct
+
+        vsm_w.use_prior_polar = true
+        wing.wing_segments = nothing
+
+        match_aero_sections_to_structure!(
+            wing, points)
+
+        @test vsm_w.n_unrefined_sections == n_struct
+        @test !isnothing(wing.wing_segments)
+        @test length(wing.wing_segments) == n_struct
+
+        R = wing.R_b_to_c
+        origin = wing.pos_cad
+        for (i, (le_idx, te_idx)) in
+                enumerate(wing.wing_segments)
+            sec = vsm_w.unrefined_sections[i]
+            le_body = R' * (points[le_idx].pos_cad -
+                            origin)
+            te_body = R' * (points[te_idx].pos_cad -
+                            origin)
+            @test isapprox(Vector(sec.LE_point),
+                           le_body; atol=1e-10)
+            @test isapprox(Vector(sec.TE_point),
+                           te_body; atol=1e-10)
+        end
+    end
+
+    @testset "use_prior_polar preserves refined polars" begin
+        sys = SymbolicAWEModels.load_sys_struct_from_yaml(
+            refine_yaml;
+            system_name="refine_polar", set, vsm_set)
+        wing = sys.wings[1]
+        points = sys.points
+        vsm_w = wing.vsm_wing
+
+        # Baseline refined polars from constructor
+        n_refined = length(vsm_w.refined_sections)
+        @test n_refined > 0
+        @test all(
+            s -> !isnothing(s.aero_data),
+            vsm_w.refined_sections)
+        baseline = [
+            deepcopy(sec.aero_data)
+            for sec in vsm_w.refined_sections
+        ]
+
+        n_struct = count(
+            p -> p.type == WING &&
+                 p.wing_idx == wing.idx,
+            points) ÷ 2
+
+        # Set all unrefined polars to a scaled value so
+        # re-interpolation would produce uniform output
+        # that differs from baseline
+        orig_ad = vsm_w.unrefined_sections[1].aero_data
+        @test !isnothing(orig_ad)
+        uniform = Tuple(v .* 2.0 for v in orig_ad)
+        for sec in vsm_w.unrefined_sections
+            sec.aero_data = deepcopy(uniform)
+        end
+
+        # Verify non-trivial: some baseline refined
+        # polars differ from the uniform unrefined data
+        @test any(
+            any(!isapprox(baseline[i][k], uniform[k])
+                for k in eachindex(uniform))
+            for i in eachindex(baseline))
+
+        # Create mismatch: add 2 extra unrefined sections
+        for _ in 1:2
+            push!(vsm_w.unrefined_sections,
+                  deepcopy(vsm_w.unrefined_sections[1]))
+        end
+        vsm_w.n_unrefined_sections =
+            Int16(length(vsm_w.unrefined_sections))
+        @test vsm_w.n_unrefined_sections != n_struct
+
+        vsm_w.use_prior_polar = true
+        wing.wing_segments = nothing
+
+        match_aero_sections_to_structure!(
+            wing, points)
+
+        # Refined count unchanged (n_panels unchanged)
+        @test length(vsm_w.refined_sections) == n_refined
+
+        # Refined polars preserved — NOT re-interpolated
+        for (i, sec) in
+                enumerate(vsm_w.refined_sections)
+            @test !isnothing(sec.aero_data)
+            for k in eachindex(baseline[i])
+                @test sec.aero_data[k] ≈ baseline[i][k]
+            end
+        end
+
+        # Preserved polars differ from rebuilt unrefined
+        # (re-interpolation would have made them equal)
+        unrefined_ad = vsm_w.unrefined_sections[1].aero_data
+        @test any(
+            any(!isapprox(sec.aero_data[k],
+                          unrefined_ad[k])
+                for k in eachindex(sec.aero_data))
+            for sec in vsm_w.refined_sections)
+
+        # Unrefined sections rebuilt to match structure
+        @test vsm_w.n_unrefined_sections == n_struct
+    end
+
+    @testset "errors when use_prior_polar=false" begin
+        sys = SymbolicAWEModels.load_sys_struct_from_yaml(
+            refine_yaml;
+            system_name="refine_err", set, vsm_set)
+        wing = sys.wings[1]
+        points = sys.points
+        vsm_w = wing.vsm_wing
+
+        # Create mismatch
+        push!(vsm_w.unrefined_sections,
+              deepcopy(vsm_w.unrefined_sections[1]))
+        vsm_w.n_unrefined_sections =
+            Int16(length(vsm_w.unrefined_sections))
+
+        vsm_w.use_prior_polar = false
+        wing.wing_segments = nothing
+
+        @test_throws ErrorException match_aero_sections_to_structure!(
+            wing, points)
+    end
+end
+
+# ─────────────────────────────────────────────────────────
+@testset "match_aero_sections — QUATERNION" begin
+
+    @testset "geometry: LE/TE match structural points" begin
+        sys_q = SymbolicAWEModels.load_sys_struct_from_yaml(
+            struc_yaml;
+            system_name="quat_geom", set, vsm_set,
+            wing_type=QUATERNION)
+        wing = sys_q.wings[1]
+        points = sys_q.points
+        vsm_w = wing.vsm_wing
+
+        # wing_segments populated by constructor
+        @test !isnothing(wing.wing_segments)
+
+        n_struct = length(wing.group_idxs)
+        @test length(wing.wing_segments) == n_struct
+        @test vsm_w.n_unrefined_sections == n_struct
+
+        # Verify LE/TE positions match
+        R = wing.R_b_to_c
+        origin = wing.pos_cad
+        for (i, (le_idx, te_idx)) in
+                enumerate(wing.wing_segments)
+            sec = vsm_w.unrefined_sections[i]
+            le_body = R' * (points[le_idx].pos_cad -
+                            origin)
+            te_body = R' * (points[te_idx].pos_cad -
+                            origin)
+            @test isapprox(Vector(sec.LE_point),
+                           le_body; atol=1e-10)
+            @test isapprox(Vector(sec.TE_point),
+                           te_body; atol=1e-10)
+        end
+    end
+
+    @testset "vsm resize after count mismatch" begin
+        sys_q = SymbolicAWEModels.load_sys_struct_from_yaml(
+            struc_yaml;
+            system_name="quat_resize", set, vsm_set,
+            wing_type=QUATERNION)
+        wing = sys_q.wings[1]
+        points = sys_q.points
+        vsm_w = wing.vsm_wing
+
+        n_struct = length(wing.group_idxs)
+
+        # Save baseline refined polars
+        n_refined = length(vsm_w.refined_sections)
+        @test n_refined > 0
+        @test all(
+            s -> !isnothing(s.aero_data),
+            vsm_w.refined_sections)
+        baseline = [
+            deepcopy(sec.aero_data)
+            for sec in vsm_w.refined_sections
+        ]
+
+        # Set all unrefined polars to a scaled value so
+        # re-interpolation would produce uniform output
+        # that differs from baseline
+        orig_ad = vsm_w.unrefined_sections[1].aero_data
+        @test !isnothing(orig_ad)
+        uniform = Tuple(v .* 2.0 for v in orig_ad)
+        for sec in vsm_w.unrefined_sections
+            sec.aero_data = deepcopy(uniform)
+        end
+
+        # Verify non-trivial: some baseline refined
+        # polars differ from the uniform unrefined data
+        @test any(
+            any(!isapprox(baseline[i][k], uniform[k])
+                for k in eachindex(uniform))
+            for i in eachindex(baseline))
+
+        # Add extra section to create mismatch
+        extra = deepcopy(vsm_w.unrefined_sections[1])
+        push!(vsm_w.unrefined_sections, extra)
+        vsm_w.n_unrefined_sections =
+            Int16(length(vsm_w.unrefined_sections))
+        @test vsm_w.n_unrefined_sections != n_struct
+
+        vsm_w.use_prior_polar = true
+        wing.wing_segments = nothing
+
+        match_aero_sections_to_structure!(
+            wing, points;
+            groups=collect(sys_q.groups))
+
+        @test vsm_w.n_unrefined_sections == n_struct
+
+        # Verify linearization vectors resized
+        ny = 3 + n_struct + 3
+        nx = 3 + 3 + n_struct
+        @test length(wing.vsm_y) == ny
+        @test length(wing.vsm_x) == nx
+        @test size(wing.vsm_jac) == (nx, ny)
+
+        # Refined polars preserved — NOT re-interpolated
+        @test length(vsm_w.refined_sections) == n_refined
+        for (i, sec) in
+                enumerate(vsm_w.refined_sections)
+            @test !isnothing(sec.aero_data)
+            for k in eachindex(baseline[i])
+                @test sec.aero_data[k] ≈ baseline[i][k]
+            end
+        end
+
+        # Preserved polars differ from rebuilt unrefined
+        # (re-interpolation would have made them equal)
+        unrefined_ad = vsm_w.unrefined_sections[1].aero_data
+        @test any(
+            any(!isapprox(sec.aero_data[k],
+                          unrefined_ad[k])
+                for k in eachindex(sec.aero_data))
+            for sec in vsm_w.refined_sections)
+    end
+end
diff --git a/test/test_quaternion_auto_groups.jl b/test/test_quaternion_auto_groups.jl
index 99650ac9..3d8a770e 100644
--- a/test/test_quaternion_auto_groups.jl
+++ b/test/test_quaternion_auto_groups.jl
@@ -3,6 +3,8 @@
 
 # Test auto-creation of groups for QUATERNION wings
 using SymbolicAWEModels
+using SymbolicAWEModels: VortexStepMethod, WING,
+    QUATERNION, REFINE
 using Test
 using LinearAlgebra
 
@@ -16,85 +18,78 @@ using LinearAlgebra
     data_path = joinpath(tmpdir, "2plate_kite")
     cp(src_data_path, data_path; force=true)
     set_data_path(data_path)
-    set = Settings("system.yaml")
 
-    # Load with REFINE (should have 0 groups)
-    println("\n=== Testing REFINE wing (no auto-groups) ===")
+    struc_yaml = joinpath(
+        data_path, "quat_struc_geometry.yaml")
     refine_yaml = joinpath(
-        data_path, "refine_struc_geometry.yaml"
-    )
+        data_path, "refine_struc_geometry.yaml")
+    aero_yaml = joinpath(
+        data_path, "aero_geometry.yaml")
+
+    SymbolicAWEModels.update_aero_yaml_from_struc_yaml!(
+        struc_yaml, aero_yaml)
+
+    set = Settings("system.yaml")
     vsm_set_path = joinpath(
-        data_path, "vsm_settings.yaml"
-    )
-    vsm_set = SymbolicAWEModels.VortexStepMethod.VSMSettings(
-        vsm_set_path; data_prefix=false
-    )
+        data_path, "vsm_settings.yaml")
+    vsm_set = VortexStepMethod.VSMSettings(
+        vsm_set_path; data_prefix=false)
+
+    # ── REFINE: should have 0 groups ──────────────
     sys_refine = load_sys_struct_from_yaml(
         refine_yaml;
-        system_name="2plate_refine", set, vsm_set
-    )
+        system_name="2plate_refine", set, vsm_set)
 
     @test length(sys_refine.wings) == 1
-    @test sys_refine.wings[1].wing_type == SymbolicAWEModels.REFINE
+    @test sys_refine.wings[1].wing_type == REFINE
     @test length(sys_refine.groups) == 0
     @test length(sys_refine.wings[1].group_idxs) == 0
-    println("✓ REFINE wing: $(length(sys_refine.groups)) groups")
-
-    # Now test manually creating a QUATERNION wing with WING points
-    println("\n=== Testing QUATERNION wing auto-group creation ===")
 
-    # Get WING points from REFINE system
-    wing_points = [p for p in sys_refine.points if p.type == SymbolicAWEModels.WING]
-    println("Found $(length(wing_points)) WING points")
-    @test length(wing_points) == 6  # 3 LE/TE pairs
-
-    # Create a QUATERNION wing with these points
-    vsm_wing = SymbolicAWEModels.Wing(
-        set, vsm_set; prn=false
-    )
-    vsm_aero = SymbolicAWEModels.BodyAerodynamics([vsm_wing])
-    vsm_solver = SymbolicAWEModels.Solver(vsm_aero;
-                                          solver_type=SymbolicAWEModels.NONLIN,
-                                          atol=2e-8, rtol=2e-8)
+    # ── QUATERNION with YAML-defined groups ───────
+    # quat_struc_geometry.yaml has 3 explicit groups
+    # and 7 WING points (6 LE/TE + kcu).
+    sys_quat = load_sys_struct_from_yaml(
+        struc_yaml;
+        system_name="2plate_quat", set, vsm_set,
+        wing_type=QUATERNION)
 
-    # Create wing with QUATERNION type and empty group_idxs
-    quat_wing = SymbolicAWEModels.VSMWing(
-        SymbolicAWEModels.BaseWing(
-            :main_wing, Int16[], I(3),
-            zeros(3), ones(3);
-            wing_type=SymbolicAWEModels.QUATERNION
-        ),
-        vsm_aero, vsm_wing, vsm_solver,
-        Float64[], Float64[], zeros(0, 0),
-        nothing, nothing,
-        nothing, nothing, nothing, nothing,
-        nothing, nothing, 0.0, 0.0
-    )
+    wing = sys_quat.wings[1]
+    @test wing.wing_type == QUATERNION
+    @test length(sys_quat.groups) == 3
+    @test length(wing.group_idxs) == 3
+    @test !isnothing(wing.wing_segments)
+    @test length(wing.wing_segments) == 3
 
-    # Create SystemStructure (should auto-create groups)
-    # collect() converts NamedCollection to Vector
-    sys_quat = SymbolicAWEModels.SystemStructure(
-        "2plate_quat", set;
-        points=collect(sys_refine.points),
-        segments=collect(sys_refine.segments),
-        pulleys=collect(sys_refine.pulleys),
-        tethers=collect(sys_refine.tethers),
-        winches=collect(sys_refine.winches),
-        wings=[quat_wing],
-        transforms=collect(sys_refine.transforms)
-    )
+    # Geometry was computed from closest VSM panel
+    for group in sys_quat.groups
+        @test !iszero(group.chord)
+        @test !iszero(group.y_airf)
+    end
 
-    # Verify groups were auto-created
-    @test length(sys_quat.groups) == 3  # One group per LE/TE pair
-    @test length(sys_quat.wings[1].group_idxs) == 3
-    @test sys_quat.wings[1].wing_type == SymbolicAWEModels.QUATERNION
+    # ── QUATERNION auto-group creation ────────────
+    # Load without explicit groups: auto-group should
+    # kick in for the 6 LE/TE WING points (kcu is
+    # excluded because it isn't in any group → the
+    # group-based path is unavailable and the fallback
+    # consecutive-pair heuristic uses only even-count
+    # subsets).
+    # Easiest approach: load REFINE YAML as QUATERNION
+    # (REFINE YAML has 6 WING points, no groups, no
+    # kcu WING point).
+    sys_auto = load_sys_struct_from_yaml(
+        refine_yaml;
+        system_name="2plate_auto", set, vsm_set,
+        wing_type=QUATERNION)
 
-    println("✓ QUATERNION wing: $(length(sys_quat.groups)) groups auto-created")
+    wing_auto = sys_auto.wings[1]
+    @test wing_auto.wing_type == QUATERNION
+    @test length(sys_auto.groups) == 3
+    @test length(wing_auto.group_idxs) == 3
+    @test !isnothing(wing_auto.wing_segments)
+    @test length(wing_auto.wing_segments) == 3
 
-    # Check that geometry was computed from closest VSM panel
-    for (i, group) in enumerate(sys_quat.groups)
-        println("  Group $i: le_pos = $(group.le_pos), chord_norm = $(norm(group.chord)), points = $(group.point_idxs)")
-        @test !iszero(group.chord)  # Chord should be computed from panel
-        @test !iszero(group.y_airf)  # y_airf should be computed from panel
+    for group in sys_auto.groups
+        @test !iszero(group.chord)
+        @test !iszero(group.y_airf)
     end
 end