diff --git a/bin/run_julia b/bin/run_julia
index 204a44c9..8408a0c7 100755
--- a/bin/run_julia
+++ b/bin/run_julia
@@ -19,7 +19,7 @@ done
 if [[ $# -gt 0 ]]; then
     JULIA_ARGS=("$@")
 else
-    JULIA_ARGS=(-i -e 'using VortexStepMethod; function menu(); include("examples/menu.jl"); end; function menu_cp(); include("examples_cp/menu_cp.jl"); end')
+    JULIA_ARGS=()
 fi
 
 julia --project "${JULIA_ARGS[@]}"
diff --git a/ext/VortexStepMethodControlPlotsExt.jl b/ext/VortexStepMethodControlPlotsExt.jl
index 6042583c..dba0be96 100644
--- a/ext/VortexStepMethodControlPlotsExt.jl
+++ b/ext/VortexStepMethodControlPlotsExt.jl
@@ -793,8 +793,9 @@ function VortexStepMethod.plot_polar_data(body_aero::BodyAerodynamics,
             ax = fig.add_subplot(1, 3, idx, projection="3d")
 
             # Create interpolation matrix
-            interp_matrix = zeros(length(alphas), length(delta_tes))
-            interp_matrix .= [interp(alpha, delta_te) for alpha in alphas, delta_te in delta_tes]
+            interp_matrix = zeros(length(delta_tes), length(alphas))
+            interp_matrix .= [interp(alpha, delta_te)
+                              for delta_te in delta_tes, alpha in alphas]
             X = collect(delta_tes) .+ zeros(length(alphas))'
             Y = collect(alphas)' .+ zeros(length(delta_tes))
 
diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl
index 75d8ae61..7f3b7a42 100644
--- a/ext/VortexStepMethodMakieExt.jl
+++ b/ext/VortexStepMethodMakieExt.jl
@@ -920,7 +920,7 @@ function VortexStepMethod.plot_polar_data(body_aero::BodyAerodynamics,
 
             # Create interpolation matrix
             interp_matrix = [interp(alpha, delta_te)
-                             for alpha in alphas, delta_te in delta_tes]
+                             for delta_te in delta_tes, alpha in alphas]
 
             # Create wireframe
             wireframe!(ax, delta_tes, alphas, interp_matrix;
diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl
index 62d96c72..da9cb9ea 100644
--- a/src/VortexStepMethod.jl
+++ b/src/VortexStepMethod.jl
@@ -29,7 +29,7 @@ using Xfoil
 export SolverSettings, VSMSettings, WingSettings
 export ObjWing, Section, Wing, refine!, reinit!
 export BodyAerodynamics
-export Solver, VSMSolution, linearize, solve, solve!, solve_base!
+export Solver, VSMSolution, linearize, solve, solve!, solve_base!, calc_forces!
 export calculate_results
 export add_section!, set_va!
 export calculate_projected_area, calculate_span
diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index 8ca7c003..7f21ed28 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -315,28 +315,29 @@ end
 )
     size(va_input) == (n_panels, 3) ||
         throw(ArgumentError("'va' must be shape (3,) or ($(n_panels), 3); got $(size(va_input))"))
-
-    T = promote_type(eltype(va_input),
-                     isnothing(panel_areas) ? Float64 : eltype(panel_areas))
-    areas = if isnothing(panel_areas)
-        ones(T, n_panels)
-    else
+    if !isnothing(panel_areas)
         length(panel_areas) == n_panels ||
             throw(ArgumentError("panel_areas must be shape ($(n_panels),), got length $(length(panel_areas))"))
-        T.(panel_areas)
     end
 
-    total_area = sum(areas)
-    total_area > 0.0 || throw(ArgumentError("Total panel area must be positive."))
-
+    T = promote_type(eltype(va_input),
+                     isnothing(panel_areas) ? Float64 : eltype(panel_areas))
+    total_area = zero(T)
     weighted_speed_sq = zero(T)
     direction = zeros(MVector{3, T})
     @inbounds for i in 1:n_panels
-        @views va_i = va_input[i, :]
-        speed_i = norm(va_i)
-        weighted_speed_sq += areas[i] * speed_i^2
-        direction .+= areas[i] .* va_i
+        area_i = isnothing(panel_areas) ? one(T) : T(panel_areas[i])
+        va1 = va_input[i, 1]
+        va2 = va_input[i, 2]
+        va3 = va_input[i, 3]
+        speed_i = sqrt(va1^2 + va2^2 + va3^2)
+        total_area += area_i
+        weighted_speed_sq += area_i * speed_i^2
+        direction[1] += area_i * va1
+        direction[2] += area_i * va2
+        direction[3] += area_i * va3
     end
+    total_area > 0.0 || throw(ArgumentError("Total panel area must be positive."))
 
     reference_speed = sqrt(weighted_speed_sq / total_area)
     direction_norm = norm(direction)
diff --git a/src/solver.jl b/src/solver.jl
index 7f8430a0..ebf021e9 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -43,6 +43,7 @@ Struct for storing the solution of the [solve!](@ref) function. Must contain all
     _chord_dist::Vector{T} = zeros(T, P)
     ### end of private vectors
     width_dist::Vector{T} = zeros(T, P)
+    panel_area_dist::Vector{T} = zeros(T, P)
     alpha_dist::Vector{T} = zeros(T, P)
     alpha_geometric_dist::Vector{T} = zeros(T, P)
     cl_dist::Vector{T} = zeros(T, P)
@@ -74,6 +75,7 @@ Struct for storing the solution of the [solve!](@ref) function. Must contain all
     va_unrefined_dist::Vector{MVector{3, T}} = [zeros(MVector{3, T}) for _ in 1:U]
     chord_unrefined_dist::MVector{U, T} = zeros(MVector{U, T})
     width_unrefined_dist::MVector{U, T} = zeros(MVector{U, T})
+    unrefined_count_dist::Vector{Int} = zeros(Int, U)
     solver_status::SolverStatus = FAILURE
 end
 
@@ -237,6 +239,20 @@ function solve!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma_dist
         solver.sol.gamma_distribution = gamma_new
     end
 
+    return calc_forces!(solver, body_aero; reference_point, moment_frac)
+end
+
+"""
+    calc_forces!(solver::Solver, body_aero::BodyAerodynamics;
+                 reference_point=solver.reference_point, moment_frac=0.1)
+
+Assemble aerodynamic forces and moments from the circulation already converged
+by [`solve_base!`](@ref) and stored in `solver`. Split out of [`solve!`](@ref).
+"""
+function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
+        reference_point=solver.reference_point, moment_frac=0.1) where {P, U, T}
+    gamma_new = solver.lr.gamma_new
+
     # Initialize arrays
     cl_dist = solver.sol.cl_dist
     cd_dist = solver.sol.cd_dist
@@ -322,7 +338,7 @@ function solve!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma_dist
 
     # Initialize result arrays
     area_all_panels = zero(T)
-    panel_areas = zeros(T, length(panels))
+    panel_areas = solver.sol.panel_area_dist
 
     # Get wing properties
     spanwise_direction = body_aero.wings[1].spanwise_direction
@@ -428,6 +444,7 @@ function solve!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma_dist
         va_unrefined_dist = solver.sol.va_unrefined_dist
         chord_unrefined_dist = solver.sol.chord_unrefined_dist
         width_unrefined_dist = solver.sol.width_unrefined_dist
+        unrefined_count_dist = solver.sol.unrefined_count_dist
 
         # Zero all unrefined arrays
         moment_unrefined_dist .= 0.0
@@ -444,13 +461,12 @@ function solve!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma_dist
         end
         chord_unrefined_dist .= 0.0
         width_unrefined_dist .= 0.0
+        fill!(unrefined_count_dist, 0)
 
         panel_idx = 1
         unrefined_idx = 1
         for wing in body_aero.wings
             if wing.n_unrefined_sections > 0
-                # Accumulate values from refined panels to unrefined sections
-                unrefined_section_counts = zeros(Int, wing.n_unrefined_sections)
                 for local_panel_idx in 1:wing.n_panels
                     panel = body_aero.panels[panel_idx]
                     original_section_idx = wing.refined_panel_mapping[local_panel_idx]
@@ -472,7 +488,7 @@ function solve!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma_dist
                     chord_unrefined_dist[target_unrefined_idx] += panel.chord
                     width_unrefined_dist[target_unrefined_idx] += panel.width
 
-                    unrefined_section_counts[original_section_idx] += 1
+                    unrefined_count_dist[target_unrefined_idx] += 1
                     panel_idx += 1
                 end
 
@@ -480,8 +496,8 @@ function solve!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma_dist
                 # moment_coeff_unrefined_dist stay summed (extensive).
                 for i in 1:wing.n_unrefined_sections
                     target_unrefined_idx = unrefined_idx + i - 1
-                    if unrefined_section_counts[i] > 0
-                        count = unrefined_section_counts[i]
+                    if unrefined_count_dist[target_unrefined_idx] > 0
+                        count = unrefined_count_dist[target_unrefined_idx]
                         moment_unrefined_dist[target_unrefined_idx] /= count
                         cl_unrefined_dist[target_unrefined_idx] /= count
                         cd_unrefined_dist[target_unrefined_idx] /= count
diff --git a/test/plotting/test_plotting.jl b/test/plotting/test_plotting.jl
index f02f52ee..b779fe15 100644
--- a/test/plotting/test_plotting.jl
+++ b/test/plotting/test_plotting.jl
@@ -245,6 +245,16 @@ end
         @test fig !== nothing
     end
 
+    fig_rect = plot_polar_data(body_aero;
+        alphas=collect(deg2rad.(-5:1.0:15)),
+        delta_tes=collect(deg2rad.(-3:1.0:5)),
+        is_show=false)
+    if backend == "Makie"
+        @test fig_rect isa Figure
+    else
+        @test fig_rect !== nothing
+    end
+
     # Tests for both backends
     body_aero_empty = create_body_aero()
     empty!(body_aero_empty.panels)
diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl
index 244dab03..d74769af 100644
--- a/test/solver/test_solver.jl
+++ b/test/solver/test_solver.jl
@@ -68,3 +68,26 @@ end
         rm(settings_file; force=true)
     end
 end
+
+calc_forces_allocs(solver, body_aero) =
+    (calc_forces!(solver, body_aero); @allocated calc_forces!(solver, body_aero))
+
+@testset "calc_forces! is zero-alloc" begin
+    settings_file = create_temp_wing_settings(
+        "solver", "solver_test_wing.yaml";
+        alpha=5.0, beta=0.0, wind_speed=10.0,
+    )
+    try
+        settings = VSMSettings(settings_file)
+        wing = Wing(settings)
+        refine!(wing)
+        body_aero = BodyAerodynamics([wing])
+        solver = Solver(body_aero, settings)
+        set_va!(body_aero, [10.0, 0.0, 0.0])
+        solve!(solver, body_aero)
+
+        @test calc_forces_allocs(solver, body_aero) == 0
+    finally
+        rm(settings_file; force=true)
+    end
+end