OpenSourceAWE · 1-Bort-1 · Sep 24, 2026 · Sep 24, 2026 · Sep 24, 2026 · Sep 28, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -4,12 +4,34 @@
 
 ### Added
 
+- `VSMSolution` holds everything the removed `solve` dictionary held, filled by every
+  `solve!`: `lift`, `drag`, `side`, `cl`, `cd`, `cs`, `cl_distribution`,
+  `cd_distribution`, `cs_distribution`, `alpha_uncorrected`, `va_ref_vec`, `q_ref`, `rey`,
+  `area_all_panels`, `projected_area`, `wing_span`, `aspect_ratio_projected`,
+  `center_of_pressure` and `panel_cp_locations`. `calc_only_f_and_gamma` skips the
+  projections, span and centers of pressure.
 - `plot_geometry`, `plot_distribution`, `plot_combined_analysis`, `plot_section_polars`,
   `plot_airfoil_fit` and `plot_airfoils` take `show_title=true`; with `false` the title
   is not drawn, and still names the saved file and window.
 - The panel plots (`plot!(ax, panel)`, `plot!(ax, body)`, `plot(panel)`, `plot(body)`)
   take `border_color`, default `:black`, for the panel outlines.
 
+### Changed
+
+- BREAKING: `solve` and `calculate_results` are removed; use `solve!`, which returns the
+  solver's `VSMSolution`, and `solve!(solver, body_aero, nothing)` to start from a fresh
+  circulation as `solve` did. Keys that became other fields: `Fx`/`Fy`/`Fz` → `force`,
+  `Mx`/`My`/`Mz` → `moment`, `cfx`… → `force_coeffs`, `cmx`… → `moment_coeffs`,
+  `F_distribution` → `f_body_3D`, `M_distribution` → `m_body_3D`, `alpha_at_ac` →
+  `alpha_dist`, `alpha_geometric` → `alpha_geometric_dist`, `Rey` → `rey`.
+- BREAKING: `plot_distribution` and `plot_combined_analysis` take `VSMSolution`s instead of
+  dictionaries. `plot_combined_analysis`, `plot_polars` and `generate_polar_data` solve
+  through `solve!`, so they leave `solver.sol` at their last angle.
+- `correct_aoa` applies to the VSM model only. With LLT, results that came from `solve`
+  move by up to 0.34 % in `cl`, 0.43 % in `cd` and 4.5 % in `cmx`.
+- `center_of_pressure` is `NaN` where the line of action crosses no panel, instead of
+  `nothing` with a warning.
+
 ### Fixed
 
 - The docs define the body frame as KiteUtils' `KA` frame: x from LE to TE, y towards

diff --git a/docs/src/examples.md b/docs/src/examples.md
@@ -91,20 +91,20 @@ julia> vsm_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic
 #### Step 6: Solve using both methods
 
 ```julia
-julia> results_llt = solve(llt_solver, body_aero)
-julia> results_vsm = solve(vsm_solver, body_aero)
+julia> results_llt = solve!(llt_solver, body_aero)
+julia> results_vsm = solve!(vsm_solver, body_aero)
 ```
 
 ##### Print results comparison
 
 ```julia
 julia> println("\nLifting Line Theory Results:")
-julia> println("CL = $(round(results_llt["cl"], digits=4))")
-julia> println("CD = $(round(results_llt["cd"], digits=4))")
+julia> println("CL = $(round(results_llt.cl, digits=4))")
+julia> println("CD = $(round(results_llt.cd, digits=4))")
 julia> println("\nVortex Step Method Results:")
-julia> println("CL = $(round(results_vsm["cl"], digits=4))")
-julia> println("CD = $(round(results_vsm["cd"], digits=4))")
-julia> println("Projected area = $(round(results_vsm["projected_area"], digits=4)) m²")
+julia> println("CL = $(round(results_vsm.cl, digits=4))")
+julia> println("CD = $(round(results_vsm.cd, digits=4))")
+julia> println("Projected area = $(round(results_vsm.projected_area, digits=4)) m²")
 ```
 
 #### Step 7: Plot combined analysis

diff --git a/docs/src/functions.md b/docs/src/functions.md
@@ -98,14 +98,12 @@ CurrentModule = VortexStepMethod
 set_va!
 apparent_wind
 section_pitch_rate
-solve
 solve!
 solve_base!
 reinit!(body_aero::BodyAerodynamics{P, W, T}) where {P, W, T}
 linearize
 stability_derivatives
 trim_angle
-calculate_results
 ```
 
 ## Main Plotting Functions

diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
@@ -65,6 +65,9 @@ flow_curvature_cm
 spanwise_flow_drag
 panel_force_directions
 prescribed_va_directions
+analysis_fields!
+find_center_of_pressure!
+compute_panel_center_of_pressures!
 panel_moment
 panel_couple_force
 panel_loads

diff --git a/docs/src/tips_and_tricks.md b/docs/src/tips_and_tricks.md
@@ -37,7 +37,7 @@ This approach is useful for:
 If running the example `ram_air_kite.jl` fails, try to run the `cleanup.jl` script and then try again. Background: this example caches the calculated polars. Reading cached polars can fail after an update.
 
 ## Output formats
-Currently, the `solve!()` function returns the results as [`VSMSolution`](@ref) struct. The function solve() returns a dictionary with the results. The `solve!()` function is faster, and the `solve()` contains many more entries, therefore the first function is good for integration in dynamic models and the second one better suited for aerodynamic analysis.
+The `solve!()` function returns the results as a [`VSMSolution`](@ref) struct, which it overwrites on every call. Keep a result from one solve past the next by solving it with a separate `Solver`.
 
 ## Performance
 Calling `reinit!(body_aero; init_aero=false)` is very fast. After calling `unrefined_deform!(wing, theta_angles, delta_angles)`, you have to run `reinit!(body_aero; init_aero=false)` to apply the deformed wing to the body aerodynamics. This is in turn necessary for the linearization from deformation to aerodynamic coefficients for RAM-air kites.

diff --git a/examples/V3_kite.jl b/examples/V3_kite.jl
@@ -76,26 +76,7 @@ sideslip_deg = settings.condition.beta
 yaw_rate = settings.condition.yaw_rate
 
 set_va!(body_aero, settings)
-results = VortexStepMethod.solve(solver, body_aero; log=true)
-
-PLOT && plot_polars(
-    solvers,
-    bodies,
-    labels,
-    literature_path_list=literature_paths,
-    angle_range=range(-5, 25, length=31),
-    angle_type="angle_of_attack",
-    angle_of_attack=angle_of_attack_deg,
-    side_slip=sideslip_deg,
-    va=va,
-    title="$(wing.n_panels)_panels_$(wing.spanwise_distribution)_from_yaml_settings",
-    save_path=OUTPUT_DIR,
-    is_save=false || SAVE_ALL,
-    is_show=true,
-    use_tex=USE_TEX,
-    show_moments=false,
-    cl_over_cd=true
-)
+results = solve!(solver, body_aero; log=true)
 
 # Plotting geometry
 PLOT && plot_geometry(
@@ -124,6 +105,25 @@ PLOT && plot_distribution(
     use_tex=USE_TEX
 )
 
+PLOT && plot_polars(
+    solvers,
+    bodies,
+    labels,
+    literature_path_list=literature_paths,
+    angle_range=range(-5, 25, length=31),
+    angle_type="angle_of_attack",
+    angle_of_attack=angle_of_attack_deg,
+    side_slip=sideslip_deg,
+    va=va,
+    title="$(wing.n_panels)_panels_$(wing.spanwise_distribution)_from_yaml_settings",
+    save_path=OUTPUT_DIR,
+    is_save=false || SAVE_ALL,
+    is_show=true,
+    use_tex=USE_TEX,
+    show_moments=false,
+    cl_over_cd=true
+)
+
 # --- Beta sweep ---
 PLOT && plot_polars(
     solvers,

diff --git a/examples/bench.jl b/examples/bench.jl
@@ -45,7 +45,6 @@ llt_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_
 vsm_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=VSM)
 
 # Step 5: Solve using both methods
-results_vsm = solve(vsm_solver, body_aero, nothing)
 sol = solve!(vsm_solver, body_aero, nothing)
 results_vsm_base = solve_base!(vsm_solver, body_aero, nothing)
 println("Rectangular wing, solve_base!:")
@@ -54,8 +53,6 @@ println("Rectangular wing, solve_base!:")
 # time Julia:   0.35 ms Ryzen 7950x
 println("Rectangular wing, solve!:")
 @time sol = solve!(vsm_solver, body_aero, nothing)
-println("Rectangular wing, solve:")
-@time solve(vsm_solver, body_aero, nothing)
 
 # Create wing geometry (convert-then-load: obj -> per-section NeuralFoil POLAR_MATRICES)
 ram_yaml = obj_to_yaml(

diff --git a/examples/billowing.jl b/examples/billowing.jl
@@ -94,15 +94,15 @@ set_va!(body_aero_flat, va_vec)
 set_va!(body_aero_bill, va_vec)
 
 # --- Solve and compare ---
-results_flat = VortexStepMethod.solve(
+results_flat = solve!(
     solver_flat, body_aero_flat; log=true)
-results_bill = VortexStepMethod.solve(
+results_bill = solve!(
     solver_bill, body_aero_bill; log=true)
 
-println("\nFlat wing: CL=$(round(results_flat["cl"]; digits=4)), " *
-        "CD=$(round(results_flat["cd"]; digits=4))")
-println("Billowed:  CL=$(round(results_bill["cl"]; digits=4)), " *
-        "CD=$(round(results_bill["cd"]; digits=4))")
+println("\nFlat wing: CL=$(round(results_flat.cl; digits=4)), " *
+        "CD=$(round(results_flat.cd; digits=4))")
+println("Billowed:  CL=$(round(results_bill.cl; digits=4)), " *
+        "CD=$(round(results_bill.cd; digits=4))")
 
 if PLOT
     # Plot geometry (flat wing)

diff --git a/examples/obj_to_yaml_kite.jl b/examples/obj_to_yaml_kite.jl
@@ -75,7 +75,7 @@ VortexStepMethod.reinit!(body_aero)
 solver = Solver(wing.n_panels, wing.n_unrefined_sections;
     aerodynamic_model_type=VSM, rtol=1e-5, solver_type=LOOP)
 set_va!(body_aero, apparent_wind(deg2rad(8), 0.0, va))
-results = VortexStepMethod.solve(solver, body_aero; log=true)
+results = solve!(solver, body_aero; log=true)
 
 if PLOT
     plot_geometry(body_aero, "Ram air kite (converted from .obj)"; is_show=true,

diff --git a/examples/pyramid_model.jl b/examples/pyramid_model.jl
@@ -26,31 +26,14 @@ sideslip_deg = vsm_settings.condition.beta
 yaw_rate = vsm_settings.condition.yaw_rate
 
 # Run the solver
-results = VortexStepMethod.solve(solver, body_aero; log=true)
+results = solve!(solver, body_aero; log=true)
 
 # Using plotting modules, to create more comprehensive plots
 PLOT = true
 SAVE_ALL = false
 USE_TEX = false
 OUTPUT_DIR = joinpath(dirname(@__DIR__), "output")
 
-# Plotting polars
-PLOT && plot_polars(
-    [solver],
-    [body_aero],
-    ["VSM Pyramid Model"],
-    angle_range=range(-5, 25, length=30),
-    angle_type="angle_of_attack",
-    angle_of_attack=angle_of_attack_deg,
-    side_slip=sideslip_deg,
-    va=va,
-    title="$(wing.n_panels)_panels_$(wing.spanwise_distribution)_pyramid_model",
-    save_path=OUTPUT_DIR,
-    is_save=false || SAVE_ALL,
-    is_show=true,
-    use_tex=USE_TEX
-)
-
 # Plotting geometry
 PLOT && plot_geometry(
     body_aero,
@@ -77,4 +60,21 @@ PLOT && plot_distribution(
     use_tex=USE_TEX
 )
 
+# Plotting polars
+PLOT && plot_polars(
+    [solver],
+    [body_aero],
+    ["VSM Pyramid Model"],
+    angle_range=range(-5, 25, length=30),
+    angle_type="angle_of_attack",
+    angle_of_attack=angle_of_attack_deg,
+    side_slip=sideslip_deg,
+    va=va,
+    title="$(wing.n_panels)_panels_$(wing.spanwise_distribution)_pyramid_model",
+    save_path=OUTPUT_DIR,
+    is_save=false || SAVE_ALL,
+    is_show=true,
+    use_tex=USE_TEX
+)
+
 nothing
diff --git a/examples/readme_figure.jl b/examples/readme_figure.jl
@@ -13,7 +13,7 @@ labels = ["VSM Julia", "CFD Re=5e5", "CFD Re=10e5", "VSM Python Re=5e5",
     "Wind tunnel Re=5e5"]
 
 set_va!(body_aero, settings)
-results = VortexStepMethod.solve(solver, body_aero)
+results = solve!(solver, body_aero)
 
 plot_combined_analysis(solver, body_aero, results;
     labels,

diff --git a/examples/rectangular_wing.jl b/examples/rectangular_wing.jl
@@ -49,19 +49,19 @@ llt_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_
 vsm_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=VSM)
 
 # Step 5: Solve using both methods
-results_llt = solve(llt_solver, body_aero)
-results_vsm = solve(vsm_solver, body_aero)
-@time solve(llt_solver, body_aero)
-@time solve(vsm_solver, body_aero)
+results_llt = solve!(llt_solver, body_aero)
+results_vsm = solve!(vsm_solver, body_aero)
+@time solve!(llt_solver, body_aero)
+@time solve!(vsm_solver, body_aero)
 
 # Print results comparison
 println("\nLifting Line Theory Results:")
-println("CL = $(round(results_llt["cl"], digits=4))")
-println("CD = $(round(results_llt["cd"], digits=4))")
+println("CL = $(round(results_llt.cl, digits=4))")
+println("CD = $(round(results_llt.cd, digits=4))")
 println("\nVortex Step Method Results:")
-println("CL = $(round(results_vsm["cl"], digits=4))")
-println("CD = $(round(results_vsm["cd"], digits=4))")
-println("Projected area = $(round(results_vsm["projected_area"], digits=4)) m²")
+println("CL = $(round(results_vsm.cl, digits=4))")
+println("CD = $(round(results_vsm.cd, digits=4))")
+println("Projected area = $(round(results_vsm.projected_area, digits=4)) m²")
 
 # Step 6: Plot geometry
 PLOT && plot_geometry(

diff --git a/examples/stall_model.jl b/examples/stall_model.jl
@@ -77,9 +77,9 @@ PLOT && plot_geometry(
 )
 
 # Solving and plotting distributions
-results = solve(vsm_solver, body_aero)
-@time results_with_stall = solve(VSM_with_stall_correction, body_aero)
-@time results_with_stall = solve(VSM_with_stall_correction, body_aero)
+results = solve!(vsm_solver, body_aero)
+@time results_with_stall = solve!(VSM_with_stall_correction, body_aero)
+@time results_with_stall = solve!(VSM_with_stall_correction, body_aero)
 
 CAD_y_coordinates = [panel.aero_center[2] for panel in body_aero.panels]