From 72922774f164dae094c1d68795fca11ccc9fae89 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 14 Nov 2025 19:42:53 +0000
Subject: [PATCH 01/11] Add grid interpolation support to Function class with
 from_grid() method

Co-authored-by: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com>
---
 rocketpy/mathutils/function.py             | 258 ++++++++++++++++++++-
 tests/unit/mathutils/test_function_grid.py | 151 ++++++++++++
 2 files changed, 407 insertions(+), 2 deletions(-)
 create mode 100644 tests/unit/mathutils/test_function_grid.py

diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py
index 9fe343687..955ef7c34 100644
--- a/rocketpy/mathutils/function.py
+++ b/rocketpy/mathutils/function.py
@@ -22,6 +22,7 @@
     LinearNDInterpolator,
     NearestNDInterpolator,
     RBFInterpolator,
+    RegularGridInterpolator,
 )
 
 from rocketpy.plots.plot_helpers import show_or_save_plot
@@ -43,6 +44,7 @@
     "spline": 3,
     "shepard": 4,
     "rbf": 5,
+    "linear_grid": 6,
 }
 EXTRAPOLATION_TYPES = {"zero": 0, "natural": 1, "constant": 2}
 
@@ -449,6 +451,37 @@ def rbf_interpolation(x, x_min, x_max, x_data, y_data, coeffs):  # pylint: disab
 
             self._interpolation_func = rbf_interpolation
 
+        elif interpolation == 6:  # linear_grid (RegularGridInterpolator)
+            # For grid interpolation, the actual interpolator is stored separately
+            # This function is a placeholder that should not be called directly
+            # since __get_value_opt_grid is used instead
+            if hasattr(self, '_grid_interpolator'):
+                def grid_interpolation(x, x_min, x_max, x_data, y_data, coeffs):  # pylint: disable=unused-argument
+                    return self._grid_interpolator(x)
+                self._interpolation_func = grid_interpolation
+            else:
+                # Fallback to shepard if grid interpolator not available
+                warnings.warn(
+                    "Grid interpolator not found, falling back to shepard interpolation"
+                )
+                def shepard_fallback(x, x_min, x_max, x_data, y_data, _):
+                    # pylint: disable=unused-argument
+                    arg_qty, arg_dim = x.shape
+                    result = np.empty(arg_qty)
+                    x = x.reshape((arg_qty, 1, arg_dim))
+                    sub_matrix = x_data - x
+                    distances_squared = np.sum(sub_matrix**2, axis=2)
+                    zero_distances = np.where(distances_squared == 0)
+                    valid_indexes = np.ones(arg_qty, dtype=bool)
+                    valid_indexes[zero_distances[0]] = False
+                    weights = distances_squared[valid_indexes] ** (-1.5)
+                    numerator_sum = np.sum(y_data * weights, axis=1)
+                    denominator_sum = np.sum(weights, axis=1)
+                    result[valid_indexes] = numerator_sum / denominator_sum
+                    result[~valid_indexes] = y_data[zero_distances[1]]
+                    return result
+                self._interpolation_func = shepard_fallback
+
         else:
             raise ValueError(f"Interpolation {interpolation} method not recognized.")
 
@@ -635,6 +668,64 @@ def __get_value_opt_nd(self, *args):
 
         return result
 
+    def __get_value_opt_grid(self, *args):
+        """Evaluate the Function using RegularGridInterpolator for structured grids.
+
+        Parameters
+        ----------
+        args : tuple
+            Values where the Function is to be evaluated. Must match the number
+            of dimensions of the grid.
+
+        Returns
+        -------
+        result : scalar or ndarray
+            Value of the Function at the specified points.
+        """
+        # Check if we have the grid interpolator
+        if not hasattr(self, '_grid_interpolator'):
+            raise RuntimeError(
+                "Grid interpolator not initialized. Use from_grid() to create "
+                "a Function with grid interpolation."
+            )
+        
+        # Convert args to appropriate format for RegularGridInterpolator
+        # RegularGridInterpolator expects points as (N, ndim) array
+        if len(args) != self.__dom_dim__:
+            raise ValueError(
+                f"Expected {self.__dom_dim__} arguments but got {len(args)}"
+            )
+        
+        # Handle single point evaluation
+        point = np.array(args).reshape(1, -1)
+        
+        # Handle extrapolation based on the extrapolation setting
+        if self.__extrapolation__ == "constant":
+            # Clamp point to grid boundaries for constant extrapolation
+            for i, axis in enumerate(self._grid_axes):
+                point[0, i] = np.clip(point[0, i], axis[0], axis[-1])
+            result = self._grid_interpolator(point)
+        elif self.__extrapolation__ == "zero":
+            # Check if point is outside bounds
+            outside_bounds = False
+            for i, axis in enumerate(self._grid_axes):
+                if point[0, i] < axis[0] or point[0, i] > axis[-1]:
+                    outside_bounds = True
+                    break
+            if outside_bounds:
+                result = np.array([0.0])
+            else:
+                result = self._grid_interpolator(point)
+        else:
+            # Natural or other extrapolation - use interpolator directly
+            result = self._grid_interpolator(point)
+        
+        # Return scalar for single evaluation
+        if result.size == 1:
+            return float(result[0])
+        
+        return result
+
     def __determine_1d_domain_bounds(self, lower, upper):
         """Determine domain bounds for 1-D function discretization.
 
@@ -3891,11 +3982,11 @@ def __validate_interpolation(self, interpolation):
         elif self.__dom_dim__ > 1:
             if interpolation is None:
                 interpolation = "shepard"
-            if interpolation.lower() not in ["shepard", "linear", "rbf"]:
+            if interpolation.lower() not in ["shepard", "linear", "rbf", "linear_grid"]:
                 warnings.warn(
                     (
                         "Interpolation method set to 'shepard'. The methods "
-                        "'linear', 'shepard' and 'rbf' are supported for "
+                        "'linear', 'shepard', 'rbf' and 'linear_grid' are supported for "
                         "multiple dimensions."
                     ),
                 )
@@ -3950,6 +4041,169 @@ def to_dict(self, **kwargs):  # pylint: disable=unused-argument
             "extrapolation": self.__extrapolation__,
         }
 
+    @classmethod
+    def from_grid(cls, grid_data, axes, inputs=None, outputs=None, 
+                  interpolation="linear_grid", extrapolation="constant", **kwargs):
+        """Creates a Function from N-dimensional grid data.
+        
+        This method is designed for structured grid data, such as CFD simulation
+        results where values are computed on a regular grid. It uses
+        scipy.interpolate.RegularGridInterpolator for efficient interpolation.
+
+        Parameters
+        ----------
+        grid_data : ndarray
+            N-dimensional array containing the function values on the grid.
+            For example, for a 3D function Cd(M, Re, α), this would be a 3D array
+            where grid_data[i, j, k] = Cd(M[i], Re[j], α[k]).
+        axes : list of ndarray
+            List of 1D arrays defining the grid points along each axis.
+            Each array should be sorted in ascending order.
+            For example: [M_axis, Re_axis, alpha_axis].
+        inputs : list of str, optional
+            Names of the input variables. If None, generic names will be used.
+            For example: ['Mach', 'Reynolds', 'Alpha'].
+        outputs : str, optional
+            Name of the output variable. For example: 'Cd'.
+        interpolation : str, optional
+            Interpolation method. Default is 'linear_grid'.
+            Currently only 'linear_grid' is supported for grid data.
+        extrapolation : str, optional
+            Extrapolation behavior. Default is 'constant', which clamps to edge values.
+            'constant': Use nearest edge value for out-of-bounds points.
+            'zero': Return zero for out-of-bounds points.
+        **kwargs : dict, optional
+            Additional arguments passed to the Function constructor.
+
+        Returns
+        -------
+        Function
+            A Function object using RegularGridInterpolator for evaluation.
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> # Create 3D drag coefficient data
+        >>> mach = np.array([0.0, 0.5, 1.0, 1.5, 2.0])
+        >>> reynolds = np.array([1e5, 5e5, 1e6])
+        >>> alpha = np.array([0.0, 2.0, 4.0, 6.0])
+        >>> # Create a simple drag coefficient function
+        >>> M, Re, A = np.meshgrid(mach, reynolds, alpha, indexing='ij')
+        >>> cd_data = 0.3 + 0.1 * M + 1e-7 * Re + 0.01 * A
+        >>> # Create Function object
+        >>> cd_func = Function.from_grid(
+        ...     cd_data, 
+        ...     [mach, reynolds, alpha],
+        ...     inputs=['Mach', 'Reynolds', 'Alpha'],
+        ...     outputs='Cd'
+        ... )
+        >>> # Evaluate at a point
+        >>> cd_func(1.2, 3e5, 3.0)
+        
+        Notes
+        -----
+        - Grid data must be on a regular (structured) grid.
+        - For unstructured data, use the regular Function constructor with
+          scattered points.
+        - Extrapolation with 'constant' mode uses the nearest edge values,
+          which is appropriate for aerodynamic coefficients where extrapolation
+          beyond the data range should be avoided.
+        """
+        # Validate inputs
+        if not isinstance(grid_data, np.ndarray):
+            grid_data = np.array(grid_data)
+        
+        if not isinstance(axes, (list, tuple)):
+            raise ValueError("axes must be a list or tuple of 1D arrays")
+        
+        # Ensure all axes are numpy arrays
+        axes = [np.array(axis) if not isinstance(axis, np.ndarray) else axis 
+                for axis in axes]
+        
+        # Check dimensions match
+        if len(axes) != grid_data.ndim:
+            raise ValueError(
+                f"Number of axes ({len(axes)}) must match grid_data dimensions "
+                f"({grid_data.ndim})"
+            )
+        
+        # Check each axis matches corresponding grid dimension
+        for i, axis in enumerate(axes):
+            if len(axis) != grid_data.shape[i]:
+                raise ValueError(
+                    f"Axis {i} has {len(axis)} points but grid dimension {i} "
+                    f"has {grid_data.shape[i]} points"
+                )
+        
+        # Set default inputs if not provided
+        if inputs is None:
+            inputs = [f"x{i}" for i in range(len(axes))]
+        elif len(inputs) != len(axes):
+            raise ValueError(
+                f"Number of inputs ({len(inputs)}) must match number of axes ({len(axes)})"
+            )
+        
+        # Create a new Function instance
+        func = cls.__new__(cls)
+        
+        # Initialize basic attributes
+        func.source = None  # Will be set to indicate grid source
+        func.__inputs__ = inputs
+        func.__outputs__ = outputs if outputs is not None else "f"
+        func.__interpolation__ = interpolation
+        func.__extrapolation__ = extrapolation
+        func.title = kwargs.get('title', None)
+        func.__img_dim__ = 1
+        func.__cropped_domain__ = (None, None)
+        func._source_type = SourceType.ARRAY
+        func.__dom_dim__ = len(axes)
+        
+        # Store grid-specific data
+        func._grid_axes = axes
+        func._grid_data = grid_data
+        
+        # Create RegularGridInterpolator
+        # We handle extrapolation manually in __get_value_opt_grid,
+        # so we set bounds_error=False and let it extrapolate linearly
+        # (which we'll override when needed)
+        func._grid_interpolator = RegularGridInterpolator(
+            axes, 
+            grid_data, 
+            method='linear',
+            bounds_error=False,
+            fill_value=None  # Linear extrapolation (will be overridden by manual handling)
+        )
+        
+        # Create placeholder domain and image for compatibility
+        # This flattens the grid for any code expecting these attributes
+        mesh = np.meshgrid(*axes, indexing='ij')
+        domain_points = np.column_stack([m.ravel() for m in mesh])
+        func._domain = domain_points
+        func._image = grid_data.ravel()
+        
+        # Set basic array attributes for compatibility
+        func.x_array = axes[0]
+        func.x_initial, func.x_final = axes[0][0], axes[0][-1]
+        func.y_array = func._image[:len(axes[0])]  # Placeholder
+        func.y_initial, func.y_final = func._image[0], func._image[-1]
+        if len(axes) > 2:
+            func.z_array = axes[2] 
+            func.z_initial, func.z_final = axes[2][0], axes[2][-1]
+        
+        # Set get_value_opt to use grid interpolation
+        func.get_value_opt = func.__get_value_opt_grid
+        
+        # Set interpolation and extrapolation functions
+        func.__set_interpolation_func()
+        func.__set_extrapolation_func()
+        
+        # Set inputs and outputs properly
+        func.set_inputs(inputs)
+        func.set_outputs(outputs)
+        func.set_title(func.title)
+        
+        return func
+
     @classmethod
     def from_dict(cls, func_dict):
         """Creates a Function instance from a dictionary.
diff --git a/tests/unit/mathutils/test_function_grid.py b/tests/unit/mathutils/test_function_grid.py
new file mode 100644
index 000000000..ffec452f8
--- /dev/null
+++ b/tests/unit/mathutils/test_function_grid.py
@@ -0,0 +1,151 @@
+"""Unit tests for Function.from_grid() method and grid interpolation."""
+
+import numpy as np
+import pytest
+
+from rocketpy import Function
+
+
+def test_from_grid_1d():
+    """Test from_grid with 1D data (edge case)."""
+    x = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
+    y_data = np.array([0.0, 1.0, 4.0, 9.0, 16.0])  # y = x^2
+    
+    func = Function.from_grid(
+        y_data, 
+        [x],
+        inputs=['x'],
+        outputs='y'
+    )
+    
+    # Test interpolation
+    assert abs(func(1.5) - 2.25) < 0.5  # Should be close to 1.5^2
+
+
+def test_from_grid_2d():
+    """Test from_grid with 2D data."""
+    x = np.array([0.0, 1.0, 2.0])
+    y = np.array([0.0, 1.0, 2.0])
+    
+    # Create grid: f(x, y) = x + 2*y
+    X, Y = np.meshgrid(x, y, indexing='ij')
+    z_data = X + 2 * Y
+    
+    func = Function.from_grid(
+        z_data,
+        [x, y],
+        inputs=['x', 'y'],
+        outputs='z'
+    )
+    
+    # Test exact points
+    assert func(0.0, 0.0) == 0.0
+    assert func(1.0, 1.0) == 3.0
+    assert func(2.0, 2.0) == 6.0
+    
+    # Test interpolation
+    result = func(1.0, 0.5)
+    expected = 1.0 + 2 * 0.5  # = 2.0
+    assert abs(result - expected) < 0.01
+
+
+def test_from_grid_3d_drag_coefficient():
+    """Test from_grid with 3D drag coefficient data (Mach, Reynolds, Alpha)."""
+    # Create sample aerodynamic data
+    mach = np.array([0.0, 0.5, 1.0, 1.5, 2.0])
+    reynolds = np.array([1e5, 5e5, 1e6])
+    alpha = np.array([0.0, 2.0, 4.0, 6.0])
+    
+    # Create a simple drag coefficient model
+    # Cd increases with Mach and alpha, slight dependency on Reynolds
+    M, Re, A = np.meshgrid(mach, reynolds, alpha, indexing='ij')
+    cd_data = 0.3 + 0.1 * M - 1e-7 * Re + 0.01 * A
+    
+    cd_func = Function.from_grid(
+        cd_data,
+        [mach, reynolds, alpha],
+        inputs=['Mach', 'Reynolds', 'Alpha'],
+        outputs='Cd'
+    )
+    
+    # Test at grid points
+    assert abs(cd_func(0.0, 1e5, 0.0) - 0.29) < 0.01  # 0.3 - 1e-7*1e5
+    assert abs(cd_func(1.0, 5e5, 0.0) - 0.35) < 0.01  # 0.3 + 0.1*1 - 1e-7*5e5
+    
+    # Test interpolation between points
+    result = cd_func(0.5, 3e5, 1.0)
+    # Expected roughly: 0.3 + 0.1*0.5 - 1e-7*3e5 + 0.01*1.0 = 0.32
+    assert 0.31 < result < 0.34
+
+
+def test_from_grid_extrapolation_constant():
+    """Test that constant extrapolation clamps to edge values."""
+    x = np.array([0.0, 1.0, 2.0])
+    y = np.array([0.0, 1.0, 4.0])  # y = x^2
+    
+    func = Function.from_grid(
+        y,
+        [x],
+        inputs=['x'],
+        outputs='y',
+        extrapolation='constant'
+    )
+    
+    # Test below lower bound - should return value at x=0
+    assert func(-1.0) == 0.0
+    
+    # Test above upper bound - should return value at x=2
+    assert func(3.0) == 4.0
+
+
+def test_from_grid_validation_errors():
+    """Test that from_grid raises appropriate errors for invalid inputs."""
+    x = np.array([0.0, 1.0, 2.0])
+    y = np.array([0.0, 1.0, 2.0])
+    
+    # Mismatched dimensions
+    X, Y = np.meshgrid(x, y, indexing='ij')
+    z_data = X + Y
+    
+    # Wrong number of axes
+    with pytest.raises(ValueError, match="Number of axes"):
+        Function.from_grid(z_data, [x], inputs=['x'], outputs='z')
+    
+    # Wrong axis length
+    with pytest.raises(ValueError, match="Axis 1 has"):
+        Function.from_grid(z_data, [x, np.array([0.0, 1.0])], inputs=['x', 'y'], outputs='z')
+    
+    # Wrong number of inputs
+    with pytest.raises(ValueError, match="Number of inputs"):
+        Function.from_grid(z_data, [x, y], inputs=['x'], outputs='z')
+
+
+def test_from_grid_default_inputs():
+    """Test that from_grid uses default input names when not provided."""
+    x = np.array([0.0, 1.0, 2.0])
+    y = np.array([0.0, 1.0, 2.0])
+    
+    X, Y = np.meshgrid(x, y, indexing='ij')
+    z_data = X + Y
+    
+    func = Function.from_grid(z_data, [x, y])
+    
+    # Should use default names
+    assert 'x0' in func.__inputs__
+    assert 'x1' in func.__inputs__
+
+
+def test_from_grid_backward_compatibility():
+    """Test that regular Function creation still works after adding from_grid."""
+    # Test 1D function from list
+    func1 = Function([[0, 0], [1, 1], [2, 4], [3, 9]])
+    assert func1(1.5) > 0  # Should interpolate
+    
+    # Test 2D function from array
+    data = np.array([[0, 0, 0], [1, 0, 1], [0, 1, 2], [1, 1, 3]])
+    func2 = Function(data)
+    assert func2(0.5, 0.5) > 0  # Should interpolate
+    
+    # Test callable function
+    func3 = Function(lambda x: x**2)
+    assert func3(2) == 4

From d1b18a156cc35ed233cf64c16cf81ca6f3770d91 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 14 Nov 2025 19:56:24 +0000
Subject: [PATCH 02/11] Add multi-dimensional drag coefficient support to
 Flight class and integration tests

Co-authored-by: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com>
---
 rocketpy/mathutils/function.py          |  28 +++--
 rocketpy/rocket/rocket.py               |  36 +++---
 rocketpy/simulation/flight.py           | 160 ++++++++++++++++++++++--
 tests/integration/test_multidim_drag.py | 122 ++++++++++++++++++
 4 files changed, 308 insertions(+), 38 deletions(-)
 create mode 100644 tests/integration/test_multidim_drag.py

diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py
index 955ef7c34..abdcfaf44 100644
--- a/rocketpy/mathutils/function.py
+++ b/rocketpy/mathutils/function.py
@@ -4146,19 +4146,7 @@ def from_grid(cls, grid_data, axes, inputs=None, outputs=None,
         # Create a new Function instance
         func = cls.__new__(cls)
         
-        # Initialize basic attributes
-        func.source = None  # Will be set to indicate grid source
-        func.__inputs__ = inputs
-        func.__outputs__ = outputs if outputs is not None else "f"
-        func.__interpolation__ = interpolation
-        func.__extrapolation__ = extrapolation
-        func.title = kwargs.get('title', None)
-        func.__img_dim__ = 1
-        func.__cropped_domain__ = (None, None)
-        func._source_type = SourceType.ARRAY
-        func.__dom_dim__ = len(axes)
-        
-        # Store grid-specific data
+        # Store grid-specific data first
         func._grid_axes = axes
         func._grid_data = grid_data
         
@@ -4181,6 +4169,20 @@ def from_grid(cls, grid_data, axes, inputs=None, outputs=None,
         func._domain = domain_points
         func._image = grid_data.ravel()
         
+        # Set source as flattened data array (for compatibility with serialization, etc.)
+        func.source = np.column_stack([domain_points, func._image])
+        
+        # Initialize basic attributes
+        func.__inputs__ = inputs
+        func.__outputs__ = outputs if outputs is not None else "f"
+        func.__interpolation__ = interpolation
+        func.__extrapolation__ = extrapolation
+        func.title = kwargs.get('title', None)
+        func.__img_dim__ = 1
+        func.__cropped_domain__ = (None, None)
+        func._source_type = SourceType.ARRAY
+        func.__dom_dim__ = len(axes)
+        
         # Set basic array attributes for compatibility
         func.x_array = axes[0]
         func.x_initial, func.x_final = axes[0][0], axes[0][-1]
diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py
index 1112a98f3..ffb546eb7 100644
--- a/rocketpy/rocket/rocket.py
+++ b/rocketpy/rocket/rocket.py
@@ -341,20 +341,28 @@ def __init__(  # pylint: disable=too-many-statements
         )
 
         # Define aerodynamic drag coefficients
-        self.power_off_drag = Function(
-            power_off_drag,
-            "Mach Number",
-            "Drag Coefficient with Power Off",
-            "linear",
-            "constant",
-        )
-        self.power_on_drag = Function(
-            power_on_drag,
-            "Mach Number",
-            "Drag Coefficient with Power On",
-            "linear",
-            "constant",
-        )
+        # If already a Function, use it directly (preserves multi-dimensional drag)
+        if isinstance(power_off_drag, Function):
+            self.power_off_drag = power_off_drag
+        else:
+            self.power_off_drag = Function(
+                power_off_drag,
+                "Mach Number",
+                "Drag Coefficient with Power Off",
+                "linear",
+                "constant",
+            )
+        
+        if isinstance(power_on_drag, Function):
+            self.power_on_drag = power_on_drag
+        else:
+            self.power_on_drag = Function(
+                power_on_drag,
+                "Mach Number",
+                "Drag Coefficient with Power On",
+                "linear",
+                "constant",
+            )
 
         # Create a, possibly, temporary empty motor
         # self.motors = Components()  # currently unused, only 1 motor is supported
diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py
index a38be7d93..83a17c500 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -1348,6 +1348,80 @@ def lateral_surface_wind(self):
 
         return -wind_u * np.cos(heading_rad) + wind_v * np.sin(heading_rad)
 
+    def __get_drag_coefficient(
+        self, drag_function, mach, z, velocity_body, freestream_velocity_body
+    ):
+        """Calculate drag coefficient, handling both 1D and multi-dimensional functions.
+
+        Parameters
+        ----------
+        drag_function : Function
+            The drag coefficient function (power_on_drag or power_off_drag)
+        mach : float
+            Mach number
+        z : float
+            Altitude in meters
+        velocity_body : Vector or array-like
+            Rocket velocity in body frame [vx_b, vy_b, vz_b]
+        freestream_velocity_body : Vector or array-like
+            Freestream velocity in body frame [stream_vx_b, stream_vy_b, stream_vz_b]
+
+        Returns
+        -------
+        float
+            Drag coefficient value
+        """
+        # Check if drag function is multi-dimensional by examining its inputs
+        if hasattr(drag_function, "__inputs__") and len(drag_function.__inputs__) > 1:
+            # Multi-dimensional drag function - calculate additional parameters
+            
+            # Calculate Reynolds number
+            # Re = rho * V * L / mu
+            # where L is characteristic length (rocket diameter)
+            rho = self.env.density.get_value_opt(z)
+            mu = self.env.dynamic_viscosity.get_value_opt(z)
+            freestream_speed = np.linalg.norm(freestream_velocity_body)
+            characteristic_length = 2 * self.rocket.radius  # Diameter
+            reynolds = rho * freestream_speed * characteristic_length / mu
+            
+            # Calculate angle of attack
+            # Angle between freestream velocity and rocket axis (z-axis in body frame)
+            # The z component of freestream velocity in body frame
+            if hasattr(freestream_velocity_body, "z"):
+                stream_vz_b = -freestream_velocity_body.z
+            else:
+                stream_vz_b = -freestream_velocity_body[2]
+            
+            # Normalize and calculate angle
+            if freestream_speed > 1e-6:
+                cos_alpha = stream_vz_b / freestream_speed
+                # Clamp to [-1, 1] to avoid numerical issues
+                cos_alpha = np.clip(cos_alpha, -1.0, 1.0)
+                alpha_rad = np.arccos(cos_alpha)
+                alpha_deg = np.rad2deg(alpha_rad)
+            else:
+                alpha_deg = 0.0
+            
+            # Determine which parameters to pass based on input names
+            input_names = [name.lower() for name in drag_function.__inputs__]
+            args = []
+            
+            for name in input_names:
+                if "mach" in name or name == "m":
+                    args.append(mach)
+                elif "reynolds" in name or name == "re":
+                    args.append(reynolds)
+                elif "alpha" in name or name == "a" or "attack" in name:
+                    args.append(alpha_deg)
+                else:
+                    # Unknown parameter, default to mach
+                    args.append(mach)
+            
+            return drag_function.get_value_opt(*args)
+        else:
+            # 1D drag function - only mach number
+            return drag_function.get_value_opt(mach)
+
     def udot_rail1(self, t, u, post_processing=False):
         """Calculates derivative of u state vector with respect to time
         when rocket is flying in 1 DOF motion in the rail.
@@ -1384,7 +1458,37 @@ def udot_rail1(self, t, u, post_processing=False):
             + (vz) ** 2
         ) ** 0.5
         free_stream_mach = free_stream_speed / self.env.speed_of_sound.get_value_opt(z)
-        drag_coeff = self.rocket.power_on_drag.get_value_opt(free_stream_mach)
+        
+        # For rail motion, rocket is constrained - velocity mostly along z-axis in body frame
+        # Calculate velocity in body frame (simplified for rail)
+        a11 = 1 - 2 * (e2**2 + e3**2)
+        a12 = 2 * (e1 * e2 - e0 * e3)
+        a13 = 2 * (e1 * e3 + e0 * e2)
+        a21 = 2 * (e1 * e2 + e0 * e3)
+        a22 = 1 - 2 * (e1**2 + e3**2)
+        a23 = 2 * (e2 * e3 - e0 * e1)
+        a31 = 2 * (e1 * e3 - e0 * e2)
+        a32 = 2 * (e2 * e3 + e0 * e1)
+        a33 = 1 - 2 * (e1**2 + e2**2)
+        
+        vx_b = a11 * vx + a21 * vy + a31 * vz
+        vy_b = a12 * vx + a22 * vy + a32 * vz
+        vz_b = a13 * vx + a23 * vy + a33 * vz
+        
+        # Freestream velocity in body frame
+        wind_vx = self.env.wind_velocity_x.get_value_opt(z)
+        wind_vy = self.env.wind_velocity_y.get_value_opt(z)
+        stream_vx_b = a11 * (wind_vx - vx) + a21 * (wind_vy - vy) + a31 * (-vz)
+        stream_vy_b = a12 * (wind_vx - vx) + a22 * (wind_vy - vy) + a32 * (-vz)
+        stream_vz_b = a13 * (wind_vx - vx) + a23 * (wind_vy - vy) + a33 * (-vz)
+        
+        drag_coeff = self.__get_drag_coefficient(
+            self.rocket.power_on_drag,
+            free_stream_mach,
+            z,
+            [vx_b, vy_b, vz_b],
+            [stream_vx_b, stream_vy_b, stream_vz_b],
+        )
 
         # Calculate Forces
         pressure = self.env.pressure.get_value_opt(z)
@@ -1552,12 +1656,34 @@ def u_dot(self, t, u, post_processing=False):  # pylint: disable=too-many-locals
         ) ** 0.5
         free_stream_mach = free_stream_speed / speed_of_sound
 
+        # Get rocket velocity in body frame (needed for drag calculation)
+        vx_b = a11 * vx + a21 * vy + a31 * vz
+        vy_b = a12 * vx + a22 * vy + a32 * vz
+        vz_b = a13 * vx + a23 * vy + a33 * vz
+        
+        # Calculate freestream velocity in body frame
+        stream_vx_b = a11 * (wind_velocity_x - vx) + a21 * (wind_velocity_y - vy) + a31 * (-vz)
+        stream_vy_b = a12 * (wind_velocity_x - vx) + a22 * (wind_velocity_y - vy) + a32 * (-vz)
+        stream_vz_b = a13 * (wind_velocity_x - vx) + a23 * (wind_velocity_y - vy) + a33 * (-vz)
+        
         # Determine aerodynamics forces
         # Determine Drag Force
         if t < self.rocket.motor.burn_out_time:
-            drag_coeff = self.rocket.power_on_drag.get_value_opt(free_stream_mach)
+            drag_coeff = self.__get_drag_coefficient(
+                self.rocket.power_on_drag,
+                free_stream_mach,
+                z,
+                [vx_b, vy_b, vz_b],
+                [stream_vx_b, stream_vy_b, stream_vz_b],
+            )
         else:
-            drag_coeff = self.rocket.power_off_drag.get_value_opt(free_stream_mach)
+            drag_coeff = self.__get_drag_coefficient(
+                self.rocket.power_off_drag,
+                free_stream_mach,
+                z,
+                [vx_b, vy_b, vz_b],
+                [stream_vx_b, stream_vy_b, stream_vz_b],
+            )
         rho = self.env.density.get_value_opt(z)
         R3 = -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff
         for air_brakes in self.rocket.air_brakes:
@@ -1579,10 +1705,6 @@ def u_dot(self, t, u, post_processing=False):  # pylint: disable=too-many-locals
         # Off center moment
         M1 += self.rocket.cp_eccentricity_y * R3
         M2 -= self.rocket.cp_eccentricity_x * R3
-        # Get rocket velocity in body frame
-        vx_b = a11 * vx + a21 * vy + a31 * vz
-        vy_b = a12 * vx + a22 * vy + a32 * vz
-        vz_b = a13 * vx + a23 * vy + a33 * vz
         # Calculate lift and moment for each component of the rocket
         velocity_in_body_frame = Vector([vx_b, vy_b, vz_b])
         w = Vector([omega1, omega2, omega3])
@@ -1831,6 +1953,12 @@ def u_dot_generalized(self, t, u, post_processing=False):  # pylint: disable=too
         speed_of_sound = self.env.speed_of_sound.get_value_opt(z)
         free_stream_mach = free_stream_speed / speed_of_sound
 
+        # Get rocket velocity in body frame (needed for drag calculation)
+        velocity_in_body_frame = Kt @ v
+        # Calculate freestream velocity in body frame
+        freestream_velocity = wind_velocity - v
+        freestream_velocity_body = Kt @ freestream_velocity
+        
         if self.rocket.motor.burn_start_time < t < self.rocket.motor.burn_out_time:
             pressure = self.env.pressure.get_value_opt(z)
             net_thrust = max(
@@ -1838,10 +1966,22 @@ def u_dot_generalized(self, t, u, post_processing=False):  # pylint: disable=too
                 + self.rocket.motor.pressure_thrust(pressure),
                 0,
             )
-            drag_coeff = self.rocket.power_on_drag.get_value_opt(free_stream_mach)
+            drag_coeff = self.__get_drag_coefficient(
+                self.rocket.power_on_drag,
+                free_stream_mach,
+                z,
+                velocity_in_body_frame,
+                freestream_velocity_body,
+            )
         else:
             net_thrust = 0
-            drag_coeff = self.rocket.power_off_drag.get_value_opt(free_stream_mach)
+            drag_coeff = self.__get_drag_coefficient(
+                self.rocket.power_off_drag,
+                free_stream_mach,
+                z,
+                velocity_in_body_frame,
+                freestream_velocity_body,
+            )
         R3 += -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff
         for air_brakes in self.rocket.air_brakes:
             if air_brakes.deployment_level > 0:
@@ -1859,8 +1999,6 @@ def u_dot_generalized(self, t, u, post_processing=False):  # pylint: disable=too
                     R3 = air_brakes_force  # Substitutes rocket drag coefficient
                 else:
                     R3 += air_brakes_force
-        # Get rocket velocity in body frame
-        velocity_in_body_frame = Kt @ v
         # Calculate lift and moment for each component of the rocket
         for aero_surface, _ in self.rocket.aerodynamic_surfaces:
             # Component cp relative to CDM in body frame
diff --git a/tests/integration/test_multidim_drag.py b/tests/integration/test_multidim_drag.py
new file mode 100644
index 000000000..90ef58467
--- /dev/null
+++ b/tests/integration/test_multidim_drag.py
@@ -0,0 +1,122 @@
+"""Integration tests for multi-dimensional drag coefficient support."""
+
+import numpy as np
+import pytest
+
+from rocketpy import Function
+
+
+def test_flight_with_1d_drag(example_plain_env, calisto):
+    """Test that flights with 1D drag curves still work (backward compatibility)."""
+    from rocketpy import Flight
+    
+    flight = Flight(
+        rocket=calisto,
+        environment=example_plain_env,
+        rail_length=5.2,
+        inclination=85,
+        heading=0,
+    )
+    
+    # Check that flight completed successfully
+    assert flight.t_final > 0
+    assert flight.apogee > 0
+    assert flight.apogee_time > 0
+
+
+def test_flight_with_3d_drag_basic():
+    """Test that a simple 3D drag function works."""
+    from rocketpy import Environment, Flight, Rocket, SolidMotor
+    
+    # Create environment
+    env = Environment(gravity=9.81)
+    env.set_atmospheric_model(type='standard_atmosphere')
+    
+    # Create motor with simple constant thrust
+    motor = SolidMotor(
+        thrust_source=lambda t: 2000 if t < 3 else 0,
+        burn_time=3.0,
+        nozzle_radius=0.033,
+        dry_mass=1.815,
+        dry_inertia=(0.125, 0.125, 0.002),
+        center_of_dry_mass_position=0.317,
+        grains_center_of_mass_position=0.397,
+        grain_number=5,
+        grain_separation=0.005,
+        grain_density=1815,
+        grain_outer_radius=0.033,
+        grain_initial_inner_radius=0.015,
+        grain_initial_height=0.120,
+        nozzle_position=0,
+        throat_radius=0.011,
+    )
+    
+    # Create 3D drag
+    mach = np.array([0.0, 0.5, 1.0, 1.5, 2.0])
+    reynolds = np.array([1e5, 5e5, 1e6])
+    alpha = np.array([0.0, 2.0, 4.0, 6.0])
+    
+    M, Re, A = np.meshgrid(mach, reynolds, alpha, indexing='ij')
+    cd_data = 0.3 + 0.1 * M - 1e-7 * Re + 0.01 * A
+    cd_data = np.clip(cd_data, 0.2, 1.0)
+    
+    power_off_drag = Function.from_grid(cd_data, [mach, reynolds, alpha], 
+                                 inputs=['Mach', 'Reynolds', 'Alpha'], outputs='Cd')
+    power_on_drag = Function.from_grid(cd_data * 1.1, [mach, reynolds, alpha], 
+                                inputs=['Mach', 'Reynolds', 'Alpha'], outputs='Cd')
+    
+    # Create rocket
+    rocket = Rocket(
+        radius=0.0635,
+        mass=16.24,
+        inertia=(6.321, 6.321, 0.034),
+        power_off_drag=power_off_drag,
+        power_on_drag=power_on_drag,
+        center_of_mass_without_motor=0,
+        coordinate_system_orientation='tail_to_nose',
+    )
+    rocket.set_rail_buttons(0.2, -0.5, 30)
+    rocket.add_motor(motor, position=-1.255)
+    
+    # Run flight
+    flight = Flight(
+        rocket=rocket,
+        environment=env,
+        rail_length=5.2,
+        inclination=85,
+        heading=0,
+    )
+    
+    # Check results - should launch and have non-zero apogee
+    assert flight.apogee > 100, f"Apogee too low: {flight.apogee}m"
+    assert flight.apogee < 5000, f"Apogee too high: {flight.apogee}m"
+    assert hasattr(flight, 'angle_of_attack')
+
+
+def test_3d_drag_with_varying_alpha():
+    """Test that 3D drag responds to angle of attack changes."""
+    # Create drag function with strong alpha dependency
+    mach = np.array([0.0, 0.5, 1.0, 1.5])
+    reynolds = np.array([1e5, 1e6])
+    alpha = np.array([0.0, 5.0, 10.0, 15.0])
+    
+    M, Re, A = np.meshgrid(mach, reynolds, alpha, indexing='ij')
+    # Strong alpha dependency: Cd increases significantly with alpha
+    cd_data = 0.3 + 0.05 * M + 0.03 * A
+    cd_data = np.clip(cd_data, 0.2, 2.0)
+    
+    drag_func = Function.from_grid(
+        cd_data,
+        [mach, reynolds, alpha],
+        inputs=['Mach', 'Reynolds', 'Alpha'],
+        outputs='Cd'
+    )
+    
+    # Test at different angles of attack
+    # At zero alpha, Cd should be lower
+    cd_0 = drag_func(0.8, 5e5, 0.0)
+    cd_10 = drag_func(0.8, 5e5, 10.0)
+    
+    # Cd should increase with alpha
+    assert cd_10 > cd_0
+    assert cd_10 - cd_0 > 0.2  # Should show significant difference

From b621a16ca154111ae93772c8f03f0744a2cddb6c Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 14 Nov 2025 19:57:51 +0000
Subject: [PATCH 03/11] Run ruff format on modified files

Co-authored-by: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com>
---
 rocketpy/mathutils/function.py             | 90 +++++++++++++---------
 rocketpy/rocket/rocket.py                  |  2 +-
 rocketpy/simulation/flight.py              | 38 +++++----
 tests/integration/test_multidim_drag.py    | 58 ++++++++------
 tests/unit/mathutils/test_function_grid.py | 88 +++++++++------------
 5 files changed, 146 insertions(+), 130 deletions(-)

diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py
index abdcfaf44..245640551 100644
--- a/rocketpy/mathutils/function.py
+++ b/rocketpy/mathutils/function.py
@@ -455,15 +455,18 @@ def rbf_interpolation(x, x_min, x_max, x_data, y_data, coeffs):  # pylint: disab
             # For grid interpolation, the actual interpolator is stored separately
             # This function is a placeholder that should not be called directly
             # since __get_value_opt_grid is used instead
-            if hasattr(self, '_grid_interpolator'):
+            if hasattr(self, "_grid_interpolator"):
+
                 def grid_interpolation(x, x_min, x_max, x_data, y_data, coeffs):  # pylint: disable=unused-argument
                     return self._grid_interpolator(x)
+
                 self._interpolation_func = grid_interpolation
             else:
                 # Fallback to shepard if grid interpolator not available
                 warnings.warn(
                     "Grid interpolator not found, falling back to shepard interpolation"
                 )
+
                 def shepard_fallback(x, x_min, x_max, x_data, y_data, _):
                     # pylint: disable=unused-argument
                     arg_qty, arg_dim = x.shape
@@ -480,6 +483,7 @@ def shepard_fallback(x, x_min, x_max, x_data, y_data, _):
                     result[valid_indexes] = numerator_sum / denominator_sum
                     result[~valid_indexes] = y_data[zero_distances[1]]
                     return result
+
                 self._interpolation_func = shepard_fallback
 
         else:
@@ -683,22 +687,22 @@ def __get_value_opt_grid(self, *args):
             Value of the Function at the specified points.
         """
         # Check if we have the grid interpolator
-        if not hasattr(self, '_grid_interpolator'):
+        if not hasattr(self, "_grid_interpolator"):
             raise RuntimeError(
                 "Grid interpolator not initialized. Use from_grid() to create "
                 "a Function with grid interpolation."
             )
-        
+
         # Convert args to appropriate format for RegularGridInterpolator
         # RegularGridInterpolator expects points as (N, ndim) array
         if len(args) != self.__dom_dim__:
             raise ValueError(
                 f"Expected {self.__dom_dim__} arguments but got {len(args)}"
             )
-        
+
         # Handle single point evaluation
         point = np.array(args).reshape(1, -1)
-        
+
         # Handle extrapolation based on the extrapolation setting
         if self.__extrapolation__ == "constant":
             # Clamp point to grid boundaries for constant extrapolation
@@ -719,11 +723,11 @@ def __get_value_opt_grid(self, *args):
         else:
             # Natural or other extrapolation - use interpolator directly
             result = self._grid_interpolator(point)
-        
+
         # Return scalar for single evaluation
         if result.size == 1:
             return float(result[0])
-        
+
         return result
 
     def __determine_1d_domain_bounds(self, lower, upper):
@@ -4042,10 +4046,18 @@ def to_dict(self, **kwargs):  # pylint: disable=unused-argument
         }
 
     @classmethod
-    def from_grid(cls, grid_data, axes, inputs=None, outputs=None, 
-                  interpolation="linear_grid", extrapolation="constant", **kwargs):
+    def from_grid(
+        cls,
+        grid_data,
+        axes,
+        inputs=None,
+        outputs=None,
+        interpolation="linear_grid",
+        extrapolation="constant",
+        **kwargs,
+    ):
         """Creates a Function from N-dimensional grid data.
-        
+
         This method is designed for structured grid data, such as CFD simulation
         results where values are computed on a regular grid. It uses
         scipy.interpolate.RegularGridInterpolator for efficient interpolation.
@@ -4092,14 +4104,14 @@ def from_grid(cls, grid_data, axes, inputs=None, outputs=None,
         >>> cd_data = 0.3 + 0.1 * M + 1e-7 * Re + 0.01 * A
         >>> # Create Function object
         >>> cd_func = Function.from_grid(
-        ...     cd_data, 
+        ...     cd_data,
         ...     [mach, reynolds, alpha],
         ...     inputs=['Mach', 'Reynolds', 'Alpha'],
         ...     outputs='Cd'
         ... )
         >>> # Evaluate at a point
         >>> cd_func(1.2, 3e5, 3.0)
-        
+
         Notes
         -----
         - Grid data must be on a regular (structured) grid.
@@ -4112,21 +4124,23 @@ def from_grid(cls, grid_data, axes, inputs=None, outputs=None,
         # Validate inputs
         if not isinstance(grid_data, np.ndarray):
             grid_data = np.array(grid_data)
-        
+
         if not isinstance(axes, (list, tuple)):
             raise ValueError("axes must be a list or tuple of 1D arrays")
-        
+
         # Ensure all axes are numpy arrays
-        axes = [np.array(axis) if not isinstance(axis, np.ndarray) else axis 
-                for axis in axes]
-        
+        axes = [
+            np.array(axis) if not isinstance(axis, np.ndarray) else axis
+            for axis in axes
+        ]
+
         # Check dimensions match
         if len(axes) != grid_data.ndim:
             raise ValueError(
                 f"Number of axes ({len(axes)}) must match grid_data dimensions "
                 f"({grid_data.ndim})"
             )
-        
+
         # Check each axis matches corresponding grid dimension
         for i, axis in enumerate(axes):
             if len(axis) != grid_data.shape[i]:
@@ -4134,7 +4148,7 @@ def from_grid(cls, grid_data, axes, inputs=None, outputs=None,
                     f"Axis {i} has {len(axis)} points but grid dimension {i} "
                     f"has {grid_data.shape[i]} points"
                 )
-        
+
         # Set default inputs if not provided
         if inputs is None:
             inputs = [f"x{i}" for i in range(len(axes))]
@@ -4142,68 +4156,68 @@ def from_grid(cls, grid_data, axes, inputs=None, outputs=None,
             raise ValueError(
                 f"Number of inputs ({len(inputs)}) must match number of axes ({len(axes)})"
             )
-        
+
         # Create a new Function instance
         func = cls.__new__(cls)
-        
+
         # Store grid-specific data first
         func._grid_axes = axes
         func._grid_data = grid_data
-        
+
         # Create RegularGridInterpolator
         # We handle extrapolation manually in __get_value_opt_grid,
         # so we set bounds_error=False and let it extrapolate linearly
         # (which we'll override when needed)
         func._grid_interpolator = RegularGridInterpolator(
-            axes, 
-            grid_data, 
-            method='linear',
+            axes,
+            grid_data,
+            method="linear",
             bounds_error=False,
-            fill_value=None  # Linear extrapolation (will be overridden by manual handling)
+            fill_value=None,  # Linear extrapolation (will be overridden by manual handling)
         )
-        
+
         # Create placeholder domain and image for compatibility
         # This flattens the grid for any code expecting these attributes
-        mesh = np.meshgrid(*axes, indexing='ij')
+        mesh = np.meshgrid(*axes, indexing="ij")
         domain_points = np.column_stack([m.ravel() for m in mesh])
         func._domain = domain_points
         func._image = grid_data.ravel()
-        
+
         # Set source as flattened data array (for compatibility with serialization, etc.)
         func.source = np.column_stack([domain_points, func._image])
-        
+
         # Initialize basic attributes
         func.__inputs__ = inputs
         func.__outputs__ = outputs if outputs is not None else "f"
         func.__interpolation__ = interpolation
         func.__extrapolation__ = extrapolation
-        func.title = kwargs.get('title', None)
+        func.title = kwargs.get("title", None)
         func.__img_dim__ = 1
         func.__cropped_domain__ = (None, None)
         func._source_type = SourceType.ARRAY
         func.__dom_dim__ = len(axes)
-        
+
         # Set basic array attributes for compatibility
         func.x_array = axes[0]
         func.x_initial, func.x_final = axes[0][0], axes[0][-1]
-        func.y_array = func._image[:len(axes[0])]  # Placeholder
+        func.y_array = func._image[: len(axes[0])]  # Placeholder
         func.y_initial, func.y_final = func._image[0], func._image[-1]
         if len(axes) > 2:
-            func.z_array = axes[2] 
+            func.z_array = axes[2]
             func.z_initial, func.z_final = axes[2][0], axes[2][-1]
-        
+
         # Set get_value_opt to use grid interpolation
         func.get_value_opt = func.__get_value_opt_grid
-        
+
         # Set interpolation and extrapolation functions
         func.__set_interpolation_func()
         func.__set_extrapolation_func()
-        
+
         # Set inputs and outputs properly
         func.set_inputs(inputs)
         func.set_outputs(outputs)
         func.set_title(func.title)
-        
+
         return func
 
     @classmethod
diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py
index ffb546eb7..f5b3f13c1 100644
--- a/rocketpy/rocket/rocket.py
+++ b/rocketpy/rocket/rocket.py
@@ -352,7 +352,7 @@ def __init__(  # pylint: disable=too-many-statements
                 "linear",
                 "constant",
             )
-        
+
         if isinstance(power_on_drag, Function):
             self.power_on_drag = power_on_drag
         else:
diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py
index 83a17c500..83059f246 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -1374,7 +1374,7 @@ def __get_drag_coefficient(
         # Check if drag function is multi-dimensional by examining its inputs
         if hasattr(drag_function, "__inputs__") and len(drag_function.__inputs__) > 1:
             # Multi-dimensional drag function - calculate additional parameters
-            
+
             # Calculate Reynolds number
             # Re = rho * V * L / mu
             # where L is characteristic length (rocket diameter)
@@ -1383,7 +1383,7 @@ def __get_drag_coefficient(
             freestream_speed = np.linalg.norm(freestream_velocity_body)
             characteristic_length = 2 * self.rocket.radius  # Diameter
             reynolds = rho * freestream_speed * characteristic_length / mu
-            
+
             # Calculate angle of attack
             # Angle between freestream velocity and rocket axis (z-axis in body frame)
             # The z component of freestream velocity in body frame
@@ -1391,7 +1391,7 @@ def __get_drag_coefficient(
                 stream_vz_b = -freestream_velocity_body.z
             else:
                 stream_vz_b = -freestream_velocity_body[2]
-            
+
             # Normalize and calculate angle
             if freestream_speed > 1e-6:
                 cos_alpha = stream_vz_b / freestream_speed
@@ -1401,11 +1401,11 @@ def __get_drag_coefficient(
                 alpha_deg = np.rad2deg(alpha_rad)
             else:
                 alpha_deg = 0.0
-            
+
             # Determine which parameters to pass based on input names
             input_names = [name.lower() for name in drag_function.__inputs__]
             args = []
-            
+
             for name in input_names:
                 if "mach" in name or name == "m":
                     args.append(mach)
@@ -1416,7 +1416,7 @@ def __get_drag_coefficient(
                 else:
                     # Unknown parameter, default to mach
                     args.append(mach)
-            
+
             return drag_function.get_value_opt(*args)
         else:
             # 1D drag function - only mach number
@@ -1458,7 +1458,7 @@ def udot_rail1(self, t, u, post_processing=False):
             + (vz) ** 2
         ) ** 0.5
         free_stream_mach = free_stream_speed / self.env.speed_of_sound.get_value_opt(z)
-        
+
         # For rail motion, rocket is constrained - velocity mostly along z-axis in body frame
         # Calculate velocity in body frame (simplified for rail)
         a11 = 1 - 2 * (e2**2 + e3**2)
@@ -1470,18 +1470,18 @@ def udot_rail1(self, t, u, post_processing=False):
         a31 = 2 * (e1 * e3 - e0 * e2)
         a32 = 2 * (e2 * e3 + e0 * e1)
         a33 = 1 - 2 * (e1**2 + e2**2)
-        
+
         vx_b = a11 * vx + a21 * vy + a31 * vz
         vy_b = a12 * vx + a22 * vy + a32 * vz
         vz_b = a13 * vx + a23 * vy + a33 * vz
-        
+
         # Freestream velocity in body frame
         wind_vx = self.env.wind_velocity_x.get_value_opt(z)
         wind_vy = self.env.wind_velocity_y.get_value_opt(z)
         stream_vx_b = a11 * (wind_vx - vx) + a21 * (wind_vy - vy) + a31 * (-vz)
         stream_vy_b = a12 * (wind_vx - vx) + a22 * (wind_vy - vy) + a32 * (-vz)
         stream_vz_b = a13 * (wind_vx - vx) + a23 * (wind_vy - vy) + a33 * (-vz)
-        
+
         drag_coeff = self.__get_drag_coefficient(
             self.rocket.power_on_drag,
             free_stream_mach,
@@ -1660,12 +1660,18 @@ def u_dot(self, t, u, post_processing=False):  # pylint: disable=too-many-locals
         vx_b = a11 * vx + a21 * vy + a31 * vz
         vy_b = a12 * vx + a22 * vy + a32 * vz
         vz_b = a13 * vx + a23 * vy + a33 * vz
-        
+
         # Calculate freestream velocity in body frame
-        stream_vx_b = a11 * (wind_velocity_x - vx) + a21 * (wind_velocity_y - vy) + a31 * (-vz)
-        stream_vy_b = a12 * (wind_velocity_x - vx) + a22 * (wind_velocity_y - vy) + a32 * (-vz)
-        stream_vz_b = a13 * (wind_velocity_x - vx) + a23 * (wind_velocity_y - vy) + a33 * (-vz)
-        
+        stream_vx_b = (
+            a11 * (wind_velocity_x - vx) + a21 * (wind_velocity_y - vy) + a31 * (-vz)
+        )
+        stream_vy_b = (
+            a12 * (wind_velocity_x - vx) + a22 * (wind_velocity_y - vy) + a32 * (-vz)
+        )
+        stream_vz_b = (
+            a13 * (wind_velocity_x - vx) + a23 * (wind_velocity_y - vy) + a33 * (-vz)
+        )
+
         # Determine aerodynamics forces
         # Determine Drag Force
         if t < self.rocket.motor.burn_out_time:
@@ -1958,7 +1964,7 @@ def u_dot_generalized(self, t, u, post_processing=False):  # pylint: disable=too
         # Calculate freestream velocity in body frame
         freestream_velocity = wind_velocity - v
         freestream_velocity_body = Kt @ freestream_velocity
-        
+
         if self.rocket.motor.burn_start_time < t < self.rocket.motor.burn_out_time:
             pressure = self.env.pressure.get_value_opt(z)
             net_thrust = max(
diff --git a/tests/integration/test_multidim_drag.py b/tests/integration/test_multidim_drag.py
index 90ef58467..accf20afa 100644
--- a/tests/integration/test_multidim_drag.py
+++ b/tests/integration/test_multidim_drag.py
@@ -9,7 +9,7 @@
 def test_flight_with_1d_drag(example_plain_env, calisto):
     """Test that flights with 1D drag curves still work (backward compatibility)."""
     from rocketpy import Flight
-    
+
     flight = Flight(
         rocket=calisto,
         environment=example_plain_env,
@@ -17,7 +17,7 @@ def test_flight_with_1d_drag(example_plain_env, calisto):
         inclination=85,
         heading=0,
     )
-    
+
     # Check that flight completed successfully
     assert flight.t_final > 0
     assert flight.apogee > 0
@@ -27,11 +27,11 @@ def test_flight_with_1d_drag(example_plain_env, calisto):
 def test_flight_with_3d_drag_basic():
     """Test that a simple 3D drag function works."""
     from rocketpy import Environment, Flight, Rocket, SolidMotor
-    
+
     # Create environment
     env = Environment(gravity=9.81)
-    env.set_atmospheric_model(type='standard_atmosphere')
-    
+    env.set_atmospheric_model(type="standard_atmosphere")
+
     # Create motor with simple constant thrust
     motor = SolidMotor(
         thrust_source=lambda t: 2000 if t < 3 else 0,
@@ -50,21 +50,29 @@ def test_flight_with_3d_drag_basic():
         nozzle_position=0,
         throat_radius=0.011,
     )
-    
+
     # Create 3D drag
     mach = np.array([0.0, 0.5, 1.0, 1.5, 2.0])
     reynolds = np.array([1e5, 5e5, 1e6])
     alpha = np.array([0.0, 2.0, 4.0, 6.0])
-    
-    M, Re, A = np.meshgrid(mach, reynolds, alpha, indexing='ij')
+
+    M, Re, A = np.meshgrid(mach, reynolds, alpha, indexing="ij")
     cd_data = 0.3 + 0.1 * M - 1e-7 * Re + 0.01 * A
     cd_data = np.clip(cd_data, 0.2, 1.0)
-    
-    power_off_drag = Function.from_grid(cd_data, [mach, reynolds, alpha], 
-                                 inputs=['Mach', 'Reynolds', 'Alpha'], outputs='Cd')
-    power_on_drag = Function.from_grid(cd_data * 1.1, [mach, reynolds, alpha], 
-                                inputs=['Mach', 'Reynolds', 'Alpha'], outputs='Cd')
-    
+
+    power_off_drag = Function.from_grid(
+        cd_data,
+        [mach, reynolds, alpha],
+        inputs=["Mach", "Reynolds", "Alpha"],
+        outputs="Cd",
+    )
+    power_on_drag = Function.from_grid(
+        cd_data * 1.1,
+        [mach, reynolds, alpha],
+        inputs=["Mach", "Reynolds", "Alpha"],
+        outputs="Cd",
+    )
+
     # Create rocket
     rocket = Rocket(
         radius=0.0635,
@@ -73,11 +81,11 @@ def test_flight_with_3d_drag_basic():
         power_off_drag=power_off_drag,
         power_on_drag=power_on_drag,
         center_of_mass_without_motor=0,
-        coordinate_system_orientation='tail_to_nose',
+        coordinate_system_orientation="tail_to_nose",
     )
     rocket.set_rail_buttons(0.2, -0.5, 30)
     rocket.add_motor(motor, position=-1.255)
-    
+
     # Run flight
     flight = Flight(
         rocket=rocket,
@@ -86,11 +94,11 @@ def test_flight_with_3d_drag_basic():
         inclination=85,
         heading=0,
     )
-    
+
     # Check results - should launch and have non-zero apogee
     assert flight.apogee > 100, f"Apogee too low: {flight.apogee}m"
     assert flight.apogee < 5000, f"Apogee too high: {flight.apogee}m"
-    assert hasattr(flight, 'angle_of_attack')
+    assert hasattr(flight, "angle_of_attack")
 
 
 def test_3d_drag_with_varying_alpha():
@@ -99,24 +107,24 @@ def test_3d_drag_with_varying_alpha():
     mach = np.array([0.0, 0.5, 1.0, 1.5])
     reynolds = np.array([1e5, 1e6])
     alpha = np.array([0.0, 5.0, 10.0, 15.0])
-    
-    M, Re, A = np.meshgrid(mach, reynolds, alpha, indexing='ij')
+
+    M, Re, A = np.meshgrid(mach, reynolds, alpha, indexing="ij")
     # Strong alpha dependency: Cd increases significantly with alpha
     cd_data = 0.3 + 0.05 * M + 0.03 * A
     cd_data = np.clip(cd_data, 0.2, 2.0)
-    
+
     drag_func = Function.from_grid(
         cd_data,
         [mach, reynolds, alpha],
-        inputs=['Mach', 'Reynolds', 'Alpha'],
-        outputs='Cd'
+        inputs=["Mach", "Reynolds", "Alpha"],
+        outputs="Cd",
     )
-    
+
     # Test at different angles of attack
     # At zero alpha, Cd should be lower
     cd_0 = drag_func(0.8, 5e5, 0.0)
     cd_10 = drag_func(0.8, 5e5, 10.0)
-    
+
     # Cd should increase with alpha
     assert cd_10 > cd_0
     assert cd_10 - cd_0 > 0.2  # Should show significant difference
diff --git a/tests/unit/mathutils/test_function_grid.py b/tests/unit/mathutils/test_function_grid.py
index ffec452f8..994689544 100644
--- a/tests/unit/mathutils/test_function_grid.py
+++ b/tests/unit/mathutils/test_function_grid.py
@@ -10,14 +10,9 @@ def test_from_grid_1d():
     """Test from_grid with 1D data (edge case)."""
     x = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
     y_data = np.array([0.0, 1.0, 4.0, 9.0, 16.0])  # y = x^2
-    
-    func = Function.from_grid(
-        y_data, 
-        [x],
-        inputs=['x'],
-        outputs='y'
-    )
-    
+
+    func = Function.from_grid(y_data, [x], inputs=["x"], outputs="y")
+
     # Test interpolation
     assert abs(func(1.5) - 2.25) < 0.5  # Should be close to 1.5^2
 
@@ -26,23 +21,18 @@ def test_from_grid_2d():
     """Test from_grid with 2D data."""
     x = np.array([0.0, 1.0, 2.0])
     y = np.array([0.0, 1.0, 2.0])
-    
+
     # Create grid: f(x, y) = x + 2*y
-    X, Y = np.meshgrid(x, y, indexing='ij')
+    X, Y = np.meshgrid(x, y, indexing="ij")
     z_data = X + 2 * Y
-    
-    func = Function.from_grid(
-        z_data,
-        [x, y],
-        inputs=['x', 'y'],
-        outputs='z'
-    )
-    
+
+    func = Function.from_grid(z_data, [x, y], inputs=["x", "y"], outputs="z")
+
     # Test exact points
     assert func(0.0, 0.0) == 0.0
     assert func(1.0, 1.0) == 3.0
     assert func(2.0, 2.0) == 6.0
-    
+
     # Test interpolation
     result = func(1.0, 0.5)
     expected = 1.0 + 2 * 0.5  # = 2.0
@@ -55,23 +45,23 @@ def test_from_grid_3d_drag_coefficient():
     mach = np.array([0.0, 0.5, 1.0, 1.5, 2.0])
     reynolds = np.array([1e5, 5e5, 1e6])
     alpha = np.array([0.0, 2.0, 4.0, 6.0])
-    
+
     # Create a simple drag coefficient model
     # Cd increases with Mach and alpha, slight dependency on Reynolds
-    M, Re, A = np.meshgrid(mach, reynolds, alpha, indexing='ij')
+    M, Re, A = np.meshgrid(mach, reynolds, alpha, indexing="ij")
     cd_data = 0.3 + 0.1 * M - 1e-7 * Re + 0.01 * A
-    
+
     cd_func = Function.from_grid(
         cd_data,
         [mach, reynolds, alpha],
-        inputs=['Mach', 'Reynolds', 'Alpha'],
-        outputs='Cd'
+        inputs=["Mach", "Reynolds", "Alpha"],
+        outputs="Cd",
     )
-    
+
     # Test at grid points
     assert abs(cd_func(0.0, 1e5, 0.0) - 0.29) < 0.01  # 0.3 - 1e-7*1e5
     assert abs(cd_func(1.0, 5e5, 0.0) - 0.35) < 0.01  # 0.3 + 0.1*1 - 1e-7*5e5
-    
+
     # Test interpolation between points
     result = cd_func(0.5, 3e5, 1.0)
     # Expected roughly: 0.3 + 0.1*0.5 - 1e-7*3e5 + 0.01*1.0 = 0.32
@@ -82,18 +72,14 @@ def test_from_grid_extrapolation_constant():
     """Test that constant extrapolation clamps to edge values."""
     x = np.array([0.0, 1.0, 2.0])
     y = np.array([0.0, 1.0, 4.0])  # y = x^2
-    
+
     func = Function.from_grid(
-        y,
-        [x],
-        inputs=['x'],
-        outputs='y',
-        extrapolation='constant'
+        y, [x], inputs=["x"], outputs="y", extrapolation="constant"
     )
-    
+
     # Test below lower bound - should return value at x=0
     assert func(-1.0) == 0.0
-    
+
     # Test above upper bound - should return value at x=2
     assert func(3.0) == 4.0
 
@@ -102,37 +88,39 @@ def test_from_grid_validation_errors():
     """Test that from_grid raises appropriate errors for invalid inputs."""
     x = np.array([0.0, 1.0, 2.0])
     y = np.array([0.0, 1.0, 2.0])
-    
+
     # Mismatched dimensions
-    X, Y = np.meshgrid(x, y, indexing='ij')
+    X, Y = np.meshgrid(x, y, indexing="ij")
     z_data = X + Y
-    
+
     # Wrong number of axes
     with pytest.raises(ValueError, match="Number of axes"):
-        Function.from_grid(z_data, [x], inputs=['x'], outputs='z')
-    
+        Function.from_grid(z_data, [x], inputs=["x"], outputs="z")
+
     # Wrong axis length
     with pytest.raises(ValueError, match="Axis 1 has"):
-        Function.from_grid(z_data, [x, np.array([0.0, 1.0])], inputs=['x', 'y'], outputs='z')
-    
+        Function.from_grid(
+            z_data, [x, np.array([0.0, 1.0])], inputs=["x", "y"], outputs="z"
+        )
+
     # Wrong number of inputs
     with pytest.raises(ValueError, match="Number of inputs"):
-        Function.from_grid(z_data, [x, y], inputs=['x'], outputs='z')
+        Function.from_grid(z_data, [x, y], inputs=["x"], outputs="z")
 
 
 def test_from_grid_default_inputs():
     """Test that from_grid uses default input names when not provided."""
     x = np.array([0.0, 1.0, 2.0])
     y = np.array([0.0, 1.0, 2.0])
-    
-    X, Y = np.meshgrid(x, y, indexing='ij')
+
+    X, Y = np.meshgrid(x, y, indexing="ij")
     z_data = X + Y
-    
+
     func = Function.from_grid(z_data, [x, y])
-    
+
     # Should use default names
-    assert 'x0' in func.__inputs__
-    assert 'x1' in func.__inputs__
+    assert "x0" in func.__inputs__
+    assert "x1" in func.__inputs__
 
 
 def test_from_grid_backward_compatibility():
@@ -140,12 +128,12 @@ def test_from_grid_backward_compatibility():
     # Test 1D function from list
     func1 = Function([[0, 0], [1, 1], [2, 4], [3, 9]])
     assert func1(1.5) > 0  # Should interpolate
-    
+
     # Test 2D function from array
     data = np.array([[0, 0, 0], [1, 0, 1], [0, 1, 2], [1, 1, 3]])
     func2 = Function(data)
     assert func2(0.5, 0.5) > 0  # Should interpolate
-    
+
     # Test callable function
     func3 = Function(lambda x: x**2)
     assert func3(2) == 4

From efe5cd72b35c4396e53bf85b6dba961c1d49b2ed Mon Sep 17 00:00:00 2001
From: Ishan <99824864+aZira371@users.noreply.github.com>
Date: Sun, 16 Nov 2025 00:46:59 +0530
Subject: [PATCH 04/11] MNt: refactoring get_drag_coefficient in flight.py

- MNT: velocity_body was not being used in get_drag_coefficient, removed it as an input
---
 rocketpy/simulation/flight.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py
index 83059f246..2735c57e0 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -1349,7 +1349,7 @@ def lateral_surface_wind(self):
         return -wind_u * np.cos(heading_rad) + wind_v * np.sin(heading_rad)
 
     def __get_drag_coefficient(
-        self, drag_function, mach, z, velocity_body, freestream_velocity_body
+        self, drag_function, mach, z, freestream_velocity_body
     ):
         """Calculate drag coefficient, handling both 1D and multi-dimensional functions.
 

From 17fc587477c058c37bdacd8666fa3cba9c693be0 Mon Sep 17 00:00:00 2001
From: Ishan <dubeyishan371@gmail.com>
Date: Thu, 20 Nov 2025 00:24:50 +0530
Subject: [PATCH 05/11] MNT: refactoring in flight.py and lint corrections to
 function.py and test_multidim_drag.py

-MNT:  removed unused velocity in body frame parameter requirement from all instances of get_drag_coefficient in flight.py

- MNT: corrected docstring for get_value_opt_grid in function.py

- MNT: shifted import of classes before the definition of functions in test_multidim_drag.py
---
 rocketpy/mathutils/function.py          |  2 ++
 rocketpy/simulation/flight.py           | 11 +----------
 tests/integration/test_multidim_drag.py |  6 +-----
 3 files changed, 4 insertions(+), 15 deletions(-)

diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py
index 245640551..a1fbc1595 100644
--- a/rocketpy/mathutils/function.py
+++ b/rocketpy/mathutils/function.py
@@ -674,6 +674,8 @@ def __get_value_opt_nd(self, *args):
 
     def __get_value_opt_grid(self, *args):
         """Evaluate the Function using RegularGridInterpolator for structured grids.
+        
+        This method is dynamically assigned in from_grid() class method.
 
         Parameters
         ----------
diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py
index 2735c57e0..47ffa9af4 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -1348,9 +1348,7 @@ def lateral_surface_wind(self):
 
         return -wind_u * np.cos(heading_rad) + wind_v * np.sin(heading_rad)
 
-    def __get_drag_coefficient(
-        self, drag_function, mach, z, freestream_velocity_body
-    ):
+    def __get_drag_coefficient(self, drag_function, mach, z, freestream_velocity_body):
         """Calculate drag coefficient, handling both 1D and multi-dimensional functions.
 
         Parameters
@@ -1361,8 +1359,6 @@ def __get_drag_coefficient(
             Mach number
         z : float
             Altitude in meters
-        velocity_body : Vector or array-like
-            Rocket velocity in body frame [vx_b, vy_b, vz_b]
         freestream_velocity_body : Vector or array-like
             Freestream velocity in body frame [stream_vx_b, stream_vy_b, stream_vz_b]
 
@@ -1486,7 +1482,6 @@ def udot_rail1(self, t, u, post_processing=False):
             self.rocket.power_on_drag,
             free_stream_mach,
             z,
-            [vx_b, vy_b, vz_b],
             [stream_vx_b, stream_vy_b, stream_vz_b],
         )
 
@@ -1679,7 +1674,6 @@ def u_dot(self, t, u, post_processing=False):  # pylint: disable=too-many-locals
                 self.rocket.power_on_drag,
                 free_stream_mach,
                 z,
-                [vx_b, vy_b, vz_b],
                 [stream_vx_b, stream_vy_b, stream_vz_b],
             )
         else:
@@ -1687,7 +1681,6 @@ def u_dot(self, t, u, post_processing=False):  # pylint: disable=too-many-locals
                 self.rocket.power_off_drag,
                 free_stream_mach,
                 z,
-                [vx_b, vy_b, vz_b],
                 [stream_vx_b, stream_vy_b, stream_vz_b],
             )
         rho = self.env.density.get_value_opt(z)
@@ -1976,7 +1969,6 @@ def u_dot_generalized(self, t, u, post_processing=False):  # pylint: disable=too
                 self.rocket.power_on_drag,
                 free_stream_mach,
                 z,
-                velocity_in_body_frame,
                 freestream_velocity_body,
             )
         else:
@@ -1985,7 +1977,6 @@ def u_dot_generalized(self, t, u, post_processing=False):  # pylint: disable=too
                 self.rocket.power_off_drag,
                 free_stream_mach,
                 z,
-                velocity_in_body_frame,
                 freestream_velocity_body,
             )
         R3 += -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff
diff --git a/tests/integration/test_multidim_drag.py b/tests/integration/test_multidim_drag.py
index accf20afa..642a2a263 100644
--- a/tests/integration/test_multidim_drag.py
+++ b/tests/integration/test_multidim_drag.py
@@ -1,14 +1,12 @@
 """Integration tests for multi-dimensional drag coefficient support."""
 
 import numpy as np
-import pytest
 
-from rocketpy import Function
+from rocketpy import Environment, Flight, Function, Rocket, SolidMotor
 
 
 def test_flight_with_1d_drag(example_plain_env, calisto):
     """Test that flights with 1D drag curves still work (backward compatibility)."""
-    from rocketpy import Flight
 
     flight = Flight(
         rocket=calisto,
@@ -26,8 +24,6 @@ def test_flight_with_1d_drag(example_plain_env, calisto):
 
 def test_flight_with_3d_drag_basic():
     """Test that a simple 3D drag function works."""
-    from rocketpy import Environment, Flight, Rocket, SolidMotor
-
     # Create environment
     env = Environment(gravity=9.81)
     env.set_atmospheric_model(type="standard_atmosphere")

From 0a416b568a831e748e8e11d9aeea4ab6fc9f1f52 Mon Sep 17 00:00:00 2001
From: Ishan <dubeyishan371@gmail.com>
Date: Tue, 25 Nov 2025 01:19:41 +0530
Subject: [PATCH 06/11] MNT: refactoring flight.py to remove unused parameters

---
 rocketpy/simulation/flight.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py
index 47ffa9af4..226f6cc66 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -1467,10 +1467,6 @@ def udot_rail1(self, t, u, post_processing=False):
         a32 = 2 * (e2 * e3 + e0 * e1)
         a33 = 1 - 2 * (e1**2 + e2**2)
 
-        vx_b = a11 * vx + a21 * vy + a31 * vz
-        vy_b = a12 * vx + a22 * vy + a32 * vz
-        vz_b = a13 * vx + a23 * vy + a33 * vz
-
         # Freestream velocity in body frame
         wind_vx = self.env.wind_velocity_x.get_value_opt(z)
         wind_vy = self.env.wind_velocity_y.get_value_opt(z)

From 9ee6a90370748e9df34e891a37e87ba6ca4065d5 Mon Sep 17 00:00:00 2001
From: Ishan <dubeyishan371@gmail.com>
Date: Tue, 25 Nov 2025 01:27:14 +0530
Subject: [PATCH 07/11] MNT: correction of docstring function.py

- MNT: rearranged the docstring of from_grid in function.py to match the expected output of doctest
---
 rocketpy/mathutils/function.py | 20 +++++++++++---------
 1 file changed, 11 insertions(+), 9 deletions(-)

diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py
index a1fbc1595..706d0237d 100644
--- a/rocketpy/mathutils/function.py
+++ b/rocketpy/mathutils/function.py
@@ -4093,7 +4093,16 @@ def from_grid(
         -------
         Function
             A Function object using RegularGridInterpolator for evaluation.
-
+        
+        Notes
+        -----
+        - Grid data must be on a regular (structured) grid.
+        - For unstructured data, use the regular Function constructor with
+          scattered points.
+        - Extrapolation with 'constant' mode uses the nearest edge values,
+          which is appropriate for aerodynamic coefficients where extrapolation
+          beyond the data range should be avoided.
+        
         Examples
         --------
         >>> import numpy as np
@@ -4113,15 +4122,8 @@ def from_grid(
         ... )
         >>> # Evaluate at a point
         >>> cd_func(1.2, 3e5, 3.0)
+        0.48000000000000004
 
-        Notes
-        -----
-        - Grid data must be on a regular (structured) grid.
-        - For unstructured data, use the regular Function constructor with
-          scattered points.
-        - Extrapolation with 'constant' mode uses the nearest edge values,
-          which is appropriate for aerodynamic coefficients where extrapolation
-          beyond the data range should be avoided.
         """
         # Validate inputs
         if not isinstance(grid_data, np.ndarray):

From 099155abf337a034625f5b0290ae3adb0b5f70c1 Mon Sep 17 00:00:00 2001
From: Ishan <dubeyishan371@gmail.com>
Date: Tue, 25 Nov 2025 01:30:47 +0530
Subject: [PATCH 08/11] MNT: make format and lint corrections to function.py

- MNT: reran make format and lint on function.py to correct after previous changes to from_grid
---
 rocketpy/mathutils/function.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py
index 706d0237d..e0635cdf0 100644
--- a/rocketpy/mathutils/function.py
+++ b/rocketpy/mathutils/function.py
@@ -674,7 +674,7 @@ def __get_value_opt_nd(self, *args):
 
     def __get_value_opt_grid(self, *args):
         """Evaluate the Function using RegularGridInterpolator for structured grids.
-        
+
         This method is dynamically assigned in from_grid() class method.
 
         Parameters
@@ -4093,7 +4093,7 @@ def from_grid(
         -------
         Function
             A Function object using RegularGridInterpolator for evaluation.
-        
+
         Notes
         -----
         - Grid data must be on a regular (structured) grid.
@@ -4102,7 +4102,7 @@ def from_grid(
         - Extrapolation with 'constant' mode uses the nearest edge values,
           which is appropriate for aerodynamic coefficients where extrapolation
           beyond the data range should be avoided.
-        
+
         Examples
         --------
         >>> import numpy as np

From 26e2e5590799e2738e6639ecfb0bf051adaaab36 Mon Sep 17 00:00:00 2001
From: Ishan <dubeyishan371@gmail.com>
Date: Tue, 25 Nov 2025 01:47:29 +0530
Subject: [PATCH 09/11] MNT: pylint adjustments for new methods in function.py

- MNT: disables pylint unused private member for get_value_opt_grid as it is called upon dynamically by from_grid

- MNT: disabled pylint too many statement for from_grid for now and added a to-do to refactor it into smaller methods/helper functions

- MNT: updated .pylintrc to record Re as good name
---
 .pylintrc                      | 1 +
 rocketpy/mathutils/function.py | 4 ++--
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/.pylintrc b/.pylintrc
index e417e0b11..991d9445f 100644
--- a/.pylintrc
+++ b/.pylintrc
@@ -229,6 +229,7 @@ good-names=FlightPhases,
            center_of_mass_without_motor_to_CDM,
            motor_center_of_dry_mass_to_CDM,
            generic_motor_cesaroni_M1520,
+           Re, # Reynolds number
 
 # Good variable names regexes, separated by a comma. If names match any regex,
 # they will always be accepted
diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py
index e0635cdf0..7f82540c3 100644
--- a/rocketpy/mathutils/function.py
+++ b/rocketpy/mathutils/function.py
@@ -672,7 +672,7 @@ def __get_value_opt_nd(self, *args):
 
         return result
 
-    def __get_value_opt_grid(self, *args):
+    def __get_value_opt_grid(self, *args): # pylint: disable=unused-private-member
         """Evaluate the Function using RegularGridInterpolator for structured grids.
 
         This method is dynamically assigned in from_grid() class method.
@@ -4057,7 +4057,7 @@ def from_grid(
         interpolation="linear_grid",
         extrapolation="constant",
         **kwargs,
-    ):
+    ): # pylint: disable=too-many-statements #TODO: Refactor this method into smaller methods
         """Creates a Function from N-dimensional grid data.
 
         This method is designed for structured grid data, such as CFD simulation

From bf6a18284fbec60d0dc6044b500e123658f72b03 Mon Sep 17 00:00:00 2001
From: Ishan <dubeyishan371@gmail.com>
Date: Tue, 25 Nov 2025 01:50:18 +0530
Subject: [PATCH 10/11] MNt: make format after previous change to function.py

---
 rocketpy/mathutils/function.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py
index 7f82540c3..c01d64c24 100644
--- a/rocketpy/mathutils/function.py
+++ b/rocketpy/mathutils/function.py
@@ -672,7 +672,7 @@ def __get_value_opt_nd(self, *args):
 
         return result
 
-    def __get_value_opt_grid(self, *args): # pylint: disable=unused-private-member
+    def __get_value_opt_grid(self, *args):  # pylint: disable=unused-private-member
         """Evaluate the Function using RegularGridInterpolator for structured grids.
 
         This method is dynamically assigned in from_grid() class method.
@@ -4057,7 +4057,7 @@ def from_grid(
         interpolation="linear_grid",
         extrapolation="constant",
         **kwargs,
-    ): # pylint: disable=too-many-statements #TODO: Refactor this method into smaller methods
+    ):  # pylint: disable=too-many-statements #TODO: Refactor this method into smaller methods
         """Creates a Function from N-dimensional grid data.
 
         This method is designed for structured grid data, such as CFD simulation

From 69caf403c9c10d1ec4eea0ff4e885a1f858c7c1e Mon Sep 17 00:00:00 2001
From: Ishan <dubeyishan371@gmail.com>
Date: Tue, 25 Nov 2025 02:08:50 +0530
Subject: [PATCH 11/11] MNT: removed Re where unused in test_multidim_drag.py

- MNT: Re variable was unused in test_3d_drag_with_varying_alpha thus replaced it
---
 tests/integration/test_multidim_drag.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/integration/test_multidim_drag.py b/tests/integration/test_multidim_drag.py
index 642a2a263..be5c79ee5 100644
--- a/tests/integration/test_multidim_drag.py
+++ b/tests/integration/test_multidim_drag.py
@@ -104,7 +104,7 @@ def test_3d_drag_with_varying_alpha():
     reynolds = np.array([1e5, 1e6])
     alpha = np.array([0.0, 5.0, 10.0, 15.0])
 
-    M, Re, A = np.meshgrid(mach, reynolds, alpha, indexing="ij")
+    M, _, A = np.meshgrid(mach, reynolds, alpha, indexing="ij")
     # Strong alpha dependency: Cd increases significantly with alpha
     cd_data = 0.3 + 0.05 * M + 0.03 * A
     cd_data = np.clip(cd_data, 0.2, 2.0)