ENH: Add multi-dimensional drag coefficient support (Cd as function of M, Re, α)

[Copilot]

Pull request type

• Code changes (bugfix, features)
• Code maintenance (refactoring, formatting, tests)

Checklist

• Tests for the changes have been added (if needed)
• Docs have been reviewed and added / updated
• Lint (black rocketpy/ tests/) has passed locally
• All tests (pytest tests -m slow --runslow) have passed locally
• CHANGELOG.md has been updated (if relevant)

Current behavior

Drag coefficients are limited to 1D functions of Mach number. CFD simulations typically produce drag data varying with Reynolds number and angle of attack, which cannot be properly represented.

New behavior

Drag coefficients now support N-dimensional structured grids via Function.from_grid(). Flight simulations automatically calculate and pass Reynolds number (from velocity, density, viscosity, rocket diameter) and angle of attack (from freestream velocity vectors in body frame) to multi-dimensional drag functions.

Implementation

Function Class (rocketpy/mathutils/function.py)

• Added from_grid() class method using scipy.interpolate.RegularGridInterpolator
• Added 'linear_grid' interpolation type with edge clamping (no extrapolation beyond bounds)
• Grid data stored as flattened array for serialization compatibility

Flight Class (rocketpy/simulation/flight.py)

• Added __get_drag_coefficient() helper detecting function dimensionality
• Calculates Re = ρVD/μ and α from body-frame velocities during integration
• Updated udot_rail1(), u_dot(), u_dot_generalized() to use helper

Rocket Class (rocketpy/rocket/rocket.py)

• Preserves Function objects when passed as drag coefficients (previously wrapped in 1D Function)

Usage Example

import numpy as np
from rocketpy import Function, Rocket

# Create 3D drag data from CFD
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])  # degrees

M, Re, A = np.meshgrid(mach, reynolds, alpha, indexing='ij')
cd_data = 0.3 + 0.1*M - 1e-7*Re + 0.01*A

# Create structured grid function
drag_func = Function.from_grid(
    cd_data, 
    [mach, reynolds, alpha],
    inputs=['Mach', 'Reynolds', 'Alpha'],
    outputs='Cd'
)

# Use directly in Rocket - Flight automatically handles M, Re, α
rocket = Rocket(..., power_off_drag=drag_func, power_on_drag=drag_func)

Testing

• 7 unit tests for Function.from_grid() covering 1D/2D/3D grids, extrapolation, validation
• 3 integration tests for 3D drag in flight simulations
• All 593 existing tests pass (full backward compatibility)

Breaking change

• No

Existing 1D drag curves (Mach only) work unchanged. Multi-dimensional support is opt-in via from_grid().

Additional information

Resolves #874

CodeQL security scan: 0 alerts

Original prompt

---

This section details on the original issue you should resolve

<issue_title>ENH: Drag curve as a function of Re and alpha as well</issue_title>
<issue_description>Is your feature request related to a problem? Please describe.

As of now drag curve is a function of speed (or mach number). Since a lot of rocketeers use CFD simulations, it is expected that drag values won't be just a function of mach number. Drag of any body moving in fluid can depend on many major parameters. Some of them are listed below:

• Reynolds Number (Re)
• Angle of Attack (alpha)
• Altitude/Air Properties
• Rocket Geometry/Shape
• Surface Roughness

Re accounts for the viscous effects of the flow (skin friction vs. pressure drag). Reynolds number depends on velocity, characteristic length, air viscosity, and density, and is crucial for modeling drag at low speeds and during laminar-turbulent transitions. Non-zero angles of attack change the effective cross-sectional area and cause flow separation or vortex formation, increasing drag. Altitude or atmospheric variation can be handled by RocketPy separately to update Re or atmospheric data rather than being built directly into the drag curve.
Drag could be input as Cd(M, α) or Cd(M, Re, α).

Describe the solution you'd like

1. A multidimensional array implementation of Cd=f(M,Re,α) in Rocket Class.
2. Extension of drag_coeff in flight class to allow multiple params as input - M, Re and alpha. Probably implement new get_value_opt like function specifically for this.
3. Check if current interpolation and extrapolation methods of the Function class are applicable in case of 2D/3D arrays. If not then update the same.
4. Ensure backwards compatibility to still allow for use of Cd=f(M).

Additional context

It does not make sense that the drag curve will be given as an input of Cd=f(Re) or Cd=f(α). Physically these are not capturing the entire flow physics at the same accuracy at which Cd=f(M) does as mach number based effects like shock wave propagation become more important in rocket flights. May be still allow for single array functionality with Re and α but provide a note that the results may not be accurate!</issue_description>

<agent_instructions>Here is a specific, technical prompt you can use to delegate this issue to the GitHub Copilot agent. This prompt directs it to the correct implementation strategy we discussed (enhancing the Function class).

---

Optional Prompt for GitHub Copilot

Here is a specific prompt to guide the implementation for Issue #874.

1. Core Objective:
Implement the enhancement described in issue #874. The goal is to allow the drag coefficient Cd to be defined as a 3D function of Mach (M), Reynolds number (Re), and Angle of Attack ($\alpha$).

2. Implementation Location:
The primary changes must be made within the rocketpy.Function class (in rocketpy/function.py). This is not just a change to the Flight or Rocket class; the Function class itself must be enhanced to handle N-dimensional structured grid data.

3. Add New Interpolation Method:
The current shepard interpolation is inefficient for structured CFD grids.

• Add a new interpolation option (e.g., interpolation='linear_grid' or interpolation='grid') to the Function class __init__.
• This new method MUST use scipy.interpolate.RegularGridInterpolator as its backend. You will need to import this from scipy.interpolate.
• The Function class state should store the initialized RegularGridInterpolator object to be called later.

4. Add New Function.from_grid() Class Method:
Manually creating a flattened list of points is impractical for users.

• Create a new class method: @classmethod def from_grid(cls, grid_data, axes, inputs, output, interpolation='linear_grid', **kwargs).
• This method will accept:
  • grid_data: An N-dimensional NumPy array (e.g., a 3D array for $C_d$).
  • axes: A list of 1D arrays representing the axes of the grid (e.g., [M_axis, Re_axis, A_axis]).
  • inputs: A list of strings, e.g., ['Mach', 'Reynolds', 'Alpha'].
  • output: A string, e.g., 'Cd'.
• Inside this method, you will initialize the scipy.interpolate.RegularGridInterpolator with the provided axes and grid_data and store it in the new Function object.

5. Handle Extrapolation (Crucial):
The new grid interpolation must clamp to the edge (i.e., use the nearest edge value) for out-of-bounds requests. It should not extrapolate.

• When initializing RegularGridInterpolator, set bounds_error=False.
• The fill_value should be set to None, which (when bounds_error=False) causes the interpolator to use the edge values for out-of-bounds points.

6. Update Flight Class:
The Flight class (likely in flight.py or related files) must be modified to use this new function.

• The flight simulation loop needs to detect the type of drag function it has.
• If the `rocket.aerodynamic_surfaces.dra...

• Fixes ENH: Drag curve as a function of Re and alpha as well #874

---

💬 We'd love your input! Share your thoughts on Copilot coding agent in our 2 minute survey.

[aZira371]

The code implementation is as I planned it to be done. There might be a slight edge case which should be addressed in regular grid interpolator unit tests or with warning: that the grid was not sorted in ascending order as intended. This can happen with repeated values in the grid. Apart from that the changes are okay to be merged! @Gui-FernandesBR @copilot

[aZira371]

@Gui-FernandesBR Is there a way to see this branch on VSC? I am unable to access it there. It is not being listed as a remote branch. I think one of the previous lint errors was related to the new grid method being too big in function.py. It should be ideally broken down into helper functions or pylint has to be disabled for it.

[codecov]

Codecov Report

❌ Patch coverage is 89.03226% with 17 lines in your changes missing coverage. Please review.
✅ Project coverage is 80.36%. Comparing base (9cf3dd4) to head (74fe825).
⚠️ Report is 6 commits behind head on develop.

Files with missing lines	Patch %	Lines
rocketpy/mathutils/function.py	82.79%	16 Missing ⚠️
rocketpy/simulation/flight.py	98.21%	1 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           develop     #875      +/-   ##
===========================================
+ Coverage    80.27%   80.36%   +0.08%     
===========================================
  Files          104      106       +2     
  Lines        12769    13122     +353     
===========================================
+ Hits         10250    10545     +295     
- Misses        2519     2577      +58     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[aZira371]

The pull request is ready for review now. Pylint is failing due to a problem in .pylintrc, I am not sure if this file should be corrected in this PR or not as this is an unrelated error. Tests did pass locally on my end before merging from develop. @RocketPy-Team/code-owners will need help from you, because different tests are failing at different places - on github test_sensor failed but on local after merging develop test_environment failed.

[Gui-FernandesBR]

@aZira371 is linear_grid really the best name here? regular_grid should be more similar to scipy, no?

[Copilot]

Pull request overview

This PR implements multi-dimensional drag coefficient support, allowing drag to vary with Mach number, Reynolds number, and angle of attack. This enhancement enables RocketPy users to utilize CFD simulation data that naturally depends on multiple flow parameters, significantly improving aerodynamic modeling accuracy.

Key changes:

• Added Function.from_grid() class method for creating N-dimensional structured grid interpolators using scipy.interpolate.RegularGridInterpolator
• Modified Flight class to automatically calculate Reynolds number and angle of attack during trajectory integration, passing these to multi-dimensional drag functions
• Updated Rocket class to preserve Function objects when provided directly as drag coefficients (backward compatible with existing 1D usage)

Reviewed changes

Copilot reviewed 6 out of 6 changed files in this pull request and generated 7 comments.

Show a summary per file

File	Description
rocketpy/mathutils/function.py	Adds from_grid() class method, linear_grid interpolation type, and __get_value_opt_grid() for structured grid evaluation with edge-clamping extrapolation
rocketpy/simulation/flight.py	Implements __get_drag_coefficient() helper to detect function dimensionality and calculate Re/α during integration; updates all u_dot methods
rocketpy/rocket/rocket.py	Preserves Function objects when passed as drag coefficients instead of wrapping them, enabling multi-dimensional drag
tests/unit/mathutils/test_function_grid.py	Comprehensive unit tests for grid interpolation (1D/2D/3D), extrapolation modes, validation, and fallback behavior
tests/integration/test_multidim_drag.py	Integration tests verifying 3D drag works in full flight simulations and maintains backward compatibility with 1D drag
.pylintrc	Adds 'Re' (Reynolds number) to accepted variable names

[Copilot]

The docstring for test_shepard_fallback_warning states that "the Function class should fall back to shepard interpolation and emit a warning," but the test implementation doesn't actually test the fallback scenario correctly. The test manually changes interpolation to 'linear_grid' on a Function that wasn't created with from_grid(), which is not a realistic scenario.

Consider renaming the test to better reflect what it actually tests (e.g., test_linear_grid_without_grid_interpolator_warns) or restructure it to test a more realistic fallback scenario if one exists.

[Copilot]

Potential division by zero: reynolds = rho * freestream_speed * characteristic_length / mu on line 1412 doesn't check if mu (dynamic viscosity) is zero or near-zero before division. While atmospheric models typically provide non-zero viscosity values, defensive programming suggests adding a check or using a safe division operation to prevent crashes in edge cases.

Consider: reynolds = rho * freestream_speed * characteristic_length / max(mu, 1e-10) or adding an explicit check.

Suggested change

	reynolds = rho * freestream_speed * characteristic_length / mu
	# Defensive: avoid division by zero in case mu is zero or near-zero
	reynolds = rho * freestream_speed * characteristic_length / max(mu, 1e-10)

[Copilot]

The test test_3d_drag_with_varying_alpha verifies that the drag function responds to different alpha values when called directly, but doesn't actually verify that the flight simulation is using angle of attack correctly during integration. Consider adding assertions that check the drag coefficient values during an actual flight simulation (e.g., by examining flight data or verifying that a rocket with high-alpha-sensitivity drag behaves differently than one without).

[Copilot]

The docstring states that extrapolation can be 'constant' or 'zero', but the __get_value_opt_grid method (line 726-727) allows other extrapolation types (like 'natural') to fall through to the interpolator's default behavior. The docstring should document what happens with 'natural' extrapolation (linear extrapolation from edge gradients) or explicitly state that only 'constant' and 'zero' are supported for grid interpolation.

Additionally, consider adding validation in from_grid() to warn or error if an unsupported extrapolation type is specified.

[Copilot]

Potential issue with y_array initialization for multi-dimensional grids. Line 4207 sets func.y_array = func._image[: len(axes[0])] which only takes the first n values from the flattened grid where n is the length of the first axis. This creates an inconsistent y_array that doesn't properly represent the function data.

For grid functions with dimensions > 1, y_array is not a meaningful concept (it's designed for 1D functions). Consider either:

1. Setting y_array to the full flattened image: func.y_array = func._image
2. Setting it to None or an empty array with a comment explaining that y_array is not applicable for grid functions
3. Documenting this limitation clearly

This could cause issues if other parts of the code expect y_array to contain all function values.

Suggested change

	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:
	# For 1D functions, y_array contains the function values.

[Copilot]

Performance consideration: Creating the full flattened domain and image arrays (lines 4185-4191) for compatibility purposes requires O(n^d) memory where n is the typical axis length and d is the number of dimensions. For large 3D grids (e.g., 100x100x100 points), this creates redundant data since the grid interpolator already has the structured data.

While this is necessary for backward compatibility with code expecting _domain and _image attributes, consider documenting this memory overhead and potentially adding a parameter to skip this flattening if compatibility with existing Function methods is not required.

[Copilot]

[nitpick] The input parameter name matching uses lowercase comparison and substring matching (e.g., "mach" in name, "reynolds" in name). While functional, this could match unintended parameter names (e.g., "machining" would match "mach").

Consider using more precise matching such as:

• Exact matches: name in ['mach', 'm']
• Or document the substring matching behavior clearly in the method docstring

This is a minor robustness improvement for edge cases.

[Gui-FernandesBR]

I'm not confortable with this comparision here.

        if hasattr(drag_function, "__inputs__") and len(drag_function.__inputs__) > 1:

This could be done within the Function class.

An attribute called is_multidimensioanl: bool should be really useful.

Then the comparision here would be:

if isinstance(foo, Function) and foo.is_multidimensional:

muche easier to read!

[Gui-FernandesBR]

• @MateusStano please review changes in the flight.py file
• @phmbressan please review changes in the function.py ile

Everything is looking good, @aZira371 great work here! I'm confortable with what has been done so far. I believe we can include this one on the next release.

---

@aZira371 could you please update changelog?

[aZira371]

Thanks for the review @Gui-FernandesBR will soon work on this and update the changelog.