Merge origin/main into agent/330-add-rigid-body-stability-derivatives-and

set_va!(body_aero, settings) keeps both sides: the inflow comes from
apparent_wind, and the body turns at settings.condition.yaw_rate about
body_aero.reference_point, as #353 made it. Both Added entries in the
changelog are kept.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Author: 1-Bort-1

CHANGELOG.md

@@ -7,10 +7,19 @@
 - `stability_derivatives` gives the force and moment coefficients and their derivatives
   with respect to angle of attack and sideslip, and `trim_angle` the angles of attack at
   which `CMy` changes sign, with the slope that says whether each trim is stable.
+- `set_va!(body_aero, va_vec, omega; reference_point)` turns the body about
+  `reference_point` [m] instead of the origin. The point is stored on
+  `BodyAerodynamics`, starts at the origin, and is kept by later `set_va!`, `reinit!`
+  and `linearize` calls until it is given again.
 - Spanwise-flow viscous drag correction (Gaunaa et al. 2024,
   doi:10.1088/1742-6596/2767/2/022068): each section gets a drag increment and a force
   along its span from the flow across it, in `solve!`, `solve` and `linearize`. Opt-in
   via `is_with_viscous_drag_correction` (default `false`) on the solver settings.
+- `plot_section_polars(body_aero; panels, alphas, delta)` draws cl, cd and cm against α
+  per panel through `calculate_cl`/`calculate_cd`/`calculate_cm`, for every aero model
+  and at flap deflection `delta`, in one figure instead of one coefficient per call.
+- `linearize` takes a `BodyAerodynamics` with more than one wing; `theta_idxs` and
+  `delta_idxs` then run over the unrefined sections of all wings in order.
 
 ### Changed
 
@@ -21,11 +30,13 @@
 
 ### Fixed
 
+- `set_va!(body_aero, settings)` applies `condition.yaw_rate` as a turn rate about the
+  body z axis; it was read from the settings file and ignored.
 - The `VSMSolution` docstring gives `lift_dist`, `drag_dist` and `panel_moment_dist` in
   the per-unit-span units they hold, [N/m] and [Nm/m], instead of [N] and [Nm].
-
-### Fixed
-
+- Inside its vortex core, `velocity_3D_trailing_vortex!` induces an azimuthal velocity
+  instead of a radial one. Only points within the millimetre-scale Oseen core of a
+  panel's chordwise trailing segment were affected.
 - With `artificial_damping` on, an iteration whose circulation is already smooth no longer
   re-applies the previous iteration's damping correction.
 

docs/src/private_functions.md

@@ -67,6 +67,7 @@ panel_loads
 ```@docs
 velocity_3D_bound_vortex!
 velocity_3D_trailing_vortex!
+velocity_3D_vortex_segment!
 velocity_3D_trailing_vortex_semiinfinite!
 calculate_velocity_induced_bound_2D!
 calculate_velocity_induced_single_ring_semiinfinite!
@@ -87,6 +88,7 @@ calculate_filaments_for_plotting
 ### Mesh refinement and billowing
 ```@docs
 unrefined_deform!
+unrefined_section_range
 deform!
 compute_refined_panel_mapping!
 compute_refined_section_interpolation!
@@ -278,6 +280,7 @@ panel_contour
 panel_normal
 plate_hinge_local
 panel_plate_geometry
+panel_polar_curves
 PLATE_FACES
 Makie.plot!(ax, panel::VortexStepMethod.Panel)
 Makie.plot!(ax, body::VortexStepMethod.BodyAerodynamics)

docs/src/settings.md

@@ -32,7 +32,7 @@ condition:
   wind_speed: 10.0                # free-stream velocity magnitude [m/s]
   alpha: 5.0                      # angle of attack [°]
   beta: 0.0                       # sideslip angle [°]
-  yaw_rate: 0.0                   # yaw rate [°/s]
+  yaw_rate: 0.0                   # turn rate about the body z axis [°/s]
 
 wings:
   - name: main_wing               # label the wing carries into plots and output

examples/obj_to_yaml_kite.jl

@@ -80,10 +80,9 @@ if PLOT
     plot_geometry(body_aero, "Ram air kite (converted from .obj)"; is_show=true,
         view_elevation=15, view_azimuth=-120, use_tex=USE_TEX)
 
-    # Airfoils and per-section polars recovered from the converted geometry
+    # Airfoils and panel polars recovered from the converted geometry
     plot_airfoils(geometry_file; symmetric=true, is_show=true)
-    plot_section_polars(body_aero, :cl; is_show=true)
-    plot_section_polars(body_aero, :cd; is_show=true)
+    plot_section_polars(body_aero; panels=[1, 10], is_show=true)
 
     plot_polars([solver], [body_aero], ["VSM (NeuralFoil polars from .obj)"];
         angle_range=range(-5, 20, length=26), v_a=va,

ext/VortexStepMethodMakieExt.jl

@@ -1,7 +1,8 @@
 module VortexStepMethodMakieExt
 using MakieControlPlots.Makie, VortexStepMethod, LinearAlgebra, Statistics, DelimitedFiles
 import MakieControlPlots
-import VortexStepMethod: calculate_filaments_for_plotting
+import VortexStepMethod: calculate_filaments_for_plotting, calculate_cl, calculate_cd,
+    calculate_cm
 import VortexStepMethod: ObjAdapter, AirfoilAero
 
 export plot_geometry, plot_distribution, plot_polars, save_plot, show_plot,
@@ -1324,13 +1325,12 @@ function VortexStepMethod.plot_combined_analysis(
             xlabel="α [°]",
             ylabel="Cm")
 
-        cl_vals = [first_body.panels[1].cl_interp(a) for a in alphas]
-        cd_vals = [first_body.panels[1].cd_interp(a) for a in alphas]
-        cm_vals = [first_body.panels[1].cm_interp(a) for a in alphas]
+        panel = first_body.panels[1]
+        cl, cd, cm = panel_polar_curves([panel], alphas, panel.delta)
 
-        lines!(ax_cl_curve, alphas_deg, cl_vals; color=:blue, linewidth=2)
-        lines!(ax_cd_curve, alphas_deg, cd_vals; color=:red, linewidth=2)
-        lines!(ax_cm_curve, alphas_deg, cm_vals; color=:green, linewidth=2)
+        lines!(ax_cl_curve, alphas_deg, only(cl); color=:blue, linewidth=2)
+        lines!(ax_cd_curve, alphas_deg, only(cd); color=:red, linewidth=2)
+        lines!(ax_cm_curve, alphas_deg, only(cm); color=:green, linewidth=2)
     end
 
     # [2,1] Spanwise Distributions (3×3 grid)
@@ -1477,48 +1477,43 @@ function VortexStepMethod.plot_combined_analysis(
 end
 
 """
-    plot_section_polars(body_aero, coefficient=:cl; is_show=true,
-                        is_save=false, save_path=nothing, data_type=".png")
+    panel_polar_curves(panels, alphas, deltas) -> (cl, cd, cm)
+
+Lift, drag and moment coefficients of each of `panels` over `alphas` [rad], each panel
+at its flap deflection in `deltas` [rad] (or one shared deflection), as one vector per
+panel per coefficient.
+"""
+function panel_polar_curves(panels, alphas, deltas)
+    cl = collect.(eachrow(calculate_cl.(panels, alphas', deltas)))
+    cd = collect.(eachrow(calculate_cd.(panels, alphas', deltas)))
+    cm = collect.(eachrow(calculate_cm.(panels, alphas', deltas)))
+    return cl, cd, cm
+end
+
+"""
+    plot_section_polars(body_aero; kwargs...)
 
 Implementation of [`plot_section_polars`](@ref); rendered through `MakieControlPlots`.
 """
-function VortexStepMethod.plot_section_polars(body_aero::BodyAerodynamics,
-    coefficient::Symbol=:cl; is_show::Bool=true, is_save::Bool=false,
-    save_path=nothing, data_type::String=".png")
-
-    coefficient in (:cl, :cd, :cm) ||
-        throw(ArgumentError("coefficient must be :cl, :cd, or :cm, got :$coefficient"))
-    idx = coefficient === :cl ? 2 : coefficient === :cd ? 3 : 4
-    label = uppercasefirst(string(coefficient))
-
-    alphas_deg = nothing
-    series = Vector{Float64}[]
-    labels = String[]
-    for wing in body_aero.wings
-        for (s, section) in enumerate(wing.unrefined_sections)
-            section.aero_model == POLAR_VECTORS || continue
-            aero = section.aero_data
-            aero === nothing && continue
-            section_alphas = rad2deg.(aero[1])
-            if isnothing(alphas_deg)
-                alphas_deg = collect(section_alphas)
-            elseif length(section_alphas) != length(alphas_deg)
-                @warn "section $s has a different α grid; plotting against the first section's α"
-            end
-            push!(series, Float64.(aero[idx]))
-            push!(labels, "section $s")
-        end
-    end
-    isempty(series) && error("No POLAR_VECTORS sections found in body")
+function VortexStepMethod.plot_section_polars(body_aero::BodyAerodynamics;
+    panels=eachindex(body_aero.panels), alphas=deg2rad.(-20:0.5:30), delta=nothing,
+    is_show::Bool=true, is_save::Bool=false, save_path=nothing,
+    data_type::String=".png")
+
+    panel_indices = vcat(panels)
+    chosen_panels = body_aero.panels[panel_indices]
+    deltas = something.(delta, getproperty.(chosen_panels, :delta))
+    cl, cd, cm = panel_polar_curves(chosen_panels, alphas, deltas)
+    labels = ["panel $i ($(panel.aero_model))"
+              for (i, panel) in zip(panel_indices, chosen_panels)]
 
-    plt = MakieControlPlots.plot(alphas_deg, series;
-        xlabel="α [deg]", ylabel=label, title="$label per section",
-        labels=labels, disp=(is_show || is_save))
+    plt = MakieControlPlots.plotx(rad2deg.(alphas), cl, cd, cm;
+        xlabel="α [deg]", ylabels=["cl", "cd", "cm"], title="Section polars",
+        labels=[labels], disp=(is_show || is_save))
 
     if is_save && !isnothing(save_path)
         isdir(save_path) || mkpath(save_path)
-        MakieControlPlots.savefig(
-            joinpath(save_path, "section_polars_$(coefficient)$(data_type)"))
+        MakieControlPlots.savefig(joinpath(save_path, "section_polars$(data_type)"))
     end
     return plt
 end

src/VortexStepMethod.jl

@@ -175,17 +175,22 @@ in sequence.
 function plot_combined_analysis end
 
 """
-    plot_section_polars(body_aero::BodyAerodynamics, coefficient=:cl; kwargs...)
+    plot_section_polars(body_aero::BodyAerodynamics; panels=eachindex(body_aero.panels),
+                        alphas=deg2rad.(-20:0.5:30), delta=nothing, kwargs...)
 
-Plot one polar coefficient (`:cl`, `:cd`, or `:cm`) against angle of attack for
-every section of a wing using stored `POLAR_VECTORS` data. Rendered through
-`MakieControlPlots`.
+Plot the lift, drag and moment coefficients against angle of attack for the chosen
+`panels`, one curve per panel, as each panel's [`calculate_cl`](@ref),
+[`calculate_cd`](@ref) and [`calculate_cm`](@ref) evaluate them for its aero model.
+Rendered through `MakieControlPlots`; returns its plot object.
 
 # Arguments
 - `body_aero`: the [`BodyAerodynamics`](@ref) to plot
-- `coefficient`: `:cl`, `:cd`, or `:cm` (default: `:cl`)
 
 # Keyword arguments
+- `panels`: index or indices into `body_aero.panels` (default: all panels)
+- `alphas`: angles of attack [rad] (default: `deg2rad.(-20:0.5:30)`)
+- `delta`: flap deflection [rad] a `POLAR_MATRICES` panel is evaluated at
+  (default: `nothing`, each panel's own `delta`)
 - `is_show`: whether to display (default: `true`)
 - `is_save`: whether to save (default: `false`)
 - `save_path`: directory to save the figure (default: `nothing`)

src/body_aerodynamics.jl

@@ -8,6 +8,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru
 - wings::Vector{W}: A vector of wings of type `W <: AbstractWing`; a body can have multiple wings
 - `va::MVec3` = zeros(MVec3):   A vector of the apparent wind speed, see: [`MVec3`](@ref)
 - `omega`::MVec3 = zeros(MVec3): A vector of the turn rates around the kite body axes
+- `reference_point`::MVec3 = zeros(MVec3): The point `omega` turns the body about [m]
 - `gamma_distribution`=zeros(Float64, P): A vector of the circulation
                         of the velocity field; Length: Number of segments. [m²/s]
 - `alpha_uncorrected`=zeros(Float64, P): angles of attack per panel
@@ -35,6 +36,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru
     _va::MVector{3, T} = zeros(MVector{3, T})
     has_distributed_va::Bool = false
     omega::MVector{3, T} = zeros(MVector{3, T})
+    reference_point::MVector{3, T} = zeros(MVector{3, T})
     gamma_distribution::MVector{P, T} = zeros(MVector{P, T})
     alpha_uncorrected::MVector{P, T} = zeros(MVector{P, T})
     alpha_corrected::MVector{P, T} = zeros(MVector{P, T})
@@ -155,6 +157,8 @@ function Base.setproperty!(obj::BodyAerodynamics, sym::Symbol, val)
         set_va!(obj, val)
     elseif sym === :omega
         set_va!(obj, obj._va, val)
+    elseif sym === :reference_point
+        set_va!(obj, obj._va, obj.omega; reference_point=val)
     else
         setfield!(obj, sym, val)
     end
@@ -252,6 +256,37 @@ function calculate_stall_angle_list!(stall_angles::AbstractVector,
     return nothing
 end
 
+"""
+    unrefined_section_range(body_aero::BodyAerodynamics, wing_idx)
+
+Indices of the unrefined sections of wing `wing_idx` in a distribution that runs over the
+unrefined sections of all wings in order, such as `moment_unrefined_dist`.
+"""
+function unrefined_section_range(body_aero::BodyAerodynamics, wing_idx)
+    offset = 0
+    for i in 1:wing_idx-1
+        offset += body_aero.wings[i].n_unrefined_sections
+    end
+    return offset .+ (1:body_aero.wings[wing_idx].n_unrefined_sections)
+end
+
+"""
+    unrefined_deform!(body_aero::BodyAerodynamics, theta_angles, delta_angles)
+
+Deform each wing of `body_aero` by its entries of `theta_angles` and `delta_angles` [rad],
+which run over the unrefined sections of all wings in order; `nothing` leaves that angle
+unchanged. Call [`reinit!`](@ref) afterwards to update the panels.
+"""
+function unrefined_deform!(body_aero::BodyAerodynamics, theta_angles, delta_angles)
+    for (wing_idx, wing) in enumerate(body_aero.wings)
+        section_range = unrefined_section_range(body_aero, wing_idx)
+        unrefined_deform!(wing,
+            isnothing(theta_angles) ? nothing : view(theta_angles, section_range),
+            isnothing(delta_angles) ? nothing : view(delta_angles, section_range))
+    end
+    return nothing
+end
+
 """
     reinit!(body_aero::BodyAerodynamics; init_aero, va, omega, refine_mesh, recompute_mapping, sort_sections)
 
@@ -1052,44 +1087,32 @@ end
 
 
 """
-    set_va!(body_aero::BodyAerodynamics, va_vec::VelVector, omega=zeros(MVec3))
+    set_va!(body_aero::BodyAerodynamics, va_vec::VelVector, omega=zeros(MVec3);
+            reference_point=body_aero.reference_point)
 
-Set velocity array and update wake filaments.
+Set a uniform apparent wind and a body turn rate, and update the wake filaments. Each
+panel sees `va_vec - omega × (control_point - reference_point)`.
 
 # Arguments
 - body_aero::BodyAerodynamics: The [`BodyAerodynamics`](@ref) struct to modify
 - `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]
+- `reference_point`: Point the body turns about, stored on `body_aero`  [m]
 
 `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_vec::AbstractVector,
-                 omega=zeros(MVector{3, T})) where {P, W, T}
-    n_panels = length(body_aero.panels)
-    va_vec_dist = zeros(T, n_panels, 3)
+                 omega=zeros(MVector{3, T});
+                 reference_point=body_aero.reference_point) where {P, W, T}
     body_aero.omega .= omega
+    body_aero.reference_point .= reference_point
     set_pitch_rate_dist!(body_aero, omega)
 
-    if all(iszero, omega)
-        va_vec_dist .= reshape(va_vec, 1, 3)
-    else
-        idx = 1
-        for wing in body_aero.wings
-            panel_end = idx + wing.n_panels - 1
-
-            # Calculate velocities for each panel in this wing slice
-            for j in idx:panel_end
-                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
-    end
-
-    # Update panel velocities
+    va_vec_dist = zeros(T, P, 3)
     for (i, panel) in enumerate(body_aero.panels)
-        panel.va .= va_vec_dist[i,:]
+        panel.va .= va_vec .- omega × (panel.control_point .- body_aero.reference_point)
+        va_vec_dist[i, :] .= panel.va
     end
 
     # Update wake elements
@@ -1147,7 +1170,8 @@ apparent_wind(alpha, beta, wind_speed) =
     set_va!(body_aero::BodyAerodynamics, settings::VSMSettings)
 
 Set the uniform inflow of `body_aero` to the [`apparent_wind`](@ref) at the `alpha` and
-`beta` [°] and `wind_speed` [m/s] of `settings.condition`.
+`beta` [°] and `wind_speed` [m/s] of `settings.condition`, turning the body about
+`body_aero.reference_point` at its `yaw_rate` [°/s] about Z_b.
 
 # Example
 ```julia
@@ -1160,5 +1184,5 @@ function set_va!(body_aero::BodyAerodynamics, settings::VSMSettings)
     condition = settings.condition
     va_vec = apparent_wind(deg2rad(condition.alpha), deg2rad(condition.beta),
         condition.wind_speed)
-    set_va!(body_aero, va_vec)
+    set_va!(body_aero, va_vec, [0.0, 0.0, deg2rad(condition.yaw_rate)])
 end