diff --git a/README.md b/README.md
index 7a555f90..42bf46c5 100644
--- a/README.md
+++ b/README.md
@@ -89,7 +89,7 @@ Three kinds of input data is needed:
   - kite wing: model of polars included, n sections to define
 
 - The airflow and turn rate:
-  - `v_app` vector and `omega` (turn rate) vector in Kite Body (KB) reference frame
+  - `va_vec` vector and `omega` (turn rate) vector in Kite Body (KB) reference frame
 
 - The configuration:
   - how many panels  
@@ -106,7 +106,7 @@ A whole run — the flight condition, each wing and the solver — is configured
 n_panels = 20          # Number of panels
 span = 20.0            # Wing span [m]
 chord = 1.0            # Chord length [m]
-v_a = 20.0             # Magnitude of inflow velocity [m/s]
+va = 20.0              # Magnitude of inflow velocity [m/s]
 density = 1.225        # Air density [kg/m³]
 alpha_deg = 30.0       # Angle of attack [degrees]
 alpha = deg2rad(alpha_deg)
@@ -131,8 +131,8 @@ refine!(wing)
 body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
-vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
-set_va!(body_aero, vel_app)
+va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+set_va!(body_aero, va_vec)
 ```
 It is possible to import the wing geometry using an `.obj` file as shown in the example `ram_air_kite.jl`. During the import the polars are calculated automatically using XFoil. This approach is valid for rigid wings and ram-air kites, but not for leading edge inflatable kites.
 
diff --git a/docs/src/examples.md b/docs/src/examples.md
index 2c37cfbe..f9d2fa43 100644
--- a/docs/src/examples.md
+++ b/docs/src/examples.md
@@ -35,7 +35,7 @@ julia> using VortexStepMethod
 julia> n_panels = 20          # Number of panels
 julia> span = 20.0            # Wing span [m]
 julia> chord = 1.0            # Chord length [m]
-julia> v_a = 20.0             # Magnitude of inflow velocity [m/s]
+julia> va = 20.0              # Magnitude of inflow velocity [m/s]
 julia> alpha_deg = 30.0       # Angle of attack [degrees]
 julia> alpha = deg2rad(alpha_deg)
 ```
@@ -77,8 +77,8 @@ multiple wings.
 ###### Set inflow conditions
 
 ```julia
-julia> vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
-julia> set_va!(body_aero, vel_app, [0, 0, 0.1])
+julia> va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+julia> set_va!(body_aero, va_vec, [0, 0, 0.1])
 ```
 
 #### Step 5: Initialize solvers for both LLT and VSM methods
@@ -118,7 +118,7 @@ julia> plot_combined_analysis(
            solver_label=["LLT", "VSM"],
            angle_range=angle_range,
            angle_type="angle_of_attack",
-           v_a=v_a,
+           v_a=va,
            title="Rectangular Wing",
            is_show=true,
        )
diff --git a/docs/src/index.md b/docs/src/index.md
index 2e594ece..fc1cfb2f 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -94,7 +94,7 @@ Three kinds of input data is needed:
   - kite wing: model of polars included, n sections to define
 
 - The airflow and turn rate:
-  - `v_app` vector and `omega` (turn rate) vector in Kite Body (KB) reference frame
+  - `va_vec` vector and `omega` (turn rate) vector in Kite Body (KB) reference frame
 
 - The configuration:
   - how many panels  
@@ -112,7 +112,7 @@ A whole run — the flight condition, each wing and the solver — is configured
 n_panels = 20          # Number of panels
 span = 20.0            # Wing span [m]
 chord = 1.0            # Chord length [m]
-v_a = 20.0             # Magnitude of inflow velocity [m/s]
+va = 20.0              # Magnitude of inflow velocity [m/s]
 density = 1.225        # Air density [kg/m³]
 alpha_deg = 30.0       # Angle of attack [degrees]
 alpha = deg2rad(alpha_deg)
@@ -137,8 +137,8 @@ refine!(wing)
 body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
-vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
-set_va!(body_aero, vel_app)
+va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+set_va!(body_aero, va_vec)
 ```
 
 It is possible to import the wing geometry using an `.obj` file as shown in the example `ram_air_kite.jl`. During the import the polars are calculated automatically, using NeuralFoil by default or XFoil as a viscous cross-check. This approach is valid for rigid wings and ram-air kites, but not for leading edge inflatable kites. See [From CAD mesh to aerodynamic model](@ref) for the full pipeline.
diff --git a/docs/src/reference_frames.md b/docs/src/reference_frames.md
index 7be95d1f..04a9f004 100644
--- a/docs/src/reference_frames.md
+++ b/docs/src/reference_frames.md
@@ -3,7 +3,7 @@
 ### Introduction
 Reference frames are needed for following purposes:
 - for creating a CAD model of the wing (or the wings)
-- for defining the apparent wind speed vector $v_a$
+- for defining the apparent wind speed vector `va_vec`
 - for calculating the lift and drag and side force coefficients
 - for calculating the resulting forces and moments
 
@@ -30,6 +30,6 @@ The turn rates $\mathrm{omega} = [\mathrm{omega_x}, \mathrm{omega_y} ,\mathrm{om
 
 ## Input and output
 - when running a simulation, the turnrate of the kite must be provided on each time step
-- the apparent wind speed vector `v_a` is defined in the **KB** reference frame
+- the apparent wind speed vector `va_vec` is defined in the **KB** reference frame
 - the resulting forces are defined in the **KB** reference frame
 - the **CL**, **CD**, **CS** and the resulting moments and moment coefficients are defined in the **KB** reference frame
\ No newline at end of file
diff --git a/examples/V3_kite.jl b/examples/V3_kite.jl
index 4a7b831c..e3edb46d 100644
--- a/examples/V3_kite.jl
+++ b/examples/V3_kite.jl
@@ -141,7 +141,7 @@ PLOT && plot_distribution(
     [body_y_coordinates],
     [results],
     ["VSM"];
-    title="CAD_spanwise_distributions_alpha_$(round(angle_of_attack_deg, digits=1))_delta_$(round(sideslip_deg, digits=1))_yaw_$(round(yaw_rate, digits=1))_v_a_$(round(wind_speed, digits=1))",
+    title="CAD_spanwise_distributions_alpha_$(round(angle_of_attack_deg, digits=1))_delta_$(round(sideslip_deg, digits=1))_yaw_$(round(yaw_rate, digits=1))_va_$(round(wind_speed, digits=1))",
     save_path=OUTPUT_DIR,
     is_save=false || SAVE_ALL,
     is_show=true,
diff --git a/examples/V3_neuralfoil.jl b/examples/V3_neuralfoil.jl
index 88064a5b..25795625 100644
--- a/examples/V3_neuralfoil.jl
+++ b/examples/V3_neuralfoil.jl
@@ -69,7 +69,7 @@ nf_yaml = obj_to_yaml(OBJ_PATH, gen_dir; n_sections=N_SLICES, Re=RE,
     rotation=ROTATION, wrap_method=WRAP, aero_solver=NF_SOLVER, verbose=true)
 
 # Flight conditions
-v_a = 10.0
+va = 10.0
 angle_range = range(-5, 25, length=31)
 
 # Load settings and create wing with CFD polars
@@ -105,7 +105,7 @@ fig = plot_polars(
      "Wind tunnel (Poland 2025)"];
     literature_path_list=literature_paths,
     angle_range,
-    v_a,
+    v_a=va,
     title="TU Delft V3 Kite: CFD vs NeuralFoil (Re=$RE)",
     is_save=false,
 )
diff --git a/examples/bench.jl b/examples/bench.jl
index ca9414ec..be28e2d3 100644
--- a/examples/bench.jl
+++ b/examples/bench.jl
@@ -12,7 +12,7 @@ using VortexStepMethod: solve_base!
 n_panels = 20          # Number of panels
 span = 20.0            # Wing span [m]
 chord = 1.0            # Chord length [m]
-v_a = 20.0             # Magnitude of inflow velocity [m/s]
+va = 20.0              # Magnitude of inflow velocity [m/s]
 density = 1.225        # Air density [kg/m³]
 alpha_deg = 30.0       # Angle of attack [degrees]
 alpha = deg2rad(alpha_deg)
@@ -37,8 +37,8 @@ refine!(wing)
 body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
-vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
-set_va!(body_aero, vel_app)
+va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+set_va!(body_aero, va_vec)
 
 # Step 4: Initialize solvers for both LLT and VSM methods
 llt_solver = Solver(body_aero; aerodynamic_model_type=LLT)
@@ -73,17 +73,17 @@ vsm_solver = Solver(
 )
 
 # Setting velocity conditions
-v_a = 15.0
+va = 15.0
 aoa = 15.0
 side_slip = 0.0
 yaw_rate = 0.0
 aoa_rad = deg2rad(aoa)
-vel_app = [
+va_vec = [
     cos(aoa_rad) * cos(side_slip),
     sin(side_slip),
     sin(aoa_rad)
-] * v_a
-set_va!(body_aero, vel_app)
+] * va
+set_va!(body_aero, va_vec)
 
 # Solving
 solve_base!(vsm_solver, body_aero, nothing)
diff --git a/examples/billowing.jl b/examples/billowing.jl
index 219b7e54..d6e956d3 100644
--- a/examples/billowing.jl
+++ b/examples/billowing.jl
@@ -124,9 +124,9 @@ sideslip_deg = condition_cfg["beta"]
 
 α0 = deg2rad(angle_of_attack_deg)
 β0 = deg2rad(sideslip_deg)
-va = wind_speed .* [cos(α0) * cos(β0), sin(β0), sin(α0) * cos(β0)]
-set_va!(body_aero_flat, va)
-set_va!(body_aero_bill, va)
+va_vec = wind_speed .* [cos(α0) * cos(β0), sin(β0), sin(α0) * cos(β0)]
+set_va!(body_aero_flat, va_vec)
+set_va!(body_aero_bill, va_vec)
 
 # --- Solve and compare ---
 results_flat = VortexStepMethod.solve(
diff --git a/examples/linearize_check.jl b/examples/linearize_check.jl
index 63614520..63791e5b 100644
--- a/examples/linearize_check.jl
+++ b/examples/linearize_check.jl
@@ -31,22 +31,22 @@ solver = Solver(body_aero;
     use_gamma_prev=false,
 )
 
-v_a       = 15.0
-aoa_deg   = 10.0
-aoa_rad   = deg2rad(aoa_deg)
-side_slip = 0.0
-va_b_0    = [
+va         = 15.0
+aoa_deg    = 10.0
+aoa_rad    = deg2rad(aoa_deg)
+side_slip  = 0.0
+va_vec_b_0 = [
     cos(aoa_rad) * cos(side_slip),
     sin(side_slip),
     sin(aoa_rad),
-] * v_a
-omega_b_0 = zeros(3)
-theta_0   = zeros(n_unrefined)
+] * va
+omega_b_0  = zeros(3)
+theta_0    = zeros(n_unrefined)
 
 theta_idxs = 1:n_unrefined
 va_idxs    = (n_unrefined + 1):(n_unrefined + 3)
 omega_idxs = (n_unrefined + 4):(n_unrefined + 6)
-y0         = [theta_0; va_b_0; omega_b_0]
+y0         = [theta_0; va_vec_b_0; omega_b_0]
 
 @info "Computing FiniteDiff Jacobian …"
 t_fd = @elapsed begin
@@ -97,7 +97,7 @@ n_outputs = length(output_labels)
 
 input_scales = [
     fill(0.05, n_unrefined)...,    # θ [rad] : ±0.05 rad ≈ ±2.9°
-    fill(1.0, 3)...,               # va: ±1 m/s
+    fill(1.0, 3)...,               # va_vec: ±1 m/s
     fill(0.05, 3)...,              # ω : ±0.05 rad/s
 ]
 n_sweep      = 11
@@ -109,14 +109,14 @@ last_theta = fill(NaN, n_unrefined)
 
 function solve_at!(y)
     theta = y[theta_idxs]
-    va    = y[va_idxs]
+    va_vec = y[va_idxs]
     omega = y[omega_idxs]
     if !all(theta .== last_theta)
         unrefined_deform!(wing, theta, nothing; smooth=false)
         reinit!(body_aero; init_aero=false)
         last_theta .= theta
     end
-    set_va!(body_aero, va, omega)
+    set_va!(body_aero, va_vec, omega)
     solve!(solver, body_aero; log=false)
     return [
         solver.sol.force_coeffs...,
diff --git a/examples/obj_to_yaml_kite.jl b/examples/obj_to_yaml_kite.jl
index c76531b8..91eff2f9 100644
--- a/examples/obj_to_yaml_kite.jl
+++ b/examples/obj_to_yaml_kite.jl
@@ -25,7 +25,7 @@ using LinearAlgebra
 
 PLOT = true
 USE_TEX = false
-v_a = 15.0
+va = 15.0
 project_dir = dirname(@__DIR__)
 obj_path = joinpath(project_dir, "data", "ram_air_kite", "ram_air_kite_body.obj")
 output_dir = joinpath(project_dir, "output", "ram_air_kite_converted")
@@ -73,7 +73,7 @@ body_aero = BodyAerodynamics([wing])
 VortexStepMethod.reinit!(body_aero)
 
 solver = Solver(body_aero; aerodynamic_model_type=VSM, rtol=1e-5, solver_type=LOOP)
-set_va!(body_aero, [cos(deg2rad(8)) * v_a, 0.0, sin(deg2rad(8)) * v_a])
+set_va!(body_aero, [cos(deg2rad(8)) * va, 0.0, sin(deg2rad(8)) * va])
 results = VortexStepMethod.solve(solver, body_aero; log=true)
 
 if PLOT
@@ -86,7 +86,7 @@ if PLOT
     plot_section_polars(body_aero, :cd; is_show=true)
 
     plot_polars([solver], [body_aero], ["VSM (NeuralFoil polars from .obj)"];
-        angle_range=range(-5, 20, length=26), v_a=v_a,
+        angle_range=range(-5, 20, length=26), v_a=va,
         title="Ram air kite: obj_to_yaml route", is_save=false, use_tex=USE_TEX)
 end
 
diff --git a/examples/pyramid_model.jl b/examples/pyramid_model.jl
index 445638c0..c6fe02c3 100644
--- a/examples/pyramid_model.jl
+++ b/examples/pyramid_model.jl
@@ -70,7 +70,7 @@ PLOT && plot_distribution(
     [body_y_coordinates],
     [results],
     ["VSM"];
-    title="pyramid_spanwise_distributions_alpha_$(round(angle_of_attack_deg, digits=1))_delta_$(round(sideslip_deg, digits=1))_yaw_$(round(yaw_rate, digits=1))_v_a_$(round(wind_speed, digits=1))",
+    title="pyramid_spanwise_distributions_alpha_$(round(angle_of_attack_deg, digits=1))_delta_$(round(sideslip_deg, digits=1))_yaw_$(round(yaw_rate, digits=1))_va_$(round(wind_speed, digits=1))",
     save_path=OUTPUT_DIR,
     is_save=false || SAVE_ALL,
     is_show=true,
diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl
index 2eabe7fe..dc09ad6f 100644
--- a/examples/ram_air_kite.jl
+++ b/examples/ram_air_kite.jl
@@ -11,7 +11,7 @@ using LinearAlgebra
 
 PLOT = true
 USE_TEX = false
-v_a = 15.0
+va = 15.0
 RE = 1e6
 
 # Shared boundary-layer transition settings — both backends use the e^N model.
@@ -97,7 +97,7 @@ if PLOT
         [body_xfoil, body_nf],
         ["XFoil", "NeuralFoil"];
         angle_range=range(-5, 25, length=31),
-        v_a=v_a,
+        v_a=va,
         title="Ram Air Kite: XFoil vs NeuralFoil",
         is_save=false,
         use_tex=USE_TEX
diff --git a/examples/rectangular_wing.jl b/examples/rectangular_wing.jl
index 069db9b1..11aba7e3 100644
--- a/examples/rectangular_wing.jl
+++ b/examples/rectangular_wing.jl
@@ -16,7 +16,7 @@ OUTPUT_DIR = joinpath(dirname(@__DIR__), "output")
 n_panels = 20          # Number of panels
 span = 20.0            # Wing span [m]
 chord = 1.0            # Chord length [m]
-v_a = 20.0             # Magnitude of inflow velocity [m/s]
+va = 20.0              # Magnitude of inflow velocity [m/s]
 density = 1.225        # Air density [kg/m³]
 alpha_deg = 30.0       # Angle of attack [degrees]
 alpha = deg2rad(alpha_deg)
@@ -41,8 +41,8 @@ refine!(wing)
 body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
-vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
-set_va!(body_aero, vel_app, [0, 0, 0.1])
+va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+set_va!(body_aero, va_vec, [0, 0, 0.1])
 
 # Step 4: Initialize solvers for both LLT and VSM methods
 llt_solver = Solver(body_aero; aerodynamic_model_type=LLT)
@@ -94,7 +94,7 @@ PLOT && plot_polars(
     ["LLT", "VSM"];
     angle_range,
     angle_type="angle_of_attack",
-    v_a,
+    v_a=va,
     title="Rectangular Wing Polars",
     save_path=OUTPUT_DIR,
     is_save=false || SAVE_ALL,
diff --git a/examples/stall_model.jl b/examples/stall_model.jl
index 1b675c3c..e504eff3 100644
--- a/examples/stall_model.jl
+++ b/examples/stall_model.jl
@@ -56,17 +56,17 @@ VSM_with_stall_correction = Solver(body_aero;
 )
 
 # Setting velocity conditions
-v_a = 15.0
+va = 15.0
 aoa = 17.0
 side_slip = 0.0
 yaw_rate = 0.0
 aoa_rad = deg2rad(aoa)
-vel_app = [
+va_vec = [
     cos(aoa_rad) * cos(side_slip),
     sin(side_slip),
     sin(aoa_rad)
-] * v_a
-set_va!(body_aero, vel_app)
+] * va
+set_va!(body_aero, va_vec)
 
 # Plotting geometry
 PLOT && plot_geometry(
@@ -91,7 +91,7 @@ PLOT && plot_distribution(
     [CAD_y_coordinates, CAD_y_coordinates],
     [results, results_with_stall],
     ["VSM", "VSM with stall correction"];
-    title="CAD_spanwise_distributions_alpha_$(round(aoa, digits=1))_delta_$(round(side_slip, digits=1))_yaw_$(round(yaw_rate, digits=1))_v_a_$(round(v_a, digits=1))",
+    title="CAD_spanwise_distributions_alpha_$(round(aoa, digits=1))_delta_$(round(side_slip, digits=1))_yaw_$(round(yaw_rate, digits=1))_va_$(round(va, digits=1))",
     save_path=OUTPUT_DIR,
     is_save=false || SAVE_ALL,
     is_show=true,
@@ -127,7 +127,7 @@ PLOT && plot_polars(
     angle_type="angle_of_attack",
     angle_of_attack=aoa,
     side_slip=side_slip,
-    v_a=v_a,
+    v_a=va,
     title="tutorial_testing_stall_model_n_panels_$(n_panels)_distribution_$(spanwise_distribution)",
     save_path=OUTPUT_DIR,
     is_save=false || SAVE_ALL,
diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl
index 01cefdaf..f1227f8c 100644
--- a/ext/VortexStepMethodMakieExt.jl
+++ b/ext/VortexStepMethodMakieExt.jl
@@ -606,7 +606,7 @@ Create a 3D Makie plot of wing geometry including panels and filaments.
 function create_geometry_plot_makie(body_aero::BodyAerodynamics, title,
     view_elevation, view_azimuth; zoom=0.5)
     panels = body_aero.panels
-    va = getfield(body_aero, :_va)
+    va_vec = getfield(body_aero, :_va)
 
     # Create figure
     fig = Figure(size=(1400, 1400))
@@ -651,9 +651,9 @@ function create_geometry_plot_makie(body_aero::BodyAerodynamics, title,
 
     # Plot velocity vector
     max_chord = maximum(panel.chord for panel in panels)
-    va_mag = norm(va)
-    va_vector_begin = -2 * max_chord * va / va_mag
-    va_vector_end = va_vector_begin + 1.5 * va / va_mag
+    va_unit = va_vec / norm(va_vec)
+    va_vector_begin = -2 * max_chord * va_unit
+    va_vector_end = va_vector_begin + 1.5 * va_unit
     plot_line_segment_makie!(ax, [va_vector_begin, va_vector_end], :lightblue, "va")
 
     # Set equal axes
@@ -851,28 +851,6 @@ function VortexStepMethod.plot_distribution(y_coordinates_list, results_list, la
     return fig
 end
 
-"""
-    generate_polar_data(solver, body_aero::BodyAerodynamics, angle_range;
-                        angle_type="angle_of_attack", angle_of_attack=0.0,
-                        side_slip=0.0, v_a=10.0)
-
-Generate polar data for aerodynamic analysis over a range of angles.
-
-# Arguments
-- `solver`: Aerodynamic solver object
-- `body_aero`: Wing aerodynamics struct
-- `angle_range`: Range of angles to analyze
-
-# Keyword arguments
-- `angle_type`: Type of angle variation ("angle_of_attack" or "side_slip")
-- `angle_of_attack`: Initial angle of attack [°]
-- `side_slip`: Initial side slip angle [°]
-- `v_a`: Norm of apparent wind speed [m/s]
-
-# Returns
-- Tuple of polar data array and Reynolds number
-"""
-
 """
     plot_polars(solver_list, body_aero_list, label_list;
                 literature_path_list=String[],
@@ -936,7 +914,7 @@ function VortexStepMethod.plot_polars(
     for (i, (solver, body_aero)) in enumerate(zip(solver_list, body_aero_list))
         result = VortexStepMethod.generate_polar_data(
             solver, body_aero, angle_range;
-            angle_type, angle_of_attack, side_slip, v_a
+            angle_type, angle_of_attack, side_slip, va=v_a
         )
         push!(polar_data_list, result.polar_data)
         push!(cm_data_list, (cmx=result.cmx, cmy=result.cmy,
@@ -1241,7 +1219,7 @@ function VortexStepMethod.plot_combined_analysis(
     # Use first body_aero for geometry and polar data display
     first_body = body_aeros[1]
     panels = first_body.panels
-    va = getfield(first_body, :_va)
+    va_vec = getfield(first_body, :_va)
 
     # Compute spanwise results for each solver
     results_spanwise_list = copy(results_list)
@@ -1249,13 +1227,13 @@ function VortexStepMethod.plot_combined_analysis(
         α_span = deg2rad(angle_of_attack_for_spanwise_distribution)
         β_span = deg2rad(side_slip)
         for (i, (s, ba)) in enumerate(zip(solvers, body_aeros))
-            va_old = copy(getfield(ba, :_va))
+            va_vec_old = copy(getfield(ba, :_va))
             omega_old = copy(ba.omega)
             set_va!(ba, [cos(α_span) * cos(β_span), sin(β_span),
                 sin(α_span)] * v_a)
             results_spanwise_list[i] = solve(s, ba,
                 s.sol.gamma_distribution)
-            set_va!(ba, va_old, omega_old)
+            set_va!(ba, va_vec_old, omega_old)
         end
     end
 
@@ -1295,9 +1273,9 @@ function VortexStepMethod.plot_combined_analysis(
     end
 
     max_chord = maximum(panel.chord for panel in panels)
-    va_mag = norm(va)
-    va_vector_begin = -2 * max_chord * va / va_mag
-    va_vector_end = va_vector_begin + 1.5 * va / va_mag
+    va_unit = va_vec / norm(va_vec)
+    va_vector_begin = -2 * max_chord * va_unit
+    va_vector_end = va_vector_begin + 1.5 * va_unit
     plot_line_segment_makie!(ax_geo, [va_vector_begin, va_vector_end],
         :lightblue, "va")
 
@@ -1428,7 +1406,7 @@ function VortexStepMethod.plot_combined_analysis(
             zip(solvers, body_aeros, solver_labels))
         result = VortexStepMethod.generate_polar_data(
             s, ba, angle_range;
-            angle_type, angle_of_attack, side_slip, v_a)
+            angle_type, angle_of_attack, side_slip, va=v_a)
         pd = result.polar_data
         label_re = "$lbl Re = $(round(Int64,
                      result.rey * 1e-5))e5"
diff --git a/mwes/mwe_01.jl b/mwes/mwe_01.jl
index cfc40ca5..c6f0990b 100644
--- a/mwes/mwe_01.jl
+++ b/mwes/mwe_01.jl
@@ -1,4 +1,4 @@
-# Replace va_norm_dist = norm.(eachrow(solver.sol._va_dist)) with a for loop
+# Replace va_dist = norm.(eachrow(solver.sol._va_dist)) with a for loop
 
 # Testcase that shows that the new function is equivalent to the old, allocating line of code.
 using Test
@@ -9,31 +9,31 @@ struct MockSolver
     sol::NamedTuple
 end
 
-function calc_norm_dist!(va_norm_dist, va_dist)
-    for i in 1:size(va_dist, 1)
-        va_norm_dist[i] = norm(view(va_dist, i, :))
+function calc_norm_dist!(va_dist, va_vec_dist)
+    for i in 1:size(va_vec_dist, 1)
+        va_dist[i] = norm(view(va_vec_dist, i, :))
     end
 end
 
-@testset "va_norm_dist calculation" begin
-    global va_norm_dist
+@testset "va_dist calculation" begin
+    global va_dist
 
     # Create a sample 2D array
-    sample_va_dist = [
+    sample_va_vec_dist = [
         1.0 2.0 3.0;
         4.0 5.0 6.0;
         7.0 8.0 9.0
     ]
 
     # Create a mock solver with the sample array
-    mock_solver = MockSolver((va_dist = sample_va_dist,))
+    mock_solver = MockSolver((va_vec_dist = sample_va_vec_dist,))
 
-    # Calculate va_norm_dist
-    n = @allocated va_norm_dist = norm.(eachrow(mock_solver.sol.va_dist))
+    # Calculate va_dist
+    n = @allocated va_dist = norm.(eachrow(mock_solver.sol.va_vec_dist))
     println(n)
 
-    va_norm_dist2 = zeros(3)
-    m = @allocated calc_norm_dist!(va_norm_dist2, sample_va_dist)
+    va_dist2 = zeros(3)
+    m = @allocated calc_norm_dist!(va_dist2, sample_va_vec_dist)
     println(m)
 
     # Expected results (calculated manually)
@@ -44,20 +44,20 @@ end
     ]
 
     # Test the results
-    @test length(va_norm_dist) == size(sample_va_dist, 1)
-    @test va_norm_dist ≈ expected_norms atol=1e-10
+    @test length(va_dist) == size(sample_va_vec_dist, 1)
+    @test va_dist ≈ expected_norms atol=1e-10
 
     # Test individual values
-    @test va_norm_dist[1] ≈ norm(sample_va_dist[1, :]) atol=1e-10
-    @test va_norm_dist[2] ≈ norm(sample_va_dist[2, :]) atol=1e-10
-    @test va_norm_dist[3] ≈ norm(sample_va_dist[3, :]) atol=1e-10
+    @test va_dist[1] ≈ norm(sample_va_vec_dist[1, :]) atol=1e-10
+    @test va_dist[2] ≈ norm(sample_va_vec_dist[2, :]) atol=1e-10
+    @test va_dist[3] ≈ norm(sample_va_vec_dist[3, :]) atol=1e-10
 
-    @test length(va_norm_dist2) == size(sample_va_dist, 1)
-    @test va_norm_dist2 ≈ expected_norms atol=1e-10
+    @test length(va_dist2) == size(sample_va_vec_dist, 1)
+    @test va_dist2 ≈ expected_norms atol=1e-10
 
     # Test individual values
-    @test va_norm_dist2[1] ≈ norm(sample_va_dist[1, :]) atol=1e-10
-    @test va_norm_dist2[2] ≈ norm(sample_va_dist[2, :]) atol=1e-10
-    @test va_norm_dist2[3] ≈ norm(sample_va_dist[3, :]) atol=1e-10
+    @test va_dist2[1] ≈ norm(sample_va_vec_dist[1, :]) atol=1e-10
+    @test va_dist2[2] ≈ norm(sample_va_vec_dist[2, :]) atol=1e-10
+    @test va_dist2[3] ≈ norm(sample_va_vec_dist[3, :]) atol=1e-10
 end
 nothing
\ No newline at end of file
diff --git a/mwes/mwe_03.jl b/mwes/mwe_03.jl
index fe4b5349..731c18e9 100644
--- a/mwes/mwe_03.jl
+++ b/mwes/mwe_03.jl
@@ -9,7 +9,7 @@ using VortexStepMethod, PreallocationTools
 n_panels = 20          # Number of panels
 span = 20.0            # Wing span [m]
 chord = 1.0            # Chord length [m]
-v_a = 20.0             # Magnitude of inflow velocity [m/s]
+va = 20.0              # Magnitude of inflow velocity [m/s]
 density = 1.225        # Air density [kg/m³]
 alpha_deg = 30.0       # Angle of attack [degrees]
 alpha = deg2rad(alpha_deg)
diff --git a/mwes/mwe_warntype.jl b/mwes/mwe_warntype.jl
index 1eef5db2..1e5f714c 100644
--- a/mwes/mwe_warntype.jl
+++ b/mwes/mwe_warntype.jl
@@ -25,11 +25,11 @@ add_section!(wing, [0.0, span/2, 0.0], [chord, span/2, 0.0], INVISCID)
 add_section!(wing, [0.0, -span/2, 0.0], [chord, -span/2, 0.0], INVISCID)
 refine!(wing)
 body_aero = BodyAerodynamics([wing])
-vel_app = [cos(alpha), 0.0, sin(alpha)] .* 20.0
-set_va!(body_aero, vel_app)
+va_vec = [cos(alpha), 0.0, sin(alpha)] .* 20.0
+set_va!(body_aero, va_vec)
 solver = Solver(body_aero)
 
-va_norm_dist = ones(n_panels)
+va_dist = ones(n_panels)
 va_unit_dist = ones(n_panels, 3)
 
 # Prepare args for individual functions
@@ -73,7 +73,7 @@ printstyled("\n$sep\n calculate_velocity_induced_bound_2D!\n$sep\n"; color=:cyan
 
 printstyled("\n$sep\n calculate_AIC_matrices!\n$sep\n"; color=:cyan)
 @code_warntype calculate_AIC_matrices!(
-    body_aero, VSM, 0.001, va_norm_dist, va_unit_dist)
+    body_aero, VSM, 0.001, va_dist, va_unit_dist)
 
 printstyled("\n$sep\n gamma_loop!\n$sep\n"; color=:cyan)
 @code_warntype gamma_loop!(
diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index f0dfad1b..edac27a7 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -337,7 +337,8 @@ the area-weighted mean direction.
     panel_areas::Union{Nothing, AbstractVector}=nothing
 )
     length(va_input) == 3 ||
-        throw(ArgumentError("'va' must be shape (3,) or ($(n_panels), 3); got length $(length(va_input))"))
+        throw(ArgumentError("va_vec must be shape (3,) or va_vec_dist ($(n_panels), 3); " *
+                            "got length $(length(va_input))"))
     T = eltype(va_input)
     return MVector{3, T}(va_input[1], va_input[2], va_input[3])
 end
@@ -348,7 +349,8 @@ end
     panel_areas::Union{Nothing, AbstractVector}=nothing
 )
     size(va_input) == (n_panels, 3) ||
-        throw(ArgumentError("'va' must be shape (3,) or ($(n_panels), 3); got $(size(va_input))"))
+        throw(ArgumentError("va_vec must be shape (3,) or va_vec_dist ($(n_panels), 3); " *
+                            "got $(size(va_input))"))
     if !isnothing(panel_areas)
         length(panel_areas) == n_panels ||
             throw(ArgumentError("panel_areas must be shape ($(n_panels),), got length $(length(panel_areas))"))
@@ -385,7 +387,7 @@ end
 
 """
     calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::Model, core_radius_fraction,
-                            va_norm_dist, va_unit_dist, target=body_aero.AIC)
+                            va_dist, va_unit_dist, target=body_aero.AIC)
 
 Calculate Aerodynamic Influence Coefficient matrices.
 
@@ -395,7 +397,7 @@ Returns: nothing
 """
 @inline function calculate_AIC_matrices!(body_aero::BodyAerodynamics{P, W, T}, model::Model,
                               core_radius_fraction,
-                              va_norm_dist::AbstractVector{T},
+                              va_dist::AbstractVector{T},
                               va_unit_dist::AbstractMatrix{T},
                               target::AbstractArray{T, 3}=body_aero.AIC) where {P, W, T}
     # Determine evaluation point based on model
@@ -410,19 +412,19 @@ Returns: nothing
 
     # Python parity: one shared area-weighted wake vector for all panels.
     panel_areas = [panel.chord * panel.width for panel in body_aero.panels]
-    va_distribution = zeros(T, length(body_aero.panels), 3)
+    va_vec_dist = zeros(T, length(body_aero.panels), 3)
     @inbounds for i in 1:length(body_aero.panels), k in 1:3
-        va_distribution[i, k] = va_unit_dist[i, k] * va_norm_dist[i]
+        va_vec_dist[i, k] = va_unit_dist[i, k] * va_dist[i]
     end
     wake_velocity = _compute_reference_velocity_from_distribution(
-        va_distribution,
+        va_vec_dist,
         length(body_aero.panels),
         panel_areas
     )
     wake_speed = norm(wake_velocity)
     wake_speed > 0.0 || throw(ArgumentError("Wake reference speed must be positive."))
     va_unit .= wake_velocity ./ wake_speed
-    va_norm = wake_speed
+    va = wake_speed
     
     # Calculate influence coefficients
     for jring in eachindex(body_aero.panels)
@@ -438,7 +440,7 @@ Returns: nothing
                 filaments,
                 ep,
                 evaluation_point_on_bound,
-                va_norm,
+                va,
                 va_unit,
                 one(T),
                 core_radius_fraction,
@@ -488,7 +490,7 @@ end
 """
     update_effective_angle_of_attack!(alpha_corrected, body_aero::BodyAerodynamics, gamma,
                                       core_radius_fraction, z_airf_dist, x_airf_dist,
-                                      va_dist, va_norm_dist, va_unit_dist)
+                                      va_vec_dist, va_dist, va_unit_dist)
 
 Update angle of attack at aerodynamic center for VSM method.
 
@@ -501,23 +503,23 @@ function update_effective_angle_of_attack!(alpha_corrected,
     core_radius_fraction,
     z_airf_dist,
     x_airf_dist,
+    va_vec_dist,
     va_dist,
-    va_norm_dist,
     va_unit_dist)
 
     # 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_dist,
+    calculate_AIC_matrices!(body_aero, LLT, core_radius_fraction, va_dist,
                             va_unit_dist, body_aero.AIC_aero_center)
 
-    induced_velocity = body_aero.cache[1][va_dist]
+    induced_velocity = body_aero.cache[1][va_vec_dist]
     for k in 1:3
         mul!(view(induced_velocity, :, k), view(body_aero.AIC_aero_center, :, :, k), gamma)
     end
 
     # In-place relative velocity calculation
-    relative_velocity = body_aero.cache[2][va_dist]
-    relative_velocity .= va_dist .+ induced_velocity
+    relative_velocity = body_aero.cache[2][va_vec_dist]
+    relative_velocity .= va_vec_dist .+ induced_velocity
 
     # Preallocate and compute dot products manually
     n = size(relative_velocity, 1)
@@ -745,8 +747,8 @@ end
 
 """
     calculate_results(body_aero::BodyAerodynamics, gamma_new, reference_point, density,
-                      core_radius_fraction, mu, alpha_dist, v_a_dist, chord_dist,
-                      x_airf_dist, z_airf_dist, va_dist, va_norm_dist, va_unit_dist,
+                      core_radius_fraction, mu, alpha_dist, v_rel_dist, chord_dist,
+                      x_airf_dist, z_airf_dist, va_vec_dist, va_dist, va_unit_dist,
                       panels::Vector{<:Panel}, is_only_f_and_gamma_output::Bool;
                       correct_aoa=false, flow_curvature=false,
                       is_with_viscous_drag_correction=false, v_span_dist=nothing)
@@ -769,12 +771,12 @@ function calculate_results(
     core_radius_fraction,
     mu,
     alpha_dist,
-    v_a_dist,
+    v_rel_dist,
     chord_dist,
     x_airf_dist,
     z_airf_dist,
+    va_vec_dist,
     va_dist,
-    va_norm_dist,
     va_unit_dist,
     panels::Vector{<:Panel},
     is_only_f_and_gamma_output::Bool;
@@ -811,20 +813,20 @@ function calculate_results(
             panel, alpha_dist[i])
         if flow_curvature
             cm_dist[i] += flow_curvature_cm(
-                body_aero.pitch_rate_dist[i], chord_dist[i], v_a_dist[i])
+                body_aero.pitch_rate_dist[i], chord_dist[i], v_rel_dist[i])
         end
         panel_width_dist[i] = panel.width
-        va_norm = va_norm_dist[i]
+        va = va_dist[i]
         x_norm = norm3(panel.x_airf)
         z_norm = norm3(panel.z_airf)
-        if va_norm == 0.0 || x_norm == 0.0 || z_norm == 0.0
+        if va == 0.0 || x_norm == 0.0 || z_norm == 0.0
             alpha_geometric[i] = NaN
         else
-            inv_va_norm = 1.0 / va_norm
+            inv_va = 1.0 / va
             v_tangential = -dot3(panel.x_airf, panel.va) *
-                           inv_va_norm / x_norm
+                           inv_va / x_norm
             v_normal = -dot3(panel.z_airf, panel.va) *
-                       inv_va_norm / z_norm
+                       inv_va / z_norm
             alpha_geometric[i] = atan(-v_normal, -v_tangential)
         end
     end
@@ -838,8 +840,8 @@ function calculate_results(
             core_radius_fraction,
             z_airf_dist,
             x_airf_dist,
+            va_vec_dist,
             va_dist,
-            va_norm_dist,
             va_unit_dist
         )
     else
@@ -853,40 +855,40 @@ function calculate_results(
 
     # Get wing properties and reference velocity
     spanwise_direction = body_aero.wings[1].spanwise_direction
-    va_ref_vector = MVec3(0.0, 0.0, 0.0)
+    va_ref_vec = MVec3(0.0, 0.0, 0.0)
     weighted_speed_sq = 0.0
     total_area = 0.0
     @inbounds for i in 1:n_panels
         area_i = chord_dist[i] * panel_width_dist[i]
         total_area += area_i
-        speed_i = va_norm_dist[i]
+        speed_i = va_dist[i]
         weighted_speed_sq += area_i * speed_i^2
-        va_ref_vector[1] += area_i * va_dist[i, 1]
-        va_ref_vector[2] += area_i * va_dist[i, 2]
-        va_ref_vector[3] += area_i * va_dist[i, 3]
+        va_ref_vec[1] += area_i * va_vec_dist[i, 1]
+        va_ref_vec[2] += area_i * va_vec_dist[i, 2]
+        va_ref_vec[3] += area_i * va_vec_dist[i, 3]
     end
     total_area > 0.0 || throw(ArgumentError(
         "Total panel area must be positive."))
     reference_speed = sqrt(weighted_speed_sq / total_area)
-    direction_norm = norm3(va_ref_vector)
+    direction_norm = norm3(va_ref_vec)
     if direction_norm <= 0.0
-        va_ref_vector .= (1.0, 0.0, 0.0)
+        va_ref_vec .= (1.0, 0.0, 0.0)
         direction_norm = 1.0
     end
     @inbounds for k in 1:3
-        va_ref_vector[k] = va_ref_vector[k] / direction_norm *
-                           reference_speed
+        va_ref_vec[k] = va_ref_vec[k] / direction_norm *
+                        reference_speed
     end
-    va_ref_mag = norm3(va_ref_vector)
-    va_ref_mag > 0.0 || throw(ArgumentError(
+    va_ref = norm3(va_ref_vec)
+    va_ref > 0.0 || throw(ArgumentError(
         "Reference freestream magnitude must be positive."))
     va_ref_unit = body_aero.work_vectors[1]
-    inv_va_ref = 1.0 / va_ref_mag
+    inv_va_ref = 1.0 / va_ref
     @inbounds for k in 1:3
-        va_ref_unit[k] = va_ref_vector[k] * inv_va_ref
+        va_ref_unit[k] = va_ref_vec[k] * inv_va_ref
     end
     dir_lift_ref = body_aero.work_vectors[2]
-    cross3!(dir_lift_ref, va_ref_vector, spanwise_direction)
+    cross3!(dir_lift_ref, va_ref_vec, spanwise_direction)
     dir_lift_ref_norm = norm3(dir_lift_ref)
     dir_lift_ref_norm > 0.0 || throw(ArgumentError(
         "Reference lift direction is undefined because " *
@@ -896,7 +898,7 @@ function calculate_results(
     end
     dir_side_ref = body_aero.work_vectors[3]
     cross3!(dir_side_ref, dir_lift_ref, va_ref_unit)
-    q_ref = 0.5 * density * va_ref_mag^2
+    q_ref = 0.5 * density * va_ref^2
 
     dir_lift_prescribed_va = body_aero.work_vectors[9]
     temp_vec = body_aero.work_vectors[10]
@@ -911,38 +913,38 @@ function calculate_results(
         dirs = panel_force_directions(axes, alpha_corrected[i], spanwise_unit)
         c_span = 0.0
         if is_with_viscous_drag_correction
-            viscous = spanwise_flow_drag(v_a_dist[i], v_span_dist[i], panel.chord,
+            viscous = spanwise_flow_drag(v_rel_dist[i], v_span_dist[i], panel.chord,
                 density, mu)
             cd_dist[i] += viscous.delta_cd
             c_span = viscous.c_span
         end
         loads = panel_loads(axes, dirs,
-            dynamic_pressure(density, density, v_a_dist[i]),
+            dynamic_pressure(density, density, v_rel_dist[i]),
             cl_dist[i], cd_dist[i], cm_dist[i]; c_span)
         force = loads.force
 
-        va_panel_mag = va_norm_dist[i]
-        va_panel_mag > 0.0 || throw(ArgumentError(
+        va_panel = va_dist[i]
+        va_panel > 0.0 || throw(ArgumentError(
             "Panel $i has non-positive apparent " *
             "velocity magnitude."))
-        q_panel = 0.5 * density * va_panel_mag^2
+        q_panel = 0.5 * density * va_panel^2
         cross3!(dir_lift_prescribed_va,
                 panel.va, spanwise_direction)
         normalize3!(dir_lift_prescribed_va)
 
         cross3!(temp_vec, dir_lift_prescribed_va, panel.va)
-        inv_vpm = 1.0 / va_panel_mag
+        inv_va_panel = 1.0 / va_panel
         @inbounds for k in 1:3
-            temp_vec[k] *= inv_vpm
+            temp_vec[k] *= inv_va_panel
         end
         lift_prescribed_va = dot(force, dir_lift_prescribed_va)
-        drag_prescribed_va = dot(force, panel.va) * inv_vpm
+        drag_prescribed_va = dot(force, panel.va) * inv_va_panel
         side_prescribed_va = dot(force, temp_vec)
 
         lift_wing_3D_sum += lift_prescribed_va *
             dot3(dir_lift_prescribed_va, dir_lift_ref)
         drag_wing_3D_sum += drag_prescribed_va *
-            (dot3(panel.va, va_ref_unit) * inv_vpm)
+            (dot3(panel.va, va_ref_unit) * inv_va_panel)
         side_wing_3D_sum += side_prescribed_va *
             dot3(temp_vec, dir_side_ref)
 
@@ -973,7 +975,7 @@ function calculate_results(
 
     # Calculate Reynolds number
     c_ref = body_aero.c_ref
-    reynolds_number = density * va_ref_mag * c_ref / mu
+    reynolds_number = density * va_ref * c_ref / mu
 
     force_total = body_aero.work_vectors[9]
     moment_total = body_aero.work_vectors[10]
@@ -1035,7 +1037,7 @@ function calculate_results(
         "aspect_ratio_projected" => aspect_ratio_projected,
         "Rey" => reynolds_number,
         "q_ref" => q_ref,
-        "va_ref" => va_ref_vector,
+        "va_ref" => va_ref_vec,
         "center_of_pressure" => center_of_pressure,
         "panel_cp_locations" => panel_cp_locations
     )
@@ -1050,26 +1052,27 @@ end
 
 
 """
-    set_va!(body_aero::BodyAerodynamics, va::VelVector, omega=zeros(MVec3))
+    set_va!(body_aero::BodyAerodynamics, va_vec::VelVector, omega=zeros(MVec3))
 
 Set velocity array and update wake filaments.
 
 # Arguments
 - body_aero::BodyAerodynamics: The [`BodyAerodynamics`](@ref) struct to modify
-- `va::VelVector`: Velocity vector of the apparent wind speed           [m/s]
+- `va_vec::VelVector`: Velocity vector of the apparent wind speed       [m/s]
 - `omega::VelVector`: Turn rate vector around x y and z axis            [rad/s]
 
 `omega` is also projected onto each panel's spanwise axis into
 `pitch_rate_dist`, which the solver reads when `flow_curvature` is enabled.
 """
-function set_va!(body_aero::BodyAerodynamics{P, W, T}, va::AbstractVector, omega=zeros(MVector{3, T})) where {P, W, T}
+function set_va!(body_aero::BodyAerodynamics{P, W, T}, va_vec::AbstractVector,
+                 omega=zeros(MVector{3, T})) where {P, W, T}
     n_panels = length(body_aero.panels)
-    va_distribution = zeros(T, n_panels, 3)
+    va_vec_dist = zeros(T, n_panels, 3)
     body_aero.omega .= omega
     set_pitch_rate_dist!(body_aero, omega)
 
     if all(iszero, omega)
-        va_distribution .= reshape(va, 1, 3)
+        va_vec_dist .= reshape(va_vec, 1, 3)
     else
         idx = 1
         for wing in body_aero.wings
@@ -1077,8 +1080,8 @@ function set_va!(body_aero::BodyAerodynamics{P, W, T}, va::AbstractVector, omega
 
             # Calculate velocities for each panel in this wing slice
             for j in idx:panel_end
-                omega_va = -omega × body_aero.panels[j].control_point
-                va_distribution[j, :] .= omega_va .+ va
+                omega_va_vec = -omega × body_aero.panels[j].control_point
+                va_vec_dist[j, :] .= omega_va_vec .+ va_vec
             end
             idx = panel_end + 1
         end
@@ -1086,18 +1089,18 @@ function set_va!(body_aero::BodyAerodynamics{P, W, T}, va::AbstractVector, omega
 
     # Update panel velocities
     for (i, panel) in enumerate(body_aero.panels)
-        panel.va .= va_distribution[i,:]
+        panel.va .= va_vec_dist[i,:]
     end
 
     # Update wake elements
-    frozen_wake!(body_aero, va_distribution)
-    body_aero._va .= va
+    frozen_wake!(body_aero, va_vec_dist)
+    body_aero._va .= va_vec
     body_aero.has_distributed_va = false
     return nothing
 end
 
 """
-    set_va!(body_aero::BodyAerodynamics, va_distribution::AbstractMatrix;
+    set_va!(body_aero::BodyAerodynamics, va_vec_dist::AbstractMatrix;
             pitch_rate_dist=nothing)
 
 Set a per-panel inflow distribution. `pitch_rate_dist` gives each panel's rotation
@@ -1107,10 +1110,11 @@ rates differ per section and no single body rate describes them. It is reset to
 zero when omitted, because this method takes no `omega` and a stale one would
 silently feed the `flow_curvature` moment.
 """
-function set_va!(body_aero::BodyAerodynamics, va_distribution::AbstractMatrix;
+function set_va!(body_aero::BodyAerodynamics, va_vec_dist::AbstractMatrix;
                  pitch_rate_dist=nothing)
-    size(va_distribution, 1) != length(body_aero.panels) &&
-        throw(ArgumentError("Number of rows in va distribution should be equal to number of panels."))
+    size(va_vec_dist, 1) != length(body_aero.panels) &&
+        throw(ArgumentError(
+            "Number of rows in va_vec_dist should be equal to number of panels."))
     if isnothing(pitch_rate_dist)
         body_aero.pitch_rate_dist .= 0
     else
@@ -1120,12 +1124,12 @@ function set_va!(body_aero::BodyAerodynamics, va_distribution::AbstractMatrix;
     end
 
     for (i, panel) in enumerate(body_aero.panels)
-        panel.va .= va_distribution[i, :]
+        panel.va .= va_vec_dist[i, :]
     end
 
     # Update wake elements
-    frozen_wake!(body_aero, va_distribution)
-    body_aero._va .= [mean(va_distribution[:,i]) for i in 1:3]
+    frozen_wake!(body_aero, va_vec_dist)
+    body_aero._va .= [mean(va_vec_dist[:,i]) for i in 1:3]
     body_aero.has_distributed_va = true
     return nothing
 end
@@ -1162,11 +1166,11 @@ function set_va!(body_aero::BodyAerodynamics, settings::VSMSettings)
     β = deg2rad(settings.condition.beta)
     wind_speed = settings.condition.wind_speed
     
-    va = wind_speed * [
+    va_vec = wind_speed * [
         cos(α)*cos(β),  # X_b (forward)
         sin(β),         # Y_b (right)
         sin(α)*cos(β)   # Z_b (down)
     ]
     
-    set_va!(body_aero, va)
+    set_va!(body_aero, va_vec)
 end
diff --git a/src/filament.jl b/src/filament.jl
index 6ea8bceb..c675a860 100644
--- a/src/filament.jl
+++ b/src/filament.jl
@@ -122,14 +122,14 @@ end
 
 """
     velocity_3D_trailing_vortex(vel, filament::BoundFilament, 
-                              XVP, gamma, v_a, work_vectors)
+                              XVP, gamma, va, work_vectors)
 
 Calculate induced velocity by a trailing vortex filament.
 
 # Arguments
 - `XVP`: Control point coordinates
 - `gamma`: Vortex strength
-- `v_a`: Inflow velocity magnitude
+- `va`: Inflow velocity magnitude
 - work_vectors: preallocated array of intermediate variables
 
 Reference: Rick Damiani et al. "A vortex step method for nonlinear airfoil polar data 
@@ -140,7 +140,7 @@ as implemented in KiteAeroDyn".
     filament::BoundFilament,
     XVP,
     gamma,
-    v_a,
+    va,
     work_vectors
 )
     r1 = work_vectors[2]
@@ -161,7 +161,7 @@ as implemented in KiteAeroDyn".
 
     # Cut-off radius. The perpendicular component has length |r1.r0|/|r0|, so the
     # vector itself is only needed inside the core.
-    epsilon = sqrt(4 * ALPHA0 * NU * abs(d_r1_r0) / nr0 / v_a)
+    epsilon = sqrt(4 * ALPHA0 * NU * abs(d_r1_r0) / nr0 / va)
 
     cross3!(r1Xr0, r1, r0)
 
@@ -233,10 +233,11 @@ Represents a semi-infinite vortex filament.
     initialized::Bool = false
 end
 
-function reinit!(filament::SemiInfiniteFilament{T}, x1::AbstractVector, direction::AbstractVector, vel_mag::Real, filament_direction::Real) where T
+function reinit!(filament::SemiInfiniteFilament{T}, x1::AbstractVector,
+                 direction::AbstractVector, va::Real, filament_direction::Real) where T
     filament.x1 .= x1
     filament.direction .= direction
-    filament.vel_mag = vel_mag
+    filament.vel_mag = va
     filament.filament_direction = filament_direction
     filament.initialized = true
     return nothing
@@ -244,7 +245,7 @@ end
 
 """
     velocity_3D_trailing_vortex_semiinfinite(filament::SemiInfiniteFilament, 
-                                             Vf, XVP, GAMMA, v_a, work_vectors)
+                                             Vf, XVP, GAMMA, va, work_vectors)
 
 Calculate induced velocity by a semi-infinite trailing vortex filament.
 """
@@ -254,7 +255,7 @@ function velocity_3D_trailing_vortex_semiinfinite!(
     Vf,
     XVP,
     GAMMA,
-    v_a,
+    va,
     work_vectors
 )
     r1 = work_vectors[1]
@@ -266,7 +267,7 @@ function velocity_3D_trailing_vortex_semiinfinite!(
     # the vector itself is only needed inside the core.
     d_r1_Vf = dot3(r1, Vf)
     nVf = norm3(Vf)
-    epsilon = sqrt(4 * ALPHA0 * NU * abs(d_r1_Vf) * nVf / v_a)
+    epsilon = sqrt(4 * ALPHA0 * NU * abs(d_r1_Vf) * nVf / va)
 
     cross3!(r1XVf, r1, Vf)
 
diff --git a/src/panel.jl b/src/panel.jl
index 175a2a5d..6957dd0c 100644
--- a/src/panel.jl
+++ b/src/panel.jl
@@ -633,7 +633,7 @@ end
         filaments,
         evaluation_point::MVec3,
         evaluation_point_on_bound::Bool,
-        va_norm::Float64,
+        va::Float64,
         va_unit::MVec3,
         gamma::Float64,
         core_radius_fraction::Float64,
@@ -648,7 +648,7 @@ Calculate the velocity induced by a vortex ring at a control point.
 - filaments
 - `evaluation_point`::MVec3:         Point where induced velocity is evaluated
 - `evaluation_point_on_bound`::Bool: Whether evaluation point is on bound vortex
-- `va_norm`::Float64:                Norm of apparent velocity
+- `va`::Float64:                     Norm of apparent velocity
 - `va_unit`::MVec3:                  Unit vector of apparent velocity
 - `gamma`::Float64:                  Circulation strength
 - `core_radius_fraction`::Float64:   Vortex core radius as fraction of panel width
@@ -663,7 +663,7 @@ Calculate the velocity induced by a vortex ring at a control point.
     filaments,
     evaluation_point::AbstractVector{T},
     evaluation_point_on_bound::Bool,
-    va_norm::T,
+    va::T,
     va_unit::AbstractVector{T},
     gamma::T,
     core_radius_fraction::Real,
@@ -691,7 +691,7 @@ Calculate the velocity induced by a vortex ring at a control point.
         "Filament not initialized: 2."))
     velocity_3D_trailing_vortex!(
         tempvel, f2, evaluation_point, gamma,
-        va_norm, work_vectors)
+        va, work_vectors)
     velind .+= tempvel
 
     # Filament 3: trailing filament (BoundFilament)
@@ -700,7 +700,7 @@ Calculate the velocity induced by a vortex ring at a control point.
         "Filament not initialized: 3."))
     velocity_3D_trailing_vortex!(
         tempvel, f3, evaluation_point, gamma,
-        va_norm, work_vectors)
+        va, work_vectors)
     velind .+= tempvel
 
     # Filament 4: semi-infinite trailing (SemiInfiniteFilament)
@@ -709,7 +709,7 @@ Calculate the velocity induced by a vortex ring at a control point.
         "Filament not initialized: 4."))
     velocity_3D_trailing_vortex_semiinfinite!(
         tempvel, f4, va_unit, evaluation_point, gamma,
-        va_norm, work_vectors)
+        va, work_vectors)
     velind .+= tempvel
 
     # Filament 5: semi-infinite trailing (SemiInfiniteFilament)
@@ -718,7 +718,7 @@ Calculate the velocity induced by a vortex ring at a control point.
         "Filament not initialized: 5."))
     velocity_3D_trailing_vortex_semiinfinite!(
         tempvel, f5, va_unit, evaluation_point, gamma,
-        va_norm, work_vectors)
+        va, work_vectors)
     velind .+= tempvel
 
     return nothing
diff --git a/src/panel_aerodynamics.jl b/src/panel_aerodynamics.jl
index c1e103f6..9cd2cd00 100644
--- a/src/panel_aerodynamics.jl
+++ b/src/panel_aerodynamics.jl
@@ -70,18 +70,18 @@ on section ordering.
 end
 
 """
-    section_pitch_rate(delta_va, z_airf, chord)
+    section_pitch_rate(delta_va_vec, z_airf, chord)
     section_pitch_rate(velocity_leading, velocity_trailing, z_airf, chord)
 
-Rate a section rotates about its own spanwise axis, positive nose-up. `delta_va`
+Rate a section rotates about its own spanwise axis, positive nose-up. `delta_va_vec`
 is the trailing minus leading edge apparent wind; apparent wind is
 `wind - velocity`, so that is the leading minus trailing edge velocity, hence the
 reversed order in the four-argument form. Chordwise wind variation enters here
 too. Builds a `pitch_rate_dist` for [`set_va!`](@ref) on a deforming structure,
 where no single body rate describes every section.
 """
-@inline section_pitch_rate(delta_va, z_airf, chord) =
-    ifelse(chord > 0, dot(delta_va, z_airf) / smooth_norm(chord), zero(chord))
+@inline section_pitch_rate(delta_va_vec, z_airf, chord) =
+    ifelse(chord > 0, dot(delta_va_vec, z_airf) / smooth_norm(chord), zero(chord))
 @inline section_pitch_rate(velocity_leading, velocity_trailing, z_airf, chord) =
     section_pitch_rate(velocity_leading .- velocity_trailing, z_airf, chord)
 
@@ -106,23 +106,23 @@ geometric angle still turns the force, so a lag shifts the coefficients only.
 @inline effective_alpha(alpha, deficiency) = alpha - deficiency
 
 """
-    panel_inflow(axes, va_1, va_2, v_ind, dva_1=nothing, dva_2=nothing,
+    panel_inflow(axes, va_vec_1, va_vec_2, v_ind, dva_vec_1=nothing, dva_vec_2=nothing,
                  deficiency=0)
 
 Flow a panel sees, as `(; v_eff, alpha, alpha_eff, v_span, pitch_rate)`, where
 `v_span` is the effective velocity across the span. `axes` is a
-[`panel_axes`](@ref) result. `dva_1`/`dva_2` are the sections' trailing minus
+[`panel_axes`](@ref) result. `dva_vec_1`/`dva_vec_2` are the sections' trailing minus
 leading edge apparent wind, giving the [`section_pitch_rate`](@ref); `nothing`
 leaves it zero. `deficiency` feeds [`effective_alpha`](@ref).
 """
-@inline function panel_inflow(axes, va_1, va_2, v_ind, dva_1=nothing,
-                              dva_2=nothing, deficiency=0)
+@inline function panel_inflow(axes, va_vec_1, va_vec_2, v_ind, dva_vec_1=nothing,
+                              dva_vec_2=nothing, deficiency=0)
     (; x_airf, y_airf, z_airf, chord) = axes
-    v_eff = 0.5 .* (va_1 .+ va_2) .+ v_ind
+    v_eff = 0.5 .* (va_vec_1 .+ va_vec_2) .+ v_ind
     alpha = atan(dot(v_eff, z_airf), dot(v_eff, x_airf))
     v_span = cross(v_eff, y_airf)
-    pitch_rate = isnothing(dva_1) ? zero(chord) :
-        section_pitch_rate(0.5 .* (dva_1 .+ dva_2), z_airf, chord)
+    pitch_rate = isnothing(dva_vec_1) ? zero(chord) :
+        section_pitch_rate(0.5 .* (dva_vec_1 .+ dva_vec_2), z_airf, chord)
     return (; v_eff, alpha, alpha_eff=effective_alpha(alpha, deficiency),
             v_span, pitch_rate)
 end
@@ -167,17 +167,17 @@ the panel normal at its leading and trailing edge.
     scale * width * cm * q_dyn * chord
 
 """
-    spanwise_flow_drag(v_a, v_span, chord, density, mu)
+    spanwise_flow_drag(v_rel, v_span, chord, density, mu)
 
 Viscous force increments from spanwise flow (Gaunaa et al. 2024,
 doi:10.1088/1742-6596/2767/2/022068), as `(; delta_cd, c_span)`: an addition to the
-section drag coefficient and a force coefficient along `y_airf`. Both refer to `v_a`,
+section drag coefficient and a force coefficient along `y_airf`. Both refer to `v_rel`,
 the speed normal to the span; `v_span` is the velocity along `y_airf`.
 """
-@inline function spanwise_flow_drag(v_a, v_span, chord, density, mu)
-    f0 = 0.062 * (density * v_a * chord / mu)^(-1 / 7)
-    skew_factor = (hypot(v_a, v_span) / v_a)^(5 / 7)
-    return (; delta_cd = f0 * (skew_factor - 1), c_span = f0 * v_span / v_a * skew_factor)
+@inline function spanwise_flow_drag(v_rel, v_span, chord, density, mu)
+    f0 = 0.062 * (density * v_rel * chord / mu)^(-1 / 7)
+    skew_factor = (hypot(v_rel, v_span) / v_rel)^(5 / 7)
+    return (; delta_cd = f0 * (skew_factor - 1), c_span = f0 * v_span / v_rel * skew_factor)
 end
 
 """
diff --git a/src/plotting_helpers.jl b/src/plotting_helpers.jl
index 06c8baba..bbc61a72 100644
--- a/src/plotting_helpers.jl
+++ b/src/plotting_helpers.jl
@@ -76,7 +76,7 @@ end
 """
     generate_polar_data(solver, body_aero, angle_range;
         angle_type="angle_of_attack", angle_of_attack=0.0,
-        side_slip=0.0, v_a=10.0)
+        side_slip=0.0, va=10.0)
 
 Sweep over `angle_range` (degrees), solving at each angle. Returns
 a named tuple `(polar_data, cmx, cmy, cmz, rey)` where `polar_data`
@@ -90,7 +90,7 @@ function generate_polar_data(
     angle_type="angle_of_attack",
     angle_of_attack=0.0,
     side_slip=0.0,
-    v_a=10.0
+    va=10.0
 )
     n_panels = length(body_aero.panels)
     n_angles = length(angle_range)
@@ -121,7 +121,7 @@ function generate_polar_data(
         end
 
         set_va!(body_aero,
-            [cos(α) * cos(β), sin(β), sin(α)] * v_a)
+            [cos(α) * cos(β), sin(β), sin(α)] * va)
 
         results = solve(solver, body_aero,
             gamma_distribution[i, :])
diff --git a/src/solver.jl b/src/solver.jl
index e43f143a..a2683b02 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -316,7 +316,7 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
     alpha_dist = solver.lr.alpha_dist
     alpha_corrected = solver.sol.alpha_dist
     alpha_geometric_dist = solver.sol.alpha_geometric_dist
-    v_a_dist = solver.lr.v_a_dist
+    v_rel_dist = solver.lr.v_a_dist
     panels = body_aero.panels
    
     width_dist = solver.sol.width_dist
@@ -334,7 +334,7 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
         if solver.flow_curvature
             cm_dist[i] += flow_curvature_cm(
                 body_aero.pitch_rate_dist[i], solver.sol._chord_dist[i],
-                v_a_dist[i])
+                v_rel_dist[i])
         end
         width_dist[i] = panel.width
 
@@ -344,7 +344,7 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
             va1 = solver.sol._va_dist[i,1]
             va2 = solver.sol._va_dist[i,2]
             va3 = solver.sol._va_dist[i,3]
-            va_norm = sqrt(va1^2 + va2^2 + va3^2)
+            va = sqrt(va1^2 + va2^2 + va3^2)
             x1 = solver.sol._x_airf_dist[i,1]
             x2 = solver.sol._x_airf_dist[i,2]
             x3 = solver.sol._x_airf_dist[i,3]
@@ -353,10 +353,10 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
             z2 = solver.sol._z_airf_dist[i,2]
             z3 = solver.sol._z_airf_dist[i,3]
             z_norm = sqrt(z1^2 + z2^2 + z3^2)
-            if va_norm == 0 || x_norm == 0 || z_norm == 0
+            if va == 0 || x_norm == 0 || z_norm == 0
                 alpha_geometric_dist[i] = NaN
             else
-                inv_va = -1.0 / va_norm
+                inv_va = -1.0 / va
                 vu1 = va1 * inv_va
                 vu2 = va2 * inv_va
                 vu3 = va3 * inv_va
@@ -412,13 +412,13 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
         dirs = panel_force_directions(axes, alpha_corrected[i], spanwise_unit)
         c_span = zero(T)
         if solver.is_with_viscous_drag_correction
-            viscous = spanwise_flow_drag(v_a_dist[i], solver.lr.v_span_dist[i],
+            viscous = spanwise_flow_drag(v_rel_dist[i], solver.lr.v_span_dist[i],
                 panel.chord, density, solver.mu)
             cd_dist[i] += viscous.delta_cd
             c_span = viscous.c_span
         end
         loads = panel_loads(axes, dirs,
-            dynamic_pressure(density, density, v_a_dist[i]),
+            dynamic_pressure(density, density, v_rel_dist[i]),
             cl_dist[i], cd_dist[i], cm_dist[i]; c_span)
         lift[i] = loads.lift
         drag[i] = loads.drag
@@ -437,14 +437,14 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
     end
 
     # Python parity: normalize with area-weighted reference velocity for distributed inflow.
-    va_ref_vector = _compute_reference_velocity_from_distribution(
+    va_ref_vec = _compute_reference_velocity_from_distribution(
         solver.sol._va_dist,
         length(panels),
         panel_areas
     )
-    va_ref_mag = norm(va_ref_vector)
-    va_ref_mag > 0.0 || throw(ArgumentError("Reference freestream magnitude must be positive."))
-    q_ref = 0.5 * density * va_ref_mag^2
+    va_ref = norm(va_ref_vec)
+    va_ref > 0.0 || throw(ArgumentError("Reference freestream magnitude must be positive."))
+    q_ref = 0.5 * density * va_ref^2
     moment_coeff_dist .= moment_dist ./ (q_ref * projected_area * c_ref)
 
     # Only compute unrefined arrays if there are unrefined sections
@@ -616,11 +616,11 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n
     return results
 end
 
-@inline @inbounds function calc_norm_dist!(va_norm_dist, va_dist)
-    for i in axes(va_dist, 1)
-        va_norm_dist[i] = sqrt(
-            va_dist[i,1]^2 + va_dist[i,2]^2 +
-            va_dist[i,3]^2)
+@inline @inbounds function calc_norm_dist!(va_dist, va_vec_dist)
+    for i in axes(va_vec_dist, 1)
+        va_dist[i] = sqrt(
+            va_vec_dist[i,1]^2 + va_vec_dist[i,2]^2 +
+            va_vec_dist[i,3]^2)
     end
 end
 
@@ -639,7 +639,8 @@ function solve_base!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma
                log=false) where {P, U, T}
     
     # check arguments
-    isnothing(body_aero.panels[1].va) && throw(ArgumentError("Inflow conditions are not set, use set_va!(body_aero, va)"))
+    isnothing(body_aero.panels[1].va) && throw(ArgumentError(
+        "Inflow conditions are not set, use set_va!(body_aero, va_vec)"))
     
     # Initialize variables
     panels = body_aero.panels
@@ -723,7 +724,7 @@ end
     velocity_view_x,
     velocity_view_y,
     velocity_view_z,
-    va_dist,
+    va_vec_dist,
     induced_velocity_all,
     relative_velocity_dist,
     y_airf_dist,
@@ -741,7 +742,7 @@ end
     mul!(velocity_view_y, AIC_y, gamma_in)
     mul!(velocity_view_z, AIC_z, gamma_in)
 
-    relative_velocity_dist .= va_dist .+ induced_velocity_all
+    relative_velocity_dist .= va_vec_dist .+ induced_velocity_all
     @inbounds for i in 1:n_panels
         ax = relative_velocity_dist[i,1]
         ay = relative_velocity_dist[i,2]
@@ -752,9 +753,9 @@ end
         relative_velocity_crossz[i,1] = ay*bz - az*by
         relative_velocity_crossz[i,2] = az*bx - ax*bz
         relative_velocity_crossz[i,3] = ax*by - ay*bx
-        ax = va_dist[i,1]
-        ay = va_dist[i,2]
-        az = va_dist[i,3]
+        ax = va_vec_dist[i,1]
+        ay = va_vec_dist[i,2]
+        az = va_vec_dist[i,3]
         v_acrossz_dist[i,1] = ay*bz - az*by
         v_acrossz_dist[i,2] = az*bx - ax*bz
         v_acrossz_dist[i,3] = ax*by - ay*bx
@@ -898,7 +899,7 @@ function gamma_loop!(
     relaxation_factor;
     log::Bool = true
 ) where {P, U, T}
-    va_dist = solver.sol._va_dist
+    va_vec_dist = solver.sol._va_dist
     chord_dist = solver.sol._chord_dist
     x_airf_dist = solver.sol._x_airf_dist
     y_airf_dist = solver.sol._y_airf_dist
@@ -911,10 +912,10 @@ function gamma_loop!(
     va_magw_dist             = solver.cache[1][solver.lr.v_a_dist]
     gamma                    = solver.cache[2][solver.lr.gamma_new]
     abs_gamma_new            = solver.cache[3][solver.lr.gamma_new]
-    induced_velocity_all     = solver.cache[4][va_dist]
-    relative_velocity_dist   = solver.cache[5][va_dist]
-    relative_velocity_crossz = solver.cache[6][va_dist]
-    v_acrossz_dist           = solver.cache[7][va_dist]
+    induced_velocity_all     = solver.cache[4][va_vec_dist]
+    relative_velocity_dist   = solver.cache[5][va_vec_dist]
+    relative_velocity_crossz = solver.cache[6][va_vec_dist]
+    v_acrossz_dist           = solver.cache[7][va_vec_dist]
     cl_dist                  = solver.cache[8][solver.lr.gamma_new]
     damp                     = solver.cache[9][solver.lr.gamma_new]
     damp                    .= zero(T)
@@ -946,7 +947,7 @@ function gamma_loop!(
             residual, gamma_iter, solver, panels, n_panels,
             AIC_x, AIC_y, AIC_z,
             velocity_view_x, velocity_view_y, velocity_view_z,
-            va_dist, induced_velocity_all, relative_velocity_dist,
+            va_vec_dist, induced_velocity_all, relative_velocity_dist,
             y_airf_dist, relative_velocity_crossz, v_acrossz_dist,
             z_airf_dist, x_airf_dist,
             v_normal_dist, v_tangential_dist,
@@ -968,7 +969,7 @@ function gamma_loop!(
                     residual_perturbed, gamma_perturbed, solver, panels, n_panels,
                     AIC_x, AIC_y, AIC_z,
                     velocity_view_x, velocity_view_y, velocity_view_z,
-                    va_dist, induced_velocity_all, relative_velocity_dist,
+                    va_vec_dist, induced_velocity_all, relative_velocity_dist,
                     y_airf_dist, relative_velocity_crossz, v_acrossz_dist,
                     z_airf_dist, x_airf_dist,
                     v_normal_dist, v_tangential_dist,
@@ -997,7 +998,7 @@ function gamma_loop!(
                     residual_perturbed, gamma_perturbed, solver, panels, n_panels,
                     AIC_x, AIC_y, AIC_z,
                     velocity_view_x, velocity_view_y, velocity_view_z,
-                    va_dist, induced_velocity_all, relative_velocity_dist,
+                    va_vec_dist, induced_velocity_all, relative_velocity_dist,
                     y_airf_dist, relative_velocity_crossz, v_acrossz_dist,
                     z_airf_dist, x_airf_dist,
                     v_normal_dist, v_tangential_dist,
@@ -1058,7 +1059,7 @@ function gamma_loop!(
                 velocity_view_x,
                 velocity_view_y,
                 velocity_view_z,
-                va_dist,
+                va_vec_dist,
                 induced_velocity_all,
                 relative_velocity_dist,
                 y_airf_dist,
@@ -1335,9 +1336,9 @@ function linearize(solver::Solver, body_aero::BodyAerodynamics, y::Vector{T};
             VortexStepMethod.reinit!(body_aero_c; init_aero=false)
         end
 
-        va = isnothing(va_idxs) ? MVector{3, TI}(body_aero_c._va) : y_in[va_idxs]
+        va_vec = isnothing(va_idxs) ? MVector{3, TI}(body_aero_c._va) : y_in[va_idxs]
         om = isnothing(omega_idxs) ? MVector{3, TI}(body_aero_c.omega) : y_in[omega_idxs]
-        set_va!(body_aero_c, va, om)
+        set_va!(body_aero_c, va_vec, om)
 
         solve!(solver_c, body_aero_c; kwargs...)
         solver_c.lr.converged || (n_failed[] += 1)
diff --git a/src/wake.jl b/src/wake.jl
index fdd8003a..971702bb 100644
--- a/src/wake.jl
+++ b/src/wake.jl
@@ -1,6 +1,6 @@
 
 """
-    frozen_wake(body_aero::BodyAerodynamics, va_distribution)
+    frozen_wake(body_aero::BodyAerodynamics, va_vec_dist)
 
 Update the filaments of the panels with frozen wake model.
 Uses one shared wake vector computed from area-weighted distributed inflow.
@@ -9,19 +9,20 @@ Replaces older filaments if present by checking length of filaments.
 
 # Arguments
 - `body_aero`::BodyAerodynamics: see: [`BodyAerodynamics`](@ref) 
-- `va_distribution::Matrix{Float64}`: Array of velocity vectors at each panel
+- `va_vec_dist::Matrix{Float64}`: Array of velocity vectors at each panel
 
 # Returns
 - nothing
 """
-function frozen_wake!(body_aero::BodyAerodynamics, va_distribution)
+function frozen_wake!(body_aero::BodyAerodynamics, va_vec_dist)
     n_panels = length(body_aero.panels)
-    size(va_distribution) == (n_panels, 3) ||
-        throw(ArgumentError("va_distribution must be shape ($(n_panels), 3), got $(size(va_distribution))"))
+    size(va_vec_dist) == (n_panels, 3) ||
+        throw(ArgumentError(
+            "va_vec_dist must be shape ($(n_panels), 3), got $(size(va_vec_dist))"))
 
     panel_areas = [panel.chord * panel.width for panel in body_aero.panels]
     wake_velocity = _compute_reference_velocity_from_distribution(
-        va_distribution,
+        va_vec_dist,
         n_panels,
         panel_areas
     )
diff --git a/test/bench.jl b/test/bench.jl
index c956add9..887f25a0 100644
--- a/test/bench.jl
+++ b/test/bench.jl
@@ -27,7 +27,7 @@ using LinearAlgebra
     n_panels = 20          # Number of panels
     span = 20.0            # Wing span [m]
     chord = 1.0            # Chord length [m]
-    v_a = 20.0             # Magnitude of inflow velocity [m/s]
+    va = 20.0              # Magnitude of inflow velocity [m/s]
     density = 1.225        # Air density [kg/m³]
     alpha_deg = 30.0       # Angle of attack [degrees]
     alpha = deg2rad(alpha_deg)
@@ -68,8 +68,8 @@ using LinearAlgebra
         @test result.allocs ≤ 50
     end
 
-    vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
-    set_va!(body_aero, vel_app)
+    va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+    set_va!(body_aero, va_vec)
 
     # Initialize solvers for both LLT and VSM methods
     solver = Solver(body_aero)
@@ -83,7 +83,7 @@ using LinearAlgebra
     AIC_z = similar(AIC_x)
     v_ind = zeros(3)
     point = rand(3)
-    va_norm_dist = ones(n_panels)
+    va_dist = ones(n_panels)
     va_unit_dist = ones(n_panels, 3)
     
     models = [VSM, LLT]
@@ -94,7 +94,7 @@ using LinearAlgebra
         for model in models
             for frac in core_radius_fractions
                 @testset "Model $model Core Radius Fraction $frac" begin
-                    result = @benchmark calculate_AIC_matrices!($body_aero, $model, $frac, $va_norm_dist, $va_unit_dist) samples=1 evals=1
+                    result = @benchmark calculate_AIC_matrices!($body_aero, $model, $frac, $va_dist, $va_unit_dist) samples=1 evals=1
                     @test result.allocs ≤ 30
                     @info "Model: $(model) \t Core radius fraction: $(frac) \t Allocations: $(result.allocs) \t Memory: $(result.memory)"
                 end
@@ -106,7 +106,7 @@ using LinearAlgebra
         @info "Gamma Loop"
         # Pre-allocate arrays
         gamma_new = zeros(n_panels)
-        va_dist = zeros(n_panels, 3)
+        va_vec_dist = zeros(n_panels, 3)
         chord_dist = zeros(n_panels)
         x_airf_dist = zeros(n_panels, 3)
         y_airf_dist = zeros(n_panels, 3)
@@ -114,7 +114,7 @@ using LinearAlgebra
         
         # Fill arrays with data
         for (i, panel) in enumerate(body_aero.panels)
-            va_dist[i, :] .= panel.va
+            va_vec_dist[i, :] .= panel.va
             chord_dist[i] = panel.chord
             x_airf_dist[i, :] .= panel.x_airf
             y_airf_dist[i, :] .= panel.y_airf
@@ -145,7 +145,7 @@ using LinearAlgebra
                 solver = Solver(body_aero;
                     aerodynamic_model_type=model
                 )
-                solver.sol._va_dist .= va_dist
+                solver.sol._va_dist .= va_vec_dist
                 solver.sol._chord_dist .= chord_dist
                 solver.sol._x_airf_dist .= x_airf_dist
                 solver.sol._y_airf_dist .= y_airf_dist
@@ -166,28 +166,28 @@ using LinearAlgebra
     @testset "Results Calculation" begin
         # Pre-allocate arrays
         alpha_dist = zeros(n_panels)
-        v_a_dist = zeros(n_panels)
+        v_rel_dist = zeros(n_panels)
         chord_dist = zeros(n_panels)
         x_airf_dist = zeros(n_panels, 3)
         y_airf_dist = zeros(n_panels, 3)
         z_airf_dist = zeros(n_panels, 3)
-        va_dist = zeros(n_panels, 3)
-        va_norm_dist = zeros(n_panels)
+        va_vec_dist = zeros(n_panels, 3)
+        va_dist = zeros(n_panels)
         va_unit_dist = zeros(n_panels, 3)
         reference_point = zeros(3)
         
 
-        set_va!(body_aero, vel_app)
+        set_va!(body_aero, va_vec)
         # Fill arrays with panel data to satisfy calculate_results preconditions.
         for (i, panel) in enumerate(body_aero.panels)
             chord_dist[i] = panel.chord
             x_airf_dist[i, :] .= panel.x_airf
             y_airf_dist[i, :] .= panel.y_airf
             z_airf_dist[i, :] .= panel.z_airf
-            va_dist[i, :] .= panel.va
-            va_norm_dist[i] = norm(panel.va)
-            va_unit_dist[i, :] .= va_norm_dist[i] > 0.0 ? panel.va ./ va_norm_dist[i] : [1.0, 0.0, 0.0]
-            v_a_dist[i] = va_norm_dist[i]
+            va_vec_dist[i, :] .= panel.va
+            va_dist[i] = norm(panel.va)
+            va_unit_dist[i, :] .= va_dist[i] > 0.0 ? panel.va ./ va_dist[i] : [1.0, 0.0, 0.0]
+            v_rel_dist[i] = va_dist[i]
         end
         results = @MVector zeros(3)
         
@@ -199,12 +199,12 @@ using LinearAlgebra
             1e-20,
             0.0,
             $alpha_dist,
-            $v_a_dist,
+            $v_rel_dist,
             $chord_dist,
             $x_airf_dist,
             $z_airf_dist,
+            $va_vec_dist,
             $va_dist,
-            $va_norm_dist,
             $va_unit_dist,
             $body_aero.panels,
             false
diff --git a/test/bench_solve.jl b/test/bench_solve.jl
index f9dcc4be..1ede63d8 100644
--- a/test/bench_solve.jl
+++ b/test/bench_solve.jl
@@ -11,7 +11,7 @@ using Test
 n_panels = 20          # Number of panels
 span = 20.0            # Wing span [m]
 chord = 1.0            # Chord length [m]
-v_a = 20.0             # Magnitude of inflow velocity [m/s]
+va = 20.0              # Magnitude of inflow velocity [m/s]
 density = 1.225        # Air density [kg/m³]
 alpha_deg = 30.0       # Angle of attack [degrees]
 alpha = deg2rad(alpha_deg)
@@ -34,8 +34,8 @@ refine!(wing)
 body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
-vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
-set_va!(body_aero, vel_app)
+va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+set_va!(body_aero, va_vec)
 
 # Step 4: Initialize solvers for both LLT and VSM methods
 vsm_solver = Solver(body_aero; aerodynamic_model_type=VSM)
diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl
index 037a1b6f..28121e1d 100644
--- a/test/body_aerodynamics/test_body_aerodynamics.jl
+++ b/test/body_aerodynamics/test_body_aerodynamics.jl
@@ -38,9 +38,9 @@ end
     @debug "N: $N"
     @debug "size(coord): $(size(coord))"
 
-    v_a = 20.0
+    va = 20.0
     aoa = 5.7106 * π / 180
-    v_a = [cos(aoa), 0.0, sin(aoa)] .* v_a
+    va_vec = [cos(aoa), 0.0, sin(aoa)] .* va
 
     # Create wing geometry
     core_radius_fraction = 1e-20
@@ -57,11 +57,11 @@ end
     
     refine!(wing)
     body_aero = BodyAerodynamics([wing])
-    set_va!(body_aero, v_a)
+    set_va!(body_aero, va_vec)
 
     # Calculate reference matrices using thesis functions
     controlpoints, rings, bladepanels, ringvec, coord_L =
-        create_geometry_general(coord, v_a, N, "5fil", LLT)
+        create_geometry_general(coord, va_vec, N, "5fil", LLT)
     
     # Test LLT matrices
     @testset "LLT Matrices" begin
@@ -70,7 +70,7 @@ end
             deepcopy(ringvec),
             deepcopy(controlpoints),
             deepcopy(rings),
-            deepcopy(v_a),
+            deepcopy(va_vec),
             zeros(N-1),
             nothing,  # data_airf not needed
             nothing,  # conv_crit not needed
@@ -78,13 +78,14 @@ end
         )
 
         # Calculate new matrices
-        va_norm_dist = fill(norm(v_a), length(body_aero.panels))
-        va_unit_dist = repeat(reshape(v_a ./ norm(v_a), 1, 3), length(body_aero.panels))
+        va_dist = fill(norm(va_vec), length(body_aero.panels))
+        va_unit_dist = repeat(reshape(va_vec ./ norm(va_vec), 1, 3),
+                              length(body_aero.panels))
         calculate_AIC_matrices!(
             body_aero,
             LLT,
             core_radius_fraction,
-            va_norm_dist,
+            va_dist,
             va_unit_dist
         )
         AIC_x, AIC_y, AIC_z = @views body_aero.AIC[:, :, 1], body_aero.AIC[:, :, 2], body_aero.AIC[:, :, 3]
@@ -99,13 +100,13 @@ end
     @testset "VSM Matrices" begin
         # Calculate reference matrices for VSM
         controlpoints, rings, bladepanels, ringvec, coord_L = 
-            create_geometry_general(coord, v_a, N, "5fil", VSM)
+            create_geometry_general(coord, va_vec, N, "5fil", VSM)
         
         MatrixU, MatrixV, MatrixW = thesis_induction_matrix_creation(
             deepcopy(ringvec),
             deepcopy(controlpoints),
             deepcopy(rings),
-            deepcopy(v_a),
+            deepcopy(va_vec),
             zeros(N-1),
             nothing,
             nothing,
@@ -113,13 +114,14 @@ end
         )
 
         # Calculate new matrices
-        va_norm_dist = fill(norm(v_a), length(body_aero.panels))
-        va_unit_dist = repeat(reshape(v_a ./ norm(v_a), 1, 3), length(body_aero.panels))
+        va_dist = fill(norm(va_vec), length(body_aero.panels))
+        va_unit_dist = repeat(reshape(va_vec ./ norm(va_vec), 1, 3),
+                              length(body_aero.panels))
         calculate_AIC_matrices!(
             body_aero,
             VSM,
             core_radius_fraction,
-            va_norm_dist,
+            va_dist,
             va_unit_dist
         )
         AIC_x, AIC_y, AIC_z = body_aero.AIC[:, :, 1], body_aero.AIC[:, :, 2], body_aero.AIC[:, :, 3]
@@ -154,9 +156,9 @@ end
         span = 17.0
         AR = span^2 / (π * span * max_chord / 4)
         @debug "AR: $AR"
-        v_a = 20.0
+        va = 20.0
         aoa = 5.7106 * π / 180
-        v_a = [cos(aoa), 0.0, sin(aoa)] .* v_a
+        va_vec = [cos(aoa), 0.0, sin(aoa)] .* va
     
         coord = if wing_type === :rectangular
             theta = range(-0.5, 0.5, length=N)
@@ -191,23 +193,23 @@ end
         end
         refine!(wing)
         body_aero = BodyAerodynamics([wing])
-        set_va!(body_aero, v_a)
+        set_va!(body_aero, va_vec)
 
-        return body_aero, coord, v_a, model
+        return body_aero, coord, va_vec, model
     end
 
     for model in [VSM, LLT]
         @debug "model: $model"
         for wing_type in [:rectangular, :curved, :elliptical]
             @debug "wing_type: $wing_type"
-            body_aero, coord, v_a, model = create_geometry(
+            body_aero, coord, va_vec, model = create_geometry(
                 model=model, wing_type=wing_type
             )
             
             # Generate geometry
             expected_controlpoints, expected_rings, expected_bladepanels, 
                 expected_ringvec, expected_coord_L = create_geometry_general(
-                coord, v_a, div(size(coord,1), 2), "5fil", model
+                coord, va_vec, div(size(coord,1), 2), "5fil", model
             )
 
             for i in 1:length(body_aero.panels)
@@ -298,10 +300,10 @@ end
     N = 40
     max_chord = 1.0
     span = 15.709  # AR = 20
-    v_a = 20.0
+    va = 20.0
     AR = span^2 / (π * span * max_chord / 4)
     aoa = deg2rad(5)
-    v_a = [cos(aoa), 0.0, sin(aoa)] .* v_a
+    va_vec = [cos(aoa), 0.0, sin(aoa)] .* va
     model = VSM
 
     # Setup wing geometry
@@ -322,7 +324,7 @@ end
     
     refine!(wing)
     body_aero = BodyAerodynamics([wing])
-    set_va!(body_aero, v_a)
+    set_va!(body_aero, va_vec)
 
     # Run analysis
     loop_solver = Solver(body_aero;
@@ -374,7 +376,7 @@ end
 
     # Calculate forces using uncorrected alpha
     alpha = results_NEW["alpha_uncorrected"]
-    dyn_visc = 0.5 * density * norm(v_a)^2
+    dyn_visc = 0.5 * density * norm(va_vec)^2
     n_panels = length(body_aero.panels)
     lift = zeros(n_panels)
     drag = zeros(n_panels)
@@ -404,7 +406,7 @@ end
     Atot = calculate_projected_area(wing)
 
     F_rel_ref, F_gl_ref, Ltot_ref, Dtot_ref, CL_ref, CD_ref, CS_ref = 
-        output_results(Fmag, aero_coeffs, ringvec, v_a, controlpoints, Atot)
+        output_results(Fmag, aero_coeffs, ringvec, va_vec, controlpoints, Atot)
 
     # Compare results
     @info "Comparing results"
@@ -441,12 +443,12 @@ end
         set_va!(body_aero, settings)
 
         α, β, wind_speed = deg2rad(10.0), deg2rad(5.0), 15.0
-        expected_va = wind_speed .* [cos(α)*cos(β), sin(β), sin(α)*cos(β)]
+        expected_va_vec = wind_speed .* [cos(α)*cos(β), sin(β), sin(α)*cos(β)]
 
         for p in body_aero.panels
-            @test p.va ≈ expected_va atol=1e-10
+            @test p.va ≈ expected_va_vec atol=1e-10
         end
-        @test body_aero._va ≈ expected_va atol=1e-10
+        @test body_aero._va ≈ expected_va_vec atol=1e-10
     finally
         isfile(settings_file) && rm(settings_file; force=true)
     end
@@ -455,11 +457,11 @@ end
 @testset "set_va! with distributed inflow blocks body_aero.va access" begin
     body_aero = BodyAerodynamics([inviscid_wing([0.0, 1.0, 2.0])])
 
-    va_distribution = [
+    va_vec_dist = [
         10.0 0.0 0.0
         9.0 0.0 1.0
     ]
-    set_va!(body_aero, va_distribution)
+    set_va!(body_aero, va_vec_dist)
 
     @test body_aero.has_distributed_va
     try
@@ -479,25 +481,25 @@ end
     body_aero = BodyAerodynamics([inviscid_wing([0.0, 1.0, 2.0]),
                                   inviscid_wing([10.0, 11.0, 12.0])])
 
-    va = [10.0, 0.0, 0.0]
+    va_vec = [10.0, 0.0, 0.0]
     omega = [0.0, 0.0, 1.0]
-    set_va!(body_aero, va, omega)
+    set_va!(body_aero, va_vec, omega)
 
     for panel in body_aero.panels
-        expected_va = va .+ (-omega × panel.control_point)
-        @test panel.va ≈ expected_va atol=1e-12
+        expected_va_vec = va_vec .+ (-omega × panel.control_point)
+        @test panel.va ≈ expected_va_vec atol=1e-12
     end
     @test body_aero.omega ≈ omega
     @test !body_aero.has_distributed_va
-    @test body_aero.va ≈ va
+    @test body_aero.va ≈ va_vec
 
     new_omega = [0.0, 0.0, 2.0]
-    @test body_aero._va ≈ va
+    @test body_aero._va ≈ va_vec
     body_aero.omega = new_omega
 
     for panel in body_aero.panels
-        expected_va = va .+ (-new_omega × panel.control_point)
-        @test panel.va ≈ expected_va atol=1e-12
+        expected_va_vec = va_vec .+ (-new_omega × panel.control_point)
+        @test panel.va ≈ expected_va_vec atol=1e-12
     end
     @test body_aero.omega ≈ new_omega
 end
diff --git a/test/body_aerodynamics/test_results.jl b/test/body_aerodynamics/test_results.jl
index 83aa8fc2..2d3667d1 100644
--- a/test/body_aerodynamics/test_results.jl
+++ b/test/body_aerodynamics/test_results.jl
@@ -23,7 +23,7 @@ end
     #      the Jacobian columns reflect true local sensitivity rather than
     #      numerical noise (a noise-driven Jacobian would not extrapolate).
 
-    va = [15.0, 1.0, 0.5]
+    va_vec = [15.0, 1.0, 0.5]
     theta = deg2rad.([2.0, 1.0, -1.0, -2.0])
     delta = deg2rad.([1.0, 0.5, -0.5, -1.0])
     omega = [0.0, 0.1, 0.0]
@@ -31,7 +31,7 @@ end
     fd_step = 1e-3
 
     VortexStepMethod.unrefined_deform!(ram_wing, theta, delta; smooth=false)
-    body_aero = BodyAerodynamics([ram_wing]; va, omega)
+    body_aero = BodyAerodynamics([ram_wing]; va=va_vec, omega)
     solver = Solver(body_aero;
         aerodynamic_model_type=VSM,
         is_with_artificial_damping=false,
@@ -40,7 +40,7 @@ end
         solver_type=NONLIN,
     )
 
-    base_inputs = [theta; va; omega; delta]
+    base_inputs = [theta; va_vec; omega; delta]
     jac, lin_res, lin_converged = VortexStepMethod.linearize(
         solver, body_aero, base_inputs;
         theta_idxs=1:4, va_idxs=5:7, omega_idxs=8:10, delta_idxs=11:14,
@@ -54,14 +54,14 @@ end
     # differ at the noise floor (~sqrt(eps)) and that floor dominates
     # the small Δ used at scale=1.
     function evaluate_at!(input_vec)
-        perturbed_theta = input_vec[1:4]
-        perturbed_va    = input_vec[5:7]
-        perturbed_omega = input_vec[8:10]
-        perturbed_delta = input_vec[11:14]
+        perturbed_theta  = input_vec[1:4]
+        perturbed_va_vec = input_vec[5:7]
+        perturbed_omega  = input_vec[8:10]
+        perturbed_delta  = input_vec[11:14]
         VortexStepMethod.unrefined_deform!(
             ram_wing, perturbed_theta, perturbed_delta; smooth=false)
         reinit!(body_aero; init_aero=false,
-            va=perturbed_va, omega=perturbed_omega)
+            va=perturbed_va_vec, omega=perturbed_omega)
         VortexStepMethod.solve!(solver, body_aero; log=false)
         return [solver.sol.force; solver.sol.moment;
                 solver.sol.moment_unrefined_dist]
diff --git a/test/filament/test_semi_infinite_filament.jl b/test/filament/test_semi_infinite_filament.jl
index 43b2149f..826ae0f9 100644
--- a/test/filament/test_semi_infinite_filament.jl
+++ b/test/filament/test_semi_infinite_filament.jl
@@ -7,13 +7,13 @@ function create_test_filament2()
     x1 = [0.0, 0.0, 0.0]
     direction = [1.0, 0.0, 0.0]
     filament_direction = 1
-    vel_mag = 1.0
+    va = 1.0
     filament = SemiInfiniteFilament{Float64}()
-    reinit!(filament, x1, direction, vel_mag, filament_direction)
+    reinit!(filament, x1, direction, va, filament_direction)
     return filament
 end
 
-function analytical_solution(control_point, gamma, x1, direction, filament_direction, vel_mag)
+function analytical_solution(control_point, gamma, x1, direction, filament_direction, va)
     gamma = -gamma  # Sign convention difference
     r1 = control_point - x1
     r1_cross_direction = cross(r1, direction)
@@ -21,7 +21,7 @@ function analytical_solution(control_point, gamma, x1, direction, filament_direc
     
     alpha0 = 1.25643
     nu = 1.48e-5
-    epsilon = sqrt(4 * alpha0 * nu * norm(r_perp) / vel_mag)
+    epsilon = sqrt(4 * alpha0 * nu * norm(r_perp) / va)
 
     if norm(r1_cross_direction) > epsilon
         K = (gamma / (4π * norm(r1_cross_direction)^2)) * 
@@ -162,16 +162,16 @@ end
 
     @testset "Constant azimuthal direction inside core" begin
         filament = create_test_filament2()
-        v_a = filament.vel_mag
+        va = filament.vel_mag
 
         d_inside = 1e-4
         v1 = zeros(3); v2 = zeros(3)
         velocity_3D_trailing_vortex_semiinfinite!(
             v1, filament, filament.direction,
-            [0.5, d_inside, 0.0], gamma, v_a, work_vectors)
+            [0.5, d_inside, 0.0], gamma, va, work_vectors)
         velocity_3D_trailing_vortex_semiinfinite!(
             v2, filament, filament.direction,
-            [0.5, 2 * d_inside, 0.0], gamma, v_a, work_vectors)
+            [0.5, 2 * d_inside, 0.0], gamma, va, work_vectors)
 
         @test isapprox(normalize(v2), normalize(v1); atol=1e-8)
     end
diff --git a/test/plotting/test_plotting.jl b/test/plotting/test_plotting.jl
index 5a0476b7..d3ffc9cd 100644
--- a/test/plotting/test_plotting.jl
+++ b/test/plotting/test_plotting.jl
@@ -32,7 +32,7 @@ function create_body_aero()
     n_panels = 20          # Number of panels
     span = 20.0            # Wing span [m]
     chord = 1.0            # Chord length [m]
-    v_a = 20.0             # Magnitude of inflow velocity [m/s]
+    va = 20.0              # Magnitude of inflow velocity [m/s]
     alpha_deg = 30.0       # Angle of attack [degrees]
     alpha = deg2rad(alpha_deg)
 
@@ -49,8 +49,8 @@ function create_body_aero()
 
     refine!(wing)
     body_aero = BodyAerodynamics([wing])
-    vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
-    set_va!(body_aero, vel_app)
+    va_vec = [cos(alpha), 0.0, sin(alpha)] .* va
+    set_va!(body_aero, va_vec)
     body_aero
 end
 
@@ -98,7 +98,7 @@ end
     @test fig isa Figure
 
     # Plot polar curves
-    v_a = 20.0
+    va = 20.0
     angle_range = range(0, 20, 20)
     fig = plot_polars(
         [llt_solver, vsm_solver],
@@ -106,7 +106,7 @@ end
         ["VSM", "LLT"],
         angle_range=angle_range,
         angle_type="angle_of_attack",
-        v_a=v_a,
+        v_a=va,
         title="Rectangular Wing Polars",
         data_type=".png",
         save_path=save_dir,
@@ -123,7 +123,7 @@ end
         ["VSM", "LLT"],
         angle_range=angle_range,
         angle_type="angle_of_attack",
-        v_a=v_a,
+        v_a=va,
         title="Polars CL vs CD",
         is_save=false,
         is_show=false,
@@ -137,7 +137,7 @@ end
         angle_range=angle_range,
         angle_type="angle_of_attack",
         angle_of_attack=30.0,
-        v_a=v_a,
+        v_a=va,
         title="Combined Analysis",
         is_save=false,
         is_show=false,
@@ -151,7 +151,7 @@ end
         angle_range=angle_range,
         angle_type="angle_of_attack",
         angle_of_attack=30.0,
-        v_a=v_a,
+        v_a=va,
         title="Combined CL vs CD",
         is_save=false,
         is_show=false,
@@ -182,8 +182,8 @@ end
 
     body_aero_distributed = create_body_aero()
     n_panels = length(body_aero_distributed.panels)
-    va_distribution = repeat([12.0 0.0 1.0], n_panels, 1)
-    set_va!(body_aero_distributed, va_distribution)
+    va_vec_dist = repeat([12.0 0.0 1.0], n_panels, 1)
+    set_va!(body_aero_distributed, va_vec_dist)
 
     @test body_aero_distributed.has_distributed_va
     fig = plot_geometry(
diff --git a/test/solver/test_flow_curvature.jl b/test/solver/test_flow_curvature.jl
index 645a7d12..7c895ca9 100644
--- a/test/solver/test_flow_curvature.jl
+++ b/test/solver/test_flow_curvature.jl
@@ -134,16 +134,16 @@ end
 
     @testset "distributed rates drive a deformation mode" begin
         n = length(body_aero.panels)
-        va_dist = repeat([V 0.0 0.0], n)
+        va_vec_dist = repeat([V 0.0 0.0], n)
 
-        set_va!(body_aero, va_dist)
+        set_va!(body_aero, va_vec_dist)
         @test all(iszero, body_aero.pitch_rate_dist)
         solve!(solver_on, body_aero)
         base = copy(solver_on.sol.cm_dist)
 
         # antisymmetric twist rate: no rigid-body omega can express this
         rates = [panel.aero_center[2] > 0 ? 1.0 : -1.0 for panel in body_aero.panels]
-        set_va!(body_aero, va_dist; pitch_rate_dist=rates)
+        set_va!(body_aero, va_vec_dist; pitch_rate_dist=rates)
         @test body_aero.pitch_rate_dist ≈ rates
         solve!(solver_on, body_aero)
 
@@ -156,7 +156,7 @@ end
         @test sign(solver_on.sol.cm_dist[1] - base[1]) ==
               -sign(solver_on.sol.cm_dist[n] - base[n])
 
-        @test_throws ArgumentError set_va!(body_aero, va_dist;
+        @test_throws ArgumentError set_va!(body_aero, va_vec_dist;
                                            pitch_rate_dist=rates[1:end-1])
     end
 end
diff --git a/test/solver/test_forwarddiff.jl b/test/solver/test_forwarddiff.jl
index c052aae4..fae465f6 100644
--- a/test/solver/test_forwarddiff.jl
+++ b/test/solver/test_forwarddiff.jl
@@ -17,9 +17,9 @@ relative_error(jac, reference) = maximum(abs.(jac .- reference)) / maximum(abs,
     refine!(wing)
     body_aero = BodyAerodynamics([wing])
 
-    va = [15.0, 1.0, 2.0]
+    va_vec = [15.0, 1.0, 2.0]
     omega = [0.0, 0.0, 0.0]
-    y0 = [va; omega]
+    y0 = [va_vec; omega]
 
     @testset "AutoForwardDiff matches AutoFiniteDiff (LOOP, INVISCID)" begin
         solver = Solver(body_aero;
@@ -65,10 +65,10 @@ relative_error(jac, reference) = maximum(abs.(jac .- reference)) / maximum(abs,
             use_gamma_prev=false,
         )
 
-        v_a = 15.0
+        va = 15.0
         aoa_rad = deg2rad(7.5)
         y_op = [zeros(4);
-                [cos(aoa_rad), 0.0, sin(aoa_rad)] * v_a;
+                [cos(aoa_rad), 0.0, sin(aoa_rad)] * va;
                 zeros(3)]
 
         jac_fwd, _, conv_fwd = VortexStepMethod.linearize(
diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl
index 56896de1..59e37487 100644
--- a/test/solver/test_solver.jl
+++ b/test/solver/test_solver.jl
@@ -25,8 +25,8 @@ end
             @test solver.density == 1.225
 
             # Test that the solver can solve
-            va = [10.0, 0.0, 0.0]
-            set_va!(body_aero, va)
+            va_vec = [10.0, 0.0, 0.0]
+            set_va!(body_aero, va_vec)
             sol = solve!(solver, body_aero)
             @test sol isa VSMSolution
 
@@ -81,16 +81,16 @@ end
         wing = Wing(settings)
         refine!(wing)
         body_aero = BodyAerodynamics([wing])
-        va = [10.0, 0.0, 5.0]   # 26.6 deg angle of attack, past stall
+        va_vec = [10.0, 0.0, 5.0]   # 26.6 deg angle of attack, past stall
         nonlin = Solver(body_aero; solver_type=NONLIN, aerodynamic_model_type=VSM,
             type_initial_gamma_distribution=ELLIPTIC)
         loop = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM,
             type_initial_gamma_distribution=ELLIPTIC)
 
-        set_va!(body_aero, va)
+        set_va!(body_aero, va_vec)
         sol_nonlin = solve!(nonlin, body_aero)
         gamma_nonlin = copy(sol_nonlin.gamma_distribution)
-        set_va!(body_aero, va)
+        set_va!(body_aero, va_vec)
         sol_loop = solve!(loop, body_aero)
 
         @test sol_nonlin.solver_status == FEASIBLE
@@ -128,8 +128,9 @@ end
         solver = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM,
             type_initial_gamma_distribution=ELLIPTIC)
 
-        for va in ([10.0, 0.0, 0.0], [10.0, 0.0, 5.0])   # 0 deg, and 26.6 deg past stall
-            set_va!(body_aero, va)
+        # 0 deg, and 26.6 deg past stall
+        for va_vec in ([10.0, 0.0, 0.0], [10.0, 0.0, 5.0])
+            set_va!(body_aero, va_vec)
             gamma = copy(solve!(solver, body_aero).gamma_distribution)
             @test solver.lr.converged
             residual = maximum(abs, unrelaxed_step(body_aero, gamma) .- gamma)
diff --git a/test/solver/test_unrefined_dist.jl b/test/solver/test_unrefined_dist.jl
index 67f0bee3..08e955e7 100644
--- a/test/solver/test_unrefined_dist.jl
+++ b/test/solver/test_unrefined_dist.jl
@@ -24,8 +24,8 @@ using Test
             solver = Solver(body_aero, settings)
 
             # Set conditions and solve
-            va = [10.0, 0.0, 0.0]
-            set_va!(body_aero, va)
+            va_vec = [10.0, 0.0, 0.0]
+            set_va!(body_aero, va_vec)
             sol = solve!(solver, body_aero)
 
             # Test 1: Unrefined arrays exist and have correct size
@@ -98,8 +98,8 @@ using Test
             body_aero = BodyAerodynamics([wing])
             solver = Solver(body_aero, settings)
 
-            va = [10.0, 0.0, 0.0]
-            set_va!(body_aero, va)
+            va_vec = [10.0, 0.0, 0.0]
+            set_va!(body_aero, va_vec)
             sol = solve!(solver, body_aero)
 
             panels = body_aero.panels
@@ -171,8 +171,8 @@ using Test
                 body_aero = BodyAerodynamics([wing])
                 solver = Solver(body_aero, settings)
 
-                va = [10.0, 0.0, 0.0]
-                set_va!(body_aero, va)
+                va_vec = [10.0, 0.0, 0.0]
+                set_va!(body_aero, va_vec)
                 sol = solve!(solver, body_aero)
 
                 # Verify arrays have correct size
diff --git a/test/solver/test_viscous_drag_correction.jl b/test/solver/test_viscous_drag_correction.jl
index 2b6a254a..8e8bf8e8 100644
--- a/test/solver/test_viscous_drag_correction.jl
+++ b/test/solver/test_viscous_drag_correction.jl
@@ -32,19 +32,20 @@ end
     solver_off = Solver(body_aero; use_gamma_prev=false)
     solver_on = Solver(body_aero; use_gamma_prev=false,
                        is_with_viscous_drag_correction=true)
-    va_sideslip = V .* [cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta)]
-    va_straight = V .* [cos(alpha), 0.0, sin(alpha)]
+    va_vec_sideslip = V .* [cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta)]
+    va_vec_straight = V .* [cos(alpha), 0.0, sin(alpha)]
 
-    function force_dist_at(solver, va)
-        set_va!(body_aero, va)
+    function force_dist_at(solver, va_vec)
+        set_va!(body_aero, va_vec)
         solve!(solver, body_aero)
         return copy(solver.sol.f_body_3D)
     end
 
     @testset "adds the model's drag and spanwise force to each panel" begin
         density, mu = solver_on.density, solver_on.mu
-        for va in (va_sideslip, va_straight)
-            delta_force = force_dist_at(solver_on, va) .- force_dist_at(solver_off, va)
+        for va_vec in (va_vec_sideslip, va_vec_straight)
+            delta_force = force_dist_at(solver_on, va_vec) .-
+                          force_dist_at(solver_off, va_vec)
             for (i, panel) in enumerate(body_aero.panels)
                 v_normal = solver_on.lr.v_a_dist[i]
                 v_span = solver_on.lr.v_span_dist[i]
@@ -61,16 +62,17 @@ end
     end
 
     @testset "sideslip drives the spanwise flow and raises the drag" begin
-        delta_force = force_dist_at(solver_on, va_sideslip) .-
-                      force_dist_at(solver_off, va_sideslip)
+        delta_force = force_dist_at(solver_on, va_vec_sideslip) .-
+                      force_dist_at(solver_off, va_vec_sideslip)
         for (i, panel) in enumerate(body_aero.panels)
-            @test solver_on.lr.v_span_dist[i] ≈ dot(va_sideslip, panel.y_airf) rtol = 0.05
-            @test dot(delta_force[:, i], va_sideslip) > 0
+            @test solver_on.lr.v_span_dist[i] ≈
+                  dot(va_vec_sideslip, panel.y_airf) rtol = 0.05
+            @test dot(delta_force[:, i], va_vec_sideslip) > 0
         end
     end
 
     @testset "solve reports the corrected forces" begin
-        set_va!(body_aero, va_sideslip)
+        set_va!(body_aero, va_vec_sideslip)
         solve!(solver_on, body_aero)
         results = solve(solver_on, body_aero)
         @test [results["Fx"], results["Fy"], results["Fz"]] ≈ solver_on.sol.force
@@ -78,16 +80,16 @@ end
     end
 
     @testset "linearize reports the corrected forces" begin
-        y = [va_sideslip; zeros(3)]
+        y = [va_vec_sideslip; zeros(3)]
         results_for(solver) = VortexStepMethod.linearize(solver, body_aero, y;
             theta_idxs=nothing, va_idxs=1:3, omega_idxs=4:6)[2]
         results_on, results_off = results_for(solver_on), results_for(solver_off)
-        @test results_on[1:3] ≈ vec(sum(force_dist_at(solver_on, va_sideslip); dims=2))
+        @test results_on[1:3] ≈ vec(sum(force_dist_at(solver_on, va_vec_sideslip); dims=2))
         @test !(results_on[1:3] ≈ results_off[1:3])
     end
 
     @testset "calc_forces! stays zero-alloc" begin
-        set_va!(body_aero, va_sideslip)
+        set_va!(body_aero, va_vec_sideslip)
         solve!(solver_on, body_aero)
         calc_forces!(solver_on, body_aero)
         @test (@allocated calc_forces!(solver_on, body_aero)) == 0
diff --git a/test/wake/test_wake.jl b/test/wake/test_wake.jl
index f46feaf9..4c21405d 100644
--- a/test/wake/test_wake.jl
+++ b/test/wake/test_wake.jl
@@ -11,16 +11,16 @@ using VortexStepMethod
 
         # Test that frozen_wake! doesn't throw errors
         n_panels = length(body_aero.panels)
-        va_distribution = ones(n_panels, 3)  # Create velocity distribution for all panels
-        va_distribution[:, 1] .= 15.0  # X velocity = 15.0 m/s
-        va_distribution[:, 2] .= 0.0   # Y velocity = 0.0
-        va_distribution[:, 3] .= 0.0   # Z velocity = 0.0
+        va_vec_dist = ones(n_panels, 3)  # Create velocity distribution for all panels
+        va_vec_dist[:, 1] .= 15.0  # X velocity = 15.0 m/s
+        va_vec_dist[:, 2] .= 0.0   # Y velocity = 0.0
+        va_vec_dist[:, 3] .= 0.0   # Z velocity = 0.0
 
         # This should not throw an error
-        @test_nowarn VortexStepMethod.frozen_wake!(body_aero, va_distribution)
+        @test_nowarn VortexStepMethod.frozen_wake!(body_aero, va_vec_dist)
 
         # Test that the function accepts the right number of velocity vectors
-        @test size(va_distribution, 1) == length(body_aero.panels)
-        @test size(va_distribution, 2) == 3  # 3D velocity vectors
+        @test size(va_vec_dist, 1) == length(body_aero.panels)
+        @test size(va_vec_dist, 2) == 3  # 3D velocity vectors
     end
 end