ENH: rework aerodynamic coefficients onto a body-frame GenericSurface

Every aerodynamic surface is now rooted in GenericSurface and stores its
force coefficients in the body frame (cN/cY/cA) plus the cm/cn/cl moments,
exposing all nine coefficients (cL/cD/cQ/cN/cY/cA/cm/cn/cl) with the wind
trio lazily derived. A force_convention argument lets users supply wind- or
body-frame coefficients. Barrowman surfaces (nose, tail, fin sets) keep the
classic geometric normal-force/moment method, report the force at the
geometric center of pressure via the classic 180-degree surface rotation,
and expose cN_alpha/cY_beta stability slopes (the old clalpha relabelled).

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

Author: MateusStano

rocketpy/mathutils/function.py

@@ -40,6 +40,32 @@
     "regular_grid": 6,
 }
 EXTRAPOLATION_TYPES = {"zero": 0, "natural": 1, "constant": 2}
+# Maps a requested interpolation name onto a scipy ``RegularGridInterpolator``
+# ``method`` for gridded (N-D Cartesian) data. The 1-D-only names ``spline`` and
+# ``akima`` fall back to their closest grid analogs (``cubic`` and the
+# shape-preserving ``pchip``); anything unrecognized defaults to ``linear``.
+REGULAR_GRID_METHODS = {
+    "linear": "linear",
+    "nearest": "nearest",
+    "slinear": "slinear",
+    "cubic": "cubic",
+    "quintic": "quintic",
+    "pchip": "pchip",
+    "spline": "cubic",
+    "akima": "pchip",
+    "polynomial": "cubic",
+}
+# Minimum points per axis required by each ``RegularGridInterpolator`` method.
+# A grid with fewer samples on any axis cannot use the higher-order methods, so
+# the caller falls back to linear rather than letting SciPy raise mid-build.
+REGULAR_GRID_MIN_POINTS = {
+    "nearest": 1,
+    "linear": 2,
+    "slinear": 2,
+    "pchip": 2,
+    "cubic": 4,
+    "quintic": 6,
+}
 
 
 class SourceType(Enum):
@@ -157,7 +183,12 @@ def __init__(
 
     @classmethod
     def from_regular_grid_csv(
-        cls, csv_source, variable_names, coeff_name, extrapolation
+        cls,
+        csv_source,
+        variable_names,
+        coeff_name,
+        extrapolation,
+        interpolation="linear",
     ):
         """Create a regular-grid Function from CSV samples when possible.
 
@@ -171,6 +202,14 @@ def from_regular_grid_csv(
             Name of the output coefficient.
         extrapolation : str
             Extrapolation method passed to the Function constructor.
+        interpolation : str, optional
+            Requested interpolation. Mapped onto a
+            :class:`scipy.interpolate.RegularGridInterpolator` ``method`` via
+            :data:`REGULAR_GRID_METHODS` (e.g. ``"spline"`` -> ``"cubic"``,
+            ``"akima"`` -> ``"pchip"``); unrecognized names fall back to
+            ``"linear"``. Smooth methods require enough points per axis
+            (``"cubic"`` needs at least 4), otherwise SciPy raises. Default
+            ``"linear"``.
 
         Returns
         -------
@@ -215,13 +254,33 @@ def from_regular_grid_csv(
             return None
 
         grid_data = sorted_values.reshape(tuple(axis.size for axis in axes))
-        return cls(
+        grid_function = cls(
             (axes, grid_data),
             inputs=variable_names,
             outputs=[coeff_name],
             interpolation="regular_grid",
             extrapolation=extrapolation,
         )
+        # Honor the requested interpolation on the grid by rebuilding the
+        # interpolator/extrapolator with the mapped scipy ``method``. The
+        # constructor above always builds the default ("linear"); only rebuild
+        # when a different method was asked for.
+        grid_method = REGULAR_GRID_METHODS.get(interpolation, "linear")
+        smallest_axis = min(axis.size for axis in axes)
+        if smallest_axis < REGULAR_GRID_MIN_POINTS.get(grid_method, 2):
+            warnings.warn(
+                f"Grid interpolation method '{grid_method}' needs at least "
+                f"{REGULAR_GRID_MIN_POINTS[grid_method]} points per axis, but the "
+                f"coarsest axis of '{coeff_name}' has {smallest_axis}; falling "
+                "back to 'linear'.",
+                UserWarning,
+            )
+            grid_method = "linear"
+        if grid_method != "linear":
+            grid_function._grid_method = grid_method
+            grid_function.set_interpolation("regular_grid")
+            grid_function.set_extrapolation(grid_function.get_extrapolation_method())
+        return grid_function
 
     # Define all set methods
     def set_inputs(self, inputs):
@@ -318,6 +377,10 @@ def set_source(self, source):  # pylint: disable=too-many-statements
             self.__dom_dim__ = source.shape[1] - 1
             self._domain = source[:, :-1]
             self._image = source[:, -1]
+            # Cache per-dimension domain bounds so the N-D hot evaluation path
+            # (``__get_value_opt_nd``) does not recompute them on every call.
+            self._domain_min = self._domain.min(axis=0)
+            self._domain_max = self._domain.max(axis=0)
 
             # set x and y. If Function is 2D, also set z
             if self.__dom_dim__ == 1:
@@ -488,11 +551,20 @@ def __process_grid_source(self, source):
                     f"{grid_data.shape[i]} points."
                 )
             if not np.all(np.diff(ax) > 0):
-                warnings.warn(
-                    f"Axis {i} is not strictly sorted in ascending order. "
-                    "RegularGridInterpolator requires sorted axes.",
-                    UserWarning,
-                )
+                # RegularGridInterpolator requires strictly ascending axes. Sort
+                # this axis (and reorder the grid data along it) so descending or
+                # shuffled inputs are accepted; repeated coordinates cannot form
+                # a regular grid and are rejected with a clear error rather than
+                # a cryptic SciPy failure.
+                order = np.argsort(ax, kind="stable")
+                ax = ax[order]
+                grid_data = np.take(grid_data, order, axis=i)
+                axes[i] = ax
+                if not np.all(np.diff(ax) > 0):
+                    raise ValueError(
+                        f"Axis {i} has repeated coordinates; a regular grid "
+                        "requires strictly increasing values along each axis."
+                    )
 
         self._grid_axes = axes
         self._grid_data = grid_data
@@ -596,7 +668,7 @@ def rbf_interpolation(x, x_min, x_max, x_data, y_data, coeffs):  # pylint: disab
                 grid_interpolator = RegularGridInterpolator(
                     self._grid_axes,
                     self._grid_data,
-                    method="linear",
+                    method=getattr(self, "_grid_method", "linear"),
                     bounds_error=True,
                 )
                 # Store so extrapolation funcs can reuse it
@@ -720,9 +792,9 @@ def natural_extrapolation(  # pylint: disable=function-redefined
                         grid_extrapolator = RegularGridInterpolator(
                             self._grid_axes,
                             self._grid_data,
-                            method="linear",
+                            method=getattr(self, "_grid_method", "linear"),
                             bounds_error=False,
-                            fill_value=None,  # linear extrapolation beyond edges
+                            fill_value=None,  # extrapolation beyond edges
                         )
 
                         def natural_extrapolation(  # pylint: disable=function-redefined
@@ -824,8 +896,15 @@ def __get_value_opt_nd(self, *args):
         arg_qty = len(args)
         result = np.empty(arg_qty)
 
-        min_domain = self._domain.T.min(axis=1)
-        max_domain = self._domain.T.max(axis=1)
+        # Domain bounds are fixed once the source is set, so they are cached in
+        # ``set_source`` (this hot path runs per integration step); fall back to
+        # computing them for any Function built without going through it.
+        min_domain = getattr(self, "_domain_min", None)
+        if min_domain is None:
+            min_domain = self._domain.min(axis=0)
+            max_domain = self._domain.max(axis=0)
+        else:
+            max_domain = self._domain_max
 
         lower, upper = args < min_domain, args > max_domain
         extrap = np.logical_or(lower.any(axis=1), upper.any(axis=1))
@@ -4162,7 +4241,7 @@ def to_dict(self, **kwargs):  # pylint: disable=unused-argument
             else:
                 source = source.__name__
 
-        return {
+        function_dict = {
             "source": source,
             "title": self.title,
             "inputs": self.__inputs__,
@@ -4171,6 +4250,20 @@ def to_dict(self, **kwargs):  # pylint: disable=unused-argument
             "extrapolation": self.__extrapolation__,
         }
 
+        # A regular-grid Function cannot be rebuilt from its flat scatter
+        # ``source``; persist the ``(axes, grid_data)`` structure (and the mapped
+        # scipy method) instead, so it round-trips through ``from_dict``.
+        if self.__interpolation__ == "regular_grid":
+            function_dict["source"] = [
+                [np.asarray(axis).tolist() for axis in self._grid_axes],
+                np.asarray(self._grid_data).tolist(),
+            ]
+            grid_method = getattr(self, "_grid_method", "linear")
+            if grid_method != "linear":
+                function_dict["grid_method"] = grid_method
+
+        return function_dict
+
     @classmethod
     def from_dict(cls, func_dict):
         """Creates a Function instance from a dictionary.
@@ -4184,7 +4277,7 @@ def from_dict(cls, func_dict):
         if func_dict["interpolation"] is None and func_dict["extrapolation"] is None:
             source = from_hex_decode(source)
 
-        return cls(
+        function = cls(
             source=source,
             interpolation=func_dict["interpolation"],
             extrapolation=func_dict["extrapolation"],
@@ -4193,6 +4286,16 @@ def from_dict(cls, func_dict):
             title=func_dict["title"],
         )
 
+        # Restore a non-default regular-grid method (the constructor above builds
+        # the "linear" default); rebuild the interpolator/extrapolator with it.
+        grid_method = func_dict.get("grid_method")
+        if grid_method and grid_method != "linear":
+            function._grid_method = grid_method
+            function.set_interpolation("regular_grid")
+            function.set_extrapolation(function.get_extrapolation_method())
+
+        return function
+
     @staticmethod
     def __make_arith_lambda(
         operator, func, other, func_dim, other_dim=0, reverse=False

rocketpy/plots/aero_surface_plots.py

@@ -114,15 +114,14 @@ class _BarrowmanSurfacePlots(_LinearGenericSurfacePlots):
     geometry drawing and the lift-coefficient surface plot."""
 
     def lift(self):
-        """Plots the lift coefficient of the aero surface as a function of Mach
-        and the angle of attack. A 3D plot is expected. See the rocketpy.Function
-        class for more information on how this plot is made.
+        """Plots the lift-curve slope (``clalpha``) of the aero surface as a
+        function of Mach number.
 
         Returns
         -------
         None
         """
-        self.aero_surface.cl()
+        self.aero_surface.clalpha()
 
     def all(self):
         """Plots the surface geometry, the lift coefficient and the
@@ -274,8 +273,7 @@ class for more information on how this plot is made. Also, this method
         -------
         None
         """
-        print("Lift coefficient:")
-        self.aero_surface.cl(filename=filename)
+        print("Lift coefficient derivative:")
         self.aero_surface.clalpha_single_fin(filename=filename)
         self.aero_surface.clalpha_multiple_fins(filename=filename)
 
@@ -356,8 +354,7 @@ class for more information on how this plot is made. Also, this method
         -------
         None
         """
-        print("Lift coefficient:")
-        self.aero_surface.cl(filename=filename)
+        print("Lift coefficient derivative:")
         self.aero_surface.clalpha_single_fin(filename=filename)
 
     def all(self, *, filename=None):

rocketpy/plots/rocket_plots.py

@@ -794,10 +794,10 @@ def _center_of_pressure_range(self, plane, max_angle=np.deg2rad(15), samples=31)
         for angle in angles:
             if plane == "yz":
                 coeffs = rocket.aerodynamic_coefficients_full(0.0, angle, 0.0)
-                force, moment = coeffs["cQ"], coeffs["cn"]
+                force, moment = coeffs["cY"], coeffs["cn"]
             else:
                 coeffs = rocket.aerodynamic_coefficients_full(angle, 0.0, 0.0)
-                force, moment = coeffs["cL"], coeffs["cm"]
+                force, moment = coeffs["cN"], coeffs["cm"]
             if force == 0:
                 continue
             position = cdm + csys * diameter * moment / force