diff --git a/.github/copilot-instructions.md b/.github/copilot-instructions.md new file mode 100644 index 000000000..f5366cb3b --- /dev/null +++ b/.github/copilot-instructions.md @@ -0,0 +1,221 @@ +# GitHub Copilot Instructions for RocketPy + +This file provides instructions for GitHub Copilot when working on the RocketPy codebase. +These guidelines help ensure consistency with the project's coding standards and development practices. + +## Project Overview + +RocketPy is a Python library for 6-DOF rocket trajectory simulation. +It's designed for high-power rocketry applications with focus on accuracy, performance, and ease of use. + +## Coding Standards + +### Naming Conventions +- **Use `snake_case` for all new code** - variables, functions, methods, and modules +- **Use descriptive names** - prefer `angle_of_attack` over `a` or `alpha` +- **Class names use PascalCase** - e.g., `SolidMotor`, `Environment`, `Flight` +- **Constants use UPPER_SNAKE_CASE** - e.g., `DEFAULT_GRAVITY`, `EARTH_RADIUS` + +### Code Style +- Follow **PEP 8** guidelines +- Line length: **88 characters** (Black's default) +- Organize imports with **isort** +- Our official formatter is the **ruff frmat** + +### Documentation +- **All public classes, methods, and functions must have docstrings** +- Use **NumPy style docstrings** +- Include **Parameters**, **Returns**, and **Examples** sections +- Document **units** for physical quantities (e.g., "in meters", "in radians") + +### Testing +- Write **unit tests** for all new features using pytest +- Follow **AAA pattern** (Arrange, Act, Assert) +- Use descriptive test names following: `test_methodname_expectedbehaviour` +- Include test docstrings explaining expected behavior +- Use **parameterization** for testing multiple scenarios +- Create pytest fixtures to avoid code repetition + +## Domain-Specific Guidelines + +### Physical Units and Conventions +- **SI units by default** - meters, kilograms, seconds, radians +- **Document coordinate systems** clearly (e.g., "tail_to_nose", "nozzle_to_combustion_chamber") +- **Position parameters** are critical - always document reference points +- Use **descriptive variable names** for physical quantities + +### Rocket Components +- **Motors**: SolidMotor, HybridMotor and LiquidMotor classes are children classes of the Motor class +- **Aerodynamic Surfaces**: They have Drag curves and lift coefficients +- **Parachutes**: Trigger functions, deployment conditions +- **Environment**: Atmospheric models, weather data, wind profiles + +### Mathematical Operations +- Use **numpy arrays** for vectorized operations (this improves performance) +- Prefer **scipy functions** for numerical integration and optimization +- **Handle edge cases** in calculations (division by zero, sqrt of negative numbers) +- **Validate input ranges** for physical parameters +- Monte Carlo simulations: sample from `numpy.random` for random number generation and creates several iterations to assess uncertainty in simulations. + +## File Structure and Organization + +### Source Code Organization + +Reminds that `rocketpy` is a Python package served as a library, and its source code is organized into several modules to facilitate maintainability and clarity. The following structure is recommended: + +``` +rocketpy/ +├── core/ # Core simulation classes +├── motors/ # Motor implementations +├── environment/ # Atmospheric and environmental models +├── plots/ # Plotting and visualization +├── tools/ # Utility functions +└── mathutils/ # Mathematical utilities +``` + +Please refer to popular Python packages like `scipy`, `numpy`, and `matplotlib` for inspiration on module organization. + +### Test Organization +``` +tests/ +├── unit/ # Unit tests +├── integration/ # Integration tests +├── acceptance/ # Acceptance tests +└── fixtures/ # Test fixtures organized by component +``` + +### Documentation Structure +``` +docs/ +├── user/ # User guides and tutorials +├── development/ # Development documentation +├── reference/ # API reference +├── examples/ # Flight examples and notebooks +└── technical/ # Technical documentation +``` + +## Common Patterns and Practices + +### Error Handling +- Use **descriptive error messages** with context +- **Validate inputs** at class initialization and method entry +- Raise **appropriate exception types** (ValueError, TypeError, etc.) +- Include **suggestions for fixes** in error messages + +### Performance Considerations +- Use **vectorized operations** where possible +- **Cache expensive computations** when appropriate (we frequently use `cached_property`) +- Keep in mind that RocketPy must be fast! + +### Backward Compatibility +- **Avoid breaking changes** in public APIs +- Use **deprecation warnings** before removing features +- **Document code changes** in docstrings and CHANGELOG + +## AI Assistant Guidelines + +### Code Generation +- **Always include docstrings** for new functions and classes +- **Follow existing patterns** in the codebase +- **Consider edge cases** and error conditions + +### Code Review and Suggestions +- **Check for consistency** with existing code style +- **Verify physical units** and coordinate systems +- **Ensure proper error handling** and input validation +- **Suggest performance improvements** when applicable +- **Recommend additional tests** for new functionality + +### Documentation Assistance +- **Use NumPy docstring format** consistently +- **Include practical examples** in docstrings +- **Document physical meanings** of parameters +- **Cross-reference related functions** and classes + +## Testing Guidelines + +### Unit Tests +- **Test individual methods** in isolation +- **Use fixtures** from the appropriate test fixture modules +- **Mock external dependencies** when necessary +- **Test both happy path and error conditions** + +### Integration Tests +- **Test interactions** between components +- **Verify end-to-end workflows** (Environment → Motor → Rocket → Flight) + +### Test Data +- **Use realistic parameters** for rocket simulations +- **Include edge cases** (very small/large rockets, extreme conditions) +- **Test with different coordinate systems** and orientations + +## Project-Specific Considerations + +### User Experience +- **Provide helpful error messages** with context and suggestions +- **Include examples** in docstrings and documentation +- **Support common use cases** with reasonable defaults + +## Examples of Good Practices + +### Function Definition +```python +def calculate_drag_force( + velocity, + air_density, + drag_coefficient, + reference_area +): + """Calculate drag force using the standard drag equation. + + Parameters + ---------- + velocity : float + Velocity magnitude in m/s. + air_density : float + Air density in kg/m³. + drag_coefficient : float + Dimensionless drag coefficient. + reference_area : float + Reference area in m². + + Returns + ------- + float + Drag force in N. + + Examples + -------- + >>> drag_force = calculate_drag_force(100, 1.225, 0.5, 0.01) + >>> print(f"Drag force: {drag_force:.2f} N") + """ + if velocity < 0: + raise ValueError("Velocity must be non-negative") + if air_density <= 0: + raise ValueError("Air density must be positive") + if reference_area <= 0: + raise ValueError("Reference area must be positive") + + return 0.5 * air_density * velocity**2 * drag_coefficient * reference_area +``` + +### Test Example +```python +def test_calculate_drag_force_returns_correct_value(): + """Test drag force calculation with known inputs.""" + # Arrange + velocity = 100.0 # m/s + air_density = 1.225 # kg/m³ + drag_coefficient = 0.5 + reference_area = 0.01 # m² + expected_force = 30.625 # N + + # Act + result = calculate_drag_force(velocity, air_density, drag_coefficient, reference_area) + + # Assert + assert abs(result - expected_force) < 1e-6 +``` + + +Remember: RocketPy prioritizes accuracy, performance, and usability. Always consider the physical meaning of calculations and provide clear, well-documented interfaces for users. diff --git a/.github/workflows/linters.yml b/.github/workflows/linters.yml index d998b5f38..3e5a3e81d 100644 --- a/.github/workflows/linters.yml +++ b/.github/workflows/linters.yml @@ -16,7 +16,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: ["3.9"] + python-version: ["3.10"] steps: - uses: actions/checkout@main - name: Set up Python ${{ matrix.python-version }} diff --git a/.github/workflows/publish-to-pypi.yml b/.github/workflows/publish-to-pypi.yml index 8724c2a94..ba2de7648 100644 --- a/.github/workflows/publish-to-pypi.yml +++ b/.github/workflows/publish-to-pypi.yml @@ -19,7 +19,7 @@ jobs: - name: Set up Python uses: actions/setup-python@main with: - python-version: "3.9" + python-version: "3.10" - name: Install dependencies run: | python -m pip install --upgrade pip diff --git a/.github/workflows/test-pytest-slow.yaml b/.github/workflows/test-pytest-slow.yaml index 76e6e7e28..fb58f62de 100644 --- a/.github/workflows/test-pytest-slow.yaml +++ b/.github/workflows/test-pytest-slow.yaml @@ -5,7 +5,7 @@ on: - cron: "0 17 * * 5" # at 05:00 PM, only on Friday push: branches: - - main + - master paths: - "**.py" - ".github/**" @@ -21,9 +21,10 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.9, 3.13] + python-version: ["3.10", "3.14"] env: PYTHON: ${{ matrix.python-version }} + MPLBACKEND: Agg steps: - uses: actions/checkout@main - name: Set up Python diff --git a/.github/workflows/test_pytest.yaml b/.github/workflows/test_pytest.yaml index 97fb068cc..f2a13fb13 100644 --- a/.github/workflows/test_pytest.yaml +++ b/.github/workflows/test_pytest.yaml @@ -19,24 +19,22 @@ jobs: strategy: matrix: os: [ubuntu-latest, macos-latest, windows-latest] - python-version: [3.9, 3.13] + python-version: ["3.10", "3.14"] env: OS: ${{ matrix.os }} PYTHON: ${{ matrix.python-version }} + MPLBACKEND: Agg steps: - uses: actions/checkout@main - name: Set up Python uses: actions/setup-python@main with: python-version: ${{ matrix.python-version }} - - - name: Cache Python dependencies - uses: actions/cache@main - with: - path: ~/.cache/pip - key: ${{ runner.os }}-pip-${{ hashFiles('**/requirements-tests.txt') }} - restore-keys: | - ${{ runner.os }}-pip- + cache: 'pip' + cache-dependency-path: | + requirements.txt + requirements-tests.txt + requirements-optional.txt - name: Install rocketpy run: pip install . diff --git a/.pylintrc b/.pylintrc index 4ae88a7c6..e417e0b11 100644 --- a/.pylintrc +++ b/.pylintrc @@ -88,7 +88,7 @@ persistent=yes # Minimum Python version to use for version dependent checks. Will default to # the version used to run pylint. -py-version=3.9 +py-version=3.10 # Discover python modules and packages in the file system subtree. recursive=no diff --git a/.vscode/settings.json b/.vscode/settings.json index 7fa0c1ff8..6399d8142 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -119,6 +119,7 @@ "filt", "firstsimulation", "flightcsys", + "flightusage", "Fluxogram", "fmax", "fmin", @@ -181,6 +182,7 @@ "Kaleb", "Karman", "Krasser", + "Kutta", "labelrotation", "linalg", "Lince", @@ -273,6 +275,7 @@ "rtol", "rtype", "rucsoundings", + "Runge", "runslow", "rwork", "savetxt", diff --git a/CHANGELOG.md b/CHANGELOG.md index 28586a1b7..144a2f1a1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,19 +32,38 @@ Attention: The newest changes should be on top --> ### Added +- ENH: Tank Fluids with Variable Density from Temperature and Pressure [#852](https://github.com/RocketPy-Team/RocketPy/pull/852) +- ENH: Controller (AirBrakes) and Sensors Encoding [#849](https://github.com/RocketPy-Team/RocketPy/pull/849) +- EHN: Addition of ensemble variable to ECMWF dictionaries [#842](https://github.com/RocketPy-Team/RocketPy/pull/842) +- ENH: Added Crop and Clip Methods to Function Class [#817](https://github.com/RocketPy-Team/RocketPy/pull/817) +- DOC: Add Flight class usage documentation and update index [#841](https://github.com/RocketPy-Team/RocketPy/pull/841) +- ENH: Discretized and No-Pickle Encoding Options [#827](https://github.com/RocketPy-Team/RocketPy/pull/827) +- ENH: Add the Coriolis Force to the Flight class [#799](https://github.com/RocketPy-Team/RocketPy/pull/799) +- ENH: Improve parachute geometric parametrization [#835](https://github.com/RocketPy-Team/RocketPy/pull/835) +- ENH: Changing ellipses plot axis label [#855](https://github.com/RocketPy-Team/RocketPy/pull/855) ### Changed +- MNT: allow for exporting of non apogee flights. [#863](https://github.com/RocketPy-Team/RocketPy/pull/863) +- TST: remove remaining files after test session. [#862](https://github.com/RocketPy-Team/RocketPy/pull/862) +- MNT: bumps min python version to 3.10 [#857](https://github.com/RocketPy-Team/RocketPy/pull/857) +- DOC: Update docs dependencies and sub dependencies [#851](https://github.com/RocketPy-Team/RocketPy/pull/851) +- MNT: extract flight data exporters [#845](https://github.com/RocketPy-Team/RocketPy/pull/845) +- ENH: _MotorPrints inheritance - issue #460 [#828](https://github.com/RocketPy-Team/RocketPy/pull/828) +- MNT: fix deprecations and warnings [#829](https://github.com/RocketPy-Team/RocketPy/pull/829) ### Fixed +- BUG: correct encoding for trapezoidal sweep length and angle. [#861](https://github.com/RocketPy-Team/RocketPy/pull/861) +- BUG: Fix no time initialization when passing initial_solution as array to Flight object [#844](https://github.com/RocketPy-Team/RocketPy/pull/844) + ## [v1.10.0] - 2025-05-16 ### Added - ENH: Support for ND arithmetic in Function class. [#810] (https://github.com/RocketPy-Team/RocketPy/pull/810) - ENH: allow users to provide custom samplers [#803](https://github.com/RocketPy-Team/RocketPy/pull/803) -- ENH: Implement Multivariate Rejection Sampling (MRS) [#738] (https://github.com/RocketPy-Team/RocketPy/pull/738) +- ENH: Implement Multivariate Rejection Sampling (MRS) [#738] (https://github.com/RocketPy-Team/RocketPy/pull/738) - ENH: Create a rocketpy file to store flight simulations [#800](https://github.com/RocketPy-Team/RocketPy/pull/800) - ENH: Support for the RSE file format has been added to the library [#798](https://github.com/RocketPy-Team/RocketPy/pull/798) - ENH: Introduce Net Thrust with pressure corrections [#789](https://github.com/RocketPy-Team/RocketPy/pull/789) @@ -303,6 +322,7 @@ You can install this version by running `pip install rocketpy==1.3.0` ### Fixed +- BUG: Fixes StochasticNoseCone powerseries issue #838 [#839](https://github.com/RocketPy-Team/RocketPy/pull/839) - MNT: Alter PYPI classifier naming. [#615](https://github.com/RocketPy-Team/RocketPy/pull/615) - DOC: Solve Dependencies Conflicts and pyproject build [#613](https://github.com/RocketPy-Team/RocketPy/pull/613) - BUG: Fixes nose cone bluffness issue #610 [#611](https://github.com/RocketPy-Team/RocketPy/pull/611) diff --git a/README.md b/README.md index 3eb63c9a6..8247ab1a3 100644 --- a/README.md +++ b/README.md @@ -262,6 +262,9 @@ main = calisto.add_parachute( sampling_rate=105, lag=1.5, noise=(0, 8.3, 0.5), + radius=1.5, + height=1.5, + porosity=0.0432, ) drogue = calisto.add_parachute( @@ -271,6 +274,9 @@ drogue = calisto.add_parachute( sampling_rate=105, lag=1.5, noise=(0, 8.3, 0.5), + radius=1.5, + height=1.5, + porosity=0.0432, ) ``` diff --git a/data/motors/liquid_motor_example/n2o_fluid_parameters.csv b/data/motors/liquid_motor_example/n2o_fluid_parameters.csv new file mode 100644 index 000000000..b6bfb2ee9 --- /dev/null +++ b/data/motors/liquid_motor_example/n2o_fluid_parameters.csv @@ -0,0 +1,101 @@ +Temperature (K),Pressure (Pa),Density (kg/m³) +193.150000,101325.000000,2.839688 +205.372222,101325.000000,2.659196 +217.594444,101325.000000,2.501930 +229.816667,101325.000000,2.363172 +242.038889,101325.000000,2.239591 +254.261111,101325.000000,2.128693 +266.483333,101325.000000,2.028539 +278.705556,101325.000000,1.937588 +290.927778,101325.000000,1.854594 +303.150000,101325.000000,1.778532 +193.150000,990733.333333,1206.532256 +205.372222,990733.333333,1168.931406 +217.594444,990733.333333,1129.335464 +229.816667,990733.333333,1087.042281 +242.038889,990733.333333,24.480299 +254.261111,990733.333333,22.807515 +266.483333,990733.333333,21.414919 +278.705556,990733.333333,20.221402 +290.927778,990733.333333,19.178742 +303.150000,990733.333333,18.255258 +193.150000,1880141.666667,1207.901606 +205.372222,1880141.666667,1170.560205 +217.594444,1880141.666667,1131.305292 +229.816667,1880141.666667,1089.479276 +242.038889,1880141.666667,1044.107836 +254.261111,1880141.666667,993.574261 +266.483333,1880141.666667,44.818526 +278.705556,1880141.666667,41.521255 +290.927778,1880141.666667,38.847302 +303.150000,1880141.666667,36.598433 +193.150000,2769550.000000,1209.256890 +205.372222,2769550.000000,1172.168610 +217.594444,2769550.000000,1133.244437 +229.816667,2769550.000000,1091.867709 +242.038889,2769550.000000,1047.140735 +254.261111,2769550.000000,997.611513 +266.483333,2769550.000000,940.542069 +278.705556,2769550.000000,67.699444 +290.927778,2769550.000000,61.910544 +303.150000,2769550.000000,57.447399 +193.150000,3658958.333333,1210.598454 +205.372222,3658958.333333,1173.757228 +217.594444,3658958.333333,1135.154031 +229.816667,3658958.333333,1094.209859 +242.038889,3658958.333333,1050.096018 +254.261111,3658958.333333,1001.503955 +266.483333,3658958.333333,946.059899 +278.705556,3658958.333333,878.329006 +290.927778,3658958.333333,90.504952 +303.150000,3658958.333333,81.893190 +193.150000,4548366.666667,1211.926631 +205.372222,4548366.666667,1175.326641 +217.594444,4548366.666667,1137.035143 +229.816667,4548366.666667,1096.507838 +242.038889,4548366.666667,1052.978332 +254.261111,4548366.666667,1005.263627 +266.483333,4548366.666667,951.292446 +278.705556,4548366.666667,886.629576 +290.927778,4548366.666667,130.554270 +303.150000,4548366.666667,112.079628 +193.150000,5437775.000000,1213.241741 +205.372222,5437775.000000,1176.877402 +217.594444,5437775.000000,1138.888779 +229.816667,5437775.000000,1098.763608 +242.038889,5437775.000000,1055.791888 +254.261111,5437775.000000,1008.901026 +266.483333,5437775.000000,956.273388 +278.705556,5437775.000000,894.237174 +290.927778,5437775.000000,812.409895 +303.150000,5437775.000000,153.463224 +193.150000,6327183.333333,1214.544092 +205.372222,6327183.333333,1178.410038 +217.594444,6327183.333333,1140.715888 +229.816667,6327183.333333,1100.978992 +242.038889,6327183.333333,1058.540519 +254.261111,6327183.333333,1012.425369 +266.483333,6327183.333333,961.030387 +278.705556,6327183.333333,901.279439 +290.927778,6327183.333333,825.421408 +303.150000,6327183.333333,689.418589 +193.150000,7216591.666667,1215.833981 +205.372222,7216591.666667,1179.925053 +217.594444,7216591.666667,1142.517369 +229.816667,7216591.666667,1103.155693 +242.038889,7216591.666667,1061.227727 +254.261111,7216591.666667,1015.844801 +266.483333,7216591.666667,965.586483 +278.705556,7216591.666667,907.849718 +290.927778,7216591.666667,836.658448 +303.150000,7216591.666667,729.094125 +193.150000,8106000.000000,1217.111692 +205.372222,8106000.000000,1181.422931 +217.594444,8106000.000000,1144.294074 +229.816667,8106000.000000,1105.295297 +242.038889,8106000.000000,1063.856717 +254.261111,8106000.000000,1019.166555 +266.483333,8106000.000000,969.961109 +278.705556,8106000.000000,914.018666 +290.927778,8106000.000000,846.619843 +303.150000,8106000.000000,753.272731 diff --git a/docs/requirements.in b/docs/requirements.in index 249681598..c38228d35 100644 --- a/docs/requirements.in +++ b/docs/requirements.in @@ -1,10 +1,10 @@ -nbsphinx==0.9.4 -pydata-sphinx-theme==0.15.2 +nbsphinx==0.9.7 +pydata-sphinx-theme==0.16.1 jupyter-sphinx==0.5.3 -sphinx_design==0.5.0 -ipykernel==6.29.4 -pandas==2.2.2 -lxml_html_clean==0.1.1 +sphinx_design==0.6.1 +ipykernel==6.30.1 +pandas==2.3.2 +lxml_html_clean==0.4.2 sphinx-copybutton==0.5.2 sphinx-tabs==3.4.7 -plotly>=5 +plotly==6.3.0 diff --git a/docs/requirements.txt b/docs/requirements.txt index 14f379cdc..0e9b5ea3d 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -2,87 +2,85 @@ # This file is autogenerated by pip-compile with Python 3.12 # by the following command: # -# pip-compile requirements.in +# pip-compile docs/requirements.in # accessible-pygments==0.0.5 # via pydata-sphinx-theme -alabaster==0.7.16 +alabaster==1.0.0 # via sphinx -asttokens==2.4.1 +asttokens==3.0.0 # via stack-data -attrs==23.2.0 +attrs==25.3.0 # via # jsonschema # referencing -babel==2.15.0 +babel==2.17.0 # via # pydata-sphinx-theme # sphinx -beautifulsoup4==4.12.3 +beautifulsoup4==4.13.5 # via # nbconvert # pydata-sphinx-theme -bleach==6.1.0 +bleach[css]==6.2.0 # via nbconvert -certifi==2024.2.2 +certifi==2025.8.3 # via requests -charset-normalizer==3.3.2 +charset-normalizer==3.4.3 # via requests -comm==0.2.2 +comm==0.2.3 # via # ipykernel # ipywidgets -debugpy==1.8.1 +debugpy==1.8.17 # via ipykernel -decorator==5.1.1 +decorator==5.2.1 # via ipython defusedxml==0.7.1 # via nbconvert -docutils==0.20.1 +docutils==0.21.2 # via # nbsphinx # pydata-sphinx-theme # sphinx # sphinx-tabs -entrypoints==0.4 - # via nbconvert -exceptiongroup==1.2.2 - # via ipython -executing==2.0.1 +executing==2.2.1 # via stack-data -fastjsonschema==2.19.1 +fastjsonschema==2.21.2 # via nbformat -idna==3.7 +idna==3.10 # via requests imagesize==1.4.1 # via sphinx -ipykernel==6.29.4 +ipykernel==6.30.1 # via - # -r requirements.in + # -r docs/requirements.in # jupyter-sphinx -ipython==8.24.0 +ipython==9.5.0 # via # ipykernel # ipywidgets # jupyter-sphinx -ipywidgets==8.1.2 +ipython-pygments-lexers==1.1.1 + # via ipython +ipywidgets==8.1.7 # via jupyter-sphinx -jedi==0.19.1 +jedi==0.19.2 # via ipython -jinja2==3.1.4 +jinja2==3.1.6 # via # nbconvert # nbsphinx # sphinx -jsonschema==4.22.0 +jsonschema==4.25.1 # via nbformat -jsonschema-specifications==2023.12.1 +jsonschema-specifications==2025.9.1 # via jsonschema -jupyter-client==8.6.2 +jupyter-client==8.6.3 # via # ipykernel # nbclient -jupyter-core==5.7.2 +jupyter-core==5.8.1 # via # ipykernel # jupyter-client @@ -90,18 +88,16 @@ jupyter-core==5.7.2 # nbconvert # nbformat jupyter-sphinx==0.5.3 - # via -r requirements.in + # via -r docs/requirements.in jupyterlab-pygments==0.3.0 # via nbconvert -jupyterlab-widgets==3.0.10 +jupyterlab-widgets==3.0.15 # via ipywidgets -lxml==5.2.2 - # via - # lxml-html-clean - # nbconvert -lxml-html-clean==0.1.1 - # via -r requirements.in -markupsafe==2.1.5 +lxml==6.0.1 + # via lxml-html-clean +lxml-html-clean==0.4.2 + # via -r docs/requirements.in +markupsafe==3.0.2 # via # jinja2 # nbconvert @@ -109,11 +105,13 @@ matplotlib-inline==0.1.7 # via # ipykernel # ipython -mistune==0.8.4 +mistune==3.1.4 # via nbconvert -nbclient==0.10.0 +narwhals==2.5.0 + # via plotly +nbclient==0.10.2 # via nbconvert -nbconvert==6.5.4 +nbconvert==7.16.6 # via # jupyter-sphinx # nbsphinx @@ -123,45 +121,45 @@ nbformat==5.10.4 # nbclient # nbconvert # nbsphinx -nbsphinx==0.9.4 - # via -r requirements.in +nbsphinx==0.9.7 + # via -r docs/requirements.in nest-asyncio==1.6.0 # via ipykernel -numpy==1.26.4 +numpy==2.3.3 # via pandas -packaging==24.0 +packaging==25.0 # via # ipykernel # nbconvert # plotly - # pydata-sphinx-theme # sphinx -pandas==2.2.2 - # via -r requirements.in +pandas==2.3.2 + # via -r docs/requirements.in pandocfilters==1.5.1 # via nbconvert -parso==0.8.4 +parso==0.8.5 # via jedi pexpect==4.9.0 # via ipython -platformdirs==4.2.2 +platformdirs==4.4.0 # via jupyter-core -plotly==5.24.1 - # via -r requirements.in -prompt-toolkit==3.0.43 +plotly==6.3.0 + # via -r docs/requirements.in +prompt-toolkit==3.0.52 # via ipython -psutil==5.9.8 +psutil==7.1.0 # via ipykernel ptyprocess==0.7.0 # via pexpect -pure-eval==0.2.2 +pure-eval==0.2.3 # via stack-data -pydata-sphinx-theme==0.15.2 - # via -r requirements.in -pygments==2.18.0 +pydata-sphinx-theme==0.16.1 + # via -r docs/requirements.in +pygments==2.19.2 # via # accessible-pygments # ipython + # ipython-pygments-lexers # nbconvert # pydata-sphinx-theme # sphinx @@ -170,32 +168,29 @@ python-dateutil==2.9.0.post0 # via # jupyter-client # pandas -pytz==2024.1 +pytz==2025.2 # via pandas -pyzmq==26.0.3 +pyzmq==27.1.0 # via # ipykernel # jupyter-client -referencing==0.35.1 +referencing==0.36.2 # via # jsonschema # jsonschema-specifications -requests==2.32.2 +requests==2.32.5 # via sphinx -rpds-py==0.18.1 +rpds-py==0.27.1 # via # jsonschema # referencing -six==1.16.0 - # via - # asttokens - # bleach - # python-dateutil -snowballstemmer==2.2.0 +six==1.17.0 + # via python-dateutil +snowballstemmer==3.0.1 # via sphinx -soupsieve==2.5 +soupsieve==2.8 # via beautifulsoup4 -sphinx==7.3.7 +sphinx==8.1.3 # via # jupyter-sphinx # nbsphinx @@ -204,38 +199,33 @@ sphinx==7.3.7 # sphinx-design # sphinx-tabs sphinx-copybutton==0.5.2 - # via -r requirements.in -sphinx-design==0.5.0 - # via -r requirements.in + # via -r docs/requirements.in +sphinx-design==0.6.1 + # via -r docs/requirements.in sphinx-tabs==3.4.7 - # via -r requirements.in -sphinxcontrib-applehelp==1.0.8 + # via -r docs/requirements.in +sphinxcontrib-applehelp==2.0.0 # via sphinx -sphinxcontrib-devhelp==1.0.6 +sphinxcontrib-devhelp==2.0.0 # via sphinx -sphinxcontrib-htmlhelp==2.0.5 +sphinxcontrib-htmlhelp==2.1.0 # via sphinx sphinxcontrib-jsmath==1.0.1 # via sphinx -sphinxcontrib-qthelp==1.0.7 +sphinxcontrib-qthelp==2.0.0 # via sphinx -sphinxcontrib-serializinghtml==1.1.10 +sphinxcontrib-serializinghtml==2.0.0 # via sphinx stack-data==0.6.3 # via ipython -tenacity==9.0.0 - # via plotly -tinycss2==1.3.0 - # via nbconvert -tomli==2.2.1 - # via sphinx -tornado==6.4 +tinycss2==1.4.0 + # via bleach +tornado==6.5.2 # via # ipykernel # jupyter-client traitlets==5.14.3 # via - # comm # ipykernel # ipython # ipywidgets @@ -246,13 +236,13 @@ traitlets==5.14.3 # nbconvert # nbformat # nbsphinx -typing-extensions==4.12.0 +typing-extensions==4.15.0 # via - # ipython + # beautifulsoup4 # pydata-sphinx-theme -tzdata==2024.1 +tzdata==2025.2 # via pandas -urllib3==2.2.1 +urllib3==2.5.0 # via requests wcwidth==0.2.13 # via prompt-toolkit @@ -260,5 +250,5 @@ webencodings==0.5.1 # via # bleach # tinycss2 -widgetsnbextension==4.0.10 +widgetsnbextension==4.0.14 # via ipywidgets diff --git a/docs/user/first_simulation.rst b/docs/user/first_simulation.rst index 5624ed926..1e341e9b1 100644 --- a/docs/user/first_simulation.rst +++ b/docs/user/first_simulation.rst @@ -276,6 +276,9 @@ Finally, we can add any number of Parachutes to the ``Rocket`` object. sampling_rate=105, lag=1.5, noise=(0, 8.3, 0.5), + radius=1.5, + height=1.5, + porosity=0.0432, ) drogue = calisto.add_parachute( @@ -285,6 +288,9 @@ Finally, we can add any number of Parachutes to the ``Rocket`` object. sampling_rate=105, lag=1.5, noise=(0, 8.3, 0.5), + radius=1.5, + height=1.5, + porosity=0.0432, ) We can then see if the rocket is stable by plotting the static margin: @@ -559,12 +565,16 @@ Visualizing the Trajectory in Google Earth We can export the trajectory to ``.kml`` to visualize it in Google Earth: +Use the dedicated exporter class: + .. jupyter-input:: - test_flight.export_kml( + from rocketpy.simulation import FlightDataExporter + + FlightDataExporter(test_flight).export_kml( file_name="trajectory.kml", extrude=True, - altitude_mode="relative_to_ground", + altitude_mode="relativetoground", ) .. note:: @@ -572,6 +582,10 @@ We can export the trajectory to ``.kml`` to visualize it in Google Earth: To learn more about the ``.kml`` format, see `KML Reference `_. +.. note:: + + The legacy method ``Flight.export_kml`` is deprecated. Use + :meth:`rocketpy.simulation.flight_data_exporter.FlightDataExporter.export_kml`. Manipulating results -------------------- @@ -604,17 +618,27 @@ In this section, we will explore how to export specific data from your RocketPy simulations to CSV files. This is particularly useful if you want to insert the data into spreadsheets or other software for further analysis. -The main method that is used to export data is the :meth:`rocketpy.Flight.export_data` method. This method exports selected flight attributes to a CSV file. In this first example, we will export the rocket angle of attack (see :meth:`rocketpy.Flight.angle_of_attack`) and the rocket mach number (see :meth:`rocketpy.Flight.mach_number`) to the file ``calisto_flight_data.csv``. +The recommended API is +:meth:`rocketpy.simulation.flight_data_exporter.FlightDataExporter.export_data`, +which exports selected flight attributes to a CSV file. In this first example, +we export the rocket angle of attack (see :meth:`rocketpy.Flight.angle_of_attack`) +and the rocket Mach number (see :meth:`rocketpy.Flight.mach_number`) to the file +``calisto_flight_data.csv``. .. jupyter-execute:: - test_flight.export_data( + from rocketpy.simulation import FlightDataExporter + + exporter = FlightDataExporter(test_flight) + exporter.export_data( "calisto_flight_data.csv", "angle_of_attack", "mach_number", ) -| As you can see, the first argument of the method is the name of the file to be created. The following arguments are the attributes to be exported. We can check the file that was created by reading it with the :func:`pandas.read_csv` function: +| As you can see, the first argument is the file name to be created. The following +arguments are the attributes to be exported. We can check the file by reading it +with :func:`pandas.read_csv`: .. jupyter-execute:: @@ -622,11 +646,13 @@ The main method that is used to export data is the :meth:`rocketpy.Flight.export pd.read_csv("calisto_flight_data.csv") -| The file header specifies the meaning of each column. The time samples are obtained from the simulation solver steps. Should you want to export the data at a different sampling rate, you can use the ``time_step`` argument of the :meth:`rocketpy.Flight.export_data` method as follows. +| The file header specifies the meaning of each column. The time samples are +obtained from the simulation solver steps. To export the data at a different +sampling rate, use the ``time_step`` argument: .. jupyter-execute:: - test_flight.export_data( + exporter.export_data( "calisto_flight_data.csv", "angle_of_attack", "mach_number", @@ -635,15 +661,16 @@ The main method that is used to export data is the :meth:`rocketpy.Flight.export pd.read_csv("calisto_flight_data.csv") -This will export the same data at a sampling rate of 1 second. The flight data will be interpolated to match the new sampling rate. +This exports the same data at a sampling rate of 1 second. The flight data is +interpolated to match the new sampling rate. -Finally, the :meth:`rocketpy.Flight.export_data` method also provides a convenient way to export the entire flight solution (see :meth:`rocketpy.Flight.solution_array`) to a CSV file. This is done by not passing any attributes names to the method. +Finally, ``FlightDataExporter.export_data`` also provides a convenient way to +export the entire flight solution (see :meth:`rocketpy.Flight.solution_array`) +by not passing any attribute names: .. jupyter-execute:: - test_flight.export_data( - "calisto_flight_data.csv", - ) + exporter.export_data("calisto_flight_data.csv") .. jupyter-execute:: :hide-code: @@ -653,6 +680,10 @@ Finally, the :meth:`rocketpy.Flight.export_data` method also provides a convenie import os os.remove("calisto_flight_data.csv") +.. note:: + + The legacy method ``Flight.export_data`` is deprecated. Use + :meth:`rocketpy.simulation.flight_data_exporter.FlightDataExporter.export_data`. Saving and Storing Plots ------------------------ @@ -664,7 +695,7 @@ For instance, we can save our rocket drawing as a ``.png`` file: calisto.draw(filename="calisto_drawing.png") -Also, if you want to save a specific rocketpy plot, every RocketPy +Also, if you want to save a specific rocketpy plot, every RocketPy attribute of type :class:`rocketpy.Function` is capable of saving its plot as an image file. For example, we can save our rocket's speed plot and the trajectory plot as ``.jpg`` files: diff --git a/docs/user/flight.rst b/docs/user/flight.rst new file mode 100644 index 000000000..17dde836b --- /dev/null +++ b/docs/user/flight.rst @@ -0,0 +1,613 @@ +.. _flightusage: + +Flight Class Usage +================== + +The :class:`rocketpy.Flight` class is the heart of RocketPy's simulation engine. +It takes a :class:`rocketpy.Rocket`, an :class:`rocketpy.Environment`, and +launch parameters to simulate the complete flight trajectory of a rocket from +launch to landing. + +.. seealso:: + + For a complete example of Flight simulation, see the + :doc:`First Simulation ` guide. + +Creating a Flight Simulation +---------------------------- + +Basic Flight Creation +~~~~~~~~~~~~~~~~~~~~~ + +The most basic way to create a Flight simulation requires three mandatory parameters: +a rocket, an environment, and a rail length. + +.. jupyter-execute:: + :hide-code: + :hide-output: + + import datetime + from rocketpy import Environment, SolidMotor, Rocket, Flight + + # Create a basic environment (Spaceport America) + env = Environment(latitude=32.990254, longitude=-106.974998, elevation=1400) + # Use standard atmosphere to avoid weather data warnings in docs + env.set_atmospheric_model(type="standard_atmosphere") + + # Create a simple solid motor + motor = SolidMotor( + thrust_source="../data/motors/cesaroni/Cesaroni_M1670.eng", + dry_mass=1.815, + dry_inertia=(0.125, 0.125, 0.002), + nozzle_radius=33 / 1000, + grain_number=5, + grain_density=1815, + grain_outer_radius=33 / 1000, + grain_initial_inner_radius=15 / 1000, + grain_initial_height=120 / 1000, + grain_separation=5 / 1000, + grains_center_of_mass_position=0.317, + center_of_dry_mass_position=0.317, + ) + + # Create a simple rocket + rocket = Rocket( + radius=127 / 2000, + mass=14.426, + inertia=(6.321, 6.321, 0.034), + power_off_drag="../data/rockets/calisto/powerOffDragCurve.csv", + power_on_drag="../data/rockets/calisto/powerOnDragCurve.csv", + center_of_mass_without_motor=0, + coordinate_system_orientation="tail_to_nose", + ) + rocket.add_motor(motor, position=-1.255) + +.. jupyter-execute:: + + from rocketpy import Environment, SolidMotor, Rocket, Flight + + # Assuming you have already defined env, motor, and rocket objects + # (See Environment, Motor, and Rocket documentation for details) + + # Create a basic flight + flight = Flight( + rocket=rocket, # Your Rocket object + environment=env, # Your Environment object + rail_length=5.2, # Length of launch rail in meters + ) + +Once created, the Flight object automatically runs the simulation and stores +all results internally. + +Launch Parameters +~~~~~~~~~~~~~~~~~ + +You can customize the launch conditions by specifying additional parameters: + +.. jupyter-execute:: + + flight_custom = Flight( + rocket=rocket, + environment=env, + rail_length=5.2, + inclination=85, # Rail inclination angle (degrees from horizontal) + heading=90, # Launch direction (degrees from North) + ) + +**Key Launch Parameters:** + +- **inclination**: Rail inclination angle relative to ground (0° = horizontal, 90° = vertical) +- **heading**: Launch direction in degrees from North (0° = North, 90° = East, 180° = South, 270° = West) + +Understanding Flight Parameters +------------------------------- + +Coordinate Systems and Positioning +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +RocketPy uses a launch-centered coordinate system: + +- **X-axis**: Points East (positive values = East direction) +- **Y-axis**: Points North (positive values = North direction) +- **Z-axis**: Points upward (positive values = altitude above launch site) + +The rocket's position is tracked relative to the launch site throughout the flight. + +Initial Conditions +~~~~~~~~~~~~~~~~~~ + +The Flight class automatically calculates appropriate initial conditions based on: + +- Rail inclination and heading angles +- Rocket orientation (affected by rail button positions if present) +- Environment conditions (wind, atmospheric properties) + +You can also specify custom initial conditions by passing an ``initial_solution`` +array or another Flight object to continue from a previous state. + +**Custom Initial Solution Vector** + +The ``initial_solution`` parameter accepts a 14-element array defining the complete +initial state of the rocket: + +.. code-block:: python + + initial_solution = [ + t_initial, # Initial time (s) + x_init, # Initial X position - East coordinate (m) + y_init, # Initial Y position - North coordinate (m) + z_init, # Initial Z position - altitude above launch site (m) + vx_init, # Initial velocity in X direction - East (m/s) + vy_init, # Initial velocity in Y direction - North (m/s) + vz_init, # Initial velocity in Z direction - upward (m/s) + e0_init, # Initial Euler parameter 0 (quaternion scalar part) + e1_init, # Initial Euler parameter 1 (quaternion i component) + e2_init, # Initial Euler parameter 2 (quaternion j component) + e3_init, # Initial Euler parameter 3 (quaternion k component) + w1_init, # Initial angular velocity about rocket's x-axis (rad/s) + w2_init, # Initial angular velocity about rocket's y-axis (rad/s) + w3_init # Initial angular velocity about rocket's z-axis (rad/s) + ] + +**Using a Previous Flight as Initial Condition** + +You can also continue a simulation from where another flight ended: + +.. code-block:: python + + # Continue from the final state of a previous flight + continued_flight = Flight( + rocket=rocket, + environment=env, + rail_length=0, # Set to 0 when continuing from free flight + initial_solution=flight # Use previous Flight object + ) + +This is particularly useful for multi-stage simulations or when analyzing +different scenarios from a specific flight condition. + +Simulation Control Parameters +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +**Time Control:** + +- ``max_time``: Maximum simulation duration (default: 600 seconds) +- ``max_time_step``: Maximum integration step size (default: infinity) +- ``min_time_step``: Minimum integration step size (default: 0) + +.. note:: + + These time control parameters can significantly help with integration stability \ + in challenging simulation cases. This is particularly useful for liquid and \ + hybrid motors, which often have more complex thrust curves and transient \ + behaviors that can cause numerical integration difficulties. + +**Accuracy Control:** + +- ``rtol``: Relative tolerance for numerical integration (default: 1e-6) +- ``atol``: Absolute tolerance for numerical integration (default: auto-calculated) + +.. note:: + + Increasing the tolerance values can speed up simulations but may reduce accuracy. + +**Simulation Behavior:** + +- ``terminate_on_apogee``: Stop simulation at apogee (default: False) +- ``time_overshoot``: Allow time step overshoot for efficiency (default: True) + +Accessing Simulation Results +---------------------------- + +Once a Flight simulation is complete, you can access a wealth of data about +the rocket's trajectory and performance. + +Trajectory Data +~~~~~~~~~~~~~~~ + +Basic position and velocity data: + +.. jupyter-execute:: + + # Position coordinates (as functions of time) + x_trajectory = flight.x # East coordinate (m) + y_trajectory = flight.y # North coordinate (m) + altitude = flight.z # Altitude above launch site (m) + + # Velocity components (as functions of time) + vx = flight.vx # East velocity (m/s) + vy = flight.vy # North velocity (m/s) + vz = flight.vz # Vertical velocity (m/s) + + # Access specific values at given times + altitude_at_10s = flight.z(10) # Altitude at t=10 seconds + max_altitude = flight.apogee # Maximum altitude reached + +Key Flight Events +~~~~~~~~~~~~~~~~~ + +Important events during the flight: + +.. jupyter-execute:: + + # Rail departure + rail_departure_time = flight.out_of_rail_time + rail_departure_velocity = flight.out_of_rail_velocity + + # Apogee + apogee_time = flight.apogee_time + apogee_altitude = flight.apogee + apogee_coordinates = (flight.apogee_x, flight.apogee_y) + + # Landing/Impact + impact_time = flight.impact_state[0] + impact_velocity = flight.impact_velocity + impact_coordinates = (flight.x_impact, flight.y_impact) + +Forces and Accelerations +~~~~~~~~~~~~~~~~~~~~~~~~ + +The Flight object provides access to all forces and accelerations acting on the rocket: + +.. jupyter-execute:: + + # Linear accelerations in inertial frame (m/s²) + ax = flight.ax # East acceleration + ay = flight.ay # North acceleration + az = flight.az # Vertical acceleration + + # Aerodynamic forces in body frame (N) + R1 = flight.R1 # X-axis aerodynamic force + R2 = flight.R2 # Y-axis aerodynamic force + R3 = flight.R3 # Z-axis aerodynamic force (drag) + + # Aerodynamic moments in body frame (N⋅m) + M1 = flight.M1 # Roll moment + M2 = flight.M2 # Pitch moment + M3 = flight.M3 # Yaw moment + +Attitude and Orientation +~~~~~~~~~~~~~~~~~~~~~~~~ + +Rocket orientation throughout the flight: + +.. jupyter-execute:: + + # Euler parameters (quaternions) + e0, e1, e2, e3 = flight.e0, flight.e1, flight.e2, flight.e3 + + # Angular velocities in body frame (rad/s) + w1 = flight.w1 # Roll rate + w2 = flight.w2 # Pitch rate + w3 = flight.w3 # Yaw rate + + # Derived attitude angles + attitude_angle = flight.attitude_angle + path_angle = flight.path_angle + +Performance Metrics +~~~~~~~~~~~~~~~~~~~ + +Key performance indicators: + +.. jupyter-execute:: + + # Velocity and speed + total_speed = flight.speed + mach_number = flight.mach_number + + # Stability indicators + static_margin = flight.static_margin + stability_margin = flight.stability_margin + +Accessing Raw Simulation Data +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +For users who need direct access to the raw numerical simulation results, +the Flight object provides the complete solution array through the ``solution`` +and ``solution_array`` attributes. + +**Flight.solution** + +The ``Flight.solution`` attribute contains the raw simulation data as a list of +state vectors, where each row represents the rocket's complete state at a specific time: + +.. jupyter-execute:: + + # Access the raw solution list + raw_solution = flight.solution + + # Each element is a 14-element state vector: + # [time, x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3] + initial_state = flight.solution[0] # First time step + final_state = flight.solution[-1] # Last time step + + print(f"Initial state: {initial_state}") + print(f"Final state: {final_state}") + +**Flight.solution_array** + +For easier numerical analysis, use ``solution_array`` which provides the same data +as a NumPy array: + +.. jupyter-execute:: + + import numpy as np + + # Get solution as NumPy array for easier manipulation + solution_array = flight.solution_array # Shape: (n_time_steps, 14) + + # Extract specific columns (state variables) + time_array = solution_array[:, 0] # Time values + position_data = solution_array[:, 1:4] # X, Y, Z positions + velocity_data = solution_array[:, 4:7] # Vx, Vy, Vz velocities + quaternions = solution_array[:, 7:11] # e0, e1, e2, e3 + angular_velocities = solution_array[:, 11:14] # w1, w2, w3 + + # Example: Calculate velocity magnitude manually + velocity_magnitude = np.sqrt(np.sum(velocity_data**2, axis=1)) + +**State Vector Format** + +Each row in the solution array follows this 14-element format: + +.. code-block:: python + + [time, x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3] + +Where: + - ``time``: Simulation time in seconds + - ``x, y, z``: Position coordinates in meters (East, North, Up) + - ``vx, vy, vz``: Velocity components in m/s (East, North, Up) + - ``e0, e1, e2, e3``: Euler parameters (quaternions) for attitude + - ``w1, w2, w3``: Angular velocities in rad/s (body frame: roll, pitch, yaw rates) + +**Getting State at Specific Time** + +You can extract the rocket's state at any specific time during the flight: + +.. jupyter-execute:: + + # Get complete state vector at t=10 seconds + state_at_10s = flight.get_solution_at_time(10.0) + + print(f"State at t=10s: {state_at_10s}") + + # Extract specific values from the state vector + time_10s = state_at_10s[0] + altitude_10s = state_at_10s[3] # Z coordinate + speed_10s = np.sqrt(state_at_10s[4]**2 + state_at_10s[5]**2 + state_at_10s[6]**2) + + print(f"At t={time_10s}s: altitude={altitude_10s:.1f}m, speed={speed_10s:.1f}m/s") + +This raw data access is particularly useful for: + +- Custom post-processing and analysis +- Exporting data to external tools +- Implementing custom flight metrics +- Monte Carlo analysis and statistical studies +- Integration with other simulation frameworks + + +Plotting Flight Data +-------------------- + +The Flight class provides comprehensive plotting capabilities through the +``plots`` attribute. + +Trajectory Visualization +~~~~~~~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + # 3D trajectory plot + flight.plots.trajectory_3d() + + # 2D trajectory (ground track) + flight.plots.linear_kinematics_data() + +Flight Data Plots +~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + # Velocity and acceleration plots + flight.plots.linear_kinematics_data() + + # Attitude and angular motion + flight.plots.attitude_data() + flight.plots.angular_kinematics_data() + + # Flight path and orientation + flight.plots.flight_path_angle_data() + +Forces and Moments +~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + # Aerodynamic forces + flight.plots.aerodynamic_forces() + + # Rail button forces (if applicable) + flight.plots.rail_buttons_forces() + +Energy Analysis +~~~~~~~~~~~~~~~ + +.. code-block:: python + + # Energy plots + flight.plots.energy_data() + + # Stability analysis + flight.plots.stability_and_control_data() + +Comprehensive Analysis +~~~~~~~~~~~~~~~~~~~~~~ + +For a complete overview of all plots: + +.. jupyter-execute:: + + # Show all available plots + flight.all_info() + +Printing Flight Information +--------------------------- + +The Flight class also provides detailed text output through the ``prints`` attribute. + +Flight Summary +~~~~~~~~~~~~~~ + +.. code-block:: python + + # Complete flight information + flight.info() + + # All detailed information + flight.all_info() + +Specific Information Sections +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + # Initial conditions + flight.prints.initial_conditions() + + # Wind and environment conditions + flight.prints.surface_wind_conditions() + + # Launch rail information + flight.prints.launch_rail_conditions() + + # Rail departure conditions + flight.prints.out_of_rail_conditions() + + # Motor burn out conditions + flight.prints.burn_out_conditions() + + # Apogee conditions + flight.prints.apogee_conditions() + + # Landing/impact conditions + flight.prints.impact_conditions() + + # Maximum values during flight + flight.prints.maximum_values() + +Advanced Features +----------------- + +Custom Equations of Motion +~~~~~~~~~~~~~~~~~~~~~~~~~~ + +RocketPy supports different sets of equations of motion: + +.. code-block:: python + + # Standard 6-DOF equations (default) + flight_6dof = Flight( + rocket=rocket, + environment=env, + rail_length=5.2, + equations_of_motion="standard" + ) + + # Simplified solid propulsion equations (legacy) + # This may run a bit faster with no accuracy loss, but only works for solid motors + flight_simple = Flight( + rocket=rocket, + environment=env, + rail_length=5.2, + equations_of_motion="solid_propulsion" + ) + +Integration Method Selection +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +You can choose different numerical integration methods using the ``ode_solver`` parameter. +RocketPy supports the following integration methods from ``scipy.integrate.solve_ivp``: + +**Available ODE Solvers:** + +- **'LSODA'** (default): Recommended for most flights. Automatically switches between stiff and non-stiff methods +- **'RK45'**: Explicit Runge-Kutta method of order 5(4). Good for non-stiff problems +- **'RK23'**: Explicit Runge-Kutta method of order 3(2). Faster but less accurate than RK45 +- **'DOP853'**: Explicit Runge-Kutta method of order 8. High accuracy for smooth problems +- **'Radau'**: Implicit Runge-Kutta method of order 5. Good for stiff problems +- **'BDF'**: Implicit multi-step variable-order method. Efficient for stiff problems + +.. code-block:: python + + # High-accuracy integration (default, recommended for most cases) + flight_default = Flight( + rocket=rocket, + environment=env, + rail_length=5.2, + ode_solver="LSODA" + ) + + # Fast integration for quick simulations + flight_fast = Flight( + rocket=rocket, + environment=env, + rail_length=5.2, + ode_solver="RK45" + ) + + # Very high accuracy for smooth problems + flight_high_accuracy = Flight( + rocket=rocket, + environment=env, + rail_length=5.2, + ode_solver="DOP853" + ) + + # For stiff problems (e.g., complex motor thrust curves) + flight_stiff = Flight( + rocket=rocket, + environment=env, + rail_length=5.2, + ode_solver="BDF" + ) + +You can also pass a custom ``scipy.integrate.OdeSolver`` object for advanced use cases. +For more information on integration methods, see the `scipy documentation +`_. + +Exporting Flight Data +--------------------- + +You can export flight data for external analysis: + +.. code-block:: python + + # Convert to dictionary format + flight_data = flight.to_dict(include_outputs=True) + # NOTE: RocketPy offers an unofficial json serializer, see rocketpy._encoders for details + + +Common Issues and Solutions +--------------------------- + +**Integration Problems:** + +- Reduce ``max_time_step`` for more accuracy +- Increase ``rtol`` for faster but less accurate simulations +- Check rocket and environment definitions for unrealistic values + +**Missing Events:** + +- Ensure ``max_time`` is sufficient for complete flight +- Verify parachute trigger conditions +- Check for premature termination conditions +- You can set ``verbose=True`` in the Flight constructor to get detailed logs during simulation + +**Performance Issues:** + +- Set ``time_overshoot=True`` for better performance +- Use simpler integration methods for quick runs +- Consider reducing the complexity of atmospheric models diff --git a/docs/user/index.rst b/docs/user/index.rst index a09872ad5..8620e3ad5 100644 --- a/docs/user/index.rst +++ b/docs/user/index.rst @@ -16,6 +16,7 @@ RocketPy's User Guide Motors Rocket Environment + Flight .. toctree:: :maxdepth: 2 diff --git a/docs/user/installation.rst b/docs/user/installation.rst index d19a42ff3..dcc188d17 100644 --- a/docs/user/installation.rst +++ b/docs/user/installation.rst @@ -85,7 +85,7 @@ Requirements Python Version ^^^^^^^^^^^^^^ -RocketPy supports Python 3.9 and above. +RocketPy supports Python 3.10 and above. Sorry, there are currently no plans to support earlier versions. If you really need to run RocketPy on Python 3.8 or earlier, feel free to submit an issue and we will see what we can do! diff --git a/docs/user/motors/tanks.rst b/docs/user/motors/tanks.rst index 49b5a0687..f1a3342b8 100644 --- a/docs/user/motors/tanks.rst +++ b/docs/user/motors/tanks.rst @@ -73,6 +73,21 @@ density as such: Fluid are then passed to tanks when they are defined. +.. note:: + + One may define the fluid density as a function of temperature (K) and + pressure (Pa). The data can be imported from an external source, such as + a dataset or external libraries. + In this case, the fluid would be defined as such: + + >>> Fluid(name="N2O", density=lambda t, p: 44 * p / (8.314 * t)) + >>> from CoolProp.CoolProp import PropsSI # external library + >>> Fluid(name="N2O", density=lambda t, p: PropsSI('D', 'T', t, 'P', p, 'N2O')) + + In fact, the density parameter can be any ``Function`` source, such as a + ``callable``, csv file or an array of points. See more on + :class:`rocketpy.Function`. + Tank Geometry ------------- diff --git a/docs/user/rocket/rocket_usage.rst b/docs/user/rocket/rocket_usage.rst index 1b4d2bbd5..3496da6fd 100644 --- a/docs/user/rocket/rocket_usage.rst +++ b/docs/user/rocket/rocket_usage.rst @@ -302,6 +302,9 @@ apogee and another that will be deployed at 800 meters above ground level: sampling_rate=105, lag=1.5, noise=(0, 8.3, 0.5), + radius=1.5, + height=1.5, + porosity=0.0432, ) drogue = calisto.add_parachute( @@ -311,6 +314,9 @@ apogee and another that will be deployed at 800 meters above ground level: sampling_rate=105, lag=1.5, noise=(0, 8.3, 0.5), + radius=1.5, + height=1.5, + porosity=0.0432, ) .. seealso:: diff --git a/pyproject.toml b/pyproject.toml index f471085ed..23d6b3998 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ version = "1.10.0" description="Advanced 6-DOF trajectory simulation for High-Power Rocketry." dynamic = ["dependencies"] readme = "README.md" -requires-python = ">=3.9" +requires-python = ">=3.10" authors = [ {name = "Giovani Hidalgo Ceotto", email = "ghceotto@gmail.com"}, {name = "Guilherme Fernandes Alves", email = "guilherme_fernandes@usp.br"}, diff --git a/requirements.txt b/requirements.txt index 68aebe503..61a594320 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,6 +1,6 @@ numpy>=1.13 scipy>=1.0 -matplotlib>=3.0 +matplotlib>=3.9.0 # Released May 15th 2024 netCDF4>=1.6.4 requests pytz diff --git a/rocketpy/_encoders.py b/rocketpy/_encoders.py index 141ff6baa..58eeae809 100644 --- a/rocketpy/_encoders.py +++ b/rocketpy/_encoders.py @@ -13,33 +13,46 @@ class RocketPyEncoder(json.JSONEncoder): """Custom JSON encoder for RocketPy objects. It defines how to encode - different types of objects to a JSON supported format.""" + different types of objects to a JSON supported format. + """ def __init__(self, *args, **kwargs): + """Initializes the encoder with parameter options. + + Parameters + ---------- + *args : tuple + Positional arguments to pass to the parent class. + **kwargs : dict + Keyword arguments to configure the encoder. The following + options are available: + - include_outputs: bool, whether to include simulation outputs. + Default is False. + - include_function_data: bool, whether to include Function + data in the encoding. If False, Functions will be encoded by their + ``__repr__``. This is useful for reducing the size of the outputs, + but it prevents full restoration of the object upon decoding. + Default is True. + - discretize: bool, whether to discretize Functions whose source + are callables. If True, the accuracy of the decoding may be reduced. + Default is False. + - allow_pickle: bool, whether to pickle callable objects. If + False, callable sources (such as user-defined functions, parachute + triggers or simulation callable outputs) will have their name + stored instead of the function itself. This is useful for + reducing the size of the outputs, but it prevents full restoration + of the object upon decoding. + Default is True. + """ self.include_outputs = kwargs.pop("include_outputs", False) self.include_function_data = kwargs.pop("include_function_data", True) + self.discretize = kwargs.pop("discretize", False) + self.allow_pickle = kwargs.pop("allow_pickle", True) super().__init__(*args, **kwargs) def default(self, o): - if isinstance( - o, - ( - np.int_, - np.intc, - np.intp, - np.int8, - np.int16, - np.int32, - np.int64, - np.uint8, - np.uint16, - np.uint32, - np.uint64, - ), - ): - return int(o) - elif isinstance(o, (np.float16, np.float32, np.float64)): - return float(o) + if isinstance(o, np.generic): + return o.item() elif isinstance(o, np.ndarray): return o.tolist() elif isinstance(o, datetime): @@ -50,11 +63,19 @@ def default(self, o): if not self.include_function_data: return str(o) else: - encoding = o.to_dict(self.include_outputs) + encoding = o.to_dict( + include_outputs=self.include_outputs, + discretize=self.discretize, + allow_pickle=self.allow_pickle, + ) encoding["signature"] = get_class_signature(o) return encoding elif hasattr(o, "to_dict"): - encoding = o.to_dict(self.include_outputs) + encoding = o.to_dict( + include_outputs=self.include_outputs, + discretize=self.discretize, + allow_pickle=self.allow_pickle, + ) encoding = remove_circular_references(encoding) encoding["signature"] = get_class_signature(o) @@ -86,15 +107,21 @@ def object_hook(self, obj): try: class_ = get_class_from_signature(signature) + hash_ = signature.get("hash", None) if class_.__name__ == "Flight" and not self.resimulate: new_flight = class_.__new__(class_) new_flight.prints = _FlightPrints(new_flight) new_flight.plots = _FlightPlots(new_flight) set_minimal_flight_attributes(new_flight, obj) + if hash_ is not None: + setattr(new_flight, "__rpy_hash", hash_) return new_flight elif hasattr(class_, "from_dict"): - return class_.from_dict(obj) + new_obj = class_.from_dict(obj) + if hash_ is not None: + setattr(new_obj, "__rpy_hash", hash_) + return new_obj else: # Filter keyword arguments kwargs = { @@ -102,8 +129,10 @@ def object_hook(self, obj): for key, value in obj.items() if key in class_.__init__.__code__.co_varnames } - - return class_(**kwargs) + new_obj = class_(**kwargs) + if hash_ is not None: + setattr(new_obj, "__rpy_hash", hash_) + return new_obj except (ImportError, AttributeError): return obj else: @@ -136,7 +165,6 @@ def set_minimal_flight_attributes(flight, obj): "x_impact", "y_impact", "t_final", - "flight_phases", "ax", "ay", "az", @@ -186,7 +214,14 @@ def get_class_signature(obj): class_ = obj.__class__ name = getattr(class_, "__qualname__", class_.__name__) - return {"module": class_.__module__, "name": name} + signature = {"module": class_.__module__, "name": name} + + try: + signature.update({"hash": hash(obj)}) + except TypeError: + pass + + return signature def get_class_from_signature(signature): diff --git a/rocketpy/control/controller.py b/rocketpy/control/controller.py index 8338e05b4..27ad62361 100644 --- a/rocketpy/control/controller.py +++ b/rocketpy/control/controller.py @@ -1,4 +1,7 @@ from inspect import signature +from typing import Iterable + +from rocketpy.tools import from_hex_decode, to_hex_encode from ..prints.controller_prints import _ControllerPrints @@ -181,3 +184,46 @@ def info(self): def all_info(self): """Prints out all information about the controller.""" self.info() + + def to_dict(self, **kwargs): + allow_pickle = kwargs.get("allow_pickle", True) + + if allow_pickle: + controller_function = to_hex_encode(self.controller_function) + else: + controller_function = self.controller_function.__name__ + + return { + "controller_function": controller_function, + "sampling_rate": self.sampling_rate, + "initial_observed_variables": self.initial_observed_variables, + "name": self.name, + "_interactive_objects_hash": hash(self.interactive_objects) + if not isinstance(self.interactive_objects, Iterable) + else [hash(obj) for obj in self.interactive_objects], + } + + @classmethod + def from_dict(cls, data): + interactive_objects = data.get("interactive_objects", []) + controller_function = data.get("controller_function") + sampling_rate = data.get("sampling_rate") + initial_observed_variables = data.get("initial_observed_variables") + name = data.get("name", "Controller") + + try: + controller_function = from_hex_decode(controller_function) + except (TypeError, ValueError): + pass + + obj = cls( + interactive_objects=interactive_objects, + controller_function=controller_function, + sampling_rate=sampling_rate, + initial_observed_variables=initial_observed_variables, + name=name, + ) + setattr( + obj, "_interactive_objects_hash", data.get("_interactive_objects_hash", []) + ) + return obj diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py index a0a8c4238..eb2eacd5a 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -27,6 +27,7 @@ find_latitude_index, find_longitude_index, find_time_index, + geodesic_to_utm, get_elevation_data_from_dataset, get_final_date_from_time_array, get_initial_date_from_time_array, @@ -34,8 +35,6 @@ get_pressure_levels_from_file, mask_and_clean_dataset, ) -from rocketpy.environment.tools import geodesic_to_utm as geodesic_to_utm_tools -from rocketpy.environment.tools import utm_to_geodesic as utm_to_geodesic_tools from rocketpy.environment.weather_model_mapping import WeatherModelMapping from rocketpy.mathutils.function import NUMERICAL_TYPES, Function, funcify_method from rocketpy.plots.environment_plots import _EnvironmentPlots @@ -248,6 +247,8 @@ class Environment: Number of ensemble members. Only defined when using Ensembles. Environment.ensemble_member : int Current selected ensemble member. Only defined when using Ensembles. + Environment.earth_rotation_vector : list[float] + Earth's angular velocity vector in the Flight Coordinate System. Notes ----- @@ -353,6 +354,7 @@ def __init__( self.set_location(latitude, longitude) self.__initialize_earth_geometry(datum) self.__initialize_utm_coordinates() + self.__set_earth_rotation_vector() # Set the gravity model self.gravity = self.set_gravity_model(gravity) @@ -451,7 +453,7 @@ def __initialize_utm_coordinates(self): self.initial_utm_letter, self.initial_hemisphere, self.initial_ew, - ) = self.geodesic_to_utm( + ) = geodesic_to_utm( lat=self.latitude, lon=self.longitude, flattening=self.ellipsoid.flattening, @@ -584,6 +586,23 @@ def __reset_wind_direction_function(self): self.wind_direction.set_outputs("Wind Direction (Deg True)") self.wind_direction.set_title("Wind Direction Profile") + def __set_earth_rotation_vector(self): + """Calculates and stores the Earth's angular velocity vector in the Flight + Coordinate System, which is essential for evaluating inertial forces. + """ + # Sidereal day + T = 86164.1 # seconds + + # Earth's angular velocity magnitude + w_earth = 2 * np.pi / T + + # Vector in the Flight Coordinate System + lat = np.radians(self.latitude) + w_local = [0, w_earth * np.cos(lat), w_earth * np.sin(lat)] + + # Store the attribute + self.earth_rotation_vector = w_local + # Validators (used to verify an attribute is being set correctly.) def __validate_dictionary(self, file, dictionary): @@ -2523,98 +2542,7 @@ def set_earth_geometry(self, datum): f"the following recognized datum: {available_datums}" ) from e - # Auxiliary functions - Geodesic Coordinates - - @staticmethod - def geodesic_to_utm( - lat, lon, semi_major_axis=6378137.0, flattening=1 / 298.257223563 - ): - """Function which converts geodetic coordinates, i.e. lat/lon, to UTM - projection coordinates. Can be used only for latitudes between -80.00° - and 84.00° - - Parameters - ---------- - lat : float - The latitude coordinates of the point of analysis, must be contained - between -80.00° and 84.00° - lon : float - The longitude coordinates of the point of analysis, must be - contained between -180.00° and 180.00° - semi_major_axis : float - The semi-major axis of the ellipsoid used to represent the Earth, - must be given in meters (default is 6,378,137.0 m, which corresponds - to the WGS84 ellipsoid) - flattening : float - The flattening of the ellipsoid used to represent the Earth, usually - between 1/250 and 1/150 (default is 1/298.257223563, which - corresponds to the WGS84 ellipsoid) - - Returns - ------- - x : float - East coordinate, always positive - y : float - North coordinate, always positive - utm_zone : int - The number of the UTM zone of the point of analysis, can vary - between 1 and 60 - utm_letter : string - The letter of the UTM zone of the point of analysis, can vary - between C and X, omitting the letters "I" and "O" - hemis : string - Returns "S" for southern hemisphere and "N" for Northern hemisphere - EW : string - Returns "W" for western hemisphere and "E" for eastern hemisphere - """ - warnings.warn( - "This function is deprecated and will be removed in v1.10.0. " - "Please use the new method `tools.geodesic_to_utm` instead.", - DeprecationWarning, - ) - return geodesic_to_utm_tools(lat, lon, semi_major_axis, flattening) - - @staticmethod - def utm_to_geodesic( - x, y, utm_zone, hemis, semi_major_axis=6378137.0, flattening=1 / 298.257223563 - ): - """Function to convert UTM coordinates to geodesic coordinates - (i.e. latitude and longitude). - - Parameters - ---------- - x : float - East UTM coordinate in meters - y : float - North UTM coordinate in meters - utm_zone : int - The number of the UTM zone of the point of analysis, can vary - between 1 and 60 - hemis : string - Equals to "S" for southern hemisphere and "N" for Northern - hemisphere - semi_major_axis : float - The semi-major axis of the ellipsoid used to represent the Earth, - must be given in meters (default is 6,378,137.0 m, which corresponds - to the WGS84 ellipsoid) - flattening : float - The flattening of the ellipsoid used to represent the Earth, usually - between 1/250 and 1/150 (default is 1/298.257223563, which - corresponds to the WGS84 ellipsoid) - - Returns - ------- - lat : float - latitude of the analyzed point - lon : float - latitude of the analyzed point - """ - warnings.warn( - "This function is deprecated and will be removed in v1.10.0. " - "Please use the new method `tools.utm_to_geodesic` instead.", - DeprecationWarning, - ) - return utm_to_geodesic_tools(x, y, utm_zone, hemis, semi_major_axis, flattening) + # Auxiliary functions @staticmethod def calculate_earth_radius( @@ -2702,7 +2630,21 @@ def decimal_degrees_to_arc_seconds(angle): arc_seconds = (remainder * 60 - arc_minutes) * 60 return degrees, arc_minutes, arc_seconds - def to_dict(self, include_outputs=False): + def to_dict(self, **kwargs): + wind_velocity_x = self.wind_velocity_x + wind_velocity_y = self.wind_velocity_y + wind_heading = self.wind_heading + wind_direction = self.wind_direction + wind_speed = self.wind_speed + density = self.density + if kwargs.get("discretize", False): + wind_velocity_x = wind_velocity_x.set_discrete(0, self.max_expected_height) + wind_velocity_y = wind_velocity_y.set_discrete(0, self.max_expected_height) + wind_heading = wind_heading.set_discrete(0, self.max_expected_height) + wind_direction = wind_direction.set_discrete(0, self.max_expected_height) + wind_speed = wind_speed.set_discrete(0, self.max_expected_height) + density = density.set_discrete(0, self.max_expected_height) + env_dict = { "gravity": self.gravity, "date": self.date, @@ -2715,15 +2657,15 @@ def to_dict(self, include_outputs=False): "atmospheric_model_type": self.atmospheric_model_type, "pressure": self.pressure, "temperature": self.temperature, - "wind_velocity_x": self.wind_velocity_x, - "wind_velocity_y": self.wind_velocity_y, - "wind_heading": self.wind_heading, - "wind_direction": self.wind_direction, - "wind_speed": self.wind_speed, + "wind_velocity_x": wind_velocity_x, + "wind_velocity_y": wind_velocity_y, + "wind_heading": wind_heading, + "wind_direction": wind_direction, + "wind_speed": wind_speed, } - if include_outputs: - env_dict["density"] = self.density + if kwargs.get("include_outputs", False): + env_dict["density"] = density env_dict["barometric_height"] = self.barometric_height env_dict["speed_of_sound"] = self.speed_of_sound env_dict["dynamic_viscosity"] = self.dynamic_viscosity diff --git a/rocketpy/environment/weather_model_mapping.py b/rocketpy/environment/weather_model_mapping.py index 59c6e2fd9..c8596d705 100644 --- a/rocketpy/environment/weather_model_mapping.py +++ b/rocketpy/environment/weather_model_mapping.py @@ -32,6 +32,7 @@ class WeatherModelMapping: "latitude": "latitude", "longitude": "longitude", "level": "level", + "ensemble": "number", "temperature": "t", "surface_geopotential_height": None, "geopotential_height": None, @@ -44,6 +45,7 @@ class WeatherModelMapping: "latitude": "latitude", "longitude": "longitude", "level": "pressure_level", + "ensemble": "number", "temperature": "t", "surface_geopotential_height": None, "geopotential_height": None, diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index e60b39286..9fe343687 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -10,10 +10,10 @@ from bisect import bisect_left from collections.abc import Iterable from copy import deepcopy +from enum import Enum from functools import cached_property from inspect import signature from pathlib import Path -from enum import Enum import matplotlib.pyplot as plt import numpy as np @@ -24,9 +24,8 @@ RBFInterpolator, ) -from rocketpy.tools import from_hex_decode, to_hex_encode - -from ..plots.plot_helpers import show_or_save_plot +from rocketpy.plots.plot_helpers import show_or_save_plot +from rocketpy.tools import deprecated, from_hex_decode, to_hex_encode # Numpy 1.x compatibility, # TODO: remove these lines when all dependencies support numpy>=2.0.0 @@ -153,6 +152,7 @@ def __init__( self.__extrapolation__ = extrapolation self.title = title self.__img_dim__ = 1 # always 1, here for backwards compatibility + self.__cropped_domain__ = (None, None) # the x interval if cropped # args must be passed from self. self.set_source(self.source) @@ -449,6 +449,9 @@ def rbf_interpolation(x, x_min, x_max, x_data, y_data, coeffs): # pylint: disab self._interpolation_func = rbf_interpolation + else: + raise ValueError(f"Interpolation {interpolation} method not recognized.") + def __set_extrapolation_func(self): # pylint: disable=too-many-statements """Defines extrapolation function used by the Function. Each extrapolation method has its own function. The function is stored in @@ -532,6 +535,11 @@ def natural_extrapolation(x, x_min, x_max, x_data, y_data, _): def natural_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): # pylint: disable=unused-argument return interpolator(x) + else: + raise ValueError( + f"Natural extrapolation not defined for {interpolation}." + ) + self._extrapolation_func = natural_extrapolation elif extrapolation == 2: # constant if self.__dom_dim__ == 1: @@ -547,6 +555,8 @@ def constant_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): return extrapolator(x) self._extrapolation_func = constant_extrapolation + else: + raise ValueError(f"Extrapolation {extrapolation} method not recognized.") def set_get_value_opt(self): """Defines a method that evaluates interpolations. @@ -625,10 +635,121 @@ def __get_value_opt_nd(self, *args): return result + def __determine_1d_domain_bounds(self, lower, upper): + """Determine domain bounds for 1-D function discretization. + + Parameters + ---------- + lower : scalar, optional + Lower bound. If None, will use cropped domain or default. + upper : scalar, optional + Upper bound. If None, will use cropped domain or default. + + Returns + ------- + tuple + (lower_bound, upper_bound) for the domain. + """ + domain = [0, 10] # default boundaries + cropped = self.__cropped_domain__ + + if cropped[0] is not None and cropped[0] > domain[0]: + domain[0] = cropped[0] + + if cropped[1] is not None and cropped[1] < domain[1]: + domain[1] = cropped[1] + + # Input bounds have preference + domain[0] = lower if lower is not None else domain[0] + domain[1] = upper if upper is not None else domain[1] + + return domain + + def __determine_2d_domain_bounds(self, lower, upper, samples): + """Determine domain bounds for 2-D function discretization. + + Parameters + ---------- + lower : scalar or list, optional + Lower bounds. If None, will use cropped domain or default. + upper : scalar or list, optional + Upper bounds. If None, will use cropped domain or default. + samples : int or list + Number of samples for each dimension. + + Returns + ------- + tuple + (lower_bounds, upper_bounds, sample_counts) for the 2D domain. + """ + default_bounds = [[0, 10], [0, 10]] + + # Apply cropped domain constraints if they exist + final_bounds = deepcopy(default_bounds) + if self.__cropped_domain__ is not None: + for dim in range(2): + cropped_limits = self.__cropped_domain__[dim] + if cropped_limits is not None: + # Use the more restrictive bounds (cropped domain takes precedence) + final_bounds[dim][0] = max( + default_bounds[dim][0], cropped_limits[0] + ) + final_bounds[dim][1] = min( + default_bounds[dim][1], cropped_limits[1] + ) + + # Convert parameters to consistent list format + lower_bounds = self.__normalize_2d_parameter( + lower, [final_bounds[0][0], final_bounds[1][0]] + ) + upper_bounds = self.__normalize_2d_parameter( + upper, [final_bounds[0][1], final_bounds[1][1]] + ) + sample_counts = self.__normalize_2d_parameter(samples, samples) + + return lower_bounds, upper_bounds, sample_counts + + def __normalize_2d_parameter(self, param, default_values): + if param is None: + return ( + default_values + if isinstance(default_values, list) + else [default_values, default_values] + ) + + if isinstance(param, NUMERICAL_TYPES): + return [param, param] + + return param + + def __discretize_1d_function( + self, func, lower, upper, samples, interpolation, extrapolation, one_by_one + ): + lower, upper = self.__determine_1d_domain_bounds(lower, upper) + xs = np.linspace(lower, upper, samples) + ys = func.get_value(xs.tolist()) if one_by_one else func.get_value(xs) + func.__interpolation__ = interpolation + func.__extrapolation__ = extrapolation + func.set_source(np.column_stack((xs, ys))) + + def __discretize_2d_function(self, func, lower, upper, samples): + lower, upper, sam = self.__determine_2d_domain_bounds(lower, upper, samples) + + # Create nodes to evaluate function + xs = np.linspace(lower[0], upper[0], sam[0]) + ys = np.linspace(lower[1], upper[1], sam[1]) + xs, ys = np.array(np.meshgrid(xs, ys)).reshape(2, xs.size * ys.size) + + # Evaluate function at all mesh nodes and convert it to matrix + zs = np.array(func.get_value(xs, ys)) + func.set_source(np.concatenate(([xs], [ys], [zs])).transpose()) + func.__interpolation__ = "shepard" + func.__extrapolation__ = "natural" + def set_discrete( self, - lower=0, - upper=10, + lower=None, + upper=None, samples=200, interpolation="spline", extrapolation="constant", @@ -647,9 +768,9 @@ def set_discrete( Parameters ---------- lower : scalar, optional - Value where sampling range will start. Default is 0. + Value where sampling range will start. Default is None. upper : scalar, optional - Value where sampling range will end. Default is 10. + Value where sampling range will end. Default is None. samples : int, optional Number of samples to be taken from inside range. Default is 200. interpolation : string @@ -689,24 +810,11 @@ def set_discrete( func = deepcopy(self) if not mutate_self else self if func.__dom_dim__ == 1: - xs = np.linspace(lower, upper, samples) - ys = func.get_value(xs.tolist()) if one_by_one else func.get_value(xs) - func.__interpolation__ = interpolation - func.__extrapolation__ = extrapolation - func.set_source(np.column_stack((xs, ys))) + self.__discretize_1d_function( + func, lower, upper, samples, interpolation, extrapolation, one_by_one + ) elif func.__dom_dim__ == 2: - lower = 2 * [lower] if isinstance(lower, NUMERICAL_TYPES) else lower - upper = 2 * [upper] if isinstance(upper, NUMERICAL_TYPES) else upper - sam = 2 * [samples] if isinstance(samples, NUMERICAL_TYPES) else samples - # Create nodes to evaluate function - xs = np.linspace(lower[0], upper[0], sam[0]) - ys = np.linspace(lower[1], upper[1], sam[1]) - xs, ys = np.array(np.meshgrid(xs, ys)).reshape(2, xs.size * ys.size) - # Evaluate function at all mesh nodes and convert it to matrix - zs = np.array(func.get_value(xs, ys)) - func.set_source(np.concatenate(([xs], [ys], [zs])).transpose()) - func.__interpolation__ = "shepard" - func.__extrapolation__ = "natural" + self.__discretize_2d_function(func, lower, upper, samples) else: raise ValueError( "Discretization is only supported for 1-D and 2-D Functions." @@ -897,6 +1005,244 @@ def reset( return self + def __crop_array_source(self, cropped_func, x_lim): + """Crop the array source of a Function based on domain limits. + + Parameters + ---------- + cropped_func : Function + The Function instance to be cropped. + x_lim : list[tuple] + Range of values with lower and upper limits for cropping. + """ + if cropped_func.__dom_dim__ == 1: + cropped_func.source = cropped_func.source[ + (cropped_func.source[:, 0] >= x_lim[0][0]) + & (cropped_func.source[:, 0] <= x_lim[0][1]) + ] + elif cropped_func.__dom_dim__ == 2: + cropped_func.source = cropped_func.source[ + (cropped_func.source[:, 0] >= x_lim[0][0]) + & (cropped_func.source[:, 0] <= x_lim[0][1]) + & (cropped_func.source[:, 1] >= x_lim[1][0]) + & (cropped_func.source[:, 1] <= x_lim[1][1]) + ] + + def __set_cropped_domain_1d(self, cropped_func, x_lim): + """Set the cropped domain for 1-D functions. + + Parameters + ---------- + cropped_func : Function + The Function instance to set the cropped domain for. + x_lim : list[tuple] + Range of values with lower and upper limits. + """ + if x_lim[0][0] < x_lim[0][1]: + cropped_func.__cropped_domain__ = x_lim[0] + + def __set_cropped_domain_2d(self, cropped_func, x_lim): + """Set the cropped domain for 2-D functions. + + Parameters + ---------- + cropped_func : Function + The Function instance to set the cropped domain for. + x_lim : list[tuple] + Range of values with lower and upper limits. + """ + if len(x_lim) < 2: + raise IndexError("x_lim must have a length of 2 for 2-D function") + + if x_lim[0] is not None and x_lim[0][0] < x_lim[0][1]: + cropped_func.__cropped_domain__ = [x_lim[0]] + else: + cropped_func.__cropped_domain__ = [None] + + if len(x_lim) >= 2 and x_lim[1] is not None and x_lim[1][0] < x_lim[1][1]: + cropped_func.__cropped_domain__.append(x_lim[1]) + else: + cropped_func.__cropped_domain__.append(None) + + def crop(self, x_lim): + """Restrict the **input** domain of the Function to specified ranges. + + This method limits the input values of the Function to the intervals + defined in `x_lim`, effectively trimming the data so that only values + within the specified ranges are retained. For multi-dimensional + functions, each dimension can be cropped independently by providing a + tuple with lower and upper bounds for each input variable. If a + dimension is set to `None`, it will not be cropped. + + Parameters + ---------- + x_lim : list[tuple] + Range of values with lower and upper limits for input values to be + cropped within. + + Returns + ------- + Function + A new Function instance with the cropped domain. + + See also + -------- + Function.clip + + Examples + -------- + >>> from rocketpy import Function + >>> import numpy as np + + Create two 2D functions: + >>> f1 = Function( + ... lambda x1, x2: np.sin(x1)*np.cos(x2), + ... inputs=['x1', 'x2'], + ... outputs='y' + ... ) + >>> f2 = Function( + ... lambda x1, x2: np.cos(x1)*np.sin(x2), + ... inputs=['x1', 'x2'], + ... outputs='y' + ... ) + + Crop their domains: + >>> f1_cropped = f1.crop([(-1, 1), (-2, 2)]) + >>> f2_cropped = f2.crop([None, (-2, 2)]) + + Compare the cropped functions using Function.compare_plots: + >>> # Function.compare_plots([ + >>> # (f1_cropped, 'sin(x1)*cos(x2), cropped'), + >>> # (f2_cropped, 'cos(x1)*sin(x2), cropped') + >>> # ]) + """ + if not isinstance(x_lim, list): + raise TypeError("x_lim must be a list of tuples.") + + if len(x_lim) > self.__dom_dim__: + raise ValueError( + "x_lim must not exceed the length of the domain dimension." + ) + + cropped_func = deepcopy(self) + + if isinstance(cropped_func.source, np.ndarray): + self.__crop_array_source(cropped_func, x_lim) + + if cropped_func.__dom_dim__ == 1: + self.__set_cropped_domain_1d(cropped_func, x_lim) + elif cropped_func.__dom_dim__ == 2: + self.__set_cropped_domain_2d(cropped_func, x_lim) + + cropped_func.set_source(cropped_func.source) + return cropped_func + + def __validate_clip_parameters(self, y_lim): + if not isinstance(y_lim, list): + raise TypeError("y_lim must be a list of tuples.") + + if len(y_lim) != len(self.__outputs__): + raise ValueError( + "y_lim must have the same length as the output dimensions." + ) + + def __clip_array_source(self, clipped_func, y_lim: list[tuple]): + clipped_func.source = clipped_func.source[ + (clipped_func.source[:, clipped_func.__dom_dim__] >= y_lim[0][0]) + & (clipped_func.source[:, clipped_func.__dom_dim__] <= y_lim[0][1]) + ] + + def __clip_numerical_source(self, clipped_func, y_lim: list[tuple]): + try: + if clipped_func.source < y_lim[0][0]: + raise ArithmeticError("Constant function outside range") + if clipped_func.source > y_lim[0][1]: + raise ArithmeticError("Constant function outside range") + except TypeError as e: + raise TypeError("y_lim must be the same type as the function source") from e + + def __clip_callable_source(self, clipped_func, y_lim: list[tuple]): + original_function = clipped_func.source + + def clipped_function(*args): + results = original_function(*args) + clipped_results = [] + + if isinstance(results, (tuple, list)): + # Multi-dimensional output + for i, (lower, upper) in enumerate(y_lim): + clipped_results.append(max(lower, min(upper, results[i]))) + else: + # Single value output + for lower, upper in y_lim: + clipped_results.append(max(lower, min(upper, results))) + + return ( + tuple(clipped_results) + if len(clipped_results) > 1 + else clipped_results[0] + ) + + clipped_func.source = clipped_function + + def clip(self, y_lim): + """Restrict the **output** values of the Function to specified ranges. + + This method limits the output values of the Function to the intervals + defined in `y_lim`, effectively removing all input-output pairs where + the output values fall outside the specified ranges. This operation + filters the data based on output constraints rather than input domain + restrictions. + + Parameters + ---------- + y_lim : list[tuple] + Range of values with lower and upper limits for output values to be + clipped within. + + Returns + ------- + Function + A new Function instance with the clipped output values. + + See also + -------- + Function.crop + + Examples + -------- + >>> from rocketpy import Function + >>> + >>> f = Function(lambda x: x**2, inputs='x', outputs='y') + >>> print(f) + Function from R1 to R1 : (x) → (y) + >>> f_clipped = f.clip([(-5, 5)]) + >>> print(f_clipped) + Function from R1 to R1 : (x) → (y) + """ + self.__validate_clip_parameters(y_lim) + + clipped_func = deepcopy(self) + + if isinstance(clipped_func.source, np.ndarray): + self.__clip_array_source(clipped_func, y_lim) + elif isinstance(clipped_func.source, NUMERICAL_TYPES): + self.__clip_numerical_source(clipped_func, y_lim) + elif callable(clipped_func.source): + self.__clip_callable_source(clipped_func, y_lim) + + try: + clipped_func.set_source(clipped_func.source) + except ValueError as e: + raise ValueError( + "Cannot clip function as function reduces to " + f"{len(clipped_func.source) if isinstance(clipped_func.source, (list, np.ndarray)) else 'unknown'} points (too few data points to define" + " a domain). Ensure that the source is array-like and has " + "sufficient data points after applying the clipping function." + ) from e + + return clipped_func + # Define all get methods def get_inputs(self): "Return tuple of inputs of the function." @@ -1459,14 +1805,13 @@ def plot(self, *args, **kwargs): else: print("Error: Only functions with 1D or 2D domains can be plotted.") + @deprecated( + reason="The `Function.plot1D` method is set to be deprecated and fully " + "removed in rocketpy v2.0.0", + alternative="Function.plot_1d", + ) def plot1D(self, *args, **kwargs): # pragma: no cover """Deprecated method, use Function.plot_1d instead.""" - warnings.warn( - "The `Function.plot1D` method is set to be deprecated and fully " - + "removed in rocketpy v2.0.0, use `Function.plot_1d` instead. " - + "This method is calling `Function.plot_1d`.", - DeprecationWarning, - ) return self.plot_1d(*args, **kwargs) def plot_1d( # pylint: disable=too-many-statements @@ -1525,8 +1870,13 @@ def plot_1d( # pylint: disable=too-many-statements ax = fig.axes if self._source_type is SourceType.CALLABLE: # Determine boundaries - lower = 0 if lower is None else lower - upper = 10 if upper is None else upper + domain = [0, 10] + if self.__cropped_domain__[0] and self.__cropped_domain__[0] > domain[0]: + domain[0] = self.__cropped_domain__[0] + if self.__cropped_domain__[1] and self.__cropped_domain__[1] < domain[1]: + domain[1] = self.__cropped_domain__[1] + lower = domain[0] if lower is None else lower + upper = domain[1] if upper is None else upper else: # Determine boundaries x_data = self.x_array @@ -1559,14 +1909,13 @@ def plot_1d( # pylint: disable=too-many-statements if return_object: return fig, ax + @deprecated( + reason="The `Function.plot2D` method is set to be deprecated and fully " + "removed in rocketpy v2.0.0", + alternative="Function.plot_2d", + ) def plot2D(self, *args, **kwargs): # pragma: no cover """Deprecated method, use Function.plot_2d instead.""" - warnings.warn( - "The `Function.plot2D` method is set to be deprecated and fully " - + "removed in rocketpy v2.0.0, use `Function.plot_2d` instead. " - + "This method is calling `Function.plot_2d`.", - DeprecationWarning, - ) return self.plot_2d(*args, **kwargs) def plot_2d( # pylint: disable=too-many-statements @@ -1637,9 +1986,17 @@ def plot_2d( # pylint: disable=too-many-statements # Define a mesh and f values at mesh nodes for plotting if self._source_type is SourceType.CALLABLE: # Determine boundaries - lower = [0, 0] if lower is None else lower + domain = [[0, 10], [0, 10]] + if self.__cropped_domain__ is not None: + for i in range(0, 2): + if self.__cropped_domain__[i] is not None: + if self.__cropped_domain__[i][0] > domain[i][0]: + domain[i][0] = self.__cropped_domain__[i][0] + if self.__cropped_domain__[i][1] < domain[i][1]: + domain[i][1] = self.__cropped_domain__[i][1] + lower = [domain[0][0], domain[1][0]] if lower is None else lower lower = 2 * [lower] if isinstance(lower, NUMERICAL_TYPES) else lower - upper = [10, 10] if upper is None else upper + upper = [domain[0][1], domain[1][1]] if upper is None else upper upper = 2 * [upper] if isinstance(upper, NUMERICAL_TYPES) else upper else: # Determine boundaries @@ -3348,8 +3705,9 @@ def savetxt( + " must be provided." ) # Generate the data points using the callable - x = np.linspace(lower, upper, samples) - data_points = np.column_stack((x, self(x))) + data_points = self.set_discrete( + lower, upper, samples, mutate_self=False + ).source else: # If the source is already an array, use it as is data_points = self.source @@ -3567,7 +3925,7 @@ def __validate_extrapolation(self, extrapolation): extrapolation = "natural" return extrapolation - def to_dict(self, include_outputs=False): # pylint: disable=unused-argument + def to_dict(self, **kwargs): # pylint: disable=unused-argument """Serializes the Function instance to a dictionary. Returns @@ -3578,7 +3936,10 @@ def to_dict(self, include_outputs=False): # pylint: disable=unused-argument source = self.source if callable(source): - source = to_hex_encode(source) + if kwargs.get("allow_pickle", True): + source = to_hex_encode(source) + else: + source = source.__name__ return { "source": source, diff --git a/rocketpy/mathutils/vector_matrix.py b/rocketpy/mathutils/vector_matrix.py index ceac9a08b..9c8bde144 100644 --- a/rocketpy/mathutils/vector_matrix.py +++ b/rocketpy/mathutils/vector_matrix.py @@ -403,7 +403,7 @@ def zeros(): """Returns the zero vector.""" return Vector([0, 0, 0]) - def to_dict(self, include_outputs=False): # pylint: disable=unused-argument + def to_dict(self, **kwargs): # pylint: disable=unused-argument """Returns the vector as a JSON compatible element.""" return list(self.components) @@ -1007,7 +1007,7 @@ def __repr__(self): + f" [{self.zx}, {self.zy}, {self.zz}])" ) - def to_dict(self, include_outputs=False): # pylint: disable=unused-argument + def to_dict(self, **kwargs): # pylint: disable=unused-argument """Returns the matrix as a JSON compatible element.""" return [list(row) for row in self.components] diff --git a/rocketpy/motors/fluid.py b/rocketpy/motors/fluid.py index e027702d0..33d882abf 100644 --- a/rocketpy/motors/fluid.py +++ b/rocketpy/motors/fluid.py @@ -1,44 +1,147 @@ -from dataclasses import dataclass +import numpy as np +from scipy.constants import atm, zero_Celsius +from ..mathutils.function import NUMERICAL_TYPES, Function from ..plots.fluid_plots import _FluidPlots from ..prints.fluid_prints import _FluidPrints -@dataclass class Fluid: - """Class that represents a fluid. + """Class that represents a fluid object and its attributes, + such as density. Attributes ---------- name : str Name of the fluid. - density : float - Density of the fluid in kg/m³. + density : int, float, callable, string, array, Function + Input density of the fluid in kg/m³. + Used internally to compute the density_function. + density_function : Function + Density of the fluid as a function of temperature and pressure. """ - name: str - density: float - - def __post_init__(self): - """Post initialization method. - - Raises - ------ - ValueError - If the name is not a string. - ValueError - If the density is not a positive number. + def __init__(self, name, density): + """Initializes a Fluid object. + + Parameters + ---------- + name : str + Name of the fluid. + density : int, float, callable, string, array, Function + Density of the fluid in kg/m³ as a function of temperature + and pressure. If a int or float is given, it is considered + a constant function. A callable, csv file or an list of points + can be given to express the density values. + The parameter is used as a ``Function`` source. + + See more on :class:`rocketpy.Function` source types. """ - if not isinstance(self.name, str): # pragma: no cover - raise ValueError("The name must be a string.") - if self.density < 0: # pragma: no cover - raise ValueError("The density must be a positive number.") + self.name = name + self.density = density + self.density_function = density - # Initialize plots and prints object self.prints = _FluidPrints(self) self.plots = _FluidPlots(self) + @property + def density(self): + """Density of the fluid from the class parameter input.""" + return self._density + + @density.setter + def density(self, value): + """Setter of the density class parameter.""" + if isinstance(value, NUMERICAL_TYPES): + # Numerical value kept for retro-compatibility + self._density = value + else: + self._density = Function( + value, + interpolation="linear", + extrapolation="natural", + inputs=["Temperature (K)", "Pressure (Pa)"], + outputs="Density (kg/m³)", + ) + + @property + def density_function(self): + """Density of the fluid as a function of temperature and pressure.""" + return self._density_function + + @density_function.setter + def density_function(self, value): + """Setter for the density function of the fluid. Numeric values + are converted to constant functions. + + Parameters + ---------- + value : int, float, callable, string, array, Function + Density of the fluid in kg/m³ as a function of temperature + and pressure. The value is used as a ``Function`` source. + """ + if isinstance(value, NUMERICAL_TYPES): + # pylint: disable=unused-argument + def density_function(temperature, pressure): + return value + + else: + density_function = value + + self._density_function = Function( + density_function, + interpolation="linear", + extrapolation="natural", + inputs=["Temperature (K)", "Pressure (Pa)"], + outputs="Density (kg/m³)", + ) + + def get_time_variable_density(self, temperature, pressure): + """Get the density of the fluid as a function of time. + + Parameters + ---------- + temperature : Function + Temperature of the fluid in Kelvin. + pressure : Function + Pressure of the fluid in Pascals. + + Returns + ------- + Function + Density of the fluid in kg/m³ as function of time. + """ + is_temperature_callable = callable(temperature.source) + is_pressure_callable = callable(pressure.source) + if is_temperature_callable and is_pressure_callable: + return Function( + lambda time: self.density_function.get_value( + temperature.source(time), pressure.source(time) + ), + inputs="Time (s)", + outputs="Density (kg/m³)", + ) + + if is_temperature_callable or is_pressure_callable: + time_scale = ( + temperature.x_array if not is_temperature_callable else pressure.x_array + ) + else: + time_scale = np.unique( + np.concatenate((temperature.x_array, pressure.x_array)) + ) + density_time = self.density_function.get_value( + temperature.get_value(time_scale), pressure.get_value(time_scale) + ) + return Function( + np.column_stack((time_scale, density_time)), + interpolation="linear", + extrapolation="constant", + inputs="Time (s)", + outputs="Density (kg/m³)", + ) + def __repr__(self): """Representation method. @@ -61,8 +164,19 @@ def __str__(self): return f"Fluid: {self.name}" - def to_dict(self, include_outputs=False): # pylint: disable=unused-argument - return {"name": self.name, "density": self.density} + def to_dict(self, **kwargs): # pylint: disable=unused-argument + discretize = kwargs.get("discretize", False) + + density = self.density + if discretize and isinstance(density, Function): + density = density.set_discrete( + lower=(100, atm * 0.9), + upper=(zero_Celsius + 30, atm * 50), + samples=(25, 25), + mutate_self=False, + ) + + return {"name": self.name, "density": density} @classmethod def from_dict(cls, data): diff --git a/rocketpy/motors/hybrid_motor.py b/rocketpy/motors/hybrid_motor.py index e227e4550..7cb28670c 100644 --- a/rocketpy/motors/hybrid_motor.py +++ b/rocketpy/motors/hybrid_motor.py @@ -641,8 +641,8 @@ def draw(self, *, filename=None): """ self.plots.draw(filename=filename) - def to_dict(self, include_outputs=False): - data = super().to_dict(include_outputs) + def to_dict(self, **kwargs): + data = super().to_dict(**kwargs) data.update( { "grain_number": self.grain_number, @@ -660,13 +660,18 @@ def to_dict(self, include_outputs=False): } ) - if include_outputs: + if kwargs.get("include_outputs", False): + burn_rate = self.solid.burn_rate + if kwargs.get("discretize", False): + burn_rate = burn_rate.set_discrete_based_on_model( + self.thrust, mutate_self=False + ) data.update( { "grain_inner_radius": self.solid.grain_inner_radius, "grain_height": self.solid.grain_height, "burn_area": self.solid.burn_area, - "burn_rate": self.solid.burn_rate, + "burn_rate": burn_rate, } ) diff --git a/rocketpy/motors/liquid_motor.py b/rocketpy/motors/liquid_motor.py index acce874af..1cd76a52f 100644 --- a/rocketpy/motors/liquid_motor.py +++ b/rocketpy/motors/liquid_motor.py @@ -497,8 +497,8 @@ def draw(self, *, filename=None): """ self.plots.draw(filename=filename) - def to_dict(self, include_outputs=False): - data = super().to_dict(include_outputs) + def to_dict(self, **kwargs): + data = super().to_dict(**kwargs) data.update( { "positioned_tanks": [ diff --git a/rocketpy/motors/motor.py b/rocketpy/motors/motor.py index 7178cdcf5..7930ed52b 100644 --- a/rocketpy/motors/motor.py +++ b/rocketpy/motors/motor.py @@ -1151,7 +1151,7 @@ def vacuum_thrust(self): Returns ------- vacuum_thrust : Function - The rocket's thrust in a vaccum. + The rocket's thrust in a vacuum. """ if self.reference_pressure is None: warnings.warn( @@ -1231,14 +1231,7 @@ def get_attr_value(obj, attr_name, multiplier=1): # Write last line file.write(f"{self.thrust.source[-1, 0]:.4f} {0:.3f}\n") - def to_dict(self, include_outputs=False): - thrust_source = self.thrust_source - - if isinstance(thrust_source, str): - thrust_source = self.thrust.source - elif callable(thrust_source) and not isinstance(thrust_source, Function): - thrust_source = Function(thrust_source) - + def to_dict(self, **kwargs): data = { "thrust_source": self.thrust, "dry_I_11": self.dry_I_11, @@ -1258,31 +1251,94 @@ def to_dict(self, include_outputs=False): "reference_pressure": self.reference_pressure, } - if include_outputs: + if kwargs.get("include_outputs", False): + total_mass = self.total_mass + propellant_mass = self.propellant_mass + mass_flow_rate = self.total_mass_flow_rate + center_of_mass = self.center_of_mass + center_of_propellant_mass = self.center_of_propellant_mass + exhaust_velocity = self.exhaust_velocity + I_11 = self.I_11 + I_22 = self.I_22 + I_33 = self.I_33 + I_12 = self.I_12 + I_13 = self.I_13 + I_23 = self.I_23 + propellant_I_11 = self.propellant_I_11 + propellant_I_22 = self.propellant_I_22 + propellant_I_33 = self.propellant_I_33 + propellant_I_12 = self.propellant_I_12 + propellant_I_13 = self.propellant_I_13 + propellant_I_23 = self.propellant_I_23 + if kwargs.get("discretize", False): + total_mass = total_mass.set_discrete_based_on_model( + self.thrust, mutate_self=False + ) + propellant_mass = propellant_mass.set_discrete_based_on_model( + self.thrust, mutate_self=False + ) + mass_flow_rate = mass_flow_rate.set_discrete_based_on_model( + self.thrust, mutate_self=False + ) + center_of_mass = center_of_mass.set_discrete_based_on_model( + self.thrust, mutate_self=False + ) + center_of_propellant_mass = ( + center_of_propellant_mass.set_discrete_based_on_model( + self.thrust, mutate_self=False + ) + ) + exhaust_velocity = exhaust_velocity.set_discrete_based_on_model( + self.thrust, mutate_self=False + ) + I_11 = I_11.set_discrete_based_on_model(self.thrust, mutate_self=False) + I_22 = I_22.set_discrete_based_on_model(self.thrust, mutate_self=False) + I_33 = I_33.set_discrete_based_on_model(self.thrust, mutate_self=False) + I_12 = I_12.set_discrete_based_on_model(self.thrust, mutate_self=False) + I_13 = I_13.set_discrete_based_on_model(self.thrust, mutate_self=False) + I_23 = I_23.set_discrete_based_on_model(self.thrust, mutate_self=False) + propellant_I_11 = propellant_I_11.set_discrete_based_on_model( + self.thrust, mutate_self=False + ) + propellant_I_22 = propellant_I_22.set_discrete_based_on_model( + self.thrust, mutate_self=False + ) + propellant_I_33 = propellant_I_33.set_discrete_based_on_model( + self.thrust, mutate_self=False + ) + propellant_I_12 = propellant_I_12.set_discrete_based_on_model( + self.thrust, mutate_self=False + ) + propellant_I_13 = propellant_I_13.set_discrete_based_on_model( + self.thrust, mutate_self=False + ) + propellant_I_23 = propellant_I_23.set_discrete_based_on_model( + self.thrust, mutate_self=False + ) data.update( { "vacuum_thrust": self.vacuum_thrust, - "total_mass": self.total_mass, - "propellant_mass": self.propellant_mass, - "mass_flow_rate": self.mass_flow_rate, - "center_of_mass": self.center_of_mass, - "center_of_propellant_mass": self.center_of_propellant_mass, + "total_mass": total_mass, + "propellant_mass": propellant_mass, + "mass_flow_rate": mass_flow_rate, + "center_of_mass": center_of_mass, + "center_of_propellant_mass": center_of_propellant_mass, "total_impulse": self.total_impulse, - "exhaust_velocity": self.exhaust_velocity, + "exhaust_velocity": exhaust_velocity, "propellant_initial_mass": self.propellant_initial_mass, "structural_mass_ratio": self.structural_mass_ratio, - "I_11": self.I_11, - "I_22": self.I_22, - "I_33": self.I_33, - "I_12": self.I_12, - "I_13": self.I_13, - "I_23": self.I_23, - "propellant_I_11": self.propellant_I_11, - "propellant_I_22": self.propellant_I_22, - "propellant_I_33": self.propellant_I_33, - "propellant_I_12": self.propellant_I_12, - "propellant_I_13": self.propellant_I_13, - "propellant_I_23": self.propellant_I_23, + "I_11": I_11, + "I_22": I_22, + "I_33": I_33, + "I_12": I_12, + "I_13": I_13, + "I_23": I_23, + "propellant_I_11": propellant_I_11, + "propellant_I_22": propellant_I_22, + "propellant_I_33": propellant_I_33, + "propellant_I_12": propellant_I_12, + "propellant_I_13": propellant_I_13, + "propellant_I_23": propellant_I_23, } ) @@ -1510,7 +1566,7 @@ def center_of_propellant_mass(self): Function Function representing the center of mass of the motor. """ - return self.chamber_position + return Function(self.chamber_position).set_discrete_based_on_model(self.thrust) @funcify_method("Time (s)", "Inertia I_11 (kg m²)") def propellant_I_11(self): @@ -1532,11 +1588,11 @@ def propellant_I_11(self): ---------- https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ - return ( + return Function( self.propellant_mass * (3 * self.chamber_radius**2 + self.chamber_height**2) / 12 - ) + ).set_discrete_based_on_model(self.thrust) @funcify_method("Time (s)", "Inertia I_22 (kg m²)") def propellant_I_22(self): @@ -1580,19 +1636,21 @@ def propellant_I_33(self): ---------- https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ - return self.propellant_mass * self.chamber_radius**2 / 2 + return Function( + self.propellant_mass * self.chamber_radius**2 / 2 + ).set_discrete_based_on_model(self.thrust) @funcify_method("Time (s)", "Inertia I_12 (kg m²)") def propellant_I_12(self): - return Function(0) + return Function(0).set_discrete_based_on_model(self.thrust) @funcify_method("Time (s)", "Inertia I_13 (kg m²)") def propellant_I_13(self): - return Function(0) + return Function(0).set_discrete_based_on_model(self.thrust) @funcify_method("Time (s)", "Inertia I_23 (kg m²)") def propellant_I_23(self): - return Function(0) + return Function(0).set_discrete_based_on_model(self.thrust) @staticmethod def load_from_eng_file( @@ -1862,8 +1920,8 @@ def all_info(self): self.prints.all() self.plots.all() - def to_dict(self, include_outputs=False): - data = super().to_dict(include_outputs) + def to_dict(self, **kwargs): + data = super().to_dict(**kwargs) data.update( { "chamber_radius": self.chamber_radius, diff --git a/rocketpy/motors/solid_motor.py b/rocketpy/motors/solid_motor.py index a8d823966..8a00eeec9 100644 --- a/rocketpy/motors/solid_motor.py +++ b/rocketpy/motors/solid_motor.py @@ -765,8 +765,8 @@ def draw(self, *, filename=None): """ self.plots.draw(filename=filename) - def to_dict(self, include_outputs=False): - data = super().to_dict(include_outputs) + def to_dict(self, **kwargs): + data = super().to_dict(**kwargs) data.update( { "nozzle_radius": self.nozzle_radius, @@ -781,13 +781,18 @@ def to_dict(self, include_outputs=False): } ) - if include_outputs: + if kwargs.get("include_outputs", False): + burn_rate = self.burn_rate + if kwargs.get("discretize", False): + burn_rate = burn_rate.set_discrete_based_on_model( + self.thrust, mutate_self=False + ) data.update( { "grain_inner_radius": self.grain_inner_radius, "grain_height": self.grain_height, "burn_area": self.burn_area, - "burn_rate": self.burn_rate, + "burn_rate": burn_rate, "Kn": self.Kn, } ) diff --git a/rocketpy/motors/tank.py b/rocketpy/motors/tank.py index 18060f1a5..6f3496b32 100644 --- a/rocketpy/motors/tank.py +++ b/rocketpy/motors/tank.py @@ -1,11 +1,12 @@ from abc import ABC, abstractmethod import numpy as np +from scipy.constants import atm, zero_Celsius from ..mathutils.function import Function, funcify_method from ..plots.tank_plots import _TankPlots from ..prints.tank_prints import _TankPrints -from ..tools import tuple_handler +from ..tools import deprecated, tuple_handler class Tank(ABC): @@ -75,7 +76,17 @@ class Tank(ABC): Tank symmetry axis. The reference point is the Tank center of mass. """ - def __init__(self, name, geometry, flux_time, liquid, gas, discretize=100): + def __init__( + self, + name, + geometry, + flux_time, + liquid, + gas, + discretize=100, + temperature=None, + pressure=None, + ): """Initialize Tank class. Parameters @@ -98,8 +109,24 @@ def __init__(self, name, geometry, flux_time, liquid, gas, discretize=100): Liquid inside the tank as a Fluid object. discretize : int, optional Number of points to discretize fluid inputs. If the input - already has a appropriate discretization, this parameter + already has a appropriate uniform discretization, this parameter must be set to None. The default is 100. + temperature : int, float, callable, string, array, Function + Temperature inside the tank as a function of time in K. If a callable + is given, it must be a function of time in seconds. An array of points + can also be given as a list or a string with the path to a .csv file + with two columns, the first one being the time in seconds and the + second one being the temperature in K. The default is None. This + parameter is only required if fluid (``liquid`` or ``gas`` parameters) + densities are functions of temperature. + pressure : int, float, callable, string, array, Function + Pressure inside the tank as a function of time in Pa. If a callable + is given, it must be a function of time in seconds. An array of points + can also be given as a list or a string with the path to a .csv file + with two columns, the first one being the time in seconds and the + second one being the pressure in Pa. The default is None. This + parameter is only required if fluid (``liquid`` or ``gas`` parameters) + densities are functions of pressure. """ self.name = name self.geometry = geometry @@ -107,6 +134,15 @@ def __init__(self, name, geometry, flux_time, liquid, gas, discretize=100): self.gas = gas self.liquid = liquid self.discretize = discretize + self.temperature = temperature + self.pressure = pressure + + self._liquid_density = self.liquid.get_time_variable_density( + self.temperature, self.pressure + ) + self._gas_density = self.gas.get_time_variable_density( + self.temperature, self.pressure + ) # Initialize plots and prints object self.prints = _TankPrints(self) @@ -134,6 +170,80 @@ def flux_time(self, flux_time): """ self._flux_time = tuple_handler(flux_time) + @property + def temperature(self): + """Returns the temperature of the tank as a function of time. + + Returns + ------- + Function + Temperature of the tank as a function of time. + """ + return self._temperature + + @temperature.setter + def temperature(self, temperature): + """Sets the temperature of the tank as a function of time. + + Parameters + ---------- + temperature : int, float, callable, string, array, Function + Temperature inside the tank as a function of time in K as + a ``Function`` source. + """ + if temperature is None: + temperature = zero_Celsius + + _temperature = Function( + temperature, + interpolation="linear", + extrapolation="constant", + inputs="Time (s)", + outputs="Temperature (K)", + ) + + if self.discretize: + _temperature = _temperature.set_discrete(*self.flux_time, self.discretize) + + self._temperature = _temperature + + @property + def pressure(self): + """Returns the pressure of the tank as a function of time. + + Returns + ------- + Function + Pressure of the tank as a function of time. + """ + return self._pressure + + @pressure.setter + def pressure(self, pressure): + """Sets the pressure of the tank as a function of time. + + Parameters + ---------- + pressure : int, float, callable, string, array, Function + Pressure inside the tank as a function of time in Pa as + a ``Function`` source. + """ + if pressure is None: + pressure = atm + + _pressure = Function( + pressure, + interpolation="linear", + extrapolation="constant", + inputs="Time (s)", + outputs="Pressure (Pa)", + ) + + if self.discretize: + _pressure = _pressure.set_discrete(*self.flux_time, self.discretize) + + self._pressure = _pressure + @property @abstractmethod def fluid_mass(self): @@ -326,7 +436,7 @@ def center_of_mass(self): # Check for zero mass bound_mass = ( - self.fluid_mass < 0.001 * self.geometry.total_volume * self.gas.density + self.fluid_mass < 0.001 * self.geometry.total_volume * self._gas_density ) if bound_mass.any(): # TODO: pending Function setter impl. @@ -360,7 +470,7 @@ def liquid_inertia(self): self.liquid_volume * (self.liquid_center_of_mass - self.center_of_mass) ** 2 ) - return self.liquid.density * Ix_volume + return self._liquid_density * Ix_volume @funcify_method("Time (s)", "Gas Inertia (kg*m^2)") def gas_inertia(self): @@ -386,7 +496,7 @@ def gas_inertia(self): self.gas_volume * (self.gas_center_of_mass - self.center_of_mass) ** 2 ) - return self.gas.density * inertia_volume + return self._gas_density * inertia_volume @funcify_method("Time (s)", "Fluid Inertia (kg*m^2)") def inertia(self): @@ -477,6 +587,10 @@ def underfill_height_exception(param_name, param): elif (height < bottom_tolerance).any(): underfill_height_exception(name, height) + @abstractmethod + def _discretize_fluid_inputs(self): + """Uniformly discretizes the parameter of inputs of fluid data .""" + def draw(self, *, filename=None): """Draws the tank geometry. @@ -505,7 +619,7 @@ def all_info(self): self.prints.all() self.plots.all() - def to_dict(self, include_outputs=False): + def to_dict(self, **kwargs): data = { "name": self.name, "geometry": self.geometry, @@ -513,8 +627,10 @@ def to_dict(self, include_outputs=False): "liquid": self.liquid, "gas": self.gas, "discretize": self.discretize, + "temperature": self.temperature, + "pressure": self.pressure, } - if include_outputs: + if kwargs.get("include_outputs", False): data.update( { "fluid_mass": self.fluid_mass, @@ -560,6 +676,8 @@ def __init__( liquid_mass_flow_rate_out, gas_mass_flow_rate_out, discretize=100, + temperature=None, + pressure=None, ): """Initializes the MassFlowRateBasedTank class. @@ -619,8 +737,26 @@ def __init__( this parameter may be set to None. Otherwise, an uniform discretization will be applied based on the discretize value. The default is 100. - """ - super().__init__(name, geometry, flux_time, liquid, gas, discretize) + temperature : int, float, callable, string, array, Function + Temperature inside the tank as a function of time in K. If a callable + is given, it must be a function of time in seconds. An array of points + can also be given as a list or a string with the path to a .csv file + with two columns, the first one being the time in seconds and the + second one being the temperature in K. The default is None. This + parameter is only required if fluid (``liquid`` or ``gas`` parameters) + densities are functions of temperature. + pressure : int, float, callable, string, array, Function + Pressure inside the tank as a function of time in Pa. If a callable + is given, it must be a function of time in seconds. An array of points + can also be given as a list or a string with the path to a .csv file + with two columns, the first one being the time in seconds and the + second one being the pressure in Pa. The default is None. This + parameter is only required if fluid (``liquid`` or ``gas`` parameters) + densities are functions of pressure. + """ + super().__init__( + name, geometry, flux_time, liquid, gas, discretize, temperature, pressure + ) self.initial_liquid_mass = initial_liquid_mass self.initial_gas_mass = initial_gas_mass @@ -654,9 +790,7 @@ def __init__( extrapolation="zero", ) - # Discretize input flow if needed - if discretize: - self.discretize_flow() + self._discretize_fluid_inputs() # Check if the tank is overfilled or underfilled self._check_volume_bounds() @@ -792,7 +926,7 @@ def liquid_volume(self): Function Volume of the liquid as a function of time. """ - return self.liquid_mass / self.liquid.density + return self.liquid_mass / self._liquid_density @funcify_method("Time (s)", "Gas Volume (m³)") def gas_volume(self): @@ -804,7 +938,7 @@ def gas_volume(self): Function Volume of the gas as a function of time. """ - return self.gas_mass / self.gas.density + return self.gas_mass / self._gas_density @funcify_method("Time (s)", "Liquid Height (m)") def liquid_height(self): @@ -869,25 +1003,41 @@ def gas_height(self): ) return gas_height + def _discretize_fluid_inputs(self): + """Uniformly discretizes the parameter of inputs of fluid data .""" + if self.discretize: + self.liquid_mass_flow_rate_in.set_discrete( + *self.flux_time, self.discretize, "linear" + ) + self.gas_mass_flow_rate_in.set_discrete( + *self.flux_time, self.discretize, "linear" + ) + self.liquid_mass_flow_rate_out.set_discrete( + *self.flux_time, self.discretize, "linear" + ) + self.gas_mass_flow_rate_out.set_discrete( + *self.flux_time, self.discretize, "linear" + ) + else: + # Discretize densities for backward compatibility + self._liquid_density.set_discrete_based_on_model( + self.liquid_mass_flow_rate_in + ) + self._gas_density.set_discrete_based_on_model(self.gas_mass_flow_rate_in) + + @deprecated( + "Should not be a public member of the class.", + "1.12.0", + "_discretize_fluid_inputs", + ) def discretize_flow(self): """Discretizes the mass flow rate inputs according to the flux time and the discretize parameter. """ - self.liquid_mass_flow_rate_in.set_discrete( - *self.flux_time, self.discretize, "linear" - ) - self.gas_mass_flow_rate_in.set_discrete( - *self.flux_time, self.discretize, "linear" - ) - self.liquid_mass_flow_rate_out.set_discrete( - *self.flux_time, self.discretize, "linear" - ) - self.gas_mass_flow_rate_out.set_discrete( - *self.flux_time, self.discretize, "linear" - ) + self._discretize_fluid_inputs() - def to_dict(self, include_outputs=False): - data = super().to_dict(include_outputs) + def to_dict(self, **kwargs): + data = super().to_dict(**kwargs) data.update( { "initial_liquid_mass": self.initial_liquid_mass, @@ -915,6 +1065,8 @@ def from_dict(cls, data): liquid_mass_flow_rate_out=data["liquid_mass_flow_rate_out"], gas_mass_flow_rate_out=data["gas_mass_flow_rate_out"], discretize=data["discretize"], + temperature=data.get("temperature"), + pressure=data.get("pressure"), ) @@ -939,6 +1091,8 @@ def __init__( gas, ullage, discretize=100, + temperature=None, + pressure=None, ): """ Parameters @@ -972,15 +1126,32 @@ def __init__( an uniform discretization will be applied based on the discretize value. The default is 100. - """ - super().__init__(name, geometry, flux_time, liquid, gas, discretize) + temperature : int, float, callable, string, array, Function + Temperature inside the tank as a function of time in K. If a callable + is given, it must be a function of time in seconds. An array of points + can also be given as a list or a string with the path to a .csv file + with two columns, the first one being the time in seconds and the + second one being the temperature in K. The default is None. This + parameter is only required if fluid (``liquid`` or ``gas`` parameters) + densities are functions of temperature. + pressure : int, float, callable, string, array, Function + Pressure inside the tank as a function of time in Pa. If a callable + is given, it must be a function of time in seconds. An array of points + can also be given as a list or a string with the path to a .csv file + with two columns, the first one being the time in seconds and the + second one being the pressure in Pa. The default is None. This + parameter is only required if fluid (``liquid`` or ``gas`` parameters) + densities are functions of pressure. + """ + super().__init__( + name, geometry, flux_time, liquid, gas, discretize, temperature, pressure + ) # Define ullage self.ullage = Function(ullage, "Time (s)", "Volume (m³)", "linear") # Discretize input if needed - if discretize: - self.discretize_ullage() + self._discretize_fluid_inputs() # Check if the tank is overfilled or underfilled self._check_volume_bounds() @@ -1065,7 +1236,7 @@ def gas_mass(self): Function Mass of the gas as a function of time. """ - return self.gas_volume * self.gas.density + return self.gas_volume * self._gas_density @funcify_method("Time (s)", "Liquid Mass (kg)") def liquid_mass(self): @@ -1077,7 +1248,7 @@ def liquid_mass(self): Function Mass of the liquid as a function of time. """ - return self.liquid_volume * self.liquid.density + return self.liquid_volume * self._liquid_density @funcify_method("Time (s)", "Liquid Height (m)") def liquid_height(self): @@ -1108,13 +1279,28 @@ def gas_height(self): """ return Function(self.geometry.top).set_discrete_based_on_model(self.gas_volume) + def _discretize_fluid_inputs(self): + """Uniformly discretizes the parameter of inputs of fluid data .""" + if self.discretize: + self.ullage.set_discrete(*self.flux_time, self.discretize, "linear") + else: + # Discretize densities for backward compatibility + self._liquid_density.set_discrete_based_on_model(self.ullage) + self._gas_density.set_discrete_based_on_model(self.ullage) + + @deprecated( + "Should not be a public member of the class.", + "1.12.0", + "_discretize_fluid_inputs", + ) def discretize_ullage(self): - """Discretizes the ullage input according to the flux time and the - discretize parameter.""" - self.ullage.set_discrete(*self.flux_time, self.discretize, "linear") + """Discretizes the ullage input according to the flux time + and the discretize parameter. + """ + self._discretize_fluid_inputs() - def to_dict(self, include_outputs=False): - data = super().to_dict(include_outputs) + def to_dict(self, **kwargs): + data = super().to_dict(**kwargs) data.update({"ullage": self.ullage}) return data @@ -1128,6 +1314,8 @@ def from_dict(cls, data): gas=data["gas"], ullage=data["ullage"], discretize=data["discretize"], + temperature=data.get("temperature"), + pressure=data.get("pressure"), ) @@ -1152,6 +1340,8 @@ def __init__( gas, liquid_height, discretize=100, + temperature=None, + pressure=None, ): """ Parameters @@ -1185,14 +1375,31 @@ def __init__( Otherwise, an uniform discretization will be applied based on the discretize value. The default is 100. - """ - super().__init__(name, geometry, flux_time, liquid, gas, discretize) + temperature : int, float, callable, string, array, Function + Temperature inside the tank as a function of time in K. If a callable + is given, it must be a function of time in seconds. An array of points + can also be given as a list or a string with the path to a .csv file + with two columns, the first one being the time in seconds and the + second one being the temperature in K. The default is None. This + parameter is only required if fluid (``liquid`` or ``gas`` parameters) + densities are functions of temperature. + pressure : int, float, callable, string, array, Function + Pressure inside the tank as a function of time in Pa. If a callable + is given, it must be a function of time in seconds. An array of points + can also be given as a list or a string with the path to a .csv file + with two columns, the first one being the time in seconds and the + second one being the pressure in Pa. The default is None. This + parameter is only required if fluid (``liquid`` or ``gas`` parameters) + densities are functions of pressure. + """ + super().__init__( + name, geometry, flux_time, liquid, gas, discretize, temperature, pressure + ) # Define liquid level function self.liquid_level = Function(liquid_height, "Time (s)", "height (m)", "linear") - if discretize: - self.discretize_liquid_height() + self._discretize_fluid_inputs() # Check if the tank is overfilled or underfilled self._check_height_bounds() @@ -1303,7 +1510,7 @@ def gas_mass(self): Function Mass of the gas as a function of time. """ - return self.gas_volume * self.gas.density + return self.gas_volume * self._gas_density @funcify_method("Time (s)", "Liquid Mass (kg)") def liquid_mass(self): @@ -1315,7 +1522,7 @@ def liquid_mass(self): Function Mass of the liquid as a function of time. """ - return self.liquid_volume * self.liquid.density + return self.liquid_volume * self._liquid_density @funcify_method("Time (s)", "Gas Height (m)", "linear") def gas_height(self): @@ -1335,14 +1542,28 @@ def gas_height(self): self.liquid_level ) + @deprecated( + "Should not be a public member of the class.", + "1.12.0", + "_discretize_fluid_inputs", + ) def discretize_liquid_height(self): """Discretizes the liquid height input according to the flux time and the discretize parameter. """ - self.liquid_level.set_discrete(*self.flux_time, self.discretize, "linear") + self._discretize_fluid_inputs() - def to_dict(self, include_outputs=False): - data = super().to_dict(include_outputs) + def _discretize_fluid_inputs(self): + """Uniformly discretizes the parameter of inputs of fluid data .""" + if self.discretize: + self.liquid_level.set_discrete(*self.flux_time, self.discretize, "linear") + else: + # Discretize densities for backward compatibility + self._liquid_density.set_discrete_based_on_model(self.liquid_level) + self._gas_density.set_discrete_based_on_model(self.liquid_level) + + def to_dict(self, **kwargs): + data = super().to_dict(**kwargs) data.update({"liquid_height": self.liquid_level}) return data @@ -1356,6 +1577,8 @@ def from_dict(cls, data): gas=data["gas"], liquid_height=data["liquid_height"], discretize=data["discretize"], + temperature=data.get("temperature"), + pressure=data.get("pressure"), ) @@ -1379,6 +1602,8 @@ def __init__( liquid_mass, gas_mass, discretize=100, + temperature=None, + pressure=None, ): """ Parameters @@ -1417,15 +1642,32 @@ def __init__( may be set to None. Otherwise, an uniform discretization will be applied based on the discretize value. The default is 100. - """ - super().__init__(name, geometry, flux_time, liquid, gas, discretize) + temperature : int, float, callable, string, array, Function + Temperature inside the tank as a function of time in K. If a callable + is given, it must be a function of time in seconds. An array of points + can also be given as a list or a string with the path to a .csv file + with two columns, the first one being the time in seconds and the + second one being the temperature in K. The default is None. This + parameter is only required if fluid (``liquid`` or ``gas`` parameters) + densities are functions of temperature. + pressure : int, float, callable, string, array, Function + Pressure inside the tank as a function of time in Pa. If a callable + is given, it must be a function of time in seconds. An array of points + can also be given as a list or a string with the path to a .csv file + with two columns, the first one being the time in seconds and the + second one being the pressure in Pa. The default is None. This + parameter is only required if fluid (``liquid`` or ``gas`` parameters) + densities are functions of pressure. + """ + super().__init__( + name, geometry, flux_time, liquid, gas, discretize, temperature, pressure + ) # Define fluid masses self.liquid_mass = Function(liquid_mass, "Time (s)", "Mass (kg)", "linear") self.gas_mass = Function(gas_mass, "Time (s)", "Mass (kg)", "linear") - if discretize: - self.discretize_masses() + self._discretize_fluid_inputs() # Check if the tank is overfilled or underfilled self._check_volume_bounds() @@ -1519,7 +1761,7 @@ def gas_volume(self): Function Volume of the gas as a function of time. """ - return self.gas_mass / self.gas.density + return self.gas_mass / self._gas_density @funcify_method("Time (s)", "Liquid Volume (m³)") def liquid_volume(self): @@ -1531,7 +1773,7 @@ def liquid_volume(self): Function Volume of the liquid as a function of time. """ - return self.liquid_mass / self.liquid.density + return self.liquid_mass / self._liquid_density @funcify_method("Time (s)", "Liquid Height (m)") def liquid_height(self): @@ -1593,15 +1835,29 @@ def gas_height(self): ) return gas_height + @deprecated( + "Should not be a public member of the class.", + "1.12.0", + "_discretize_fluid_inputs", + ) def discretize_masses(self): """Discretizes the fluid mass inputs according to the flux time and the discretize parameter. """ - self.liquid_mass.set_discrete(*self.flux_time, self.discretize, "linear") - self.gas_mass.set_discrete(*self.flux_time, self.discretize, "linear") + self._discretize_fluid_inputs() + + def _discretize_fluid_inputs(self): + """Uniformly discretizes the parameter of inputs of fluid data .""" + if self.discretize: + self.liquid_mass.set_discrete(*self.flux_time, self.discretize, "linear") + self.gas_mass.set_discrete(*self.flux_time, self.discretize, "linear") + else: + # Discretize densities for backward compatibility + self._liquid_density.set_discrete_based_on_model(self.liquid_mass) + self._gas_density.set_discrete_based_on_model(self.gas_mass) - def to_dict(self, include_outputs=False): - data = super().to_dict(include_outputs) + def to_dict(self, **kwargs): + data = super().to_dict(**kwargs) data.update( { "liquid_mass": self.liquid_mass, @@ -1621,4 +1877,6 @@ def from_dict(cls, data): liquid_mass=data["liquid_mass"], gas_mass=data["gas_mass"], discretize=data["discretize"], + temperature=data.get("temperature"), + pressure=data.get("pressure"), ) diff --git a/rocketpy/motors/tank_geometry.py b/rocketpy/motors/tank_geometry.py index 59644a537..485f57b09 100644 --- a/rocketpy/motors/tank_geometry.py +++ b/rocketpy/motors/tank_geometry.py @@ -346,14 +346,17 @@ def add_geometry(self, domain, radius_function): self._geometry[domain] = Function(radius_function) self.radius = PiecewiseFunction(self._geometry, "Height (m)", "radius (m)") - def to_dict(self, include_outputs=False): + def to_dict(self, **kwargs): data = { "geometry": { - str(domain): function for domain, function in self._geometry.items() + str(domain): function.set_discrete(*domain, 50, mutate_self=False) + if kwargs.get("discretize", False) + else function + for domain, function in self._geometry.items() } } - if include_outputs: + if kwargs.get("include_outputs", False): data["outputs"] = { "average_radius": self.average_radius, "bottom": self.bottom, @@ -442,15 +445,15 @@ def upper_cap_radius(h): else: raise ValueError("Tank already has caps.") - def to_dict(self, include_outputs=False): + def to_dict(self, **kwargs): data = { "radius": self.__input_radius, "height": self.height, "spherical_caps": self.has_caps, } - if include_outputs: - data.update(super().to_dict(include_outputs)) + if kwargs.get("include_outputs", False): + data.update(super().to_dict(**kwargs)) return data @@ -482,11 +485,11 @@ def __init__(self, radius, geometry_dict=None): self.__input_radius = radius self.add_geometry((-radius, radius), lambda h: (radius**2 - h**2) ** 0.5) - def to_dict(self, include_outputs=False): + def to_dict(self, **kwargs): data = {"radius": self.__input_radius} - if include_outputs: - data.update(super().to_dict(include_outputs)) + if kwargs.get("include_outputs", False): + data.update(super().to_dict(**kwargs)) return data diff --git a/rocketpy/plots/fluid_plots.py b/rocketpy/plots/fluid_plots.py index b0292a8c6..8528ce4f8 100644 --- a/rocketpy/plots/fluid_plots.py +++ b/rocketpy/plots/fluid_plots.py @@ -1,3 +1,8 @@ +import warnings + +from scipy.constants import atm, zero_Celsius + + class _FluidPlots: """Class that holds plot methods for Fluid class. @@ -23,6 +28,28 @@ def __init__(self, fluid): self.fluid = fluid + def density_function(self, lower=None, upper=None): + """Plots the density as a function of temperature in Kelvin + and Pressure in Pascal. + + Parameters + ---------- + lower: tuple + Lower range of the temperature and pressure interval. If None + default values are used. + upper: tuple + Upper range of the temperature and pressure interval. If None + default values are used. + """ + if lower is None: + lower = (100, atm) + if upper is None: + upper = (zero_Celsius + 30, atm * 50) + try: + self.fluid.density_function.plot(lower, upper) + except ValueError: + warnings.warn("Invalid value while attempting density plot.") + def all(self): """Prints out all graphs available about the Fluid. It simply calls all the other plotter methods in this class. @@ -31,3 +58,4 @@ def all(self): ------ None """ + self.density_function() diff --git a/rocketpy/plots/monte_carlo_plots.py b/rocketpy/plots/monte_carlo_plots.py index cfb865c5f..22d320b6b 100644 --- a/rocketpy/plots/monte_carlo_plots.py +++ b/rocketpy/plots/monte_carlo_plots.py @@ -1,5 +1,6 @@ import matplotlib.pyplot as plt import numpy as np +from matplotlib.transforms import offset_copy from ..tools import generate_monte_carlo_ellipses, import_optional_dependency @@ -114,9 +115,20 @@ def ellipses( ) plt.legend() + ax.set_title("1$\\sigma$, 2$\\sigma$ and 3$\\sigma$ Monte Carlo Ellipses") - ax.set_ylabel("North (m)") - ax.set_xlabel("East (m)") + north_south_offset = offset_copy( + ax.transAxes, fig=plt.gcf(), x=-72, y=0, units="points" + ) + east_west_offset = offset_copy( + ax.transAxes, fig=plt.gcf(), x=0, y=-30, units="points" + ) + ax.text(0, 0, "West", va="bottom", ha="center", transform=east_west_offset) + ax.text(1, 0, "East", va="bottom", ha="center", transform=east_west_offset) + ax.text(0, 0, "South", va="bottom", ha="left", transform=north_south_offset) + ax.text(0, 1, "North", va="top", ha="left", transform=north_south_offset) + ax.set_ylabel("Y (m)") + ax.set_xlabel("X (m)") # Add background image to plot # TODO: In the future, integrate with other libraries to plot the map (e.g. cartopy, ee, etc.) @@ -183,6 +195,7 @@ def all(self, keys=None): ax2 = fig.add_subplot(gs[1]) # Plot boxplot + # TODO: changes vert to orientation="horizontal" when support for Py3.9 ends ax1.boxplot(self.monte_carlo.results[key], vert=False) ax1.set_title(f"Box Plot of {key}") ax1.set_yticks([]) diff --git a/rocketpy/prints/fluid_prints.py b/rocketpy/prints/fluid_prints.py index 7b2d7cac8..f1aac5fe2 100644 --- a/rocketpy/prints/fluid_prints.py +++ b/rocketpy/prints/fluid_prints.py @@ -1,3 +1,6 @@ +from ..mathutils.function import NUMERICAL_TYPES + + class _FluidPrints: """Class that holds prints methods for Fluid class. @@ -33,4 +36,7 @@ def all(self): None """ print(f"Name: {self.fluid.name}") - print(f"Density: {self.fluid.density:.4f} kg/m^3") + if isinstance(self.fluid.density, NUMERICAL_TYPES): + print(f"Density: {self.fluid.density:.4f} kg/m^3") + else: + print(f"Density: {self.fluid.density_function}") diff --git a/rocketpy/prints/hybrid_motor_prints.py b/rocketpy/prints/hybrid_motor_prints.py index d527042a3..6f086a162 100644 --- a/rocketpy/prints/hybrid_motor_prints.py +++ b/rocketpy/prints/hybrid_motor_prints.py @@ -1,7 +1,9 @@ import numpy as np +from .motor_prints import _MotorPrints -class _HybridMotorPrints: + +class _HybridMotorPrints(_MotorPrints): """Class that holds prints methods for HybridMotor class. Attributes @@ -26,6 +28,7 @@ def __init__( ------- None """ + super().__init__(hybrid_motor) self.hybrid_motor = hybrid_motor def nozzle_details(self): @@ -63,28 +66,6 @@ def grain_details(self): print(f"Grain Volume: {self.hybrid_motor.solid.grain_initial_volume:.3f} m3") print(f"Grain Mass: {self.hybrid_motor.solid.grain_initial_mass:.3f} kg\n") - def motor_details(self): - """Prints out all data available about the HybridMotor. - - Returns - ------- - None - """ - print("Motor Details") - print(f"Total Burning Time: {self.hybrid_motor.burn_duration} s") - print( - f"Total Propellant Mass: {self.hybrid_motor.propellant_initial_mass:.3f} kg" - ) - print(f"Structural Mass Ratio: {self.hybrid_motor.structural_mass_ratio:.3f}") - avg = self.hybrid_motor.exhaust_velocity.average(*self.hybrid_motor.burn_time) - print(f"Average Propellant Exhaust Velocity: {avg:.3f} m/s") - print(f"Average Thrust: {self.hybrid_motor.average_thrust:.3f} N") - print( - f"Maximum Thrust: {self.hybrid_motor.max_thrust} N at " - f"{self.hybrid_motor.max_thrust_time} s after ignition." - ) - print(f"Total Impulse: {self.hybrid_motor.total_impulse:.3f} Ns\n") - def all(self): """Prints out all data available about the HybridMotor. diff --git a/rocketpy/prints/liquid_motor_prints.py b/rocketpy/prints/liquid_motor_prints.py index 4c80326ab..7411ca59a 100644 --- a/rocketpy/prints/liquid_motor_prints.py +++ b/rocketpy/prints/liquid_motor_prints.py @@ -1,4 +1,7 @@ -class _LiquidMotorPrints: +from .motor_prints import _MotorPrints + + +class _LiquidMotorPrints(_MotorPrints): """Class that holds prints methods for LiquidMotor class. Attributes @@ -23,6 +26,7 @@ def __init__( ------- None """ + super().__init__(liquid_motor) self.liquid_motor = liquid_motor def nozzle_details(self): @@ -35,28 +39,6 @@ def nozzle_details(self): print("Nozzle Details") print("Nozzle Radius: " + str(self.liquid_motor.nozzle_radius) + " m\n") - def motor_details(self): - """Prints out all data available about the motor. - - Returns - ------- - None - """ - print("Motor Details") - print(f"Total Burning Time: {self.liquid_motor.burn_duration} s") - print( - f"Total Propellant Mass: {self.liquid_motor.propellant_initial_mass:.3f} kg" - ) - print(f"Structural Mass Ratio: {self.liquid_motor.structural_mass_ratio:.3f}") - avg = self.liquid_motor.exhaust_velocity.average(*self.liquid_motor.burn_time) - print(f"Average Propellant Exhaust Velocity: {avg:.3f} m/s") - print(f"Average Thrust: {self.liquid_motor.average_thrust:.3f} N") - print( - f"Maximum Thrust: {self.liquid_motor.max_thrust} N at " - f"{self.liquid_motor.max_thrust_time} s after ignition." - ) - print(f"Total Impulse: {self.liquid_motor.total_impulse:.3f} Ns\n") - def all(self): """Prints out all data available about the LiquidMotor. diff --git a/rocketpy/prints/solid_motor_prints.py b/rocketpy/prints/solid_motor_prints.py index 6f4c28d5b..16ef5dc38 100644 --- a/rocketpy/prints/solid_motor_prints.py +++ b/rocketpy/prints/solid_motor_prints.py @@ -1,4 +1,7 @@ -class _SolidMotorPrints: +from .motor_prints import _MotorPrints + + +class _SolidMotorPrints(_MotorPrints): """Class that holds prints methods for SolidMotor class. Attributes @@ -23,6 +26,7 @@ def __init__( ------- None """ + super().__init__(solid_motor) self.solid_motor = solid_motor def nozzle_details(self): @@ -53,28 +57,6 @@ def grain_details(self): print(f"Grain Volume: {self.solid_motor.grain_initial_volume:.3f} m3") print(f"Grain Mass: {self.solid_motor.grain_initial_mass:.3f} kg\n") - def motor_details(self): - """Prints out all data available about the SolidMotor. - - Returns - ------- - None - """ - print("Motor Details") - print("Total Burning Time: " + str(self.solid_motor.burn_duration) + " s") - print( - f"Total Propellant Mass: {self.solid_motor.propellant_initial_mass:.3f} kg" - ) - print(f"Structural Mass Ratio: {self.solid_motor.structural_mass_ratio:.3f}") - average = self.solid_motor.exhaust_velocity.average(*self.solid_motor.burn_time) - print(f"Average Propellant Exhaust Velocity: {average:.3f} m/s") - print(f"Average Thrust: {self.solid_motor.average_thrust:.3f} N") - print( - f"Maximum Thrust: {self.solid_motor.max_thrust} N " - f"at {self.solid_motor.max_thrust_time} s after ignition." - ) - print(f"Total Impulse: {self.solid_motor.total_impulse:.3f} Ns\n") - def all(self): """Prints out all data available about the SolidMotor. diff --git a/rocketpy/rocket/aero_surface/air_brakes.py b/rocketpy/rocket/aero_surface/air_brakes.py index ee4830808..2fa5c782c 100644 --- a/rocketpy/rocket/aero_surface/air_brakes.py +++ b/rocketpy/rocket/aero_surface/air_brakes.py @@ -131,7 +131,8 @@ def deployment_level(self, value): warnings.warn( f"Deployment level of {self.name} is smaller than 0 or " + "larger than 1. Extrapolation for the drag coefficient " - + "curve will be used." + + "curve will be used.", + UserWarning, ) self._deployment_level = value @@ -205,3 +206,24 @@ def all_info(self): """ self.info() self.plots.drag_coefficient_curve() + + def to_dict(self, **kwargs): # pylint: disable=unused-argument + return { + "drag_coefficient_curve": self.drag_coefficient, + "reference_area": self.reference_area, + "clamp": self.clamp, + "override_rocket_drag": self.override_rocket_drag, + "deployment_level": self.initial_deployment_level, + "name": self.name, + } + + @classmethod + def from_dict(cls, data): + return cls( + drag_coefficient_curve=data.get("drag_coefficient_curve"), + reference_area=data.get("reference_area"), + clamp=data.get("clamp"), + override_rocket_drag=data.get("override_rocket_drag"), + deployment_level=data.get("deployment_level"), + name=data.get("name"), + ) diff --git a/rocketpy/rocket/aero_surface/fins/elliptical_fins.py b/rocketpy/rocket/aero_surface/fins/elliptical_fins.py index 08c51ab64..61d98bc0d 100644 --- a/rocketpy/rocket/aero_surface/fins/elliptical_fins.py +++ b/rocketpy/rocket/aero_surface/fins/elliptical_fins.py @@ -317,9 +317,9 @@ def all_info(self): self.prints.all() self.plots.all() - def to_dict(self, include_outputs=False): - data = super().to_dict(include_outputs) - if include_outputs: + def to_dict(self, **kwargs): + data = super().to_dict(**kwargs) + if kwargs.get("include_outputs", False): data.update( { "Af": self.Af, diff --git a/rocketpy/rocket/aero_surface/fins/fins.py b/rocketpy/rocket/aero_surface/fins/fins.py index 26a4a7111..b401164a0 100644 --- a/rocketpy/rocket/aero_surface/fins/fins.py +++ b/rocketpy/rocket/aero_surface/fins/fins.py @@ -426,22 +426,40 @@ def compute_forces_and_moments( M3 = M3_forcing - M3_damping return R1, R2, R3, M1, M2, M3 - def to_dict(self, include_outputs=False): + def to_dict(self, **kwargs): + if self.airfoil: + if kwargs.get("discretize", False): + lower = -np.pi / 6 if self.airfoil[1] == "radians" else -30 + upper = np.pi / 6 if self.airfoil[1] == "radians" else 30 + airfoil = ( + self.airfoil_cl.set_discrete(lower, upper, 50, mutate_self=False), + self.airfoil[1], + ) + else: + airfoil = (self.airfoil_cl, self.airfoil[1]) if self.airfoil else None + else: + airfoil = None data = { "n": self.n, "root_chord": self.root_chord, "span": self.span, "rocket_radius": self.rocket_radius, "cant_angle": self.cant_angle, - "airfoil": self.airfoil, + "airfoil": airfoil, "name": self.name, } - if include_outputs: + if kwargs.get("include_outputs", False): + cl = self.cl + if kwargs.get("discretize", False): + cl = cl.set_discrete( + (-np.pi / 6, 0), (np.pi / 6, 2), (10, 10), mutate_self=False + ) + data.update( { "cp": self.cp, - "cl": self.cl, + "cl": cl, "roll_parameters": self.roll_parameters, "d": self.d, "ref_area": self.ref_area, diff --git a/rocketpy/rocket/aero_surface/fins/free_form_fins.py b/rocketpy/rocket/aero_surface/fins/free_form_fins.py index 7cae4e556..72758171e 100644 --- a/rocketpy/rocket/aero_surface/fins/free_form_fins.py +++ b/rocketpy/rocket/aero_surface/fins/free_form_fins.py @@ -359,11 +359,11 @@ def evaluate_shape(self): x_array, y_array = zip(*self.shape_points) self.shape_vec = [np.array(x_array), np.array(y_array)] - def to_dict(self, include_outputs=False): - data = super().to_dict(include_outputs) + def to_dict(self, **kwargs): + data = super().to_dict(**kwargs) data["shape_points"] = self.shape_points - if include_outputs: + if kwargs.get("include_outputs", False): data.update( { "Af": self.Af, diff --git a/rocketpy/rocket/aero_surface/fins/trapezoidal_fins.py b/rocketpy/rocket/aero_surface/fins/trapezoidal_fins.py index 61e2b78fc..c6b4ea633 100644 --- a/rocketpy/rocket/aero_surface/fins/trapezoidal_fins.py +++ b/rocketpy/rocket/aero_surface/fins/trapezoidal_fins.py @@ -176,9 +176,6 @@ def __init__( sweep_length = np.tan(sweep_angle * np.pi / 180) * span elif sweep_length is None: sweep_length = root_chord - tip_chord - else: - # Sweep length is given - pass self._tip_chord = tip_chord self._sweep_length = sweep_length @@ -348,15 +345,15 @@ def all_info(self): self.prints.all() self.plots.all() - def to_dict(self, include_outputs=False): - data = super().to_dict(include_outputs) + def to_dict(self, **kwargs): + data = super().to_dict(**kwargs) data["tip_chord"] = self.tip_chord + data["sweep_length"] = self.sweep_length + data["sweep_angle"] = self.sweep_angle - if include_outputs: + if kwargs.get("include_outputs", False): data.update( { - "sweep_length": self.sweep_length, - "sweep_angle": self.sweep_angle, "shape_vec": self.shape_vec, "Af": self.Af, "AR": self.AR, @@ -382,4 +379,5 @@ def from_dict(cls, data): cant_angle=data["cant_angle"], airfoil=data["airfoil"], name=data["name"], + sweep_length=data.get("sweep_length"), ) diff --git a/rocketpy/rocket/aero_surface/nose_cone.py b/rocketpy/rocket/aero_surface/nose_cone.py index 954ce1ef8..4e6ea8e95 100644 --- a/rocketpy/rocket/aero_surface/nose_cone.py +++ b/rocketpy/rocket/aero_surface/nose_cone.py @@ -523,7 +523,7 @@ def all_info(self): self.prints.all() self.plots.all() - def to_dict(self, include_outputs=False): + def to_dict(self, **kwargs): data = { "_length": self._length, "_kind": self._kind, @@ -533,10 +533,17 @@ def to_dict(self, include_outputs=False): "_power": self._power, "name": self.name, } - if include_outputs: + if kwargs.get("include_outputs", False): + clalpha = self.clalpha + cl = self.cl + if kwargs.get("discretize", False): + clalpha = clalpha.set_discrete(0, 4, 50) + cl = cl.set_discrete( + (-np.pi / 6, 0), (np.pi / 6, 2), (10, 10), mutate_self=False + ) data["cp"] = self.cp - data["clalpha"] = self.clalpha - data["cl"] = self.cl + data["clalpha"] = clalpha + data["cl"] = cl return data diff --git a/rocketpy/rocket/aero_surface/rail_buttons.py b/rocketpy/rocket/aero_surface/rail_buttons.py index 879f013d7..c27e2949a 100644 --- a/rocketpy/rocket/aero_surface/rail_buttons.py +++ b/rocketpy/rocket/aero_surface/rail_buttons.py @@ -100,7 +100,7 @@ def evaluate_geometrical_parameters(self): None """ - def to_dict(self, include_outputs=False): # pylint: disable=unused-argument + def to_dict(self, **kwargs): # pylint: disable=unused-argument return { "buttons_distance": self.buttons_distance, "angular_position": self.angular_position, diff --git a/rocketpy/rocket/aero_surface/tail.py b/rocketpy/rocket/aero_surface/tail.py index d32836e18..b5e4b878a 100644 --- a/rocketpy/rocket/aero_surface/tail.py +++ b/rocketpy/rocket/aero_surface/tail.py @@ -205,7 +205,7 @@ def all_info(self): self.prints.all() self.plots.all() - def to_dict(self, include_outputs=False): + def to_dict(self, **kwargs): data = { "top_radius": self._top_radius, "bottom_radius": self._bottom_radius, @@ -214,11 +214,20 @@ def to_dict(self, include_outputs=False): "name": self.name, } - if include_outputs: + if kwargs.get("include_outputs", False): + clalpha = self.clalpha + cl = self.cl + if kwargs.get("discretize", False): + clalpha = clalpha.set_discrete(0, 4, 50) + cl = cl.set_discrete( + (-np.pi / 6, 0), (np.pi / 6, 2), (10, 10), mutate_self=False + ) + data.update( { "cp": self.cp, - "cl": self.clalpha, + "clalpha": clalpha, + "cl": cl, "slant_length": self.slant_length, "surface_area": self.surface_area, } diff --git a/rocketpy/rocket/components.py b/rocketpy/rocket/components.py index 66448eb69..9ce8595b3 100644 --- a/rocketpy/rocket/components.py +++ b/rocketpy/rocket/components.py @@ -200,7 +200,7 @@ def sort_by_position(self, reverse=False): """ self._components.sort(key=lambda x: x.position.z, reverse=reverse) - def to_dict(self, include_outputs=False): # pylint: disable=unused-argument + def to_dict(self, **kwargs): # pylint: disable=unused-argument return { "components": [ {"component": c.component, "position": c.position} diff --git a/rocketpy/rocket/parachute.py b/rocketpy/rocket/parachute.py index a18cb4793..e27216e26 100644 --- a/rocketpy/rocket/parachute.py +++ b/rocketpy/rocket/parachute.py @@ -9,7 +9,7 @@ class Parachute: - """Keeps parachute information. + """Keeps information of the parachute, which is modeled as a hemispheroid. Attributes ---------- @@ -92,6 +92,20 @@ class Parachute: Function of noisy_pressure_signal. Parachute.clean_pressure_signal_function : Function Function of clean_pressure_signal. + Parachute.radius : float + Length of the non-unique semi-axis (radius) of the inflated hemispheroid + parachute in meters. + Parachute.height : float, None + Length of the unique semi-axis (height) of the inflated hemispheroid + parachute in meters. + Parachute.porosity : float + Geometric porosity of the canopy (ratio of open area to total canopy area), + in [0, 1]. Affects only the added-mass scaling during descent; it does + not change ``cd_s`` (drag). The default, 0.0432, yields an added-mass + of 1.0 (“neutral” behavior). + Parachute.added_mass_coefficient : float + Coefficient used to calculate the added-mass due to dragged air. It is + calculated from the porosity of the parachute. """ def __init__( @@ -102,6 +116,9 @@ def __init__( sampling_rate, lag=0, noise=(0, 0, 0), + radius=1.5, + height=None, + porosity=0.0432, ): """Initializes Parachute class. @@ -154,6 +171,19 @@ def __init__( The values are used to add noise to the pressure signal which is passed to the trigger function. Default value is ``(0, 0, 0)``. Units are in Pa. + radius : float, optional + Length of the non-unique semi-axis (radius) of the inflated hemispheroid + parachute. Default value is 1.5. + Units are in meters. + height : float, optional + Length of the unique semi-axis (height) of the inflated hemispheroid + parachute. Default value is the radius of the parachute. + Units are in meters. + porosity : float, optional + Geometric porosity of the canopy (ratio of open area to total canopy area), + in [0, 1]. Affects only the added-mass scaling during descent; it does + not change ``cd_s`` (drag). The default, 0.0432, yields an added-mass + of 1.0 (“neutral” behavior). """ self.name = name self.cd_s = cd_s @@ -170,6 +200,15 @@ def __init__( self.clean_pressure_signal_function = Function(0) self.noisy_pressure_signal_function = Function(0) self.noise_signal_function = Function(0) + self.radius = radius + self.height = height or radius + self.porosity = porosity + self.added_mass_coefficient = 1.068 * ( + 1 + - 1.465 * self.porosity + - 0.25975 * self.porosity**2 + + 1.2626 * self.porosity**3 + ) alpha, beta = self.noise_corr self.noise_function = lambda: alpha * self.noise_signal[-1][ @@ -251,11 +290,15 @@ def all_info(self): self.info() # self.plots.all() # Parachutes still doesn't have plots - def to_dict(self, include_outputs=False): + def to_dict(self, **kwargs): + allow_pickle = kwargs.get("allow_pickle", True) trigger = self.trigger if callable(self.trigger) and not isinstance(self.trigger, Function): - trigger = to_hex_encode(trigger) + if allow_pickle: + trigger = to_hex_encode(trigger) + else: + trigger = trigger.__name__ data = { "name": self.name, @@ -264,11 +307,18 @@ def to_dict(self, include_outputs=False): "sampling_rate": self.sampling_rate, "lag": self.lag, "noise": self.noise, + "radius": self.radius, + "height": self.height, + "porosity": self.porosity, } - if include_outputs: + if kwargs.get("include_outputs", False): data["noise_signal"] = self.noise_signal - data["noise_function"] = to_hex_encode(self.noise_function) + data["noise_function"] = ( + to_hex_encode(self.noise_function) + if allow_pickle + else self.noise_function.__name__ + ) data["noisy_pressure_signal"] = self.noisy_pressure_signal data["clean_pressure_signal"] = self.clean_pressure_signal @@ -290,6 +340,9 @@ def from_dict(cls, data): sampling_rate=data["sampling_rate"], lag=data["lag"], noise=data["noise"], + radius=data.get("radius", 1.5), + height=data.get("height", None), + porosity=data.get("porosity", 0.0432), ) return parachute diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index bf938d4be..1112a98f3 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -1,5 +1,6 @@ import math import warnings +from typing import Iterable import numpy as np @@ -22,7 +23,11 @@ from rocketpy.rocket.aero_surface.generic_surface import GenericSurface from rocketpy.rocket.components import Components from rocketpy.rocket.parachute import Parachute -from rocketpy.tools import parallel_axis_theorem_from_com +from rocketpy.tools import ( + deprecated, + find_obj_from_hash, + parallel_axis_theorem_from_com, +) # pylint: disable=too-many-instance-attributes, too-many-public-methods, too-many-instance-attributes @@ -1173,16 +1178,16 @@ def add_nose( self.add_surfaces(nose, position) return nose + @deprecated( + reason="This method is set to be deprecated in version 1.0.0 and fully " + "removed by version 2.0.0", + alternative="Rocket.add_trapezoidal_fins", + ) def add_fins(self, *args, **kwargs): # pragma: no cover """See Rocket.add_trapezoidal_fins for documentation. This method is set to be deprecated in version 1.0.0 and fully removed by version 2.0.0. Use Rocket.add_trapezoidal_fins instead. It keeps the same arguments and signature.""" - warnings.warn( - "This method is set to be deprecated in version 1.0.0 and fully " - "removed by version 2.0.0. Use Rocket.add_trapezoidal_fins instead", - DeprecationWarning, - ) return self.add_trapezoidal_fins(*args, **kwargs) def add_trapezoidal_fins( @@ -1434,7 +1439,16 @@ def add_free_form_fins( return fin_set def add_parachute( - self, name, cd_s, trigger, sampling_rate=100, lag=0, noise=(0, 0, 0) + self, + name, + cd_s, + trigger, + sampling_rate=100, + lag=0, + noise=(0, 0, 0), + radius=1.5, + height=None, + porosity=0.0432, ): """Creates a new parachute, storing its parameters such as opening delay, drag coefficients and trigger function. @@ -1493,16 +1507,39 @@ def add_parachute( The values are used to add noise to the pressure signal which is passed to the trigger function. Default value is (0, 0, 0). Units are in pascal. + radius : float, optional + Length of the non-unique semi-axis (radius) of the inflated hemispheroid + parachute. Default value is 1.5. + Units are in meters. + height : float, optional + Length of the unique semi-axis (height) of the inflated hemispheroid + parachute. Default value is the radius of the parachute. + Units are in meters. + porosity : float, optional + Geometric porosity of the canopy (ratio of open area to total canopy area), + in [0, 1]. Affects only the added-mass scaling during descent; it does + not change ``cd_s`` (drag). The default, 0.0432, yields an added-mass + of 1.0 (“neutral” behavior). Returns ------- parachute : Parachute - Parachute containing trigger, sampling_rate, lag, cd_s, noise - and name. Furthermore, it stores clean_pressure_signal, + Parachute containing trigger, sampling_rate, lag, cd_s, noise, radius, + height, porosity and name. Furthermore, it stores clean_pressure_signal, noise_signal and noisyPressureSignal which are filled in during Flight simulation. """ - parachute = Parachute(name, cd_s, trigger, sampling_rate, lag, noise) + parachute = Parachute( + name, + cd_s, + trigger, + sampling_rate, + lag, + noise, + radius, + height, + porosity, + ) self.parachutes.append(parachute) return self.parachutes[-1] @@ -1892,7 +1929,16 @@ def all_info(self): self.info() self.plots.all() - def to_dict(self, include_outputs=False): + # pylint: disable=too-many-statements + def to_dict(self, **kwargs): + discretize = kwargs.get("discretize", False) + + power_off_drag = self.power_off_drag + power_on_drag = self.power_on_drag + if discretize: + power_off_drag = power_off_drag.set_discrete(0, 4, 50, mutate_self=False) + power_on_drag = power_on_drag.set_discrete(0, 4, 50, mutate_self=False) + rocket_dict = { "radius": self.radius, "mass": self.mass, @@ -1902,8 +1948,8 @@ def to_dict(self, include_outputs=False): "I_12_without_motor": self.I_12_without_motor, "I_13_without_motor": self.I_13_without_motor, "I_23_without_motor": self.I_23_without_motor, - "power_off_drag": self.power_off_drag, - "power_on_drag": self.power_on_drag, + "power_off_drag": power_off_drag, + "power_on_drag": power_on_drag, "center_of_mass_without_motor": self.center_of_mass_without_motor, "coordinate_system_orientation": self.coordinate_system_orientation, "motor": self.motor, @@ -1916,7 +1962,51 @@ def to_dict(self, include_outputs=False): "sensors": self.sensors, } - if include_outputs: + if kwargs.get("include_outputs", False): + thrust_to_weight = self.thrust_to_weight + cp_position = self.cp_position + stability_margin = self.stability_margin + center_of_mass = self.center_of_mass + motor_center_of_mass_position = self.motor_center_of_mass_position + reduced_mass = self.reduced_mass + total_mass = self.total_mass + total_mass_flow_rate = self.total_mass_flow_rate + center_of_propellant_position = self.center_of_propellant_position + + if discretize: + thrust_to_weight = thrust_to_weight.set_discrete_based_on_model( + self.motor.thrust, mutate_self=False + ) + cp_position = cp_position.set_discrete(0, 4, 25, mutate_self=False) + stability_margin = stability_margin.set_discrete( + (0, self.motor.burn_time[0]), + (2, self.motor.burn_time[1]), + (10, 10), + mutate_self=False, + ) + center_of_mass = center_of_mass.set_discrete_based_on_model( + self.motor.thrust, mutate_self=False + ) + motor_center_of_mass_position = ( + motor_center_of_mass_position.set_discrete_based_on_model( + self.motor.thrust, mutate_self=False + ) + ) + reduced_mass = reduced_mass.set_discrete_based_on_model( + self.motor.thrust, mutate_self=False + ) + total_mass = total_mass.set_discrete_based_on_model( + self.motor.thrust, mutate_self=False + ) + total_mass_flow_rate = total_mass_flow_rate.set_discrete_based_on_model( + self.motor.thrust, mutate_self=False + ) + center_of_propellant_position = ( + center_of_propellant_position.set_discrete_based_on_model( + self.motor.thrust, mutate_self=False + ) + ) + rocket_dict["area"] = self.area rocket_dict["center_of_dry_mass_position"] = ( self.center_of_dry_mass_position @@ -1924,30 +2014,26 @@ def to_dict(self, include_outputs=False): rocket_dict["center_of_mass_without_motor"] = ( self.center_of_mass_without_motor ) - rocket_dict["motor_center_of_mass_position"] = ( - self.motor_center_of_mass_position - ) + rocket_dict["motor_center_of_mass_position"] = motor_center_of_mass_position rocket_dict["motor_center_of_dry_mass_position"] = ( self.motor_center_of_dry_mass_position ) - rocket_dict["center_of_mass"] = self.center_of_mass - rocket_dict["reduced_mass"] = self.reduced_mass - rocket_dict["total_mass"] = self.total_mass - rocket_dict["total_mass_flow_rate"] = self.total_mass_flow_rate - rocket_dict["thrust_to_weight"] = self.thrust_to_weight + rocket_dict["center_of_mass"] = center_of_mass + rocket_dict["reduced_mass"] = reduced_mass + rocket_dict["total_mass"] = total_mass + rocket_dict["total_mass_flow_rate"] = total_mass_flow_rate + rocket_dict["thrust_to_weight"] = thrust_to_weight rocket_dict["cp_eccentricity_x"] = self.cp_eccentricity_x rocket_dict["cp_eccentricity_y"] = self.cp_eccentricity_y rocket_dict["thrust_eccentricity_x"] = self.thrust_eccentricity_x rocket_dict["thrust_eccentricity_y"] = self.thrust_eccentricity_y - rocket_dict["cp_position"] = self.cp_position - rocket_dict["stability_margin"] = self.stability_margin + rocket_dict["cp_position"] = cp_position + rocket_dict["stability_margin"] = stability_margin rocket_dict["static_margin"] = self.static_margin rocket_dict["nozzle_position"] = self.nozzle_position rocket_dict["nozzle_to_cdm"] = self.nozzle_to_cdm rocket_dict["nozzle_gyration_tensor"] = self.nozzle_gyration_tensor - rocket_dict["center_of_propellant_position"] = ( - self.center_of_propellant_position - ) + rocket_dict["center_of_propellant_position"] = center_of_propellant_position return rocket_dict @@ -1990,17 +2076,29 @@ def from_dict(cls, data): for parachute in data["parachutes"]: rocket.parachutes.append(parachute) - for air_brakes in data["air_brakes"]: - rocket.add_air_brakes( - drag_coefficient_curve=air_brakes["drag_coefficient_curve"], - controller_function=air_brakes["controller_function"], - sampling_rate=air_brakes["sampling_rate"], - clamp=air_brakes["clamp"], - reference_area=air_brakes["reference_area"], - initial_observed_variables=air_brakes["initial_observed_variables"], - override_rocket_drag=air_brakes["override_rocket_drag"], - name=air_brakes["name"], - controller_name=air_brakes["controller_name"], - ) + for sensor, position in data["sensors"]: + rocket.add_sensor(sensor, position) + + for air_brake in data["air_brakes"]: + rocket.air_brakes.append(air_brake) + + for controller in data["_controllers"]: + interactive_objects_hash = getattr(controller, "_interactive_objects_hash") + if interactive_objects_hash is not None: + is_iterable = isinstance(interactive_objects_hash, Iterable) + if not is_iterable: + interactive_objects_hash = [interactive_objects_hash] + for hash_ in interactive_objects_hash: + if (hashed_obj := find_obj_from_hash(data, hash_)) is not None: + if not is_iterable: + controller.interactive_objects = hashed_obj + else: + controller.interactive_objects.append(hashed_obj) + else: + warnings.warn( + "Could not find controller interactive objects." + "Deserialization will proceed, results may not be accurate." + ) + rocket._add_controllers(controller) return rocket diff --git a/rocketpy/sensors/accelerometer.py b/rocketpy/sensors/accelerometer.py index de249c173..1037c0e9a 100644 --- a/rocketpy/sensors/accelerometer.py +++ b/rocketpy/sensors/accelerometer.py @@ -275,3 +275,28 @@ def export_measured_data(self, filename, file_format="csv"): file_format=file_format, data_labels=("t", "ax", "ay", "az"), ) + + def to_dict(self, **kwargs): + data = super().to_dict(**kwargs) + data.update({"consider_gravity": self.consider_gravity}) + return data + + @classmethod + def from_dict(cls, data): + return cls( + sampling_rate=data["sampling_rate"], + orientation=data["orientation"], + measurement_range=data["measurement_range"], + resolution=data["resolution"], + noise_density=data["noise_density"], + noise_variance=data["noise_variance"], + random_walk_density=data["random_walk_density"], + random_walk_variance=data["random_walk_variance"], + constant_bias=data["constant_bias"], + operating_temperature=data["operating_temperature"], + temperature_bias=data["temperature_bias"], + temperature_scale_factor=data["temperature_scale_factor"], + cross_axis_sensitivity=data["cross_axis_sensitivity"], + consider_gravity=data["consider_gravity"], + name=data["name"], + ) diff --git a/rocketpy/sensors/barometer.py b/rocketpy/sensors/barometer.py index bf5b538e2..cc26bb056 100644 --- a/rocketpy/sensors/barometer.py +++ b/rocketpy/sensors/barometer.py @@ -190,3 +190,20 @@ def export_measured_data(self, filename, file_format="csv"): file_format=file_format, data_labels=("t", "pressure"), ) + + @classmethod + def from_dict(cls, data): + return cls( + sampling_rate=data["sampling_rate"], + measurement_range=data["measurement_range"], + resolution=data["resolution"], + noise_density=data["noise_density"], + noise_variance=data["noise_variance"], + random_walk_density=data["random_walk_density"], + random_walk_variance=data["random_walk_variance"], + constant_bias=data["constant_bias"], + operating_temperature=data["operating_temperature"], + temperature_bias=data["temperature_bias"], + temperature_scale_factor=data["temperature_scale_factor"], + name=data["name"], + ) diff --git a/rocketpy/sensors/gnss_receiver.py b/rocketpy/sensors/gnss_receiver.py index 686cd29e5..548f9c879 100644 --- a/rocketpy/sensors/gnss_receiver.py +++ b/rocketpy/sensors/gnss_receiver.py @@ -124,3 +124,20 @@ def export_measured_data(self, filename, file_format="csv"): file_format=file_format, data_labels=("t", "latitude", "longitude", "altitude"), ) + + def to_dict(self, **kwargs): + return { + "sampling_rate": self.sampling_rate, + "position_accuracy": self.position_accuracy, + "altitude_accuracy": self.altitude_accuracy, + "name": self.name, + } + + @classmethod + def from_dict(cls, data): + return cls( + sampling_rate=data["sampling_rate"], + position_accuracy=data["position_accuracy"], + altitude_accuracy=data["altitude_accuracy"], + name=data["name"], + ) diff --git a/rocketpy/sensors/gyroscope.py b/rocketpy/sensors/gyroscope.py index 5acbb6f50..6a18af601 100644 --- a/rocketpy/sensors/gyroscope.py +++ b/rocketpy/sensors/gyroscope.py @@ -296,3 +296,28 @@ def export_measured_data(self, filename, file_format="csv"): file_format=file_format, data_labels=("t", "wx", "wy", "wz"), ) + + def to_dict(self, **kwargs): + data = super().to_dict(**kwargs) + data.update({"acceleration_sensitivity": self.acceleration_sensitivity}) + return data + + @classmethod + def from_dict(cls, data): + return cls( + sampling_rate=data["sampling_rate"], + orientation=data["orientation"], + measurement_range=data["measurement_range"], + resolution=data["resolution"], + noise_density=data["noise_density"], + noise_variance=data["noise_variance"], + random_walk_density=data["random_walk_density"], + random_walk_variance=data["random_walk_variance"], + constant_bias=data["constant_bias"], + operating_temperature=data["operating_temperature"], + temperature_bias=data["temperature_bias"], + temperature_scale_factor=data["temperature_scale_factor"], + cross_axis_sensitivity=data["cross_axis_sensitivity"], + acceleration_sensitivity=data["acceleration_sensitivity"], + name=data["name"], + ) diff --git a/rocketpy/sensors/sensor.py b/rocketpy/sensors/sensor.py index 416201019..0b44aeb18 100644 --- a/rocketpy/sensors/sensor.py +++ b/rocketpy/sensors/sensor.py @@ -271,6 +271,23 @@ def _generic_export_measured_data(self, filename, file_format, data_labels): print(f"Data saved to {filename}") return + # pylint: disable=unused-argument + def to_dict(self, **kwargs): + return { + "sampling_rate": self.sampling_rate, + "measurement_range": self.measurement_range, + "resolution": self.resolution, + "operating_temperature": self.operating_temperature, + "noise_density": self.noise_density, + "noise_variance": self.noise_variance, + "random_walk_density": self.random_walk_density, + "random_walk_variance": self.random_walk_variance, + "constant_bias": self.constant_bias, + "temperature_bias": self.temperature_bias, + "temperature_scale_factor": self.temperature_scale_factor, + "name": self.name, + } + class InertialSensor(Sensor): """Model of an inertial sensor (accelerometer, gyroscope, magnetometer). @@ -574,6 +591,16 @@ def apply_temperature_drift(self, value): ) return value & scale_factor + def to_dict(self, **kwargs): + data = super().to_dict(**kwargs) + data.update( + { + "orientation": self.orientation, + "cross_axis_sensitivity": self.cross_axis_sensitivity, + } + ) + return data + class ScalarSensor(Sensor): """Model of a scalar sensor (e.g. Barometer). Scalar sensors are used diff --git a/rocketpy/simulation/__init__.py b/rocketpy/simulation/__init__.py index 6b98fdcf4..1ade0f16f 100644 --- a/rocketpy/simulation/__init__.py +++ b/rocketpy/simulation/__init__.py @@ -1,4 +1,5 @@ from .flight import Flight +from .flight_data_exporter import FlightDataExporter from .flight_data_importer import FlightDataImporter from .monte_carlo import MonteCarlo from .multivariate_rejection_sampler import MultivariateRejectionSampler diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index ce728dafe..a38be7d93 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1,20 +1,21 @@ # pylint: disable=too-many-lines -import json import math import warnings from copy import deepcopy from functools import cached_property import numpy as np -import simplekml from scipy.integrate import BDF, DOP853, LSODA, RK23, RK45, OdeSolver, Radau +from rocketpy.simulation.flight_data_exporter import FlightDataExporter + from ..mathutils.function import Function, funcify_method from ..mathutils.vector_matrix import Matrix, Vector from ..plots.flight_plots import _FlightPlots from ..prints.flight_prints import _FlightPrints from ..tools import ( calculate_cubic_hermite_coefficients, + deprecated, euler313_to_quaternions, find_closest, find_root_linear_interpolation, @@ -143,7 +144,7 @@ class Flight: Flight.heading : int, float Launch heading angle relative to north given in degrees. Flight.initial_solution : list - List defines initial condition - [tInit, x_init, + List defines initial condition - [t_initial, x_init, y_init, z_init, vx_init, vy_init, vz_init, e0_init, e1_init, e2_init, e3_init, w1_init, w2_init, w3_init] Flight.t_initial : int, float @@ -155,12 +156,6 @@ class Flight: Current integration time. Flight.y : list Current integration state vector u. - Flight.post_processed : bool - Defines if solution data has been post processed. - Flight.initial_solution : list - List defines initial condition - [tInit, x_init, - y_init, z_init, vx_init, vy_init, vz_init, e0_init, e1_init, - e2_init, e3_init, w1_init, w2_init, w3_init] Flight.out_of_rail_time : int, float Time, in seconds, in which the rocket completely leaves the rail. @@ -546,7 +541,7 @@ def __init__( # pylint: disable=too-many-arguments,too-many-statements rtol : float, array, optional Maximum relative error tolerance to be tolerated in the integration scheme. Can be given as array for each - state space variable. Default is 1e-3. + state space variable. Default is 1e-6. atol : float, optional Maximum absolute error tolerance to be tolerated in the integration scheme. Can be given as array for each @@ -767,7 +762,22 @@ def __simulate(self, verbose): callbacks = [ lambda self, parachute_cd_s=parachute.cd_s: setattr( self, "parachute_cd_s", parachute_cd_s - ) + ), + lambda self, parachute_radius=parachute.radius: setattr( + self, "parachute_radius", parachute_radius + ), + lambda self, parachute_height=parachute.height: setattr( + self, "parachute_height", parachute_height + ), + lambda self, parachute_porosity=parachute.porosity: setattr( + self, "parachute_porosity", parachute_porosity + ), + lambda self, + added_mass_coefficient=parachute.added_mass_coefficient: setattr( + self, + "parachute_added_mass_coefficient", + added_mass_coefficient, + ), ] self.flight_phases.add_phase( node.t + parachute.lag, @@ -1013,7 +1023,31 @@ def __simulate(self, verbose): lambda self, parachute_cd_s=parachute.cd_s: setattr( self, "parachute_cd_s", parachute_cd_s - ) + ), + lambda self, + parachute_radius=parachute.radius: setattr( + self, + "parachute_radius", + parachute_radius, + ), + lambda self, + parachute_height=parachute.height: setattr( + self, + "parachute_height", + parachute_height, + ), + lambda self, + parachute_porosity=parachute.porosity: setattr( + self, + "parachute_porosity", + parachute_porosity, + ), + lambda self, + added_mass_coefficient=parachute.added_mass_coefficient: setattr( + self, + "parachute_added_mass_coefficient", + added_mass_coefficient, + ), ] self.flight_phases.add_phase( overshootable_node.t + parachute.lag, @@ -1096,13 +1130,13 @@ def __init_solution_monitors(self): self.out_of_rail_time_index = 0 self.out_of_rail_state = np.array([0]) self.apogee_state = np.array([0]) + self.apogee = 0 self.apogee_time = 0 self.x_impact = 0 self.y_impact = 0 self.impact_velocity = 0 self.impact_state = np.array([0]) self.parachute_events = [] - self.post_processed = False self.__post_processed_variables = [] def __init_flight_state(self): @@ -1171,6 +1205,7 @@ def __init_flight_state(self): self.out_of_rail_state = self.initial_solution[1:] self.out_of_rail_time = self.initial_solution[0] self.out_of_rail_time_index = 0 + self.t_initial = self.initial_solution[0] self.initial_derivative = self.u_dot_generalized if self._controllers or self.sensors: # Handle post process during simulation, get initial accel/forces @@ -1679,6 +1714,12 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals ax, ay, az = K @ Vector(L) az -= self.env.gravity.get_value_opt(z) # Include gravity + # Coriolis acceleration + _, w_earth_y, w_earth_z = self.env.earth_rotation_vector + ax -= 2 * (vz * w_earth_y - vy * w_earth_z) + ay -= 2 * (vx * w_earth_z) + az -= 2 * (-vx * w_earth_y) + # Create u_dot u_dot = [ vx, @@ -1745,7 +1786,7 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too _, _, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u # Create necessary vectors - # r = Vector([x, y, z]) # CDM position vector + # r = Vector([x, y, z]) # CDM position vector v = Vector([vx, vy, vz]) # CDM velocity vector e = [e0, e1, e2, e3] # Euler parameters/quaternions w = Vector([omega1, omega2, omega3]) # Angular velocity vector @@ -1904,9 +1945,6 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too # Angular velocity derivative w_dot = I_CM.inverse @ (T21 + (T20 ^ r_CM)) - # Velocity vector derivative - v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - # Euler parameters derivative e_dot = [ 0.5 * (-omega1 * e1 - omega2 * e2 - omega3 * e3), @@ -1915,6 +1953,10 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too 0.5 * (omega3 * e0 + omega2 * e1 - omega1 * e2), ] + # Velocity vector derivative + Coriolis acceleration + w_earth = Vector(self.env.earth_rotation_vector) + v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - 2 * (w_earth ^ v) + # Position vector derivative r_dot = [vx, vy, vz] @@ -1959,15 +2001,9 @@ def u_dot_parachute(self, t, u, post_processing=False): wind_velocity_x = self.env.wind_velocity_x.get_value_opt(z) wind_velocity_y = self.env.wind_velocity_y.get_value_opt(z) - # Get Parachute data - cd_s = self.parachute_cd_s - # Get the mass of the rocket mp = self.rocket.dry_mass - # Define constants - ka = 1 # Added mass coefficient (depends on parachute's porosity) - R = 1.5 # Parachute radius # to = 1.2 # eta = 1 # Rdot = (6 * R * (1 - eta) / (1.2**6)) * ( @@ -1975,8 +2011,17 @@ def u_dot_parachute(self, t, u, post_processing=False): # ) # Rdot = 0 + # tf = 8 * nominal diameter / velocity at line stretch + # Calculate added mass - ma = ka * rho * (4 / 3) * np.pi * R**3 + ma = ( + self.parachute_added_mass_coefficient + * rho + * (2 / 3) + * np.pi + * self.parachute_radius**2 + * self.parachute_height + ) # Calculate freestream speed freestream_x = vx - wind_velocity_x @@ -1985,14 +2030,20 @@ def u_dot_parachute(self, t, u, post_processing=False): free_stream_speed = (freestream_x**2 + freestream_y**2 + freestream_z**2) ** 0.5 # Determine drag force - pseudo_drag = -0.5 * rho * cd_s * free_stream_speed + pseudo_drag = -0.5 * rho * self.parachute_cd_s * free_stream_speed # pseudo_drag = pseudo_drag - ka * rho * 4 * np.pi * (R**2) * Rdot - Dx = pseudo_drag * freestream_x + Dx = pseudo_drag * freestream_x # add eta efficiency for wake Dy = pseudo_drag * freestream_y Dz = pseudo_drag * freestream_z ax = Dx / (mp + ma) ay = Dy / (mp + ma) - az = (Dz - 9.8 * mp) / (mp + ma) + az = (Dz - mp * self.env.gravity.get_value_opt(z)) / (mp + ma) + + # Add coriolis acceleration + _, w_earth_y, w_earth_z = self.env.earth_rotation_vector + ax -= 2 * (vz * w_earth_y - vy * w_earth_z) + ay -= 2 * (vx * w_earth_z) + az -= 2 * (-vx * w_earth_y) if post_processing: self.__post_processed_variables.append( @@ -2070,7 +2121,7 @@ def x(self): @funcify_method("Time (s)", "Y (m)", "spline", "constant") def y(self): - """Rocket y position relative to the lauch pad as a Function of + """Rocket y position relative to the launch pad as a Function of time.""" return self.solution_array[:, [0, 2]] @@ -2966,7 +3017,7 @@ def max_rail_button2_shear_force(self): "Time (s)", "Horizontal Distance to Launch Point (m)", "spline", "constant" ) def drift(self): - """Rocket horizontal distance to tha launch point, in meters, as a + """Rocket horizontal distance to the launch point, in meters, as a Function of time.""" return np.column_stack( (self.time, (self.x[:, 1] ** 2 + self.y[:, 1] ** 2) ** 0.5) @@ -3055,14 +3106,16 @@ def __calculate_rail_button_forces(self): # TODO: complex method. null_force = Function(0) if self.out_of_rail_time_index == 0: # No rail phase, no rail button forces warnings.warn( - "Trying to calculate rail button forces without a rail phase defined." - + "The rail button forces will be set to zero." + "Trying to calculate rail button forces without a rail phase defined. " + + "The rail button forces will be set to zero.", + UserWarning, ) return null_force, null_force, null_force, null_force if len(self.rocket.rail_buttons) == 0: warnings.warn( - "Trying to calculate rail button forces without rail buttons defined." - + "The rail button forces will be set to zero." + "Trying to calculate rail button forces without rail buttons defined. " + + "The rail button forces will be set to zero.", + UserWarning, ) return null_force, null_force, null_force, null_force @@ -3174,28 +3227,6 @@ def __evaluate_post_process(self): return np.array(self.__post_processed_variables) - def post_process(self, interpolation="spline", extrapolation="natural"): - """This method is **deprecated** and is only kept here for backwards - compatibility. All attributes that need to be post processed are - computed just in time. - - Post-process all Flight information produced during - simulation. Includes the calculation of maximum values, - calculation of secondary values such as energy and conversion - of lists to Function objects to facilitate plotting. - - Returns - ------- - None - """ - # pylint: disable=unused-argument - warnings.warn( - "The method post_process is deprecated and will be removed in v1.10. " - "All attributes that need to be post processed are computed just in time.", - DeprecationWarning, - ) - self.post_processed = True - def calculate_stall_wind_velocity(self, stall_angle): # TODO: move to utilities """Function to calculate the maximum wind velocity before the angle of attack exceeds a desired angle, at the instant of departing rail launch. @@ -3235,191 +3266,53 @@ def calculate_stall_wind_velocity(self, stall_angle): # TODO: move to utilities + f" of attack exceeds {stall_angle:.3f}°: {w_v:.3f} m/s" ) - def export_pressures(self, file_name, time_step): # TODO: move out - """Exports the pressure experienced by the rocket during the flight to - an external file, the '.csv' format is recommended, as the columns will - be separated by commas. It can handle flights with or without - parachutes, although it is not possible to get a noisy pressure signal - if no parachute is added. - - If a parachute is added, the file will contain 3 columns: time in - seconds, clean pressure in Pascals and noisy pressure in Pascals. - For flights without parachutes, the third column will be discarded - - This function was created especially for the 'Projeto Jupiter' - Electronics Subsystems team and aims to help in configuring - micro-controllers. - - Parameters - ---------- - file_name : string - The final file name, - time_step : float - Time step desired for the final file - - Return - ------ - None + @deprecated( + reason="Moved to FlightDataExporter.export_pressures()", + version="v1.12.0", + alternative="rocketpy.simulation.flight_data_exporter.FlightDataExporter.export_pressures", + ) + def export_pressures(self, file_name, time_step): """ - time_points = np.arange(0, self.t_final, time_step) - # pylint: disable=W1514, E1121 - with open(file_name, "w") as file: - if len(self.rocket.parachutes) == 0: - print("No parachutes in the rocket, saving static pressure.") - for t in time_points: - file.write(f"{t:f}, {self.pressure.get_value_opt(t):.5f}\n") - else: - for parachute in self.rocket.parachutes: - for t in time_points: - p_cl = parachute.clean_pressure_signal_function.get_value_opt(t) - p_ns = parachute.noisy_pressure_signal_function.get_value_opt(t) - file.write(f"{t:f}, {p_cl:.5f}, {p_ns:.5f}\n") - # We need to save only 1 parachute data - break + .. deprecated:: 1.11 + Use :class:`rocketpy.simulation.flight_data_exporter.FlightDataExporter` + and call ``.export_pressures(...)``. + """ + return FlightDataExporter(self).export_pressures(file_name, time_step) + @deprecated( + reason="Moved to FlightDataExporter.export_data()", + version="v1.12.0", + alternative="rocketpy.simulation.flight_data_exporter.FlightDataExporter.export_data", + ) def export_data(self, file_name, *variables, time_step=None): - """Exports flight data to a comma separated value file (.csv). - - Data is exported in columns, with the first column representing time - steps. The first line of the file is a header line, specifying the - meaning of each column and its units. - - Parameters - ---------- - file_name : string - The file name or path of the exported file. Example: flight_data.csv - Do not use forbidden characters, such as / in Linux/Unix and - `<, >, :, ", /, \\, | ?, *` in Windows. - variables : strings, optional - Names of the data variables which shall be exported. Must be Flight - class attributes which are instances of the Function class. Usage - example: test_flight.export_data('test.csv', 'z', 'angle_of_attack', - 'mach_number'). - time_step : float, optional - Time step desired for the data. If None, all integration time steps - will be exported. Otherwise, linear interpolation is carried out to - calculate values at the desired time steps. Example: 0.001. """ - # TODO: we should move this method to outside of class. - - # Fast evaluation for the most basic scenario - if time_step is None and len(variables) == 0: - np.savetxt( - file_name, - self.solution, - fmt="%.6f", - delimiter=",", - header="" - "Time (s)," - "X (m)," - "Y (m)," - "Z (m)," - "E0," - "E1," - "E2," - "E3," - "W1 (rad/s)," - "W2 (rad/s)," - "W3 (rad/s)", - ) - return - - # Not so fast evaluation for general case - if variables is None: - variables = [ - "x", - "y", - "z", - "vx", - "vy", - "vz", - "e0", - "e1", - "e2", - "e3", - "w1", - "w2", - "w3", - ] - - if time_step is None: - time_points = self.time - else: - time_points = np.arange(self.t_initial, self.t_final, time_step) - - exported_matrix = [time_points] - exported_header = "Time (s)" - - # Loop through variables, get points and names (for the header) - for variable in variables: - if variable in self.__dict__: - variable_function = self.__dict__[variable] - # Deal with decorated Flight methods - else: - try: - obj = getattr(self.__class__, variable) - variable_function = obj.__get__(self, self.__class__) - except AttributeError as exc: - raise AttributeError( - f"Variable '{variable}' not found in Flight class" - ) from exc - variable_points = variable_function(time_points) - exported_matrix += [variable_points] - exported_header += f", {variable_function.__outputs__[0]}" - - exported_matrix = np.array(exported_matrix).T # Fix matrix orientation - - np.savetxt( - file_name, - exported_matrix, - fmt="%.6f", - delimiter=",", - header=exported_header, - encoding="utf-8", + .. deprecated:: 1.11 + Use :class:`rocketpy.simulation.flight_data_exporter.FlightDataExporter` + and call ``.export_data(...)``. + """ + return FlightDataExporter(self).export_data( + file_name, *variables, time_step=time_step ) + @deprecated( + reason="Moved to FlightDataExporter.export_sensor_data()", + version="v1.12.0", + alternative="rocketpy.simulation.flight_data_exporter.FlightDataExporter.export_sensor_data", + ) def export_sensor_data(self, file_name, sensor=None): - """Exports sensors data to a file. The file format can be either .csv or - .json. - - Parameters - ---------- - file_name : str - The file name or path of the exported file. Example: flight_data.csv - Do not use forbidden characters, such as / in Linux/Unix and - `<, >, :, ", /, \\, | ?, *` in Windows. - sensor : Sensor, string, optional - The sensor to export data from. Can be given as a Sensor object or - as a string with the sensor name. If None, all sensors data will be - exported. Default is None. """ - if sensor is None: - data_dict = {} - for used_sensor, measured_data in self.sensor_data.items(): - data_dict[used_sensor.name] = measured_data - else: - # export data of only that sensor - data_dict = {} - - if not isinstance(sensor, str): - data_dict[sensor.name] = self.sensor_data[sensor] - else: # sensor is a string - matching_sensors = [s for s in self.sensor_data if s.name == sensor] - - if len(matching_sensors) > 1: - data_dict[sensor] = [] - for s in matching_sensors: - data_dict[s.name].append(self.sensor_data[s]) - elif len(matching_sensors) == 1: - data_dict[sensor] = self.sensor_data[matching_sensors[0]] - else: - raise ValueError("Sensor not found in the Flight.sensor_data.") - - with open(file_name, "w") as file: - json.dump(data_dict, file) - print("Sensor data exported to: ", file_name) + .. deprecated:: 1.11 + Use :class:`rocketpy.simulation.flight_data_exporter.FlightDataExporter` + and call ``.export_sensor_data(...)``. + """ + return FlightDataExporter(self).export_sensor_data(file_name, sensor=sensor) - def export_kml( # TODO: should be moved out of this class. + @deprecated( + reason="Moved to FlightDataExporter.export_kml()", + version="v1.12.0", + alternative="rocketpy.simulation.flight_data_exporter.FlightDataExporter.export_kml", + ) + def export_kml( self, file_name="trajectory.kml", time_step=None, @@ -3427,78 +3320,18 @@ def export_kml( # TODO: should be moved out of this class. color="641400F0", altitude_mode="absolute", ): - """Exports flight data to a .kml file, which can be opened with Google - Earth to display the rocket's trajectory. - - Parameters - ---------- - file_name : string - The file name or path of the exported file. Example: flight_data.csv - time_step : float, optional - Time step desired for the data. If None, all integration time steps - will be exported. Otherwise, linear interpolation is carried out to - calculate values at the desired time steps. Example: 0.001. - extrude: bool, optional - To be used if you want to project the path over ground by using an - extruded polygon. In case False only the linestring containing the - flight path will be created. Default is True. - color : str, optional - Color of your trajectory path, need to be used in specific kml - format. Refer to http://www.zonums.com/gmaps/kml_color/ for more - info. - altitude_mode: str - Select elevation values format to be used on the kml file. Use - 'relativetoground' if you want use Above Ground Level elevation, or - 'absolute' if you want to parse elevation using Above Sea Level. - Default is 'relativetoground'. Only works properly if the ground - level is flat. Change to 'absolute' if the terrain is to irregular - or contains mountains. """ - # Define time points vector - if time_step is None: - time_points = self.time - else: - time_points = np.arange(self.t_initial, self.t_final + time_step, time_step) - # Open kml file with simplekml library - kml = simplekml.Kml(open=1) - trajectory = kml.newlinestring(name="Rocket Trajectory - Powered by RocketPy") - - if altitude_mode == "relativetoground": - # In this mode the elevation data will be the Above Ground Level - # elevation. Only works properly if the ground level is similar to - # a plane, i.e. it might not work well if the terrain has mountains - coords = [ - ( - self.longitude.get_value_opt(t), - self.latitude.get_value_opt(t), - self.altitude.get_value_opt(t), - ) - for t in time_points - ] - trajectory.coords = coords - trajectory.altitudemode = simplekml.AltitudeMode.relativetoground - else: # altitude_mode == 'absolute' - # In this case the elevation data will be the Above Sea Level elevation - # Ensure you use the correct value on self.env.elevation, otherwise - # the trajectory path can be offset from ground - coords = [ - ( - self.longitude.get_value_opt(t), - self.latitude.get_value_opt(t), - self.z.get_value_opt(t), - ) - for t in time_points - ] - trajectory.coords = coords - trajectory.altitudemode = simplekml.AltitudeMode.absolute - # Modify style of trajectory linestring - trajectory.style.linestyle.color = color - trajectory.style.polystyle.color = color - if extrude: - trajectory.extrude = 1 - # Save the KML - kml.save(file_name) - print("File ", file_name, " saved with success!") + .. deprecated:: 1.11 + Use :class:`rocketpy.simulation.flight_data_exporter.FlightDataExporter` + and call ``.export_kml(...)``. + """ + return FlightDataExporter(self).export_kml( + file_name=file_name, + time_step=time_step, + extrude=extrude, + color=color, + altitude_mode=altitude_mode, + ) def info(self): """Prints out a summary of the data available about the Flight.""" @@ -3515,7 +3348,7 @@ def time_iterator(self, node_list): yield i, node_list[i] i += 1 - def to_dict(self, include_outputs=False): + def to_dict(self, **kwargs): data = { "rocket": self.rocket, "env": self.env, @@ -3544,7 +3377,6 @@ def to_dict(self, include_outputs=False): "x_impact": self.x_impact, "y_impact": self.y_impact, "t_final": self.t_final, - "flight_phases": self.flight_phases, "function_evaluations": self.function_evaluations, "ax": self.ax, "ay": self.ay, @@ -3558,9 +3390,10 @@ def to_dict(self, include_outputs=False): "M1": self.M1, "M2": self.M2, "M3": self.M3, + "net_thrust": self.net_thrust, } - if include_outputs: + if kwargs.get("include_outputs", False): data.update( { "time": self.time, diff --git a/rocketpy/simulation/flight_data_exporter.py b/rocketpy/simulation/flight_data_exporter.py new file mode 100644 index 000000000..d73c0897a --- /dev/null +++ b/rocketpy/simulation/flight_data_exporter.py @@ -0,0 +1,298 @@ +""" +Exports a rocketpy.Flight object's data to external files. +""" + +import json + +import numpy as np +import simplekml + + +class FlightDataExporter: + """Export data from a rocketpy.Flight object to various formats.""" + + def __init__(self, flight, name="Flight Data"): + """ + Parameters + ---------- + flight : rocketpy.simulation.flight.Flight + The Flight instance to export from. + name : str, optional + A label for this exporter instance. + """ + self.name = name + self._flight = flight + + def __repr__(self): + return f"FlightDataExporter(name='{self.name}', flight='{type(self._flight).__name__}')" + + def export_pressures(self, file_name, time_step): + """Exports the pressure experienced by the rocket during the flight to + an external file, the '.csv' format is recommended, as the columns will + be separated by commas. It can handle flights with or without + parachutes, although it is not possible to get a noisy pressure signal + if no parachute is added. + + If a parachute is added, the file will contain 3 columns: time in + seconds, clean pressure in Pascals and noisy pressure in Pascals. + For flights without parachutes, the third column will be discarded + + This function was created especially for the 'Projeto Jupiter' + Electronics Subsystems team and aims to help in configuring + micro-controllers. + + Parameters + ---------- + file_name : string + The final file name, + time_step : float + Time step desired for the final file + + Return + ------ + None + """ + f = self._flight + time_points = np.arange(0, f.t_final, time_step) + # pylint: disable=W1514, E1121 + with open(file_name, "w") as file: + if len(f.rocket.parachutes) == 0: + print("No parachutes in the rocket, saving static pressure.") + for t in time_points: + file.write(f"{t:f}, {f.pressure.get_value_opt(t):.5f}\n") + else: + for parachute in f.rocket.parachutes: + for t in time_points: + p_cl = parachute.clean_pressure_signal_function.get_value_opt(t) + p_ns = parachute.noisy_pressure_signal_function.get_value_opt(t) + file.write(f"{t:f}, {p_cl:.5f}, {p_ns:.5f}\n") + # We need to save only 1 parachute data + break + + def export_data(self, file_name, *variables, time_step=None): + """Exports flight data to a comma separated value file (.csv). + + Data is exported in columns, with the first column representing time + steps. The first line of the file is a header line, specifying the + meaning of each column and its units. + + Parameters + ---------- + file_name : string + The file name or path of the exported file. Example: flight_data.csv + Do not use forbidden characters, such as / in Linux/Unix and + `<, >, :, ", /, \\, | ?, *` in Windows. + variables : strings, optional + Names of the data variables which shall be exported. Must be Flight + class attributes which are instances of the Function class. Usage + example: test_flight.export_data('test.csv', 'z', 'angle_of_attack', + 'mach_number'). + time_step : float, optional + Time step desired for the data. If None, all integration time steps + will be exported. Otherwise, linear interpolation is carried out to + calculate values at the desired time steps. Example: 0.001. + """ + f = self._flight + + # Fast evaluation for the most basic scenario + if time_step is None and len(variables) == 0: + np.savetxt( + file_name, + f.solution, + fmt="%.6f", + delimiter=",", + header="" + "Time (s)," + "X (m)," + "Y (m)," + "Z (m)," + "E0," + "E1," + "E2," + "E3," + "W1 (rad/s)," + "W2 (rad/s)," + "W3 (rad/s)", + ) + return + + # Not so fast evaluation for general case + if variables is None: + variables = [ + "x", + "y", + "z", + "vx", + "vy", + "vz", + "e0", + "e1", + "e2", + "e3", + "w1", + "w2", + "w3", + ] + + if time_step is None: + time_points = f.time + else: + time_points = np.arange(f.t_initial, f.t_final, time_step) + + exported_matrix = [time_points] + exported_header = "Time (s)" + + # Loop through variables, get points and names (for the header) + for variable in variables: + if variable in f.__dict__: + variable_function = f.__dict__[variable] + # Deal with decorated Flight methods + else: + try: + variable_function = getattr(f, variable) + except AttributeError as exc: + raise AttributeError( + f"Variable '{variable}' not found in Flight class" + ) from exc + variable_points = variable_function(time_points) + exported_matrix += [variable_points] + exported_header += f", {variable_function.__outputs__[0]}" + + exported_matrix = np.array(exported_matrix).T # Fix matrix orientation + + np.savetxt( + file_name, + exported_matrix, + fmt="%.6f", + delimiter=",", + header=exported_header, + encoding="utf-8", + ) + + def export_sensor_data(self, file_name, sensor=None): + """Exports sensors data to a file. The file format can be either .csv or + .json. + + Parameters + ---------- + file_name : str + The file name or path of the exported file. Example: flight_data.csv + Do not use forbidden characters, such as / in Linux/Unix and + `<, >, :, ", /, \\, | ?, *` in Windows. + sensor : Sensor, string, optional + The sensor to export data from. Can be given as a Sensor object or + as a string with the sensor name. If None, all sensors data will be + exported. Default is None. + """ + f = self._flight + + if sensor is None: + data_dict = {} + for used_sensor, measured_data in f.sensor_data.items(): + data_dict[used_sensor.name] = measured_data + else: + # export data of only that sensor + data_dict = {} + + if not isinstance(sensor, str): + data_dict[sensor.name] = f.sensor_data[sensor] + else: # sensor is a string + matching_sensors = [s for s in f.sensor_data if s.name == sensor] + + if len(matching_sensors) > 1: + data_dict[sensor] = [] + for s in matching_sensors: + data_dict[s.name].append(f.sensor_data[s]) + elif len(matching_sensors) == 1: + data_dict[sensor] = f.sensor_data[matching_sensors[0]] + else: + raise ValueError("Sensor not found in the Flight.sensor_data.") + + with open(file_name, "w") as file: + json.dump(data_dict, file) + print("Sensor data exported to: ", file_name) + + def export_kml( + self, + file_name="trajectory.kml", + time_step=None, + extrude=True, + color="641400F0", + altitude_mode="absolute", + ): + """Exports flight data to a .kml file, which can be opened with Google + Earth to display the rocket's trajectory. + + Parameters + ---------- + file_name : string + The file name or path of the exported file. Example: flight_data.csv + time_step : float, optional + Time step desired for the data. If None, all integration time steps + will be exported. Otherwise, linear interpolation is carried out to + calculate values at the desired time steps. Example: 0.001. + extrude: bool, optional + To be used if you want to project the path over ground by using an + extruded polygon. In case False only the linestring containing the + flight path will be created. Default is True. + color : str, optional + Color of your trajectory path, need to be used in specific kml + format. Refer to http://www.zonums.com/gmaps/kml_color/ for more + info. + altitude_mode: str + Select elevation values format to be used on the kml file. Use + 'relativetoground' if you want use Above Ground Level elevation, or + 'absolute' if you want to parse elevation using Above Sea Level. + Default is 'relativetoground'. Only works properly if the ground + level is flat. Change to 'absolute' if the terrain is to irregular + or contains mountains. + """ + f = self._flight + + # Define time points vector + if time_step is None: + time_points = f.time + else: + time_points = np.arange(f.t_initial, f.t_final + time_step, time_step) + + kml = simplekml.Kml(open=1) + trajectory = kml.newlinestring(name="Rocket Trajectory - Powered by RocketPy") + + if altitude_mode == "relativetoground": + # In this mode the elevation data will be the Above Ground Level + # elevation. Only works properly if the ground level is similar to + # a plane, i.e. it might not work well if the terrain has mountains + coords = [ + ( + f.longitude.get_value_opt(t), + f.latitude.get_value_opt(t), + f.altitude.get_value_opt(t), + ) + for t in time_points + ] + trajectory.coords = coords + trajectory.altitudemode = simplekml.AltitudeMode.relativetoground + else: # altitude_mode == 'absolute' + # In this case the elevation data will be the Above Sea Level elevation + # Ensure you use the correct value on self.env.elevation, otherwise + # the trajectory path can be offset from ground + coords = [ + ( + f.longitude.get_value_opt(t), + f.latitude.get_value_opt(t), + f.z.get_value_opt(t), + ) + for t in time_points + ] + trajectory.coords = coords + trajectory.altitudemode = simplekml.AltitudeMode.absolute + + # Modify style of trajectory linestring + trajectory.style.linestyle.color = color + trajectory.style.polystyle.color = color + if extrude: + trajectory.extrude = 1 + + # Save the KML + kml.save(file_name) + print("File ", file_name, " saved with success!") diff --git a/rocketpy/stochastic/stochastic_aero_surfaces.py b/rocketpy/stochastic/stochastic_aero_surfaces.py index 5e3750f48..07c50f8ee 100644 --- a/rocketpy/stochastic/stochastic_aero_surfaces.py +++ b/rocketpy/stochastic/stochastic_aero_surfaces.py @@ -54,6 +54,7 @@ def __init__( base_radius=None, bluffness=None, rocket_radius=None, + power=None, ): """Initializes the Stochastic Nose Cone class. @@ -84,6 +85,7 @@ def __init__( base_radius=base_radius, bluffness=bluffness, rocket_radius=rocket_radius, + power=power, name=None, ) diff --git a/rocketpy/stochastic/stochastic_parachute.py b/rocketpy/stochastic/stochastic_parachute.py index 4cf7746d7..dea8a077d 100644 --- a/rocketpy/stochastic/stochastic_parachute.py +++ b/rocketpy/stochastic/stochastic_parachute.py @@ -29,6 +29,12 @@ class StochasticParachute(StochasticModel): time-correlation). name : list[str] List with the name of the parachute object. This cannot be randomized. + radius : tuple, list, int, float + Radius of the parachute in meters. + height : tuple, list, int, float + Height of the parachute in meters. + porosity : tuple, list, int, float + Porosity of the parachute. """ def __init__( @@ -39,6 +45,9 @@ def __init__( sampling_rate=None, lag=None, noise=None, + radius=None, + height=None, + porosity=None, ): """Initializes the Stochastic Parachute class. @@ -63,6 +72,12 @@ def __init__( noise : list List of tuples in the form of (mean, standard deviation, time-correlation). + radius : tuple, list, int, float + Radius of the parachute in meters. + height : tuple, list, int, float + Height of the parachute in meters. + porosity : tuple, list, int, float + Porosity of the parachute. """ self.parachute = parachute self.cd_s = cd_s @@ -70,6 +85,9 @@ def __init__( self.sampling_rate = sampling_rate self.lag = lag self.noise = noise + self.radius = radius + self.height = height + self.porosity = porosity self._validate_trigger(trigger) self._validate_noise(noise) @@ -81,6 +99,9 @@ def __init__( lag=lag, noise=noise, name=None, + radius=radius, + height=height, + porosity=porosity, ) def _validate_trigger(self, trigger): diff --git a/rocketpy/tools.py b/rocketpy/tools.py index 692d5e224..0e80152a7 100644 --- a/rocketpy/tools.py +++ b/rocketpy/tools.py @@ -14,6 +14,7 @@ import math import re import time +import warnings from bisect import bisect_left import dill @@ -28,6 +29,66 @@ INSTALL_MAPPING = {"IPython": "ipython"} +def deprecated(reason=None, version=None, alternative=None): + """ + Decorator to mark functions or methods as deprecated. + + This decorator issues a DeprecationWarning when the decorated function + is called, indicating that it will be removed in future versions. + + Parameters + ---------- + reason : str, optional + Custom deprecation message. If not provided, a default message will be used. + version : str, optional + Version when the function will be removed. If provided, it will be + included in the warning message. + alternative : str, optional + Name of the alternative function/method that should be used instead. + If provided, it will be included in the warning message. + + Returns + ------- + callable + The decorated function with deprecation warning functionality. + + Examples + -------- + >>> @deprecated(reason="This function is obsolete", version="v2.0.0", + ... alternative="new_function") + ... def old_function(): + ... return "old result" + + >>> @deprecated() + ... def another_old_function(): + ... return "result" + """ + + def decorator(func): + @functools.wraps(func) + def wrapper(*args, **kwargs): + # Build the deprecation message + if reason: + message = reason + else: + message = f"The function `{func.__name__}` is deprecated" + + if version: + message += f" and will be removed in {version}" + + if alternative: + message += f". Use `{alternative}` instead" + + message += "." + + warnings.warn(message, DeprecationWarning, stacklevel=2) + return func(*args, **kwargs) + + return wrapper + + return decorator + + def tuple_handler(value): """Transforms the input value into a tuple that represents a range. If the input is an int or float, the output is a tuple from zero to the input @@ -121,12 +182,13 @@ def find_roots_cubic_function(a, b, c, d): Examples -------- >>> from rocketpy.tools import find_roots_cubic_function + >>> import cmath First we define the coefficients of the function ax**3 + bx**2 + cx + d >>> a, b, c, d = 1, -3, -1, 3 >>> x1, x2, x3 = find_roots_cubic_function(a, b, c, d) - >>> x1 - (-1+0j) + >>> cmath.isclose(x1, (-1+0j)) + True To get the real part of the roots, use the real attribute of the complex number. @@ -1234,6 +1296,43 @@ def from_hex_decode(obj_bytes, decoder=base64.b85decode): return dill.loads(decoder(bytes.fromhex(obj_bytes))) +def find_obj_from_hash(obj, hash_, depth_limit=None): + """Searches the object (and its children) for + an object whose '__rpy_hash' field has a particular hash value. + + Parameters + ---------- + obj : object + Object to search. + hash_ : int + Hash value to search for in the '__rpy_hash' field. + depth_limit : int, optional + Maximum depth to search recursively. If None, no limit. + + Returns + ------- + object + The object whose '__rpy_hash' matches hash_, or None if not found. + """ + + stack = [(obj, 0)] + while stack: + current_obj, current_depth = stack.pop() + if depth_limit is not None and current_depth > depth_limit: + continue + + if getattr(current_obj, "__rpy_hash", None) == hash_: + return current_obj + + if isinstance(current_obj, dict): + stack.extend((v, current_depth + 1) for v in current_obj.values()) + + elif isinstance(current_obj, (list, tuple, set)): + stack.extend((item, current_depth + 1) for item in current_obj) + + return None + + if __name__ == "__main__": # pragma: no cover import doctest diff --git a/rocketpy/utilities.py b/rocketpy/utilities.py index 3ed81df41..9ec707ffe 100644 --- a/rocketpy/utilities.py +++ b/rocketpy/utilities.py @@ -1,8 +1,6 @@ -import ast import inspect import json import os -import traceback import warnings from datetime import date from importlib.metadata import version @@ -444,96 +442,6 @@ def _flutter_prints( print(f"Altitude of minimum Safety Factor: {altitude_min_sf:.3f} m (AGL)\n") -def create_dispersion_dictionary(filename): # pragma: no cover - """Creates a dictionary with the rocket data provided by a .csv file. - File should be organized in four columns: attribute_class, parameter_name, - mean_value, standard_deviation. The first row should be the header. - It is advised to use ";" as separator, but "," should work on most of cases. - The "," separator might cause problems if the data set contains lists where - the items are separated by commas. - - Parameters - ---------- - filename : string - String with the path to the .csv file. The file should follow the - following structure: - - .. code-block:: - - attribute_class; parameter_name; mean_value; standard_deviation; - - environment; ensemble_member; [0, 1, 2, 3, 4, 5, 6, 7, 8, 9];; - - motor; impulse; 1415.15; 35.3; - - motor; burn_time; 5.274; 1; - - motor; nozzle_radius; 0.021642; 0.0005; - - motor; throat_radius; 0.008; 0.0005; - - motor; grain_separation; 0.006; 0.001; - - motor; grain_density; 1707; 50; - - Returns - ------- - dictionary - Dictionary with all rocket data to be used in dispersion analysis. The - dictionary will follow the following structure: - - .. code-block:: python - - analysis_parameters = { - 'environment': { - 'ensemble_member': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] - }, - 'motor': { - 'impulse': (1415.15, 35.3), - 'burn_time': (5.274, 1), - 'nozzle_radius': (0.021642, 0.0005), - 'throat_radius': (0.008, 0.0005), - 'grain_separation': (0.006, 0.001), - 'grain_density': (1707, 50), - } - } - """ - warnings.warn( - "This function is deprecated and will be removed in v1.10.0.", - DeprecationWarning, - ) - try: - file = np.genfromtxt( - filename, usecols=(1, 2, 3), skip_header=1, delimiter=";", dtype=str - ) - except ValueError: - warnings.warn( - "Error caught: the recommended delimiter is ';'. If using ',' " - "instead, be aware that some resources might not work as " - "expected if your data set contains lists where the items are " - "separated by commas. Please consider changing the delimiter to " - "';' if that is the case." - ) - warnings.warn(traceback.format_exc()) - file = np.genfromtxt( - filename, usecols=(1, 2, 3), skip_header=1, delimiter=",", dtype=str - ) - analysis_parameters = {} - for row in file: - if row[0] != "": - if row[2] == "": - try: - analysis_parameters[row[0].strip()] = float(row[1]) - except ValueError: - analysis_parameters[row[0].strip()] = ast.literal_eval(row[1]) - else: - try: - analysis_parameters[row[0].strip()] = (float(row[1]), float(row[2])) - except ValueError: - analysis_parameters[row[0].strip()] = "" - return analysis_parameters - - def apogee_by_mass(flight, min_mass, max_mass, points=10, plot=True): """Returns a Function object that estimates the apogee of a rocket given its mass (no motor). The function will use the rocket's mass as the diff --git a/tests/acceptance/test_bella_lui_rocket.py b/tests/acceptance/test_bella_lui_rocket.py index 76f2b4fde..a67547780 100644 --- a/tests/acceptance/test_bella_lui_rocket.py +++ b/tests/acceptance/test_bella_lui_rocket.py @@ -182,7 +182,6 @@ def drogue_trigger(p, h, y): inclination=parameters.get("inclination")[0], heading=parameters.get("heading")[0], ) - test_flight.post_process() # Comparison with Real Data flight_data = np.loadtxt( diff --git a/tests/conftest.py b/tests/conftest.py index 980e0b6ac..2f6124900 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -10,6 +10,7 @@ "tests.fixtures.motor.solid_motor_fixtures", "tests.fixtures.motor.empty_motor_fixtures", "tests.fixtures.motor.tanks_fixtures", + "tests.fixtures.motor.fluid_fixtures", "tests.fixtures.motor.tank_geometry_fixtures", "tests.fixtures.motor.generic_motor_fixtures", "tests.fixtures.parachutes.parachute_fixtures", diff --git a/tests/fixtures/motor/fluid_fixtures.py b/tests/fixtures/motor/fluid_fixtures.py new file mode 100644 index 000000000..733600873 --- /dev/null +++ b/tests/fixtures/motor/fluid_fixtures.py @@ -0,0 +1,18 @@ +import pytest + +from rocketpy import Fluid + + +@pytest.fixture +def nitrous_oxide_non_constant_fluid(): + """A nitrous_oxide fluid whose density varies with temperature + and pressure. + + Returns + ------- + rocketpy.Fluid + """ + return Fluid( + name="N2O", + density="./data/motors/liquid_motor_example/n2o_fluid_parameters.csv", + ) diff --git a/tests/fixtures/motor/tanks_fixtures.py b/tests/fixtures/motor/tanks_fixtures.py index f30ce49c6..88721966d 100644 --- a/tests/fixtures/motor/tanks_fixtures.py +++ b/tests/fixtures/motor/tanks_fixtures.py @@ -1,5 +1,8 @@ +from math import exp + import numpy as np import pytest +from scipy.constants import atm, zero_Celsius from rocketpy import ( CylindricalTank, @@ -457,3 +460,52 @@ def spherical_oxidizer_tank(oxidizer_fluid, oxidizer_pressurant): ) return oxidizer_tank + + +@pytest.fixture +def cylindrical_variable_density_oxidizer_tank(nitrous_oxide_non_constant_fluid): + """Fixture for creating a cylindrical variable density oxidizer + tank. This fixture returns a `MassBasedTank` object representing + a cylindrical oxidizer tank with variable density properties. + + Parameters + ---------- + nitrous_oxide_non_constant_fluid : rocketpy.Fluid + The fluid object representing nitrous oxide. + + Returns + ------- + MassBasedTank + The configured MassBasedTank fixture. + """ + + def liquid_mass(t): + if t < 4: + return 4.2 - 3.7 / 4 * t + else: + return 0.5 * exp(4 - t) + + def pressure(t): + if t < 4: + return (60 - 25 / 4 * t) * atm + else: + return (34 * exp(4 - t)) * atm + atm + + def temperature(t): + if t < 4: + return (22 - 5 * t) + zero_Celsius + else: + return -40 * (t - 4) + zero_Celsius + 2 + + return MassBasedTank( + name="Variable Density N2O Tank", + geometry=CylindricalTank(height=0.8, radius=0.06, spherical_caps=True), + flux_time=7, + liquid=nitrous_oxide_non_constant_fluid, + gas=nitrous_oxide_non_constant_fluid, + discretize=50, + liquid_mass=liquid_mass, + gas_mass=0, + pressure=pressure, + temperature=temperature, + ) diff --git a/tests/fixtures/surfaces/surface_fixtures.py b/tests/fixtures/surfaces/surface_fixtures.py index bf6e384c4..c48627478 100644 --- a/tests/fixtures/surfaces/surface_fixtures.py +++ b/tests/fixtures/surfaces/surface_fixtures.py @@ -69,6 +69,50 @@ def calisto_trapezoidal_fins(): ) +@pytest.fixture +def calisto_trapezoidal_fins_custom_sweep_length(): + """The trapezoidal fins of the Calisto rocket with + a custom sweep length. + + Returns + ------- + rocketpy.TrapezoidalFins + The trapezoidal fins of the Calisto rocket. + """ + return TrapezoidalFins( + n=4, + span=0.100, + root_chord=0.120, + tip_chord=0.040, + rocket_radius=0.0635, + name="calisto_trapezoidal_fins_custom_sweep_length", + cant_angle=0, + sweep_length=0.1, + ) + + +@pytest.fixture +def calisto_trapezoidal_fins_custom_sweep_angle(): + """The trapezoidal fins of the Calisto rocket with + a custom sweep angle. + + Returns + ------- + rocketpy.TrapezoidalFins + The trapezoidal fins of the Calisto rocket. + """ + return TrapezoidalFins( + n=4, + span=0.100, + root_chord=0.120, + tip_chord=0.040, + rocket_radius=0.0635, + name="calisto_trapezoidal_fins_custom_sweep_angle", + cant_angle=0, + sweep_angle=30, + ) + + @pytest.fixture def calisto_free_form_fins(): """The free form fins of the Calisto rocket. diff --git a/tests/integration/environment/__init__.py b/tests/integration/environment/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/integration/test_environment.py b/tests/integration/environment/test_environment.py similarity index 100% rename from tests/integration/test_environment.py rename to tests/integration/environment/test_environment.py diff --git a/tests/integration/test_environment_analysis.py b/tests/integration/environment/test_environment_analysis.py similarity index 99% rename from tests/integration/test_environment_analysis.py rename to tests/integration/environment/test_environment_analysis.py index e6043c85a..2b12a2057 100644 --- a/tests/integration/test_environment_analysis.py +++ b/tests/integration/environment/test_environment_analysis.py @@ -4,6 +4,7 @@ import matplotlib as plt import pytest + from rocketpy import Environment plt.rcParams.update({"figure.max_open_warning": 0}) diff --git a/tests/integration/motors/__init__.py b/tests/integration/motors/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/integration/test_empty_motor.py b/tests/integration/motors/test_empty_motor.py similarity index 100% rename from tests/integration/test_empty_motor.py rename to tests/integration/motors/test_empty_motor.py diff --git a/tests/integration/test_genericmotor.py b/tests/integration/motors/test_genericmotor.py similarity index 100% rename from tests/integration/test_genericmotor.py rename to tests/integration/motors/test_genericmotor.py diff --git a/tests/integration/test_hybridmotor.py b/tests/integration/motors/test_hybridmotor.py similarity index 100% rename from tests/integration/test_hybridmotor.py rename to tests/integration/motors/test_hybridmotor.py diff --git a/tests/integration/test_liquidmotor.py b/tests/integration/motors/test_liquidmotor.py similarity index 100% rename from tests/integration/test_liquidmotor.py rename to tests/integration/motors/test_liquidmotor.py diff --git a/tests/integration/test_tank.py b/tests/integration/motors/test_tank.py similarity index 94% rename from tests/integration/test_tank.py rename to tests/integration/motors/test_tank.py index f79d670ab..1e3c9366a 100644 --- a/tests/integration/test_tank.py +++ b/tests/integration/motors/test_tank.py @@ -22,6 +22,7 @@ "fuel_tank", "oxidizer_tank", "spherical_oxidizer_tank", + "cylindrical_variable_density_oxidizer_tank", ], ) def test_tank_all_info(mock_show, fixture_name, request): diff --git a/tests/integration/simulation/__init__.py b/tests/integration/simulation/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/integration/test_flight.py b/tests/integration/simulation/test_flight.py similarity index 70% rename from tests/integration/test_flight.py rename to tests/integration/simulation/test_flight.py index 8fddb6486..f40eb6b27 100644 --- a/tests/integration/test_flight.py +++ b/tests/integration/simulation/test_flight.py @@ -1,4 +1,3 @@ -import os from unittest.mock import patch import matplotlib as plt @@ -69,131 +68,6 @@ def test_all_info_different_solvers( assert test_flight.all_info() is None -class TestExportData: - """Tests the export_data method of the Flight class.""" - - def test_basic_export(self, flight_calisto): - """Tests basic export functionality""" - file_name = "test_export_data_1.csv" - flight_calisto.export_data(file_name) - self.validate_basic_export(flight_calisto, file_name) - os.remove(file_name) - - def test_custom_export(self, flight_calisto): - """Tests custom export functionality""" - file_name = "test_export_data_2.csv" - flight_calisto.export_data( - file_name, - "z", - "vz", - "e1", - "w3", - "angle_of_attack", - time_step=0.1, - ) - self.validate_custom_export(flight_calisto, file_name) - os.remove(file_name) - - def validate_basic_export(self, flight_calisto, file_name): - """Validates the basic export file content""" - test_data = np.loadtxt(file_name, delimiter=",") - assert np.allclose(flight_calisto.x[:, 0], test_data[:, 0], atol=1e-5) - assert np.allclose(flight_calisto.x[:, 1], test_data[:, 1], atol=1e-5) - assert np.allclose(flight_calisto.y[:, 1], test_data[:, 2], atol=1e-5) - assert np.allclose(flight_calisto.z[:, 1], test_data[:, 3], atol=1e-5) - assert np.allclose(flight_calisto.vx[:, 1], test_data[:, 4], atol=1e-5) - assert np.allclose(flight_calisto.vy[:, 1], test_data[:, 5], atol=1e-5) - assert np.allclose(flight_calisto.vz[:, 1], test_data[:, 6], atol=1e-5) - assert np.allclose(flight_calisto.e0[:, 1], test_data[:, 7], atol=1e-5) - assert np.allclose(flight_calisto.e1[:, 1], test_data[:, 8], atol=1e-5) - assert np.allclose(flight_calisto.e2[:, 1], test_data[:, 9], atol=1e-5) - assert np.allclose(flight_calisto.e3[:, 1], test_data[:, 10], atol=1e-5) - assert np.allclose(flight_calisto.w1[:, 1], test_data[:, 11], atol=1e-5) - assert np.allclose(flight_calisto.w2[:, 1], test_data[:, 12], atol=1e-5) - assert np.allclose(flight_calisto.w3[:, 1], test_data[:, 13], atol=1e-5) - - def validate_custom_export(self, flight_calisto, file_name): - """Validates the custom export file content""" - test_data = np.loadtxt(file_name, delimiter=",") - time_points = np.arange(flight_calisto.t_initial, flight_calisto.t_final, 0.1) - assert np.allclose(time_points, test_data[:, 0], atol=1e-5) - assert np.allclose(flight_calisto.z(time_points), test_data[:, 1], atol=1e-5) - assert np.allclose(flight_calisto.vz(time_points), test_data[:, 2], atol=1e-5) - assert np.allclose(flight_calisto.e1(time_points), test_data[:, 3], atol=1e-5) - assert np.allclose(flight_calisto.w3(time_points), test_data[:, 4], atol=1e-5) - assert np.allclose( - flight_calisto.angle_of_attack(time_points), test_data[:, 5], atol=1e-5 - ) - - -def test_export_kml(flight_calisto_robust): - """Tests weather the method Flight.export_kml is working as intended. - - Parameters: - ----------- - flight_calisto_robust : rocketpy.Flight - Flight object to be tested. See the conftest.py file for more info - regarding this pytest fixture. - """ - - test_flight = flight_calisto_robust - - # Basic export - test_flight.export_kml( - "test_export_data_1.kml", time_step=None, extrude=True, altitude_mode="absolute" - ) - - # Load exported files and fixtures and compare them - with open("test_export_data_1.kml", "r") as test_1: - for row in test_1: - if row[:29] == " ": - r = row[29:-15] - r = r.split(",") - for i, j in enumerate(r): - r[i] = j.split(" ") - lon, lat, z, coords = [], [], [], [] - for i in r: - for j in i: - coords.append(j) - for i in range(0, len(coords), 3): - lon.append(float(coords[i])) - lat.append(float(coords[i + 1])) - z.append(float(coords[i + 2])) - os.remove("test_export_data_1.kml") - - assert np.allclose(test_flight.latitude[:, 1], lat, atol=1e-3) - assert np.allclose(test_flight.longitude[:, 1], lon, atol=1e-3) - assert np.allclose(test_flight.z[:, 1], z, atol=1e-3) - - -def test_export_pressures(flight_calisto_robust): - """Tests if the method Flight.export_pressures is working as intended. - - Parameters - ---------- - flight_calisto_robust : Flight - Flight object to be tested. See the conftest.py file for more info - regarding this pytest fixture. - """ - file_name = "pressures.csv" - time_step = 0.5 - parachute = flight_calisto_robust.rocket.parachutes[0] - - flight_calisto_robust.export_pressures(file_name, time_step) - - with open(file_name, "r") as file: - contents = file.read() - - expected_data = "" - for t in np.arange(0, flight_calisto_robust.t_final, time_step): - p_cl = parachute.clean_pressure_signal_function(t) - p_ns = parachute.noisy_pressure_signal_function(t) - expected_data += f"{t:f}, {p_cl:.5f}, {p_ns:.5f}\n" - - assert contents == expected_data - os.remove(file_name) - - @patch("matplotlib.pyplot.show") def test_hybrid_motor_flight(mock_show, flight_calisto_hybrid_modded): # pylint: disable=unused-argument """Test the flight of a rocket with a hybrid motor. This test only validates @@ -414,6 +288,7 @@ def test_air_brakes_flight(mock_show, flight_calisto_air_brakes): # pylint: dis """ test_flight = flight_calisto_air_brakes air_brakes = test_flight.rocket.air_brakes[0] + assert air_brakes.plots.all() is None assert air_brakes.prints.all() is None @@ -491,8 +366,9 @@ def test_empty_motor_flight(mock_show, example_plain_env, calisto_motorless): # def test_freestream_speed_at_apogee(example_plain_env, calisto): """ - Asserts that a rocket at apogee has a free stream speed of 0.0 m/s in all - directions given that the environment doesn't have any wind. + Asserts that a rocket at apogee has a free stream speed of near 0.0 m/s + in all directions given that the environment doesn't have any wind. Any + speed values comes from coriolis effect. """ # NOTE: this rocket doesn't move in x or z direction. There's no wind. hard_atol = 1e-12 @@ -508,7 +384,9 @@ def test_freestream_speed_at_apogee(example_plain_env, calisto): ) assert np.isclose( - test_flight.stream_velocity_x(test_flight.apogee_time), 0.0, atol=hard_atol + test_flight.stream_velocity_x(test_flight.apogee_time), + 0.4641492104717301, + atol=hard_atol, ) assert np.isclose( test_flight.stream_velocity_y(test_flight.apogee_time), 0.0, atol=hard_atol @@ -518,9 +396,13 @@ def test_freestream_speed_at_apogee(example_plain_env, calisto): test_flight.stream_velocity_z(test_flight.apogee_time), 0.0, atol=soft_atol ) assert np.isclose( - test_flight.free_stream_speed(test_flight.apogee_time), 0.0, atol=soft_atol + test_flight.free_stream_speed(test_flight.apogee_time), + 0.4641492104717798, + atol=hard_atol, + ) + assert np.isclose( + test_flight.apogee_freestream_speed, 0.4641492104717798, atol=hard_atol ) - assert np.isclose(test_flight.apogee_freestream_speed, 0.0, atol=soft_atol) def test_rocket_csys_equivalence( @@ -546,6 +428,7 @@ def test_rocket_csys_equivalence( assert np.isclose( flight_calisto_robust.x_impact, flight_calisto_nose_to_tail_robust.x_impact, + atol=1e-3, ) assert np.isclose( flight_calisto_robust.y_impact, diff --git a/tests/integration/test_flight_data_importer.py b/tests/integration/simulation/test_flight_data_importer.py similarity index 100% rename from tests/integration/test_flight_data_importer.py rename to tests/integration/simulation/test_flight_data_importer.py diff --git a/tests/integration/simulation/test_monte_carlo.py b/tests/integration/simulation/test_monte_carlo.py new file mode 100644 index 000000000..50ebdaa5e --- /dev/null +++ b/tests/integration/simulation/test_monte_carlo.py @@ -0,0 +1,209 @@ +# pylint: disable=unused-argument +import os +from unittest.mock import patch + +import matplotlib as plt +import numpy as np +import pytest + +plt.rcParams.update({"figure.max_open_warning": 0}) + + +def _post_test_file_cleanup(): + """Clean monte carlo files after test session if they exist.""" + if os.path.exists("monte_carlo_class_example.kml"): + os.remove("monte_carlo_class_example.kml") + if os.path.exists("monte_carlo_test.errors.txt"): + os.remove("monte_carlo_test.errors.txt") + if os.path.exists("monte_carlo_test.inputs.txt"): + os.remove("monte_carlo_test.inputs.txt") + if os.path.exists("monte_carlo_test.outputs.txt"): + os.remove("monte_carlo_test.outputs.txt") + + +@pytest.mark.slow +@pytest.mark.parametrize("parallel", [False, True]) +def test_monte_carlo_simulate(monte_carlo_calisto, parallel): + """Tests the simulate method of the MonteCarlo class. + + Parameters + ---------- + monte_carlo_calisto : MonteCarlo + The MonteCarlo object, this is a pytest fixture. + """ + try: + # NOTE: this is really slow, it runs 10 flight simulations + monte_carlo_calisto.simulate( + number_of_simulations=10, append=False, parallel=parallel + ) + + assert monte_carlo_calisto.num_of_loaded_sims == 10 + assert monte_carlo_calisto.number_of_simulations == 10 + assert str(monte_carlo_calisto.filename.name) == "monte_carlo_test" + assert str(monte_carlo_calisto.error_file.name) == "monte_carlo_test.errors.txt" + assert ( + str(monte_carlo_calisto.output_file.name) == "monte_carlo_test.outputs.txt" + ) + assert np.isclose( + monte_carlo_calisto.processed_results["apogee"][0], 4711, rtol=0.2 + ) + assert np.isclose( + monte_carlo_calisto.processed_results["impact_velocity"][0], + -5.234, + rtol=0.2, + ) + finally: + _post_test_file_cleanup() + + +def test_monte_carlo_set_inputs_log(monte_carlo_calisto): + """Tests the set_inputs_log method of the MonteCarlo class. + + Parameters + ---------- + monte_carlo_calisto : MonteCarlo + The MonteCarlo object, this is a pytest fixture. + """ + try: + monte_carlo_calisto.input_file = "tests/fixtures/monte_carlo/example.inputs.txt" + monte_carlo_calisto.set_inputs_log() + assert len(monte_carlo_calisto.inputs_log) == 100 + assert all(isinstance(item, dict) for item in monte_carlo_calisto.inputs_log) + assert all( + "gravity" in item and "elevation" in item + for item in monte_carlo_calisto.inputs_log + ) + finally: + _post_test_file_cleanup() + + +def test_monte_carlo_set_outputs_log(monte_carlo_calisto): + """Tests the set_outputs_log method of the MonteCarlo class. + + Parameters + ---------- + monte_carlo_calisto : MonteCarlo + The MonteCarlo object, this is a pytest fixture. + """ + try: + monte_carlo_calisto.output_file = ( + "tests/fixtures/monte_carlo/example.outputs.txt" + ) + monte_carlo_calisto.set_outputs_log() + assert len(monte_carlo_calisto.outputs_log) == 100 + assert all(isinstance(item, dict) for item in monte_carlo_calisto.outputs_log) + assert all( + "apogee" in item and "impact_velocity" in item + for item in monte_carlo_calisto.outputs_log + ) + finally: + _post_test_file_cleanup() + + +# def test_monte_carlo_set_errors_log(monte_carlo_calisto): +# monte_carlo_calisto.error_file = "tests/fixtures/monte_carlo/example.errors.txt" +# monte_carlo_calisto.set_errors_log() +# assert + + +def test_monte_carlo_prints(monte_carlo_calisto): + """Tests the prints methods of the MonteCarlo class.""" + try: + monte_carlo_calisto.info() + monte_carlo_calisto.compare_info(monte_carlo_calisto) + finally: + _post_test_file_cleanup() + + +@patch("matplotlib.pyplot.show") # pylint: disable=unused-argument +def test_monte_carlo_plots(mock_show, monte_carlo_calisto_pre_loaded): + """Tests the plots methods of the MonteCarlo class.""" + try: + assert monte_carlo_calisto_pre_loaded.all_info() is None + assert ( + monte_carlo_calisto_pre_loaded.compare_plots(monte_carlo_calisto_pre_loaded) + is None + ) + assert ( + monte_carlo_calisto_pre_loaded.compare_ellipses( + monte_carlo_calisto_pre_loaded + ) + is None + ) + finally: + _post_test_file_cleanup() + + +def test_monte_carlo_export_ellipses_to_kml(monte_carlo_calisto_pre_loaded): + """Tests the export_ellipses_to_kml method of the MonteCarlo class. + + Parameters + ---------- + monte_carlo_calisto_pre_loaded : MonteCarlo + The MonteCarlo object, this is a pytest fixture. + """ + try: + assert ( + monte_carlo_calisto_pre_loaded.export_ellipses_to_kml( + filename="monte_carlo_class_example.kml", + origin_lat=32.990254, + origin_lon=-106.974998, + type="all", + ) + is None + ) + finally: + _post_test_file_cleanup() + + +@pytest.mark.slow +def test_monte_carlo_callback(monte_carlo_calisto): + """Tests the data_collector argument of the MonteCarlo class. + + Parameters + ---------- + monte_carlo_calisto : MonteCarlo + The MonteCarlo object, this is a pytest fixture. + """ + try: + # define valid data collector + valid_data_collector = { + "name": lambda flight: flight.name, + "density_t0": lambda flight: flight.env.density(0), + } + + monte_carlo_calisto.data_collector = valid_data_collector + # NOTE: this is really slow, it runs 10 flight simulations + monte_carlo_calisto.simulate(number_of_simulations=10, append=False) + + # tests if print works when we have None in summary + monte_carlo_calisto.info() + + ## tests if an error is raised for invalid data_collector definitions + # invalid type + def invalid_data_collector(flight): + return flight.name + + with pytest.raises(ValueError): + monte_carlo_calisto._check_data_collector(invalid_data_collector) + + # invalid key overwrite + invalid_data_collector = {"apogee": lambda flight: flight.apogee} + with pytest.raises(ValueError): + monte_carlo_calisto._check_data_collector(invalid_data_collector) + + # invalid callback definition + invalid_data_collector = {"name": "Calisto"} # callbacks must be callables! + with pytest.raises(ValueError): + monte_carlo_calisto._check_data_collector(invalid_data_collector) + + # invalid logic (division by zero) + invalid_data_collector = { + "density_t0": lambda flight: flight.env.density(0) / "0", + } + monte_carlo_calisto.data_collector = invalid_data_collector + # NOTE: this is really slow, it runs 10 flight simulations + with pytest.raises(ValueError): + monte_carlo_calisto.simulate(number_of_simulations=10, append=False) + finally: + _post_test_file_cleanup() diff --git a/tests/integration/test_multivariate_rejection_sampler.py b/tests/integration/simulation/test_multivariate_rejection_sampler.py similarity index 100% rename from tests/integration/test_multivariate_rejection_sampler.py rename to tests/integration/simulation/test_multivariate_rejection_sampler.py diff --git a/tests/integration/test_encoding.py b/tests/integration/test_encoding.py index 3ef71252c..38d5aef41 100644 --- a/tests/integration/test_encoding.py +++ b/tests/integration/test_encoding.py @@ -1,12 +1,16 @@ import json import os +from unittest.mock import patch import numpy as np import pytest from rocketpy._encoders import RocketPyDecoder, RocketPyEncoder +from rocketpy.tools import from_hex_decode +# pylint: disable=unused-argument +@patch("matplotlib.pyplot.show") @pytest.mark.parametrize( ["flight_name", "include_outputs"], [ @@ -15,9 +19,13 @@ ("flight_calisto_robust", True), ("flight_calisto_liquid_modded", False), ("flight_calisto_hybrid_modded", False), + ("flight_calisto_air_brakes", False), + ("flight_calisto_with_sensors", False), ], ) -def test_flight_save_load_no_resimulate(flight_name, include_outputs, request): +def test_flight_save_load_no_resimulate( + mock_show, flight_name, include_outputs, request +): """Test encoding a ``rocketpy.Flight``. Parameters @@ -48,6 +56,8 @@ def test_flight_save_load_no_resimulate(flight_name, include_outputs, request): # Higher tolerance due to random parachute trigger assert np.isclose(flight_to_save.t_final, flight_loaded.t_final, rtol=1e-3) + flight_loaded.all_info() + os.remove("flight.json") @@ -60,6 +70,8 @@ def test_flight_save_load_no_resimulate(flight_name, include_outputs, request): ("flight_calisto_robust", True), ("flight_calisto_liquid_modded", False), ("flight_calisto_hybrid_modded", False), + ("flight_calisto_air_brakes", False), + ("flight_calisto_with_sensors", False), ], ) def test_flight_save_load_resimulate(flight_name, include_outputs, request): @@ -91,7 +103,7 @@ def test_flight_save_load_resimulate(flight_name, include_outputs, request): assert np.isclose(flight_to_save.apogee_time, flight_loaded.apogee_time) # Higher tolerance due to random parachute trigger - assert np.isclose(flight_to_save.t_final, flight_loaded.t_final, rtol=1e-3) + assert np.isclose(flight_to_save.t_final, flight_loaded.t_final, rtol=5e-3) os.remove("flight.json") @@ -228,3 +240,90 @@ def test_rocket_encoder(rocket_name, request): rocket_to_encode.static_margin(sample_times), rocket_loaded.static_margin(sample_times), ) + + +@pytest.mark.parametrize( + "fin_name", + [ + "calisto_trapezoidal_fins", + "calisto_trapezoidal_fins_custom_sweep_length", + "calisto_trapezoidal_fins_custom_sweep_angle", + ], +) +def test_trapezoidal_fins_encoder(fin_name, request): + """Test encoding a ``rocketpy.TrapezoidalFins``. + + Parameters + ---------- + fin_name : str + Name of the fin fixture to encode. + request : pytest.FixtureRequest + Pytest request object. + """ + fin_to_encode = request.getfixturevalue(fin_name) + + json_encoded = json.dumps(fin_to_encode, cls=RocketPyEncoder) + + fin_loaded = json.loads(json_encoded, cls=RocketPyDecoder) + + assert isinstance(fin_loaded, type(fin_to_encode)) + assert fin_to_encode.n == fin_loaded.n + assert np.isclose(fin_to_encode.span, fin_loaded.span) + assert np.isclose(fin_to_encode.root_chord, fin_loaded.root_chord) + assert np.isclose(fin_to_encode.tip_chord, fin_loaded.tip_chord) + assert np.isclose(fin_to_encode.rocket_radius, fin_loaded.rocket_radius) + assert np.isclose(fin_to_encode.sweep_length, fin_loaded.sweep_length) + if fin_to_encode._sweep_angle is not None and fin_loaded._sweep_angle is not None: + assert np.isclose(fin_to_encode.sweep_angle, fin_loaded.sweep_angle) + + +@pytest.mark.parametrize("rocket_name", ["calisto_robust", "calisto_hybrid_modded"]) +def test_encoder_discretize(rocket_name, request): + """Test encoding the total mass of ``rocketpy.Rocket`` with + discretized encoding. + + Parameters + ---------- + rocket_name : str + Name of the rocket fixture to encode. + request : pytest.FixtureRequest + Pytest request object. + """ + rocket_to_encode = request.getfixturevalue(rocket_name) + + json_encoded = json.dumps( + rocket_to_encode, cls=RocketPyEncoder, discretize=True, include_outputs=True + ) + + mass_loaded = json.loads( + json.dumps(json.loads(json_encoded)["total_mass"]), cls=RocketPyDecoder + ) + + sample_times = np.linspace(*rocket_to_encode.motor.burn_time, 100) + + np.testing.assert_allclose( + mass_loaded(sample_times), + rocket_to_encode.total_mass(sample_times), + rtol=1e-3, + atol=1e-1, + ) + assert isinstance(mass_loaded.source, np.ndarray) + + +@pytest.mark.parametrize("parachute_name", ["calisto_main_chute"]) +def test_encoder_no_pickle(parachute_name, request): + """Test encoding of a ``rocketpy.Parachute`` disallowing + pickle usage. + """ + parachute_to_encode = request.getfixturevalue(parachute_name) + + json_encoded = json.dumps( + parachute_to_encode, + cls=RocketPyEncoder, + allow_pickle=False, + ) + + trigger_loaded = json.loads(json_encoded)["trigger"] + + with pytest.raises(ValueError, match=r"non-hexadecimal number found"): + from_hex_decode(trigger_loaded) diff --git a/tests/integration/test_monte_carlo.py b/tests/integration/test_monte_carlo.py deleted file mode 100644 index 5da5066a2..000000000 --- a/tests/integration/test_monte_carlo.py +++ /dev/null @@ -1,180 +0,0 @@ -# pylint: disable=unused-argument -import os -from unittest.mock import patch - -import matplotlib as plt -import numpy as np -import pytest - -plt.rcParams.update({"figure.max_open_warning": 0}) - - -@pytest.mark.slow -@pytest.mark.parametrize("parallel", [False, True]) -def test_monte_carlo_simulate(monte_carlo_calisto, parallel): - """Tests the simulate method of the MonteCarlo class. - - Parameters - ---------- - monte_carlo_calisto : MonteCarlo - The MonteCarlo object, this is a pytest fixture. - """ - # NOTE: this is really slow, it runs 10 flight simulations - monte_carlo_calisto.simulate( - number_of_simulations=10, append=False, parallel=parallel - ) - - assert monte_carlo_calisto.num_of_loaded_sims == 10 - assert monte_carlo_calisto.number_of_simulations == 10 - assert str(monte_carlo_calisto.filename.name) == "monte_carlo_test" - assert str(monte_carlo_calisto.error_file.name) == "monte_carlo_test.errors.txt" - assert str(monte_carlo_calisto.output_file.name) == "monte_carlo_test.outputs.txt" - assert np.isclose( - monte_carlo_calisto.processed_results["apogee"][0], 4711, rtol=0.2 - ) - assert np.isclose( - monte_carlo_calisto.processed_results["impact_velocity"][0], - -5.234, - rtol=0.2, - ) - os.remove("monte_carlo_test.errors.txt") - os.remove("monte_carlo_test.outputs.txt") - os.remove("monte_carlo_test.inputs.txt") - - -def test_monte_carlo_set_inputs_log(monte_carlo_calisto): - """Tests the set_inputs_log method of the MonteCarlo class. - - Parameters - ---------- - monte_carlo_calisto : MonteCarlo - The MonteCarlo object, this is a pytest fixture. - """ - monte_carlo_calisto.input_file = "tests/fixtures/monte_carlo/example.inputs.txt" - monte_carlo_calisto.set_inputs_log() - assert len(monte_carlo_calisto.inputs_log) == 100 - assert all(isinstance(item, dict) for item in monte_carlo_calisto.inputs_log) - assert all( - "gravity" in item and "elevation" in item - for item in monte_carlo_calisto.inputs_log - ) - - -def test_monte_carlo_set_outputs_log(monte_carlo_calisto): - """Tests the set_outputs_log method of the MonteCarlo class. - - Parameters - ---------- - monte_carlo_calisto : MonteCarlo - The MonteCarlo object, this is a pytest fixture. - """ - monte_carlo_calisto.output_file = "tests/fixtures/monte_carlo/example.outputs.txt" - monte_carlo_calisto.set_outputs_log() - assert len(monte_carlo_calisto.outputs_log) == 100 - assert all(isinstance(item, dict) for item in monte_carlo_calisto.outputs_log) - assert all( - "apogee" in item and "impact_velocity" in item - for item in monte_carlo_calisto.outputs_log - ) - - -# def test_monte_carlo_set_errors_log(monte_carlo_calisto): -# monte_carlo_calisto.error_file = "tests/fixtures/monte_carlo/example.errors.txt" -# monte_carlo_calisto.set_errors_log() -# assert - - -def test_monte_carlo_prints(monte_carlo_calisto): - """Tests the prints methods of the MonteCarlo class.""" - monte_carlo_calisto.info() - monte_carlo_calisto.compare_info(monte_carlo_calisto) - - -@patch("matplotlib.pyplot.show") # pylint: disable=unused-argument -def test_monte_carlo_plots(mock_show, monte_carlo_calisto_pre_loaded): - """Tests the plots methods of the MonteCarlo class.""" - assert monte_carlo_calisto_pre_loaded.all_info() is None - assert ( - monte_carlo_calisto_pre_loaded.compare_plots(monte_carlo_calisto_pre_loaded) - is None - ) - assert ( - monte_carlo_calisto_pre_loaded.compare_ellipses(monte_carlo_calisto_pre_loaded) - is None - ) - - -def test_monte_carlo_export_ellipses_to_kml(monte_carlo_calisto_pre_loaded): - """Tests the export_ellipses_to_kml method of the MonteCarlo class. - - Parameters - ---------- - monte_carlo_calisto_pre_loaded : MonteCarlo - The MonteCarlo object, this is a pytest fixture. - """ - assert ( - monte_carlo_calisto_pre_loaded.export_ellipses_to_kml( - filename="monte_carlo_class_example.kml", - origin_lat=32.990254, - origin_lon=-106.974998, - type="all", - ) - is None - ) - - os.remove("monte_carlo_class_example.kml") - - -@pytest.mark.slow -def test_monte_carlo_callback(monte_carlo_calisto): - """Tests the data_collector argument of the MonteCarlo class. - - Parameters - ---------- - monte_carlo_calisto : MonteCarlo - The MonteCarlo object, this is a pytest fixture. - """ - - # define valid data collector - valid_data_collector = { - "name": lambda flight: flight.name, - "density_t0": lambda flight: flight.env.density(0), - } - - monte_carlo_calisto.data_collector = valid_data_collector - # NOTE: this is really slow, it runs 10 flight simulations - monte_carlo_calisto.simulate(number_of_simulations=10, append=False) - - # tests if print works when we have None in summary - monte_carlo_calisto.info() - - ## tests if an error is raised for invalid data_collector definitions - # invalid type - def invalid_data_collector(flight): - return flight.name - - with pytest.raises(ValueError): - monte_carlo_calisto._check_data_collector(invalid_data_collector) - - # invalid key overwrite - invalid_data_collector = {"apogee": lambda flight: flight.apogee} - with pytest.raises(ValueError): - monte_carlo_calisto._check_data_collector(invalid_data_collector) - - # invalid callback definition - invalid_data_collector = {"name": "Calisto"} # callbacks must be callables! - with pytest.raises(ValueError): - monte_carlo_calisto._check_data_collector(invalid_data_collector) - - # invalid logic (division by zero) - invalid_data_collector = { - "density_t0": lambda flight: flight.env.density(0) / "0", - } - monte_carlo_calisto.data_collector = invalid_data_collector - # NOTE: this is really slow, it runs 10 flight simulations - with pytest.raises(ValueError): - monte_carlo_calisto.simulate(number_of_simulations=10, append=False) - - os.remove("monte_carlo_test.errors.txt") - os.remove("monte_carlo_test.outputs.txt") - os.remove("monte_carlo_test.inputs.txt") diff --git a/tests/unit/environment/__init__.py b/tests/unit/environment/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/unit/test_environment.py b/tests/unit/environment/test_environment.py similarity index 98% rename from tests/unit/test_environment.py rename to tests/unit/environment/test_environment.py index 714470a40..6ad3e51db 100644 --- a/tests/unit/test_environment.py +++ b/tests/unit/environment/test_environment.py @@ -6,6 +6,7 @@ import pytz from rocketpy import Environment +from rocketpy.environment.tools import geodesic_to_utm, utm_to_geodesic @pytest.mark.parametrize( @@ -78,7 +79,7 @@ def test_geodesic_coordinate_geodesic_to_utm_converts_coordinate(): utm_letter, north_south_hemis, east_west_hemis, - ) = Environment.geodesic_to_utm( + ) = geodesic_to_utm( lat=32.990254, lon=-106.974998, semi_major_axis=6378137.0, # WGS84 @@ -98,7 +99,7 @@ class and checks the conversion results from UTM to geodesic coordinates. """ - lat, lon = Environment.utm_to_geodesic( + lat, lon = utm_to_geodesic( x=315468.64, y=3651938.65, utm_zone=13, diff --git a/tests/unit/test_environment_analysis.py b/tests/unit/environment/test_environment_analysis.py similarity index 100% rename from tests/unit/test_environment_analysis.py rename to tests/unit/environment/test_environment_analysis.py diff --git a/tests/unit/mathutils/__init__.py b/tests/unit/mathutils/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/unit/test_function.py b/tests/unit/mathutils/test_function.py similarity index 92% rename from tests/unit/test_function.py rename to tests/unit/mathutils/test_function.py index 77f5916f4..96acf45f5 100644 --- a/tests/unit/test_function.py +++ b/tests/unit/mathutils/test_function.py @@ -250,6 +250,113 @@ def test_set_discrete_based_on_model_non_mutator(linear_func): assert callable(func.source) +source_array = np.array( + [ + [-2, -4, -6], + [-0.75, -1.5, -2.25], + [0, 0, 0], + [0, 1, 1], + [0.5, 1, 1.5], + [1.5, 1, 2.5], + [2, 4, 6], + ] +) +cropped_array = np.array([[-0.75, -1.5, -2.25], [0, 0, 0], [0, 1, 1], [0.5, 1, 1.5]]) +clipped_array = np.array([[0, 0, 0], [0, 1, 1], [0.5, 1, 1.5]]) + + +@pytest.mark.parametrize( + "array3dsource, array3dcropped", + [ + (source_array, cropped_array), + ], +) +def test_crop_ndarray(array3dsource, array3dcropped): # pylint: disable=unused-argument + """Tests the functionality of crop method of the Function class. + The source is initialized as a ndarray before cropping. + """ + func = Function(array3dsource, inputs=["x1", "x2"], outputs="y") + cropped_func = func.crop([(-1, 1), (-2, 2)]) + + assert isinstance(func, Function) + assert isinstance(cropped_func, Function) + assert np.array_equal(cropped_func.source, array3dcropped) + assert isinstance(cropped_func.source, type(func.source)) + + +def test_crop_function(): + """Tests the functionality of crop method of the Function class. + The source is initialized as a function before cropping. + """ + func = Function( + lambda x1, x2: np.sin(x1) * np.cos(x2), inputs=["x1", "x2"], outputs="y" + ) + cropped_func = func.crop([(-1, 1), (-2, 2)]) + + assert isinstance(func, Function) + assert isinstance(cropped_func, Function) + assert callable(func.source) + assert callable(cropped_func.source) + + +def test_crop_constant(): + """Tests the functionality of crop method of the Function class. + The source is initialized as a single integer constant before cropping. + """ + func = Function(13) + cropped_func = func.crop([(-1, 1)]) + + assert isinstance(func, Function) + assert isinstance(cropped_func, Function) + assert callable(func.source) + assert callable(cropped_func.source) + + +@pytest.mark.parametrize( + "array3dsource, array3dclipped", + [ + (source_array, clipped_array), + ], +) +def test_clip_ndarray(array3dsource, array3dclipped): # pylint: disable=unused-argument + """Tests the functionality of clip method of the Function class. + The source is initialized as a ndarray before clipping. + """ + func = Function(array3dsource, inputs=["x1", "x2"], outputs="y") + clipped_func = func.clip([(-2, 2)]) + + assert isinstance(func, Function) + assert isinstance(clipped_func, Function) + assert np.array_equal(clipped_func.source, array3dclipped) + assert isinstance(clipped_func.source, type(func.source)) + + +def test_clip_function(): + """Tests the functionality of clip method of the Function class. + The source is initialized as a function before clipping. + """ + func = Function(lambda x: x**2, inputs="x", outputs="y") + clipped_func = func.clip([(-1, 1)]) + + assert isinstance(func, Function) + assert isinstance(clipped_func, Function) + assert callable(func.source) + assert callable(clipped_func.source) + + +def test_clip_constant(): + """Tests the functionality of clip method of the Function class. + The source is initialized as a single integer constant before clipping. + """ + func = Function(1) + clipped_func = func.clip([(-2, 2)]) + + assert isinstance(func, Function) + assert isinstance(clipped_func, Function) + assert callable(func.source) + assert callable(clipped_func.source) + + @pytest.mark.parametrize( "x, y, expected_x, expected_y", [ diff --git a/tests/unit/test_piecewise_function.py b/tests/unit/mathutils/test_piecewise_function.py similarity index 100% rename from tests/unit/test_piecewise_function.py rename to tests/unit/mathutils/test_piecewise_function.py diff --git a/tests/unit/test_tools_matrix.py b/tests/unit/mathutils/test_tools_matrix.py similarity index 100% rename from tests/unit/test_tools_matrix.py rename to tests/unit/mathutils/test_tools_matrix.py diff --git a/tests/unit/test_tools_vector.py b/tests/unit/mathutils/test_tools_vector.py similarity index 100% rename from tests/unit/test_tools_vector.py rename to tests/unit/mathutils/test_tools_vector.py diff --git a/tests/unit/motors/__init__.py b/tests/unit/motors/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/unit/motors/test_fluid.py b/tests/unit/motors/test_fluid.py new file mode 100644 index 000000000..1c2cbff4b --- /dev/null +++ b/tests/unit/motors/test_fluid.py @@ -0,0 +1,57 @@ +import numpy as np +from scipy.constants import atm, zero_Celsius + +from rocketpy.mathutils.function import Function +from rocketpy.motors.fluid import Fluid + + +def test_constant_density_fluid_properties(): + """Test Fluid with constant density has correct properties.""" + water = Fluid(name="Water", density=1000.0) + assert water.name == "Water" + assert water.density == 1000.0 + assert isinstance(water.density_function, Function) + assert np.isclose(water.density_function.get_value(zero_Celsius, atm), 1000.0) + + +def test_variable_density_fluid_properties(): + """Test Fluid with variable density function.""" + test_fluid = Fluid(name="TestFluid", density=lambda t, p: 44 * p / (8.314 * t)) + assert test_fluid.name == "TestFluid" + assert callable(test_fluid.density) + assert isinstance(test_fluid.density, Function) + assert isinstance(test_fluid.density_function, Function) + expected_density = 44 * atm / (8.314 * zero_Celsius) + assert np.isclose(test_fluid.density_function(zero_Celsius, atm), expected_density) + + +def test_get_time_variable_density_with_callable_sources(): + """Test get_time_variable_density with callable temperature and pressure.""" + test_fluid = Fluid("TestFluid", lambda t, p: t + p) + temperature = Function(lambda t: 300 + t) + pressure = Function(lambda t: 100_000 + 10 * t) + density_time_func = test_fluid.get_time_variable_density(temperature, pressure) + assert np.isclose(density_time_func.get_value(0), 300 + 100_000) + assert np.isclose(density_time_func.get_value(10), 310 + 100_100) + + +def test_get_time_variable_density_with_arrays(): + """Test get_time_variable_density with array sources.""" + fluid = Fluid("TestFluid", lambda t, p: t + p) + times = np.array([0, 10]) + temperature = Function(np.column_stack((times, [300, 310]))) + pressure = Function(np.column_stack((times, [100_000, 100_100]))) + density_time_func = fluid.get_time_variable_density(temperature, pressure) + expected = [300 + 100_000, 310 + 100_100] + np.testing.assert_allclose(density_time_func(times), expected) + + +def test_get_time_variable_density_with_callable_and_arrays(): + """Test get_time_variable_density with mixed sources.""" + fluid = Fluid("TestFluid", lambda t, p: t + p) + times = np.array([0, 10]) + temperature = Function(np.column_stack((times, [300, 310]))) + pressure = Function(lambda t: 100_000 + 10 * t) + density_time_func = fluid.get_time_variable_density(temperature, pressure) + expected = [300 + 100_000, 310 + 100_100] + np.testing.assert_allclose(density_time_func(times), expected) diff --git a/tests/unit/test_genericmotor.py b/tests/unit/motors/test_genericmotor.py similarity index 100% rename from tests/unit/test_genericmotor.py rename to tests/unit/motors/test_genericmotor.py diff --git a/tests/unit/test_hybridmotor.py b/tests/unit/motors/test_hybridmotor.py similarity index 100% rename from tests/unit/test_hybridmotor.py rename to tests/unit/motors/test_hybridmotor.py diff --git a/tests/unit/test_liquidmotor.py b/tests/unit/motors/test_liquidmotor.py similarity index 100% rename from tests/unit/test_liquidmotor.py rename to tests/unit/motors/test_liquidmotor.py diff --git a/tests/unit/test_solidmotor.py b/tests/unit/motors/test_solidmotor.py similarity index 100% rename from tests/unit/test_solidmotor.py rename to tests/unit/motors/test_solidmotor.py diff --git a/tests/unit/test_tank.py b/tests/unit/motors/test_tank.py similarity index 91% rename from tests/unit/test_tank.py rename to tests/unit/motors/test_tank.py index 95179ccfa..87059be0b 100644 --- a/tests/unit/test_tank.py +++ b/tests/unit/motors/test_tank.py @@ -1,13 +1,10 @@ from math import isclose from pathlib import Path -from unittest.mock import patch import numpy as np import pytest import scipy.integrate as spi -from rocketpy.motors import TankGeometry - BASE_PATH = Path("./data/rockets/berkeley/") @@ -360,6 +357,34 @@ def expected_gas_inertia(t): ) -@patch("matplotlib.pyplot.show") -def test_tank_geometry_plots_info(mock_show): # pylint: disable=unused-argument - assert TankGeometry({(0, 5): 1}).plots.all() is None +def test_variable_density_mass_tank(cylindrical_variable_density_oxidizer_tank): + """Tests a cylindrical tank with variable density fluids + from its temperature and pressure values. + + Parameters + ---------- + cylindrical_variable_density_oxidizer_tank: MassBasedTank + The tank to be tested. + """ + tank = cylindrical_variable_density_oxidizer_tank + time_steps = np.linspace(*tank.flux_time, 75) + + assert (tank._liquid_density(time_steps) > 0).all() + assert (tank._gas_density(time_steps) > 0).all() + assert (tank._liquid_density(time_steps) < 1e5).all() + assert (tank._gas_density(time_steps) < 1e5).all() + np.testing.assert_allclose( + tank.liquid_mass(time_steps), + tank.liquid_volume(time_steps) * tank._liquid_density(time_steps), + atol=1e-2, + ) + np.testing.assert_allclose( + tank.gas_mass(time_steps), + tank.gas_volume(time_steps) * tank._gas_density(time_steps), + atol=1e-2, + ) + np.testing.assert_allclose( + tank.gas_mass(time_steps), + 0, + atol=1e-4, + ) diff --git a/tests/unit/test_tank_geometry.py b/tests/unit/motors/test_tank_geometry.py similarity index 93% rename from tests/unit/test_tank_geometry.py rename to tests/unit/motors/test_tank_geometry.py index 54e3671c2..ff4a525ba 100644 --- a/tests/unit/test_tank_geometry.py +++ b/tests/unit/motors/test_tank_geometry.py @@ -1,8 +1,11 @@ from pathlib import Path +from unittest.mock import patch import numpy as np import pytest +from rocketpy.motors import TankGeometry + PRESSURANT_PARAMS = (0.135 / 2, 0.981) PROPELLANT_PARAMS = (0.0744, 0.8068) SPHERICAL_PARAMS = (0.05, 0.1) @@ -124,3 +127,8 @@ def test_tank_inertia(params, request): rtol=1e-5, atol=1e-9, ) + + +@patch("matplotlib.pyplot.show") +def test_tank_geometry_plots_info(mock_show): # pylint: disable=unused-argument + assert TankGeometry({(0, 5): 1}).plots.all() is None diff --git a/tests/unit/rocket/__init__.py b/tests/unit/rocket/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/unit/rocket/aero_surface/__init__.py b/tests/unit/rocket/aero_surface/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/unit/test_aero_surfaces.py b/tests/unit/rocket/aero_surface/test_aero_surfaces.py similarity index 100% rename from tests/unit/test_aero_surfaces.py rename to tests/unit/rocket/aero_surface/test_aero_surfaces.py diff --git a/tests/unit/test_generic_surfaces.py b/tests/unit/rocket/aero_surface/test_generic_surfaces.py similarity index 100% rename from tests/unit/test_generic_surfaces.py rename to tests/unit/rocket/aero_surface/test_generic_surfaces.py diff --git a/tests/unit/test_linear_generic_surfaces.py b/tests/unit/rocket/aero_surface/test_linear_generic_surfaces.py similarity index 100% rename from tests/unit/test_linear_generic_surfaces.py rename to tests/unit/rocket/aero_surface/test_linear_generic_surfaces.py diff --git a/tests/unit/test_rocket.py b/tests/unit/rocket/test_rocket.py similarity index 100% rename from tests/unit/test_rocket.py rename to tests/unit/rocket/test_rocket.py diff --git a/tests/unit/sensors/__init__.py b/tests/unit/sensors/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/unit/test_sensor.py b/tests/unit/sensors/test_sensor.py similarity index 100% rename from tests/unit/test_sensor.py rename to tests/unit/sensors/test_sensor.py diff --git a/tests/unit/simulation/__init__.py b/tests/unit/simulation/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/unit/test_flight.py b/tests/unit/simulation/test_flight.py similarity index 97% rename from tests/unit/test_flight.py rename to tests/unit/simulation/test_flight.py index 260a2b138..1bd9384e2 100644 --- a/tests/unit/test_flight.py +++ b/tests/unit/simulation/test_flight.py @@ -109,15 +109,15 @@ def test_get_solution_at_time(flight_calisto): flight_calisto.get_solution_at_time(flight_calisto.t_final), np.array( [ - 48.4313533, - 0.0, - 985.7665845, - -0.00000229951048, - 0.0, - 11.2223284, - -341.028803, - 0.999048222, - -0.0436193874, + 48.43719482805657, + -14.836008075478597, + 985.9858934483618, + -3.4415459237894554e-05, + 0.0007572309307800201, + 11.21695000766671, + -341.1460775169661, + 0.9990482215818578, + -0.043619387365336, 0.0, 0.0, 0.0, @@ -212,7 +212,7 @@ def test_export_sensor_data(flight_calisto_with_sensors): [ ("t_initial", (0.25886, -0.649623, 0)), ("out_of_rail_time", (0.792028, -1.987634, 0)), - ("apogee_time", (-0.522875, -0.741825, 0)), + ("apogee_time", (-0.509420, -0.732933, -2.089120e-14)), ("t_final", (0, 0, 0)), ], ) @@ -251,8 +251,8 @@ def test_aerodynamic_moments(flight_calisto_custom_wind, flight_time, expected_v [ ("t_initial", (1.654150, 0.659142, -0.067103)), ("out_of_rail_time", (5.052628, 2.013361, -1.75370)), - ("apogee_time", (2.339424, -1.648934, -0.938867)), - ("t_final", (0, 0, 159.2210)), + ("apogee_time", (2.321838, -1.613641, -0.962108)), + ("t_final", (-0.019802, 0.012030, 159.051604)), ], ) def test_aerodynamic_forces(flight_calisto_custom_wind, flight_time, expected_values): @@ -292,7 +292,7 @@ def test_aerodynamic_forces(flight_calisto_custom_wind, flight_time, expected_va ("out_of_rail_time", (0, 2.248540, 25.700928)), ( "apogee_time", - (-14.488364, 15.638049, -0.000191), + (-14.826350, 15.670022, -0.000264), ), ("t_final", (5, 2, -5.660155)), ], @@ -540,7 +540,7 @@ def test_lat_lon_conversion_from_origin(mock_show, example_plain_env, calisto_ro heading=0, ) - assert abs(test_flight.longitude(test_flight.t_final) - 0) < 1e-12 + assert abs(test_flight.longitude(test_flight.t_final)) < 1e-4 assert test_flight.latitude(test_flight.t_final) > 0 diff --git a/tests/unit/test_flight_time_nodes.py b/tests/unit/simulation/test_flight_time_nodes.py similarity index 100% rename from tests/unit/test_flight_time_nodes.py rename to tests/unit/simulation/test_flight_time_nodes.py diff --git a/tests/unit/test_monte_carlo.py b/tests/unit/simulation/test_monte_carlo.py similarity index 100% rename from tests/unit/test_monte_carlo.py rename to tests/unit/simulation/test_monte_carlo.py diff --git a/tests/unit/test_multivariate_rejection_sampler.py b/tests/unit/simulation/test_multivariate_rejection_sampler.py similarity index 100% rename from tests/unit/test_multivariate_rejection_sampler.py rename to tests/unit/simulation/test_multivariate_rejection_sampler.py diff --git a/tests/unit/test_stochastic_model.py b/tests/unit/stochastic/test_stochastic_model.py similarity index 100% rename from tests/unit/test_stochastic_model.py rename to tests/unit/stochastic/test_stochastic_model.py diff --git a/tests/unit/test_flight_data_exporter.py b/tests/unit/test_flight_data_exporter.py new file mode 100644 index 000000000..b7f8f4e9d --- /dev/null +++ b/tests/unit/test_flight_data_exporter.py @@ -0,0 +1,194 @@ +"""Unit tests for FlightDataExporter class. + +This module tests the data export functionality of the FlightDataExporter class, +which exports flight simulation data to various formats (CSV, JSON, KML). +""" + +import json + +import numpy as np + +from rocketpy.simulation import FlightDataExporter + + +def test_export_pressures_writes_csv_rows(flight_calisto_robust, tmp_path): + """Test that export_pressures writes CSV rows with pressure data. + + Validates that the exported file contains multiple data rows in CSV format + with 2-3 columns (time and pressure values). + + Parameters + ---------- + flight_calisto_robust : rocketpy.Flight + Flight object with parachutes configured. + tmp_path : pathlib.Path + Pytest fixture for temporary directories. + """ + out = tmp_path / "pressures.csv" + FlightDataExporter(flight_calisto_robust).export_pressures(str(out), time_step=0.2) + lines = out.read_text(encoding="utf-8").strip().splitlines() + assert len(lines) > 5 + # Basic CSV shape "t, value" or "t, clean, noisy" + parts = lines[0].split(",") + assert len(parts) in (2, 3) + + +def test_export_sensor_data_writes_json(flight_calisto, tmp_path, monkeypatch): + """Test that export_sensor_data writes JSON with sensor data. + + Validates that sensor data is exported as JSON with sensor names as keys + and measurement arrays as values. + + Parameters + ---------- + flight_calisto : rocketpy.Flight + Flight object to be tested. + tmp_path : pathlib.Path + Pytest fixture for temporary directories. + monkeypatch : pytest.MonkeyPatch + Pytest fixture for modifying attributes. + """ + + class DummySensor: + """Dummy sensor with name attribute for testing.""" + + def __init__(self, name): + self.name = name + + s1 = DummySensor("DummySensor") + monkeypatch.setattr( + flight_calisto, "sensor_data", {s1: [1.0, 2.0, 3.0]}, raising=False + ) + out = tmp_path / "sensors.json" + + FlightDataExporter(flight_calisto).export_sensor_data(str(out)) + + data = json.loads(out.read_text(encoding="utf-8")) + assert data["DummySensor"] == [1.0, 2.0, 3.0] + + +def test_export_data_default_variables(flight_calisto, tmp_path): + """Test export_data with default variables (full solution matrix). + + Validates that all state variables are exported correctly when no specific + variables are requested: position (x, y, z), velocity (vx, vy, vz), + quaternions (e0, e1, e2, e3), and angular velocities (w1, w2, w3). + + Parameters + ---------- + flight_calisto : rocketpy.Flight + Flight object to be tested. + tmp_path : pathlib.Path + Pytest fixture for temporary directories. + """ + file_name = tmp_path / "flight_data.csv" + FlightDataExporter(flight_calisto).export_data(str(file_name)) + + test_data = np.loadtxt(file_name, delimiter=",", skiprows=1) + + # Verify time column + assert np.allclose(flight_calisto.x[:, 0], test_data[:, 0], atol=1e-5) + + # Verify position + assert np.allclose(flight_calisto.x[:, 1], test_data[:, 1], atol=1e-5) + assert np.allclose(flight_calisto.y[:, 1], test_data[:, 2], atol=1e-5) + assert np.allclose(flight_calisto.z[:, 1], test_data[:, 3], atol=1e-5) + + # Verify velocity + assert np.allclose(flight_calisto.vx[:, 1], test_data[:, 4], atol=1e-5) + assert np.allclose(flight_calisto.vy[:, 1], test_data[:, 5], atol=1e-5) + assert np.allclose(flight_calisto.vz[:, 1], test_data[:, 6], atol=1e-5) + + # Verify quaternions + assert np.allclose(flight_calisto.e0[:, 1], test_data[:, 7], atol=1e-5) + assert np.allclose(flight_calisto.e1[:, 1], test_data[:, 8], atol=1e-5) + assert np.allclose(flight_calisto.e2[:, 1], test_data[:, 9], atol=1e-5) + assert np.allclose(flight_calisto.e3[:, 1], test_data[:, 10], atol=1e-5) + + # Verify angular velocities + assert np.allclose(flight_calisto.w1[:, 1], test_data[:, 11], atol=1e-5) + assert np.allclose(flight_calisto.w2[:, 1], test_data[:, 12], atol=1e-5) + assert np.allclose(flight_calisto.w3[:, 1], test_data[:, 13], atol=1e-5) + + +def test_export_data_custom_variables_and_time_step(flight_calisto, tmp_path): + """Test export_data with custom variables and time step. + + Validates that specific variables can be exported with custom time intervals, + including derived quantities like angle of attack. + + Parameters + ---------- + flight_calisto : rocketpy.Flight + Flight object to be tested. + tmp_path : pathlib.Path + Pytest fixture for temporary directories. + """ + file_name = tmp_path / "custom_flight_data.csv" + time_step = 0.1 + + FlightDataExporter(flight_calisto).export_data( + str(file_name), + "z", + "vz", + "e1", + "w3", + "angle_of_attack", + time_step=time_step, + ) + + test_data = np.loadtxt(file_name, delimiter=",", skiprows=1) + time_points = np.arange(flight_calisto.t_initial, flight_calisto.t_final, time_step) + + # Verify time column + assert np.allclose(time_points, test_data[:, 0], atol=1e-5) + + # Verify custom variables + assert np.allclose(flight_calisto.z(time_points), test_data[:, 1], atol=1e-5) + assert np.allclose(flight_calisto.vz(time_points), test_data[:, 2], atol=1e-5) + assert np.allclose(flight_calisto.e1(time_points), test_data[:, 3], atol=1e-5) + assert np.allclose(flight_calisto.w3(time_points), test_data[:, 4], atol=1e-5) + assert np.allclose( + flight_calisto.angle_of_attack(time_points), test_data[:, 5], atol=1e-5 + ) + + +def test_export_kml_trajectory(flight_calisto_robust, tmp_path): + """Test export_kml creates valid KML file with trajectory coordinates. + + Validates that the KML export contains correct latitude, longitude, and + altitude coordinates for the flight trajectory in absolute altitude mode. + + Parameters + ---------- + flight_calisto_robust : rocketpy.Flight + Flight object to be tested. + tmp_path : pathlib.Path + Pytest fixture for temporary directories. + """ + file_name = tmp_path / "trajectory.kml" + FlightDataExporter(flight_calisto_robust).export_kml( + str(file_name), time_step=None, extrude=True, altitude_mode="absolute" + ) + + # Parse KML coordinates + with open(file_name, "r") as kml_file: + for row in kml_file: + if row.strip().startswith(""): + coords_str = ( + row.strip() + .replace("", "") + .replace("", "") + ) + coords_list = coords_str.strip().split(" ") + + # Extract lon, lat, z from coordinates + parsed_coords = [c.split(",") for c in coords_list] + lon = [float(point[0]) for point in parsed_coords] + lat = [float(point[1]) for point in parsed_coords] + z = [float(point[2]) for point in parsed_coords] + + # Verify coordinates match flight data + assert np.allclose(flight_calisto_robust.latitude[:, 1], lat, atol=1e-3) + assert np.allclose(flight_calisto_robust.longitude[:, 1], lon, atol=1e-3) + assert np.allclose(flight_calisto_robust.z[:, 1], z, atol=1e-3) diff --git a/tests/unit/test_flight_export_deprecation.py b/tests/unit/test_flight_export_deprecation.py new file mode 100644 index 000000000..6fb6952b4 --- /dev/null +++ b/tests/unit/test_flight_export_deprecation.py @@ -0,0 +1,40 @@ +from unittest.mock import patch + +import pytest + +# TODO: these tests should be deleted after the deprecated methods are removed + + +def test_export_data_deprecated_emits_warning_and_delegates(flight_calisto, tmp_path): + """Expect: calling Flight.export_data emits DeprecationWarning and delegates to FlightDataExporter.export_data.""" + out = tmp_path / "out.csv" + with patch( + "rocketpy.simulation.flight_data_exporter.FlightDataExporter.export_data" + ) as spy: + with pytest.warns(DeprecationWarning): + flight_calisto.export_data(str(out)) + spy.assert_called_once() + + +def test_export_pressures_deprecated_emits_warning_and_delegates( + flight_calisto_robust, tmp_path +): + """Expect: calling Flight.export_pressures emits DeprecationWarning and delegates to FlightDataExporter.export_pressures.""" + out = tmp_path / "p.csv" + with patch( + "rocketpy.simulation.flight_data_exporter.FlightDataExporter.export_pressures" + ) as spy: + with pytest.warns(DeprecationWarning): + flight_calisto_robust.export_pressures(str(out), time_step=0.1) + spy.assert_called_once() + + +def test_export_kml_deprecated_emits_warning_and_delegates(flight_calisto, tmp_path): + """Expect: calling Flight.export_kml emits DeprecationWarning and delegates to FlightDataExporter.export_kml.""" + out = tmp_path / "traj.kml" + with patch( + "rocketpy.simulation.flight_data_exporter.FlightDataExporter.export_kml" + ) as spy: + with pytest.warns(DeprecationWarning): + flight_calisto.export_kml(str(out), time_step=0.5) + spy.assert_called_once() diff --git a/tests/unit/test_utilities.py b/tests/unit/test_utilities.py index 55d45de95..2bb8608b5 100644 --- a/tests/unit/test_utilities.py +++ b/tests/unit/test_utilities.py @@ -114,7 +114,7 @@ def test_fin_flutter_analysis(flight_calisto_custom_wind): assert np.isclose(flutter_mach(np.inf), 1.0048188594647927, atol=5e-3) assert np.isclose(safety_factor(0), 64.78797, atol=5e-3) assert np.isclose(safety_factor(10), 2.1948620401502072, atol=5e-3) - assert np.isclose(safety_factor(np.inf), 61.64222220697017, atol=5e-3) + assert np.isclose(safety_factor(np.inf), 61.669562809629035, atol=5e-3) def test_flutter_prints(flight_calisto_custom_wind):