airfoil_settings(data) -> AirfoilSettings

Read a wing's airfoil: block, checking that solver and table_format name backends that exist. Every key is optional and falls back to the default; delta_range: null (or an absent key) means no flap sweep.

settings_range(triple) -> StepRangeLen

A settings file's [first, step, last] as a range.

n_unrefined_sections(wing_settings::WingSettings) -> Int

Number of unrefined sections the Wing built from wing_settings carries: the rows of its geometry_file's wing_sections.

n_unrefined_sections(body_aero::BodyAerodynamics) -> Int

Number of unrefined sections summed over the wings of body_aero.

solver_kwargs(solver_settings::SolverSettings) -> NamedTuple

The Solver fields that solver_settings sets, as keyword arguments.

check_dimensions(solver::Solver, body_aero::BodyAerodynamics)

Throw a DimensionMismatch unless solver was built for as many panels and unrefined sections as body_aero has.

calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::Model, core_radius_fraction,
                        va_dist, va_unit_dist, target=body_aero.AIC)

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.

Both solvers converge on the fixed-point residual F(gamma) - gamma, measured relative to the largest circulation: the LOOP solver takes under-relaxed steps towards F(gamma), the NONLIN solver a Newton step with a finite-difference Jacobian, backtracked until it reduces the residual.

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_dist, 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_vec_dist)

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_vec_dist::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!.

finite_full(x) -> Bool

true if x is finite. A ForwardDiff.Dual needs every partial finite too, so a non-finite derivative of a solve is caught as well as a non-finite value.

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).

set_pitch_rate_dist!(body_aero, omega)

Fill body_aero.pitch_rate_dist from a rigid-body turn rate by projecting it onto each panel's own spanwise axis. Panels with different dihedral see different rates from the same omega.

coeffs_at_angles(solver, body_aero, alpha, beta, va)

Aerodynamic coefficients [CFx, CFy, CFz, CMx, CMy, CMz] of body_aero solved at angle of attack alpha [rad], sideslip beta [rad] and apparent wind speed va [m/s], at the rotation rate body_aero.omega, which it turns about and stores as body_aero.reference_point. Throws a SolveFailure if the solve misses the solver's tolerances.

nose_down(solver, body_aero, alpha, beta, va)

Whether CMy from coeffs_at_angles is negative.

bisect_sign_change(predicate, low, high, tol)

Bisect [low, high], across which the boolean predicate flips, to a width of tol and return the midpoint.

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!(alpha_corrected, body_aero::BodyAerodynamics, gamma,
                                  core_radius_fraction, va_vec_dist, va_dist,
                                  va_unit_dist)

Write each panel's inflow_angle at its aerodynamic centre into alpha_corrected, the induced velocity taken from gamma on the LLT influence matrix.

induced_velocity_at(body_aero::BodyAerodynamics, point, gamma, core_radius_fraction)

Velocity [m/s] induced at point by the horseshoe vortices of all panels, panel j carrying circulation gamma[j] [m²/s] and its frozen wake.

attached_trailed_loads(body_aero::BodyAerodynamics, i, gamma, density,
                       core_radius_fraction, reference_point)

Kutta–Joukowski force [N] on the two chordwise trailed vortex segments of panel i, from its quarter-chord bound points to its trailing edge, and its moment [N·m] about reference_point, as (; force, moment). Each segment carries gamma[i] and sees, at its three-quarter-chord point, the inflow turned by body_aero.omega plus the induced_velocity_at that point.

panel_force_moment(body_aero::BodyAerodynamics, i, loads, y_airf, gamma, density,
                   core_radius_fraction, reference_point,
                   is_with_attached_trailed_force)

Force [N] and moment [N·m] about reference_point of panel i: its section loads from panel_loads acting at the aerodynamic centre, plus attached_trailed_loads when is_with_attached_trailed_force. Returns (; force, moment).

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

wing_span_flip(wing) -> Int8

-1 when wing's sections run against its spanwise_direction, +1 otherwise: the flip every panel of the wing is reinitialized with (reinit!). One answer per wing, so neighbouring panels cannot disagree and invert a single normal by 180°.

spanwise_extent(wing::AbstractWing)
spanwise_extent(wings, spanwise_direction)

Lowest and highest projection of the unrefined sections' LE and TE points on wing's spanwise_direction, or of all wings on spanwise_direction, as (lo, hi) [m].

calculate_circulation_distribution_elliptical_wing(gamma_i, body_aero::BodyAerodynamics,
                                                   gamma_0=1.0)

Write into gamma_i an elliptic circulation of peak gamma_0 over each wing, its control points measured along that wing's spanwise_direction from the wing's mid-span.

_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_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.

Row i of an n×3 matrix as a 3-vector.

Length smooth_norm adds in quadrature to keep a norm positive at zero.

smooth_norm(vec, floor=SMOOTH_FLOOR)

Norm floored in quadrature: positive and differentiable at the origin, and branch-free, which is what makes the functions below traceable.

panel_span_vector(le_1, te_1, le_2, te_2)

Quarter-chord vector from section 2 to section 1: the panel's spanwise axis, its length the span width.

panel_chord_weight(width_prev, width_own, width_next)

Section 1's share of the edge blend setting a panel's chord direction, weighted by panel spacing so a panel leans towards a narrower neighbour. nothing for a missing neighbour at a tip; a lone panel gets 0.5.

panel_axes(le_1, te_1, le_2, te_2, chord_weight=0.5, orient=1)

Airfoil frame and size of the panel between two sections, as (; x_airf, y_airf, z_airf, chord, width). x_airf is chordwise, y_airf runs along the quarter-chord line the bound vortex sits on — not square to x_airf on a swept panel — and z_airf is x_airf × y_airf normalised. chord_weight (panel_chord_weight) is section 1's share of the chord-direction blend; orient is ±1 and flips y_airf and z_airf together.

panel_axes(panel::Panel)

The panel's stored airfoil frame and size in the shape panel_axes returns, so a panel built by the geometry pass feeds the same functions as one built straight from section points.

inflow_angle(velocity, y_airf, z_airf)

Angle of attack of velocity in the plane square to y_airf: from the chord's projection on that plane, y_airf × z_airf, towards z_airf.

effective_alpha(alpha, deficiency)

Angle the polars are read at: the geometric inflow angle less an unsteady lag. The geometric angle still turns the force, so a lag shifts the coefficients only.

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 result. dva_vec_1/dva_vec_2 are the sections' trailing minus leading edge apparent wind, giving the section_pitch_rate; nothing leaves it zero. deficiency feeds effective_alpha.

dynamic_pressure(rho_1, rho_2, v_span)

Panel dynamic pressure from its two section densities and the panel_inflow span-wise velocity, as that vector or its magnitude.

flow_curvature_cm(pitch_rate, chord, v_rel)

Thin-airfoil quarter-chord moment increment of a section pitching about its own spanwise axis: Δcm = -(π/4)·q̂ with q̂ = q·c/(2·v_rel), q positive nose-up. Pivot-independent, and lift needs no matching correction because the inflow is already sampled at three-quarter chord.

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_rel, the speed normal to the span; v_span is the velocity along y_airf.

panel_force_directions(axes, alpha_dir, spanwise)

Lift and drag unit vectors, (; dir_lift, dir_drag). alpha_dir is the angle that turns the force, spanwise the wing's spanwise direction.

prescribed_va_directions(va, spanwise)

Lift and side unit vectors, (; dir_lift, dir_side), of inflow va on a wing along spanwise: lift normal to both, side normal to lift and va.

Panel pitching moment per unit span about y_airf, positive nose-up.

panel_couple_force(cm, q_dyn, chord, width, scale=1)

panel_moment as a force couple: the magnitude two opposed forces one chord apart need to produce it. A particle model places the moment this way, along the panel normal at its leading and trailing edge.

panel_loads(axes, dirs, q_dyn, cl, cd, cm, scale=1; c_span=0)

Panel load from its polar coefficients, panel_axes and panel_force_directions, as (; lift, drag, moment, force, pitching_moment). The first three are per unit span; force and pitching_moment are the whole panel's, scale included. c_span adds a force along y_airf, as spanwise_flow_drag gives it.

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, with a core radius of core_radius_fraction times the filament length.

velocity_3D_trailing_vortex!(vel, filament::BoundFilament,
                             XVP, gamma, va, work_vectors)

Calculate induced velocity by a trailing vortex filament, with a Lamb–Oseen core radius grown over the axial distance of XVP from the filament start.

Arguments

• XVP: Control point coordinates
• gamma: Vortex strength
• 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 as implemented in KiteAeroDyn".

velocity_3D_vortex_segment!(vel, filament::BoundFilament, XVP,
                            gamma, epsilon, work_vectors)

Calculate the Biot–Savart velocity induced by a straight vortex segment at XVP. Inside the core radius epsilon the velocity is evaluated on the core boundary and scaled linearly with the distance to the axis.

velocity_3D_trailing_vortex_semiinfinite(filament::SemiInfiniteFilament, 
                                         Vf, XVP, GAMMA, va, 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::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::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; kwargs...)

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

flip reverses the section order so y_airf points along +spanwise_direction and z_airf to the airfoil upper surface, making the aero independent of section ordering. The caller owns it and must not derive it from the live geometry, which would invert normals mid-run; reinit!(::BodyAerodynamics) decides it once per wing.

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.

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! afterwards to update the panels.

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.

panel_range(body_aero::BodyAerodynamics, wing_idx)

Indices of the panels of wing wing_idx in body_aero.panels.

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). Sections within 1e-9 of the closest distance, relative, tie, and a tie goes to the one farther along spanwise_direction from the middle of the span, so a mirror-symmetric wing maps mirror panels to mirror sections. Works after refinement is complete.

section_centre_dist2(section, x, y, z)

Squared distance from the centre of section (the midpoint of its LE and TE) to (x, y, z).

section_span_position(wing, section)

Position of the centre of section along wing.spanwise_direction.

compute_refined_section_interpolation!(wing::AbstractWing;
                                       reuse_aero_data=false)

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).

reuse_aero_data keeps the surface tables the refined sections already hold instead of reblending them from the unrefined ones, the same preservation refine! applies to aero_data under use_prior_polar. A remesh that replaces the unrefined sections with a coarser set would otherwise resample the surface tables down to that set even while the polars stay at full resolution.

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.

The rotation is a full Rodrigues rotation about the spanwise axis. A swept or dihedral section has a chord component along that axis, and dropping the axial term would scale it by cos(theta), shortening the chord and tilting it.

_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.

span_order_key(section) -> Float64

Spanwise coordinate that normalize_span_order! and refine!'s sort_sections order sections on.

normalize_span_order!(sections) -> sections

Reverse sections if they do not already run +y to -y, the order panel normals are built from. Use refine!'s sort_sections for a scrambled list.

can_reuse_prior_refined_surface_tables(wing) -> Bool

Whether every refined section already carries a SectionAero, so a remesh can keep them instead of reblending from the unrefined sections. False on a wing that has none, where the blend is what fills them in the first 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.

set_polar!(panel, alphas, cl, cd, cm; shape=nothing)

Rewrite a panel's polar table with values at alphas [rad], ascending, and rebuild its interpolations. The table keeps the shape the panel was built with: a POLAR_VECTORS panel takes the angles as they are, and a POLAR_MATRICES panel takes them at every delta it already spans, since a regenerated polar carries its deflection as shape rather than as a flap angle and so says the same thing at each. The panel's alpha_ref and alpha_window are taken from the angles, so a table covering only a window around one angle is held at its ends rather than extrapolated past them (see window_alpha).

Written in place: the knots and the three value vectors reuse the panel's own alpha_knots and cl_coeffs/cd_coeffs/cm_coeffs storage whenever the sample count is unchanged (polar_column!), and the interpolations read that storage rather than a copy of it, so writing it is the whole update and the interpolation objects stay put (refresh_polar!). Only a table that changes shape is rebuilt. Values may be views into a caller's own buffers; nothing here holds on to them.

shape is the KulfanParameters the values were generated from, stored on the panel as live_shape so a panel's polar and the shape behind it are set together and cannot drift apart. The panel holds the shape itself, not a copy, so a source that rewrites one in place has already updated every panel flying it.

refresh_polar!(old, knots, values) -> Extrapolation

An interpolation reading knots and values themselves, carrying old's extrapolation and knot container so it keeps the type the panel's field was parameterised with. interpolate! takes both arrays by reference, so once an interpolation reads a panel's own storage, writing that storage is the whole update and there is nothing to rebuild. A table that changes shape is rebuilt once, and reads in place from then on. The two dimensional form spreads the angles over every delta the panel spans.

polar_column!(stored, values) -> Vector{Float64}

One column of a rewritten polar table in the panel's own stored vector, which is returned, or in a fresh vector when the sample count has changed — the one case a rewritten table has to grow, and the one that costs the panel's interpolations a rebuild.

polar_model(interp) -> AeroModel

The model a rewritten table leaves the panel on: the one its interpolations were built with, since their shape is fixed by the panel's type.

Whether an interpolation's knot container is a view of knots rather than a copy.

The angles in the container the panel's interpolations were parameterised with.

polar_knots(alphas) -> AbstractVector

The angles of a 1D polar in the container whose knot search fits its length, see ScanKnots.

window_alpha(panel, alpha)

alpha [rad] clamped into the range the panel's polar was built over. A table generated over a window around one angle of attack says nothing past its ends, so it is held at its end values rather than extrapolated out of them; alpha_window of 0 means unbounded and leaves alpha alone, which is what a full-range polar wants.

assemble_polar_matrix(alpha, delta, coefficients) -> (alphas, deltas, grids...)

Scatter long-format polar rows (radian alpha/delta per row, one coefficient per column of coefficients) onto the sorted unique alphas × deltas grid, one matrix per coefficient. Grid points without a row stay NaN.

load_polar_data(path::String) -> (aero_data, model)

Load an airfoil polar table, CSV or Arrow as the suffix of path says (see read_node_table), with columns alpha, cl, cd and cm (case-insensitive, any order, alpha in degrees). A delta column makes it a long-format POLAR_MATRICES table, returned as (alpha, delta, cl, cd, cm) on its sorted grid; without one it is POLAR_VECTORS (alpha, cl, cd, cm). Angles come back in radians. A missing, empty or malformed file warns and returns (nothing, INVISCID).

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_coordinates(path) -> (x, y)

Read Selig .dat airfoil coordinates (two whitespace-separated columns), skipping the name/header line, comments, and any row that is not a finite pair. The single .dat reader; write_dat is its writer.

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

Read an aero table into a radian alpha vector (one entry per row), a radian delta vector or nothing when the table has no delta column, the nrow × ncol matrix of the value columns and their names. The file suffix picks the format: .arrow (angle columns, a per-row list column values, the names in the columns metadata entry) or anything else CSV (a header line naming every column, found case-insensitively in any order). Angles are stored in degrees.

read_arrow_columns(path) -> (names, data)

Read an Arrow aero table into its column names and one Float64 matrix: the angle columns, then the list column values spread under the names its columns metadata entry holds, n0, n1, … where it has none.

read_csv_columns(path) -> (names, data)

Read a CSV table with one header line into its column names and one Float64 matrix. Throws on a row whose field count differs from the header's.

write_node_rows(path, alpha, delta, values; columns) -> path

Write an aero table from a radian alpha vector (one entry per row), a radian delta vector or nothing for a table without a delta column, and a nrow × ncol value matrix whose columns are named columns (default n0, n1, …). The file suffix picks the format, matching read_node_table: .arrow, or anything else CSV at 16 significant digits. Angles are written in degrees at 16 significant digits either way.

csv_fields(values) -> String

Join values into one comma-separated CSV line at 16 significant digits.

convert_node_table(src, dst) -> dst

Rewrite an aero table, a per-node Cp/cf table or a polar, in the format dst's suffix names, keeping its column names. Use to move an existing dataset between CSV and Arrow without re-running the (slow) airfoil solver that produced it.

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, 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 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).

deform_kulfan!(out, basis, base, upper_deflection, lower_deflection, residual)
    -> out

deform_kulfan written into the shape out already is, with residual as the scratch chord_residual! fills — the form a live polar refreshes a panel's shape with, since it neither allocates nor replaces the object the panel points at.

chord_residual!(residual, basis, deflection) -> residual

chord_residual written into storage the caller owns.

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.

truncated_least_squares(A, b)

Minimum-norm least-squares solution of A * coeffs = b with the singular values of A below KULFAN_FIT_RTOL times the largest dropped. Returns (coeffs, n_dropped).

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)

Heights of the lower and upper surfaces of the contour (x, y) at the chord station x = crease_frac: the lowest and highest crossing of that vertical line with the segments between consecutive points, not including the one from the last point back to the first. Throws an ArgumentError when the contour crosses it fewer than twice.

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.

point_buckets(x, y, reach) -> Dict

Sparse uniform grid of the point indices, one bucket per reach-sized cell, so every point within reach of a query lies in one of the nine buckets around it.

turn_measure(cross, dot) -> Float64

Pseudo-angle in [0, 4) of the direction with the given cross and dot against a reference, rising monotonically with the counterclockwise angle it stands for. Ordering directions by it is ordering them by angle, without the atan that dominates the pivot's inner loop.

pivot_step(x, y, buckets, r, p, centre) -> (q, centre) or nothing

Rotate the empty disk of radius r counterclockwise about point p, starting from centre, until its rim reaches a second point. Returns that point and the disk centre touching both, or nothing when nothing lies within 2r of p.

push_arc!(px, py, cx, cy, radius, from, sweep)

Append the arc about (cx, cy) that starts at angle from and turns by the signed angle sweep, endpoints included, stepped fine enough to stay within 1e-6 of the true arc. A zero radius appends the centre once.

edge_normal(ax, ay, bx, by) -> Float64

Angle of the outward normal of the wrap edge running from (ax, ay) to (bx, by), which pivot_contour traces with the cloud on its left.

pivot_contour(x, y, r, clearance) -> (px, py)

Trace the alpha shape of the cloud (x, y) for alpha r, offset outward by clearance, as a closed polyline. A disk of radius r is pivoted around the outside of the cloud from its leftmost point; the points it touches are the polygon's vertices in order, the ones inside a concavity narrower than the disk having been skipped and bridged straight. The offset rounds each convex vertex with an arc of radius clearance and chamfers each reflex one across the bisector, so it holds that distance from every contact, and the loops the offset makes where contacts lie closer than clearance are cut out (cut_loops).

cut_loops(px, py) -> (x, y)

The closed polyline (px, py) walked from its first node with every loop it makes by crossing itself cut out, the crossing point taking the loop's place, so what comes back is a simple closed curve.

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.

densify_contour(contour, max_edge) -> Vector

Insert evenly spaced points along the edges of the point sequence contour longer than max_edge, at most 100 pieces per edge.

enforce_min_spacing!(stations, floor_length) -> stations

Push the sampling stations apart until no gap is shorter than floor_length, taking the room out of the gaps that have it to spare so the first and last station stay put. floor_length is capped at the mean gap, which a uniform spacing already meets.

resample_arc(ax, ay, n, curvature_weight; min_panel_fraction=0.25) -> (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). The turning is taken as a curvature density smoothed over a few panel lengths (smoothed_curvature), so the measure grows smoothly and neighbouring panels stay comparable in length. No panel is left shorter than min_panel_fraction of the mean (enforce_min_spacing!). Endpoints are preserved.

smoothed_curvature(node_arclength, turn, band) -> Vector{Float64}

Turning per unit arclength at each node, averaged over a window of length band centred on it. Averaging over a length leaves the measure independent of how densely the trace happens to be stepped, which varies by a factor of several between the clearance arcs and the straight edges between them.

crossing_panels(x, y) -> Tuple{Int,Int} or nothing

The first pair of non-neighbouring panels of the contour (x, y) that cross, or nothing when the contour is a simple closed curve. A contour that does not end on its first node is closed by a panel back to it, which carries the highest panel index.

segments_cross(p, q, r, s) -> Bool

Whether the segments p-q and r-s cross properly, each strictly separating the other's endpoints. Touching at an endpoint or lying along each other does not count, nor do segments whose coordinate extents do not overlap.

extents_overlap(a1, a2, b1, b2) -> Bool

Whether the intervals spanned by a1, a2 and by b1, b2 overlap.

side_of_line(a, b, p)

Twice the signed area of the triangle a, b, p: positive with p left of the line from a to b, negative right of it, zero on it.

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, ue_upper, ue_lower), with x of length N and the four matrices sized N × n_alpha. A reconstruction should interpolate ue and square afterwards: it is linear in arc length through a stagnation point, where Cp is quadratic.

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).

fused_output(x, model) -> Matrix

Symmetry-fused forward pass over a prepared input matrix, see neuralfoil_fused_output. Takes the inputs already built so a caller that assembles its own batch does not go back through prepare_inputs.

fused_output!(work, x, model) -> AbstractMatrix

fused_output run entirely inside work, allocating nothing for a batch the workspace has room for. The result is a view of the workspace's own storage and stays valid until the next pass through it.

decode_surface_velocity(y) -> (station_x, ue_upper, ue_lower)

The edge-velocity ratios a fused output matrix carries, at NeuralFoil's own N fixed stations. Each ue matrix is N × n_cases. Interpolate these rather than Cp = 1 - ue²: ue is linear in arc length through a stagnation point, where Cp is quadratic.

decode_surface_velocity!(ue_upper, ue_lower, y) -> (ue_upper, ue_lower)

decode_surface_velocity into matrices the caller owns. The station axis is fixed by the network and is carried by NeuralFoilWorkspace instead, so a live pressure refresh reads a pass out without allocating anything for it.

surface_velocity_rows(y) -> (n_stations, upper, lower)

The number of boundary-layer stations a fused output matrix carries and the row ranges its upper and lower edge-velocity ratios sit in. Viewing those rows is how a caller that already holds the output reads the velocities without copying them out (decode_surface_velocity is the copying form).

decode_coefficients(y) -> (cl, cd, cm, confidence)

Turn a fused network output matrix into the integrated coefficients per case, undoing NeuralFoil's output scaling. The single place that scaling is written down.

decode_coefficients!(cl, cd, cm, confidence, y) -> (cl, cd, cm, confidence)

decode_coefficients into vectors the caller owns, so a live polar reads a refresh out of the network without allocating four vectors for it.

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)

nn_forward!(layers, x, model, n_cases=size(x, 2)) -> AbstractMatrix

Forward pass writing each layer's activations into layers, returning a view of the last one over the n_cases columns used. The matrix products go through BLAS in place, so a pass over a batch the buffers were sized for allocates nothing.

add_bias!(out, bias, activate) -> out

Add a layer's bias to its activations in place, passing them through swish unless this is the output layer, which carries no activation.

layer_buffers(model, n_cases) -> Vector{Matrix{Float32}}

One activation matrix per network layer, n_cases wide: the storage a forward pass writes through. NeuralFoilWorkspace holds two sets of them, one per symmetry.

How many cases a workspace was sized for.

prepare_inputs(params::KulfanParameters, alpha, Re;
               n_crit=9.0, xtr_upper=1.0, xtr_lower=1.0)
prepare_inputs(params::AbstractVector{KulfanParameters}, alpha, Re; kwargs...)

Prepare neural network inputs from Kulfan parameters and flow conditions. The vector form takes one shape per case, so a whole wing's panels go through the network in a single forward pass.

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)

fill_case_input!(x, case, params::KulfanParameters, alpha_deg, Re, n_crit,
                 xtr_upper, xtr_lower)

Write one NeuralFoil input column: rows 1–18 the Kulfan shape, rows 19–25 the flow condition (alpha_deg in degrees). The single point where the network's input layout is defined, shared by every prepare_inputs method.

flip_inputs(x::AbstractMatrix)

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

flip_inputs!(x_flip, x) -> x_flip

flip_inputs written into storage the caller owns: the upper and lower weights swap with a sign change, the leading-edge weight and sin(2·alpha) negate, and the two transition locations swap. Everything else carries over.

flip_outputs(y::AbstractMatrix)

Flip outputs back after evaluating with flipped inputs.

flipped_row(row, n_stations) -> (source, sign)

Which row of a network output a top/bottom-flipped case reads row from, and with which sign. The single place the mirror symmetry of the output layout is written down, shared by flip_outputs and fuse_flipped!.

CL and CM change sign, the two transition locations swap, and of the six n_stations-long boundary-layer blocks the upper and lower halves swap — the momentum thickness and shape factor as they are, the edge velocity mirrored, so it changes sign.

fuse_flipped!(y, y_flip) -> y

Average a direct output with its flipped counterpart in place, undoing the flip on the way (flipped_row). Reading y_flip while writing y is what lets the fusion land in storage that is already there instead of a third matrix.

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

Compute squared Mahalanobis distance from training distribution, one per case.

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

mahalanobis_case(x, case, model, centered) -> Float64

One case's squared Mahalanobis distance, taking centered as the scratch the centred input column is written into rather than allocating one per case.

penalize_confidence!(y, x, model, centered) -> y

Subtract every case's Mahalanobis penalty from the confidence logit, row 1 of the network output. This is what makes a shape far from the training distribution report a low confidence, and it is applied to each symmetry before the two are fused.

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, table_format=:csv)
    -> (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 table {id}.{table_format} (:csv or :arrow, 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.

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

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

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

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

write_polar(filepath, result::NeuralFoilResult)

Write a NeuralFoil result to a POLAR_VECTORS table (alpha, Cd, Cs, Cl, Cm), CSV or Arrow as the suffix of filepath says (see write_node_rows).

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

Write an (alpha × delta) sweep to a long-format POLAR_MATRICES table with columns alpha, delta, Cl, Cd, Cm (both angles in degrees), one row per grid point, CSV or Arrow as the suffix of filepath says. 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), one row per (alpha, delta) with angles in degrees. The file suffix picks the format: .arrow (columns alpha, delta and a per-row list column values, an order of magnitude faster to load), anything else a human-readable CSV with header alpha, delta, n0, n1, … (node columns in the contour node order). read_node_table reads both.

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, arc, cf) -> SectionSolution

Assemble a full-contour SectionSolution for case i of a NeuralFoil sweep res: reconstruct Cp on the contour nodes from the predicted edge velocities (contour_pressure), carrying the precomputed cf.

contour_arc(x, y, le) -> Vector

Signed arc length of every contour node from the leading edge node le: positive along the upper surface toward the trailing edge, negative along the lower. The natural coordinate across a blunt nose, where a small chordwise step is a long one along the skin.

contour_arc!(arc, x, y, le) -> arc

contour_arc written into arc, which is resized to the contour. A live pressure refresh walks panel after panel through one buffer this way.

arc_at_chord(x, arc, indices, fractions) -> Vector

Signed arc length at each chord fraction, read off the contour nodes indices — the one surface the fractions belong to.

arc_at_chord!(out, scratch, x, arc, indices, fractions) -> out

arc_at_chord written into out, taking the sorted copy of the surface from scratch rather than allocating one per call.

sorted_surface!(scratch, x, arc, indices) -> (surface_x, surface_arc)

One surface's nodes as chord fraction and signed arc length, ascending in chord, in the scratch's own storage. Sorted by insertion, which a surface listed from one end to the other is either already in or exactly reversed from, and which needs no scratch of its own.

sort_pairs!(sorted, carried)

Insertion-sort sorted ascending, carrying carried along with it. Stable, in place, and linear on input that is already ordered, which both of its callers hand it.

trailing_edge_speed(ue_upper, ue_lower, upper_arc, lower_arc, upper_te, lower_te)
    -> Float64

One edge speed for both surfaces at the trailing edge: each surface's last two stations extrapolated to its own trailing-edge arc length, the two magnitudes averaged, so a sharp edge carries a single pressure (Kutta).

velocity_knots(station_x, ue_upper, ue_lower, x, arc, le) -> (knots, values)

The whole contour's edge velocity as one curve of (signed arc length, signed ue/u∞), strictly increasing in arc length.

NeuralFoil reports nothing over the first and last 1/2N of chord. Signing the lower surface negative joins both surfaces into one curve through the stagnation point, so the unsampled nose is interpolated between the innermost station on each side rather than extrapolated off the end of one. The stagnation point enters as its own knot, where ue interpolated linearly between those two stations reaches zero — the stagnation-point-flow result, ue being linear in arc length there.

velocity_knots!(scratch, station_x, ue_upper, ue_lower, x, arc, le)
    -> (knots, values)

velocity_knots assembled in scratch. The two vectors returned are the scratch's own and are overwritten by the next call.

contour_pressure(station_x, ue_upper, ue_lower, x, arc, le) -> Vector{Float64}

Cp at every node of a closed contour x with signed arc length arc and nose at le, from one case's edge-velocity ratios at NeuralFoil's station_x. Builds the one-curve edge velocity (velocity_knots), interpolates it with a shape-preserving monotone cubic in arc length, and squares it into Cp = 1 - ue².

contour_pressure!(cp, scratch, station_x, ue_upper, ue_lower, x, arc, le) -> cp

contour_pressure written into cp, with scratch carrying the edge-velocity curve. Only the monotone interpolation itself still allocates, once per contour.

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.

live_xfoil_solver(live::LivePolars) -> XFoilSolver

XFoil settings matching the ones the live polars were evaluated at: the same critical amplification factor, and free transition on both surfaces, which is what refresh_live_polars! feeds the network. A comparison run at other settings measures the settings rather than the network.

xfoil_ramp(alpha, step) -> Vector{Float64}

Angles [rad] from zero out to alpha in step increments, alpha itself last. The viscous march refuses a blunt inflated section jumped straight to its angle, and analyze_sweep reinitialises at zero and walks outward, so handing it the whole ramp is what gets the reference solved at all.

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(march, n; wingtip_distance=0.0, min_chord_frac=0.01) -> Vector{Int}

Indices of the march_edges stations nearest n targets spread evenly over the spanwise length of the quarter-chord line: its arc length with the chordwise x component dropped. Stations whose chord has closed to less than min_chord_frac of the longest one are left out of that range first, so a wing tapering to a point puts its outermost sections on the last stations that still have an airfoil to slice. The remaining first and last targets sit a further wingtip_distance [m] inboard.

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.

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)

Cut the mesh at mid-span and march the leading edge outward from there in both directions with march_stations. Returns, ordered along the span, the LE/TE points, each cut's plane origin, and the LE tangent: the central difference over the neighbouring stations, where the end stations reuse their inner neighbour's. Build the airfoil for a chosen station with build_section.

cut_station(vertices, faces, point, tangent) -> (; le, te) or nothing

Cut the mesh with the vertical spanwise plane through point whose normal is the spanwise component of tangent, and return the crossings with the least and greatest x as the leading and trailing edge; nothing where the plane misses the mesh.

tip_station(vertices, faces, prev_le, tangent, step) -> (; le, te, point) or nothing

Bisect the step of length step from prev_le along tangent for the outermost cut that still yields an airfoil section whose leading edge stays within step of prev_le chordwise; nothing where no such cut exists.

march_stations(vertices, faces, start_le, step, direction) -> Vector{NamedTuple}

March the leading edge from start_le in steps of arc length step towards +y for direction = 1.0 or -y for -1.0, returning one (; le, te, point) per station, the tip station last. Stops when the leading edge stops advancing spanwise.

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, polar_file_path, and optionally raw_dat_file, cp_file, cf_file.

resolve_aero_geometry(yaml_in, out_dir; table_format=:csv, 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 table, others pass through). The generated polars go to out_dir/polars as table_format (:csv or :arrow), the resolved YAML to out_dir/geometry.yaml. wing_sections (incl. any VUP up-vectors) pass through unchanged. Load the result with Wing(yaml_out).

table_path_prefix(geometry_path, output_dir) -> String

Path from the geometry YAML's directory to the table directory, empty when they are the same. Table references resolve against the YAML's own directory, so a YAML written outside output_dir has to carry this hop.

prefix_table_paths!(airfoil_rows, prefix) -> airfoil_rows

Prepend prefix to every relative table reference in the info_dict of each row, so a geometry YAML written outside the table directory still resolves them. A no-op on an empty prefix.

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.

migrate_node_tables(yaml_path, output_dir, table_format; verbose=true)

Rewrite a generated dataset's per-node Cp/cf tables and polars in table_format and point geometry.yaml at them, when they are not in that format already; a legacy csv_file_path key is written back as polar_file_path. This is what lets an existing directory change format without re-running the airfoil solver that produced it. The source tables are left in place.

output_dir is what the YAML's relative table references resolve against, which is its own directory — the rule the geometry loader follows.

display_named(fig, name) -> fig

Show fig in the window registered under name, opening a new titled window for a name that has none yet: plotting the same title again replaces its predecessor, while a new title gets its own window. Backends without titled screens, such as CairoMakie, fall back to a plain display.

span_axis(position, title, ylabel) -> Axis

Axis for a spanwise distribution, +y on the left, matching the kite seen from the front and the +y to -y order its sections are stored in.

create_geometry_plot_makie(body_aero::BodyAerodynamics, title,
                           view_elevation, view_azimuth; zoom=0.5, show_title=true)

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: 0.5)
• show_title: draw title above the axis (default: true)

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)

polar_style(idx, n_solvers)

Marker, marker size and line style of the idx-th polar: solver output for idx <= n_solvers, literature otherwise.

widen_flat_ylims!(ax, series; min_span=0.1)

Set the y-limits of ax to span at least min_span around the non-NaN values of all vectors in series.

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 contour is fitted between the panel's corner_points by a 2D similarity, TE pinned so a deflection bulges the fore body up.

The contour is the panel's own live_shape when it has one — the very KulfanParameters object the live polar source deformed and handed to the airfoil solver, not a re-derivation of it, so a deformation bug shows up in the picture instead of being papered over. Otherwise it is section_surface at the panel's delta, which 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_contour(panel, section) -> (x, y)

The airfoil contour to draw a panel with: coordinates of its live_shape when a live polar source has put one there, else the section's tabulated contour at the panel's delta. The live branch reads the stored shape object itself, so what is drawn is what was flown.

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.

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.

PLATE_FACES

The 4 triangles that mesh the 6 panel_plate_geometry vertices of a panel's flat-plate skin, two per side of the hinge.

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, outlined in border_color at border_linewidth.

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 outlined in border_color; 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.