Write the panel aerodynamics once, as traceable functions

CHANGELOG.md

@@ -1,5 +1,31 @@
 # Changelog
 
+## Unreleased
+
+### Added
+- Shared panel aerodynamics (`src/panel_aerodynamics.jl`): the per-panel physics
+  written once as pure, branch-free, number-type-generic functions of the section
+  geometry and the flow — `panel_axes`, `panel_inflow`, `panel_force_directions`,
+  `panel_loads` and the small helpers around them. `update_panel_properties!`,
+  `init_pos!`, `calc_forces!` and `calculate_results` now call them instead of
+  spelling the algebra out three times, and `SymbolicAWEModels` traces the same
+  functions with symbolic arguments to build its equations, so the two packages
+  can no longer drift apart. Panel geometry is unchanged bit for bit; forces,
+  moments and coefficients agree to within 1 ulp, the products having been
+  reassociated. `calc_forces!` stays zero-allocation.
+- `effective_alpha` and the `deficiency` argument of `panel_inflow` carry an
+  unsteady lag (a Wagner indicial deficiency) into the angle the polars are read
+  at, while the geometric angle still turns the force. Unused by the solver,
+  which has no unsteady state; it is shared so a symbolic consumer that does have
+  one reads the same definition.
+
+### Changed
+- `section_pitch_rate` gained a three-argument form taking the trailing minus
+  leading edge apparent wind directly. The four-argument form is unchanged.
+- Norms inside the shared panel aerodynamics are floored (`smooth_norm`) rather
+  than guarded by branches, so the expressions are differentiable and traceable.
+  The floor is `1e-12` m, which moves no reported quantity.
+
 ## VortexStepMethod v4.1.2 2026-08-22
 
 ### Fixed

Project.toml

@@ -1,7 +1,7 @@
 name = "VortexStepMethod"
 uuid = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9"
 authors = ["1-Bart-1 <bart@vandelint.net>", "Oriol Cayon and contributors"]
-version = "4.1.2"
+version = "4.2.0"
 
 [workspace]
 projects = ["examples", "docs", "test"]

docs/src/private_functions.md

@@ -17,7 +17,6 @@ calculate_cl
 calculate_cd
 calculate_cm
 calculate_cd_cm
-flow_curvature_cm
 set_pitch_rate_dist!
 calculate_relative_alpha_and_velocity
 calculate_relative_alpha_and_relative_velocity
@@ -31,6 +30,27 @@ smooth_distribution!
 make_dual_shadow
 ```
 
+### Panel aerodynamics
+The per-panel aerodynamics, written once as pure, branch-free functions of the
+section geometry and the flow. `SymbolicAWEModels` traces the same functions with
+symbolic arguments to build its equations, so both packages evaluate one
+definition of the physics.
+```@docs
+SMOOTH_FLOOR
+smooth_norm
+panel_span_vector
+panel_chord_weight
+panel_axes
+effective_alpha
+panel_inflow
+dynamic_pressure
+flow_curvature_cm
+panel_force_directions
+panel_moment
+panel_couple_force
+panel_loads
+```
+
 ### Induced velocities
 ```@docs
 velocity_3D_bound_vortex!

src/VortexStepMethod.jl

@@ -414,6 +414,7 @@ function help(url)
 end
 
 # Include core functionality
+include("panel_aerodynamics.jl")
 include("settings.jl")
 include("section_aero.jl")
 include("wing_geometry.jl")

src/body_aerodynamics.jl

@@ -727,38 +727,6 @@ function compute_panel_center_of_pressures(
     return panel_cp_locations
 end
 
-"""
-    flow_curvature_cm(pitch_rate, chord, v_rel)
-
-Quarter-chord moment increment of a section rotating about its own spanwise axis,
-from thin airfoil theory. The rotation makes the local incidence vary linearly
-along the chord, which is equivalent to parabolic camber and yields
-`Δcm = -(π/4) q̂` with `q̂ = q c / (2 v_rel)` and `q` positive nose-up.
-Independent of the pivot location; the lift response to `q` needs no correction
-because the inflow is already sampled at the three-quarter-chord control point.
-"""
-@inline function flow_curvature_cm(pitch_rate, chord, v_rel)
-    v_rel > 0 || return zero(chord)
-    return -0.25π * pitch_rate * chord / (2v_rel)
-end
-
-"""
-    section_pitch_rate(velocity_leading, velocity_trailing, z_airf, chord)
-
-Rate at which a section rotates about its own spanwise axis, from the velocities
-of its leading and trailing edge. Positive nose-up, matching
-[`flow_curvature_cm`](@ref). Use this to build a `pitch_rate_dist` for
-[`set_va!`](@ref) from a deforming structure, where twist and flapping rates
-differ per section and no single body rate describes them.
-"""
-@inline function section_pitch_rate(velocity_leading, velocity_trailing,
-                                    z_airf, chord)
-    chord > 0 || return zero(chord)
-    normal_rate = dot3(velocity_trailing, z_airf) -
-                  dot3(velocity_leading, z_airf)
-    return -normal_rate / chord
-end
-
 """
     set_pitch_rate_dist!(body_aero, omega)
 
@@ -927,43 +895,26 @@ function calculate_results(
     cross3!(dir_side_ref, dir_lift_ref, va_ref_unit)
     q_ref = 0.5 * density * va_ref_mag^2
 
-    induced_va_airfoil = body_aero.work_vectors[4]
-    dir_lift_induced_va = body_aero.work_vectors[5]
-    dir_drag_induced_va = body_aero.work_vectors[6]
     lift_induced_va = body_aero.work_vectors[7]
     drag_induced_va = body_aero.work_vectors[8]
     dir_lift_prescribed_va = body_aero.work_vectors[9]
     temp_vec = body_aero.work_vectors[10]
+    spanwise_unit = SVector{3}(spanwise_direction)
 
     # Main calculation loop
     for (i, panel) in enumerate(panels)
         panel_area = panel.chord * panel.width
         area_all_panels += panel_area
 
-        alpha_corrected_i = alpha_corrected[i]
-        c_alpha = cos(alpha_corrected_i)
-        s_alpha = sin(alpha_corrected_i)
-        @inbounds for k in 1:3
-            induced_va_airfoil[k] = c_alpha * panel.x_airf[k] +
-                                    s_alpha * panel.z_airf[k]
-        end
-        normalize3!(induced_va_airfoil)
-
-        cross3!(dir_lift_induced_va,
-                induced_va_airfoil, panel.y_airf)
-        normalize3!(dir_lift_induced_va)
-        cross3!(dir_drag_induced_va,
-                spanwise_direction, dir_lift_induced_va)
-        normalize3!(dir_drag_induced_va)
-
-        q_lift = 0.5 * density * v_a_dist[i]^2
-        lift_i = cl_array[i] * q_lift * chord_array[i]
-        drag_i = cd_array[i] * q_lift * chord_array[i]
-        moment_i = cm_array[i] * q_lift * chord_array[i]^2
-
+        axes = panel_axes(panel)
+        dirs = panel_force_directions(axes, alpha_corrected[i], spanwise_unit)
+        loads = panel_loads(axes, dirs,
+            dynamic_pressure(density, density, v_a_dist[i]),
+            cl_array[i], cd_array[i], cm_array[i])
+        moment_i = loads.moment
         @inbounds for k in 1:3
-            lift_induced_va[k] = lift_i * dir_lift_induced_va[k]
-            drag_induced_va[k] = drag_i * dir_drag_induced_va[k]
+            lift_induced_va[k] = loads.lift * dirs.dir_lift[k]
+            drag_induced_va[k] = loads.drag * dirs.dir_drag[k]
         end
 
         va_panel_mag = va_norm_array[i]

src/panel.jl

@@ -100,16 +100,15 @@ function init_pos!(
     panel.LE_point_1 .= section_1.LE_point
     panel.TE_point_2 .= section_2.TE_point
     panel.LE_point_2 .= section_2.LE_point
-    panel.chord = (
-        norm(panel.TE_point_1 - panel.LE_point_1) +
-        norm(panel.TE_point_2 - panel.LE_point_2)
-    ) / 2
+    panel.chord = panel_chord(
+        SVector{3}(section_1.LE_point), SVector{3}(section_1.TE_point),
+        SVector{3}(section_2.LE_point), SVector{3}(section_2.TE_point))
     panel.corner_points[:, 1] = panel.LE_point_1
     panel.corner_points[:, 2] = panel.TE_point_1
     panel.corner_points[:, 3] = panel.TE_point_2
     panel.corner_points[:, 4] = panel.LE_point_2
     vec .= bound_point_2 .- bound_point_1
-    panel.width = norm(vec)
+    panel.width = smooth_norm(vec)
     reinit!(panel.filaments[1], bound_point_2, bound_point_1, vec)
     reinit!(panel.filaments[2], bound_point_1, panel.TE_point_1, vec)
     reinit!(panel.filaments[3], panel.TE_point_2, bound_point_2, vec)
@@ -268,6 +267,18 @@ function reinit!(
 end
 
 
+"""
+    panel_axes(panel::Panel)
+
+The panel's stored airfoil frame and size in the shape [`panel_axes`](@ref)
+returns, so a panel built by the geometry pass feeds the same functions as one built
+straight from section points.
+"""
+@inline panel_axes(panel::Panel) =
+    (x_airf = SVector{3}(panel.x_airf), y_airf = SVector{3}(panel.y_airf),
+     z_airf = SVector{3}(panel.z_airf), chord = panel.chord,
+     width = panel.width)
+
 """
     calculate_relative_alpha_and_relative_velocity(panel::Panel, induced_velocity::Vector{Float64})
 
@@ -286,14 +297,8 @@ function calculate_relative_alpha_and_relative_velocity(
     panel::Panel{T},
     induced_velocity::AbstractVector{T}
 ) where T
-    # Calculate relative velocity and angle of attack
-    # Constants throughout iterations: panel.va, panel.x_airf, panel.y_airf
-    relative_velocity = panel.va .+ induced_velocity
-    v_normal = dot(panel.z_airf, relative_velocity)
-    v_tangential = dot(panel.x_airf, relative_velocity)
-    alpha = atan(v_normal, v_tangential)
-
-    return alpha, relative_velocity
+    flow = panel_inflow(panel_axes(panel), panel.va, panel.va, induced_velocity)
+    return flow.alpha, flow.v_eff
 end
 
 """
@@ -302,11 +307,8 @@ end
 Calculate relative angle of attack and relative velocity of the panel.
 """
 function calculate_relative_alpha_and_velocity(panel::Panel, induced_velocity)
-    relative_velocity = panel.va + induced_velocity
-    v_normal = dot(panel.z_airf, relative_velocity)
-    v_tangential = dot(panel.x_airf, relative_velocity)
-    alpha = atan(v_normal, v_tangential)
-    return alpha, relative_velocity
+    flow = panel_inflow(panel_axes(panel), panel.va, panel.va, induced_velocity)
+    return flow.alpha, flow.v_eff
 end
 
 """