MNT: minor fixes and docstring update.

Author: phmbressan

rocketpy/__init__.py

@@ -1,3 +1,4 @@
+from . import utilities
 from .control import _Controller
 from .environment import Environment, EnvironmentAnalysis
 from .exceptions import (
@@ -68,8 +69,3 @@
     StochasticTail,
     StochasticTrapezoidalFins,
 )
-
-# Imported last: utilities pulls in Environment/Rocket/encoders, which are only
-# fully available once the imports above have run. Exposes
-# ``rocketpy.utilities`` (including ``enable_logging``) on ``import rocketpy``.
-from . import utilities

rocketpy/mathutils/_calc/evaluator.py

@@ -6,9 +6,20 @@
 
 from rocketpy.mathutils._calc.polation_base import PolationBase
 
+_COMPLEX_ND_ERROR = (
+    "Complex coordinates are not supported for N-D array-based Function interpolation."
+)
+
 
 class PolationEvaluator1D(PolationBase):
-    """Route 1D evaluation to interpolation or extrapolation."""
+    """Route 1-D evaluation to interpolation or extrapolation.
+
+    Complex queries are propagated only to support complex-step
+    differentiation of a real 1-D sampled Function. Routing uses the real
+    component so the imaginary perturbation cannot change the selected real
+    interval. This does not define general interpolation over the complex
+    plane.
+    """
 
     __slots__ = (
         "_interpolator",
@@ -38,6 +49,9 @@ def expose(self):
         """Flattens the evaluator into a fast closure for the
         simulation loop.
         """
+        if hasattr(self, "_exposed_fn"):
+            return self._exposed_fn
+
         scalar_eval = self._scalar_fn
         vector_eval = self._vector_fn
 
@@ -53,6 +67,9 @@ def _eval(x, _is_iterable=None):
 
     def expose_scalar(self):
         """Expose a scalar-only evaluator with no iterable checks."""
+        if hasattr(self, "_scalar_fn"):
+            return self._scalar_fn
+
         x_min = self._x_min
         x_max = self._x_max
         interp_eval = self._interpolator.evaluate
@@ -74,6 +91,9 @@ def _scalar(x):
 
     def expose_vector(self):
         """Expose a vector-only evaluator with no iterable checks."""
+        if hasattr(self, "_vector_fn"):
+            return self._vector_fn
+
         x_min = self._x_min
         x_max = self._x_max
         interp_eval = self._interpolator.evaluate
@@ -241,6 +261,9 @@ def evaluate(self, *args, _is_iterable=None):
         return self._exposed_fn(*args, _is_iterable=_is_iterable)
 
     def expose(self):
+        if hasattr(self, "_exposed_fn"):
+            return self._exposed_fn
+
         scalar_eval = self._scalar_fn
         vector_eval = self._vector_fn
 
@@ -256,6 +279,9 @@ def _eval(*args, _is_iterable=None):
         return _eval
 
     def expose_scalar(self):
+        if hasattr(self, "_scalar_fn"):
+            return self._scalar_fn
+
         min_domain = self._min_domain
         max_domain = self._max_domain
         extrapolation_mask = self._extrapolation_mask
@@ -264,7 +290,10 @@ def expose_scalar(self):
         np_local = np
 
         def _scalar(*args):
-            points = np_local.array([args], dtype=float)
+            points = np_local.array([args])
+            if np_local.iscomplexobj(points):
+                raise TypeError(_COMPLEX_ND_ERROR)
+            points = points.astype(float, copy=False)
             point = points[0]
             outside_bounds = ((point < min_domain) | (point > max_domain)).any()
             outside_hull = (
@@ -283,6 +312,9 @@ def _scalar(*args):
         return _scalar
 
     def expose_vector(self):  # pylint: disable=too-many-statements
+        if hasattr(self, "_vector_fn"):
+            return self._vector_fn
+
         min_domain = self._min_domain
         max_domain = self._max_domain
         extrapolation_mask = self._extrapolation_mask
@@ -292,11 +324,11 @@ def expose_vector(self):  # pylint: disable=too-many-statements
 
         def _vector(*args):
             points = np_local.column_stack(args)
-            out_dtype = complex if np_local.iscomplexobj(points) else float
-            result = np_local.empty(len(points), dtype=out_dtype)
-            points_real = points.real
-            lower = points_real < min_domain
-            upper = points_real > max_domain
+            if np_local.iscomplexobj(points):
+                raise TypeError(_COMPLEX_ND_ERROR)
+            result = np_local.empty(len(points), dtype=float)
+            lower = points < min_domain
+            upper = points > max_domain
             extrap_mask = lower.any(axis=1) | upper.any(axis=1)
             if extrapolation_mask is not None:
                 interp_candidates = ~extrap_mask
@@ -352,6 +384,9 @@ def evaluate(self, *args, _is_iterable=None):
         return self._exposed_fn(*args, _is_iterable=_is_iterable)
 
     def expose(self):
+        if hasattr(self, "_exposed_fn"):
+            return self._exposed_fn
+
         scalar_eval = self._scalar_fn
         vector_eval = self._vector_fn
 
@@ -367,6 +402,9 @@ def _eval(*args, _is_iterable=None):
         return _eval
 
     def expose_scalar(self):
+        if hasattr(self, "_scalar_fn"):
+            return self._scalar_fn
+
         min_domain = self._min_domain
         max_domain = self._max_domain
         interp_eval = self._interpolator.evaluate
@@ -376,14 +414,20 @@ def expose_scalar(self):
         if self._interpolator is self._extrapolator:
 
             def _scalar_same(*args):
-                points = np_local.array([args], dtype=float)
+                points = np_local.array([args])
+                if np_local.iscomplexobj(points):
+                    raise TypeError(_COMPLEX_ND_ERROR)
+                points = points.astype(float, copy=False)
                 res = interp_eval(points)
                 return float(res[0])
 
             return _scalar_same
 
         def _scalar(*args):
-            points = np_local.array([args], dtype=float)
+            points = np_local.array([args])
+            if np_local.iscomplexobj(points):
+                raise TypeError(_COMPLEX_ND_ERROR)
+            points = points.astype(float, copy=False)
             point = points[0]
             if ((point < min_domain) | (point > max_domain)).any():
                 return float(extrap_eval(points)[0])
@@ -392,6 +436,9 @@ def _scalar(*args):
         return _scalar
 
     def expose_vector(self):
+        if hasattr(self, "_vector_fn"):
+            return self._vector_fn
+
         min_domain = self._min_domain
         max_domain = self._max_domain
         interp_eval = self._interpolator.evaluate
@@ -402,12 +449,16 @@ def expose_vector(self):
 
             def _vector_same(*args):
                 points = np_local.column_stack(np_local.broadcast_arrays(*args))
+                if np_local.iscomplexobj(points):
+                    raise TypeError(_COMPLEX_ND_ERROR)
                 return interp_eval(points)
 
             return _vector_same
 
         def _vector(*args):
             points = np_local.column_stack(np_local.broadcast_arrays(*args))
+            if np_local.iscomplexobj(points):
+                raise TypeError(_COMPLEX_ND_ERROR)
             result = np_local.empty(len(points), dtype=float)
             lower = points < min_domain
             upper = points > max_domain

rocketpy/mathutils/_calc/polation_1d.py

@@ -1,4 +1,10 @@
-"""1D interpolation and extrapolation strategies."""
+"""1-D interpolation and extrapolation strategies.
+
+The sampled domain and image are real. Complex query values are propagated
+only for complex-step differentiation: their real component selects the
+piecewise interval and the full value is passed to the selected polynomial.
+These strategies do not provide general complex-plane interpolation.
+"""
 
 from __future__ import annotations
 
@@ -30,8 +36,10 @@ def _find_index(
     ----------
     x_arr : np.ndarray
         Sorted 1D array of x-coordinates.
-    xq : float or np.ndarray or complex
-        The query coordinate(s).
+    xq : float, complex or np.ndarray
+        The query coordinate(s). A complex component is meaningful only as an
+        imaginary perturbation for complex-step differentiation. Interval
+        ordering is determined exclusively from the real component.
     n : int
         The size of x_arr.
     _is_iterable : bool, optional
@@ -41,6 +49,12 @@ def _find_index(
     -------
     int or np.ndarray
         The index or indices representing the interval.
+
+    Notes
+    -----
+    Complex numbers have no ordering compatible with the real sampled domain.
+    Using the real component keeps ``x + i*h`` in the same interval as ``x``,
+    which is required for complex-step differentiation.
     """
     if _is_iterable is None:
         _is_iterable = hasattr(xq, "__iter__") and np.ndim(xq) > 0
@@ -50,7 +64,7 @@ def _find_index(
         idx = bisect_left(x_arr, xq.real)
         return 1 if idx < 1 else (idx if idx < n else n - 1)
     else:
-        idx = np.searchsorted(x_arr, xq, side="left")
+        idx = np.searchsorted(x_arr, np.real(xq), side="left")
         return np.clip(idx, 1, n - 1)
 
 
@@ -61,7 +75,7 @@ def _cubic_eval_vec(
     c: float | NDArray[np.float64],
     d: float | NDArray[np.float64],
 ) -> float | NDArray[np.float64]:
-    """Evaluate a cubic polynomial: a + t * (b + t * (c + t * d)).
+    """Evaluate a cubic polynomial: a + b*t + c*t**2 + d*t**3.
 
     Parameters
     ----------

rocketpy/mathutils/_calc/polation_nd.py

@@ -1,4 +1,8 @@
-"""ND interpolation and extrapolation strategies."""
+"""N-D interpolation and extrapolation strategies.
+
+Sampled N-D coordinates are real-valued. Complex query coordinates are not
+supported and are rejected by the evaluator before reaching these strategies.
+"""
 
 from __future__ import annotations