API Reference

This page provides the automatically generated API documentation for the lightaero library.

lightaero

aerodynamics

Aerodynamics discipline subpackage for lightaero.

Provides VLM and aeroelastic solvers. See: docs/theory/aerodynamics.md

AerostructDiscipline

Bases: DisciplineBase

Static aeroelastic discipline: VLM coupled to Euler-Bernoulli beam.

See: docs/theory/aerodynamics.md

Source code in src/lightaero/aerodynamics/aeroelastic.py

@register("aerostruct", "vlm_euler_beam")
class AerostructDiscipline(DisciplineBase):
    """Static aeroelastic discipline: VLM coupled to Euler-Bernoulli beam.

    See: docs/theory/aerodynamics.md
    """

    def __init__(self) -> None:
        self._vlm = VLMDiscipline()
        self._struct = StructuralDiscipline()

    def __call__(self, **inputs: Any) -> AerostructOutput:
        """Evaluate the coupled aerostructural discipline.

        Args:
            **inputs: Named physical quantities in SI units.
                wing (WingGeometry): Undeformed wing geometry.
                alpha_rad (float): Angle of attack in radians.
                altitude_m (float): Altitude in metres.
                MTOW_kg (float): Maximum take-off weight in kg.
                V_ms (float) | M (float): Freestream speed or Mach number (not both).

        Returns:
            AerostructOutput: Combined results including aero, struct, n_iter, and converged status.
        """
        wing: WingGeometry = inputs["wing"]
        MTOW_kg: float = float(inputs["MTOW_kg"])

        # Forward aero inputs to VLM (alpha, altitude, speed)
        aero_inputs = {k: inputs[k] for k in ("alpha_rad", "altitude_m") if k in inputs}
        if "V_ms" in inputs:
            aero_inputs["V_ms"] = inputs["V_ms"]
        elif "M" in inputs:
            aero_inputs["M"] = inputs["M"]

        # State at start of iteration
        deflection_z = np.zeros(wing.n_panels + 1, dtype=float)
        delta_twist = np.zeros_like(wing.twist)  # shape (n_stations,)

        # y_nodes used to map structural twist (at VLM nodes) to y_stations
        y_nodes = np.concatenate([[0.0], wing.y_cp])

        aero: AeroOutput | None = None
        struct: StructuralOutput | None = None
        converged = False

        for _n_iter in range(1, _MAX_ITER + 1):
            wing_def = _rebuild_wing(wing, deflection_z, delta_twist)

            aero = self._vlm(wing=wing_def, **aero_inputs)
            struct = self._struct(
                wing=wing_def,
                span_load=aero.span_load,
                span_drag=aero.span_drag,
                MTOW_kg=MTOW_kg,
            )

            new_deflection = np.asarray(struct.span_deflection, dtype=float)
            change = float(np.max(np.abs(new_deflection - deflection_z)))

            # Interpolate elastic twist from y_nodes to y_stations for geometry rebuild
            new_delta_twist = np.interp(wing.y_stations, y_nodes, np.asarray(struct.span_twist, dtype=float))

            deflection_z = new_deflection
            delta_twist = new_delta_twist

            if change < _CONV_TOL:
                converged = True
                break

        return AerostructOutput(
            aero=aero,
            struct=struct,
            n_iter=_n_iter,
            converged=converged,
        )

__call__(**inputs)

Evaluate the coupled aerostructural discipline.

Parameters:

Name	Type	Description	Default
**inputs	Any	Named physical quantities in SI units. wing (WingGeometry): Undeformed wing geometry. alpha_rad (float): Angle of attack in radians. altitude_m (float): Altitude in metres. MTOW_kg (float): Maximum take-off weight in kg. V_ms (float) | M (float): Freestream speed or Mach number (not both).	{}

Returns:

Name	Type	Description
AerostructOutput	AerostructOutput	Combined results including aero, struct, n_iter, and converged status.

Source code in src/lightaero/aerodynamics/aeroelastic.py

def __call__(self, **inputs: Any) -> AerostructOutput:
    """Evaluate the coupled aerostructural discipline.

    Args:
        **inputs: Named physical quantities in SI units.
            wing (WingGeometry): Undeformed wing geometry.
            alpha_rad (float): Angle of attack in radians.
            altitude_m (float): Altitude in metres.
            MTOW_kg (float): Maximum take-off weight in kg.
            V_ms (float) | M (float): Freestream speed or Mach number (not both).

    Returns:
        AerostructOutput: Combined results including aero, struct, n_iter, and converged status.
    """
    wing: WingGeometry = inputs["wing"]
    MTOW_kg: float = float(inputs["MTOW_kg"])

    # Forward aero inputs to VLM (alpha, altitude, speed)
    aero_inputs = {k: inputs[k] for k in ("alpha_rad", "altitude_m") if k in inputs}
    if "V_ms" in inputs:
        aero_inputs["V_ms"] = inputs["V_ms"]
    elif "M" in inputs:
        aero_inputs["M"] = inputs["M"]

    # State at start of iteration
    deflection_z = np.zeros(wing.n_panels + 1, dtype=float)
    delta_twist = np.zeros_like(wing.twist)  # shape (n_stations,)

    # y_nodes used to map structural twist (at VLM nodes) to y_stations
    y_nodes = np.concatenate([[0.0], wing.y_cp])

    aero: AeroOutput | None = None
    struct: StructuralOutput | None = None
    converged = False

    for _n_iter in range(1, _MAX_ITER + 1):
        wing_def = _rebuild_wing(wing, deflection_z, delta_twist)

        aero = self._vlm(wing=wing_def, **aero_inputs)
        struct = self._struct(
            wing=wing_def,
            span_load=aero.span_load,
            span_drag=aero.span_drag,
            MTOW_kg=MTOW_kg,
        )

        new_deflection = np.asarray(struct.span_deflection, dtype=float)
        change = float(np.max(np.abs(new_deflection - deflection_z)))

        # Interpolate elastic twist from y_nodes to y_stations for geometry rebuild
        new_delta_twist = np.interp(wing.y_stations, y_nodes, np.asarray(struct.span_twist, dtype=float))

        deflection_z = new_deflection
        delta_twist = new_delta_twist

        if change < _CONV_TOL:
            converged = True
            break

    return AerostructOutput(
        aero=aero,
        struct=struct,
        n_iter=_n_iter,
        converged=converged,
    )

AerostructOutput dataclass

Combined aerostructural result from AerostructDiscipline.

Attributes:

Name	Type	Description
aero	AeroOutput	Final AeroOutput (CL, CDi, CM, lift_N, drag_N, span_cl, span_load, span_drag) from the last VLM evaluation.
struct	StructuralOutput	Final StructuralOutput (tip_deflection_m, root_stress_Pa, wing_mass_kg, oew_kg, span_deflection, tip_twist_rad, span_twist) from the last structural evaluation.
n_iter	int	Number of coupling iterations performed.
converged	bool	True if convergence criterion was met before _MAX_ITER.

Source code in src/lightaero/aerodynamics/aeroelastic.py

@dataclass
class AerostructOutput:
    """Combined aerostructural result from AerostructDiscipline.

    Attributes:
        aero: Final AeroOutput (CL, CDi, CM, lift_N, drag_N, span_cl,
              span_load, span_drag) from the last VLM evaluation.
        struct: Final StructuralOutput (tip_deflection_m, root_stress_Pa,
                wing_mass_kg, oew_kg, span_deflection, tip_twist_rad,
                span_twist) from the last structural evaluation.
        n_iter: Number of coupling iterations performed.
        converged: True if convergence criterion was met before _MAX_ITER.
    """

    aero: AeroOutput
    struct: StructuralOutput
    n_iter: int
    converged: bool

VLMDiscipline

Bases: DisciplineBase

Low-fidelity VLM aerodynamics discipline.

See: docs/theory/aerodynamics.md

Note

Inputs: wing (WingGeometry), alpha_rad (float), V_ms (float), altitude_m (float). Returns: AeroOutput TypedDict. AIC matrix is cached based on geometry hash.

Source code in src/lightaero/aerodynamics/solver.py

@register("aero", "vlm")
class VLMDiscipline(DisciplineBase):
    """Low-fidelity VLM aerodynamics discipline.

    See: docs/theory/aerodynamics.md

    Note:
        Inputs: wing (WingGeometry), alpha_rad (float), V_ms (float), altitude_m (float).
        Returns: AeroOutput TypedDict.
        AIC matrix is cached based on geometry hash.
    """

    def __init__(self) -> None:
        self._aic_cache: np.ndarray | None = None
        self._aic_key: tuple | None = None
        self._results: AeroOutput | None = None

    def __call__(self, **inputs: Any) -> AeroOutput:
        """Evaluate VLM aerodynamics at the given flight condition.

        See: docs/theory/aerodynamics.md

        Args:
            **inputs (Any): wing, alpha_rad, V_ms (or M), altitude_m.

        Returns:
            AeroOutput: Validated results.
        """

        wing: WingGeometry = inputs["wing"]
        altitude_m: float = float(inputs["altitude_m"])
        rho, T_K, _, a = isa(altitude_m)

        if inputs.get("M", None) is not None and inputs.get("V_ms", None) is not None:
            raise ValueError("VLMDiscipline only accepts M or V_ms but not both")
        elif inputs.get("M", None) is None and inputs.get("V_ms", None) is None:
            raise ValueError("VLMDiscipline requires either M or V_ms as input.")
        elif inputs.get("M", None) is None and inputs.get("V_ms", None) is not None:
            V_ms: float = float(inputs["V_ms"])
            M = V_ms / a
        else:
            M: float = float(inputs["M"])
            V_ms = M * a

        if M >= 1.0:
            raise ValueError(f"Mach number {M:.2f} >= 1.0. Prandtl-Glauert correction is valid only for subsonic flow.")

        # Prandtl-Glauert compressibility factor
        beta = math.sqrt(1.0 - M**2)

        mesh = build_panel_mesh(wing)

        # Handle spanwise-varying alpha_rad (scalar, station-wise, or panel-wise)
        alpha_in = inputs["alpha_rad"]
        if isinstance(alpha_in, (list, np.ndarray)):
            alpha_in = np.asarray(alpha_in, dtype=np.float64)
            if len(alpha_in) == len(wing.y_stations):
                # Interpolate from stations to right semi-span control points
                alpha_cp_right = np.interp(wing.y_cp, wing.y_stations, alpha_in)
            elif len(alpha_in) == wing.n_panels:
                alpha_cp_right = alpha_in
            else:
                raise ValueError(
                    f"alpha_rad array length {len(alpha_in)} must match "
                    f"n_stations ({len(wing.y_stations)}) or n_panels ({wing.n_panels})"
                )
            # Full-span mirror: left semi-span (-tip -> root) + right semi-span (root -> tip)
            alpha_rad = np.concatenate([alpha_cp_right[::-1], alpha_cp_right])
        else:
            alpha_rad = float(alpha_in)

        # Safety Guards: Check for validity regime (Mach < 0.3, AoA < 10 deg)
        check_regime_validity(M, float(np.max(np.abs(alpha_rad))))

        epsilon = 1e-6 * 2.0 * wing.half_span

        cache_key = (
            wing.chord.tobytes(),
            wing.sweep_le.tobytes(),
            wing.y_stations.tobytes(),
            wing.n_panels,
            wing.half_span,
        )
        if self._aic_key != cache_key or self._aic_cache is None:
            self._aic_cache = assemble_aic(mesh, epsilon=epsilon)
            self._aic_key = cache_key
        AIC = self._aic_cache

        RHS = build_rhs(mesh, alpha_rad=alpha_rad, V_ms=V_ms)
        gamma = np.linalg.solve(AIC, RHS)

        CL, CDi_kj, CM, lift_N, drag_i_N, span_cl, span_load, span_drag = near_field_forces(
            mesh,
            gamma,
            rho,
            V_ms,
            mesh.S_ref,
            alpha_rad,
        )
        CDi_tp = trefftz_cdi(mesh, gamma, V_ms, mesh.S_ref)

        # Panel-local cos-sweep Prandtl-Glauert correction (A3).
        # Each strip sees an effective Mach M_eff = M * cos(local LE sweep).
        # This replaces the global beta = sqrt(1 - M^2) for span quantities.
        sweep_cp = np.interp(wing.y_cp, wing.y_stations, wing.sweep_le)
        M_eff_cp = np.clip(M * np.cos(sweep_cp), 0.0, 0.99)
        beta_cp = np.sqrt(1.0 - M_eff_cp**2)
        beta_cp = np.maximum(beta_cp, 0.01)  # numerical floor

        # Full-span mirror: left semi-span (-tip -> root) + right semi-span (root -> tip)
        beta_full = np.concatenate([beta_cp[::-1], beta_cp])

        span_cl = span_cl / beta_full
        span_load = span_load / beta_full
        span_drag = span_drag / beta_full

        # Reintegrate CL and lift_N from corrected spanwise distributions (right semi-span)
        q_inf = 0.5 * rho * V_ms**2
        CL = 2.0 * float(np.trapezoid(span_cl[wing.n_panels :] * wing.chord_cp, wing.y_cp)) / mesh.S_ref
        lift_N = CL * q_inf * mesh.S_ref

        # CDi, CM, drag_N: global Prandtl-Glauert (less sensitive to span variation)
        CDi_tp /= beta
        CM /= beta

        # Profile drag via strip theory (right semi-span only; symmetric)
        tc_stations = np.array([float(sec.thickness_at(np.array([0.35]))[0]) for sec in wing.airfoil_sections])
        tc_cp = np.interp(wing.y_cp, wing.y_stations, tc_stations)

        CD0 = profile_drag_cd0(
            rho=rho,
            V_ms=V_ms,
            T_K=T_K,
            S_ref=mesh.S_ref,
            y_cp_half=wing.y_cp,
            chord_cp_half=wing.chord_cp,
            tc_cp_half=tc_cp,
        )
        CD_total = float(CDi_tp) + CD0

        # Total Drag Force: Induced (corrected) + Viscous Profile
        drag_i_total_N = CDi_tp * q_inf * mesh.S_ref
        drag_v_total_N = CD0 * q_inf * mesh.S_ref
        drag_total_N = drag_i_total_N + drag_v_total_N

        self._results: AeroOutput = AeroOutput(
            CL=float(CL),
            CDi=float(CDi_tp),
            CD0=CD0,
            CD=CD_total,
            CM=float(CM),
            lift_N=float(lift_N),
            drag_N=float(drag_total_N),
            span_cl=span_cl,
            span_load=span_load,
            span_drag=span_drag,
            circulation=gamma,
        )

        if inputs.get("validate_output", True):
            validate_discipline_output(self._results, AERO_OUTPUT_SPEC)
        return self._results

    def __str__(self) -> str:
        if self._results is None:
            return "VLMDiscipline: No results yet. Call the discipline with valid inputs to compute aerodynamics."
        return "\n".join(
            [
                f"Lift Coefficient           (CL) : {self._results.CL:.4e}",
                f"Drag Coefficient           (CD) : {self._results.CD:.4e}",
                f"Induced Drag Coefficient   (CDi): {self._results.CDi:.4e}",
                f"Zero-Lift Drag Coefficient (CD0): {self._results.CD0:.4e}",
                f"Coefficient of Moment      (CM) : {self._results.CM:.4e}",
            ]
        )

__call__(**inputs)

Evaluate VLM aerodynamics at the given flight condition.

See: docs/theory/aerodynamics.md

Parameters:

Name	Type	Description	Default
**inputs	Any	wing, alpha_rad, V_ms (or M), altitude_m.	{}

Returns:

Name	Type	Description
AeroOutput	AeroOutput	Validated results.

Source code in src/lightaero/aerodynamics/solver.py

def __call__(self, **inputs: Any) -> AeroOutput:
    """Evaluate VLM aerodynamics at the given flight condition.

    See: docs/theory/aerodynamics.md

    Args:
        **inputs (Any): wing, alpha_rad, V_ms (or M), altitude_m.

    Returns:
        AeroOutput: Validated results.
    """

    wing: WingGeometry = inputs["wing"]
    altitude_m: float = float(inputs["altitude_m"])
    rho, T_K, _, a = isa(altitude_m)

    if inputs.get("M", None) is not None and inputs.get("V_ms", None) is not None:
        raise ValueError("VLMDiscipline only accepts M or V_ms but not both")
    elif inputs.get("M", None) is None and inputs.get("V_ms", None) is None:
        raise ValueError("VLMDiscipline requires either M or V_ms as input.")
    elif inputs.get("M", None) is None and inputs.get("V_ms", None) is not None:
        V_ms: float = float(inputs["V_ms"])
        M = V_ms / a
    else:
        M: float = float(inputs["M"])
        V_ms = M * a

    if M >= 1.0:
        raise ValueError(f"Mach number {M:.2f} >= 1.0. Prandtl-Glauert correction is valid only for subsonic flow.")

    # Prandtl-Glauert compressibility factor
    beta = math.sqrt(1.0 - M**2)

    mesh = build_panel_mesh(wing)

    # Handle spanwise-varying alpha_rad (scalar, station-wise, or panel-wise)
    alpha_in = inputs["alpha_rad"]
    if isinstance(alpha_in, (list, np.ndarray)):
        alpha_in = np.asarray(alpha_in, dtype=np.float64)
        if len(alpha_in) == len(wing.y_stations):
            # Interpolate from stations to right semi-span control points
            alpha_cp_right = np.interp(wing.y_cp, wing.y_stations, alpha_in)
        elif len(alpha_in) == wing.n_panels:
            alpha_cp_right = alpha_in
        else:
            raise ValueError(
                f"alpha_rad array length {len(alpha_in)} must match "
                f"n_stations ({len(wing.y_stations)}) or n_panels ({wing.n_panels})"
            )
        # Full-span mirror: left semi-span (-tip -> root) + right semi-span (root -> tip)
        alpha_rad = np.concatenate([alpha_cp_right[::-1], alpha_cp_right])
    else:
        alpha_rad = float(alpha_in)

    # Safety Guards: Check for validity regime (Mach < 0.3, AoA < 10 deg)
    check_regime_validity(M, float(np.max(np.abs(alpha_rad))))

    epsilon = 1e-6 * 2.0 * wing.half_span

    cache_key = (
        wing.chord.tobytes(),
        wing.sweep_le.tobytes(),
        wing.y_stations.tobytes(),
        wing.n_panels,
        wing.half_span,
    )
    if self._aic_key != cache_key or self._aic_cache is None:
        self._aic_cache = assemble_aic(mesh, epsilon=epsilon)
        self._aic_key = cache_key
    AIC = self._aic_cache

    RHS = build_rhs(mesh, alpha_rad=alpha_rad, V_ms=V_ms)
    gamma = np.linalg.solve(AIC, RHS)

    CL, CDi_kj, CM, lift_N, drag_i_N, span_cl, span_load, span_drag = near_field_forces(
        mesh,
        gamma,
        rho,
        V_ms,
        mesh.S_ref,
        alpha_rad,
    )
    CDi_tp = trefftz_cdi(mesh, gamma, V_ms, mesh.S_ref)

    # Panel-local cos-sweep Prandtl-Glauert correction (A3).
    # Each strip sees an effective Mach M_eff = M * cos(local LE sweep).
    # This replaces the global beta = sqrt(1 - M^2) for span quantities.
    sweep_cp = np.interp(wing.y_cp, wing.y_stations, wing.sweep_le)
    M_eff_cp = np.clip(M * np.cos(sweep_cp), 0.0, 0.99)
    beta_cp = np.sqrt(1.0 - M_eff_cp**2)
    beta_cp = np.maximum(beta_cp, 0.01)  # numerical floor

    # Full-span mirror: left semi-span (-tip -> root) + right semi-span (root -> tip)
    beta_full = np.concatenate([beta_cp[::-1], beta_cp])

    span_cl = span_cl / beta_full
    span_load = span_load / beta_full
    span_drag = span_drag / beta_full

    # Reintegrate CL and lift_N from corrected spanwise distributions (right semi-span)
    q_inf = 0.5 * rho * V_ms**2
    CL = 2.0 * float(np.trapezoid(span_cl[wing.n_panels :] * wing.chord_cp, wing.y_cp)) / mesh.S_ref
    lift_N = CL * q_inf * mesh.S_ref

    # CDi, CM, drag_N: global Prandtl-Glauert (less sensitive to span variation)
    CDi_tp /= beta
    CM /= beta

    # Profile drag via strip theory (right semi-span only; symmetric)
    tc_stations = np.array([float(sec.thickness_at(np.array([0.35]))[0]) for sec in wing.airfoil_sections])
    tc_cp = np.interp(wing.y_cp, wing.y_stations, tc_stations)

    CD0 = profile_drag_cd0(
        rho=rho,
        V_ms=V_ms,
        T_K=T_K,
        S_ref=mesh.S_ref,
        y_cp_half=wing.y_cp,
        chord_cp_half=wing.chord_cp,
        tc_cp_half=tc_cp,
    )
    CD_total = float(CDi_tp) + CD0

    # Total Drag Force: Induced (corrected) + Viscous Profile
    drag_i_total_N = CDi_tp * q_inf * mesh.S_ref
    drag_v_total_N = CD0 * q_inf * mesh.S_ref
    drag_total_N = drag_i_total_N + drag_v_total_N

    self._results: AeroOutput = AeroOutput(
        CL=float(CL),
        CDi=float(CDi_tp),
        CD0=CD0,
        CD=CD_total,
        CM=float(CM),
        lift_N=float(lift_N),
        drag_N=float(drag_total_N),
        span_cl=span_cl,
        span_load=span_load,
        span_drag=span_drag,
        circulation=gamma,
    )

    if inputs.get("validate_output", True):
        validate_discipline_output(self._results, AERO_OUTPUT_SPEC)
    return self._results

aeroelastic

Aeroelastic discipline: coupled VLM + Euler-Bernoulli beam iteration.

Iterates VLM and structural solvers to predict spanwise deflection and twist. See: docs/theory/aerodynamics.md

AerostructDiscipline

Bases: DisciplineBase

Static aeroelastic discipline: VLM coupled to Euler-Bernoulli beam.

See: docs/theory/aerodynamics.md

Source code in src/lightaero/aerodynamics/aeroelastic.py

@register("aerostruct", "vlm_euler_beam")
class AerostructDiscipline(DisciplineBase):
    """Static aeroelastic discipline: VLM coupled to Euler-Bernoulli beam.

    See: docs/theory/aerodynamics.md
    """

    def __init__(self) -> None:
        self._vlm = VLMDiscipline()
        self._struct = StructuralDiscipline()

    def __call__(self, **inputs: Any) -> AerostructOutput:
        """Evaluate the coupled aerostructural discipline.

        Args:
            **inputs: Named physical quantities in SI units.
                wing (WingGeometry): Undeformed wing geometry.
                alpha_rad (float): Angle of attack in radians.
                altitude_m (float): Altitude in metres.
                MTOW_kg (float): Maximum take-off weight in kg.
                V_ms (float) | M (float): Freestream speed or Mach number (not both).

        Returns:
            AerostructOutput: Combined results including aero, struct, n_iter, and converged status.
        """
        wing: WingGeometry = inputs["wing"]
        MTOW_kg: float = float(inputs["MTOW_kg"])

        # Forward aero inputs to VLM (alpha, altitude, speed)
        aero_inputs = {k: inputs[k] for k in ("alpha_rad", "altitude_m") if k in inputs}
        if "V_ms" in inputs:
            aero_inputs["V_ms"] = inputs["V_ms"]
        elif "M" in inputs:
            aero_inputs["M"] = inputs["M"]

        # State at start of iteration
        deflection_z = np.zeros(wing.n_panels + 1, dtype=float)
        delta_twist = np.zeros_like(wing.twist)  # shape (n_stations,)

        # y_nodes used to map structural twist (at VLM nodes) to y_stations
        y_nodes = np.concatenate([[0.0], wing.y_cp])

        aero: AeroOutput | None = None
        struct: StructuralOutput | None = None
        converged = False

        for _n_iter in range(1, _MAX_ITER + 1):
            wing_def = _rebuild_wing(wing, deflection_z, delta_twist)

            aero = self._vlm(wing=wing_def, **aero_inputs)
            struct = self._struct(
                wing=wing_def,
                span_load=aero.span_load,
                span_drag=aero.span_drag,
                MTOW_kg=MTOW_kg,
            )

            new_deflection = np.asarray(struct.span_deflection, dtype=float)
            change = float(np.max(np.abs(new_deflection - deflection_z)))

            # Interpolate elastic twist from y_nodes to y_stations for geometry rebuild
            new_delta_twist = np.interp(wing.y_stations, y_nodes, np.asarray(struct.span_twist, dtype=float))

            deflection_z = new_deflection
            delta_twist = new_delta_twist

            if change < _CONV_TOL:
                converged = True
                break

        return AerostructOutput(
            aero=aero,
            struct=struct,
            n_iter=_n_iter,
            converged=converged,
        )

__call__(**inputs)

Evaluate the coupled aerostructural discipline.

Parameters:

Name	Type	Description	Default
**inputs	Any	Named physical quantities in SI units. wing (WingGeometry): Undeformed wing geometry. alpha_rad (float): Angle of attack in radians. altitude_m (float): Altitude in metres. MTOW_kg (float): Maximum take-off weight in kg. V_ms (float) | M (float): Freestream speed or Mach number (not both).	{}

Returns:

Name	Type	Description
AerostructOutput	AerostructOutput	Combined results including aero, struct, n_iter, and converged status.

Source code in src/lightaero/aerodynamics/aeroelastic.py

def __call__(self, **inputs: Any) -> AerostructOutput:
    """Evaluate the coupled aerostructural discipline.

    Args:
        **inputs: Named physical quantities in SI units.
            wing (WingGeometry): Undeformed wing geometry.
            alpha_rad (float): Angle of attack in radians.
            altitude_m (float): Altitude in metres.
            MTOW_kg (float): Maximum take-off weight in kg.
            V_ms (float) | M (float): Freestream speed or Mach number (not both).

    Returns:
        AerostructOutput: Combined results including aero, struct, n_iter, and converged status.
    """
    wing: WingGeometry = inputs["wing"]
    MTOW_kg: float = float(inputs["MTOW_kg"])

    # Forward aero inputs to VLM (alpha, altitude, speed)
    aero_inputs = {k: inputs[k] for k in ("alpha_rad", "altitude_m") if k in inputs}
    if "V_ms" in inputs:
        aero_inputs["V_ms"] = inputs["V_ms"]
    elif "M" in inputs:
        aero_inputs["M"] = inputs["M"]

    # State at start of iteration
    deflection_z = np.zeros(wing.n_panels + 1, dtype=float)
    delta_twist = np.zeros_like(wing.twist)  # shape (n_stations,)

    # y_nodes used to map structural twist (at VLM nodes) to y_stations
    y_nodes = np.concatenate([[0.0], wing.y_cp])

    aero: AeroOutput | None = None
    struct: StructuralOutput | None = None
    converged = False

    for _n_iter in range(1, _MAX_ITER + 1):
        wing_def = _rebuild_wing(wing, deflection_z, delta_twist)

        aero = self._vlm(wing=wing_def, **aero_inputs)
        struct = self._struct(
            wing=wing_def,
            span_load=aero.span_load,
            span_drag=aero.span_drag,
            MTOW_kg=MTOW_kg,
        )

        new_deflection = np.asarray(struct.span_deflection, dtype=float)
        change = float(np.max(np.abs(new_deflection - deflection_z)))

        # Interpolate elastic twist from y_nodes to y_stations for geometry rebuild
        new_delta_twist = np.interp(wing.y_stations, y_nodes, np.asarray(struct.span_twist, dtype=float))

        deflection_z = new_deflection
        delta_twist = new_delta_twist

        if change < _CONV_TOL:
            converged = True
            break

    return AerostructOutput(
        aero=aero,
        struct=struct,
        n_iter=_n_iter,
        converged=converged,
    )

AerostructOutput dataclass

Combined aerostructural result from AerostructDiscipline.

Attributes:

Name	Type	Description
aero	AeroOutput	Final AeroOutput (CL, CDi, CM, lift_N, drag_N, span_cl, span_load, span_drag) from the last VLM evaluation.
struct	StructuralOutput	Final StructuralOutput (tip_deflection_m, root_stress_Pa, wing_mass_kg, oew_kg, span_deflection, tip_twist_rad, span_twist) from the last structural evaluation.
n_iter	int	Number of coupling iterations performed.
converged	bool	True if convergence criterion was met before _MAX_ITER.

Source code in src/lightaero/aerodynamics/aeroelastic.py

@dataclass
class AerostructOutput:
    """Combined aerostructural result from AerostructDiscipline.

    Attributes:
        aero: Final AeroOutput (CL, CDi, CM, lift_N, drag_N, span_cl,
              span_load, span_drag) from the last VLM evaluation.
        struct: Final StructuralOutput (tip_deflection_m, root_stress_Pa,
                wing_mass_kg, oew_kg, span_deflection, tip_twist_rad,
                span_twist) from the last structural evaluation.
        n_iter: Number of coupling iterations performed.
        converged: True if convergence criterion was met before _MAX_ITER.
    """

    aero: AeroOutput
    struct: StructuralOutput
    n_iter: int
    converged: bool

aic

AIC matrix assembly for VLM with Rankine core regularisation.

See: docs/theory/aerodynamics.md

assemble_aic(mesh, epsilon)

Assemble the AIC matrix for a PanelMesh.

See: docs/theory/aerodynamics.md

Parameters:

Name	Type	Description	Default
mesh	PanelMesh	PanelMesh object.	required
epsilon	float	Rankine core fraction.	required

Returns:

Name	Type	Description
AIC	ndarray	Aerodynamic Influence Coefficient matrix.

Raises:

Type	Description
AssertionError	If condition number of AIC exceeds 1e8.

Source code in src/lightaero/aerodynamics/aic.py

def assemble_aic(mesh: PanelMesh, epsilon: float) -> np.ndarray:
    """Assemble the AIC matrix for a PanelMesh.

    See: docs/theory/aerodynamics.md

    Args:
        mesh: PanelMesh object.
        epsilon: Rankine core fraction.

    Returns:
        AIC: Aerodynamic Influence Coefficient matrix.

    Raises:
        AssertionError: If condition number of AIC exceeds 1e8.
    """
    AIC = assemble_aic_kernel(
        mesh.xA.astype(np.float64),
        mesh.yA.astype(np.float64),
        mesh.zA.astype(np.float64),
        mesh.xB.astype(np.float64),
        mesh.yB.astype(np.float64),
        mesh.zB.astype(np.float64),
        mesh.xC.astype(np.float64),
        mesh.yC.astype(np.float64),
        mesh.zC.astype(np.float64),
        mesh.nx.astype(np.float64),
        mesh.ny.astype(np.float64),
        mesh.nz.astype(np.float64),
        float(epsilon),
        1.0,
        0.0,
        0.0,  # freestream unit direction
    )

    cond = np.linalg.cond(AIC)
    assert cond < 1e8, (
        f"AIC condition number {cond:.3e} >= 1e8: mesh may be degenerate (n_panels={mesh.n_panels}, epsilon={epsilon})"
    )

    return AIC

assemble_aic_kernel(xA, yA, zA, xB, yB, zB, xC, yC, zC, nx, ny, nz, epsilon, u_inf_x, u_inf_y, u_inf_z)

Assemble AIC matrix kernel via explicit double loop.

Parameters:

Name	Type	Description	Default
xA	ndarray	Inboard bound vortex x-endpoints, shape (n,).	required
yA	ndarray	Inboard bound vortex y-endpoints, shape (n,).	required
zA	ndarray	Inboard bound vortex z-endpoints, shape (n,).	required
xB	ndarray	Outboard bound vortex x-endpoints, shape (n,).	required
yB	ndarray	Outboard bound vortex y-endpoints, shape (n,).	required
zB	ndarray	Outboard bound vortex z-endpoints, shape (n,).	required
xC	ndarray	Control point x-coordinates, shape (n,).	required
yC	ndarray	Control point y-coordinates, shape (n,).	required
zC	ndarray	Control point z-coordinates, shape (n,).	required
nx	ndarray	Unit normal x-components, shape (n,).	required
ny	ndarray	Unit normal y-components, shape (n,).	required
nz	ndarray	Unit normal z-components, shape (n,).	required
epsilon	float	Rankine core fraction.	required
u_inf_x	float	Freestream unit direction x-component.	required
u_inf_y	float	Freestream unit direction y-component.	required
u_inf_z	float	Freestream unit direction z-component.	required

Returns:

Type	Description
ndarray	np.ndarray: AIC matrix of shape (n, n).

Source code in src/lightaero/aerodynamics/aic.py

@numba.njit(cache=True)
def assemble_aic_kernel(
    xA: np.ndarray,
    yA: np.ndarray,
    zA: np.ndarray,
    xB: np.ndarray,
    yB: np.ndarray,
    zB: np.ndarray,
    xC: np.ndarray,
    yC: np.ndarray,
    zC: np.ndarray,
    nx: np.ndarray,
    ny: np.ndarray,
    nz: np.ndarray,
    epsilon: float,
    u_inf_x: float,
    u_inf_y: float,
    u_inf_z: float,
) -> np.ndarray:
    """Assemble AIC matrix kernel via explicit double loop.

    Args:
        xA (np.ndarray): Inboard bound vortex x-endpoints, shape (n,).
        yA (np.ndarray): Inboard bound vortex y-endpoints, shape (n,).
        zA (np.ndarray): Inboard bound vortex z-endpoints, shape (n,).
        xB (np.ndarray): Outboard bound vortex x-endpoints, shape (n,).
        yB (np.ndarray): Outboard bound vortex y-endpoints, shape (n,).
        zB (np.ndarray): Outboard bound vortex z-endpoints, shape (n,).
        xC (np.ndarray): Control point x-coordinates, shape (n,).
        yC (np.ndarray): Control point y-coordinates, shape (n,).
        zC (np.ndarray): Control point z-coordinates, shape (n,).
        nx (np.ndarray): Unit normal x-components, shape (n,).
        ny (np.ndarray): Unit normal y-components, shape (n,).
        nz (np.ndarray): Unit normal z-components, shape (n,).
        epsilon (float): Rankine core fraction.
        u_inf_x (float): Freestream unit direction x-component.
        u_inf_y (float): Freestream unit direction y-component.
        u_inf_z (float): Freestream unit direction z-component.

    Returns:
        np.ndarray: AIC matrix of shape (n, n).
    """
    n = xA.shape[0]
    AIC = np.zeros((n, n))

    for i in range(n):
        Px = xC[i]
        Py = yC[i]
        Pz = zC[i]
        nxi = nx[i]
        nyi = ny[i]
        nzi = nz[i]

        for j in range(n):
            # Bound segment A_j -> B_j, gamma = +1
            ux_b, uy_b, uz_b = _biot_savart_finite(
                xA[j],
                yA[j],
                zA[j],
                xB[j],
                yB[j],
                zB[j],
                Px,
                Py,
                Pz,
                1.0,
                epsilon,
            )

            # Right trailing leg: B_j -> infinity, gamma = +1
            ux_r, uy_r, uz_r = _biot_savart_semi_infinite(
                xB[j],
                yB[j],
                zB[j],
                u_inf_x,
                u_inf_y,
                u_inf_z,
                Px,
                Py,
                Pz,
                1.0,
                epsilon,
            )

            # Left trailing leg: A_j -> infinity, gamma = -1
            ux_l, uy_l, uz_l = _biot_savart_semi_infinite(
                xA[j],
                yA[j],
                zA[j],
                u_inf_x,
                u_inf_y,
                u_inf_z,
                Px,
                Py,
                Pz,
                -1.0,
                epsilon,
            )

            # Sum all contributions and project onto control-point normal
            ux = ux_b + ux_r + ux_l
            uy = uy_b + uy_r + uy_l
            uz = uz_b + uz_r + uz_l

            AIC[i, j] = ux * nxi + uy * nyi + uz * nzi

    return AIC

build_rhs(mesh, alpha_rad, V_ms=1.0)

Build the right-hand side vector for the VLM linear system.

See: docs/theory/aerodynamics.md

Parameters:

Name	Type	Description	Default
mesh	PanelMesh	PanelMesh with unit normals.	required
alpha_rad	float | ndarray	Angle of attack (radians).	required
V_ms	float	Freestream speed (m/s).	1.0

Returns:

Name	Type	Description
RHS	ndarray	shape (n_panels,).

Source code in src/lightaero/aerodynamics/aic.py

def build_rhs(mesh: PanelMesh, alpha_rad: float | np.ndarray, V_ms: float = 1.0) -> np.ndarray:
    """Build the right-hand side vector for the VLM linear system.

    See: docs/theory/aerodynamics.md

    Args:
        mesh: PanelMesh with unit normals.
        alpha_rad: Angle of attack (radians).
        V_ms: Freestream speed (m/s).

    Returns:
        RHS: shape (n_panels,).
    """
    if isinstance(alpha_rad, (np.ndarray, list)):
        alpha_rad = np.asarray(alpha_rad, dtype=np.float64)
        vx = V_ms * np.cos(alpha_rad)
        vy = np.zeros_like(alpha_rad)
        vz = V_ms * np.sin(alpha_rad)
    else:
        vx = V_ms * math.cos(alpha_rad)
        vy = 0.0
        vz = V_ms * math.sin(alpha_rad)

    # RHS[i] = -(V_inf . n_hat_i)
    RHS = -(vx * mesh.nx + vy * mesh.ny + vz * mesh.nz)
    return RHS.astype(np.float64)

force

Aerodynamic force integration for VLM.

See: docs/theory/aerodynamics.md

near_field_forces(mesh, gamma, rho, V_ms, S_ref, alpha_rad)

Compute near-field aerodynamic forces and coefficients.

See: docs/theory/aerodynamics.md

Parameters:

Name	Type	Description	Default
mesh	PanelMesh	Panel mesh object.	required
gamma	ndarray	Circulation strengths.	required
rho	float	Air density.	required
V_ms	float	Freestream speed.	required
S_ref	float	Reference area.	required
alpha_rad	float | ndarray	Angle of attack in radians.	required

Returns:

Name	Type	Description
tuple	tuple	(CL, CDi_kj, CM, lift_N, drag_N, span_cl, span_load, span_drag)

Source code in src/lightaero/aerodynamics/force.py

def near_field_forces(
    mesh: PanelMesh,
    gamma: np.ndarray,
    rho: float,
    V_ms: float,
    S_ref: float,
    alpha_rad: float | np.ndarray,
) -> tuple:
    """Compute near-field aerodynamic forces and coefficients.

    See: docs/theory/aerodynamics.md

    Args:
        mesh: Panel mesh object.
        gamma: Circulation strengths.
        rho: Air density.
        V_ms: Freestream speed.
        S_ref: Reference area.
        alpha_rad: Angle of attack in radians.

    Returns:
        tuple: (CL, CDi_kj, CM, lift_N, drag_N, span_cl, span_load, span_drag)
    """
    n = mesh.n_panels
    q_inf = 0.5 * rho * V_ms**2

    if isinstance(alpha_rad, np.ndarray):
        # V_inf = V_ms * [cos(alpha), 0, sin(alpha)] for each panel
        V_inf = np.stack([V_ms * np.cos(alpha_rad), np.zeros(n), V_ms * np.sin(alpha_rad)], axis=1)
        # lift_dir = [-sin(alpha), 0, cos(alpha)] for each panel
        lift_dir = np.stack([-np.sin(alpha_rad), np.zeros(n), np.cos(alpha_rad)], axis=1)
        # drag_dir = [cos(alpha), 0, sin(alpha)] for each panel
        drag_dir = np.stack([np.cos(alpha_rad), np.zeros(n), np.sin(alpha_rad)], axis=1)
    else:
        V_inf = V_ms * np.array([np.cos(alpha_rad), 0.0, np.sin(alpha_rad)])
        lift_dir = np.array([-np.sin(alpha_rad), 0.0, np.cos(alpha_rad)])
        drag_dir = np.array([np.cos(alpha_rad), 0.0, np.sin(alpha_rad)])

    epsilon = 1e-6 * 2.0 * float(np.max(np.abs(mesh.y_cp)))

    bound_mid_x = 0.5 * (mesh.xA + mesh.xB)
    bound_mid_y = 0.5 * (mesh.yA + mesh.yB)
    bound_mid_z = 0.5 * (mesh.zA + mesh.zB)

    v_ind = _induced_velocity_kernel(
        mesh.xA,
        mesh.yA,
        mesh.zA,
        mesh.xB,
        mesh.yB,
        mesh.zB,
        bound_mid_x,
        bound_mid_y,
        bound_mid_z,
        gamma,
        1.0,
        0.0,
        0.0,
        float(epsilon),
    )

    if isinstance(alpha_rad, np.ndarray):
        V_local = V_inf + v_ind
    else:
        V_local = V_inf[np.newaxis, :] + v_ind

    dl = np.stack(
        [
            mesh.xB - mesh.xA,
            mesh.yB - mesh.yA,
            mesh.zB - mesh.zA,
        ],
        axis=1,
    )

    F = rho * gamma[:, np.newaxis] * np.cross(V_local, dl)

    if isinstance(alpha_rad, np.ndarray):
        # Row-wise dot product: sum(F * dir, axis=1)
        lift_j = np.sum(F * lift_dir, axis=1)
        drag_j = np.sum(F * drag_dir, axis=1)
    else:
        lift_j = F @ lift_dir
        drag_j = F @ drag_dir

    lift_N = float(np.sum(lift_j))
    drag_N = float(np.sum(drag_j))

    CL = lift_N / (q_inf * S_ref)
    CDi_kj = drag_N / (q_inf * S_ref)

    span_cl = 2.0 * gamma / (V_ms * mesh.chord_cp)
    safe_dy = np.where(
        np.abs(mesh.yB - mesh.yA) > 1e-14,
        np.abs(mesh.yB - mesh.yA),
        1e-14,
    )

    span_load = lift_j / safe_dy
    span_drag = drag_j / safe_dy

    n_half = n // 2
    # Use right semi-span for MAC and reference point calculations (indices n_half:)
    y_right = mesh.y_cp[n_half:]
    chord_right = mesh.chord_cp[n_half:]

    MAC = (2.0 / S_ref) * float(np.trapezoid(chord_right**2, y_right))
    x_le_panel = mesh.xC[n_half:] - 0.75 * chord_right
    x_le_mac = float((2.0 / S_ref) * np.trapezoid(x_le_panel * chord_right, y_right))
    x_ref = x_le_mac + 0.25 * MAC

    moment_arm = mesh.xC - x_ref
    CM = float(-np.sum(lift_j * moment_arm) / (q_inf * S_ref * MAC))

    return CL, CDi_kj, CM, lift_N, drag_N, span_cl, span_load, span_drag

profile_drag_cd0(rho, V_ms, T_K, S_ref, y_cp_half, chord_cp_half, tc_cp_half)

Profile drag coefficient via turbulent flat-plate strip theory.

See: docs/theory/aerodynamics.md

Parameters:

Name	Type	Description	Default
rho	float	Air density (kg/m³).	required
V_ms	float	Freestream speed (m/s).	required
T_K	float	Static temperature (K) for Sutherland viscosity.	required
S_ref	float	Reference wing area (m²).	required
y_cp_half	ndarray	Spanwise CP positions for the right semi-span, shape (n_half,).	required
chord_cp_half	ndarray	Chord at each CP, shape (n_half,).	required
tc_cp_half	ndarray	Thickness-to-chord ratio at each CP, shape (n_half,).	required

Returns:

Name	Type	Description
CD0	float	Total zero-lift profile drag coefficient (dimensionless).

Source code in src/lightaero/aerodynamics/force.py

def profile_drag_cd0(
    rho: float,
    V_ms: float,
    T_K: float,
    S_ref: float,
    y_cp_half: np.ndarray,
    chord_cp_half: np.ndarray,
    tc_cp_half: np.ndarray,
) -> float:
    """Profile drag coefficient via turbulent flat-plate strip theory.

    See: docs/theory/aerodynamics.md

    Args:
        rho: Air density (kg/m³).
        V_ms: Freestream speed (m/s).
        T_K: Static temperature (K) for Sutherland viscosity.
        S_ref: Reference wing area (m²).
        y_cp_half: Spanwise CP positions for the right semi-span, shape (n_half,).
        chord_cp_half: Chord at each CP, shape (n_half,).
        tc_cp_half: Thickness-to-chord ratio at each CP, shape (n_half,).

    Returns:
        CD0: Total zero-lift profile drag coefficient (dimensionless).
    """
    # Sutherland's law for dynamic viscosity
    _mu_ref = 1.716e-5  # Pa·s at T_ref
    _T_ref = 273.15  # K
    _S = 110.4  # K  (Sutherland constant for air)
    mu = _mu_ref * (T_K / _T_ref) ** 1.5 * (_T_ref + _S) / (T_K + _S)

    Re_c = rho * V_ms * chord_cp_half / mu
    Re_c = np.maximum(Re_c, 1e4)  # floor to avoid division instability at near-zero V

    Cf = 0.074 / Re_c**0.2
    FF = 1.0 + 2.0 * tc_cp_half
    cd0_strip = 2.0 * Cf * FF  # both surfaces

    # Integrate dD0 = cd0 * c * dy over right semi-span, multiply by 2 for both sides
    q_inf = 0.5 * rho * V_ms**2
    dD0 = q_inf * cd0_strip * chord_cp_half  # force per unit span
    D0_semi = float(np.trapezoid(dD0, y_cp_half))
    CD0 = 2.0 * D0_semi / (q_inf * S_ref)
    return CD0

trefftz_cdi(mesh, gamma, V_ms, S_ref)

Compute induced drag in the Trefftz plane.

See: docs/theory/aerodynamics.md

Parameters:

Name	Type	Description	Default
mesh	PanelMesh	Panel mesh object.	required
gamma	ndarray	Circulation strengths.	required
V_ms	float	Freestream speed.	required
S_ref	float	Reference area.	required

Returns:

Name	Type	Description
CDi_tp	float	Induced drag coefficient.

Source code in src/lightaero/aerodynamics/force.py

def trefftz_cdi(
    mesh: PanelMesh,
    gamma: np.ndarray,
    V_ms: float,
    S_ref: float,
) -> float:
    """Compute induced drag in the Trefftz plane.

    See: docs/theory/aerodynamics.md

    Args:
        mesh: Panel mesh object.
        gamma: Circulation strengths.
        V_ms: Freestream speed.
        S_ref: Reference area.

    Returns:
        CDi_tp: Induced drag coefficient.
    """
    n = len(gamma)
    y_j, z_j = mesh.yC, mesh.zC
    yA_k, zA_k, yB_k, zB_k = mesh.yA, mesh.zA, mesh.yB, mesh.zB

    eps2 = (1e-6 * 2.0 * float(np.max(np.abs(y_j)))) ** 2
    w_tp = np.zeros(n)
    v_tp = np.zeros(n)

    for j in range(n):
        Pj_y, Pj_z = y_j[j], z_j[j]
        for k in range(n):
            gk = gamma[k]
            dyB, dzB = Pj_y - yB_k[k], Pj_z - zB_k[k]
            dyA, dzA = Pj_y - yA_k[k], Pj_z - zA_k[k]
            rB2, rA2 = max(dyB**2 + dzB**2, eps2), max(dyA**2 + dzA**2, eps2)

            w_tp[j] += (gk / (2.0 * np.pi)) * (dyB / rB2 - dyA / rA2)
            v_tp[j] += (-gk / (2.0 * np.pi)) * (dzB / rB2 - dzA / rA2)

    CDi_tp = -np.sum(gamma * (w_tp * (mesh.yB - mesh.yA) - v_tp * (mesh.zB - mesh.zA))) / (V_ms**2 * S_ref)
    return float(CDi_tp)

panel_mesh

Panel mesh generation for VLM aerodynamic analysis.

See: docs/theory/aerodynamics.md

PanelMesh dataclass

3D panel mesh for VLM analysis.

See: docs/theory/aerodynamics.md

Attributes:

Name	Type	Description
n_panels	int	Total panel count (int).
xA,	(yA, zA)	Bound vortex A-endpoint (inboard node) at quarter-chord.
xB,	(yB, zB)	Bound vortex B-endpoint (outboard node) at quarter-chord.
xC,	(yC, zC)	Control point (collocation point) at three-quarter-chord.
nx,	(ny, nz)	Unit normal at each control point.
chord_cp	ndarray	Chord length at each control-point.
y_cp	ndarray	Spanwise CP coordinate.
S_ref	float	Full planform reference area.

Note

All computation here is pure Python/NumPy. No @numba.njit functions are used in panel geometry. Extract array fields explicitly if passing to compiled functions.

Source code in src/lightaero/aerodynamics/panel_mesh.py

@dataclass
class PanelMesh:
    """3D panel mesh for VLM analysis.

    See: docs/theory/aerodynamics.md

    Attributes:
        n_panels: Total panel count (int).
        xA, yA, zA: Bound vortex A-endpoint (inboard node) at quarter-chord.
        xB, yB, zB: Bound vortex B-endpoint (outboard node) at quarter-chord.
        xC, yC, zC: Control point (collocation point) at three-quarter-chord.
        nx, ny, nz: Unit normal at each control point.
        chord_cp: Chord length at each control-point.
        y_cp: Spanwise CP coordinate.
        S_ref: Full planform reference area.

    Note:
        All computation here is pure Python/NumPy. No @numba.njit functions
        are used in panel geometry. Extract array fields explicitly if
        passing to compiled functions.
    """

    n_panels: int
    xA: np.ndarray
    yA: np.ndarray
    zA: np.ndarray
    xB: np.ndarray
    yB: np.ndarray
    zB: np.ndarray
    xC: np.ndarray
    yC: np.ndarray
    zC: np.ndarray
    nx: np.ndarray
    ny: np.ndarray
    nz: np.ndarray
    chord_cp: np.ndarray
    y_cp: np.ndarray
    S_ref: float

build_panel_mesh(wing)

Convert WingGeometry into a full-span 3D PanelMesh.

See: docs/theory/aerodynamics.md

Parameters:

Name	Type	Description	Default
wing	WingGeometry	Frozen WingGeometry from build_wing_geometry().	required

Returns:

Name	Type	Description
PanelMesh	PanelMesh	The generated mesh.

Source code in src/lightaero/aerodynamics/panel_mesh.py

def build_panel_mesh(wing: WingGeometry) -> PanelMesh:
    """Convert WingGeometry into a full-span 3D PanelMesh.

    See: docs/theory/aerodynamics.md

    Args:
        wing: Frozen WingGeometry from build_wing_geometry().

    Returns:
        PanelMesh: The generated mesh.
    """
    n = wing.n_panels  # number of panels per semi-span

    # ------------------------------------------------------------------
    # Step 1: Spanwise nodes for the right semi-span (n+1 nodes)
    # Using the same cosine-spacing formula as WingGeometry._cosine_y_cp
    # but for nodes (not CPs) so that CPs from wing.y_cp are consistent.
    # ------------------------------------------------------------------
    theta = np.linspace(0.0, np.pi, n + 1)
    y_nodes = (wing.half_span / 2.0) * (1.0 - np.cos(theta))  # shape (n+1,)

    # ------------------------------------------------------------------
    # Step 2: Leading-edge x-offset and chord at spanwise nodes
    # x_le(y) = cumulative integral of tan(sweep_le) over y_stations
    # ------------------------------------------------------------------
    # x_le at each y_station break (cumulative sum of dx per segment)
    dx_le = np.diff(wing.y_stations) * np.tan(wing.sweep_le[:-1])
    x_le_at_stations = np.concatenate([[0.0], np.cumsum(dx_le)])  # shape (n_sections,)

    # Interpolate to nodes and to CPs
    x_le_nodes = np.interp(y_nodes, wing.y_stations, x_le_at_stations)  # shape (n+1,)
    x_le_cp = np.interp(wing.y_cp, wing.y_stations, x_le_at_stations)  # shape (n,)

    chord_nodes = np.interp(y_nodes, wing.y_stations, wing.chord)  # shape (n+1,)

    # ------------------------------------------------------------------
    # Step 3: Camber heights at nodes and CPs for z-coordinates
    # zA/zB at nodes (x=0.25c); zC at CPs (x=0.75c)
    # ------------------------------------------------------------------
    sections_list = list(wing.airfoil_sections)
    x_qc = np.array([0.25])
    x_tqc = np.array([0.75])

    # z at bound vortex endpoints (quarter-chord of each panel node)
    z_at_nodes = np.array(
        [
            chord_nodes[j] * float(interpolate_camber(y_nodes[j], wing.y_stations, sections_list, x_qc)[0])
            for j in range(n + 1)
        ]
    )  # shape (n+1,)
    z_at_nodes += wing.deflection_z

    # z at control points (three-quarter-chord)
    z_cp = np.array(
        [
            wing.chord_cp[j] * float(interpolate_camber(wing.y_cp[j], wing.y_stations, sections_list, x_tqc)[0])
            for j in range(n)
        ]
    )  # shape (n,)
    deflection_cp = 0.5 * (wing.deflection_z[:-1] + wing.deflection_z[1:])
    z_cp += deflection_cp

    # ------------------------------------------------------------------
    # Step 4: Bound vortex endpoints (right semi-span, panels 0..n-1)
    # Panel j runs from node j (inboard) to node j+1 (outboard)
    # ------------------------------------------------------------------
    xA_r = x_le_nodes[:-1] + 0.25 * chord_nodes[:-1]  # shape (n,)
    yA_r = y_nodes[:-1]
    zA_r = z_at_nodes[:-1]

    xB_r = x_le_nodes[1:] + 0.25 * chord_nodes[1:]  # shape (n,)
    yB_r = y_nodes[1:]
    zB_r = z_at_nodes[1:]

    # ------------------------------------------------------------------
    # Step 5: Control points (right semi-span, from wing.y_cp)
    # ------------------------------------------------------------------
    xC_r = x_le_cp + 0.75 * wing.chord_cp  # shape (n,)
    yC_r = wing.y_cp.copy()
    zC_r = z_cp

    # ------------------------------------------------------------------
    # Step 6: Unit normals at control points
    # The effective angle alpha_eff = twist + camberline slope at 3qc.
    # Normal hat = (-sin(alpha_eff), 0, cos(alpha_eff)) -- then normalised.
    # For NACA sections with ny=0 (no sweep-induced spanwise normal component
    # in strip theory), the normal lies in the xz-plane.
    # ------------------------------------------------------------------
    # Camberline slope dz_c/dx at x=0.75c for each CP
    slopes = np.zeros(n)
    for j in range(n):
        y_j = wing.y_cp[j]
        # Determine the two bracketing airfoil sections
        idx = int(
            np.clip(
                np.searchsorted(wing.y_stations, y_j, side="right") - 1,
                0,
                len(wing.y_stations) - 2,
            )
        )
        y0 = float(wing.y_stations[idx])
        y1 = float(wing.y_stations[idx + 1])
        alpha = (y_j - y0) / (y1 - y0) if y1 > y0 else 0.0

        # Slope at x=0.75c from section idx
        slope_0 = _camber_slope_at_3qc(wing.airfoil_sections[idx])
        # Slope at x=0.75c from section idx+1
        slope_1 = _camber_slope_at_3qc(wing.airfoil_sections[idx + 1])

        slopes[j] = (1.0 - alpha) * slope_0 + alpha * slope_1

    alpha_eff = wing.twist_cp + slopes  # shape (n,)

    nx_r = -np.sin(alpha_eff)
    ny_r = np.zeros(n)
    nz_r = np.cos(alpha_eff)

    # Normalise to unit length (already ~1 for small angles, but be exact)
    mag = np.sqrt(nx_r**2 + ny_r**2 + nz_r**2)
    nx_r = nx_r / mag
    ny_r = ny_r / mag
    nz_r = nz_r / mag

    # ------------------------------------------------------------------
    # Step 7: Mirror for left semi-span (panels n..2n-1)
    # The left semi-span uses reversed panel ordering so that the
    # inboard-to-outboard convention is preserved (outboard = more negative y).
    # Endpoint A on the right becomes B on the left (and vice versa).
    # ------------------------------------------------------------------
    xA_l = xB_r[::-1]
    yA_l = -yB_r[::-1]
    zA_l = zB_r[::-1]

    xB_l = xA_r[::-1]
    yB_l = -yA_r[::-1]
    zB_l = zA_r[::-1]

    xC_l = xC_r[::-1]
    yC_l = -yC_r[::-1]
    zC_l = zC_r[::-1]

    nx_l = nx_r[::-1]
    ny_l = -ny_r[::-1]  # Invert spanwise normal component only
    nz_l = nz_r[::-1]

    chord_cp_l = wing.chord_cp[::-1]

    # ------------------------------------------------------------------
    # Assemble full-span arrays (left then right)
    # The left semi-span occupies indices 0 .. n-1.
    # The right semi-span occupies indices n .. 2n-1.
    # This results in a continuous -tip to +tip spanwise ordering.
    # ------------------------------------------------------------------
    xA = np.concatenate([xA_l, xA_r])
    yA = np.concatenate([yA_l, yA_r])
    zA = np.concatenate([zA_l, zA_r])

    xB = np.concatenate([xB_l, xB_r])
    yB = np.concatenate([yB_l, yB_r])
    zB = np.concatenate([zB_l, zB_r])

    xC = np.concatenate([xC_l, xC_r])
    yC = np.concatenate([yC_l, yC_r])
    zC = np.concatenate([zC_l, zC_r])

    nx = np.concatenate([nx_l, nx_r])
    ny = np.concatenate([ny_l, ny_r])
    nz = np.concatenate([nz_l, nz_r])

    chord_cp_full = np.concatenate([chord_cp_l, wing.chord_cp])
    y_cp_full = np.concatenate([yC_l, yC_r])

    # ------------------------------------------------------------------
    # Step 8: Reference area (full span, trapz over right semi-span * 2)
    # ------------------------------------------------------------------
    S_ref = float(np.trapezoid(wing.chord_cp, wing.y_cp)) * 2.0

    return PanelMesh(
        n_panels=2 * n,
        xA=xA,
        yA=yA,
        zA=zA,
        xB=xB,
        yB=yB,
        zB=zB,
        xC=xC,
        yC=yC,
        zC=zC,
        nx=nx,
        ny=ny,
        nz=nz,
        chord_cp=chord_cp_full,
        y_cp=y_cp_full,
        S_ref=S_ref,
    )

solver

VLM discipline wrapper: VLMDiscipline.

See: docs/theory/aerodynamics.md

VLMDiscipline

Bases: DisciplineBase

Low-fidelity VLM aerodynamics discipline.

See: docs/theory/aerodynamics.md

Note

Inputs: wing (WingGeometry), alpha_rad (float), V_ms (float), altitude_m (float). Returns: AeroOutput TypedDict. AIC matrix is cached based on geometry hash.

Source code in src/lightaero/aerodynamics/solver.py

@register("aero", "vlm")
class VLMDiscipline(DisciplineBase):
    """Low-fidelity VLM aerodynamics discipline.

    See: docs/theory/aerodynamics.md

    Note:
        Inputs: wing (WingGeometry), alpha_rad (float), V_ms (float), altitude_m (float).
        Returns: AeroOutput TypedDict.
        AIC matrix is cached based on geometry hash.
    """

    def __init__(self) -> None:
        self._aic_cache: np.ndarray | None = None
        self._aic_key: tuple | None = None
        self._results: AeroOutput | None = None

    def __call__(self, **inputs: Any) -> AeroOutput:
        """Evaluate VLM aerodynamics at the given flight condition.

        See: docs/theory/aerodynamics.md

        Args:
            **inputs (Any): wing, alpha_rad, V_ms (or M), altitude_m.

        Returns:
            AeroOutput: Validated results.
        """

        wing: WingGeometry = inputs["wing"]
        altitude_m: float = float(inputs["altitude_m"])
        rho, T_K, _, a = isa(altitude_m)

        if inputs.get("M", None) is not None and inputs.get("V_ms", None) is not None:
            raise ValueError("VLMDiscipline only accepts M or V_ms but not both")
        elif inputs.get("M", None) is None and inputs.get("V_ms", None) is None:
            raise ValueError("VLMDiscipline requires either M or V_ms as input.")
        elif inputs.get("M", None) is None and inputs.get("V_ms", None) is not None:
            V_ms: float = float(inputs["V_ms"])
            M = V_ms / a
        else:
            M: float = float(inputs["M"])
            V_ms = M * a

        if M >= 1.0:
            raise ValueError(f"Mach number {M:.2f} >= 1.0. Prandtl-Glauert correction is valid only for subsonic flow.")

        # Prandtl-Glauert compressibility factor
        beta = math.sqrt(1.0 - M**2)

        mesh = build_panel_mesh(wing)

        # Handle spanwise-varying alpha_rad (scalar, station-wise, or panel-wise)
        alpha_in = inputs["alpha_rad"]
        if isinstance(alpha_in, (list, np.ndarray)):
            alpha_in = np.asarray(alpha_in, dtype=np.float64)
            if len(alpha_in) == len(wing.y_stations):
                # Interpolate from stations to right semi-span control points
                alpha_cp_right = np.interp(wing.y_cp, wing.y_stations, alpha_in)
            elif len(alpha_in) == wing.n_panels:
                alpha_cp_right = alpha_in
            else:
                raise ValueError(
                    f"alpha_rad array length {len(alpha_in)} must match "
                    f"n_stations ({len(wing.y_stations)}) or n_panels ({wing.n_panels})"
                )
            # Full-span mirror: left semi-span (-tip -> root) + right semi-span (root -> tip)
            alpha_rad = np.concatenate([alpha_cp_right[::-1], alpha_cp_right])
        else:
            alpha_rad = float(alpha_in)

        # Safety Guards: Check for validity regime (Mach < 0.3, AoA < 10 deg)
        check_regime_validity(M, float(np.max(np.abs(alpha_rad))))

        epsilon = 1e-6 * 2.0 * wing.half_span

        cache_key = (
            wing.chord.tobytes(),
            wing.sweep_le.tobytes(),
            wing.y_stations.tobytes(),
            wing.n_panels,
            wing.half_span,
        )
        if self._aic_key != cache_key or self._aic_cache is None:
            self._aic_cache = assemble_aic(mesh, epsilon=epsilon)
            self._aic_key = cache_key
        AIC = self._aic_cache

        RHS = build_rhs(mesh, alpha_rad=alpha_rad, V_ms=V_ms)
        gamma = np.linalg.solve(AIC, RHS)

        CL, CDi_kj, CM, lift_N, drag_i_N, span_cl, span_load, span_drag = near_field_forces(
            mesh,
            gamma,
            rho,
            V_ms,
            mesh.S_ref,
            alpha_rad,
        )
        CDi_tp = trefftz_cdi(mesh, gamma, V_ms, mesh.S_ref)

        # Panel-local cos-sweep Prandtl-Glauert correction (A3).
        # Each strip sees an effective Mach M_eff = M * cos(local LE sweep).
        # This replaces the global beta = sqrt(1 - M^2) for span quantities.
        sweep_cp = np.interp(wing.y_cp, wing.y_stations, wing.sweep_le)
        M_eff_cp = np.clip(M * np.cos(sweep_cp), 0.0, 0.99)
        beta_cp = np.sqrt(1.0 - M_eff_cp**2)
        beta_cp = np.maximum(beta_cp, 0.01)  # numerical floor

        # Full-span mirror: left semi-span (-tip -> root) + right semi-span (root -> tip)
        beta_full = np.concatenate([beta_cp[::-1], beta_cp])

        span_cl = span_cl / beta_full
        span_load = span_load / beta_full
        span_drag = span_drag / beta_full

        # Reintegrate CL and lift_N from corrected spanwise distributions (right semi-span)
        q_inf = 0.5 * rho * V_ms**2
        CL = 2.0 * float(np.trapezoid(span_cl[wing.n_panels :] * wing.chord_cp, wing.y_cp)) / mesh.S_ref
        lift_N = CL * q_inf * mesh.S_ref

        # CDi, CM, drag_N: global Prandtl-Glauert (less sensitive to span variation)
        CDi_tp /= beta
        CM /= beta

        # Profile drag via strip theory (right semi-span only; symmetric)
        tc_stations = np.array([float(sec.thickness_at(np.array([0.35]))[0]) for sec in wing.airfoil_sections])
        tc_cp = np.interp(wing.y_cp, wing.y_stations, tc_stations)

        CD0 = profile_drag_cd0(
            rho=rho,
            V_ms=V_ms,
            T_K=T_K,
            S_ref=mesh.S_ref,
            y_cp_half=wing.y_cp,
            chord_cp_half=wing.chord_cp,
            tc_cp_half=tc_cp,
        )
        CD_total = float(CDi_tp) + CD0

        # Total Drag Force: Induced (corrected) + Viscous Profile
        drag_i_total_N = CDi_tp * q_inf * mesh.S_ref
        drag_v_total_N = CD0 * q_inf * mesh.S_ref
        drag_total_N = drag_i_total_N + drag_v_total_N

        self._results: AeroOutput = AeroOutput(
            CL=float(CL),
            CDi=float(CDi_tp),
            CD0=CD0,
            CD=CD_total,
            CM=float(CM),
            lift_N=float(lift_N),
            drag_N=float(drag_total_N),
            span_cl=span_cl,
            span_load=span_load,
            span_drag=span_drag,
            circulation=gamma,
        )

        if inputs.get("validate_output", True):
            validate_discipline_output(self._results, AERO_OUTPUT_SPEC)
        return self._results

    def __str__(self) -> str:
        if self._results is None:
            return "VLMDiscipline: No results yet. Call the discipline with valid inputs to compute aerodynamics."
        return "\n".join(
            [
                f"Lift Coefficient           (CL) : {self._results.CL:.4e}",
                f"Drag Coefficient           (CD) : {self._results.CD:.4e}",
                f"Induced Drag Coefficient   (CDi): {self._results.CDi:.4e}",
                f"Zero-Lift Drag Coefficient (CD0): {self._results.CD0:.4e}",
                f"Coefficient of Moment      (CM) : {self._results.CM:.4e}",
            ]
        )

__call__(**inputs)

Evaluate VLM aerodynamics at the given flight condition.

See: docs/theory/aerodynamics.md

Parameters:

Name	Type	Description	Default
**inputs	Any	wing, alpha_rad, V_ms (or M), altitude_m.	{}

Returns:

Name	Type	Description
AeroOutput	AeroOutput	Validated results.

Source code in src/lightaero/aerodynamics/solver.py

def __call__(self, **inputs: Any) -> AeroOutput:
    """Evaluate VLM aerodynamics at the given flight condition.

    See: docs/theory/aerodynamics.md

    Args:
        **inputs (Any): wing, alpha_rad, V_ms (or M), altitude_m.

    Returns:
        AeroOutput: Validated results.
    """

    wing: WingGeometry = inputs["wing"]
    altitude_m: float = float(inputs["altitude_m"])
    rho, T_K, _, a = isa(altitude_m)

    if inputs.get("M", None) is not None and inputs.get("V_ms", None) is not None:
        raise ValueError("VLMDiscipline only accepts M or V_ms but not both")
    elif inputs.get("M", None) is None and inputs.get("V_ms", None) is None:
        raise ValueError("VLMDiscipline requires either M or V_ms as input.")
    elif inputs.get("M", None) is None and inputs.get("V_ms", None) is not None:
        V_ms: float = float(inputs["V_ms"])
        M = V_ms / a
    else:
        M: float = float(inputs["M"])
        V_ms = M * a

    if M >= 1.0:
        raise ValueError(f"Mach number {M:.2f} >= 1.0. Prandtl-Glauert correction is valid only for subsonic flow.")

    # Prandtl-Glauert compressibility factor
    beta = math.sqrt(1.0 - M**2)

    mesh = build_panel_mesh(wing)

    # Handle spanwise-varying alpha_rad (scalar, station-wise, or panel-wise)
    alpha_in = inputs["alpha_rad"]
    if isinstance(alpha_in, (list, np.ndarray)):
        alpha_in = np.asarray(alpha_in, dtype=np.float64)
        if len(alpha_in) == len(wing.y_stations):
            # Interpolate from stations to right semi-span control points
            alpha_cp_right = np.interp(wing.y_cp, wing.y_stations, alpha_in)
        elif len(alpha_in) == wing.n_panels:
            alpha_cp_right = alpha_in
        else:
            raise ValueError(
                f"alpha_rad array length {len(alpha_in)} must match "
                f"n_stations ({len(wing.y_stations)}) or n_panels ({wing.n_panels})"
            )
        # Full-span mirror: left semi-span (-tip -> root) + right semi-span (root -> tip)
        alpha_rad = np.concatenate([alpha_cp_right[::-1], alpha_cp_right])
    else:
        alpha_rad = float(alpha_in)

    # Safety Guards: Check for validity regime (Mach < 0.3, AoA < 10 deg)
    check_regime_validity(M, float(np.max(np.abs(alpha_rad))))

    epsilon = 1e-6 * 2.0 * wing.half_span

    cache_key = (
        wing.chord.tobytes(),
        wing.sweep_le.tobytes(),
        wing.y_stations.tobytes(),
        wing.n_panels,
        wing.half_span,
    )
    if self._aic_key != cache_key or self._aic_cache is None:
        self._aic_cache = assemble_aic(mesh, epsilon=epsilon)
        self._aic_key = cache_key
    AIC = self._aic_cache

    RHS = build_rhs(mesh, alpha_rad=alpha_rad, V_ms=V_ms)
    gamma = np.linalg.solve(AIC, RHS)

    CL, CDi_kj, CM, lift_N, drag_i_N, span_cl, span_load, span_drag = near_field_forces(
        mesh,
        gamma,
        rho,
        V_ms,
        mesh.S_ref,
        alpha_rad,
    )
    CDi_tp = trefftz_cdi(mesh, gamma, V_ms, mesh.S_ref)

    # Panel-local cos-sweep Prandtl-Glauert correction (A3).
    # Each strip sees an effective Mach M_eff = M * cos(local LE sweep).
    # This replaces the global beta = sqrt(1 - M^2) for span quantities.
    sweep_cp = np.interp(wing.y_cp, wing.y_stations, wing.sweep_le)
    M_eff_cp = np.clip(M * np.cos(sweep_cp), 0.0, 0.99)
    beta_cp = np.sqrt(1.0 - M_eff_cp**2)
    beta_cp = np.maximum(beta_cp, 0.01)  # numerical floor

    # Full-span mirror: left semi-span (-tip -> root) + right semi-span (root -> tip)
    beta_full = np.concatenate([beta_cp[::-1], beta_cp])

    span_cl = span_cl / beta_full
    span_load = span_load / beta_full
    span_drag = span_drag / beta_full

    # Reintegrate CL and lift_N from corrected spanwise distributions (right semi-span)
    q_inf = 0.5 * rho * V_ms**2
    CL = 2.0 * float(np.trapezoid(span_cl[wing.n_panels :] * wing.chord_cp, wing.y_cp)) / mesh.S_ref
    lift_N = CL * q_inf * mesh.S_ref

    # CDi, CM, drag_N: global Prandtl-Glauert (less sensitive to span variation)
    CDi_tp /= beta
    CM /= beta

    # Profile drag via strip theory (right semi-span only; symmetric)
    tc_stations = np.array([float(sec.thickness_at(np.array([0.35]))[0]) for sec in wing.airfoil_sections])
    tc_cp = np.interp(wing.y_cp, wing.y_stations, tc_stations)

    CD0 = profile_drag_cd0(
        rho=rho,
        V_ms=V_ms,
        T_K=T_K,
        S_ref=mesh.S_ref,
        y_cp_half=wing.y_cp,
        chord_cp_half=wing.chord_cp,
        tc_cp_half=tc_cp,
    )
    CD_total = float(CDi_tp) + CD0

    # Total Drag Force: Induced (corrected) + Viscous Profile
    drag_i_total_N = CDi_tp * q_inf * mesh.S_ref
    drag_v_total_N = CD0 * q_inf * mesh.S_ref
    drag_total_N = drag_i_total_N + drag_v_total_N

    self._results: AeroOutput = AeroOutput(
        CL=float(CL),
        CDi=float(CDi_tp),
        CD0=CD0,
        CD=CD_total,
        CM=float(CM),
        lift_N=float(lift_N),
        drag_N=float(drag_total_N),
        span_cl=span_cl,
        span_load=span_load,
        span_drag=span_drag,
        circulation=gamma,
    )

    if inputs.get("validate_output", True):
        validate_discipline_output(self._results, AERO_OUTPUT_SPEC)
    return self._results

atmosphere

isa

ISA atmosphere model (ISO 2533 / ICAO Doc 7488).

See: docs/theory/atmosphere.md

isa(altitude_m)

ISA atmosphere at the given altitude above MSL.

See: docs/theory/atmosphere.md

Parameters:

Name	Type	Description	Default
altitude_m	float	Geometric altitude in metres.	required

Returns:

Name	Type	Description
tuple	tuple[float, float, float, float]	(rho, T, p, a)

Source code in src/lightaero/atmosphere/isa.py

def isa(altitude_m: float) -> tuple[float, float, float, float]:
    """ISA atmosphere at the given altitude above MSL.

    See: docs/theory/atmosphere.md

    Args:
        altitude_m: Geometric altitude in metres.

    Returns:
        tuple: (rho, T, p, a)
    """
    return _isa_scalar(float(altitude_m))

geometry

build_crm_geometry(n_panels=40)

Build NASA CRM wing geometry (trapezoidal approximation).

Source code in src/lightaero/geometry/crm.py

def build_crm_geometry(n_panels: int = 40) -> WingGeometry:
    """Build NASA CRM wing geometry (trapezoidal approximation)."""
    airfoil = AirfoilSection.from_naca4("0012")
    return build_wing_geometry(
        half_span=_CRM_HALF_SPAN,
        n_panels=n_panels,
        y_stations=np.array([0.0, _CRM_HALF_SPAN]),
        chord=np.array([_CRM_ROOT_CHORD, _CRM_TIP_CHORD]),
        sweep_le=np.array([_CRM_SWEEP_LE, _CRM_SWEEP_LE]),
        taper=np.array([1.0, _CRM_TIP_CHORD / _CRM_ROOT_CHORD]),
        twist=np.zeros(2),
        airfoil_sections=(airfoil, airfoil),
    )

build_dlrf4_geometry(n_panels=40)

Return the DLR-F4 isolated wing as a WingGeometry.

Parameters:

Name	Type	Description	Default
n_panels	int	Number of semi-span VLM panels.	40

Returns:

Type	Description
WingGeometry	Frozen WingGeometry with DLR-F4 trapezoidal planform.

Source code in src/lightaero/geometry/dlrf4.py

def build_dlrf4_geometry(n_panels: int = 40) -> WingGeometry:
    """Return the DLR-F4 isolated wing as a WingGeometry.

    Args:
        n_panels: Number of semi-span VLM panels.

    Returns:
        Frozen WingGeometry with DLR-F4 trapezoidal planform.
    """
    root = AirfoilSection.from_naca4("0012")
    tip = AirfoilSection.from_naca4("0012")
    return build_wing_geometry(
        half_span=_HALF_SPAN,
        n_panels=n_panels,
        y_stations=np.array([0.0, _HALF_SPAN]),
        chord=np.array([_ROOT_CHORD, _TIP_CHORD]),
        sweep_le=np.array([_SWEEP_LE, _SWEEP_LE]),
        taper=np.array([1.0, _TIP_CHORD / _ROOT_CHORD]),
        twist=np.array([0.0, 0.0]),
        airfoil_sections=(root, tip),
    )

build_dlrf6_geometry(n_panels=40)

Return the DLR-F6 isolated wing as a WingGeometry.

Parameters:

Name	Type	Description	Default
n_panels	int	Number of semi-span VLM panels.	40

Returns:

Type	Description
WingGeometry	Frozen WingGeometry with DLR-F6 cranked trapezoidal planform.

Source code in src/lightaero/geometry/dlrf6.py

def build_dlrf6_geometry(n_panels: int = 40) -> WingGeometry:
    """Return the DLR-F6 isolated wing as a WingGeometry.

    Args:
        n_panels: Number of semi-span VLM panels.

    Returns:
        Frozen WingGeometry with DLR-F6 cranked trapezoidal planform.
    """
    root = AirfoilSection.from_naca4("0012")
    brk = AirfoilSection.from_naca4("0012")
    tip = AirfoilSection.from_naca4("0012")
    return build_wing_geometry(
        half_span=_HALF_SPAN,
        n_panels=n_panels,
        y_stations=np.array([0.0, _BREAK_Y, _HALF_SPAN]),
        chord=np.array([_ROOT_CHORD, _BREAK_CHORD, _TIP_CHORD]),
        sweep_le=np.array([_SWEEP_LE, _SWEEP_LE, _SWEEP_LE]),
        taper=np.array([1.0, _BREAK_CHORD / _ROOT_CHORD, _TIP_CHORD / _ROOT_CHORD]),
        twist=np.array([0.0, 0.0, 0.0]),
        airfoil_sections=(root, brk, tip),
    )

build_ucrm_geometry(n_panels=40)

Build uCRM geometry (jig/unmodified).

Source code in src/lightaero/geometry/ucrm.py

def build_ucrm_geometry(n_panels: int = 40) -> WingGeometry:
    """Build uCRM geometry (jig/unmodified)."""
    airfoil = AirfoilSection.from_naca4("0015")
    return build_wing_geometry(
        half_span=_UCRM_GEOMETRY["yle"].values[-1] * _INCH_TO_METERS,
        n_panels=n_panels,
        y_stations=_UCRM_GEOMETRY["yle"].values * _INCH_TO_METERS,
        chord=_UCRM_GEOMETRY["c-plan"].values * _INCH_TO_METERS,
        sweep_le=np.full(len(_UCRM_GEOMETRY), UCRM_SWEEP_LE),
        taper=_UCRM_GEOMETRY["c-plan"].values / _UCRM_GEOMETRY["c-plan"].values[0],
        twist=np.deg2rad(_UCRM_GEOMETRY["twist"].values),
        airfoil_sections=tuple([airfoil] * len(_UCRM_GEOMETRY)),
    )

airfoil

Airfoil section primitives for Foundation.

See: docs/theory/index.md for detailed NACA equations and interpolation logic.

AirfoilSection dataclass

Airfoil section with camber and thickness query methods.

See: docs/theory/index.md for theoretical background and NACA parsing details.

Source code in src/lightaero/geometry/airfoil.py

@dataclass
class AirfoilSection:
    """Airfoil section with camber and thickness query methods.

    See: docs/theory/index.md for theoretical background and NACA parsing details.
    """

    _M: float | None = field(default=None, repr=False)
    _P: float | None = field(default=None, repr=False)
    _T: float | None = field(default=None, repr=False)

    _x_upper: np.ndarray | None = field(default=None, repr=False)
    _y_upper: np.ndarray | None = field(default=None, repr=False)
    _x_lower: np.ndarray | None = field(default=None, repr=False)
    _y_lower: np.ndarray | None = field(default=None, repr=False)

    @classmethod
    def from_naca4(cls, designation: str) -> AirfoilSection:
        """Construct from NACA 4-digit string, e.g. '2412'."""
        d = designation.strip().upper().replace("NACA", "").strip()
        if len(d) != 4:
            raise ValueError(f"Expected 4-digit NACA designation, got {designation!r}")
        M = int(d[0]) / 100.0
        P = int(d[1]) / 10.0
        T = int(d[2:]) / 100.0
        return cls(_M=M, _P=P if P > 0.0 else 0.4, _T=T)

    def camber_at(self, x_over_c: np.ndarray) -> np.ndarray:
        """Camberline ordinate z_c/c at chordwise positions x/c."""
        x = np.asarray(x_over_c, dtype=float)
        if x.min() < 0.0 or x.max() > 1.0:
            raise ValueError("x_over_c values must be in [0, 1]")
        if self._M is not None:
            return naca4_camber(x, self._M, self._P)
        yu = np.interp(x, self._x_upper, self._y_upper)
        yl = np.interp(x, self._x_lower, self._y_lower)
        return 0.5 * (yu + yl)

    def thickness_at(self, x_over_c: np.ndarray) -> np.ndarray:
        """Full thickness t/c at chordwise positions x/c."""
        x = np.asarray(x_over_c, dtype=float)
        if self._T is not None:
            return 2.0 * naca4_thickness(x, self._T)
        yu = np.interp(x, self._x_upper, self._y_upper)
        yl = np.interp(x, self._x_lower, self._y_lower)
        return yu - yl

    def get_profile(self, n_points: int = 50) -> tuple[np.ndarray, np.ndarray]:
        """Return (x, z) coordinates for a closed profile, normalized to chord=1."""
        x = np.linspace(0, 1, n_points)
        if self._T is not None:
            zc = naca4_camber(x, self._M, self._P)
            zt = naca4_thickness(x, self._T)
            z_upper = zc + zt
            z_lower = zc - zt
        else:
            yu = np.interp(x, self._x_upper, self._y_upper)
            yl = np.interp(x, self._x_lower, self._y_lower)
            z_upper = yu
            z_lower = yl

        x_full = np.concatenate([x, x[::-1]])
        z_full = np.concatenate([z_upper, z_lower[::-1]])
        return x_full, z_full

camber_at(x_over_c)

Camberline ordinate z_c/c at chordwise positions x/c.

Source code in src/lightaero/geometry/airfoil.py

def camber_at(self, x_over_c: np.ndarray) -> np.ndarray:
    """Camberline ordinate z_c/c at chordwise positions x/c."""
    x = np.asarray(x_over_c, dtype=float)
    if x.min() < 0.0 or x.max() > 1.0:
        raise ValueError("x_over_c values must be in [0, 1]")
    if self._M is not None:
        return naca4_camber(x, self._M, self._P)
    yu = np.interp(x, self._x_upper, self._y_upper)
    yl = np.interp(x, self._x_lower, self._y_lower)
    return 0.5 * (yu + yl)

from_naca4(designation) classmethod

Construct from NACA 4-digit string, e.g. '2412'.

Source code in src/lightaero/geometry/airfoil.py

@classmethod
def from_naca4(cls, designation: str) -> AirfoilSection:
    """Construct from NACA 4-digit string, e.g. '2412'."""
    d = designation.strip().upper().replace("NACA", "").strip()
    if len(d) != 4:
        raise ValueError(f"Expected 4-digit NACA designation, got {designation!r}")
    M = int(d[0]) / 100.0
    P = int(d[1]) / 10.0
    T = int(d[2:]) / 100.0
    return cls(_M=M, _P=P if P > 0.0 else 0.4, _T=T)

get_profile(n_points=50)

Return (x, z) coordinates for a closed profile, normalized to chord=1.

Source code in src/lightaero/geometry/airfoil.py

def get_profile(self, n_points: int = 50) -> tuple[np.ndarray, np.ndarray]:
    """Return (x, z) coordinates for a closed profile, normalized to chord=1."""
    x = np.linspace(0, 1, n_points)
    if self._T is not None:
        zc = naca4_camber(x, self._M, self._P)
        zt = naca4_thickness(x, self._T)
        z_upper = zc + zt
        z_lower = zc - zt
    else:
        yu = np.interp(x, self._x_upper, self._y_upper)
        yl = np.interp(x, self._x_lower, self._y_lower)
        z_upper = yu
        z_lower = yl

    x_full = np.concatenate([x, x[::-1]])
    z_full = np.concatenate([z_upper, z_lower[::-1]])
    return x_full, z_full

thickness_at(x_over_c)

Full thickness t/c at chordwise positions x/c.

Source code in src/lightaero/geometry/airfoil.py

def thickness_at(self, x_over_c: np.ndarray) -> np.ndarray:
    """Full thickness t/c at chordwise positions x/c."""
    x = np.asarray(x_over_c, dtype=float)
    if self._T is not None:
        return 2.0 * naca4_thickness(x, self._T)
    yu = np.interp(x, self._x_upper, self._y_upper)
    yl = np.interp(x, self._x_lower, self._y_lower)
    return yu - yl

blended_profile(y_query, y_stations, sections, n_points)

Linearly blend normalized airfoil profiles between two spanwise stations.

Source code in src/lightaero/geometry/airfoil.py

def blended_profile(
    y_query: float,
    y_stations: np.ndarray,
    sections: list,
    n_points: int,
) -> tuple[np.ndarray, np.ndarray]:
    """Linearly blend normalized airfoil profiles between two spanwise stations."""
    y_stations = np.asarray(y_stations, dtype=float)
    y_query = float(np.clip(y_query, y_stations[0], y_stations[-1]))

    idx = int(
        np.clip(
            np.searchsorted(y_stations, y_query, side="right") - 1,
            0,
            len(y_stations) - 2,
        )
    )
    y0, y1 = float(y_stations[idx]), float(y_stations[idx + 1])
    alpha = 0.0 if y1 == y0 else (y_query - y0) / (y1 - y0)

    x0, z0 = sections[idx].get_profile(n_points)
    x1, z1 = sections[idx + 1].get_profile(n_points)

    return (1.0 - alpha) * x0 + alpha * x1, (1.0 - alpha) * z0 + alpha * z1

interpolate_camber(y_query, y_stations, sections, x_over_c)

Linear interpolation of camberline ordinates between spanwise stations.

See: docs/theory/index.md for the blending formula.

Source code in src/lightaero/geometry/airfoil.py

def interpolate_camber(
    y_query: float,
    y_stations: np.ndarray,
    sections: list[AirfoilSection],
    x_over_c: np.ndarray,
) -> np.ndarray:
    """Linear interpolation of camberline ordinates between spanwise stations.

    See: docs/theory/index.md for the blending formula.
    """
    y_stations = np.asarray(y_stations, dtype=float)
    x_over_c = np.asarray(x_over_c, dtype=float)
    y_query = float(np.clip(y_query, y_stations[0], y_stations[-1]))

    idx = np.searchsorted(y_stations, y_query, side="right") - 1
    idx = int(np.clip(idx, 0, len(y_stations) - 2))

    y0 = float(y_stations[idx])
    y1 = float(y_stations[idx + 1])

    if y1 == y0:
        return sections[idx].camber_at(x_over_c)

    alpha = (y_query - y0) / (y1 - y0)
    camber_0 = sections[idx].camber_at(x_over_c)
    camber_1 = sections[idx + 1].camber_at(x_over_c)

    return (1.0 - alpha) * camber_0 + alpha * camber_1

load_uiuc_dat(path)

Load airfoil coordinates from a UIUC .dat file (Selig or Lednicer format).

See: docs/theory/index.md for format details.

Source code in src/lightaero/geometry/airfoil.py

def load_uiuc_dat(path: Path) -> AirfoilSection:
    """Load airfoil coordinates from a UIUC .dat file (Selig or Lednicer format).

    See: docs/theory/index.md for format details.
    """
    path = Path(path)
    lines = path.read_text().splitlines()

    data_lines = []
    for line in lines:
        stripped = line.strip()
        if not stripped or stripped.startswith("#"):
            continue
        try:
            parts = stripped.split()
            floats = [float(p) for p in parts]
            data_lines.append(floats)
        except ValueError:
            continue

    if not data_lines:
        raise ValueError(f"No coordinate data found in {path}")

    first_line = data_lines[0]
    is_lednicer = len(first_line) == 2 and 1.0 < first_line[0] < 500.0 and 1.0 < first_line[1] < 500.0

    if is_lednicer:
        n_upper = int(round(first_line[0]))
        n_lower = int(round(first_line[1]))
        coords = data_lines[1:]
        upper = np.array(coords[:n_upper])
        lower = np.array(coords[n_upper : n_upper + n_lower])
        x_upper, y_upper = upper[:, 0], upper[:, 1]
        x_lower, y_lower = lower[:, 0], lower[:, 1]
    else:
        coords = np.array(data_lines)
        x_all = coords[:, 0]
        y_all = coords[:, 1]
        le_idx = int(np.argmin(x_all))
        x_upper = x_all[: le_idx + 1][::-1]
        y_upper = y_all[: le_idx + 1][::-1]
        x_lower = x_all[le_idx:]
        y_lower = y_all[le_idx:]

    x_max = max(float(x_upper.max()), float(x_lower.max()))
    if x_max > 1.01:
        x_upper = x_upper / x_max
        y_upper = y_upper / x_max
        x_lower = x_lower / x_max
        y_lower = y_lower / x_max

    return AirfoilSection(
        _x_upper=x_upper,
        _y_upper=y_upper,
        _x_lower=x_lower,
        _y_lower=y_lower,
    )

naca4_camber(x, M, P)

NACA 4-digit camberline ordinate z_c/c.

Parameters:

Name	Type	Description	Default
x	ndarray	Chordwise stations x/c in [0, 1].	required
M	float	Maximum camber as fraction of chord.	required
P	float	Chordwise position of maximum camber.	required

Returns:

Type	Description
ndarray	Camberline ordinate z_c/c.

Source code in src/lightaero/geometry/airfoil.py

def naca4_camber(x: np.ndarray, M: float, P: float) -> np.ndarray:
    """NACA 4-digit camberline ordinate z_c/c.

    Args:
        x: Chordwise stations x/c in [0, 1].
        M: Maximum camber as fraction of chord.
        P: Chordwise position of maximum camber.

    Returns:
        Camberline ordinate z_c/c.
    """
    x = np.asarray(x, dtype=float)
    if M <= 0.0:
        return np.zeros_like(x)
    yc = np.where(
        x < P,
        (M / P**2) * (2.0 * P * x - x**2),
        (M / (1.0 - P) ** 2) * (1.0 - 2.0 * P + 2.0 * P * x - x**2),
    )
    return yc

naca4_thickness(x, T)

NACA 4-digit half-thickness distribution t/c.

Parameters:

Name	Type	Description	Default
x	ndarray	Chordwise stations x/c in [0, 1].	required
T	float	Maximum thickness as fraction of chord.	required

Returns:

Type	Description
ndarray	Half-thickness t/c. Full thickness = 2 * returned value.

Source code in src/lightaero/geometry/airfoil.py

def naca4_thickness(x: np.ndarray, T: float) -> np.ndarray:
    """NACA 4-digit half-thickness distribution t/c.

    Args:
        x: Chordwise stations x/c in [0, 1].
        T: Maximum thickness as fraction of chord.

    Returns:
        Half-thickness t/c. Full thickness = 2 * returned value.
    """
    x = np.asarray(x, dtype=float)
    yt = (T / 0.2) * (0.2969 * np.sqrt(np.maximum(x, 0.0)) - 0.1260 * x - 0.3516 * x**2 + 0.2843 * x**3 - 0.1015 * x**4)
    return yt

crm

NASA Common Research Model (CRM) isolated wing geometry builder.

See: docs/theory/index.md for CRM planform parameters and references.

build_crm_geometry(n_panels=40)

Build NASA CRM wing geometry (trapezoidal approximation).

Source code in src/lightaero/geometry/crm.py

def build_crm_geometry(n_panels: int = 40) -> WingGeometry:
    """Build NASA CRM wing geometry (trapezoidal approximation)."""
    airfoil = AirfoilSection.from_naca4("0012")
    return build_wing_geometry(
        half_span=_CRM_HALF_SPAN,
        n_panels=n_panels,
        y_stations=np.array([0.0, _CRM_HALF_SPAN]),
        chord=np.array([_CRM_ROOT_CHORD, _CRM_TIP_CHORD]),
        sweep_le=np.array([_CRM_SWEEP_LE, _CRM_SWEEP_LE]),
        taper=np.array([1.0, _CRM_TIP_CHORD / _CRM_ROOT_CHORD]),
        twist=np.zeros(2),
        airfoil_sections=(airfoil, airfoil),
    )

dlrf4

DLR-F4 isolated wing geometry builder.

See: docs/theory/index.md for DLR-F4 planform parameters and references.

build_dlrf4_geometry(n_panels=40)

Return the DLR-F4 isolated wing as a WingGeometry.

Parameters:

Name	Type	Description	Default
n_panels	int	Number of semi-span VLM panels.	40

Returns:

Type	Description
WingGeometry	Frozen WingGeometry with DLR-F4 trapezoidal planform.

Source code in src/lightaero/geometry/dlrf4.py

def build_dlrf4_geometry(n_panels: int = 40) -> WingGeometry:
    """Return the DLR-F4 isolated wing as a WingGeometry.

    Args:
        n_panels: Number of semi-span VLM panels.

    Returns:
        Frozen WingGeometry with DLR-F4 trapezoidal planform.
    """
    root = AirfoilSection.from_naca4("0012")
    tip = AirfoilSection.from_naca4("0012")
    return build_wing_geometry(
        half_span=_HALF_SPAN,
        n_panels=n_panels,
        y_stations=np.array([0.0, _HALF_SPAN]),
        chord=np.array([_ROOT_CHORD, _TIP_CHORD]),
        sweep_le=np.array([_SWEEP_LE, _SWEEP_LE]),
        taper=np.array([1.0, _TIP_CHORD / _ROOT_CHORD]),
        twist=np.array([0.0, 0.0]),
        airfoil_sections=(root, tip),
    )

dlrf6

DLR-F6 isolated wing geometry builder.

See: docs/theory/index.md for DLR-F6 planform parameters and references.

build_dlrf6_geometry(n_panels=40)

Return the DLR-F6 isolated wing as a WingGeometry.

Parameters:

Name	Type	Description	Default
n_panels	int	Number of semi-span VLM panels.	40

Returns:

Type	Description
WingGeometry	Frozen WingGeometry with DLR-F6 cranked trapezoidal planform.

Source code in src/lightaero/geometry/dlrf6.py

def build_dlrf6_geometry(n_panels: int = 40) -> WingGeometry:
    """Return the DLR-F6 isolated wing as a WingGeometry.

    Args:
        n_panels: Number of semi-span VLM panels.

    Returns:
        Frozen WingGeometry with DLR-F6 cranked trapezoidal planform.
    """
    root = AirfoilSection.from_naca4("0012")
    brk = AirfoilSection.from_naca4("0012")
    tip = AirfoilSection.from_naca4("0012")
    return build_wing_geometry(
        half_span=_HALF_SPAN,
        n_panels=n_panels,
        y_stations=np.array([0.0, _BREAK_Y, _HALF_SPAN]),
        chord=np.array([_ROOT_CHORD, _BREAK_CHORD, _TIP_CHORD]),
        sweep_le=np.array([_SWEEP_LE, _SWEEP_LE, _SWEEP_LE]),
        taper=np.array([1.0, _BREAK_CHORD / _ROOT_CHORD, _TIP_CHORD / _ROOT_CHORD]),
        twist=np.array([0.0, 0.0, 0.0]),
        airfoil_sections=(root, brk, tip),
    )

ucrm

NASA uCRM isolated wing geometry builder.

See: docs/theory/index.md for uCRM planform parameters and references.

build_ucrm_geometry(n_panels=40)

Build uCRM geometry (jig/unmodified).

Source code in src/lightaero/geometry/ucrm.py

def build_ucrm_geometry(n_panels: int = 40) -> WingGeometry:
    """Build uCRM geometry (jig/unmodified)."""
    airfoil = AirfoilSection.from_naca4("0015")
    return build_wing_geometry(
        half_span=_UCRM_GEOMETRY["yle"].values[-1] * _INCH_TO_METERS,
        n_panels=n_panels,
        y_stations=_UCRM_GEOMETRY["yle"].values * _INCH_TO_METERS,
        chord=_UCRM_GEOMETRY["c-plan"].values * _INCH_TO_METERS,
        sweep_le=np.full(len(_UCRM_GEOMETRY), UCRM_SWEEP_LE),
        taper=_UCRM_GEOMETRY["c-plan"].values / _UCRM_GEOMETRY["c-plan"].values[0],
        twist=np.deg2rad(_UCRM_GEOMETRY["twist"].values),
        airfoil_sections=tuple([airfoil] * len(_UCRM_GEOMETRY)),
    )

wing

Parametric wing geometry.

See: docs/theory/index.md for parametric wing construction and cosine spacing details.

WingGeometry dataclass

Parametric wing geometry with cosine-spaced panel arrays.

See: docs/theory/index.md for detailed field descriptions and mathematical foundations.

Attributes:

Name	Type	Description
half_span	float	Semi-span in metres.
n_panels	int	Number of spanwise panels.
y_stations	ndarray	Spanwise break positions.
chord	ndarray	Chord lengths at each station.
sweep_le	ndarray	Leading-edge sweep angles (rad).
taper	ndarray	Local taper ratios.
twist	ndarray	Geometric twist angles (rad).
airfoil_sections	tuple	Tuple of AirfoilSection objects.
deflection_z	ndarray	Transverse nodal displacements (m).

Source code in src/lightaero/geometry/wing.py

@dataclass(frozen=True)
class WingGeometry:
    """Parametric wing geometry with cosine-spaced panel arrays.

    See: docs/theory/index.md for detailed field descriptions and mathematical foundations.

    Attributes:
        half_span: Semi-span in metres.
        n_panels: Number of spanwise panels.
        y_stations: Spanwise break positions.
        chord: Chord lengths at each station.
        sweep_le: Leading-edge sweep angles (rad).
        taper: Local taper ratios.
        twist: Geometric twist angles (rad).
        airfoil_sections: Tuple of AirfoilSection objects.
        deflection_z: Transverse nodal displacements (m).
    """

    # --- Planform inputs ---
    half_span: float
    n_panels: int
    y_stations: np.ndarray
    chord: np.ndarray
    sweep_le: np.ndarray
    taper: np.ndarray
    twist: np.ndarray
    airfoil_sections: tuple
    deflection_z: np.ndarray = None

    # --- Derived arrays ---
    y_cp: np.ndarray = field(init=False)
    chord_cp: np.ndarray = field(init=False)
    camber_cp: np.ndarray = field(init=False)
    twist_cp: np.ndarray = field(init=False)

    def __post_init__(self) -> None:
        """Compute derived arrays and lock all arrays as read-only."""
        y_cp = _cosine_y_cp(self.half_span, self.n_panels)

        if self.deflection_z is None:
            deflection_z = np.zeros(self.n_panels + 1)
        else:
            deflection_z = np.asarray(self.deflection_z, dtype=float)
            if deflection_z.shape != (self.n_panels + 1,):
                raise ValueError(f"deflection_z must have shape ({self.n_panels + 1},), got {deflection_z.shape}")

        chord_cp = np.interp(y_cp, self.y_stations, self.chord)
        twist_cp = np.interp(y_cp, self.y_stations, self.twist)

        x_mid = np.array([0.5])
        camber_cp = np.array(
            [
                float(
                    interpolate_camber(
                        y,
                        self.y_stations,
                        list(self.airfoil_sections),
                        x_mid,
                    )[0]
                )
                for y in y_cp
            ]
        )

        for arr_name, arr in [
            ("y_stations", self.y_stations),
            ("chord", self.chord),
            ("sweep_le", self.sweep_le),
            ("taper", self.taper),
            ("twist", self.twist),
        ]:
            arr_ro = np.array(arr)
            arr_ro.flags.writeable = False
            object.__setattr__(self, arr_name, arr_ro)

        for arr_name, arr in [
            ("y_cp", y_cp),
            ("chord_cp", chord_cp),
            ("twist_cp", twist_cp),
            ("camber_cp", camber_cp),
            ("deflection_z", deflection_z),
        ]:
            arr.flags.writeable = False
            object.__setattr__(self, arr_name, arr)

    def _section_at_y(self, y: float, n_profile_points: int) -> dict:
        """Compute scaled, swept, twisted airfoil profile at position y."""
        y_nodes = _cosine_y_nodes(self.half_span, self.n_panels)
        deflection = float(np.interp(y, y_nodes, self.deflection_z))

        x_norm, z_norm = blended_profile(y, self.y_stations, list(self.airfoil_sections), n_profile_points)
        c = float(np.interp(y, self.y_stations, self.chord))
        theta = float(np.interp(y, self.y_stations, self.twist))
        sweep = float(np.interp(y, self.y_stations, self.sweep_le))

        x_scaled = x_norm * c
        z_scaled = z_norm * c

        cos_t, sin_t = np.cos(theta), np.sin(theta)
        x_rot = cos_t * x_scaled + sin_t * z_scaled
        z_rot = -sin_t * x_scaled + cos_t * z_scaled

        x_rot = x_rot + y * np.tan(sweep)
        z_rot = z_rot + deflection

        return {"y": float(y), "x": x_rot, "z": z_rot}

    def get_wing_coordinates(self, n_profile_points: int = 50) -> list[dict]:
        """Return scaled and twisted airfoil coordinates at each control point."""
        return [self._section_at_y(float(y), n_profile_points) for y in self.y_cp]

    def plot_wing_3d(
        self,
        n_profile_points: int = 50,
        alpha_deg: float = 0.0,
        colorscale: str = "Blues",
        show_sections: bool = True,
        show_panels: bool = False,
        save_dir: str | None = None,
        show_grid=False,
        pov: str | None = None,
    ):
        """Render a 3D visualization of the wing using Plotly."""
        import plotly.graph_objects as go

        _POV = {
            "top": dict(eye=dict(x=0, y=0, z=+1), up=dict(x=-1, y=0, z=0)),
            "bottom": dict(eye=dict(x=0, y=0, z=-1), up=dict(x=-1, y=0, z=0)),
            "front": dict(eye=dict(x=-1, y=0, z=0), up=dict(x=0, y=0, z=1)),
            "rear": dict(eye=dict(x=+1.35, y=0, z=0), up=dict(x=0, y=0, z=1)),
            "left": dict(eye=dict(x=0, y=-1, z=0), up=dict(x=0, y=0, z=1)),
            "right": dict(eye=dict(x=0, y=+1, z=0), up=dict(x=0, y=0, z=1)),
            "iso": dict(eye=dict(x=1.5, y=-1.5, z=0.75), up=dict(x=0, y=0, z=1)),
        }
        if pov is not None and pov not in _POV:
            raise ValueError(f"pov must be one of {list(_POV)} or None, got {pov!r}")

        def mirror(sections):
            return [{"y": -s["y"], "x": s["x"], "z": s["z"]} for s in reversed(sections)]

        cp_sections = self.get_wing_coordinates(n_profile_points)
        all_cp = mirror(cp_sections) + cp_sections

        y_nodes = _cosine_y_nodes(self.half_span, self.n_panels)
        node_sections = [self._section_at_y(float(y), n_profile_points) for y in y_nodes]
        all_nodes = mirror(node_sections) + node_sections

        n_pts = 2 * n_profile_points
        n_sec = len(all_nodes)
        X = np.empty((n_pts, n_sec))
        Y = np.empty((n_pts, n_sec))
        Z = np.empty((n_pts, n_sec))
        for j, sec in enumerate(all_nodes):
            if alpha_deg != 0.0:
                from scipy.spatial.transform import Rotation as R

                r = R.from_euler("y", alpha_deg, degrees=True)
                sec["x"], sec["y"], sec["z"] = r.apply(np.matrix([sec["x"], [sec["y"]] * len(sec["x"]), sec["z"]]).T).T

            X[:, j] = sec["x"]
            Y[:, j] = sec["y"]
            Z[:, j] = sec["z"]

        traces = [
            go.Surface(
                x=X,
                y=Y,
                z=Z,
                colorscale=colorscale,
                showscale=False,
                lighting=dict(ambient=0.5, diffuse=0.8, specular=0.2),
                opacity=1.0,
            )
        ]

        if show_sections:
            for sec in all_cp:
                x_closed = np.append(sec["x"], sec["x"][0])
                z_closed = np.append(sec["z"], sec["z"][0])
                traces.append(
                    go.Scatter3d(
                        x=x_closed,
                        y=np.full_like(x_closed, sec["y"]),
                        z=z_closed,
                        mode="lines",
                        line=dict(color="royalblue", width=2),
                        showlegend=False,
                        hoverinfo="skip",
                    )
                )

        if show_panels:
            for sec in all_nodes:
                x_closed = np.append(sec["x"], sec["x"][0])
                z_closed = np.append(sec["z"], sec["z"][0])
                traces.append(
                    go.Scatter3d(
                        x=x_closed,
                        y=np.full_like(x_closed, sec["y"]),
                        z=z_closed,
                        mode="lines",
                        line=dict(color="crimson", width=1),
                        showlegend=False,
                        hoverinfo="skip",
                    )
                )

        all_x, all_y, all_z = X.ravel(), Y.ravel(), Z.ravel()
        max_range = max(all_x.max() - all_x.min(), all_y.max() - all_y.min(), all_z.max() - all_z.min())
        cx, cy, cz = (
            0.5 * (all_x.max() + all_x.min()),
            0.5 * (all_y.max() + all_y.min()),
            0.5 * (all_z.max() + all_z.min()),
        )
        half = 0.5 * max_range

        camera = _POV[pov] if pov is not None else None
        dragmode = False if pov is not None else "orbit"

        scene = dict(
            xaxis=dict(
                title="x (m)" if show_grid else "",
                range=[cx - half, cx + half],
                showbackground=False,
                showgrid=show_grid,
                showticklabels=show_grid,
                gridcolor="grey",
            ),
            yaxis=dict(
                title="y (m)" if show_grid else "",
                range=[cy - half, cy + half],
                showbackground=False,
                showgrid=show_grid,
                showticklabels=show_grid,
                gridcolor="grey",
            ),
            zaxis=dict(
                title="z (m)" if show_grid else "",
                range=[cz - half, cz + half],
                showbackground=False,
                showgrid=show_grid,
                showticklabels=show_grid,
                gridcolor="grey",
            ),
            aspectmode="cube",
            bgcolor="rgba(0,0,0,0)",
        )
        if camera is not None:
            scene["camera"] = camera

        fig = go.Figure(data=traces)
        fig.update_layout(
            scene=scene,
            paper_bgcolor="rgba(0,0,0,0)",
            plot_bgcolor="rgba(0,0,0,0)",
            margin=dict(l=0, r=0, t=0, b=0),
            dragmode=dragmode,
        )

        if save_dir:
            fig.write_html(save_dir)

        return fig

__post_init__()

Compute derived arrays and lock all arrays as read-only.

Source code in src/lightaero/geometry/wing.py

def __post_init__(self) -> None:
    """Compute derived arrays and lock all arrays as read-only."""
    y_cp = _cosine_y_cp(self.half_span, self.n_panels)

    if self.deflection_z is None:
        deflection_z = np.zeros(self.n_panels + 1)
    else:
        deflection_z = np.asarray(self.deflection_z, dtype=float)
        if deflection_z.shape != (self.n_panels + 1,):
            raise ValueError(f"deflection_z must have shape ({self.n_panels + 1},), got {deflection_z.shape}")

    chord_cp = np.interp(y_cp, self.y_stations, self.chord)
    twist_cp = np.interp(y_cp, self.y_stations, self.twist)

    x_mid = np.array([0.5])
    camber_cp = np.array(
        [
            float(
                interpolate_camber(
                    y,
                    self.y_stations,
                    list(self.airfoil_sections),
                    x_mid,
                )[0]
            )
            for y in y_cp
        ]
    )

    for arr_name, arr in [
        ("y_stations", self.y_stations),
        ("chord", self.chord),
        ("sweep_le", self.sweep_le),
        ("taper", self.taper),
        ("twist", self.twist),
    ]:
        arr_ro = np.array(arr)
        arr_ro.flags.writeable = False
        object.__setattr__(self, arr_name, arr_ro)

    for arr_name, arr in [
        ("y_cp", y_cp),
        ("chord_cp", chord_cp),
        ("twist_cp", twist_cp),
        ("camber_cp", camber_cp),
        ("deflection_z", deflection_z),
    ]:
        arr.flags.writeable = False
        object.__setattr__(self, arr_name, arr)

get_wing_coordinates(n_profile_points=50)

Return scaled and twisted airfoil coordinates at each control point.

Source code in src/lightaero/geometry/wing.py

def get_wing_coordinates(self, n_profile_points: int = 50) -> list[dict]:
    """Return scaled and twisted airfoil coordinates at each control point."""
    return [self._section_at_y(float(y), n_profile_points) for y in self.y_cp]

plot_wing_3d(n_profile_points=50, alpha_deg=0.0, colorscale='Blues', show_sections=True, show_panels=False, save_dir=None, show_grid=False, pov=None)

Render a 3D visualization of the wing using Plotly.

Source code in src/lightaero/geometry/wing.py

def plot_wing_3d(
    self,
    n_profile_points: int = 50,
    alpha_deg: float = 0.0,
    colorscale: str = "Blues",
    show_sections: bool = True,
    show_panels: bool = False,
    save_dir: str | None = None,
    show_grid=False,
    pov: str | None = None,
):
    """Render a 3D visualization of the wing using Plotly."""
    import plotly.graph_objects as go

    _POV = {
        "top": dict(eye=dict(x=0, y=0, z=+1), up=dict(x=-1, y=0, z=0)),
        "bottom": dict(eye=dict(x=0, y=0, z=-1), up=dict(x=-1, y=0, z=0)),
        "front": dict(eye=dict(x=-1, y=0, z=0), up=dict(x=0, y=0, z=1)),
        "rear": dict(eye=dict(x=+1.35, y=0, z=0), up=dict(x=0, y=0, z=1)),
        "left": dict(eye=dict(x=0, y=-1, z=0), up=dict(x=0, y=0, z=1)),
        "right": dict(eye=dict(x=0, y=+1, z=0), up=dict(x=0, y=0, z=1)),
        "iso": dict(eye=dict(x=1.5, y=-1.5, z=0.75), up=dict(x=0, y=0, z=1)),
    }
    if pov is not None and pov not in _POV:
        raise ValueError(f"pov must be one of {list(_POV)} or None, got {pov!r}")

    def mirror(sections):
        return [{"y": -s["y"], "x": s["x"], "z": s["z"]} for s in reversed(sections)]

    cp_sections = self.get_wing_coordinates(n_profile_points)
    all_cp = mirror(cp_sections) + cp_sections

    y_nodes = _cosine_y_nodes(self.half_span, self.n_panels)
    node_sections = [self._section_at_y(float(y), n_profile_points) for y in y_nodes]
    all_nodes = mirror(node_sections) + node_sections

    n_pts = 2 * n_profile_points
    n_sec = len(all_nodes)
    X = np.empty((n_pts, n_sec))
    Y = np.empty((n_pts, n_sec))
    Z = np.empty((n_pts, n_sec))
    for j, sec in enumerate(all_nodes):
        if alpha_deg != 0.0:
            from scipy.spatial.transform import Rotation as R

            r = R.from_euler("y", alpha_deg, degrees=True)
            sec["x"], sec["y"], sec["z"] = r.apply(np.matrix([sec["x"], [sec["y"]] * len(sec["x"]), sec["z"]]).T).T

        X[:, j] = sec["x"]
        Y[:, j] = sec["y"]
        Z[:, j] = sec["z"]

    traces = [
        go.Surface(
            x=X,
            y=Y,
            z=Z,
            colorscale=colorscale,
            showscale=False,
            lighting=dict(ambient=0.5, diffuse=0.8, specular=0.2),
            opacity=1.0,
        )
    ]

    if show_sections:
        for sec in all_cp:
            x_closed = np.append(sec["x"], sec["x"][0])
            z_closed = np.append(sec["z"], sec["z"][0])
            traces.append(
                go.Scatter3d(
                    x=x_closed,
                    y=np.full_like(x_closed, sec["y"]),
                    z=z_closed,
                    mode="lines",
                    line=dict(color="royalblue", width=2),
                    showlegend=False,
                    hoverinfo="skip",
                )
            )

    if show_panels:
        for sec in all_nodes:
            x_closed = np.append(sec["x"], sec["x"][0])
            z_closed = np.append(sec["z"], sec["z"][0])
            traces.append(
                go.Scatter3d(
                    x=x_closed,
                    y=np.full_like(x_closed, sec["y"]),
                    z=z_closed,
                    mode="lines",
                    line=dict(color="crimson", width=1),
                    showlegend=False,
                    hoverinfo="skip",
                )
            )

    all_x, all_y, all_z = X.ravel(), Y.ravel(), Z.ravel()
    max_range = max(all_x.max() - all_x.min(), all_y.max() - all_y.min(), all_z.max() - all_z.min())
    cx, cy, cz = (
        0.5 * (all_x.max() + all_x.min()),
        0.5 * (all_y.max() + all_y.min()),
        0.5 * (all_z.max() + all_z.min()),
    )
    half = 0.5 * max_range

    camera = _POV[pov] if pov is not None else None
    dragmode = False if pov is not None else "orbit"

    scene = dict(
        xaxis=dict(
            title="x (m)" if show_grid else "",
            range=[cx - half, cx + half],
            showbackground=False,
            showgrid=show_grid,
            showticklabels=show_grid,
            gridcolor="grey",
        ),
        yaxis=dict(
            title="y (m)" if show_grid else "",
            range=[cy - half, cy + half],
            showbackground=False,
            showgrid=show_grid,
            showticklabels=show_grid,
            gridcolor="grey",
        ),
        zaxis=dict(
            title="z (m)" if show_grid else "",
            range=[cz - half, cz + half],
            showbackground=False,
            showgrid=show_grid,
            showticklabels=show_grid,
            gridcolor="grey",
        ),
        aspectmode="cube",
        bgcolor="rgba(0,0,0,0)",
    )
    if camera is not None:
        scene["camera"] = camera

    fig = go.Figure(data=traces)
    fig.update_layout(
        scene=scene,
        paper_bgcolor="rgba(0,0,0,0)",
        plot_bgcolor="rgba(0,0,0,0)",
        margin=dict(l=0, r=0, t=0, b=0),
        dragmode=dragmode,
    )

    if save_dir:
        fig.write_html(save_dir)

    return fig

build_wing_geometry(half_span, n_panels, y_stations, chord, sweep_le, taper, twist, airfoil_sections, deflection_z=None)

Factory for WingGeometry with input validation.

See: docs/theory/index.md for geometric constraints and validation rules.

Source code in src/lightaero/geometry/wing.py

def build_wing_geometry(
    half_span: float,
    n_panels: int,
    y_stations: np.ndarray,
    chord: np.ndarray,
    sweep_le: np.ndarray,
    taper: np.ndarray,
    twist: np.ndarray,
    airfoil_sections: tuple,
    deflection_z: np.ndarray = None,
) -> WingGeometry:
    """Factory for WingGeometry with input validation.

    See: docs/theory/index.md for geometric constraints and validation rules.
    """
    y_stations = np.asarray(y_stations, dtype=float)
    chord = np.asarray(chord, dtype=float)
    sweep_le = np.asarray(sweep_le, dtype=float)
    taper = np.asarray(taper, dtype=float)
    twist = np.asarray(twist, dtype=float)

    if half_span <= 0.0:
        raise ValueError(f"half_span must be positive, got {half_span}")
    if n_panels < 4:
        raise ValueError(f"n_panels must be >= 4, got {n_panels}")
    n_sections = len(y_stations)
    if n_sections < 2:
        raise ValueError("y_stations must have at least 2 entries (root and tip)")
    if not np.isclose(y_stations[0], 0.0):
        raise ValueError(f"y_stations[0] must be 0.0 (root), got {y_stations[0]}")
    if not np.isclose(y_stations[-1], half_span):
        raise ValueError(f"y_stations[-1] must equal half_span={half_span}, got {y_stations[-1]}")
    if not np.all(np.diff(y_stations) > 0):
        raise ValueError("y_stations must be strictly monotonically increasing")
    for arr_name, arr in [("chord", chord), ("sweep_le", sweep_le), ("taper", taper), ("twist", twist)]:
        if arr.shape != (n_sections,):
            raise ValueError(f"{arr_name} must have shape ({n_sections},), got {arr.shape}")
    if np.any(chord <= 0.0):
        raise ValueError("All chord values must be positive")
    if len(airfoil_sections) != n_sections:
        raise ValueError(f"airfoil_sections must have {n_sections} entries, got {len(airfoil_sections)}")

    return WingGeometry(
        half_span=float(half_span),
        n_panels=int(n_panels),
        y_stations=y_stations,
        chord=chord,
        sweep_le=sweep_le,
        taper=taper,
        twist=twist,
        airfoil_sections=tuple(airfoil_sections),
        deflection_z=deflection_z,
    )

registry

Fidelity plugin registry subpackage for lightaero.

Exposes DISCIPLINE_REGISTRY, register, get_discipline, and DisciplineBase so that discipline modules can import from lightaero.registry directly:

from lightaero.registry import register, DisciplineBase, DISCIPLINE_REGISTRY

DisciplineBase

Bases: ABC

Abstract base class for all lightaero discipline implementations.

Enforces the black-box callable interface: discipline(inputs: dict) -> TypedDict

See: docs/theory/index.md

Source code in src/lightaero/registry/registry.py

class DisciplineBase(ABC):
    """Abstract base class for all lightaero discipline implementations.

    Enforces the black-box callable interface: discipline(inputs: dict) -> TypedDict

    See: docs/theory/index.md
    """

    @abstractmethod
    def __call__(self, **inputs: Any) -> Any:
        """Evaluate the discipline at the given inputs.

        Args:
            **inputs: Named physical quantities in SI units.

        Returns:
            Any: Results of the evaluation (e.g. AeroOutput, StructuralOutput).
        """
        ...

__call__(**inputs) abstractmethod

Evaluate the discipline at the given inputs.

Parameters:

Name	Type	Description	Default
**inputs	Any	Named physical quantities in SI units.	{}

Returns:

Name	Type	Description
Any	Any	Results of the evaluation (e.g. AeroOutput, StructuralOutput).

Source code in src/lightaero/registry/registry.py

@abstractmethod
def __call__(self, **inputs: Any) -> Any:
    """Evaluate the discipline at the given inputs.

    Args:
        **inputs: Named physical quantities in SI units.

    Returns:
        Any: Results of the evaluation (e.g. AeroOutput, StructuralOutput).
    """
    ...

get_discipline(family, key)

Retrieve a discipline class from the registry.

Parameters:

Name	Type	Description	Default
family	str	Discipline family (e.g. 'aero').	required
key	str	Implementation key (e.g. 'vlm').	required

Returns:

Type	Description
type	The registered class.

Raises:

Type	Description
KeyError	If family or key is not found.

Source code in src/lightaero/registry/registry.py

def get_discipline(family: str, key: str) -> type:
    """Retrieve a discipline class from the registry.

    Args:
        family: Discipline family (e.g. 'aero').
        key: Implementation key (e.g. 'vlm').

    Returns:
        The registered class.

    Raises:
        KeyError: If family or key is not found.
    """
    if family not in DISCIPLINE_REGISTRY:
        available = list(DISCIPLINE_REGISTRY.keys())
        raise KeyError(
            f"Unknown discipline family: {family!r}. "
            f"Available families: {available}. "
            f"Ensure the discipline module has been imported."
        )
    if key not in DISCIPLINE_REGISTRY[family]:
        available = list(DISCIPLINE_REGISTRY[family].keys())
        raise KeyError(f"Unknown key {key!r} for family {family!r}. Available keys: {available}.")
    return DISCIPLINE_REGISTRY[family][key]

register(family, key)

Decorator factory for two-level namespace registration.

Registers the decorated class into DISCIPLINE_REGISTRY[family][key].

Parameters:

Name	Type	Description	Default
family	str	Discipline family string (e.g. 'aero', 'structures').	required
key	str	Fidelity/implementation key (e.g. 'vlm', 'low_fi').	required

Returns:

Type	Description
Callable[[type], type]	Callable[[type], type]: Decorator that registers the class and returns it unchanged.

Source code in src/lightaero/registry/registry.py

def register(family: str, key: str) -> Callable[[type], type]:
    """Decorator factory for two-level namespace registration.

    Registers the decorated class into `DISCIPLINE_REGISTRY[family][key]`.

    Args:
        family: Discipline family string (e.g. 'aero', 'structures').
        key: Fidelity/implementation key (e.g. 'vlm', 'low_fi').

    Returns:
        Callable[[type], type]: Decorator that registers the class and returns it unchanged.
    """

    def decorator(cls: type) -> type:
        if family not in DISCIPLINE_REGISTRY:
            DISCIPLINE_REGISTRY[family] = {}
        DISCIPLINE_REGISTRY[family][key] = cls
        return cls

    return decorator

registry

Fidelity plugin registry for lightaero disciplines.

Provides a two-level namespace registry (family -> key -> class) so that discipline implementations can be swapped by string key without modifying coupling logic.

See: docs/theory/index.md

DisciplineBase

Bases: ABC

Abstract base class for all lightaero discipline implementations.

Enforces the black-box callable interface: discipline(inputs: dict) -> TypedDict

See: docs/theory/index.md

Source code in src/lightaero/registry/registry.py

class DisciplineBase(ABC):
    """Abstract base class for all lightaero discipline implementations.

    Enforces the black-box callable interface: discipline(inputs: dict) -> TypedDict

    See: docs/theory/index.md
    """

    @abstractmethod
    def __call__(self, **inputs: Any) -> Any:
        """Evaluate the discipline at the given inputs.

        Args:
            **inputs: Named physical quantities in SI units.

        Returns:
            Any: Results of the evaluation (e.g. AeroOutput, StructuralOutput).
        """
        ...

__call__(**inputs) abstractmethod

Evaluate the discipline at the given inputs.

Parameters:

Name	Type	Description	Default
**inputs	Any	Named physical quantities in SI units.	{}

Returns:

Name	Type	Description
Any	Any	Results of the evaluation (e.g. AeroOutput, StructuralOutput).

Source code in src/lightaero/registry/registry.py

@abstractmethod
def __call__(self, **inputs: Any) -> Any:
    """Evaluate the discipline at the given inputs.

    Args:
        **inputs: Named physical quantities in SI units.

    Returns:
        Any: Results of the evaluation (e.g. AeroOutput, StructuralOutput).
    """
    ...

get_discipline(family, key)

Retrieve a discipline class from the registry.

Parameters:

Name	Type	Description	Default
family	str	Discipline family (e.g. 'aero').	required
key	str	Implementation key (e.g. 'vlm').	required

Returns:

Type	Description
type	The registered class.

Raises:

Type	Description
KeyError	If family or key is not found.

Source code in src/lightaero/registry/registry.py

def get_discipline(family: str, key: str) -> type:
    """Retrieve a discipline class from the registry.

    Args:
        family: Discipline family (e.g. 'aero').
        key: Implementation key (e.g. 'vlm').

    Returns:
        The registered class.

    Raises:
        KeyError: If family or key is not found.
    """
    if family not in DISCIPLINE_REGISTRY:
        available = list(DISCIPLINE_REGISTRY.keys())
        raise KeyError(
            f"Unknown discipline family: {family!r}. "
            f"Available families: {available}. "
            f"Ensure the discipline module has been imported."
        )
    if key not in DISCIPLINE_REGISTRY[family]:
        available = list(DISCIPLINE_REGISTRY[family].keys())
        raise KeyError(f"Unknown key {key!r} for family {family!r}. Available keys: {available}.")
    return DISCIPLINE_REGISTRY[family][key]

register(family, key)

Decorator factory for two-level namespace registration.

Registers the decorated class into DISCIPLINE_REGISTRY[family][key].

Parameters:

Name	Type	Description	Default
family	str	Discipline family string (e.g. 'aero', 'structures').	required
key	str	Fidelity/implementation key (e.g. 'vlm', 'low_fi').	required

Returns:

Type	Description
Callable[[type], type]	Callable[[type], type]: Decorator that registers the class and returns it unchanged.

Source code in src/lightaero/registry/registry.py

def register(family: str, key: str) -> Callable[[type], type]:
    """Decorator factory for two-level namespace registration.

    Registers the decorated class into `DISCIPLINE_REGISTRY[family][key]`.

    Args:
        family: Discipline family string (e.g. 'aero', 'structures').
        key: Fidelity/implementation key (e.g. 'vlm', 'low_fi').

    Returns:
        Callable[[type], type]: Decorator that registers the class and returns it unchanged.
    """

    def decorator(cls: type) -> type:
        if family not in DISCIPLINE_REGISTRY:
            DISCIPLINE_REGISTRY[family] = {}
        DISCIPLINE_REGISTRY[family][key] = cls
        return cls

    return decorator

schemas

types

Output schemas and companion runtime spec dicts.

These schemas define the inter-discipline data contract for lightaero.

See: docs/theory/index.md

validation

Runtime validator for discipline output dicts.

See: docs/theory/index.md

check_regime_validity(mach, aoa_rad)

Check if flight condition is within standard VLM validity bounds.

Issues UserWarning if: - Mach > 0.3 (incompressible limit) - |AoA| > 10 degrees (linear lift limit)

Parameters:

Name	Type	Description	Default
mach	float	Mach number.	required
aoa_rad	float	Angle of attack in radians.	required

Source code in src/lightaero/schemas/validation.py

def check_regime_validity(mach: float, aoa_rad: float) -> None:
    """Check if flight condition is within standard VLM validity bounds.

    Issues UserWarning if:
    - Mach > 0.3 (incompressible limit)
    - |AoA| > 10 degrees (linear lift limit)

    Args:
        mach: Mach number.
        aoa_rad: Angle of attack in radians.
    """
    if mach > 0.3:
        msg = f"Mach {mach:.2f} exceeds incompressible limit (0.3). Results may be inaccurate."
        warnings.warn(msg, UserWarning, stacklevel=2)
        logger.warning(msg)

    aoa_deg = math.degrees(aoa_rad)
    if abs(aoa_deg) > 10.0:
        msg = f"AoA {aoa_deg:.1f} deg exceeds typical VLM linear limit (10 deg). Results may be inaccurate."
        warnings.warn(msg, UserWarning, stacklevel=2)
        logger.warning(msg)

validate_discipline_output(output, spec)

Runtime validation of a discipline output dict against a spec.

Checks key presence, type correctness, array shape, and plausibility bounds.

Parameters:

Name	Type	Description	Default
output	Any	Object (dataclass) returned by a discipline's call method.	required
spec	dict	Companion spec dict (e.g. AERO_OUTPUT_SPEC).	required

Raises:

Type	Description
ValueError	Missing required key, wrong array shape, or out-of-bounds value.
TypeError	Wrong Python type for a field.

Source code in src/lightaero/schemas/validation.py

def validate_discipline_output(output: Any, spec: dict) -> None:
    """Runtime validation of a discipline output dict against a spec.

    Checks key presence, type correctness, array shape, and plausibility bounds.

    Args:
        output: Object (dataclass) returned by a discipline's __call__ method.
        spec: Companion spec dict (e.g. AERO_OUTPUT_SPEC).

    Raises:
        ValueError: Missing required key, wrong array shape, or out-of-bounds value.
        TypeError: Wrong Python type for a field.
    """
    # --- 1. Key presence ---
    missing = set(spec.keys()) - set([f.name for f in fields(output)])
    if missing:
        # Report all missing keys in alphabetical order for deterministic messages
        missing_str = ", ".join(sorted(missing))
        raise ValueError(f"Missing required output keys: {missing_str}")

    # --- 2. Type, shape, and bounds per field ---
    for key, rules in spec.items():
        val: Any = getattr(output, key)
        expected_type = rules["type"]

        # Type check - use isinstance
        if not isinstance(val, expected_type):
            raise TypeError(f"Output field '{key}': expected {expected_type.__name__}, got {type(val).__name__}")

        # Array shape check (applies only when type is np.ndarray)
        if expected_type is np.ndarray and "shape_dim" in rules:
            expected_ndim: int = rules["shape_dim"]
            if val.ndim != expected_ndim:
                raise ValueError(
                    f"Output field '{key}': expected {expected_ndim}D array, "
                    f"got {val.ndim}D array with shape {val.shape}"
                )

        # Plausibility bounds check (applies only to scalar types with bounds)
        if "bounds" in rules and expected_type is not np.ndarray:
            lo, hi = rules["bounds"]
            scalar_val = float(val)
            if lo is not None and scalar_val < lo:
                raise ValueError(
                    f"Output field '{key}' = {scalar_val} is below minimum "
                    f"plausible value {lo}. Check SI units (possible unit mismatch)."
                )
            if hi is not None and scalar_val > hi:
                raise ValueError(
                    f"Output field '{key}' = {scalar_val} is above maximum "
                    f"plausible value {hi}. Check SI units (possible unit mismatch)."
                )

    return None

structures

Structural discipline subpackage for lightaero.

See: docs/theory/structures.md

beam

Euler-Bernoulli beam FEM kernel.

See: docs/theory/structures.md

beam_section_properties(y_nodes, chord_nodes, tc_ratio_nodes, E_Pa=73100000000.0, G_Pa=28000000000.0, t_sk_array=None)

Thin-walled wingbox section properties via numerical integration.

Parameters:

Name	Type	Description	Default
y_nodes	ndarray	Spanwise node positions (m).	required
chord_nodes	ndarray	Chord at each node (m).	required
tc_ratio_nodes	ndarray	t/c ratio at each node.	required
E_Pa	float	Young's modulus (Pa).	73100000000.0
G_Pa	float	Shear modulus (Pa).	28000000000.0
t_sk_array	ndarray | None	Optional skin thickness at each node (m).	None

Returns:

Type	Description
(EI_xx, EI_zz, GJ, h_max, w_box, I_xx, I_zz, A_e, x_sc)	each shape (n,).

Source code in src/lightaero/structures/beam.py

def beam_section_properties(
    y_nodes: np.ndarray,
    chord_nodes: np.ndarray,
    tc_ratio_nodes: np.ndarray,
    E_Pa: float = 73.1e9,
    G_Pa: float = 28.0e9,
    t_sk_array: np.ndarray | None = None,
) -> tuple:
    """Thin-walled wingbox section properties via numerical integration.

    Args:
        y_nodes: Spanwise node positions (m).
        chord_nodes: Chord at each node (m).
        tc_ratio_nodes: t/c ratio at each node.
        E_Pa: Young's modulus (Pa).
        G_Pa: Shear modulus (Pa).
        t_sk_array: Optional skin thickness at each node (m).

    Returns:
        (EI_xx, EI_zz, GJ, h_max, w_box, I_xx, I_zz, A_e, x_sc): each shape (n,).
    """
    n_nodes = len(chord_nodes)
    EI_xx = np.zeros(n_nodes)
    EI_zz = np.zeros(n_nodes)
    GJ = np.zeros(n_nodes)
    h_max = np.zeros(n_nodes)
    w_box = np.zeros(n_nodes)
    I_xx_out = np.zeros(n_nodes)
    I_zz_out = np.zeros(n_nodes)
    A_e_out = np.zeros(n_nodes)
    x_sc_out = np.zeros(n_nodes)

    x_norm = np.linspace(_X_START, _X_END, _N_PTS)

    for i in range(n_nodes):
        c = chord_nodes[i]
        tc = tc_ratio_nodes[i]
        t_sk = t_sk_array[i] if t_sk_array is not None else (tc * c) / 100.0

        y_norm = (
            5.0
            * tc
            * (
                0.2969 * np.sqrt(x_norm)
                - 0.1260 * x_norm
                - 0.3516 * x_norm**2
                + 0.2843 * x_norm**3
                - 0.1015 * x_norm**4
            )
        )

        x = x_norm * c
        y = y_norm * c
        h = 2.0 * y

        h_f = h[0]
        h_r = h[-1]

        h_max[i] = np.max(h)
        w_box[i] = c * (_X_END - _X_START)

        dx = np.diff(x)
        dy = np.diff(y)
        ds = np.sqrt(dx**2 + dy**2)

        y_mid = 0.5 * (y[:-1] + y[1:])
        x_mid = 0.5 * (x[:-1] + x[1:])

        I_xx_upper = np.sum(t_sk * y_mid**2 * ds)
        I_xx_fs = (1.0 / 12.0) * t_sk * h_f**3
        I_xx_rs = (1.0 / 12.0) * t_sk * h_r**3
        I_xx = 2.0 * I_xx_upper + I_xx_fs + I_xx_rs

        A_skin = np.sum(t_sk * ds)
        A_fs = t_sk * h_f
        A_rs = t_sk * h_r
        A_tot = 2.0 * A_skin + A_fs + A_rs

        x_centroid = (2.0 * np.sum(x_mid * t_sk * ds) + x[0] * A_fs + x[-1] * A_rs) / A_tot

        I_zz_upper = np.sum(t_sk * (x_mid - x_centroid) ** 2 * ds)
        I_zz_fs = A_fs * (x[0] - x_centroid) ** 2
        I_zz_rs = A_rs * (x[-1] - x_centroid) ** 2
        I_zz = 2.0 * I_zz_upper + I_zz_fs + I_zz_rs

        A_e = np.trapezoid(h, x)
        perimeter = 2.0 * np.sum(ds) + h_f + h_r
        J = (4.0 * A_e**2 * t_sk) / perimeter

        EI_xx[i] = E_Pa * I_xx
        EI_zz[i] = E_Pa * I_zz
        GJ[i] = G_Pa * J
        I_xx_out[i] = I_xx
        I_zz_out[i] = I_zz
        A_e_out[i] = A_e
        x_sc_out[i] = _compute_shear_center_x(c, tc, t_sk, I_xx, x_norm)

    return EI_xx, EI_zz, GJ, h_max, w_box, I_xx_out, I_zz_out, A_e_out, x_sc_out

build_bending_loads(y_nodes, q_nodes)

Build FEM load vector for bending.

Parameters:

Name	Type	Description	Default
y_nodes	ndarray	Spanwise node positions (m).	required
q_nodes	ndarray	Distributed transverse load (N/m).	required

Returns:

Name	Type	Description
F	ndarray	Load vector (N and N·m).

Source code in src/lightaero/structures/beam.py

def build_bending_loads(y_nodes: np.ndarray, q_nodes: np.ndarray) -> np.ndarray:
    """Build FEM load vector for bending.

    Args:
        y_nodes: Spanwise node positions (m).
        q_nodes: Distributed transverse load (N/m).

    Returns:
        F: Load vector (N and N·m).
    """
    n_nodes = len(y_nodes)
    F = np.zeros(2 * n_nodes)
    for i in range(n_nodes - 1):
        L_e = y_nodes[i + 1] - y_nodes[i]
        q_avg = 0.5 * (q_nodes[i] + q_nodes[i + 1])
        F[2 * i] += q_avg * L_e / 2.0
        F[2 * i + 1] += q_avg * L_e**2 / 12.0
        F[2 * i + 2] += q_avg * L_e / 2.0
        F[2 * i + 3] -= q_avg * L_e**2 / 12.0
    return F

build_torsion_loads(y_nodes, t_nodes)

Build FEM load vector for torsion.

Parameters:

Name	Type	Description	Default
y_nodes	ndarray	Spanwise node positions (m).	required
t_nodes	ndarray	Distributed torque (N·m/m).	required

Returns:

Name	Type	Description
T	ndarray	Load vector (N·m).

Source code in src/lightaero/structures/beam.py

def build_torsion_loads(y_nodes: np.ndarray, t_nodes: np.ndarray) -> np.ndarray:
    """Build FEM load vector for torsion.

    Args:
        y_nodes: Spanwise node positions (m).
        t_nodes: Distributed torque (N·m/m).

    Returns:
        T: Load vector (N·m).
    """
    n_nodes = len(y_nodes)
    T = np.zeros(n_nodes)
    for i in range(n_nodes - 1):
        L_e = y_nodes[i + 1] - y_nodes[i]
        t_avg = 0.5 * (t_nodes[i] + t_nodes[i + 1])
        T[i] += t_avg * L_e / 2.0
        T[i + 1] += t_avg * L_e / 2.0
    return T

compute_combined_stress(y_nodes, w_full, u_full, EI_xx, EI_zz, h, w, I_xx, I_zz, phi_full=None, GJ=None, A_e=None, t_sk=None)

Von Mises stress at each spanwise node.

See: docs/theory/structures.md

Parameters:

Name	Type	Description	Default
y_nodes	ndarray	Spanwise node positions (m).	required
w_full	ndarray	Bending DOF vector (vertical).	required
u_full	ndarray	Bending DOF vector (lateral).	required
EI_xx	ndarray	Bending stiffness array about xx.	required
EI_zz	ndarray	Bending stiffness array about zz.	required
h	ndarray	Section height at each node (m).	required
w	ndarray	Section width at each node (m).	required
I_xx	ndarray	Second moment of area about xx (m^4).	required
I_zz	ndarray	Second moment of area about zz (m^4).	required
phi_full	ndarray | None	Torsion DOF vector. Optional.	None
GJ	ndarray | None	Torsional stiffness array. Optional.	None
A_e	ndarray | None	Enclosed wingbox area at each node (m²). Optional.	None
t_sk	ndarray | None	Skin thickness at each node (m). Optional.	None

Returns:

Type	Description
ndarray	np.ndarray: Von Mises stress at each node (Pa).

Source code in src/lightaero/structures/beam.py

def compute_combined_stress(
    y_nodes: np.ndarray,
    w_full: np.ndarray,
    u_full: np.ndarray,
    EI_xx: np.ndarray,
    EI_zz: np.ndarray,
    h: np.ndarray,
    w: np.ndarray,
    I_xx: np.ndarray,
    I_zz: np.ndarray,
    phi_full: np.ndarray | None = None,
    GJ: np.ndarray | None = None,
    A_e: np.ndarray | None = None,
    t_sk: np.ndarray | None = None,
) -> np.ndarray:
    """Von Mises stress at each spanwise node.

    See: docs/theory/structures.md

    Args:
        y_nodes (np.ndarray): Spanwise node positions (m).
        w_full (np.ndarray): Bending DOF vector (vertical).
        u_full (np.ndarray): Bending DOF vector (lateral).
        EI_xx (np.ndarray): Bending stiffness array about xx.
        EI_zz (np.ndarray): Bending stiffness array about zz.
        h (np.ndarray): Section height at each node (m).
        w (np.ndarray): Section width at each node (m).
        I_xx (np.ndarray): Second moment of area about xx (m^4).
        I_zz (np.ndarray): Second moment of area about zz (m^4).
        phi_full (np.ndarray | None): Torsion DOF vector. Optional.
        GJ (np.ndarray | None): Torsional stiffness array. Optional.
        A_e (np.ndarray | None): Enclosed wingbox area at each node (m²). Optional.
        t_sk (np.ndarray | None): Skin thickness at each node (m). Optional.

    Returns:
        np.ndarray: Von Mises stress at each node (Pa).
    """
    n_nodes = len(y_nodes)
    stress = np.zeros(n_nodes)
    use_torsion = phi_full is not None and GJ is not None and A_e is not None and t_sk is not None

    for i in range(n_nodes - 1):
        L_e = y_nodes[i + 1] - y_nodes[i]

        EI_xx_e = 0.5 * (EI_xx[i] + EI_xx[i + 1])
        w_A, th_A, w_B, th_B = w_full[2 * i : 2 * i + 4]
        kappa_x = (6.0 * (w_A - w_B) / L_e**2) + (4.0 * th_A + 2.0 * th_B) / L_e
        M_x = EI_xx_e * kappa_x

        EI_zz_e = 0.5 * (EI_zz[i] + EI_zz[i + 1])
        u_A, thz_A, u_B, thz_B = u_full[2 * i : 2 * i + 4]
        kappa_z = (6.0 * (u_A - u_B) / L_e**2) + (4.0 * thz_A + 2.0 * thz_B) / L_e
        M_z = EI_zz_e * kappa_z

        sigma_b = abs(M_x) * (h[i] / 2.0) / max(I_xx[i], 1e-30) + abs(M_z) * (w[i] / 2.0) / max(I_zz[i], 1e-30)

        if use_torsion:
            GJ_e = 0.5 * (GJ[i] + GJ[i + 1])
            T_elem = GJ_e * (phi_full[i + 1] - phi_full[i]) / L_e
            denom = 2.0 * max(A_e[i], 1e-30) * max(t_sk[i], 1e-9)
            tau = abs(T_elem) / denom
            sigma_vm = (sigma_b**2 + 3.0 * tau**2) ** 0.5
        else:
            sigma_vm = sigma_b

        if sigma_vm > stress[i]:
            stress[i] = sigma_vm

    return stress

size_skin_thickness(y_nodes, chord_nodes, tc_ratio_nodes, q_z_nodes, sigma_yield=324000000.0, SF=1.5, t_sk_min=0.001)

Size wing skin thickness to satisfy bending stress allowable.

Parameters:

Name	Type	Description	Default
y_nodes	ndarray	Spanwise node positions (m).	required
chord_nodes	ndarray	Chord length at each node (m).	required
tc_ratio_nodes	ndarray	t/c ratio at each node.	required
q_z_nodes	ndarray	Distributed vertical load (N/m) at each node.	required
sigma_yield	float	Tensile yield stress (Pa).	324000000.0
SF	float	Safety factor.	1.5
t_sk_min	float	Minimum skin thickness (m).	0.001

Returns:

Name	Type	Description
t_sk	ndarray	Sized skin thickness at each node (m).

Source code in src/lightaero/structures/beam.py

def size_skin_thickness(
    y_nodes: np.ndarray,
    chord_nodes: np.ndarray,
    tc_ratio_nodes: np.ndarray,
    q_z_nodes: np.ndarray,
    sigma_yield: float = 324e6,
    SF: float = 1.5,
    t_sk_min: float = 0.001,
) -> np.ndarray:
    """Size wing skin thickness to satisfy bending stress allowable.

    Args:
        y_nodes: Spanwise node positions (m).
        chord_nodes: Chord length at each node (m).
        tc_ratio_nodes: t/c ratio at each node.
        q_z_nodes: Distributed vertical load (N/m) at each node.
        sigma_yield: Tensile yield stress (Pa).
        SF: Safety factor.
        t_sk_min: Minimum skin thickness (m).

    Returns:
        t_sk: Sized skin thickness at each node (m).
    """
    sigma_allow = sigma_yield / SF
    M = _bending_moment_from_loads(y_nodes, q_z_nodes)
    x_norm = np.linspace(_X_START, _X_END, _N_PTS)

    n_nodes = len(chord_nodes)
    t_sk = np.full(n_nodes, t_sk_min)

    for i in range(n_nodes):
        c = chord_nodes[i]
        tc = tc_ratio_nodes[i]

        y_norm = (
            5.0
            * tc
            * (
                0.2969 * np.sqrt(x_norm)
                - 0.1260 * x_norm
                - 0.3516 * x_norm**2
                + 0.2843 * x_norm**3
                - 0.1015 * x_norm**4
            )
        )
        x = x_norm * c
        y = y_norm * c
        h = 2.0 * y
        h_max_i = float(np.max(h))
        h_f = float(h[0])
        h_r = float(h[-1])

        ds = np.sqrt(np.diff(x) ** 2 + np.diff(y) ** 2)
        y_mid = 0.5 * (y[:-1] + y[1:])

        I_xx_unit = 2.0 * float(np.sum(y_mid**2 * ds)) + (1.0 / 12.0) * h_f**3 + (1.0 / 12.0) * h_r**3

        if I_xx_unit > 0.0 and h_max_i > 0.0:
            t_req = abs(M[i]) * (h_max_i / 2.0) / (I_xx_unit * sigma_allow)
            t_sk[i] = max(t_sk_min, t_req)

    return t_sk

solve_bending_beam(y_nodes, EI, F)

Solve FEM beam bending equations.

Parameters:

Name	Type	Description	Default
y_nodes	ndarray	Spanwise node positions (m).	required
EI	ndarray	Bending stiffness at nodes (N·m²).	required
F	ndarray	Load vector.	required

Returns:

Name	Type	Description
d	ndarray	Displacement vector (m and rad).

Source code in src/lightaero/structures/beam.py

def solve_bending_beam(y_nodes: np.ndarray, EI: np.ndarray, F: np.ndarray) -> np.ndarray:
    """Solve FEM beam bending equations.

    Args:
        y_nodes: Spanwise node positions (m).
        EI: Bending stiffness at nodes (N·m²).
        F: Load vector.

    Returns:
        d: Displacement vector (m and rad).
    """
    n_nodes = len(y_nodes)
    K = np.zeros((2 * n_nodes, 2 * n_nodes))

    for i in range(n_nodes - 1):
        L_e = y_nodes[i + 1] - y_nodes[i]
        EI_e = 0.5 * (EI[i] + EI[i + 1])
        c = EI_e / L_e**3
        k = c * np.array(
            [
                [12.0, 6.0 * L_e, -12.0, 6.0 * L_e],
                [6.0 * L_e, 4.0 * L_e**2, -6.0 * L_e, 2.0 * L_e**2],
                [-12.0, -6.0 * L_e, 12.0, -6.0 * L_e],
                [6.0 * L_e, 2.0 * L_e**2, -6.0 * L_e, 4.0 * L_e**2],
            ]
        )
        idx = [2 * i, 2 * i + 1, 2 * i + 2, 2 * i + 3]
        for r in range(4):
            for col in range(4):
                K[idx[r], idx[col]] += k[r, col]

    K_ff = K[2:, 2:]
    F_f = F[2:]
    d_free = np.linalg.solve(K_ff, F_f)
    return np.concatenate(([0.0, 0.0], d_free))

solve_torsion_beam(y_nodes, GJ, T)

Solve FEM beam torsion equations.

Parameters:

Name	Type	Description	Default
y_nodes	ndarray	Spanwise node positions (m).	required
GJ	ndarray	Torsional stiffness at nodes (N·m²).	required
T	ndarray	Torque load vector (N·m).	required

Returns:

Name	Type	Description
phi	ndarray	Twist angle vector (rad).

Source code in src/lightaero/structures/beam.py

def solve_torsion_beam(y_nodes: np.ndarray, GJ: np.ndarray, T: np.ndarray) -> np.ndarray:
    """Solve FEM beam torsion equations.

    Args:
        y_nodes: Spanwise node positions (m).
        GJ: Torsional stiffness at nodes (N·m²).
        T: Torque load vector (N·m).

    Returns:
        phi: Twist angle vector (rad).
    """
    n_nodes = len(y_nodes)
    K = np.zeros((n_nodes, n_nodes))
    for i in range(n_nodes - 1):
        L_e = y_nodes[i + 1] - y_nodes[i]
        GJ_e = 0.5 * (GJ[i] + GJ[i + 1])
        k = (GJ_e / L_e) * np.array([[1.0, -1.0], [-1.0, 1.0]])
        K[i : i + 2, i : i + 2] += k

    K_ff = K[1:, 1:]
    T_f = T[1:]
    phi_free = np.linalg.solve(K_ff, T_f)
    return np.concatenate(([0.0], phi_free))

mass

Wing structural mass integration and OEW empirical estimation.

See: docs/theory/structures.md

estimate_oew(MTOW_kg, wing_mass_kg, S_ref_m2=None)

Estimate operating empty weight using Roskam log-linear regression.

Parameters:

Name	Type	Description	Default
MTOW_kg	float	Maximum take-off weight (kg).	required
wing_mass_kg	float	Wing structural mass (kg).	required
S_ref_m2	float | None	Reference wing area (m²).	None

Returns:

Name	Type	Description
OEW	float	Estimated operating empty weight (kg).

Source code in src/lightaero/structures/mass.py

def estimate_oew(
    MTOW_kg: float,
    wing_mass_kg: float,
    S_ref_m2: float | None = None,
) -> float:
    """Estimate operating empty weight using Roskam log-linear regression.

    Args:
        MTOW_kg: Maximum take-off weight (kg).
        wing_mass_kg: Wing structural mass (kg).
        S_ref_m2: Reference wing area (m²).

    Returns:
        OEW: Estimated operating empty weight (kg).
    """

    MTOW_safe = max(MTOW_kg, 1.0)
    oew_roskam = 10.0 ** (-0.0833 + 0.9647 * math.log10(MTOW_safe))
    oew_floor = wing_mass_kg * 1.3
    return max(oew_roskam, oew_floor)

wing_structural_mass(y_nodes, t_sk, h, w, rho_mat=2780.0, k_struct=1.5)

Integrated wing structural mass.

Parameters:

Name	Type	Description	Default
y_nodes	ndarray	Spanwise node positions (m).	required
t_sk	ndarray	Skin thickness at each node (m).	required
h	ndarray	Section height at each node (m).	required
w	ndarray	Section width at each node (m).	required
rho_mat	float	Material density (kg/m³).	2780.0
k_struct	float	Structural mass factor.	1.5

Returns:

Name	Type	Description
mass	float	Total wing mass (kg).

Source code in src/lightaero/structures/mass.py

def wing_structural_mass(
    y_nodes: np.ndarray,
    t_sk: np.ndarray,
    h: np.ndarray,
    w: np.ndarray,
    rho_mat: float = 2780.0,
    k_struct: float = 1.5,
) -> float:
    """Integrated wing structural mass.

    Args:
        y_nodes: Spanwise node positions (m).
        t_sk: Skin thickness at each node (m).
        h: Section height at each node (m).
        w: Section width at each node (m).
        rho_mat: Material density (kg/m³).
        k_struct: Structural mass factor.

    Returns:
        mass: Total wing mass (kg).
    """
    # Rectangular thin-wall perimeter approximation
    A_cross = 2.0 * t_sk * (w + h)
    m_per_span = rho_mat * A_cross * k_struct
    m_semi = np.trapezoid(m_per_span, y_nodes)
    return 2.0 * m_semi

solver

StructuralDiscipline: compose beam.py and mass.py into a single discipline call.