diff --git a/.gitignore b/.gitignore
index a2266a9a..d7c105ab 100644
--- a/.gitignore
+++ b/.gitignore
@@ -46,6 +46,7 @@ LocalPreferences.toml
 # Pipeline-generated polars/dat/cp (regenerate via AirfoilAero/ObjAdapter)
 data/**/polars_neuralfoil/
 data/**/polars_xfoil/
+data/**/generated_neuralfoil/
 data/TUDELFT_V3_KITE/aero_geometry_neuralfoil.yaml
 
 # Local planning docs (not tracked)
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2e73b25a..646c9351 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,12 @@
   and `surfplan_to_aero_yaml`: `:csv` (default, readable) or `:arrow` (binary, ~40×
   faster to load and 2.6× smaller). `read_section_aero` detects the format from the
   file suffix, so a geometry YAML can reference either.
+- `geometry_path` keyword on `obj_to_yaml`, naming the geometry YAML itself
+  instead of always writing `output_dir/geometry.yaml`. Point it outside the
+  table directory and the emitted table references carry the path from the YAML's
+  directory to `output_dir`, which is what the geometry loader resolves them
+  against — so a generated dataset can keep its bulk in a subdirectory while the
+  geometry sits with the hand-written ones.
 - `convert_node_table` and `write_node_rows` rewrite a per-node table in the format
   the destination suffix names. `obj_to_yaml` migrates an existing dataset with
   them when `table_format` differs from what the directory holds, so a dataset
@@ -18,6 +24,32 @@
   instead of its own (`ZEROS`, `0.05`), and never said what `core_radius_fraction`
   measures. It now documents the `Solver` defaults and cites Damiani et al. (2019) for
   the 0.05 cut-off.
+- A remesh under `use_prior_polar` no longer resamples the refined sections'
+  `SectionAero` surface tables down to whatever unrefined sections survive it.
+  `compute_refined_section_interpolation!` reblended contour, `cp` and `cf` from the
+  unrefined sections unconditionally while `aero_data` was preserved, so a wing rebuilt
+  onto fewer structural stations kept full-resolution polars but lost the surface tables
+  pressure integration reads. The reblend is now skipped when the polars are preserved
+  and the refined sections already carry tables.
+
+- `obj_to_yaml` no longer places sections on a wingtip that has closed to a point.
+  `station_indices` spreads its targets over the stations that still have a chord,
+  so the outermost section lands on the last sliceable one. A V3 mesh sliced with
+  the default `wingtip_distance` used to put a zero-chord section at each tip,
+  whose polar was `NaN` and took the whole solve with it; working around it meant
+  guessing a `wingtip_distance` large enough to skip past the tip. That workaround
+  is no longer the default: `wingtip_distance` is now `0.0`, an inset on top of the
+  trim for meshes whose slices just short of the tip are still too thin to analyse.
+
+- Wing sections are normalized to `+y` to `-y` order on load (`normalize_span_order!`),
+  and by `refine!` for wings built through `add_section!`. Panel `y_airf` and `z_airf`
+  follow the order sections are stored in, so a geometry file written the other way
+  round inverted every panel normal, and a wing whose sections were replaced after its
+  panels were built (a structural remesh) inverted them mid-run, one panel at a time as
+  each crossed `spanwise_direction`. `obj_to_yaml` and `surfplan_to_aero_yaml` emit that
+  order too; files of either order keep loading the same.
+- Spanwise distribution plots put `+y` on the left, matching that order and the kite
+  seen from the front.
 
 ### Changed
 - `SolverSettings` now defaults to the same values as `Solver`: `core_radius_fraction`
@@ -40,6 +72,10 @@
 - `read_node_table` parses into a preallocated matrix instead of `reduce(vcat, …)`
   over a generator, which was quadratic in the row count: ~21× faster on a 16 MB
   surface table (2.49 s → 0.12 s), benefiting every existing dataset.
+- `is_show=true` draws into a window named after the plot title instead of into
+  whichever window the backend last used, so a script showing several plots gets
+  one window each and re-running it redraws them in place. `show_plot` takes the
+  window `name` as a keyword.
 
 ## VortexStepMethod v4.0.0 2026-08-03
 
diff --git a/docs/src/airfoil_pipeline.md b/docs/src/airfoil_pipeline.md
index 0ad05192..63fe42d8 100644
--- a/docs/src/airfoil_pipeline.md
+++ b/docs/src/airfoil_pipeline.md
@@ -100,9 +100,12 @@ For each unique airfoil id `j`, `obj_to_yaml` writes into `output_dir`:
 - `geometry.yaml` — `wing_sections` (leading/trailing-edge points) plus `wing_airfoils`
   (each section's `type` and the `.dat`/`.csv` paths above)
 
-A near-vanishing wingtip slice can shrink-wrap to an implausibly thick blob; such a
-degenerate section reuses its nearest valid neighbour's airfoil and polar while keeping its
-own edge positions, and a warning lists the reuse. All floats are rounded to millimetre
+A tip that tapers to a point has no airfoil to slice, so the outermost stations stop at
+the last slice that still has a chord; `wingtip_distance` moves them a further arc length
+inboard when the slices just short of the tip are still too thin to analyse. A
+near-vanishing slice that does get through can shrink-wrap to an implausibly thick blob;
+such a degenerate section reuses its nearest valid neighbour's airfoil and polar while
+keeping its own edge positions, and a warning lists the reuse. All floats are rounded to millimetre
 precision by the single [`write_yaml`](@ref VortexStepMethod.ObjAdapter.write_yaml) writer,
 so generated geometry files stay diff-friendly and consistent.
 
diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index 94c0f3e1..8841f5e7 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -21,6 +21,7 @@ calculate_relative_alpha_and_velocity
 calculate_relative_alpha_and_relative_velocity
 update_effective_angle_of_attack!
 calculate_stall_angle_list
+wing_span_flip
 calculate_circulation_distribution_elliptical_wing
 _compute_reference_velocity_from_distribution
 smooth_circulation!
@@ -44,7 +45,7 @@ update_panel_properties!
 build_interps
 panel_interp_types
 reinit!(wing::AbstractWing)
-reinit!(panel::Panel, section_1::Section, section_2::Section, aero_center, control_point, bound_point_1, bound_point_2, x_airf, y_airf, z_airf, delta, vec, spanwise_direction)
+reinit!(panel::Panel, section_1::Section, section_2::Section, aero_center, control_point, bound_point_1, bound_point_2, x_airf, y_airf, z_airf, delta, vec)
 rotated_te
 calculate_filaments_for_plotting
 ```
@@ -60,7 +61,9 @@ copy_sections_to_refined!
 _apply_refined_section_thetas!
 _panel_thetas_to_section_thetas!
 _interpolate_unrefined_to_refined
-_section_sort_key
+span_order_key
+normalize_span_order!
+can_reuse_prior_refined_surface_tables
 refine_mesh_for_linear_cosine_distribution!
 refine_mesh_by_splitting_provided_sections!
 refine_mesh_with_billowing!
@@ -176,6 +179,8 @@ center_to_com!
 airfoils_from_yaml
 write_geometry_yaml
 resolve_aero_geometry
+table_path_prefix
+prefix_table_paths!
 plot_airfoil_fit
 migrate_node_tables
 ```
@@ -185,6 +190,8 @@ migrate_node_tables
 CurrentModule = Base.get_extension(VortexStepMethod, :VortexStepMethodMakieExt)
 ```
 ```@docs
+display_named
+span_axis
 create_geometry_plot_makie
 plot_line_segment_makie!
 set_axes_equal_makie!
diff --git a/examples/V3_kite.jl b/examples/V3_kite.jl
index 2c09c54b..4a7b831c 100644
--- a/examples/V3_kite.jl
+++ b/examples/V3_kite.jl
@@ -2,16 +2,21 @@ using Pkg
 if Base.active_project() != joinpath(@__DIR__, "Project.toml")
     Pkg.activate(@__DIR__)
 end
-using LinearAlgebra
 using GLMakie
 using MakieControlPlots
 using VortexStepMethod
+using VortexStepMethod.ObjAdapter
+using VortexStepMethod.AirfoilAero: ShrinkWrap, NeuralFoilSolver
 
 PLOT = true
 SAVE_ALL = false
 USE_TEX = false
 DEFORM = false
+NEURALFOIL = true
+# Rolling-ball radius of the shrink wrap; fillets the concave tube-canopy junction.
+MIN_CONCAVE_RADIUS = 0.4
 OUTPUT_DIR = joinpath(dirname(@__DIR__), "output")
+REFERENCE_POINT = [0.422646, 0.0, 9.3667]
 
 project_dir = dirname(@__DIR__)
 literature_paths = [
@@ -24,43 +29,21 @@ literature_paths = [
     joinpath(project_dir, "data", "TUDELFT_V3_KITE", "literature_results",
         "windtunnel_alpha_sweep_beta_00_0_Poland_2025_Rey_5e5.csv"),
 ]
-labels = [
-    "VSM Julia Re=5e5",
-    "CFD Re=5e5",
-    "CFD Re=10e5", #with struts
-    "VSM Python Re=5e5",
-    "WindTunnel Re=5e5" #with struts
-]
 beta_literature_paths = [
     joinpath(project_dir, "data", "TUDELFT_V3_KITE", "literature_results",
         "windtunnel_beta_sweep_alpha_07_4_Poland_2025_Rey_5e5.csv"),
 ]
-beta_labels = [
-    labels[1],
-    "Wind Tunnel Re=5e5 beta sweep alpha=7.4",
-]
-
-# Load YAML settings directly
-settings_path = joinpath(
-    project_dir, "data", "TUDELFT_V3_KITE", "vsm_settings.yaml")
-settings_data = VortexStepMethod.YAML.load_file(settings_path)
-condition_cfg = settings_data["condition"]
-wing_cfg = settings_data["wings"][1]
-solver_cfg = settings_data["solver_settings"]
 
-# Create wing, body_aero, and solver objects using settings
-wing = Wing(
-    joinpath(project_dir, wing_cfg["geometry_file"]);
-    n_panels=wing_cfg["n_panels"],
-    spanwise_distribution=getproperty(
-        VortexStepMethod,
-        Symbol(wing_cfg["spanwise_panel_distribution"])),
-    spanwise_direction=Float64.(wing_cfg["spanwise_direction"]),
-    remove_nan=wing_cfg["remove_nan"],
-)
+settings = VSMSettings(joinpath(project_dir, "data", "TUDELFT_V3_KITE",
+    "vsm_settings.yaml"); data_prefix=false)
+settings.wings[1].geometry_file = joinpath(project_dir,
+    settings.wings[1].geometry_file)
+wing = Wing(settings)
 refine!(wing)
 body_aero = BodyAerodynamics([wing])
 VortexStepMethod.reinit!(body_aero)
+solver = Solver(body_aero, settings)
+solver.reference_point .= REFERENCE_POINT
 
 if DEFORM
     VortexStepMethod.unrefined_deform!(
@@ -72,52 +55,56 @@ if DEFORM
     VortexStepMethod.reinit!(body_aero; init_aero=false)
 end
 
-# Construct Solver using keyword arguments from solver settings
-solver = 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]),
-)
-
-# Extract values for plotting
-wind_speed = condition_cfg["wind_speed"]
-angle_of_attack_deg = condition_cfg["alpha"]
-sideslip_deg = condition_cfg["beta"]
-yaw_rate = condition_cfg["yaw_rate"]
-
-# Set flight conditions from settings
-α0 = deg2rad(angle_of_attack_deg)
-β0 = deg2rad(sideslip_deg)
-set_va!(body_aero,
-    wind_speed .* [cos(α0) * cos(β0), sin(β0), sin(α0) * cos(β0)])
+# Second sweep on generated polars: slice V3_25.obj, shrink-wrap every section into a
+# closed airfoil and sweep it with NeuralFoil, so the same kite flies on polars derived
+# from its own CAD surface instead of the checked-in CFD tables.
+if NEURALFOIL
+    obj_file = joinpath(project_dir, "data", "TUDELFT_V3_KITE", "V3_25.obj")
+    generated_dir = joinpath(project_dir, "data", "TUDELFT_V3_KITE",
+        "generated_neuralfoil")
+    nf_yaml = obj_to_yaml(obj_file, generated_dir;
+        n_sections=settings.wings[1].n_panels, Re=1e6, force=false,
+        aero_solver=NeuralFoilSolver(model_size="large", n_crit=4.0,
+            xtr_upper=0.05, xtr_lower=0.05),
+        wrap_method=ShrinkWrap(clearance=0.0,
+            min_concave_radius=MIN_CONCAVE_RADIUS),
+    )
+    settings_nf = deepcopy(settings)
+    settings_nf.wings[1].geometry_file = nf_yaml
+    wing_nf = Wing(settings_nf)
+    refine!(wing_nf)
+    body_nf = BodyAerodynamics([wing_nf])
+    VortexStepMethod.reinit!(body_nf)
+    solver_nf = Solver(body_nf, settings_nf)
+    solver_nf.reference_point .= REFERENCE_POINT
+
+    # Reading the generated directory instead of the OBJ shows the airfoils the polar
+    # pipeline actually analysed. Hover a slice to inspect its 2D fit.
+    PLOT && plot_slices_3d(generated_dir; obj_path=obj_file)
+end
 
-# Solve
+solvers = NEURALFOIL ? [solver, solver_nf] : [solver]
+bodies = NEURALFOIL ? [body_aero, body_nf] : [body_aero]
+solver_labels = NEURALFOIL ? ["VSM Julia CFD", "VSM Julia NeuralFoil"] :
+                ["VSM Julia CFD"]
+labels = [solver_labels;
+    ["CFD Re=5e5",
+     "CFD Re=10e5", #with struts
+     "VSM Python Re=5e5",
+     "WindTunnel Re=5e5"]] #with struts
+beta_labels = [solver_labels; ["Wind Tunnel Re=5e5 beta sweep alpha=7.4"]]
+
+wind_speed = settings.condition.wind_speed
+angle_of_attack_deg = settings.condition.alpha
+sideslip_deg = settings.condition.beta
+yaw_rate = settings.condition.yaw_rate
+
+set_va!(body_aero, settings)
 results = VortexStepMethod.solve(solver, body_aero; log=true)
 
-# Plotting polars with moment coefficients
 PLOT && plot_polars(
-    [solver],
-    [body_aero],
+    solvers,
+    bodies,
     labels,
     literature_path_list=literature_paths,
     angle_range=range(-5, 25, length=31),
@@ -130,7 +117,8 @@ PLOT && plot_polars(
     is_save=false || SAVE_ALL,
     is_show=true,
     use_tex=USE_TEX,
-    show_moments=true
+    show_moments=false,
+    cl_over_cd=true
 )
 
 # Plotting geometry
@@ -160,56 +148,10 @@ PLOT && plot_distribution(
     use_tex=USE_TEX
 )
 
-# --- Dual solver comparison: NONLIN vs LOOP ---
-solver_cfg["solver_type"] = "LOOP"
-solver_loop = Solver(body_aero;
-    solver_type=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]),
-)
-
-PLOT && plot_polars(
-    [solver_loop],
-    [body_aero],
-    labels;
-    literature_path_list=literature_paths,
-    angle_range=range(-5, 20, step=1),
-    angle_type="angle_of_attack",
-    angle_of_attack=angle_of_attack_deg,
-    side_slip=sideslip_deg,
-    v_a=wind_speed,
-    title="LOOP solver",
-    show_moments=true,
-    save_path=OUTPUT_DIR,
-    is_save=false || SAVE_ALL,
-    is_show=true,
-    use_tex=USE_TEX
-)
-
 # --- Beta sweep ---
 PLOT && plot_polars(
-    [solver_loop],
-    [body_aero],
+    solvers,
+    bodies,
     beta_labels;
     literature_path_list=beta_literature_paths,
     angle_range=range(0, 12, step=1),
@@ -225,4 +167,4 @@ PLOT && plot_polars(
     use_tex=USE_TEX
 )
 
-nothing
\ No newline at end of file
+nothing
diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl
index 635b4fee..fe2705c5 100644
--- a/examples/ram_air_kite.jl
+++ b/examples/ram_air_kite.jl
@@ -44,8 +44,7 @@ alpha_range = -8:2:26
 delta_range = 0:1
 N_SECTIONS = 10
 
-# Inset the tip stations 10 cm along the leading edge so they avoid the near-zero-chord
-# wingtips, which slice to degenerate airfoils that XFoil cannot analyse.
+# At 0.0 XFoil only converges above 16 deg on the near-tip slice this mesh still has.
 WINGTIP_DISTANCE = 0.1
 
 # Shrink-wrap each raw slice into a clean closed airfoil: the distance-field wrap
diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl
index dc6686f5..1356df73 100644
--- a/ext/VortexStepMethodMakieExt.jl
+++ b/ext/VortexStepMethodMakieExt.jl
@@ -12,6 +12,32 @@ export plot_geometry, plot_distribution, plot_polars, save_plot, show_plot,
 const PANEL_MESH_OBSERVABLES = Ref{Union{Nothing,Dict}}(nothing)
 # Global storage for airfoil-skin observables, keyed by body objectid.
 const AIRFOIL_SKIN_OBSERVABLES = Ref{Union{Nothing,Dict}}(nothing)
+const SCREENS = Dict{String,Any}()
+
+"""
+    display_named(fig, name) -> fig
+
+Show `fig` in the window registered under `name`, opening a new titled window for
+a name that has none yet: plotting the same title again replaces its predecessor,
+while a new title gets its own window. Backends without titled screens, such as
+CairoMakie, fall back to a plain `display`.
+"""
+function display_named(fig::Makie.Figure, name::AbstractString)
+    screen = get(SCREENS, name, nothing)
+    if !isnothing(screen) && isopen(screen)
+        display(screen, fig)
+        return fig
+    end
+    screen = try
+        Makie.current_backend().Screen(; title=String(name))
+    catch
+        display(fig)
+        return fig
+    end
+    SCREENS[name] = screen
+    display(screen, fig)
+    return fig
+end
 
 """
     PLATE_FACES
@@ -465,18 +491,19 @@ function VortexStepMethod.save_plot(fig::Makie.Figure, save_path, title; data_ty
 end
 
 """
-    show_plot(fig; dpi=130)
+    show_plot(fig; name="", dpi=130)
 
-Display a Makie figure.
+Display a Makie figure in an interactive session, in the window named `name`.
 
 # Arguments
 - `fig`: Makie Figure object
 
 # Keyword arguments
+- `name`: Window to draw in; reusing a name reuses its window (default: "")
 - `dpi`: Dots per inch for the figure (default: 130) - currently unused in Makie
 """
-function VortexStepMethod.show_plot(fig::Makie.Figure; dpi=130)
-    isinteractive() && display(fig)
+function VortexStepMethod.show_plot(fig::Makie.Figure; name="", dpi=130)
+    isinteractive() && display_named(fig, name)
 end
 
 """
@@ -671,13 +698,20 @@ function VortexStepMethod.plot_geometry(body_aero::BodyAerodynamics, title;
 
     fig = create_geometry_plot_makie(body_aero, title, view_elevation, view_azimuth)
 
-    if is_show && isinteractive()
-        display(fig)
-    end
+    is_show && show_plot(fig; name=title)
 
     return fig
 end
 
+"""
+    span_axis(position, title, ylabel) -> Axis
+
+Axis for a spanwise distribution, `+y` on the left, matching the kite seen from the
+front and the `+y` to `-y` order its sections are stored in.
+"""
+span_axis(position, title, ylabel) =
+    Axis(position; title, xlabel="Spanwise Position y/b", ylabel, xreversed=true)
+
 """
     plot_distribution(y_coordinates_list, results_list, label_list;
                       title="spanwise_distribution", data_type=nothing,
@@ -715,28 +749,22 @@ function VortexStepMethod.plot_distribution(y_coordinates_list, results_list, la
     Label(fig[0, :], title, fontsize=20)
 
     # Row 1: CL, CD, Gamma
-    ax_cl = Axis(fig[1, 1], title="CL Distribution",
-        xlabel="Spanwise Position y/b", ylabel="Lift Coefficient CL")
-    ax_cd = Axis(fig[1, 2], title="CD Distribution",
-        xlabel="Spanwise Position y/b", ylabel="Drag Coefficient CD")
-    ax_gamma = Axis(fig[1, 3], title="Γ Distribution",
-        xlabel="Spanwise Position y/b", ylabel="Circulation Γ")
+    ax_cl = span_axis(fig[1, 1], "CL Distribution", "Lift Coefficient CL")
+    ax_cd = span_axis(fig[1, 2], "CD Distribution", "Drag Coefficient CD")
+    ax_gamma = span_axis(fig[1, 3], "Γ Distribution", "Circulation Γ")
 
     # Row 2: Alpha geometric, alpha at ac, alpha uncorrected
-    ax_alpha_geo = Axis(fig[2, 1], title="α Geometric",
-        xlabel="Spanwise Position y/b", ylabel="Angle of Attack α (deg)")
-    ax_alpha_ac = Axis(fig[2, 2], title="α result (corrected to aerodynamic center)",
-        xlabel="Spanwise Position y/b", ylabel="Angle of Attack α (deg)")
-    ax_alpha_unc = Axis(fig[2, 3], title="α Uncorrected (if VSM, at control point)",
-        xlabel="Spanwise Position y/b", ylabel="Angle of Attack α (deg)")
+    alpha_label = "Angle of Attack α (deg)"
+    ax_alpha_geo = span_axis(fig[2, 1], "α Geometric", alpha_label)
+    ax_alpha_ac = span_axis(fig[2, 2],
+        "α result (corrected to aerodynamic center)", alpha_label)
+    ax_alpha_unc = span_axis(fig[2, 3],
+        "α Uncorrected (if VSM, at control point)", alpha_label)
 
     # Row 3: Force components
-    ax_fx = Axis(fig[3, 1], title="Force in x direction",
-        xlabel="Spanwise Position y/b", ylabel="Fx")
-    ax_fy = Axis(fig[3, 2], title="Force in y direction",
-        xlabel="Spanwise Position y/b", ylabel="Fy")
-    ax_fz = Axis(fig[3, 3], title="Force in z direction",
-        xlabel="Spanwise Position y/b", ylabel="Fz")
+    ax_fx = span_axis(fig[3, 1], "Force in x direction", "Fx")
+    ax_fy = span_axis(fig[3, 2], "Force in y direction", "Fy")
+    ax_fz = span_axis(fig[3, 3], "Force in z direction", "Fz")
 
     # Plot CL
     for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list)
@@ -801,9 +829,7 @@ function VortexStepMethod.plot_distribution(y_coordinates_list, results_list, la
         save_plot(fig, save_path, title, data_type=data_type)
     end
 
-    if is_show && isinteractive()
-        display(fig)
-    end
+    is_show && show_plot(fig; name=title)
 
     return fig
 end
@@ -1023,9 +1049,7 @@ function VortexStepMethod.plot_polars(
         save_plot(fig, save_path, main_title; data_type)
     end
 
-    if is_show && isinteractive()
-        display(fig)
-    end
+    is_show && show_plot(fig; name=title)
 
     return fig
 end
@@ -1082,9 +1106,7 @@ function VortexStepMethod.plot_polar_data(body_aero::BodyAerodynamics;
                 color=:blue, linewidth=0.5, transparency=true)
         end
 
-        if is_show && isinteractive()
-            display(fig)
-        end
+        is_show && show_plot(fig; name="polar data")
         return fig
     else
         throw(ArgumentError(
@@ -1454,9 +1476,7 @@ function VortexStepMethod.plot_combined_analysis(
     colsize!(fig.layout, 1, Relative(0.6))
     colsize!(fig.layout, 2, Relative(0.4))
 
-    if is_show && isinteractive()
-        display(fig)
-    end
+    is_show && show_plot(fig; name=title)
 
     return fig
 end
@@ -1651,7 +1671,7 @@ function ObjAdapter.plot_slices_3d(path::String; n_slices::Int=10, rotation=I,
         m = ObjAdapter.march_edges(vertices, faces; step=span / n_bins)
         le = reduce(hcat, m.le)
         te = reduce(hcat, m.te)
-        idx = ObjAdapter.station_indices(m.arclen, n_slices; wingtip_distance)
+        idx = ObjAdapter.station_indices(m, n_slices; wingtip_distance)
         secs = filter(!isnothing,
                       [ObjAdapter.build_section(vertices, faces, m.le[i], m.te[i],
                                                 m.point[i], m.tangent[i]) for i in idx])
@@ -1726,7 +1746,7 @@ function ObjAdapter.plot_slices_3d(path::String; n_slices::Int=10, rotation=I,
         sel[] = best
     end
 
-    is_show && display(fig)
+    is_show && display_named(fig, "slices: $(basename(rstrip(path, '/')))")
     return fig
 end
 
@@ -1780,7 +1800,7 @@ function ObjAdapter.plot_airfoil_fit(x::Vector, y::Vector; title::String="Airfoi
 
     axislegend(ax; position=:rt)
 
-    is_show && display(fig)
+    is_show && display_named(fig, title)
     return fig, params
 end
 
@@ -1845,7 +1865,7 @@ function ObjAdapter.plot_airfoils(geometry_file::String;
     if is_save && !isnothing(save_path)
         VortexStepMethod.save_plot(fig, save_path, "airfoils"; data_type)
     end
-    is_show && display(fig)
+    is_show && display_named(fig, title)
     return fig
 end
 
diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index 87b514ec..b9a56ca5 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -16,8 +16,11 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru
 - `alpha_dist::MVector{P, Float64}` = zeros(Float64, P)
 - `v_a_dist::MVector{P, Float64}` = zeros(Float64, P)
 - `work_vectors`::NTuple{10, MVec3} = ntuple(_ -> zeros(MVec3), 10)
-- `AIC::Array{Float64, 3}` = zeros(P, P, 3): influence coefficients, component last so
-                        that each `AIC[:, :, k]` slice is a contiguous BLAS matrix
+- `AIC::Array{Float64, 3}` = zeros(P, P, 3): control-point influence coefficients, the
+                        matrix the circulation is solved against; component last so that
+                        each `AIC[:, :, k]` slice is a contiguous BLAS matrix
+- `AIC_aero_center::Array{Float64, 3}` = zeros(P, P, 3): aerodynamic-centre (LLT)
+                        influence coefficients, used only for the corrected angle of attack
 - `projected_area::Float64` = 1.0: The area projected onto the xy-plane of the kite body reference frame [m²]
 - `c_ref::Float64` = 1.0: Reference chord length (max panel chord) [m]
 - `y::MVector{P, Float64}` = MVector{P,Float64}(zeros(P))
@@ -37,6 +40,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru
     v_a_dist::MVector{P, T} = zeros(MVector{P, T})
     work_vectors::NTuple{10, MVector{3, T}} = ntuple(_ -> zeros(MVector{3, T}), 10)
     AIC::Array{T, 3} = zeros(T, P, P, 3)
+    AIC_aero_center::Array{T, 3} = zeros(T, P, P, 3)
     projected_area::T = one(T)
     c_ref::T = one(T)
     y::MVector{P, T} = zeros(MVector{P, T})
@@ -118,6 +122,17 @@ function BodyAerodynamics(
     return body_aero
 end
 
+"""
+    wing_span_flip(wing) -> Int8
+
+`-1` when `wing`'s sections run against its `spanwise_direction`, `+1` otherwise: the
+`flip` every panel of the wing is reinitialized with ([`reinit!`](@ref)). One answer per
+wing, so neighbouring panels cannot disagree and invert a single normal by 180°.
+"""
+wing_span_flip(wing) =
+    dot(first(wing.refined_sections).LE_point - last(wing.refined_sections).LE_point,
+        wing.spanwise_direction) < 0 ? Int8(-1) : Int8(1)
+
 function Base.getproperty(obj::BodyAerodynamics, sym::Symbol)
     if sym === :va
         if getfield(obj, :has_distributed_va)
@@ -248,12 +263,13 @@ function reinit!(body_aero::BodyAerodynamics{P, W, T};
 ) where {P, W, T}
     idx = 1
     vec = zeros(MVector{3, T})
-    for wing in body_aero.wings
+    for (wing_idx, wing) in enumerate(body_aero.wings)
         reinit!(wing)
         validate_section_aero(wing.refined_sections)
         panel_props = wing.panel_props
         wing_init_aero = init_aero && !_can_skip_panel_aero_reinit(wing, body_aero.panels, idx)
-        
+        wing_flip = wing_span_flip(wing) == -1
+
         # Create panels
         for i in 1:wing.n_panels
             if length(wing.delta_dist) > 0
@@ -263,7 +279,7 @@ function reinit!(body_aero::BodyAerodynamics{P, W, T};
                 delta = zero(T)
             end
             @views reinit!(
-                body_aero.panels[idx], 
+                body_aero.panels[idx],
                 wing.refined_sections[i],
                 wing.refined_sections[i+1],
                 panel_props.aero_centers[i, :],
@@ -274,10 +290,10 @@ function reinit!(body_aero::BodyAerodynamics{P, W, T};
                 panel_props.y_airf[i, :],
                 panel_props.z_airf[i, :],
                 delta,
-                vec,
-                wing.spanwise_direction;
+                vec;
                 remove_nan=wing.remove_nan,
-                init_aero=wing_init_aero
+                init_aero=wing_init_aero,
+                flip=wing_flip
             )
             body_aero.panels[idx].crease_frac = wing.crease_frac
             idx += 1
@@ -370,7 +386,8 @@ Returns: nothing
 @inline function calculate_AIC_matrices!(body_aero::BodyAerodynamics{P, W, T}, model::Model,
                               core_radius_fraction,
                               va_norm_array::AbstractVector{T},
-                              va_unit_array::AbstractMatrix{T}) where {P, W, T}
+                              va_unit_array::AbstractMatrix{T},
+                              target::AbstractArray{T, 3}=body_aero.AIC) where {P, W, T}
     # Determine evaluation point based on model
     evaluation_point = model == VSM ? :control_point : :aero_center
     evaluation_point_on_bound = model == LLT
@@ -424,7 +441,7 @@ Returns: nothing
                 velocity_induced .-= U_2D
             end
             @inbounds for k in 1:3
-                body_aero.AIC[icp, jring, k] = velocity_induced[k]
+                target[icp, jring, k] = velocity_induced[k]
             end
         end
     end
@@ -482,11 +499,14 @@ function update_effective_angle_of_attack!(alpha_corrected,
     va_norm_array,
     va_unit_array)
 
-    calculate_AIC_matrices!(body_aero, LLT, core_radius_fraction, va_norm_array, va_unit_array)
+    # Its own buffer: `AIC` holds the control-point matrix the circulation was solved
+    # against, so overwriting it here would leave post-solve readers on the LLT one.
+    calculate_AIC_matrices!(body_aero, LLT, core_radius_fraction, va_norm_array,
+                            va_unit_array, body_aero.AIC_aero_center)
 
     induced_velocity = body_aero.cache[1][va_array]
     for k in 1:3
-        mul!(view(induced_velocity, :, k), view(body_aero.AIC, :, :, k), gamma)
+        mul!(view(induced_velocity, :, k), view(body_aero.AIC_aero_center, :, :, k), gamma)
     end
 
     # In-place relative velocity calculation
diff --git a/src/obj_adapter/obj_slice.jl b/src/obj_adapter/obj_slice.jl
index 592f958a..729895f7 100644
--- a/src/obj_adapter/obj_slice.jl
+++ b/src/obj_adapter/obj_slice.jl
@@ -433,6 +433,9 @@ built (`build_section`) at the marched station nearest each equal
 **leading-edge arc-length** target. Each section is
 `(; LE_point, TE_point, span_dir, contour3d, x_airfoil, y_airfoil)`.
 
+Stations closed to a point at the tips are skipped ([`station_indices`](@ref)), and
+`wingtip_distance` insets the outermost sections a further arc length.
+
 The slicer assumes `x` = chordwise, `y` = spanwise, `z` = up. Pass a `3×3` rotation
 matrix to reorient a mesh stored in another convention before slicing.
 """
@@ -442,7 +445,7 @@ function perpendicular_sections(vertices, faces, n_sections; n_bins=60, rotation
     ys = [v[2] for v in vertices]
     m = march_edges(vertices, faces; step=(maximum(ys) - minimum(ys)) / n_bins)
     out = NamedTuple[]
-    for i in station_indices(m.arclen, n_sections; wingtip_distance)
+    for i in station_indices(m, n_sections; wingtip_distance)
         sec = build_section(vertices, faces, m.le[i], m.te[i], m.point[i], m.tangent[i])
         sec === nothing || push!(out, sec)
     end
@@ -450,17 +453,21 @@ function perpendicular_sections(vertices, faces, n_sections; n_bins=60, rotation
 end
 
 """
-    station_indices(arclen, n; wingtip_distance=0.0) -> Vector{Int}
-
-Indices of the marched stations nearest `n` targets spread over the leading-edge arc
-length. The first and last targets sit `wingtip_distance` (arc length) in from the
-tips; with the default `0.0` they land exactly on the tips. Inset the tips a little
-(e.g. `0.1` m) to avoid degenerate near-zero-chord tip sections that some solvers
-(XFoil) cannot analyse.
+    station_indices(march, n; wingtip_distance=0.0, min_chord_frac=0.01) -> Vector{Int}
+
+Indices of the [`march_edges`](@ref) stations nearest `n` targets spread over the
+leading-edge arc length. Stations whose chord has closed to less than
+`min_chord_frac` of the longest one are left out of that range first, so a wing
+tapering to a point puts its outermost sections on the last stations that still
+have an airfoil to slice rather than on the point itself. The remaining first and
+last targets sit a further `wingtip_distance` (arc length) inboard.
 """
-function station_indices(arclen, n; wingtip_distance=0.0)
-    total = arclen[end]
-    d = clamp(wingtip_distance, 0.0, total / 2)
-    n == 1 && return [argmin(abs.(arclen .- total / 2))]
-    return [argmin(abs.(arclen .- t)) for t in range(d, total - d, n)]
+function station_indices(march, n; wingtip_distance=0.0, min_chord_frac=0.01)
+    chords = [norm(te .- le) for (le, te) in zip(march.le, march.te)]
+    usable = findall(≥(min_chord_frac * maximum(chords)), chords)
+    arclen = march.arclen
+    inner, outer = arclen[first(usable)], arclen[last(usable)]
+    d = clamp(wingtip_distance, 0.0, (outer - inner) / 2)
+    n == 1 && return [argmin(abs.(arclen .- (inner + outer) / 2))]
+    return [argmin(abs.(arclen .- t)) for t in range(inner + d, outer - d, n)]
 end
diff --git a/src/obj_adapter/obj_to_yaml.jl b/src/obj_adapter/obj_to_yaml.jl
index 78425319..bc05908f 100644
--- a/src/obj_adapter/obj_to_yaml.jl
+++ b/src/obj_adapter/obj_to_yaml.jl
@@ -45,7 +45,10 @@ Convert a 3D wing `.obj` mesh to the native YAML geometry route.
 Stations are placed at equal leading-edge arc-length intervals and sliced
 perpendicular to the local span (see [`perpendicular_sections`](@ref)), which
 keeps the airfoil undistorted near curved tips; each shape is then shrink-wrapped
-into a clean airfoil and evaluated with `aero_solver`.
+into a clean airfoil and evaluated with `aero_solver`. A tip that tapers to a
+point carries no airfoil, so the outermost stations stop at the last slice that
+still has a chord ([`station_indices`](@ref)); `wingtip_distance` moves them a
+further arc length inboard.
 
 `aero_solver` selects the 2D-airfoil backend: [`NeuralFoilSolver`](@ref) (default,
 fast) or [`XFoilSolver`](@ref) (viscous panel code); pass `aero_solver=XFoilSolver()`
@@ -54,8 +57,13 @@ to use XFoil instead. Each section's polar is written as `POLAR_VECTORS`.
 `crease_frac` is the chordwise hinge location (0–1) about which each `delta_range`
 trailing-edge deflection pivots.
 
-With `force=false` (default) an existing `geometry.yaml` in `output_dir` is reused;
-`force=true` regenerates it (e.g. after changing `delta_range` or the mesh).
+With `force=false` (default) an existing geometry YAML is reused; `force=true`
+regenerates it (e.g. after changing `delta_range` or the mesh).
+
+`geometry_path` names the YAML itself, `output_dir/geometry.yaml` by default. Point it
+elsewhere to keep the YAML out of the table directory — the emitted table references
+then carry the path from the YAML's directory to `output_dir`, which is what the
+geometry loader resolves them against.
 
 `wrap_method` ([`ShrinkWrap`](@ref)) wraps each slice's point cloud into a clean closed
 airfoil, robust both to the noisy interior-structure points (ribs, spars) of a ram-air
@@ -99,6 +107,9 @@ Rewrite a generated dataset's per-node `Cp`/`cf` tables in `table_format` and po
 an existing directory change format without re-running the airfoil solver that
 produced it — the polars are the slow part and they are untouched. The source tables
 are left in place.
+
+`output_dir` is what the YAML's relative table references resolve against, which is
+its own directory — the rule the geometry loader follows.
 """
 function migrate_node_tables(yaml_path::String, output_dir::String,
                              table_format::Symbol; verbose::Bool=true)
@@ -138,6 +149,41 @@ function migrate_node_tables(yaml_path::String, output_dir::String,
     return yaml_path
 end
 
+"""
+    table_path_prefix(geometry_path, output_dir) -> String
+
+Path from the geometry YAML's directory to the table directory, empty when they are
+the same. Table references resolve against the YAML's own directory, so a YAML
+written outside `output_dir` has to carry this hop.
+"""
+function table_path_prefix(geometry_path::String, output_dir::String)
+    yaml_dir = dirname(abspath(geometry_path))
+    tables = abspath(output_dir)
+    yaml_dir == tables && return ""
+    return relpath(tables, yaml_dir)
+end
+
+"""
+    prefix_table_paths!(airfoil_rows, prefix) -> airfoil_rows
+
+Prepend `prefix` to every relative table reference in the `info_dict` of each row, so
+a geometry YAML written outside the table directory still resolves them. A no-op on
+an empty prefix.
+"""
+function prefix_table_paths!(airfoil_rows, prefix::String)
+    isempty(prefix) && return airfoil_rows
+    for row in airfoil_rows
+        info = row[end]
+        info isa AbstractDict || continue
+        for (key, value) in info
+            (endswith(String(key), "_file") || key == "csv_file_path") || continue
+            value isa AbstractString && !isabspath(value) &&
+                (info[key] = joinpath(prefix, value))
+        end
+    end
+    return airfoil_rows
+end
+
 function obj_to_yaml(obj_path::String, output_dir::String;
                      n_sections::Int, Re::Real,
                      alpha_range=-180:1:180, delta_range=nothing,
@@ -145,17 +191,20 @@ function obj_to_yaml(obj_path::String, output_dir::String;
                      wrap_method::ShrinkWrap=ShrinkWrap(),
                      reuse_valid_airfoils::Bool=true, max_thickness_ratio::Real=2.0,
                      spanwise_direction=[0.0, 1.0, 0.0], rotation=I,
-                     wingtip_distance=0.05, crease_frac=0.75, force::Bool=false,
-                     verbose::Bool=true, table_format::Symbol=:csv)
+                     wingtip_distance=0.0, crease_frac=0.75, force::Bool=false,
+                     verbose::Bool=true, table_format::Symbol=:csv,
+                     geometry_path::String=joinpath(output_dir, "geometry.yaml"))
     (!endswith(obj_path, ".obj")) && (obj_path *= ".obj")
     isfile(obj_path) || error("OBJ file not found: $obj_path")
     !isapprox(spanwise_direction, [0.0, 1.0, 0.0]) &&
         throw(ArgumentError("Spanwise direction has to be [0.0, 1.0, 0.0]"))
 
-    yaml_path = joinpath(output_dir, "geometry.yaml")
+    yaml_path = geometry_path
+    mkpath(dirname(abspath(yaml_path)))
     if !force && isfile(yaml_path)
         verbose && @info "Reusing existing geometry (force=true to regenerate)" yaml_path
-        migrate_node_tables(yaml_path, output_dir, table_format; verbose)
+        migrate_node_tables(yaml_path, dirname(abspath(yaml_path)), table_format;
+                            verbose)
         return yaml_path
     end
 
@@ -188,6 +237,7 @@ function obj_to_yaml(obj_path::String, output_dir::String;
         delta_range, aero_solver, reuse_valid_airfoils, crease_frac, verbose,
         table_format)
     isempty(ok) && error("No section produced a valid polar in $obj_path")
+    prefix_table_paths!(airfoil_rows, table_path_prefix(yaml_path, output_dir))
 
     # Each section uses its nearest airfoil that actually produced a polar — covering
     # both too-thick (degenerate) fits and sections the solver could not converge.
@@ -202,7 +252,7 @@ function obj_to_yaml(obj_path::String, output_dir::String;
                                 s.TE_point[1], s.TE_point[2], s.TE_point[3]])
     end
 
-    sort!(section_rows; by = row -> row[3])          # clean spanwise order (by LE_y)
+    sort!(section_rows; by = row -> row[3], rev = true)   # +y to -y, by LE_y
     sort!(airfoil_rows; by = row -> row[1])          # airfoils by id
     write_geometry_yaml(yaml_path, section_rows, airfoil_rows)
     verbose && @info "Wrote geometry to $yaml_path ($(length(section_rows)) sections)"
diff --git a/src/panel.jl b/src/panel.jl
index 7dd980f7..782ef330 100644
--- a/src/panel.jl
+++ b/src/panel.jl
@@ -227,14 +227,14 @@ end
 
 """
     reinit!(panel, section_1, section_2, aero_center, control_point, bound_point_1,
-            bound_point_2, x_airf, y_airf, z_airf, delta, vec, spanwise_direction; kwargs...)
+            bound_point_2, x_airf, y_airf, z_airf, delta, vec; kwargs...)
 
 Reinitialize a panel's geometry, horseshoe filaments and aerodynamic interpolations.
 
-The panel is oriented so its `y_airf` (and the bound vortex `bound_2 -> bound_1`) points
-along `+spanwise_direction`, with `z_airf` pointing to the airfoil upper surface. This
-makes the aero independent of section ordering: a reversed order would otherwise flip the
-normal and make the panel look up its polar at a negated angle of attack.
+`flip` reverses the section order so `y_airf` points along `+spanwise_direction` and
+`z_airf` to the airfoil upper surface, making the aero independent of section ordering.
+The caller owns it and must not derive it from the live geometry, which would invert
+normals mid-run; [`reinit!(::BodyAerodynamics)`](@ref) decides it once per wing.
 """
 function reinit!(
     panel::Panel,
@@ -248,12 +248,11 @@ function reinit!(
     y_airf,
     z_airf,
     delta,
-    vec,
-    spanwise_direction;
+    vec;
     init_aero = true,
-    remove_nan = true
+    remove_nan = true,
+    flip::Bool = false
 )
-    flip = dot(y_airf, spanwise_direction) < 0
     if flip
         section_1, section_2 = section_2, section_1
         bound_point_1, bound_point_2 = bound_point_2, bound_point_1
diff --git a/src/surfplan_adapter/SurfplanAdapter.jl b/src/surfplan_adapter/SurfplanAdapter.jl
index e1cbe40c..a723daab 100644
--- a/src/surfplan_adapter/SurfplanAdapter.jl
+++ b/src/surfplan_adapter/SurfplanAdapter.jl
@@ -83,7 +83,7 @@ function surfplan_to_aero_yaml(adapter_dir::AbstractString, output_dir::Abstract
         push!(section_rows, Any[fid, s.LE_point[1], s.LE_point[2], s.LE_point[3],
                                 s.TE_point[1], s.TE_point[2], s.TE_point[3]])
     end
-    sort!(section_rows; by = row -> row[3])
+    sort!(section_rows; by = row -> row[3], rev = true)   # +y to -y, by LE_y
     sort!(airfoil_rows; by = row -> row[1])
     write_geometry_yaml(yaml_path, section_rows, airfoil_rows)
     verbose && @info "Wrote pressure geometry to $yaml_path " *
diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl
index 24a3556f..be6d5ebd 100644
--- a/src/wing_geometry.jl
+++ b/src/wing_geometry.jl
@@ -50,11 +50,24 @@ function Section(LE_point, TE_point, aero_model, aero_data, section_aero=nothing
 end
 
 """
-    reinit!(section::Section, LE_point::PosVector, TE_point::PosVector, aero_model=nothing, aero_data=nothing)
+    span_order_key(section) -> Float64
 
-Function to update a [Section](@ref) in place.
+Spanwise coordinate that [`normalize_span_order!`](@ref) and `refine!`'s
+`sort_sections` order sections on.
 """
-@inline _section_sort_key(s::Section) = s.LE_point[2]
+@inline span_order_key(s::Section) = s.LE_point[2]
+
+"""
+    normalize_span_order!(sections) -> sections
+
+Reverse `sections` if they do not already run `+y` to `-y`, the order panel normals
+are built from. Use `refine!`'s `sort_sections` for a scrambled list.
+"""
+function normalize_span_order!(sections)
+    length(sections) > 1 && span_order_key(last(sections)) >
+        span_order_key(first(sections)) && reverse!(sections)
+    return sections
+end
 
 function reinit!(section::Section, LE_point, TE_point, aero_model=nothing,
                 aero_data=nothing, section_aero=nothing)
@@ -812,6 +825,18 @@ end
     return all(_has_initialized_section_aero_data, wing.refined_sections)
 end
 
+"""
+    can_reuse_prior_refined_surface_tables(wing) -> Bool
+
+Whether every refined section already carries a [`SectionAero`](@ref), so a
+remesh can keep them instead of reblending from the unrefined sections. False on
+a wing that has none, where the blend is what fills them in the first place.
+"""
+@inline function can_reuse_prior_refined_surface_tables(wing::AbstractWing)
+    isempty(wing.refined_sections) && return false
+    return all(s -> !isnothing(s.section_aero), wing.refined_sections)
+end
+
 """
     copy_sections_to_refined!(wing; reuse_aero_data=false)
 
@@ -905,14 +930,16 @@ function refine!(wing::AbstractWing{T}; recompute_mapping=true, sort_sections=tr
     if sort_sections
         sorted = true
         for i in 1:(length(wing.unrefined_sections) - 1)
-            if _section_sort_key(wing.unrefined_sections[i]) <
-               _section_sort_key(wing.unrefined_sections[i+1])
+            if span_order_key(wing.unrefined_sections[i]) <
+               span_order_key(wing.unrefined_sections[i+1])
                 sorted = false
                 break
             end
         end
         sorted || sort!(wing.unrefined_sections;
-            by=_section_sort_key, rev=true)
+            by=span_order_key, rev=true)
+    elseif recompute_mapping
+        normalize_span_order!(wing.unrefined_sections)
     end
 
     n_sections = wing.n_panels + 1
@@ -924,7 +951,7 @@ function refine!(wing::AbstractWing{T}; recompute_mapping=true, sort_sections=tr
             copy_sections_to_refined!(wing; reuse_aero_data)
             if recompute_mapping
                 compute_refined_panel_mapping!(wing)
-                compute_refined_section_interpolation!(wing)
+                compute_refined_section_interpolation!(wing; reuse_aero_data)
             end
             update_non_deformed_sections!(wing)
             return nothing
@@ -939,7 +966,7 @@ function refine!(wing::AbstractWing{T}; recompute_mapping=true, sort_sections=tr
         copy_sections_to_refined!(wing; reuse_aero_data)
         if recompute_mapping
             compute_refined_panel_mapping!(wing)
-            compute_refined_section_interpolation!(wing)
+            compute_refined_section_interpolation!(wing; reuse_aero_data)
         end
         update_non_deformed_sections!(wing)
         return nothing
@@ -959,7 +986,7 @@ function refine!(wing::AbstractWing{T}; recompute_mapping=true, sort_sections=tr
             reuse_aero_data ? nothing : s2.aero_data)
         if recompute_mapping
             compute_refined_panel_mapping!(wing)
-            compute_refined_section_interpolation!(wing)
+            compute_refined_section_interpolation!(wing; reuse_aero_data)
         end
         update_non_deformed_sections!(wing)
         return nothing
@@ -983,7 +1010,7 @@ function refine!(wing::AbstractWing{T}; recompute_mapping=true, sort_sections=tr
     # Compute panel mapping by finding closest unrefined section for each refined panel
     if recompute_mapping
         compute_refined_panel_mapping!(wing)
-        compute_refined_section_interpolation!(wing)
+        compute_refined_section_interpolation!(wing; reuse_aero_data)
     end
 
     # Update n_unrefined_sections based on actual sections
@@ -1066,7 +1093,8 @@ function compute_refined_panel_mapping!(wing::AbstractWing)
 end
 
 """
-    compute_refined_section_interpolation!(wing::AbstractWing)
+    compute_refined_section_interpolation!(wing::AbstractWing;
+                                           reuse_aero_data=false)
 
 Compute per-refined-section linear-interpolation weights from the unrefined
 sections. For refined section i, the interpolated value is:
@@ -1078,8 +1106,15 @@ Positions are quarter-chord arc-length along the unrefined and refined sections.
 The first refined section is pinned to `left_idx == 1`, `weight == 1` (returns
 `unrefined[1]` exactly) and the last refined section to `left_idx == n_unref - 1`,
 `weight == 0` (returns `unrefined[end]` exactly).
+
+`reuse_aero_data` keeps the surface tables the refined sections already hold
+instead of reblending them from the unrefined ones, the same preservation
+[`refine!`](@ref) applies to `aero_data` under `use_prior_polar`. A remesh that
+replaces the unrefined sections with a coarser set would otherwise resample the
+surface tables down to that set even while the polars stay at full resolution.
 """
-function compute_refined_section_interpolation!(wing::AbstractWing{T}) where {T}
+function compute_refined_section_interpolation!(wing::AbstractWing{T};
+        reuse_aero_data::Bool=false) where {T}
     n_unref = length(wing.unrefined_sections)
     n_sections = wing.n_panels + 1
 
@@ -1147,7 +1182,8 @@ function compute_refined_section_interpolation!(wing::AbstractWing{T}) where {T}
     wing.refined_section_left_idx[n_sections] = Int16(n_unref - 1)
     wing.refined_section_weight[n_sections] = zero(T)
 
-    interpolate_section_aero_to_refined!(wing)
+    keep = reuse_aero_data && can_reuse_prior_refined_surface_tables(wing)
+    keep || interpolate_section_aero_to_refined!(wing)
     return nothing
 end
 
diff --git a/src/yaml_geometry.jl b/src/yaml_geometry.jl
index 94e101df..11da48ac 100644
--- a/src/yaml_geometry.jl
+++ b/src/yaml_geometry.jl
@@ -19,6 +19,8 @@ end
     TE_z::Float64
 end
 
+span_order_key(section::WingSectionData) = section.LE_y
+
 @with_kw struct WingAirfoilData
     airfoil_id::Int64
     type::String
@@ -247,7 +249,8 @@ function Wing(
             TE_z = section_dict["TE_z"]
         ))
     end
-    
+    normalize_span_order!(sections)
+
     # Convert wing airfoils
     wing_airfoils_data = data["wing_airfoils"]
     airfoils = WingAirfoilData[]
diff --git a/test/obj_adapter/test_obj_adapter.jl b/test/obj_adapter/test_obj_adapter.jl
index 322113ed..e2f0ae91 100644
--- a/test/obj_adapter/test_obj_adapter.jl
+++ b/test/obj_adapter/test_obj_adapter.jl
@@ -90,6 +90,29 @@ obj_path = normpath(joinpath(@__DIR__, "..", "..",
                                                                   :csv)
     end
 
+    @testset "geometry_path writes the YAML outside the table directory" begin
+        root = mktempdir()
+        tables = joinpath(root, "tables")
+        yaml_path = joinpath(root, "nf_aero_geometry.yaml")
+        written = obj_to_yaml(obj_path, tables; n_sections=3, Re=5e5,
+            verbose=false, geometry_path=yaml_path)
+        @test written == yaml_path
+        @test isfile(yaml_path)
+        @test !isfile(joinpath(tables, "geometry.yaml"))
+
+        info = Dict(YAML.load_file(yaml_path)["wing_airfoils"]["data"][1][3])
+        for key in ("csv_file_path", "dat_file")
+            @test startswith(info[key], "tables/")
+            @test isfile(joinpath(dirname(yaml_path), info[key]))
+        end
+        # The loader resolves references against the YAML's own directory, so the
+        # prefixed paths have to be what makes this work.
+        @test Wing(yaml_path; n_panels=4) isa Wing
+
+        @test ObjAdapter.table_path_prefix(joinpath(tables, "geometry.yaml"),
+                                           tables) == ""
+    end
+
     @testset "center_to_com! rejects non-triangular faces" begin
         verts = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.0]]
         @test_throws ArgumentError center_to_com!(verts, [[1, 2, 3, 4]]; prn=false)
diff --git a/test/panel/test_panel.jl b/test/panel/test_panel.jl
index 7dc014cb..eaeb6f8a 100644
--- a/test/panel/test_panel.jl
+++ b/test/panel/test_panel.jl
@@ -48,8 +48,7 @@ function create_panel(section1::Section, section2::Section)
         y_airf,
         z_airf,
         0.0,
-        zeros(MVec3),
-        MVec3([0.0, 1.0, 0.0])
+        zeros(MVec3)
     )
     return panel
 end
diff --git a/test/wing_geometry/test_wing_geometry.jl b/test/wing_geometry/test_wing_geometry.jl
index 14f92408..c78a5dba 100644
--- a/test/wing_geometry/test_wing_geometry.jl
+++ b/test/wing_geometry/test_wing_geometry.jl
@@ -433,6 +433,60 @@ end
         @test wing_no_reuse.refined_sections[3].aero_data != no_reuse_baseline
     end
 
+    @testset "Surface tables across a coarsening remesh" begin
+        surface_table(scale) = VortexStepMethod.SectionAero(
+            deg2rad.([-5.0, 0.0, 5.0]), [0.0],
+            reshape(collect(range(0.0, 1.0; length=4)), 4, 1),
+            fill(0.05, 4, 1), fill(scale, 4, 3, 1), fill(0.1 * scale, 4, 3, 1))
+
+        alpha = deg2rad.([-5.0, 0.0, 5.0])
+        polar = (collect(alpha), [1.0, 1.2, 1.4], [0.1, 0.11, 0.12],
+                 [0.01, 0.02, 0.03])
+        span_ys = [1.0, 0.5, 0.0, -0.5, -1.0]
+        scales = [1.0, 2.0, 10.0, 4.0, 5.0]
+
+        function fine_wing(; use_prior_polar, tables=true)
+            wing = Wing(8; spanwise_distribution=LINEAR, use_prior_polar)
+            for (y, scale) in zip(span_ys, scales)
+                add_section!(wing, [0.0, y, 0.0], [1.0, y, 0.0], POLAR_VECTORS,
+                             polar, tables ? surface_table(scale) : nothing)
+            end
+            refine!(wing)
+            return wing
+        end
+
+        function keep_endpoints!(wing)
+            wing.unrefined_sections = [wing.unrefined_sections[1],
+                                       wing.unrefined_sections[end]]
+            wing.n_unrefined_sections = Int16(2)
+            refine!(wing; recompute_mapping=true, sort_sections=false)
+            return wing
+        end
+
+        peak(wing) = maximum(maximum(s.section_aero.cp)
+                             for s in wing.refined_sections)
+
+        wing = fine_wing(use_prior_polar=true)
+        baseline = [copy(s.section_aero.cp) for s in wing.refined_sections]
+        @test peak(wing) > 5.0
+
+        keep_endpoints!(wing)
+        @test [s.section_aero.cp for s in wing.refined_sections] == baseline
+        # Only the dropped mid section carries a cp above 5; blending the two
+        # surviving endpoints cannot reach one.
+        @test peak(wing) > 5.0
+
+        @test peak(keep_endpoints!(fine_wing(use_prior_polar=false))) <= 5.0 + 1e-9
+
+        bare = fine_wing(use_prior_polar=true, tables=false)
+        @test all(isnothing(s.section_aero) for s in bare.refined_sections)
+        for (section, scale) in zip(bare.unrefined_sections, scales)
+            section.section_aero = surface_table(scale)
+        end
+        refine!(bare; recompute_mapping=true, sort_sections=false)
+        @test all(!isnothing(s.section_aero) for s in bare.refined_sections)
+        @test peak(bare) > 5.0
+    end
 
     @testset "Refined panel mapping" begin
         # Test that refined panel mapping actually maps each panel to its closest unrefined panel