calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::Model, 
                     core_radius_fraction,
                     va_norm_array, 
                     va_unit_array)

Calculate Aerodynamic Influence Coefficient matrices.

See also: BodyAerodynamics, Model

Returns: nothing

gamma_loop!(solver::Solver, AIC_x::Matrix{Float64},
          AIC_y::Matrix{Float64}, AIC_z::Matrix{Float64},
          panels::AbstractVector{<:Panel}, relaxation_factor::Float64; log=true)

Main iteration loop for calculating circulation distribution.

When solver.is_with_artificial_viscosity is set, the LOOP solver replaces the explicit target F(gamma) with the implicit Li/Gaunaa solution (I - diag(mu) L) gamma = F(gamma) before relaxation, stabilizing post-stall distributions. The Python reference gates this on precomputed per-panel stall angles; here the polar is an interpolation object, so the expensive linear solve is gated on any(mu > 0) instead, which is behaviourally identical.

build_spanwise_laplacian!(laplacian, n_panels)

Fill the n_panels × n_panels matrix laplacian with the discrete spanwise Laplacian used by the Li/Gaunaa post-stall artificial-viscosity regularization. Interior rows use the three-point stencil gamma[i-1] - 2 gamma[i] + gamma[i+1]; the tip rows use the second-order closures of Li, Gaunaa, Pirrung & Lønbæk (TORQUE 2026, Eq. 15) that enforce gamma -> 0 at the wing tips. Panels are assumed ordered consecutively along the span and approximately uniformly spaced. Returns the matrix unchanged (all zeros) when n_panels < 3.

local_lift_slope!(slopes, panels, alpha_dist, delta=deg2rad(0.5))

Per-panel local lift-curve slope dCl/dalpha, evaluated from each panel's own 2-D polar at its current effective angle of attack by central differences. The slope turns negative in post-stall, which is what activates the artificial-viscosity regularization in gamma_loop!.

apply_artificial_viscosity!(gamma, panels, alpha_dist, laplacian, viscosity_matrix,
                            lift_slope, mu_array, gamma_target, planform_area, factor)

Apply one implicit Li/Gaunaa artificial-viscosity step to gamma in place and return true when it fired. The per-panel viscosity is mu_i = max(0, -factor * planform_area * Cl'_i / width_i^2), with the local lift slope Cl'_i from local_lift_slope! evaluated at alpha_dist. When any panel is post-stall (mu_i > 0), gamma is replaced by the solution of (I - diag(mu) L) gamma = gamma, with L the spanwise Laplacian in laplacian (see build_spanwise_laplacian!); otherwise gamma is left unchanged. The remaining arguments are preallocated work buffers reused across iterations.

frozen_wake(body_aero::BodyAerodynamics, va_distribution)

Update the filaments of the panels with frozen wake model. Uses one shared wake vector computed from area-weighted distributed inflow.

Replaces older filaments if present by checking length of filaments.

Arguments

• body_aero::BodyAerodynamics: see: BodyAerodynamics
• va_distribution::Matrix{Float64}: Array of velocity vectors at each panel

Returns

• nothing

calc_forces!(solver::Solver, body_aero::BodyAerodynamics;
             reference_point=solver.reference_point, moment_frac=0.1)

Assemble aerodynamic forces and moments from the circulation already converged by solve_base! and stored in solver. Split out of solve!.

calculate_cl(panel::Panel, alpha)
calculate_cl(panel::Panel, alpha, delta)

Calculate lift coefficient for given angle of attack alpha [rad]. The 3-arg form evaluates the (α, δ) polar (POLAR_MATRICES) at the passed flap deflection delta [rad] instead of the panel's stored delta; the 2-arg form forwards with panel.delta. Other aero models ignore delta.

Returns

• Float64: Lift coefficient (Cl)

calculate_cd(panel::Panel, alpha)
calculate_cd(panel::Panel, alpha, delta)

Calculate the drag coefficient for the given angle of attack. The 3-arg form evaluates the (α, δ) polar at the passed flap deflection delta; see calculate_cl.

calculate_cm(panel::Panel, alpha)
calculate_cm(panel::Panel, alpha, delta)

Calculate the pitching-moment coefficient for the given angle of attack. The 3-arg form evaluates the (α, δ) polar at the passed flap deflection delta; see calculate_cl.

calculate_cd_cm(panel::Panel, alpha)

Calculate drag and moment coefficients for the given angle of attack (calculate_cd and calculate_cm).

calculate_relative_alpha_and_velocity(panel::Panel, induced_velocity)

Calculate relative angle of attack and relative velocity of the panel.

calculate_relative_alpha_and_relative_velocity(panel::Panel, induced_velocity::Vector{Float64})

Calculate the relative angle of attack and relative velocity of the panel.

Arguments

• panel::Panel: The panel object
• induced_velocity::Vector{Float64}: Induced velocity at the control point

Returns

• Tuple{Float64,Vector{Float64}}: Tuple containing:
  • alpha: Relative angle of attack of the panel (in radians)
  • relative_velocity: Relative velocity vector of the panel

update_effective_angle_of_attack_if_VSM(body_aero::BodyAerodynamics, gamma,
                                      core_radius_fraction,
                                      z_airf_array,
                                      x_airf_array,
                                      va_array,
                                      va_norm_array,
                                      va_unit_array)

Update angle of attack at aerodynamic center for VSM method.

Returns: nothing

calculate_stall_angle_list(panels::Vector{<:Panel};
                         begin_aoa=9.0,
                         end_aoa=22.0,
                         step_aoa=1.0,
                         stall_angle_if_none_detected=50.0,
                         cl_initial=-10.0)

Calculate stall angles for each panel.

Returns: Vector{Float64}: Stall angles in radians

calculate_circulation_distribution_elliptical_wing(body_aero::BodyAerodynamics, gamma_0=1.0)

Calculate circulation distribution for an elliptical wing.

Returns: nothing

_compute_reference_velocity_from_distribution(va_input, n_panels, panel_areas=nothing)

Return a single reference velocity vector from uniform or distributed inflow. For distributed inflow, the speed is area-weighted RMS and the direction is the area-weighted mean direction.

smooth_circulation!(damp, circulation, 
                  smoothness_factor::Float64, 
                  damping_factor::Float64)

Smooth circulation distribution if needed.

Returns:

• Tuple of smoothed circulation and boolean indicating if smoothing was applied

smooth_distribution!(dist, window_size)

Apply moving average smoothing to a distribution in-place. Uses a centered window of size window_size (must be odd).

Arguments

• dist::Vector{Float64}: Distribution to smooth (modified in-place)
• window_size::Int: Size of smoothing window (must be odd)

make_dual_shadow(solver, body_aero, ::Type{TD}) where TD

Build a BodyAerodynamics/Solver pair whose mutating buffers carry element type TD (typically a ForwardDiff.Dual). Wing geometry, polar coefficients and interpolators are reused from the Float64 originals; circulation/AIC/scratch buffers are freshly allocated as TD-typed.

velocity_3D_bound_vortex(vel, filament::BoundFilament, XVP, 
                       gamma, core_radius_fraction, work_vectors)

Calculate induced velocity by a bound vortex filament at a point in space.

velocity_3D_trailing_vortex(vel, filament::BoundFilament, 
                          XVP, gamma, v_a, work_vectors)

Calculate induced velocity by a trailing vortex filament.

Arguments

• XVP: Control point coordinates
• gamma: Vortex strength
• v_a: Inflow velocity magnitude
• work_vectors: preallocated array of intermediate variables

Reference: Rick Damiani et al. "A vortex step method for nonlinear airfoil polar data as implemented in KiteAeroDyn".

velocity_3D_trailing_vortex_semiinfinite(filament::SemiInfiniteFilament, 
                                         Vf, XVP, GAMMA, v_a, work_vectors)

Calculate induced velocity by a semi-infinite trailing vortex filament.

calculate_velocity_induced_bound_2D!(U2D, panel::Panel, evaluation_point, work_vectors)

Calculate velocity induced by bound vortex filaments at the control point. Only needed for VSM, as LLT bound and filament align, thus no induced velocity.

Arguments

• U_2D: Resulting 2D velocity vector
• panel::Panel: Panel object
• evaluation_point: Point where induced velocity is evaluated
• work_vectors: Pre-allocated temporary variables

calculate_velocity_induced_single_ring_semiinfinite!(
    velind::MVec3,
    tempvel::MVec3,
    filaments,
    evaluation_point::MVec3,
    evaluation_point_on_bound::Bool,
    va_norm::Float64,
    va_unit::MVec3,
    gamma::Float64,
    core_radius_fraction::Float64,
    work_vectors::NTuple{10, MVec3}
)

Calculate the velocity induced by a vortex ring at a control point.

Arguments

• velind
• tempvel
• 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_unit::MVec3: Unit vector of apparent velocity
• gamma::Float64: Circulation strength
• core_radius_fraction::Float64: Vortex core radius as fraction of panel width
• work_vectors::NTuple{10, MVec3} Pre-allocated temporary variables

Returns

• nothing

cross3!(result::AbstractVector{T}, a::AbstractVector{T}, b::AbstractVector{T}) where T

Compute cross product of 3D vectors in-place.

update_panel_properties!(panel_props::PanelProperties, section_list::Vector{Section}, n_panels::Int)

Update geometric properties for each panel.

Arguments

• section_list::Vector{Section}: List of Sections
• n_panels::Int: Number of Panels

Returns:

nothing, updates the PanelProperties in-place

build_interps(section_1, section_2, remove_nan) -> (cl, cd, cm, section_aero)

Build the averaged aerodynamic interpolations for the panel between two sections. Returns (cl_interp, cd_interp, cm_interp, section_aero), each nothing for models that do not use it (INVISCID, POLY). cl/cm clamp (Flat) past the alpha range but extrapolate linearly (Line) over delta; cd extrapolates linearly in both. The concrete return types parameterise Panel — see panel_interp_types.

panel_interp_types(section, remove_nan) -> (CL, CD, CM, CP)

Field types for a Panel whose aero model matches section. Since a wing has one aero model, every panel shares these, so Vector{Panel{...}} is concretely typed.

reinit!(wing::AbstractWing)

Reinitialize wing panel properties based on current refined_sections geometry.

This function only updates panel properties (chord, area, etc.) from the existing refinedsections. It does NOT refine the mesh - call refineaerodynamic_mesh!(wing) first if needed.

Note

After deformation via unrefined_deform!() or deform!(), call reinit! to update panel properties while preserving the deformed geometry.

reinit!(panel, section_1, section_2, aero_center, control_point, bound_point_1,
        bound_point_2, x_airf, y_airf, z_airf, delta, vec, spanwise_direction; kwargs...)

Reinitialize a panel's geometry, horseshoe filaments and aerodynamic interpolations.

The panel is oriented so its y_airf (and the bound vortex bound_2 -> bound_1) points along +spanwise_direction, with z_airf pointing to the airfoil upper surface. This makes the aero independent of section ordering: a reversed order would otherwise flip the normal and make the panel look up its polar at a negated angle of attack.

rotated_te(le, te, y_hat, θ)

Compute the trailing-edge position after rotating the chord vector te - le around y_hat by θ radians (Rodrigues). Returns an SVector{3} (stack-allocated, zero heap allocs).

calculate_filaments_for_plotting(panel::Panel)

Calculate filaments for plotting with their positions and colors.

Returns

• Vector{Tuple{Vector{Float64}, Vector{Float64}, String}}: List of tuples containing:
  • First point (x1)
  • Second point (x2)
  • Color string

unrefined_deform!(wing::Wing, theta_angles=nothing, delta_angles=nothing)

Apply deformation angles defined per unrefined section.

Refined-section twist and TE-deflection values are computed by linear interpolation from the unrefined-section inputs using the precomputed refined_section_left_idx / refined_section_weight cache built at refinement time. Endpoint refined sections take the unrefined endpoint values exactly. The panel-level theta_dist / delta_dist arrays are then filled by averaging adjacent refined-section values, so downstream consumers (solver, body aerodynamics) see a per-panel value.

Arguments

• wing::Wing: Wing to deform (must have nondeformedsections, populated by refine! for manual/YAML wings or by OBJ refinement for OBJ-based wings).
• theta_angles::AbstractVector: Twist angles in radians, one per unrefined section. Pass nothing to leave twist unchanged.
• delta_angles::AbstractVector: TE deflection angles in radians, one per unrefined section. Pass nothing to leave deflection unchanged.

Keyword arguments

• smooth, smooth_window: accepted for backwards compatibility with callers of deform!, but ignored here — the linear interpolation between unrefined sections is already smooth, so no post-hoc smoothing is applied.

deform!(wing::Wing, theta_dist::AbstractVector, delta_dist::AbstractVector;
        smooth=false, smooth_window=nothing)

Deform wing by applying theta and delta distributions at the panel level.

Arguments

• wing::Wing: Wing to deform (must support deformation)
• theta_dist::AbstractVector: Twist angles for each panel (length = n_panels)
• delta_dist::AbstractVector: TE deflections for each panel (length = n_panels)
• smooth::Bool: Whether to apply smoothing (default: false)
• smooth_window::Union{Nothing, Int}: Smoothing window size (default: auto-calculated)

Effects

Updates wing.refinedsections with deformed geometry based on wing.nondeformed_sections

deform!(wing::Wing; smooth=false, smooth_window=nothing)

Apply stored thetadist and deltadist to deform the wing geometry. Converts panel angles (npanels) to section angles (npanels+1) by averaging adjacent panels.

Arguments

• wing::Wing: Wing to deform (must have nondeformedsections)
• smooth::Bool: Whether to apply smoothing to thetadist and deltadist (default: false)
• smooth_window::Union{Nothing, Int}: Smoothing window size (default: auto-calculated)

Effects

Updates wing.refinedsections based on wing.nondeformed_sections and stored distributions

compute_refined_panel_mapping!(wing::AbstractWing)

Compute the mapping from refined panels to unrefined sections by finding the closest unrefined section for each refined panel (based on section center distance). Maps each refined panel index to its corresponding unrefined section index (1 to nunrefinedsections). Works after refinement is complete.

compute_refined_section_interpolation!(wing::AbstractWing)

Compute per-refined-section linear-interpolation weights from the unrefined sections. For refined section i, the interpolated value is:

out[i] = weight[i] * unrefined[left_idx[i]] +
         (1 - weight[i]) * unrefined[left_idx[i] + 1]

Positions are quarter-chord arc-length along the unrefined and refined sections. The first refined section is pinned to left_idx == 1, weight == 1 (returns unrefined[1] exactly) and the last refined section to left_idx == n_unref - 1, weight == 0 (returns unrefined[end] exactly).

copy_sections(sections) -> Vector{Section}

Copy a vector of Sections into fresh objects with their own LE_point/TE_point storage; the read-only aero_data/section_aero tables are shared by reference.

copy_sections_to_refined!(wing; reuse_aero_data=false)

Copy unrefined sections to refined sections 1:1 (no interpolation). If refined_sections is empty, allocates via copy; otherwise reinitialises in-place. Warns if billowing was requested but there are no intermediate sections to billow.

_apply_refined_section_thetas!(wing, section_thetas)

Rotate each refined section's TE point around its LE point by the given per-section twist angle, using wing.non_deformed_sections as the reference geometry.

_panel_thetas_to_section_thetas!(section_thetas, panel_thetas)

Convert panel-level twist angles (length n) to refined-section angles (length n+1). Interior sections take the average of adjacent panels. Boundary sections use linear extrapolation from the two nearest panels, so a linear input produces a linear output across the full refined-section range.

_interpolate_unrefined_to_refined(wing, unrefined_values) -> Vector

Linearly interpolate unrefined_values (length n_unrefined_sections) to refined sections (length n_panels + 1) using the cached weights. Endpoints match exactly.

reinit!(section::Section, LE_point::PosVector, TE_point::PosVector, aero_model=nothing, aero_data=nothing)

Function to update a Section in place.

refine_mesh_for_linear_cosine_distribution!(wing, idx, dist,
    n_sections, sections; endpoints, reuse_aero_data)

Refine wing mesh using linear or cosine spacing. Reads LE/TE directly from a Vector{Section} (zero matrix allocations).

refine_mesh_by_splitting_provided_sections!(wing; reuse_aero_data,
    billowing_percentage)

Refine mesh by splitting provided sections into desired number of panels.

When billowing_percentage > 0, rotates chord vectors around the leading edge with a sinusoidal profile to simulate fabric billowing between ribs.

refine_mesh_with_billowing!(wing; reuse_aero_data)

Refine wing mesh using SPLIT_PROVIDED spacing with TE billowing.

Between each pair of unrefined (rib) sections, chord vectors are rotated around the leading edge to simulate fabric billowing. The rotation amplitude is found iteratively so that the TE arc length matches the wing's billowing_percentage.

Delegates to refine_mesh_by_splitting_provided_sections!.

apply_billowing_to_pair!(sections, start_si, end_si, y_hat,
                         span_len, le_ref, te_left, te_right,
                         percentage)

Apply billowing to refined sections between two ribs by rotating each section's chord around y_hat. Uses Newton iteration to find the rotation amplitude that matches the target TE arc-length percentage, then applies the rotations in-place.

The angle profile is angle_max * sin(π t) (zero at ribs, maximum at centre). All angles are in radians.

Non-allocating: modifies sections[start_si:end_si].TE_point in-place using scalar arithmetic only.

billowing_arc_length(sections, start_si, end_si, y_hat,
                     span_len, le_ref, te_left, te_right,
                     angle_max)

Compute the TE arc length that would result from rotating each section's chord around y_hat by angle_max * sin(π t) radians, where t is the normalised spanwise position within the rib pair.

Non-allocating (uses scalar and SVector arithmetic only).

flip_created_coord_in_pairs_if_needed!(coord::Matrix{Float64})

Ensure coordinates are ordered from positive to negative along y-axis.

update_non_deformed_sections!(wing::AbstractWing)

Create nondeformedsections to match refinedsections. This enables deformation support for all wings (YAML and OBJ). Should be called after refinedsections are populated. Once populated, nondeformedsections serves as the undeformed reference geometry.

calculate_new_aero_data(sections, section_index, left_weight, right_weight)

Interpolate aerodynamic input between two adjacent sections (zero-copy variant).

calculate_new_aero_data(aero_model, aero_data, section_index,
                        left_weight, right_weight)

Interpolate aerodynamic input between two sections.

load_polar_data(csv_file_path::String) -> Tuple{Union{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}}}, Symbol}

Load aerodynamic polar data from a CSV file using only readlines.

The CSV file must contain a header row with columns for alpha, cl, cd, and cm (case-insensitive, order arbitrary). Each subsequent row should contain numeric values for these columns.

Arguments

• csv_file_path::String: Path to the CSV file containing polar data.

Returns

• A tuple (aero_data, model_type) where:
  • aero_data: A tuple of vectors (alpha, cl, cd, cm) if the file is valid, or nothing if invalid or missing.
    • alpha: Angle of attack in degrees (converted to radians internally).
    • cl: Lift coefficient.
    • cd: Drag coefficient.
    • cm: Moment coefficient.
  • model_type: POLAR_VECTORS if data is loaded, or INVISCID if not.

Behavior

• If the file is missing, empty, or invalid, a warning is issued and (nothing, INVISCID) is returned.

Example

# Create a YAML-based wing from configuration file
wing = Wing(
    "path/to/wing_config.yaml";
    n_panels=40,
    n_groups=4
)

load_matrix_polar_data(cl_path, cd_path, cm_path) -> (aero_data, POLAR_MATRICES)

Read the three (alpha × delta) coefficient matrices (see read_aero_matrix) and assemble the POLAR_MATRICES aero_data = (alpha, delta, cl, cd, cm).

read_aero_matrix(filepath::AbstractString) -> (Matrix{Float64}, Vector{Float64}, Vector{Float64})

Read an aerodynamic coefficient matrix from CSV with angle labels. Returns the coefficient matrix and corresponding angle ranges.

Returns

• matrix: Matrix of aerodynamic coefficients
• alpha_range: Vector of angle of attack values in radians
• delta_range: Vector of flap deflection angles in radians

read_dat(path) -> (x, y)

Read Selig .dat airfoil coordinates (two whitespace-separated columns), skipping the name/header line and any non-numeric lines.

read_node_table(path) -> (alpha, delta, values)

Read a per-node aero CSV (header alpha, delta, n0, n1, …; angles in degrees) into radian alpha/delta vectors (one entry per row) and a nrow × n_node value matrix.

delta_suffix(delta) -> String

Filename tag for a non-zero trailing-edge deflection delta [rad], e.g. d5, dm3 (−3°), d2p5 (2.5°). Uses millidegree precision (negatives as m, the decimal point as p) so sub-degree deflections get distinct .dat files instead of colliding. Shared by write_section_aero and read_section_aero.

interpolate_matrix_nans!(matrix::Matrix{Float64}; prn=true)

Replace NaN values in a matrix by interpolating from nearest non-NaN neighbors. Uses an expanding search radius until valid neighbors are found.

Arguments

• matrix: Matrix containing NaN values to be interpolated

remove_vector_nans(aero_data)

Remove the indices from aero_data where a NaN is found.

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)

Sweep over angle_range (degrees), solving at each angle. Returns a named tuple (polar_data, cmx, cmy, cmz, rey) where polar_data is [angle, cl, cd, cs, gamma_dist, cl_dist, cd_dist, cs_dist, reynolds].

extract_literature_polar_data(raw_data, path;
    angle_type="angle_of_attack")

Extract polar data from literature CSV data (as returned by readdlm). Returns a named tuple with polar_data (array of [angle, cl, cd, cs] vectors) and cmx, cmy, cmz vectors.

Supports both tuple format (table, header) from readdlm(...; header=true) and matrix format where the first row contains headers.

parse_literature_column(col)

Parse a column of mixed-type data (Real or String) into Float64, returning NaN for unparseable values.

interpolate_section_aero_to_refined!(wing)

Set each refined section's SectionAero by spanwise-interpolating the unrefined sections' surface tables (contour, Cp, cf) using the refined→unrefined mapping (refined_section_left_idx, refined_section_weight). No-op when the unrefined sections carry no surface aero.

validate_section_aero(sections)

Enforce the all-or-none surface-aero rule and a uniform node resolution: either every section in sections carries SectionAero or none does, and all present tables share the same node count. Throws ArgumentError on a violation.

copy_examples()

Copy all example scripts to the folder "examples" (it will be created if it doesn't exist).

bernstein_basis(x::AbstractVector, n::Int)

Compute Bernstein polynomial basis matrix.

Returns matrix of shape (length(x), n+1) where each column is a Bernstein polynomial B_i(x) = binomial(n,i) * x^i * (1-x)^(n-i)

class_function(x; N1=0.5, N2=1.0)

CST class function for airfoils: C(x) = x^N1 * (1-x)^N2

Standard values N1=0.5, N2=1.0 give round leading edge and pointed trailing edge.

leading_edge_basis(x::AbstractVector, n_weights::Int)

Leading-edge-modification basis used by AeroSandbox/NeuralFoil: x * (1 - x)^(n_weights + 0.5). Added identically to both surfaces.

compute_optimal_x_points(n) -> Vector{Float64}

NeuralFoil's boundary-layer station x/c: midpoints of a uniform [0,1] grid.

normalize_airfoil(x::Vector, y::Vector)

Normalize airfoil coordinates to unit chord with LE at origin.

Returns normalized (x, y) and transformation parameters.

get_lower_upper(x, y, crease_frac) -> (lower, upper)

Find y-coordinates where upper/lower surfaces intersect the hinge line.

turn_trailing_edge!(angle, x, y, lower_turn, upper_turn, crease_frac; thickness_frac=nothing)

Deflect airfoil trailing edge by rotating coordinates behind the crease line. Positive angle deflects downward.

The pivot's through-thickness location is set by thickness_frac (0 = bottom surface, 0.5 = mid, 1 = top): y_pivot = lower_turn + thickness_frac*(upper_turn - lower_turn). When thickness_frac is given, points are only rotated (no crease cleanup) — intended for a following Kulfan refit that wraps the crease robustly. When thickness_frac === nothing (legacy), the pivot is the top surface for downward deflection and the bottom for upward, and folded-over crease points are removed and the band averaged so the coordinates stay XFoil-paneable.

distance_parabolas!(d, f, n, v, z)

One pass of the Felzenszwalb–Huttenlocher distance transform: writes into d the lower envelope of the parabolas (q - p)^2 + f[p] over p, for q in 1:n. v and z are scratch (length n and n + 1).

squared_distance_transform!(field) -> field

In-place 2D squared Euclidean distance transform in grid-index units. On input field holds 0 (or a local squared offset) at seeds and a large finite value elsewhere; on output every node holds its squared distance to the nearest seed.

grid_sampler(field, x0, y0, cell) -> f(px, py)

Bilinear interpolant of field at physical points; node (i, j) sits at (x0 + (i-1)*cell, y0 + (j-1)*cell). Queries are clamped to the grid.

flood_outside(blocked) -> BitMatrix

Mark the nodes reachable from the grid border without entering blocked (4-connected flood fill); unreached nodes are the solid plus its enclosed holes.

trace_level_set(field, level) -> Vector{Vector{NTuple{2,Float64}}}

Marching-squares contours of field .== level, each returned as a closed loop of grid-frame vertices (unit = one cell, node 1 at 1.0). Saddle cells are resolved by the cell-center average. Loops touching the grid border are dropped.

largest_linking_gap(x, y) -> Float64

Largest edge of the Euclidean minimum spanning tree of the points (Prim, O(n²)) — the longest hop needed to keep the cloud connected. The rolling ball must not fall through it.

resample_arc(ax, ay, n, curvature_weight) -> (x, y)

Resample the polyline (ax, ay) at n stations cosine-clustered in a measure that blends arclength with curvature_weight extra length per radian of turning, so both ends attract points and sharp features are resolved by several panels regardless of their size (like XFoil's curvature-attracted paneling). Endpoints are preserved.

smooth_turning!(turn) -> turn

Diffuse the per-node turning-angle density along the contour (in place). A sharp corner — e.g. a deflected section's hinge, rounded to clearance — otherwise dumps its whole turn into one node, and the curvature-weighted resampling then collapses a few panels to near-zero length, which XFoil's viscous solver cannot handle. Spreading the turn over a short band refines a group of panels gradually instead, the way XFoil's PANGEN bunches panels on a smoothed curvature. The total turn is conserved, so leading-edge clustering is preserved.

load_neuralfoil_model(model_size::String="xlarge"; weights_dir=nothing)

Load NeuralFoil neural network weights from .npz files.

Arguments

• model_size: One of "xxsmall", "xsmall", "small", "medium", "large", "xlarge", "xxlarge", "xxxlarge"
• weights_dir: Directory containing the .npz files (defaults to package data dir)

Returns

• NeuralFoilModel: Loaded model ready for evaluation

neuralfoil_section(params, alpha, Re; kwargs...) -> NamedTuple

Full NeuralFoil evaluation returning integrated coefficients and the surface pressure distribution reconstructed from the predicted edge-velocity ratios (Cp = 1 - (ue/vinf)^2) at NeuralFoil's N fixed station x/c.

Returns (; alpha, cl, cd, cm, confidence, x, cp_upper, cp_lower), with x of length N and cp_upper/cp_lower sized N × n_alpha.

neuralfoil_fused_output(params, alpha, Re; model_size, weights_dir,
                        n_crit, xtr_upper, xtr_lower) -> Matrix

Forward pass with the symmetry embedding: evaluate the network, evaluate the top/bottom-flipped case, flip its outputs back, and average. Returns the fused output matrix (n_outputs × n_cases), with the Mahalanobis penalty already applied to the confidence logit (row 1).

nn_forward(x::AbstractMatrix, model::NeuralFoilModel)

Forward pass through the neural network.

Arguments

• x: Input matrix of shape (ninputs, ncases)
• model: Loaded NeuralFoil model

Returns

• Output matrix of shape (noutputs, ncases)

prepare_inputs(params::KulfanParameters, alpha, Re;
               n_crit=9.0, xtr_upper=1.0, xtr_lower=1.0)

Prepare neural network inputs from Kulfan parameters and flow conditions.

Arguments

• params: Kulfan CST parameters
• alpha: Angle of attack in degrees (scalar or vector)
• Re: Reynolds number (scalar or vector)
• n_crit: Critical amplification factor (default 9.0)
• xtr_upper: Forced transition location on upper surface (0-1)
• xtr_lower: Forced transition location on lower surface (0-1)

Returns

• Input matrix of shape (25, n_cases)

flip_inputs(x::AbstractMatrix)

Flip inputs for symmetry embedding (swap upper/lower, negate alpha).

flip_outputs(y::AbstractMatrix)

Flip outputs back after evaluating with flipped inputs.

squared_mahalanobis_distance(x::AbstractMatrix, model::NeuralFoilModel)

Compute squared Mahalanobis distance from training distribution.

This is used to penalize predictions far from the training data.

swish(x)

SiLU/Swish activation function: x * sigmoid(x) = x / (1 + exp(-x))

sigmoid(x)

Sigmoid activation function.

create_2d_polars(; dat_path, cl_polar_path, cd_polar_path, cm_polar_path, wind_vel,
                 area, width, crease_frac, alpha_range, delta_range,
                 solver=XFoilSolver(), remove_nan=true)

Generate 2D airfoil-section (alpha, delta) POLAR_MATRICES CSVs for an airfoil .dat, deflecting the trailing edge over delta_range and sweeping alpha_range (both radians) with solver. Pass solver=NeuralFoilSolver() to use NeuralFoil instead of XFoil — both emit the same three CSV files, so VSM loads them identically.

The Reynolds number is wind_vel * (area / width) / ν (kinematic viscosity of air). crease_frac is the chordwise hinge location (0–1). Writes lift, drag, and moment coefficient matrices to the three output paths.

lei_poly_coeffs(tube_diameter, camber) -> (cl_coeffs, cd_coeffs, cm_coeffs)

Breukels leading-edge-inflatable α-polynomial coefficients (α in degrees) for a section with normalized tube_diameter and camber. cl_coeffs is a cubic (4 coefficients), cd_coeffs/cm_coeffs are quadratic (3), each in ascending order (constant term first), ready to pass straight to evalpoly. Feed the result to a core POLY section as aero_data = (cl_coeffs, cd_coeffs, cm_coeffs).

resolve_airfoil(type, info, out_dir, id; Re, alpha_range) -> (new_type, new_info)

Resolve one awesIO wing_airfoils entry to a core-loadable form. breukels_regression (t, kappa) → poly coeffs (via lei_poly_coeffs); neuralfoil (dat_file_path, …) → a polars CSV (via generate_polar_from_dat) written under out_dir; polars/poly/inviscid pass through. masure_regression is not yet supported. info file paths should already be absolute.

read_dat_coordinates(path) -> (x, y)

Read airfoil coordinates from a Selig-format .dat file, skipping header and comment lines.

write_dat(filepath, name, x, y) -> filepath

Write airfoil coordinates to a Selig-format .dat file.

write_polar_csv(filepath, sols::Vector{SectionSolution})

Write a solver sweep (from any AbstractAirfoilSolver) to a POLAR_VECTORS CSV (alpha, Cd, Cs, Cl, Cm; alpha in degrees). Non-converged angles (NaN) are skipped, so this works for both NeuralFoil and XFoil sweeps.

write_polar_csv(filepath, result::NeuralFoilResult)

Write a NeuralFoil result to a POLAR_VECTORS CSV (alpha, Cd, Cs, Cl, Cm).

write_polar_matrix_csv(filepath, alpha_range, delta_range, cl, cd, cm)

Write an (alpha × delta) sweep to a long-format POLAR_MATRICES CSV with columns alpha, delta, Cl, Cd, Cm (both angles in degrees), one row per grid point. The delta column is what marks the file as a matrix polar to the loader. alpha_range and delta_range are in radians; cl/cd/cm are length(alpha) × length(delta).

write_aero_matrix(filepath, matrix, alpha_range, delta_range, label) -> filepath

Write an (alpha × delta) coefficient matrix to a labelled CSV: the header row holds the flap deflections (δ=…°), the first column the angles of attack (α=…°), both in degrees. alpha_range/delta_range are radians; label (e.g. "C_l") names the coefficient. The submodule-side writer; read_aero_matrix (main package) reads it back.

write_node_table(path, aero, values) -> path

Write a per-node aero table (Cp or cf, shaped n_node × n_alpha × n_delta) as a human-readable CSV: header alpha, delta, n0, n1, … (node columns in the contour node order), one row per (alpha, delta) with angles in degrees.

flat_plate_cf(xc, Re) -> Float64

Approximate local skin-friction coefficient at chord fraction xc (0..1) for Reynolds number Re (Re_x = Re·xc), taking the larger of the two standard flat-plate correlations: the Blasius laminar solution cf = 0.664·Re_x^(-1/2) and Prandtl's one-seventh-power-law turbulent estimate cf = 0.027·Re_x^(-1/7) (see e.g. White, Viscous Fluid Flow; Schlichting & Gersten, Boundary-Layer Theory). A stand-in per-node cf for backends that do not expose one (NeuralFoil); XFoil returns the exact distribution via bldump.

neuralfoil_contour_solution(alpha, res, i, x, y, le, cf) -> SectionSolution

Assemble a full-contour SectionSolution for case i of a NeuralFoil sweep res: interpolate res.cp_upper/res.cp_lower (at res.x) onto the contour nodes (x, y) split at the leading edge le, carrying the precomputed cf.

fill_node_nans!(grid, i)

Fill NaNs in node i's (alpha, delta) matrix grid[i, :, :] with interpolate_matrix_nans!; an all-NaN node is left untouched.

read_faces(filename)

Read vertices and faces from an OBJ file.

Arguments

• filename::String: Path to .obj file

Returns

• Tuple of (vertices, faces) where:
  • vertices: Vector of 3D coordinates [x,y,z]
  • faces: Vector of triangle vertex indices

slice_mesh_at_plane(vertices, faces, point, normal; tol=1e-6)

Slice a mesh with the plane through point with unit normal, returning the 3D segments (p1, p2) where the surface crosses it. The plane may have any orientation, so the cut can follow a curved or swept span.

order_segments_to_contour(segments; tol=1e-4)

Order line segments into a continuous contour.

Arguments

• segments: Vector of (p1, p2) tuples representing line segments

Returns

• Vector of [x, z] points forming the ordered contour

station_indices(arclen, n; wingtip_distance=0.0) -> Vector{Int}

Indices of the marched stations nearest n targets spread over the leading-edge arc length. The first and last targets sit wingtip_distance (arc length) in from the tips; with the default 0.0 they land exactly on the tips. Inset the tips a little (e.g. 0.1 m) to avoid degenerate near-zero-chord tip sections that some solvers (XFoil) cannot analyse.

airfoil_frame(LE_point, TE_point, span_tangent) -> (x_af, y_af, z_af)

Local airfoil axes for a slice, built so the slice plane contains the chord. With x̂ = [1,0,0] and the LE span_tangent: ẑ = x̂ × span (up), ŷ = ẑ × x̂ (spanwise slice normal — the span projected into the y-z plane), x̂_c = ŷ × ẑ (chord, oriented LE → TE). Slice the mesh with y_af. Returns nothing if the span is parallel to global-x.

build_section(vertices, faces, le, te, point, tangent) -> section or nothing

Slice the mesh at one marched station (airfoil_frame drops the chordwise tilt) and project it, producing the full (; LE_point, TE_point, span_dir, contour3d, x_airfoil, y_airfoil), or nothing.

contour_to_airfoil(contour::Vector{Vector{Float64}})

Convert a 2D contour to normalized airfoil coordinates.

Assumes contour is in [x, z] format where x is chordwise, z is thickness direction.

Returns

• (x, y): Normalized airfoil coordinates with chord = 1, LE at origin

plane_contour_to_airfoil(contour3d, LE_point, x_af, z_af)

Project a 3D slice contour into the local airfoil frame (chord axis x_af, up axis z_af) and return normalized Selig coordinates (x, y), or nothing if degenerate.

reorder_airfoil_selig(x::Vector, y::Vector)

Reorder airfoil coordinates to Selig format (TE upper -> LE -> TE lower).

This is the standard format expected by most airfoil tools.

densify_contour(contour, max_edge) -> Vector{Vector{Float64}}

Insert evenly spaced points along contour edges longer than max_edge. Coarse mesh triangles otherwise leave large hops in the slice cloud, which force the shrink wrap's auto-raised rolling ball far up and over-smooth the wrapped airfoil.

create_interpolations(vertices, circle_center_z, radius, gamma_tip)

Create interpolation functions for leading/trailing edges and area.

Arguments

• vertices: Vector of 3D point coordinates
• circle_center_z: Z-coordinate of circle center
• radius: Circle radius
• gamma_tip: Maximum angular extent

Returns

• Tuple of (leinterp, teinterp, area_interp) interpolation functions
• Where leinterp and teinterp are tuples themselves, containing the x, y and z interpolations

find_circle_center_and_radius(vertices)

Find the center and radius of the kite's curvature circle.

Arguments

• vertices: Vector of 3D point coordinates

Returns

• Tuple of (zcenter, radius, gammatip) where:
  • z_center: Z-coordinate of circle center
  • radius: Circle radius
  • gamma_tip: Angle of the kite tip from z-axis

march_edges(vertices, faces; step) -> (; le, te, point, tangent, arclen)

March the leading edge outward from mid-span in both directions in steps of arc length step. Each cut is a vertical spanwise plane (both the chordwise and vertical components of the running LE tangent dropped from the normal) so a tip that curls downward can't tilt the plane toward horizontal, where its min-chord "LE" pick would jump across the wing. Marching stops when the leading edge stops advancing spanwise. Cuts sample mesh edges, so the picks are robust to vertex density. Returns, ordered along the span, the LE/TE points, each cut's plane origin and tangent, and the cumulative LE arc length. Build the airfoil for a chosen station with build_section.

calculate_inertia_tensor(vertices, faces, mass, com)

Calculate the inertia tensor for a triangulated surface mesh, assuming a thin shell with uniform surface density.

Arguments

• vertices: Vector of 3D point coordinates representing mesh vertices
• faces: Vector of triangle indices, each defining a face of the mesh
• mass: Total mass of the shell in kg
• com: Center of mass coordinates [x,y,z]

Method

Uses the thin shell approximation where:

1. Mass is distributed uniformly over the surface area
2. Each triangle contributes to the inertia based on its area and position
3. For each triangle vertex p, contribution to diagonal terms is: area * (sum(p²) - p_i²)
4. For off-diagonal terms: area * (-p_i * p_j)
5. Final tensor is scaled by mass/(3*total_area) to get correct units

Returns

• 3×3 matrix representing the inertia tensor in kg⋅m²

center_to_com!(vertices, faces)

Calculate center of mass of a mesh and translate vertices so that COM is at origin.

Arguments

• vertices: Vector of 3D point coordinates
• faces: Vector of vertex indices for each face (can be triangular or non-triangular)

Returns

• Vector representing the original center of mass before translation

Notes

• Non-triangular faces are automatically triangulated into triangles
• Assumes uniform surface density

airfoils_from_yaml(geometry_file) -> Vector of (; id, x, y, x_raw, y_raw)

Read each airfoil's coordinates from the .dat files referenced by a YAML geometry's wing_airfoils. x, y come from dat_file (the fitted airfoil); x_raw, y_raw come from raw_dat_file if present (the raw sliced points), and are empty otherwise. Relative paths resolve against the geometry file's directory.

write_geometry_yaml(path, section_rows, airfoil_rows)

Write a geometry YAML via write_yaml, one line per section/airfoil row. section_rows are [airfoil_id, LE_x, LE_y, LE_z, TE_x, TE_y, TE_z]; airfoil_rows are [airfoil_id, type, info_dict] where info_dict holds dat_file, csv_file_path, and optionally raw_dat_file, cp_file, cf_file.

resolve_aero_geometry(yaml_in, out_dir; verbose=true) -> yaml_out

Read an awesIO-style geometry YAML and resolve every wing_airfoils entry to a core-loadable form via resolve_airfoil (breukels_regression → poly, neuralfoil → polars CSV, others pass through), writing generated CSVs and a resolved geometry.yaml under out_dir. wing_sections (incl. any VUP up-vectors) pass through unchanged. Load the result with Wing(yaml_out).

plot_airfoil_fit(x, y; kwargs...)

Plot a single airfoil with its Kulfan CST fit. Implemented in ObjAdapterMakieExt; load Makie (or GLMakie) to use it.

create_geometry_plot_makie(body_aero::BodyAerodynamics, title,
                           view_elevation, view_azimuth; zoom=1.8)

Create a 3D Makie plot of wing geometry including panels and filaments.

Arguments

• body_aero: struct of type BodyAerodynamics
• title: plot title
• view_elevation: initial view elevation angle [°]
• view_azimuth: initial view azimuth angle [°]

Keyword arguments

• zoom: zoom factor (default: 1.8)

plot_line_segment_makie!(ax, segment, color, label; width=3)

Plot a line segment in 3D with arrow using Makie.

Arguments

• ax: Makie Axis3
• segment: Array of two points defining the segment
• color: Color of the segment
• label: Label for the legend

Keyword Arguments

• width: Line width (default: 3)

set_axes_equal_makie!(ax, panels; zoom=1.8)

Set 3D Makie axis to equal scale based on panel data.

Arguments

• ax: Makie Axis3
• panels: Array of panels
• zoom: zoom factor (default: 1.8)

map_airfoil_3d(le, te, tangent, x, y) -> 3×N or nothing

Map normalized airfoil coordinates into 3D through the airfoil frame of the station given by its LE/TE points and leading-edge tangent (chord scale = |TE - LE|).

fitted_airfoil_3d(section, wrap_method; delta=0.0, crease_frac=0.75) -> 3×N or nothing

Shrink-wrap a section's sliced contour with wrap_method and map it back into 3D through the section's local airfoil frame, for overlaying on the 3D slice diagnostic. A nonzero delta (degrees) deflects the trailing edge and re-wraps (deform_section), showing the geometry the solvers consume. Returns nothing for a degenerate slice.

generated_slices(out_dir, delta, fit_pts) -> (slices, le, te)

Read the stations of a generated obj_to_yaml output directory and their written .dat airfoils — raw slice, wrap, and the delta-degree deformed wrap when it was generated — assembled for plot_slices_3d. Nothing is re-sliced or re-wrapped; only the Kulfan fits of the stored coordinates are recomputed (via fit_pts), exactly as the polar pipeline fits them.

airfoil_skin_geometry(body; R_b_w=nothing, T_b_w=nothing) -> (vertices, faces, ribs)

Lofted airfoil skin of a BodyAerodynamics: each section's deflected contour (section_surface at the panel's delta) is fitted between the panel's corner_points by a 2D similarity, TE pinned so a deflection bulges the fore body up. The skin reflects delta only when the geometry carries per-delta slices (obj_to_yaml with a delta_range); with δ=0-only data it renders undeflected, a deliberate cue that the deflected slices are missing. Transformed to world by R_b_w/T_b_w. vertices/faces triangulate the skin between consecutive equal-node sections; ribs is one closed contour polyline per section. Sections without contour data are skipped.

panel_normal(panel) -> Point3f

Body-frame airfoil-upper-surface unit normal of a panel, from its corner_points and sign-aligned to z_airf (both body-frame, so the sign relation survives the kite's attitude changes). Falls back to +z for a degenerate quad.

plate_hinge_local(crease_frac, delta) -> (lx, ly)

Chordwise/thickness fractions (chord = 1, LE at 0) of the flap hinge for a plate whose trailing edge is deflected by delta [rad] about crease_frac and then re-pinned so the TE stays at the panel TE corner. A downward deflection therefore lifts the hinge (ly > 0, the "bulge up"). crease_frac outside (0, 1) disables the kink (hinge stays on the chord line).

panel_plate_geometry(panel; R_b_w=nothing, T_b_w=nothing) -> Vector{Point3f}

The 6 vertices [LE_1, hinge_1, TE_1, TE_2, hinge_2, LE_2] of a panel's flat-plate skin, kinked at plate_hinge_local by the panel's delta/crease_frac (TE pinned to the corners, hinge bulging up). Triangulated by PLATE_FACES; with delta == 0 the two quads are coplanar (the original flat quad). Transformed to world by R_b_w/T_b_w.

plot!(ax, panel::VortexStepMethod.Panel; use_observables=false, kwargs...)

Plot a single Panel as a flat-plate mesh, kinked at the flap hinge by the panel's delta/crease_frac (see panel_plate_geometry); with delta == 0 this is the flat quad LE1-TE1-TE2-LE2.

If use_observables=true, creates observables for dynamic updates.

plot!(ax, body::VortexStepMethod.BodyAerodynamics; use_observables=false,
      airfoils=false, kwargs...)

Plot a BodyAerodynamics object. By default draws each panel as a flat quad; with airfoils=true instead draws the lofted airfoil skin (one see-through wing-shaped mesh with a contour rib line per section, see airfoil_skin_geometry).

If use_observables=true, creates observables for dynamic updates keyed by (bodyid, panelindex). Otherwise, creates static plots (original behavior).

plot!(body::VortexStepMethod.BodyAerodynamics; R_b_w=nothing, T_b_w=nothing)

Update existing body aerodynamics plot observables with current geometry. This updates all panels in the body using their current corner_points.

Requires that plot(body; use_observables=true) or plot!(ax, body; use_observables=true) was called first to create the observables.