"""Orthogonal airfoil implementation built from camber and thickness."""
from __future__ import annotations
from collections.abc import Callable
from typing import Literal, override
import numpy as np
from buffalo_core.numeric import as_float_array, as_float_scalar
from buffalo_core.typing import FloatArray, FloatInput, FloatScalar
from scipy.optimize import root_scalar
from buffalo_wings.airfoil.internal.base import Airfoil
from buffalo_wings.airfoil.internal.parameter_validation import surface_is_upper
from buffalo_wings.airfoil.internal.runtime_common import (
LOWER_SURFACE_BRACKET,
ROOT_ABS_TOLERANCE,
UPPER_SURFACE_BRACKET,
)
from buffalo_wings.airfoil.internal.runtime_types import (
AirfoilCamberResult,
AirfoilSurface,
SurfaceMappedValues,
)
from .camber_base import Camber
from .thickness_base import Thickness
_CAMBER_TOLERANCE = 1e-7
_ROOT_SOLVER_MAXITER = 200
_ROOT_SOLVER_TOLERANCE = 1e-12
_BREAKPOINT_ZERO_TOLERANCE = 1e-15
[docs]
class OrthogonalAirfoil(Airfoil):
"""
Airfoils that can be decomposed to camber and thickness.
This class represents airfoils that are naturally described by a camber
curve and a thickness normal to the camber curve, both above and below
the camber curve.
The parameterization of the camber representation and the thickness
representation must be based on the same transformation.
Notes
-----
This family uses two related public airfoil parameterizations.
The native airfoil parameter ``u`` is a signed smooth surface parameter in
``[-1, 1]``.
Negative values lie on the lower surface and positive values lie on the
upper surface.
This parameterization is used for full-airfoil curve queries and is chosen
to remain smooth through the leading edge.
The surface-local airfoil parameter ``xi`` is a chord-like coordinate in
``[0, 1]`` measured from the leading edge to the trailing edge on one
selected surface.
For this family, ``xi = u**2`` and therefore ``u = +/- sqrt(xi)``, with
the sign determined by the selected surface.
This convention is intentional.
It removes the leading-edge square-root singularity from the thickness
representation while preserving a smooth native airfoil parameterization
for coordinate and derivative evaluation.
Breakpoint ownership is split across the component models.
Camber breakpoints come from :meth:`Camber.breakpoints`, thickness-side
discontinuities come from :meth:`Thickness.discontinuities`, and the
airfoil breakpoint helpers compose those one-sided component values into
exact airfoil-side derivatives when available.
The ordinary :meth:`xy_u` and :meth:`xy_uu` evaluators return the
``minus``-side value when a query lands exactly on an airfoil breakpoint.
"""
def __init__(self, *, camber: Camber, thickness: Thickness) -> None:
"""
Store the camber and thickness curves for the airfoil.
Parameters
----------
camber : Camber
Camber representation used by the airfoil model.
thickness : Thickness
Thickness distribution paired with ``camber``.
"""
super().__init__()
self._camber = camber
self._thickness = thickness
self._u_xmin: FloatScalar | None = None
self._u_xmax: FloatScalar | None = None
#
# Orthogonal components
#
@property
def camber(self) -> Camber:
"""
Return the camber function for airfoil.
This property exposes the camber-curve model associated with the
airfoil.
"""
return self._camber
@property
def thickness(self) -> Thickness:
"""
Return the thickness function for airfoil.
This property exposes the thickness model associated with the
airfoil.
"""
return self._thickness
[docs]
@override
def camber_curve(
self,
*,
num_points: int = 81,
spacing: Literal["uniform", "cosine"] = "cosine",
) -> AirfoilCamberResult:
"""
Return the exact native camber curve for this orthogonal airfoil.
Parameters
----------
num_points : int, default=81
Unused compatibility argument accepted for API consistency with
the approximate extraction path on other airfoil families.
spacing : {"uniform", "cosine"}, default="cosine"
Unused compatibility argument accepted for API consistency with
the approximate extraction path on other airfoil families.
Returns
-------
AirfoilCamberResult
Exact camber result exposing the native camber runtime.
"""
_ = (num_points, spacing)
return AirfoilCamberResult(
curve=self.camber,
mode="exact",
source_family=type(self).__name__,
)
#
# Cached geometric properties
#
@property
def xmin_parameter(self) -> FloatScalar:
"""
Return the parameter of the smallest x-coordinate for the airfoil.
This property caches the native airfoil parameter at the minimum
x-location.
Raises
------
RuntimeError
If the internal scalar root solve for the x-minimum does not
converge.
"""
if self._u_xmin is None:
max_camber_parameter = self._camber.max_camber_parameter()
_, max_camber_value = self._camber.xy_from_u(max_camber_parameter)
max_camber = as_float_scalar(max_camber_value)
if abs(max_camber) < _CAMBER_TOLERANCE:
self._u_xmin = 0.0
else:
if max_camber > 0:
u0 = 1e-7
u1 = 0.1
else:
u0 = -1e-7
u1 = -0.1
def x_derivative_residual(u: FloatScalar) -> FloatScalar:
return as_float_scalar(self.xy_u(u)[0])
def _solve_airfoil_root(
function: Callable[[FloatScalar], FloatScalar],
*,
bracket: tuple[FloatScalar, FloatScalar],
description: str,
) -> FloatScalar:
"""
Solve one scalar orthogonal-airfoil root robustly.
Raises
------
RuntimeError
If the scalar root solver fails to converge.
"""
result = root_scalar(
function,
bracket=bracket,
method="brentq",
xtol=_ROOT_SOLVER_TOLERANCE,
rtol=_ROOT_SOLVER_TOLERANCE,
maxiter=_ROOT_SOLVER_MAXITER,
)
if not result.converged: # pragma: no cover
msg = f"Failed to solve {description}: {result.flag}."
raise RuntimeError(msg)
return as_float_scalar(result.root)
self._u_xmin = _solve_airfoil_root(
x_derivative_residual,
bracket=(u0, u1),
description="orthogonal airfoil x-minimum parameter",
)
return self._u_xmin
@property
def xmax_parameter(self) -> FloatScalar:
"""
Return the parameter of the largest x-coordinate for the airfoil.
This property caches the native airfoil parameter at the maximum
x-location.
"""
if self._u_xmax is None:
if self.xy_from_u(-1.0)[0] >= self.xy_from_u(1.0)[0]:
self._u_xmax = -1.0
else:
self._u_xmax = 1.0
return self._u_xmax
#
# Inverse parameter queries
#
def _surface_x_bounds(
self, *, surface: AirfoilSurface
) -> tuple[FloatScalar, FloatScalar]:
"""
Return the reachable x-range for one orthogonal airfoil surface.
Notes
-----
Cambered orthogonal airfoils can reach their minimum x-value on one
surface slightly ahead of the nominal leading-edge point, so the
selected surface must include the relevant interior extremum.
"""
upper = surface_is_upper(surface)
bracket = UPPER_SURFACE_BRACKET if upper else LOWER_SURFACE_BRACKET
candidate_parameters: list[FloatScalar] = [bracket[0], bracket[1]]
if bracket[0] <= self.xmin_parameter <= bracket[1]:
candidate_parameters.append(self.xmin_parameter)
if bracket[0] <= self.xmax_parameter <= bracket[1]:
candidate_parameters.append(self.xmax_parameter)
x_values = [
as_float_scalar(self.xy_from_u(parameter)[0])
for parameter in candidate_parameters
]
return min(x_values), max(x_values)
#
# Curve geometry overrides
#
[docs]
@override
def xy_from_u(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Calculate the airfoil coordinates at one parameter location.
Parameters
----------
u : FloatInput
Signed surface parameter values in ``[-1, 1]``.
Returns
-------
tuple[FloatArray, FloatArray]
Tuple ``(x, y)`` of ``float64`` arrays matching the normalized
shape of ``u``.
"""
uc, sgn = self._split_u(u)
xi = self._xi_from_u(uc)[0]
_, eta = self.camber.xy_from_u(xi)
nx, ny = self.camber.normal(xi)
v = self._v_from_u(uc)[0]
delta = sgn * self.thickness.delta(v)
x = xi + delta * nx
y = eta + delta * ny
return x, y
[docs]
@override
def xy_u(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Calculate first derivatives with respect to the airfoil parameter.
Parameters
----------
u : FloatInput
Parameter for desired locations.
Returns
-------
tuple[FloatArray, FloatArray]
Tuple ``(dx/du, dy/du)`` of ``float64`` arrays matching the
normalized shape of ``u``.
Notes
-----
If ``u`` matches one of :meth:`breakpoints` exactly, this method
returns the ``minus``-side derivative limit.
"""
u_array = as_float_array(u)
x_u = as_float_array(np.empty_like(u_array))
y_u = as_float_array(np.empty_like(u_array))
analytic_mask = np.ones_like(u_array, dtype=bool)
breakpoint_values = self.breakpoints()
for breakpoint_index, breakpoint_value in enumerate(breakpoint_values):
breakpoint_mask = np.isclose(
u_array,
breakpoint_value,
atol=ROOT_ABS_TOLERANCE,
rtol=0.0,
)
if not np.any(breakpoint_mask):
continue
minus, _ = self.xy_u_breakpoint(index=breakpoint_index)
x_u[breakpoint_mask] = minus[0]
y_u[breakpoint_mask] = minus[1]
analytic_mask &= ~breakpoint_mask
if np.any(analytic_mask):
x_u[analytic_mask], y_u[analytic_mask] = self._xy_u_analytic(
u_array[analytic_mask]
)
return x_u, y_u
[docs]
@override
def xy_uu(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Return second derivatives of the airfoil coordinates.
Parameters
----------
u : FloatInput
Parameter for desired locations.
Returns
-------
tuple[FloatArray, FloatArray]
Tuple ``(d^2x/du^2, d^2y/du^2)`` of ``float64`` arrays matching the
normalized shape of ``u``.
Notes
-----
If ``u`` matches one of :meth:`breakpoints` exactly, this method
returns the ``minus``-side derivative limit.
"""
u_array = as_float_array(u)
x_uu = as_float_array(np.empty_like(u_array))
y_uu = as_float_array(np.empty_like(u_array))
analytic_mask = np.ones_like(u_array, dtype=bool)
breakpoint_values = self.breakpoints()
for breakpoint_index, breakpoint_value in enumerate(breakpoint_values):
breakpoint_mask = np.isclose(
u_array,
breakpoint_value,
atol=ROOT_ABS_TOLERANCE,
rtol=0.0,
)
if not np.any(breakpoint_mask):
continue
minus, _ = self.xy_uu_breakpoint(index=breakpoint_index)
x_uu[breakpoint_mask] = minus[0]
y_uu[breakpoint_mask] = minus[1]
analytic_mask &= ~breakpoint_mask
if np.any(analytic_mask):
x_uu[analytic_mask], y_uu[analytic_mask] = self._xy_uu_analytic(
u_array[analytic_mask]
)
return x_uu, y_uu
#
# Surface parameter conversion API
#
[docs]
@override
def u_from_xi(
self, xi: FloatInput, *, surface: AirfoilSurface
) -> FloatArray:
"""
Convert one-surface ``xi`` coordinates to curve parameters.
Parameters
----------
xi : FloatInput
Surface-local coordinates in ``[0, 1]`` measured from the
leading edge to the trailing edge.
surface : {"lower", "upper"}
Surface to evaluate.
Returns
-------
FloatArray
Curve parameters matching ``xi`` on the selected surface.
"""
u = self._u_from_xi(xi)
return u if surface_is_upper(surface) else -u
[docs]
@override
def xi_from_u(self, u: FloatInput) -> SurfaceMappedValues:
"""
Convert curve airfoil parameters to surface-local ``xi`` values.
Parameters
----------
u : FloatInput
Curve airfoil parameters in ``[-1, 1]``.
Returns
-------
SurfaceMappedValues
Surface-local ``xi`` values and upper-surface membership flags.
Notes
-----
Concrete airfoil families define this mapping because ``xi`` need not
equal ``|u|`` for every airfoil parameterization.
"""
u_array = as_float_array(u)
xi = self._xi_from_u(u)[0]
return SurfaceMappedValues(value=xi, is_upper=u_array >= 0.0)
#
# Orthogonal model queries
#
[docs]
def camber_location(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Return the camber-line location associated with ``u``.
Parameters
----------
u : FloatInput
Signed or unsigned airfoil parameter values.
Returns
-------
tuple[FloatArray, FloatArray]
Tuple ``(x, y)`` of camber-line coordinates evaluated at the
mapped camber parameter.
"""
return self.camber.xy_from_u(self._xi_from_u(u)[0])
[docs]
def thickness_value(self, u: FloatInput) -> FloatArray:
"""
Return the thickness offset associated with ``u``.
Parameters
----------
u : FloatInput
Signed or unsigned airfoil parameter values.
Returns
-------
FloatArray
Thickness magnitude evaluated at the mapped thickness
parameter.
"""
return self.thickness.delta(self._v_from_u(u)[0])
#
# Breakpoint interface
#
[docs]
@override
def breakpoints(self) -> list[FloatScalar]:
"""
Return the locations of breakpoints or discontinuities.
Returns
-------
list[float]
Sorted unique parameter values induced by camber-line joins,
thickness discontinuities, and their mirrored surface partners.
Notes
-----
This list always includes ``-1``, ``0``, and ``1``.
Additional interior values are induced by mapped camber breakpoints or
mapped thickness discontinuities that require one-sided derivative
information at the airfoil level.
"""
camber_joints = [
as_float_scalar(self._u_from_xi(j))
for j in self.camber.breakpoints()
]
thickness_joints = [
as_float_scalar(self._u_from_v(d))
for d in self.thickness.discontinuities()
]
joints = camber_joints + thickness_joints
mirrored_joints = joints + [-value for value in joints]
unique_joints = np.unique(as_float_array(mirrored_joints))
return [as_float_scalar(value) for value in unique_joints]
[docs]
@override
def xy_u_breakpoint(
self,
*,
index: int,
) -> tuple[
tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]
]:
"""
Return one-sided first derivatives at one airfoil breakpoint.
Notes
-----
This method composes the one-sided camber and thickness values that
apply on each side of the requested airfoil breakpoint.
It is the exact-breakpoint contract paired with :meth:`xy_u`.
"""
u_breakpoint = self.breakpoints()[index]
u_parameter_minus, u_parameter_plus = self.breakpoint_parameter_limits(
index=index
)
return (
self._xy_u_breakpoint_parameter_side(
u_breakpoint=u_breakpoint,
u_parameter_side=u_parameter_minus,
),
self._xy_u_breakpoint_parameter_side(
u_breakpoint=u_breakpoint,
u_parameter_side=u_parameter_plus,
),
)
[docs]
@override
def xy_uu_breakpoint(
self,
*,
index: int,
) -> tuple[
tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]
]:
"""
Return one-sided second derivatives at one airfoil breakpoint.
Notes
-----
This method composes the one-sided camber and thickness values that
apply on each side of the requested airfoil breakpoint.
It is the exact-breakpoint contract paired with :meth:`xy_uu`.
"""
u_breakpoint = self.breakpoints()[index]
u_parameter_minus, u_parameter_plus = self.breakpoint_parameter_limits(
index=index
)
return (
self._xy_uu_breakpoint_parameter_side(
u_breakpoint=u_breakpoint,
u_parameter_side=u_parameter_minus,
),
self._xy_uu_breakpoint_parameter_side(
u_breakpoint=u_breakpoint,
u_parameter_side=u_parameter_plus,
),
)
[docs]
@override
def xy_s_breakpoint(
self,
*,
index: int,
) -> tuple[
tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]
]:
"""
Return one-sided arc-length derivatives at one airfoil breakpoint.
Notes
-----
This method composes the exact arc-length tangent values from the
exact native breakpoint derivatives returned by
:meth:`xy_u_breakpoint`.
When the native derivative vanishes at a cusp-like breakpoint, such as
the orthogonal leading edge, this method uses the one-sided native
second derivative to recover the exact one-sided tangent direction.
"""
minus_u, plus_u = self.xy_u_breakpoint(index=index)
minus_uu, plus_uu = self.xy_uu_breakpoint(index=index)
if np.isclose(np.hypot(*minus_u), 0.0, atol=0.0, rtol=0.0):
minus_s = self._xy_s_from_zero_speed_breakpoint(
native_uu=minus_uu,
side_sign=as_float_scalar(-1.0),
)
else:
minus_s = self._xy_s_from_native_breakpoint(minus_u)
if np.isclose(np.hypot(*plus_u), 0.0, atol=0.0, rtol=0.0):
plus_s = self._xy_s_from_zero_speed_breakpoint(
native_uu=plus_uu,
side_sign=as_float_scalar(1.0),
)
else:
plus_s = self._xy_s_from_native_breakpoint(plus_u)
return minus_s, plus_s
[docs]
@override
def xy_ss_breakpoint(
self,
*,
index: int,
) -> tuple[
tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]
]:
"""
Return one-sided arc-length second derivatives at one breakpoint.
Notes
-----
This method composes the exact arc-length curvature-vector values from
the exact native breakpoint derivatives returned by
:meth:`xy_u_breakpoint` and :meth:`xy_uu_breakpoint`.
When the native derivative vanishes at a cusp-like breakpoint, such as
the orthogonal leading edge, the generic sampled fallback is retained
because the exact arc-length second derivative depends on higher-order
native data, such as one-sided ``xy_uuu`` values, that this runtime
does not expose.
"""
minus_u, plus_u = self.xy_u_breakpoint(index=index)
if np.isclose(np.hypot(*minus_u), 0.0) or np.isclose(
np.hypot(*plus_u),
0.0,
):
return super().xy_ss_breakpoint(index=index)
minus_uu, plus_uu = self.xy_uu_breakpoint(index=index)
return (
self._xy_ss_from_native_breakpoint(minus_u, minus_uu),
self._xy_ss_from_native_breakpoint(plus_u, plus_uu),
)
@staticmethod
def _xy_s_from_native_breakpoint(
native_u: tuple[FloatScalar, FloatScalar],
) -> tuple[FloatScalar, FloatScalar]:
"""
Return one exact arc-length tangent from native derivatives.
Parameters
----------
native_u : tuple[FloatScalar, FloatScalar]
Exact one-sided native derivative tuple ``(x_u, y_u)``.
Returns
-------
tuple[FloatScalar, FloatScalar]
Exact arc-length tangent ``(x_s, y_s)`` on the same breakpoint
side.
"""
x_u, y_u = native_u
speed = as_float_scalar(np.hypot(x_u, y_u))
return (
as_float_scalar(x_u / speed),
as_float_scalar(y_u / speed),
)
@staticmethod
def _xy_s_from_zero_speed_breakpoint(
*,
native_uu: tuple[FloatScalar, FloatScalar],
side_sign: FloatScalar,
) -> tuple[FloatScalar, FloatScalar]:
"""
Return one exact arc-length tangent at a zero-speed breakpoint.
Parameters
----------
native_uu : tuple[FloatScalar, FloatScalar]
Exact one-sided native second-derivative tuple
``(x_uu, y_uu)``.
side_sign : FloatScalar
Sign of the one-sided parameter offset used to approach the
breakpoint, typically ``-1`` on the minus side and ``+1`` on the
plus side.
Returns
-------
tuple[FloatScalar, FloatScalar]
Exact one-sided arc-length tangent recovered from the directional
limit of ``xy_u / ||xy_u||`` as the native speed tends to zero.
"""
x_uu, y_uu = native_uu
speed = as_float_scalar(np.hypot(x_uu, y_uu))
return (
as_float_scalar(side_sign * x_uu / speed),
as_float_scalar(side_sign * y_uu / speed),
)
@staticmethod
def _xy_ss_from_native_breakpoint(
native_u: tuple[FloatScalar, FloatScalar],
native_uu: tuple[FloatScalar, FloatScalar],
) -> tuple[FloatScalar, FloatScalar]:
"""
Return one exact arc-length second derivative from native data.
Parameters
----------
native_u : tuple[FloatScalar, FloatScalar]
Exact one-sided native derivative tuple ``(x_u, y_u)``.
native_uu : tuple[FloatScalar, FloatScalar]
Exact one-sided native second-derivative tuple
``(x_uu, y_uu)``.
Returns
-------
tuple[FloatScalar, FloatScalar]
Exact arc-length second derivative ``(x_ss, y_ss)`` on the same
breakpoint side.
"""
x_u, y_u = native_u
x_uu, y_uu = native_uu
speed_sq = as_float_scalar(x_u**2 + y_u**2)
speed_pow4 = as_float_scalar(speed_sq**2)
projection = as_float_scalar(x_u * x_uu + y_u * y_uu)
x_ss = as_float_scalar(x_uu / speed_sq - x_u * projection / speed_pow4)
y_ss = as_float_scalar(y_uu / speed_sq - y_u * projection / speed_pow4)
if np.isclose(x_ss, 0.0, atol=_BREAKPOINT_ZERO_TOLERANCE, rtol=0.0):
x_ss = as_float_scalar(0.0)
if np.isclose(y_ss, 0.0, atol=_BREAKPOINT_ZERO_TOLERANCE, rtol=0.0):
y_ss = as_float_scalar(0.0)
return x_ss, y_ss
@staticmethod
def _breakpoint_index(
value: FloatScalar,
*,
breakpoints: list[FloatScalar],
) -> int | None:
"""
Return the matching breakpoint index for one orthogonal parameter.
Parameters
----------
value : FloatScalar
Candidate breakpoint parameter.
breakpoints : list[FloatScalar]
Ascending list of valid breakpoint parameters.
Returns
-------
int | None
Matching breakpoint index when ``value`` lies within the root
tolerance of one breakpoint, otherwise ``None``.
"""
matches = np.flatnonzero(
np.isclose(
as_float_array(breakpoints),
value,
atol=ROOT_ABS_TOLERANCE,
rtol=0.0,
)
)
if matches.size == 0:
return None
return int(matches[0])
#
# Cache invalidation
#
def _airfoil_changed(self) -> None:
"""
Clear cached values after a geometry change.
Returns
-------
None
This method updates internal caches in place.
"""
self._u_xmin = None
self._u_xmax = None
super()._airfoil_changed()
def _xy_u_analytic(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Return native first derivatives away from exact breakpoint dispatch.
Parameters
----------
u : FloatInput
Signed orthogonal airfoil parameter.
Returns
-------
tuple[FloatArray, FloatArray]
``(x_u, y_u)`` evaluated directly from the smooth camber and
thickness relations without consulting the exact breakpoint-side
helper APIs.
"""
uc, sgn = self._split_u(u)
xi, xi_u = self._xi_from_u(uc)[0:2]
xc_xi, yc_xi = self.camber.xy_u(xi)
k = self.camber.k(xi)
nx, ny = self.camber.normal(xi)
nx_xi = -xc_xi * k
ny_xi = -yc_xi * k
v, v_u = self._v_from_u(uc)[0:2]
delta = sgn * self.thickness.delta(v)
delta_v = sgn * self.thickness.delta_v(v)
x_u = sgn * (xi_u * xc_xi + v_u * delta_v * nx + delta * xi_u * nx_xi)
y_u = sgn * (xi_u * yc_xi + v_u * delta_v * ny + delta * xi_u * ny_xi)
return x_u, y_u
def _xy_u_breakpoint_parameter_side(
self,
*,
u_breakpoint: FloatScalar,
u_parameter_side: FloatScalar,
) -> tuple[FloatScalar, FloatScalar]:
"""
Return one exact native first derivative on one side of a breakpoint.
Parameters
----------
u_breakpoint : FloatScalar
Signed airfoil breakpoint parameter being evaluated.
u_parameter_side : FloatScalar
Nearby signed parameter used only to select the ``minus`` or
``plus`` side for any camber or thickness component breakpoints
that map to ``u_breakpoint``.
Returns
-------
tuple[FloatScalar, FloatScalar]
Exact one-sided ``(x_u, y_u)`` value composed from the matching
component breakpoint-side data.
"""
u_magnitude = (
as_float_scalar(-u_breakpoint)
if u_breakpoint < 0.0
else as_float_scalar(u_breakpoint)
)
side_magnitude = (
as_float_scalar(-u_parameter_side)
if u_parameter_side < 0.0
else as_float_scalar(u_parameter_side)
)
surface_sign = self._surface_sign_from_u(u_parameter_side)
xi_breakpoint, xi_u, _ = self._xi_from_u(u_magnitude)
xi_parameter_side = as_float_scalar(self._xi_from_u(side_magnitude)[0])
v_breakpoint, v_u, _ = self._v_from_u(u_magnitude)
v_parameter_side = as_float_scalar(self._v_from_u(side_magnitude)[0])
xc_xi, yc_xi = self._camber_xy_u_parameter_side(
xi_breakpoint=as_float_scalar(xi_breakpoint),
xi_parameter_side=xi_parameter_side,
)
k = self._camber_k_parameter_side(
xi_breakpoint=as_float_scalar(xi_breakpoint),
xi_parameter_side=xi_parameter_side,
)
nx, ny = self.camber.normal(xi_breakpoint)
nx_xi = -xc_xi * k
ny_xi = -yc_xi * k
delta = surface_sign * self._thickness_delta_parameter_side(
v_breakpoint=as_float_scalar(v_breakpoint),
v_parameter_side=v_parameter_side,
)
delta_v = surface_sign * self._thickness_delta_v_parameter_side(
v_breakpoint=as_float_scalar(v_breakpoint),
v_parameter_side=v_parameter_side,
)
x_u = surface_sign * (
as_float_scalar(xi_u) * xc_xi
+ as_float_scalar(v_u) * delta_v * as_float_scalar(nx)
+ delta * as_float_scalar(xi_u) * nx_xi
)
y_u = surface_sign * (
as_float_scalar(xi_u) * yc_xi
+ as_float_scalar(v_u) * delta_v * as_float_scalar(ny)
+ delta * as_float_scalar(xi_u) * ny_xi
)
return as_float_scalar(x_u), as_float_scalar(y_u)
def _xy_uu_analytic(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Return native second derivatives away from exact breakpoint dispatch.
Parameters
----------
u : FloatInput
Signed orthogonal airfoil parameter.
Returns
-------
tuple[FloatArray, FloatArray]
``(x_uu, y_uu)`` evaluated directly from the smooth camber and
thickness relations without consulting the exact breakpoint-side
helper APIs.
"""
uc, sgn = self._split_u(u)
xi, xi_u, xi_uu = self._xi_from_u(uc)
xc_xi, yc_xi = self.camber.xy_u(xi)
xc_xixi, yc_xixi = self.camber.xy_uu(xi)
k = self.camber.k(xi)
k_xi = self.camber.k_u(xi)
nx, ny = self.camber.normal(xi)
nx_xi = -xc_xi * k
nx_xixi = -(xc_xixi * k + xc_xi * k_xi)
ny_xi = -yc_xi * k
ny_xixi = -(yc_xixi * k + yc_xi * k_xi)
v, v_u, v_uu = self._v_from_u(uc)
delta = sgn * self.thickness.delta(v)
delta_v = sgn * self.thickness.delta_v(v)
delta_vv = sgn * self.thickness.delta_vv(v)
x_uu = (
xi_uu * xc_xi
+ xi_u**2 * xc_xixi
+ v_u**2 * nx * delta_vv
+ xi_u**2 * delta * nx_xixi
+ 2 * xi_u * v_u * delta_v * nx_xi
+ v_uu * nx * delta_v
+ xi_uu * delta * nx_xi
)
y_uu = (
xi_uu * yc_xi
+ xi_u**2 * yc_xixi
+ v_u**2 * ny * delta_vv
+ xi_u**2 * delta * ny_xixi
+ 2 * xi_u * v_u * delta_v * ny_xi
+ v_uu * ny * delta_v
+ xi_uu * delta * ny_xi
)
return x_uu, y_uu
def _xy_uu_breakpoint_parameter_side(
self,
*,
u_breakpoint: FloatScalar,
u_parameter_side: FloatScalar,
) -> tuple[FloatScalar, FloatScalar]:
"""
Return one exact native second derivative on one breakpoint side.
Parameters
----------
u_breakpoint : FloatScalar
Signed airfoil breakpoint parameter being evaluated.
u_parameter_side : FloatScalar
Nearby signed parameter used only to select the ``minus`` or
``plus`` side for any mapped component breakpoints.
Returns
-------
tuple[FloatScalar, FloatScalar]
Exact one-sided ``(x_uu, y_uu)`` value composed from the matching
component breakpoint-side data.
"""
u_magnitude = (
as_float_scalar(-u_breakpoint)
if u_breakpoint < 0.0
else as_float_scalar(u_breakpoint)
)
side_magnitude = (
as_float_scalar(-u_parameter_side)
if u_parameter_side < 0.0
else as_float_scalar(u_parameter_side)
)
surface_sign = self._surface_sign_from_u(u_parameter_side)
xi_breakpoint, xi_u, xi_uu = self._xi_from_u(u_magnitude)
xi_parameter_side = as_float_scalar(self._xi_from_u(side_magnitude)[0])
v_breakpoint, v_u, v_uu = self._v_from_u(u_magnitude)
v_parameter_side = as_float_scalar(self._v_from_u(side_magnitude)[0])
xc_xi, yc_xi = self._camber_xy_u_parameter_side(
xi_breakpoint=as_float_scalar(xi_breakpoint),
xi_parameter_side=xi_parameter_side,
)
xc_xixi, yc_xixi = self._camber_xy_uu_parameter_side(
xi_breakpoint=as_float_scalar(xi_breakpoint),
xi_parameter_side=xi_parameter_side,
)
k = self._camber_k_parameter_side(
xi_breakpoint=as_float_scalar(xi_breakpoint),
xi_parameter_side=xi_parameter_side,
)
k_xi = self._camber_k_u_parameter_side(
xi_breakpoint=as_float_scalar(xi_breakpoint),
xi_parameter_side=xi_parameter_side,
)
nx, ny = self.camber.normal(xi_breakpoint)
nx_xi = -xc_xi * k
nx_xixi = -(xc_xixi * k + xc_xi * k_xi)
ny_xi = -yc_xi * k
ny_xixi = -(yc_xixi * k + yc_xi * k_xi)
delta = surface_sign * self._thickness_delta_parameter_side(
v_breakpoint=as_float_scalar(v_breakpoint),
v_parameter_side=v_parameter_side,
)
delta_v = surface_sign * self._thickness_delta_v_parameter_side(
v_breakpoint=as_float_scalar(v_breakpoint),
v_parameter_side=v_parameter_side,
)
delta_vv = surface_sign * self._thickness_delta_vv_parameter_side(
v_breakpoint=as_float_scalar(v_breakpoint),
v_parameter_side=v_parameter_side,
)
x_uu = (
as_float_scalar(xi_uu) * xc_xi
+ as_float_scalar(xi_u) ** 2 * xc_xixi
+ as_float_scalar(v_u) ** 2 * as_float_scalar(nx) * delta_vv
+ as_float_scalar(xi_u) ** 2 * delta * nx_xixi
+ 2 * as_float_scalar(xi_u) * as_float_scalar(v_u) * delta_v * nx_xi
+ as_float_scalar(v_uu) * as_float_scalar(nx) * delta_v
+ as_float_scalar(xi_uu) * delta * nx_xi
)
y_uu = (
as_float_scalar(xi_uu) * yc_xi
+ as_float_scalar(xi_u) ** 2 * yc_xixi
+ as_float_scalar(v_u) ** 2 * as_float_scalar(ny) * delta_vv
+ as_float_scalar(xi_u) ** 2 * delta * ny_xixi
+ 2 * as_float_scalar(xi_u) * as_float_scalar(v_u) * delta_v * ny_xi
+ as_float_scalar(v_uu) * as_float_scalar(ny) * delta_v
+ as_float_scalar(xi_uu) * delta * ny_xi
)
return as_float_scalar(x_uu), as_float_scalar(y_uu)
def _camber_xy_u_parameter_side(
self,
*,
xi_breakpoint: FloatScalar,
xi_parameter_side: FloatScalar,
) -> tuple[FloatScalar, FloatScalar]:
"""
Return camber first derivatives on the selected side of ``xi``.
Parameters
----------
xi_breakpoint : FloatScalar
Chordwise parameter at which a camber breakpoint may occur.
xi_parameter_side : FloatScalar
Nearby chordwise parameter used only to choose the ``minus`` or
``plus`` side when ``xi_breakpoint`` is a real camber breakpoint.
Returns
-------
tuple[FloatScalar, FloatScalar]
Exact one-sided camber ``(x_xi, y_xi)`` values.
"""
breakpoint_index = self._breakpoint_index(
xi_breakpoint,
breakpoints=self.camber.breakpoints(),
)
if breakpoint_index is None:
x_u, y_u = self.camber.xy_u(xi_breakpoint)
return as_float_scalar(x_u), as_float_scalar(y_u)
minus_values, plus_values = self.camber.xy_u_breakpoint(
index=breakpoint_index
)
return self._select_parameter_side_tuple(
minus_values,
plus_values,
breakpoint_parameter=xi_breakpoint,
parameter_side=xi_parameter_side,
)
def _camber_xy_uu_parameter_side(
self,
*,
xi_breakpoint: FloatScalar,
xi_parameter_side: FloatScalar,
) -> tuple[FloatScalar, FloatScalar]:
"""
Return camber second derivatives on the selected side of ``xi``.
Parameters
----------
xi_breakpoint : FloatScalar
Chordwise parameter at which a camber breakpoint may occur.
xi_parameter_side : FloatScalar
Nearby chordwise parameter used only to choose the ``minus`` or
``plus`` side when ``xi_breakpoint`` is a real camber breakpoint.
Returns
-------
tuple[FloatScalar, FloatScalar]
Exact one-sided camber ``(x_xixi, y_xixi)`` values.
"""
breakpoint_index = self._breakpoint_index(
xi_breakpoint,
breakpoints=self.camber.breakpoints(),
)
if breakpoint_index is None:
x_uu, y_uu = self.camber.xy_uu(xi_breakpoint)
return as_float_scalar(x_uu), as_float_scalar(y_uu)
minus_values, plus_values = self.camber.xy_uu_breakpoint(
index=breakpoint_index
)
return self._select_parameter_side_tuple(
minus_values,
plus_values,
breakpoint_parameter=xi_breakpoint,
parameter_side=xi_parameter_side,
)
def _camber_k_parameter_side(
self,
*,
xi_breakpoint: FloatScalar,
xi_parameter_side: FloatScalar,
) -> FloatScalar:
"""
Return camber curvature on the selected side of ``xi``.
Parameters
----------
xi_breakpoint : FloatScalar
Chordwise parameter at which a camber breakpoint may occur.
xi_parameter_side : FloatScalar
Nearby chordwise parameter used only to choose the ``minus`` or
``plus`` side when ``xi_breakpoint`` is a real camber breakpoint.
Returns
-------
FloatScalar
Exact one-sided camber curvature.
"""
breakpoint_index = self._breakpoint_index(
xi_breakpoint,
breakpoints=self.camber.breakpoints(),
)
if breakpoint_index is None:
return as_float_scalar(self.camber.k(xi_breakpoint))
minus_value, plus_value = self.camber.k_breakpoint(
index=breakpoint_index
)
return self._select_parameter_side_scalar(
minus_value,
plus_value,
breakpoint_parameter=xi_breakpoint,
parameter_side=xi_parameter_side,
)
def _camber_k_u_parameter_side(
self,
*,
xi_breakpoint: FloatScalar,
xi_parameter_side: FloatScalar,
) -> FloatScalar:
"""
Return camber curvature derivative on the selected side of ``xi``.
Parameters
----------
xi_breakpoint : FloatScalar
Chordwise parameter at which a camber breakpoint may occur.
xi_parameter_side : FloatScalar
Nearby chordwise parameter used only to choose the ``minus`` or
``plus`` side when ``xi_breakpoint`` is a real camber breakpoint.
Returns
-------
FloatScalar
Exact one-sided derivative of camber curvature with respect to
``xi``.
"""
breakpoint_index = self._breakpoint_index(
xi_breakpoint,
breakpoints=self.camber.breakpoints(),
)
if breakpoint_index is None:
return as_float_scalar(self.camber.k_u(xi_breakpoint))
minus_value, plus_value = self.camber.k_u_breakpoint(
index=breakpoint_index
)
return self._select_parameter_side_scalar(
minus_value,
plus_value,
breakpoint_parameter=xi_breakpoint,
parameter_side=xi_parameter_side,
)
def _thickness_delta_parameter_side(
self,
*,
v_breakpoint: FloatScalar,
v_parameter_side: FloatScalar,
) -> FloatScalar:
"""
Return thickness on the selected side of one ``v`` location.
Parameters
----------
v_breakpoint : FloatScalar
Thickness parameter at which a true thickness discontinuity may
occur.
v_parameter_side : FloatScalar
Nearby thickness parameter used only to choose the ``minus`` or
``plus`` side when ``v_breakpoint`` is reported by
:meth:`Thickness.discontinuities`.
Returns
-------
FloatScalar
Exact one-sided thickness value.
"""
breakpoint_index = self._breakpoint_index(
v_breakpoint,
breakpoints=self.thickness.discontinuities(),
)
if breakpoint_index is None:
return as_float_scalar(self.thickness.delta(v_breakpoint))
minus_value, plus_value = self.thickness.delta_breakpoint(
index=breakpoint_index
)
return self._select_parameter_side_scalar(
minus_value,
plus_value,
breakpoint_parameter=v_breakpoint,
parameter_side=v_parameter_side,
)
def _thickness_delta_v_parameter_side(
self,
*,
v_breakpoint: FloatScalar,
v_parameter_side: FloatScalar,
) -> FloatScalar:
"""
Return thickness first derivative on the selected side of ``v``.
Parameters
----------
v_breakpoint : FloatScalar
Thickness parameter at which a true thickness discontinuity may
occur.
v_parameter_side : FloatScalar
Nearby thickness parameter used only to choose the ``minus`` or
``plus`` side when ``v_breakpoint`` is reported by
:meth:`Thickness.discontinuities`.
Returns
-------
FloatScalar
Exact one-sided first derivative of thickness with respect to
``v``.
"""
breakpoint_index = self._breakpoint_index(
v_breakpoint,
breakpoints=self.thickness.discontinuities(),
)
if breakpoint_index is None:
return as_float_scalar(self.thickness.delta_v(v_breakpoint))
minus_value, plus_value = self.thickness.delta_v_breakpoint(
index=breakpoint_index
)
return self._select_parameter_side_scalar(
minus_value,
plus_value,
breakpoint_parameter=v_breakpoint,
parameter_side=v_parameter_side,
)
def _thickness_delta_vv_parameter_side(
self,
*,
v_breakpoint: FloatScalar,
v_parameter_side: FloatScalar,
) -> FloatScalar:
"""
Return thickness second derivative on the selected side of ``v``.
Parameters
----------
v_breakpoint : FloatScalar
Thickness parameter at which a true thickness discontinuity may
occur.
v_parameter_side : FloatScalar
Nearby thickness parameter used only to choose the ``minus`` or
``plus`` side when ``v_breakpoint`` is reported by
:meth:`Thickness.discontinuities`.
Returns
-------
FloatScalar
Exact one-sided second derivative of thickness with respect to
``v``.
"""
breakpoint_index = self._breakpoint_index(
v_breakpoint,
breakpoints=self.thickness.discontinuities(),
)
if breakpoint_index is None:
return as_float_scalar(self.thickness.delta_vv(v_breakpoint))
minus_value, plus_value = self.thickness.delta_vv_breakpoint(
index=breakpoint_index
)
return self._select_parameter_side_scalar(
minus_value,
plus_value,
breakpoint_parameter=v_breakpoint,
parameter_side=v_parameter_side,
)
@staticmethod
def _select_parameter_side_scalar(
minus_value: FloatScalar,
plus_value: FloatScalar,
*,
breakpoint_parameter: FloatScalar,
parameter_side: FloatScalar,
) -> FloatScalar:
"""
Select one scalar breakpoint-side value by parameter ordering.
Parameters
----------
minus_value : FloatScalar
Value approached from parameters less than or equal to the
breakpoint.
plus_value : FloatScalar
Value approached from parameters greater than the breakpoint.
breakpoint_parameter : FloatScalar
Breakpoint parameter separating the two sides.
parameter_side : FloatScalar
Nearby parameter used to choose the side. Exact matches default to
the ``minus`` side.
"""
if parameter_side > breakpoint_parameter and not np.isclose(
parameter_side,
breakpoint_parameter,
atol=ROOT_ABS_TOLERANCE,
rtol=0.0,
):
return plus_value
return minus_value
@staticmethod
def _select_parameter_side_tuple(
minus_values: tuple[FloatScalar, FloatScalar],
plus_values: tuple[FloatScalar, FloatScalar],
*,
breakpoint_parameter: FloatScalar,
parameter_side: FloatScalar,
) -> tuple[FloatScalar, FloatScalar]:
"""
Select one vector-valued breakpoint-side value by ordering.
Parameters
----------
minus_values : tuple[FloatScalar, FloatScalar]
Tuple approached from parameters less than or equal to the
breakpoint.
plus_values : tuple[FloatScalar, FloatScalar]
Tuple approached from parameters greater than the breakpoint.
breakpoint_parameter : FloatScalar
Breakpoint parameter separating the two sides.
parameter_side : FloatScalar
Nearby parameter used to choose the side. Exact matches default to
the ``minus`` side.
"""
if parameter_side > breakpoint_parameter and not np.isclose(
parameter_side,
breakpoint_parameter,
atol=ROOT_ABS_TOLERANCE,
rtol=0.0,
):
return plus_values
return minus_values
@staticmethod
def _surface_sign_from_u(u: FloatScalar) -> FloatScalar:
"""
Return the signed surface orientation associated with ``u``.
Parameters
----------
u : FloatScalar
Signed orthogonal airfoil parameter.
Returns
-------
FloatScalar
``-1`` on the lower surface and ``+1`` on the upper surface.
"""
return as_float_scalar(-1.0 if u < 0.0 else 1.0)
#
# Internal parameter conversion helpers
#
@staticmethod
def _xi_from_u(
u: FloatInput,
) -> tuple[FloatArray, FloatArray, FloatArray]:
"""
Convert unsigned curve magnitude to chord parameterization.
Parameters
----------
u : FloatInput
Unsigned curve magnitude in ``[0, 1]``.
Returns
-------
tuple[FloatArray, FloatArray, FloatArray]
Chord parameter ``xi`` plus its first and second derivatives with
respect to the unsigned curve magnitude ``u``.
"""
u_array = as_float_array(u)
return u_array**2, 2 * u_array, 2 * np.ones_like(u_array)
@staticmethod
def _u_from_xi(xi: FloatInput) -> FloatArray:
"""
Convert chord parameterization to unsigned curve magnitude.
Parameters
----------
xi : FloatInput
Chord parameter in ``[0, 1]``.
Returns
-------
FloatArray
Unsigned curve magnitude in ``[0, 1]``.
"""
xi_array = as_float_array(xi)
return as_float_array(np.sqrt(xi_array))
@staticmethod
def _v_from_u(
u: FloatInput,
) -> tuple[FloatArray, FloatArray, FloatArray]:
"""
Convert unsigned curve magnitude to thickness parameterization.
Parameters
----------
u : FloatInput
Unsigned curve magnitude in ``[0, 1]``.
Returns
-------
tuple[FloatArray, FloatArray, FloatArray]
Thickness parameter ``v`` plus its first and second derivatives
with respect to the unsigned curve magnitude ``u``.
"""
u_array = as_float_array(u)
return u_array, np.ones_like(u_array), np.zeros_like(u_array)
@staticmethod
def _u_from_v(v: FloatInput) -> FloatArray:
"""
Convert thickness parameterization to unsigned curve magnitude.
Parameters
----------
v : FloatInput
Thickness parameter in ``[0, 1]``.
Returns
-------
FloatArray
Unsigned curve magnitude in ``[0, 1]``.
"""
return as_float_array(v)
@staticmethod
def _split_u(u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Split the curve parameter into magnitude and sign.
Parameters
----------
u : FloatInput
Curve parameter.
Returns
-------
tuple[FloatArray, FloatArray]
Corresponding magnitude of curve parameter and sign.
"""
u_array = as_float_array(u)
sign = np.where(u_array < 0.0, -1.0, 1.0)
return as_float_array(np.abs(u_array)), as_float_array(sign)