diff --git a/.github/agents/rocketpy-reviewer.agent.md b/.github/agents/rocketpy-reviewer.agent.md new file mode 100644 index 000000000..be1b64b13 --- /dev/null +++ b/.github/agents/rocketpy-reviewer.agent.md @@ -0,0 +1,62 @@ +--- +description: "Physics-safe RocketPy code review agent. Use for pull request review, unit consistency checks, coordinate-frame validation, cached-property risk detection, and regression-focused test-gap analysis." +name: "RocketPy Reviewer" +tools: [read, search, execute] +argument-hint: "Review these changes for physics correctness and regression risk: " +user-invocable: true +--- +You are a RocketPy-focused reviewer for physics safety and regression risk. + +## Goals + +- Detect behavioral regressions and numerical/physics risks before merge. +- Validate unit consistency and coordinate/reference-frame correctness. +- Identify stale-cache risks when `@cached_property` interacts with mutable state. +- Check test coverage quality for changed behavior. +- Verify alignment with RocketPy workflow and contributor conventions. + +## Review Priorities + +1. Correctness and safety issues (highest severity). +2. Behavioral regressions and API compatibility. +3. Numerical stability and tolerance correctness. +4. Missing tests or weak assertions. +5. Documentation mismatches affecting users. +6. Workflow violations (test placement, branch/PR conventions, or missing validation evidence). + +## RocketPy-Specific Checks + +- SI units are explicit and consistent. +- Orientation conventions are unambiguous (`tail_to_nose`, `nozzle_to_combustion_chamber`, etc.). +- New/changed simulation logic does not silently invalidate cached values. +- Floating-point assertions use `pytest.approx` where needed. +- New fixtures are wired through `tests/conftest.py` when applicable. +- Test type is appropriate for scope (`unit`, `integration`, `acceptance`) and `all_info()`-style tests + are not misclassified. +- New behavior includes at least one regression-oriented test and relevant edge-case checks. +- For docs-affecting changes, references and paths remain valid and build warnings are addressed. +- Tooling recommendations match current repository setup (prefer Makefile plus `pyproject.toml` + settings when docs are outdated). + +## Validation Expectations + +- Prefer focused test runs first, then broader relevant suites. +- Recommend `make format` and `make lint` when style/lint risks are present. +- Recommend `make build-docs` when `.rst` files or API docs are changed. + +## Output Format + +Provide findings first, ordered by severity. +For each finding include: +- Severity: Critical, High, Medium, or Low +- Location: file path and line +- Why it matters: behavioral or physics risk +- Suggested fix: concrete, minimal change + +After findings, include: +- Open questions or assumptions +- Residual risks or testing gaps +- Brief change summary +- Suggested validation commands (only when useful) + +If no findings are identified, state that explicitly and still report residual risks/testing gaps. diff --git a/.github/copilot-instructions.md b/.github/copilot-instructions.md index f5366cb3b..382aa15e0 100644 --- a/.github/copilot-instructions.md +++ b/.github/copilot-instructions.md @@ -1,221 +1,80 @@ -# 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. +# RocketPy Workspace Instructions + +## Code Style +- Use snake_case for variables, functions, methods, and modules. Use descriptive names. +- Use PascalCase for classes and UPPER_SNAKE_CASE for constants. +- Keep lines at 88 characters and follow PEP 8 unless existing code in the target file differs. +- Run Ruff as the source of truth for formatting/import organization: + - `make format` + - `make lint` +- Use NumPy-style docstrings for public classes, methods, and functions, including units. +- In case of tooling drift between docs and config, prefer current repository tooling in `Makefile` + and `pyproject.toml`. + +## Architecture +- RocketPy is a modular Python library; keep feature logic in the correct package boundary: + - `rocketpy/simulation`: flight simulation and Monte Carlo orchestration. + - `rocketpy/rocket`, `rocketpy/motors`, `rocketpy/environment`: domain models. + - `rocketpy/mathutils`: numerical primitives and interpolation utilities. + - `rocketpy/plots`, `rocketpy/prints`: output and visualization layers. +- Prefer extending existing classes/patterns over introducing new top-level abstractions. +- Preserve public API stability in `rocketpy/__init__.py` exports. + +## Build and Test +- Use Makefile targets for OS-agnostic workflows: + - `make install` + - `make pytest` + - `make pytest-slow` + - `make coverage` + - `make coverage-report` + - `make build-docs` +- Before finishing code changes, run focused tests first, then broader relevant suites. +- When running Python directly in this workspace, prefer `.venv/Scripts/python.exe`. +- Slow tests are explicitly marked with `@pytest.mark.slow` and are run with `make pytest-slow`. +- For docs changes, check `make build-docs` output and resolve warnings/errors when practical. + +## Development Workflow +- Target pull requests to `develop` by default; `master` is the stable branch. +- Use branch names in `type/description` format, such as: + - `bug/` + - `doc/` + - `enh/` + - `mnt/` + - `tst/` +- Prefer rebasing feature branches on top of `develop` to keep history linear. +- Keep commit and PR titles explicit and prefixed with project acronyms when possible: + - `BUG`, `DOC`, `ENH`, `MNT`, `TST`, `BLD`, `REL`, `REV`, `STY`, `DEV`. + +## Conventions +- SI units are the default. Document units and coordinate-system references explicitly. +- Position/reference-frame arguments are critical in this codebase. Be explicit about orientation + (for example, `tail_to_nose`, `nozzle_to_combustion_chamber`). +- Include unit tests for new behavior. Follow AAA structure and clear test names. +- Use fixtures from `tests/fixtures`; if adding a new fixture module, update `tests/conftest.py`. +- Use `pytest.approx` for floating-point checks where appropriate. +- Use `@cached_property` for expensive computations when helpful, and be careful with stale-cache + behavior when underlying mutable state changes. +- Keep behavior backward compatible across the public API exported via `rocketpy/__init__.py`. +- Prefer extending existing module patterns over creating new top-level package structure. + +## Testing Taxonomy +- Unit tests are mandatory for new behavior. +- Unit tests in RocketPy can be sociable (real collaborators allowed) but should still be fast and + method-focused. +- Treat tests as integration tests when they are strongly I/O-oriented or broad across many methods, + including `all_info()` convention cases. +- Acceptance tests represent realistic user/flight scenarios and may compare simulation thresholds to + known flight data. + +## Documentation Links +- Contributor workflow and setup: `docs/development/setting_up.rst` +- Style and naming details: `docs/development/style_guide.rst` +- Testing philosophy and structure: `docs/development/testing.rst` +- API reference conventions: `docs/reference/index.rst` +- Domain/physics background: `docs/technical/index.rst` + +## Scoped Customizations +- Simulation-specific rules: `.github/instructions/simulation-safety.instructions.md` +- Test-authoring rules: `.github/instructions/tests-python.instructions.md` +- RST/Sphinx documentation rules: `.github/instructions/sphinx-docs.instructions.md` +- Specialized review persona: `.github/agents/rocketpy-reviewer.agent.md` diff --git a/.github/instructions/simulation-safety.instructions.md b/.github/instructions/simulation-safety.instructions.md new file mode 100644 index 000000000..cc2af5d27 --- /dev/null +++ b/.github/instructions/simulation-safety.instructions.md @@ -0,0 +1,41 @@ +--- +description: "Use when editing rocketpy/simulation code, including Flight state updates, Monte Carlo orchestration, post-processing, or cached computations. Covers simulation state safety, unit/reference-frame clarity, and regression checks." +name: "Simulation Safety" +applyTo: "rocketpy/simulation/**/*.py" +--- +# Simulation Safety Guidelines + +- Keep simulation logic inside `rocketpy/simulation` and avoid leaking domain behavior that belongs in + `rocketpy/rocket`, `rocketpy/motors`, or `rocketpy/environment`. +- Preserve public API behavior and exported names used by `rocketpy/__init__.py`. +- Prefer extending existing simulation components before creating new abstractions: + - `flight.py`: simulation state, integration flow, and post-processing. + - `monte_carlo.py`: orchestration and statistical execution workflows. + - `flight_data_exporter.py` and `flight_data_importer.py`: persistence and interchange. + - `flight_comparator.py`: comparative analysis outputs. +- Be explicit with physical units and reference frames in new parameters, attributes, and docstrings. +- For position/orientation-sensitive behavior, use explicit conventions (for example + `tail_to_nose`, `nozzle_to_combustion_chamber`) and avoid implicit assumptions. +- Treat state mutation carefully when cached values exist. +- If changes can invalidate `@cached_property` values, either avoid post-computation mutation or + explicitly invalidate affected caches in a controlled, documented way. +- Keep numerical behavior deterministic unless stochastic behavior is intentional and documented. +- For Monte Carlo and stochastic code paths, make randomness controllable and reproducible when tests + rely on it. +- Prefer vectorized NumPy operations for hot paths and avoid introducing Python loops in + performance-critical sections without justification. +- Guard against numerical edge cases (zero/near-zero denominators, interpolation limits, and boundary + conditions). +- Do not change default numerical tolerances or integration behavior without documenting motivation and + validating regression impact. +- Add focused regression tests for changed behavior, including edge cases and orientation-dependent + behavior. +- For floating-point expectations, use `pytest.approx` with meaningful tolerances. +- Run focused tests first, then broader relevant tests (`make pytest` and `make pytest-slow` when + applicable). + +See: +- `docs/development/testing.rst` +- `docs/development/style_guide.rst` +- `docs/development/setting_up.rst` +- `docs/technical/index.rst` diff --git a/.github/instructions/sphinx-docs.instructions.md b/.github/instructions/sphinx-docs.instructions.md new file mode 100644 index 000000000..8c24cac53 --- /dev/null +++ b/.github/instructions/sphinx-docs.instructions.md @@ -0,0 +1,32 @@ +--- +description: "Use when writing or editing docs/**/*.rst. Covers Sphinx/reStructuredText conventions, cross-references, toctree hygiene, and RocketPy unit/reference-frame documentation requirements." +name: "Sphinx RST Conventions" +applyTo: "docs/**/*.rst" +--- +# Sphinx and RST Guidelines + +- Follow existing heading hierarchy and style in the target document. +- Prefer linking to existing documentation pages instead of duplicating content. +- Use Sphinx cross-references where appropriate (`:class:`, `:func:`, `:mod:`, `:doc:`, `:ref:`). +- Keep API names and module paths consistent with current code exports. +- Document physical units and coordinate/reference-frame conventions explicitly. +- Include concise, practical examples when introducing new user-facing behavior. +- Keep prose clear and technical; avoid marketing language in development/reference docs. +- When adding a new page, update the relevant `toctree` so it appears in navigation. +- Use RocketPy docs build workflow: + - `make build-docs` from repository root for normal validation. + - If stale artifacts appear, clean docs build outputs via `cd docs && make clean`, then rebuild. +- Treat new Sphinx warnings/errors as issues to fix or explicitly call out in review notes. +- Keep `docs/index.rst` section structure coherent with user, development, reference, technical, and + examples navigation. +- Do not edit Sphinx-generated scaffolding files unless explicitly requested: + - `docs/Makefile` + - `docs/make.bat` +- For API docs, ensure references remain aligned with exported/public objects and current module paths. + +See: +- `docs/index.rst` +- `docs/development/build_docs.rst` +- `docs/development/style_guide.rst` +- `docs/reference/index.rst` +- `docs/technical/index.rst` diff --git a/.github/instructions/tests-python.instructions.md b/.github/instructions/tests-python.instructions.md new file mode 100644 index 000000000..1e9626142 --- /dev/null +++ b/.github/instructions/tests-python.instructions.md @@ -0,0 +1,36 @@ +--- +description: "Use when creating or editing pytest files in tests/. Enforces AAA structure, naming conventions, fixture usage, parameterization, slow-test marking, and numerical assertion practices for RocketPy." +name: "RocketPy Pytest Standards" +applyTo: "tests/**/*.py" +--- +# RocketPy Test Authoring Guidelines + +- Unit tests are mandatory for new behavior. +- Follow AAA structure in each test: Arrange, Act, Assert. +- Use descriptive test names matching project convention: + - `test_methodname` + - `test_methodname_stateundertest` + - `test_methodname_expectedbehaviour` +- Include docstrings that clearly state expected behavior and context. +- Prefer parameterization for scenario matrices instead of duplicated tests. +- Classify tests correctly: + - `tests/unit`: fast, method-focused tests (sociable unit tests are acceptable in RocketPy). + - `tests/integration`: broad multi-method/component interactions and strongly I/O-oriented cases. + - `tests/acceptance`: realistic end-user/flight scenarios with threshold-based expectations. +- By RocketPy convention, tests centered on `all_info()` behavior are integration tests. +- Reuse fixtures from `tests/fixtures` whenever possible. +- Keep fixture organization aligned with existing categories under `tests/fixtures` + (environment, flight, motor, rockets, surfaces, units, etc.). +- If you add a new fixture module, update `tests/conftest.py` so fixtures are discoverable. +- Keep tests deterministic: set seeds when randomness is involved and avoid unstable external + dependencies unless integration behavior explicitly requires them. +- Use `pytest.approx` for floating-point comparisons with realistic tolerances. +- Mark expensive tests with `@pytest.mark.slow` and ensure they can run under the project slow-test + workflow. +- Include at least one negative or edge-case assertion for new behaviors. +- When adding a bug fix, include a regression test that fails before the fix and passes after it. + +See: +- `docs/development/testing.rst` +- `docs/development/style_guide.rst` +- `docs/development/setting_up.rst` diff --git a/.github/workflows/changelog.yml b/.github/workflows/changelog.yml new file mode 100644 index 000000000..528ca1729 --- /dev/null +++ b/.github/workflows/changelog.yml @@ -0,0 +1,62 @@ +name: Populate Changelog +on: + pull_request: + types: [closed] + branches: + - develop + +permissions: + contents: write + +jobs: + Changelog: + if: github.event.pull_request.merged == true + runs-on: ubuntu-latest + + steps: + - name: Clone RocketPy + uses: actions/checkout@main + with: + repository: RocketPy-Team/RocketPy + ref: develop + token: ${{ secrets.RELEASE_TOKEN }} + + - name: Update Changelog + env: + PR_TITLE: ${{ github.event.pull_request.title }} + PR_NUMBER: ${{ github.event.pull_request.number }} + PR_LABELS: ${{ join(github.event.pull_request.labels.*.name, ',') }} + run: | + SECTION="### Added" + PREFIX="ENH" + + if [[ "$PR_LABELS" == *"Bug"* ]]; then + SECTION="### Fixed" + PREFIX="BUG" + elif [[ "$PR_LABELS" == *"Refactor"* ]]; then + SECTION="### Changed" + PREFIX="MNT" + elif [[ "$PR_LABELS" == *"Docs"* ]] && [[ "$PR_LABELS" == *"Git housekeeping"* ]]; then + SECTION="### Changed" + PREFIX="DOC" + elif [[ "$PR_LABELS" == *"Tests"* ]]; then + SECTION="### Changed" + PREFIX="TST" + elif [[ "$PR_LABELS" == *"Docs"* ]]; then + # Only documentation -> Added section + SECTION="### Added" + PREFIX="DOC" + fi + + ENTRY="- $PREFIX: $PR_TITLE [#$PR_NUMBER](https://github.com/RocketPy-Team/RocketPy/pull/$PR_NUMBER)" + SECTION_LINE=$(grep -n "^$SECTION$" CHANGELOG.md | head -1 | cut -d: -f1) + + sed -i "$((SECTION_LINE + 1))a\\$ENTRY" CHANGELOG.md + + - name: Push Changes + run: | + git config user.name "github-actions[bot]" + git config user.email "github-actions[bot]@users.noreply.github.com" + git add CHANGELOG.md + git commit -m "DOC: Update Changelog for PR #${{ github.event.pull_request.number }}" + git push diff --git a/.pylintrc b/.pylintrc index b0aaf655c..709250978 100644 --- a/.pylintrc +++ b/.pylintrc @@ -225,6 +225,11 @@ good-names=FlightPhases, center_of_mass_without_motor_to_CDM, motor_center_of_dry_mass_to_CDM, generic_motor_cesaroni_M1520, + R_phi, + R_delta, + R_pi, + R_uncanted, + R_body_to_fin, Re, # Reynolds number # Good variable names regexes, separated by a comma. If names match any regex, diff --git a/CHANGELOG.md b/CHANGELOG.md index 183479c4b..1cdd8479c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,7 +32,10 @@ Attention: The newest changes should be on top --> ### Added -- +- ENH: MNT: introduce pressure unit conversion when using forecast/reanalysis/ensemble data [#955](https://github.com/RocketPy-Team/RocketPy/pull/955) +- ENH: Auto Populate Changelog [#919](https://github.com/RocketPy-Team/RocketPy/pull/919) +- ENH: Adaptive Monte Carlo via Convergence Criteria [#922](https://github.com/RocketPy-Team/RocketPy/pull/922) +- TST: Add acceptance tests for 3DOF flight simulation based on Bella Lui rocket [#914](https://github.com/RocketPy-Team/RocketPy/pull/914) ### Changed @@ -40,7 +43,7 @@ Attention: The newest changes should be on top --> ### Fixed -- +- BUG: Add wraparound logic for wind direction in environment plots [#939](https://github.com/RocketPy-Team/RocketPy/pull/939) ## [v1.12.1] - 2026-04-03 diff --git a/docs/development/first_pr.rst b/docs/development/first_pr.rst index 5f3eb876c..c0b07e7f1 100644 --- a/docs/development/first_pr.rst +++ b/docs/development/first_pr.rst @@ -96,14 +96,13 @@ The CHANGELOG file ------------------ We keep track of the changes in the ``CHANGELOG.md`` file. -When you open a PR, you should add a new entry to the "Unreleased" section of the file. -This entry should simply be the title of your PR. +When you open a PR, you should see the "Unreleased" section of the file. +An entry will simply contain the title of your PR if merged. .. note:: - In the future we would like to automate the CHANGELOG update, but for now \ - it is a manual process, unfortunately. - + The CHANGELOG is auto-updated once a PR is merged based on the associated labels, \ + which are assigned by the maintainers. The review process ------------------ diff --git a/docs/development/style_guide.rst b/docs/development/style_guide.rst index 510a49560..15a80e5e4 100644 --- a/docs/development/style_guide.rst +++ b/docs/development/style_guide.rst @@ -162,15 +162,16 @@ Standard acronyms to start the commit message with are:: Pull Requests ^^^^^^^^^^^^^ -When opening a Pull Request, the name of the PR should be clear and concise. -Similarly to the commit messages, the PR name should start with an acronym indicating the type of PR -and then a brief description of the changes. +When opening a Pull Request, the title should be clear and concise. +It should contain only a brief desctiption of the changes without the acronym (e.g. ENH:, BUG:). +The maintainers will label your PR accordingly, which will add a prefix via a workflow to indicate +type of the PR in the CHANGELOG file. Here is an example of a good PR name: -- ``BUG: fix the Frequency Response plot of the Flight class`` +- ``fix the Frequency Response plot of the Flight class`` -The PR description explain the changes and motivation behind them. There is a template \ +The PR description explains the changes and motivation behind them. There is a template \ available when opening a PR that can be used to guide you through the process of both \ describing the changes and making sure all the necessary steps were taken. Of course, \ you can always modify the template or add more information if you think it is necessary. diff --git a/docs/notebooks/monte_carlo_analysis/monte_carlo_class_usage.ipynb b/docs/notebooks/monte_carlo_analysis/monte_carlo_class_usage.ipynb index 2fb46fa86..8181c03ba 100644 --- a/docs/notebooks/monte_carlo_analysis/monte_carlo_class_usage.ipynb +++ b/docs/notebooks/monte_carlo_analysis/monte_carlo_class_usage.ipynb @@ -800,6 +800,28 @@ ")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Alternatively, we can target an attribute using the method `MonteCarlo.simulate_convergence()` such that when the tolerance is met, the flight simulations would terminate early." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "test_dispersion.simulate_convergence(\n", + " target_attribute=\"apogee_time\",\n", + " target_confidence=0.95,\n", + " tolerance=0.5, # in seconds\n", + " max_simulations=1000,\n", + " batch_size=50,\n", + ")" + ] + }, { "attachments": {}, "cell_type": "markdown", diff --git a/docs/reference/classes/aero_surfaces/EllipticalFin.rst b/docs/reference/classes/aero_surfaces/EllipticalFin.rst new file mode 100644 index 000000000..3f1403e8e --- /dev/null +++ b/docs/reference/classes/aero_surfaces/EllipticalFin.rst @@ -0,0 +1,5 @@ +EllipticalFin Class +------------------- + +.. autoclass:: rocketpy.EllipticalFin + :members: \ No newline at end of file diff --git a/docs/reference/classes/aero_surfaces/Fin.rst b/docs/reference/classes/aero_surfaces/Fin.rst new file mode 100644 index 000000000..b9c70c25e --- /dev/null +++ b/docs/reference/classes/aero_surfaces/Fin.rst @@ -0,0 +1,5 @@ +Fin Class +--------- + +.. autoclass:: rocketpy.Fin + :members: \ No newline at end of file diff --git a/docs/reference/classes/aero_surfaces/FreeFormFin.rst b/docs/reference/classes/aero_surfaces/FreeFormFin.rst new file mode 100644 index 000000000..636be2c9e --- /dev/null +++ b/docs/reference/classes/aero_surfaces/FreeFormFin.rst @@ -0,0 +1,5 @@ +FreeFormFin Class +----------------- + +.. autoclass:: rocketpy.FreeFormFin + :members: \ No newline at end of file diff --git a/docs/reference/classes/aero_surfaces/FreeFormFins.rst b/docs/reference/classes/aero_surfaces/FreeFormFins.rst new file mode 100644 index 000000000..c641ce156 --- /dev/null +++ b/docs/reference/classes/aero_surfaces/FreeFormFins.rst @@ -0,0 +1,5 @@ +FreeFormFins Class +------------------ + +.. autoclass:: rocketpy.FreeFormFins + :members: \ No newline at end of file diff --git a/docs/reference/classes/aero_surfaces/TrapezoidalFin.rst b/docs/reference/classes/aero_surfaces/TrapezoidalFin.rst new file mode 100644 index 000000000..34fa631df --- /dev/null +++ b/docs/reference/classes/aero_surfaces/TrapezoidalFin.rst @@ -0,0 +1,5 @@ +TrapezoidalFin Class +-------------------- + +.. autoclass:: rocketpy.TrapezoidalFin + :members: \ No newline at end of file diff --git a/docs/reference/classes/aero_surfaces/_BaseFin.rst b/docs/reference/classes/aero_surfaces/_BaseFin.rst new file mode 100644 index 000000000..a04b02ee5 --- /dev/null +++ b/docs/reference/classes/aero_surfaces/_BaseFin.rst @@ -0,0 +1,5 @@ +_BaseFin Class +-------------- + +.. autoclass:: rocketpy.rocket.aero_surface.fins._base_fin._BaseFin + :members: \ No newline at end of file diff --git a/docs/reference/classes/aero_surfaces/index.rst b/docs/reference/classes/aero_surfaces/index.rst index c6e6a3efe..a3dad0417 100644 --- a/docs/reference/classes/aero_surfaces/index.rst +++ b/docs/reference/classes/aero_surfaces/index.rst @@ -11,7 +11,12 @@ AeroSurface Classes Fins TrapezoidalFins EllipticalFins + FreeFormFins + Fin + TrapezoidalFin + EllipticalFin + FreeFormFin RailButtons AirBrakes GenericSurface - LinearGenericSurface + LinearGenericSurface \ No newline at end of file diff --git a/docs/static/rocket/fin_forces.png b/docs/static/rocket/fin_forces.png new file mode 100644 index 000000000..41be9bfa2 Binary files /dev/null and b/docs/static/rocket/fin_forces.png differ diff --git a/docs/static/rocket/individual_fin_frame.png b/docs/static/rocket/individual_fin_frame.png new file mode 100644 index 000000000..67d417cdb Binary files /dev/null and b/docs/static/rocket/individual_fin_frame.png differ diff --git a/docs/technical/aerodynamics/individual_fins.rst b/docs/technical/aerodynamics/individual_fins.rst new file mode 100644 index 000000000..408352e16 --- /dev/null +++ b/docs/technical/aerodynamics/individual_fins.rst @@ -0,0 +1,498 @@ +.. _individual_fins: + +==================== +Individual Fin Model +==================== + +:Author: Mateus Stano Junqueira +:Date: March 2025 + +Introduction +============ + +There are currently three ways to model the aerodynamic effects of fins in RocketPy: + +1. By use of the :class:`rocketpy.GenericSurface` or :class:`rocketpy.LinearGenericSurface` classes, + which are defined by receiving the aerodynamic coefficients directly; +2. By use of the :class:`rocketpy.Fins` classes (:class:`rocketpy.TrapezoidalFins`, :class:`rocketpy.EllipticalFins`, :class:`rocketpy.FreeFormFins`), + which are defined by receiving the geometric parameters of the **fin set** and + calculating the aerodynamic coefficients internally with the Barrowman method; +3. By use of the :class:`rocketpy.Fin` classes (:class:`rocketpy.TrapezoidalFin`, :class:`rocketpy.EllipticalFin`, :class:`rocketpy.FreeFormFin`) + (which are the base classes for the fin set classes), which are defined by receiving the + geometric parameters of **one single fin** and calculating the aerodynamic + coefficients internally with the Barrowman method, with exception of the + moment coefficients, which have been derived in this document. + +This document will focus on all the mathematical derivations and implementations +of the :class:`rocketpy.Fin` classes. + +Fin Coordinate Frame +==================== + +The fin coordinate frame is defined as follows: + +- The origin is at the leading edge of the fin, at the intersection of the + fin's root chord and the fin's leading edge; +- The z-axis is aligned with the fin's root chord, pointing towards the fin's + trailing edge; +- The y-axis is aligned with the fin's span, pointing towards the fin's tip chord; +- The x-axis completes the right-handed coordinate system. + +The fin is postioned at the rocket's body, given a ``rocket_radius`` :math:`r`, +an ``angular_position`` :math:`\Gamma` and a position :math:`p` along the +rocket's body. The angular position is given according to +:ref:`Angular Position Inputs `. + + +The fin can also be rotated by a cant angle :math:`\delta`. + +The image below shows the fin coordinate frame: + +.. image:: ../../static/rocket/individual_fin_frame.png + :width: 800 + :align: center + +.. note:: + A positive cant angle :math:`\delta` produces a negative body axis rolling + moment at zero angle of attack. + +The rotation matrix from the fin coordinate frame to the rocket's body frame is +define by, first a rotation around the y-axis by 180 degrees: + +.. math:: + \mathbf{R}_{y(\pi)} = \begin{bmatrix} + -1 & 0 & 0 \\ + 0 & 1 & 0 \\ + 0 & 0 & -1 + \end{bmatrix} + +Then a rotation around the z-axis by the angle :math:`\Gamma`: + +.. math:: + \mathbf{R}_{z(\Gamma)} = \begin{bmatrix} + \cos(\Gamma) & -\sin(\Gamma) & 0 \\ + \sin(\Gamma) & \cos(\Gamma) & 0 \\ + 0 & 0 & 1 + \end{bmatrix} + +Then a rotation around the y-axis by the cant angle :math:`\delta`: + +.. math:: + \mathbf{R}_{y(\delta)} = \begin{bmatrix} + \cos(\delta) & 0 & \sin(\delta) \\ + 0 & 1 & 0 \\ + -\sin(\delta) & 0 & \cos(\delta) + \end{bmatrix} + +The final rotation matrix is given by: + +.. math:: + \mathbf{R} = \mathbf{R}_{y(\delta)} \cdot \mathbf{R}_{z(\Gamma)} \cdot \mathbf{R}_{y(\pi)} + + +The position of the fin's coordinate frame origin in the rocket's body frame +is calculated by first assuming no a fin frame with no cant angle, then +calculating the position of the fin's leading edge (with cant angle) in this +frame, and finally translating this position to the fin's position in the +rocket's body frame. The position of the fin's real leading edge in this no cant +angle fin frame is given by the point :math:`\mathbf{P}^{\delta}_{le_f}`: + +.. math:: + \mathbf{P}^{\delta}_{le_f} = \begin{bmatrix} + -\frac{Cr}{2} \sin(\delta) \\ + 0 \\ + \frac{Cr}{2} (1 - \cos(\delta)) + \end{bmatrix} + +Then, describing this point to the rocket's body frame orientation (no +translation): + +.. math:: + \mathbf{P}^{\delta}_{le_b} = (\mathbf{R}_{z(\Gamma)} \cdot \mathbf{R}_{y(\pi)}) \cdot \mathbf{P}^{\delta}_{le_f} + +The position of the fin's leading edge with no cant angle in the rocket's body +frame is given by: + +.. math:: + \mathbf{P}^{\overline{\delta}}_{le_b} = \begin{bmatrix} + -r \sin(\Gamma) \\ + r \cos(\Gamma) \\ + p + \end{bmatrix} + +Finally, we add the position of the fin's leading edge with no cant angle to the +position of the fin's leading edge with cant angle in the rocket's body frame: + +.. math:: + \mathbf{P}_{le_b} = \mathbf{P}^{\overline{\delta}}_{le_b} + \mathbf{P}^{\delta}_{le_b} + + +Center of Pressure Position +=========================== + +In the Fin Coordinate Frame, the center of pressure is given by the Barrowman +method, and will here only be defined symbolically: + +.. math:: + \mathbf{cp}_f = \begin{bmatrix} + cp_x \\ + cp_y \\ + cp_z + \end{bmatrix} + +The center of pressure position in the rocket's body frame is given by: + +.. math:: + \mathbf{cp}_{rocket} = \mathbf{R} \cdot \mathbf{cp}_f + \mathbf{P}_{le_b} + +Aerodynamic Forces +================== + +.. note:: + The aerodynamic coefficients are defined according the Barrowman method. + +Given a stream velocity in the fin frame :math:`\mathbf{v}_{0f} = [v_{0x}, v_{0y}, v_{0z}]^{T}`, +the effective angle of attack of the fin is given by: + +.. math:: + \alpha_f = \arctan\left(\frac{v_{0x}}{v_{0z}}\right) + +This can also be seen as the angle between the fin's root chord and the stream +velocity vector in the fin frame. + +The aerodynamic force in the x-direction of the fin is given by: + +.. math:: + F_{x} = \frac{1}{2} \cdot \rho \cdot \|\mathbf{v}_{0f}\|^2 \cdot A_{r} \cdot C_{N}(\alpha_f, Ma) + +Where :math:`A_{r}` is the reference area of the fin, and :math:`C_{N}` is the +normal force coefficient, which is a function of the angle of attack and the +Mach number :math:`Ma`. +This force is then transformed to the rocket's body frame by the rotation matrix: + +.. math:: + \begin{bmatrix} + F_{x} \\ + F_{y} \\ + F_{z} + \end{bmatrix}_{rocket} = \mathbf{R} \cdot \begin{bmatrix} + F_{x} \\ + 0 \\ + 0 + \end{bmatrix}_{fin} + +Then, the moments are calculated by the cross product of the center of pressure +and the aerodynamic force: + +.. math:: + \begin{bmatrix} + M_{x} \\ + M_{y} \\ + M_{z} + \end{bmatrix}_{rocket} = \mathbf{cp}_{rocket} \times \begin{bmatrix} + F_{x} \\ + F_{y} \\ + F_{z} + \end{bmatrix}_{rocket} + +From the Barrowman method, the moment along the center axis of the rocket +(:math:`M_{z}`) is still missing the damping term, which is given by: + +.. math:: + M_{damp} = \frac{1}{2} \cdot \rho \cdot \|v_{0}\| \cdot A_{r} \cdot L_{r}^2 \cdot C_{ld\omega}(Ma) \cdot \frac{1}{2} \cdot \omega_z + +.. math:: + M_{z \, \text{final}} = M_{z} + M_{damp} + +Where :math:`C_{ld}` is the roll moment damping coefficient, :math:`L_{r}` +is the reference length, which is equal to the rocket diameter, and +:math:`\omega_z` is the angular velocity of the rocket around the z-axis. + +Adding Individual Fins to the Rocket +==================================== +.. jupyter-execute:: + :hide-code: + :hide-output: + + from rocketpy import * + env = Environment(latitude=32.990254, longitude=-106.974998, elevation=1400) + Pro75M1670 = 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.397, + center_of_dry_mass_position=0.317, + nozzle_position=0, + burn_time=3.9, + throat_radius=11 / 1000, + coordinate_system_orientation="nozzle_to_combustion_chamber", + ) + # IMPORTANT: modify the file paths below to match your own system + + example_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", + ) + + rail_buttons = example_rocket.set_rail_buttons( + upper_button_position=0.0818, + lower_button_position=-0.618, + angular_position=45, + ) + example_rocket.add_motor(Pro75M1670, position=-1.255) + nose_cone = example_rocket.add_nose(length=0.55829, kind="vonKarman", position=1.278) + tail = example_rocket.add_tail( + top_radius=0.0635, bottom_radius=0.0435, length=0.060, position=-1.194656 + ) + example_rocket.add_trapezoidal_fins( + n=4, + root_chord=0.120, + tip_chord=0.060, + span=0.110, + cant_angle=0.0, + position=-1.04956, + airfoil=("../data/airfoils/NACA0012-radians.txt", "radians"), + ) + +Given a defined ``Rocket`` object, we can add individual fins to the rocket by +using the ``add_surfaces`` method. Here is an example of adding two canards +in the Calisto rocket from the :ref:`First Simulation ` example: + +.. jupyter-execute:: + + canard1 = TrapezoidalFin( + angular_position=0, + root_chord=0.060, + tip_chord=0.020, + span=0.03, + rocket_radius=example_rocket.radius, + cant_angle=0.5, + airfoil=("../data/airfoils/NACA0012-radians.txt", "radians"), + ) + canard2 = TrapezoidalFin( + angular_position=180, + root_chord=0.060, + tip_chord=0.020, + span=0.03, + rocket_radius=example_rocket.radius, + cant_angle=0.5, + airfoil=("../data/airfoils/NACA0012-radians.txt", "radians"), + ) + + # Position along the center axis of the rocket is specified here. + # If different positions are desired, the position can be specified as a list. + example_rocket.add_surfaces([canard1, canard2], positions = 0.35) + + example_rocket.draw(plane="yz") + +.. seealso:: + + There are three classes for defining fins in RocketPy given their geometry: + + - :class:`rocketpy.TrapezoidalFin` - For how to define a trapezoidal fin + - :class:`rocketpy.EllipticalFin` - For how to define an elliptical fin + - :class:`rocketpy.FreeFormFin` - For how to define a free form fin + + + +Fin Force Conventions +===================== + +.. - Explain positive cant angle resultant force +.. - if all fins have positive cant angle then they will all generate a force that will rotate the rocket in the same direction +.. which is negative roll +.. - show what a positive and negative cant angle on oposing fins looks like. (generate pitch moment -> pitch control) +.. - show what a positive and negative cant angle on oposing fins looks like. (generate yaw moment -> yaw control) +.. - example for 4 fins, show how to count the number of fins and how to access each of them + + + +Here we exemplify the fin force conventions relating the cant angle +(deflection angle) of the fins to the pitch, yaw and roll moments. We will +consider a rocket with four fins, to illustrate the concepts. The image below +show the sign convention for the forces acting on the fins, given positive cant +angles: + +.. image:: ../../static/rocket/fin_forces.png + :width: 800 + :align: center + + +Roll +^^^^ + +.. jupyter-execute:: + :hide-code: + :hide-output: + + example_rocket.aerodynamic_surfaces.pop() + example_rocket.aerodynamic_surfaces.pop() + + +A positive cant angle :math:`\delta` produces a negative roll moment at zero +angle of attack. Any fin with a positive cant angle will produce a negative roll +moment, and any fin with a negative cant angle will produce a positive roll +moment. + +Here is a flight of the calisto with canards defined with a positive cant angle: + +.. jupyter-execute:: + + + canard1 = TrapezoidalFin( + angular_position=0, + root_chord=0.060, + tip_chord=0.020, + span=0.03, + rocket_radius=example_rocket.radius, + cant_angle=0.5, + airfoil=("../data/airfoils/NACA0012-radians.txt", "radians"), + ) + canard2 = TrapezoidalFin( + angular_position=180, + root_chord=0.060, + tip_chord=0.020, + span=0.03, + rocket_radius=example_rocket.radius, + cant_angle=0.5, + airfoil=("../data/airfoils/NACA0012-radians.txt", "radians"), + ) + + example_rocket.add_surfaces([canard1, canard2], positions = 0.35) + + test_flight = Flight( + rocket=example_rocket, + environment=env, rail_length=5.2, inclination=85, heading=0, + terminate_on_apogee=True, + ) + + # Rolling Moment + test_flight.M3() + + # Rolling Speed + test_flight.w3() + + # Angle of attack + test_flight.partial_angle_of_attack.plot(test_flight.out_of_rail_time, 5) + + # Angle of sideslip + test_flight.angle_of_sideslip.plot(test_flight.out_of_rail_time, 5) + +Pitch +^^^^^ + +.. jupyter-execute:: + :hide-code: + :hide-output: + + example_rocket.aerodynamic_surfaces.pop() + example_rocket.aerodynamic_surfaces.pop() + +Given canards fins at 90 degrees and 270 degrees, having opposite cant angles, +a positive pitch moment will be generated. The following example shows the +effect of this configuration in the non-zero angle of attack flight of the +rocket: + +.. jupyter-execute:: + + + canard1 = TrapezoidalFin( + angular_position=90, + root_chord=0.060, + tip_chord=0.020, + span=0.03, + rocket_radius=example_rocket.radius, + cant_angle=0.5, + airfoil=("../data/airfoils/NACA0012-radians.txt", "radians"), + ) + canard2 = TrapezoidalFin( + angular_position=270, + root_chord=0.060, + tip_chord=0.020, + span=0.03, + rocket_radius=example_rocket.radius, + cant_angle=-0.5, + airfoil=("../data/airfoils/NACA0012-radians.txt", "radians"), + ) + + example_rocket.add_surfaces([canard1, canard2], positions = 0.35) + + test_flight = Flight( + rocket=example_rocket, + environment=env, rail_length=5.2, inclination=85, heading=0, + terminate_on_apogee=True, + ) + + # Angle of attack + test_flight.partial_angle_of_attack.plot(test_flight.out_of_rail_time, 5) + + # Angle of sideslip + test_flight.angle_of_sideslip.plot(test_flight.out_of_rail_time, 5) + +Yaw +^^^ + +.. jupyter-execute:: + :hide-code: + :hide-output: + + example_rocket.aerodynamic_surfaces.pop() + example_rocket.aerodynamic_surfaces.pop() + +Given opposing canards at 0 degrees and 180 degrees, having opposite cant angles, +a positive yaw moment will be generated. The following example shows the +effect of this configuration in the non-zero angle of attack flight of the +rocket: + +.. jupyter-execute:: + + + canard1 = TrapezoidalFin( + angular_position=0, + root_chord=0.060, + tip_chord=0.020, + span=0.03, + rocket_radius=example_rocket.radius, + cant_angle=0.5, + airfoil=("../data/airfoils/NACA0012-radians.txt", "radians"), + ) + canard2 = TrapezoidalFin( + angular_position=180, + root_chord=0.060, + tip_chord=0.020, + span=0.03, + rocket_radius=example_rocket.radius, + cant_angle=-0.5, + airfoil=("../data/airfoils/NACA0012-radians.txt", "radians"), + ) + + example_rocket.add_surfaces([canard1, canard2], positions = 0.35) + + test_flight = Flight( + rocket=example_rocket, + environment=env, rail_length=5.2, inclination=85, heading=0, + terminate_on_apogee=True, + ) + + # Angle of attack + test_flight.partial_angle_of_attack.plot(test_flight.out_of_rail_time, 5) + + # Angle of sideslip + test_flight.angle_of_sideslip.plot(test_flight.out_of_rail_time, 5) + + + + + diff --git a/docs/technical/index.rst b/docs/technical/index.rst index b96e611ed..73583eba9 100644 --- a/docs/technical/index.rst +++ b/docs/technical/index.rst @@ -13,6 +13,7 @@ in their code. Equations of Motion v0 Equations of Motion v1 Elliptical Fins + Individual Fin Roll Moment Sensitivity Analysis References diff --git a/docs/user/environment/1-atm-models/ensemble.rst b/docs/user/environment/1-atm-models/ensemble.rst index 504cbfe60..8dffaac00 100644 --- a/docs/user/environment/1-atm-models/ensemble.rst +++ b/docs/user/environment/1-atm-models/ensemble.rst @@ -25,8 +25,8 @@ Global Ensemble Forecast System (GEFS) .. danger:: - **GEFS shortcut unavailable**: ``file="GEFS"`` is currently disabled in - RocketPy because NOMADS OPeNDAP is deactivated for this endpoint. + **GEFS unavailable**: ``file="GEFS"`` is currently disabled in + RocketPy because NOMADS OPeNDAP has been deactivated. .. note:: diff --git a/docs/user/environment/1-atm-models/forecast.rst b/docs/user/environment/1-atm-models/forecast.rst index ac91504e0..92f5b6b9e 100644 --- a/docs/user/environment/1-atm-models/forecast.rst +++ b/docs/user/environment/1-atm-models/forecast.rst @@ -6,8 +6,8 @@ Forecasts Weather forecasts can be used to set the atmospheric model in RocketPy. Here, we will showcase how to import global forecasts such as GFS, as well as -local forecasts like NAM and RAP for North America, all available through -OPeNDAP on the `NOAA's NCEP NOMADS `_ website. +local forecasts like NAM, RAP and HRRR for North America, all available through +OPeNDAP on the `UCAR THREDDS `_ server. Other generic forecasts can also be imported. .. important:: @@ -22,6 +22,10 @@ Other generic forecasts can also be imported. Global Forecast System (GFS) ---------------------------- +GFS is NOAA's global numerical weather prediction model. It provides worldwide +atmospheric forecasts and is usually a good default choice when you need broad +coverage, consistent availability, and launch planning several days ahead. + Using the latest forecast from GFS is simple. Set the atmospheric model to ``forecast`` and specify that GFS is the file you want. Note that since data is downloaded from a remote OPeNDAP server, this line of code can @@ -48,9 +52,34 @@ take longer than usual. `GFS overview page `_. +Artificial Intelligence Global Forecast System (AIGFS) +------------------------------------------------------ + +AIGFS is a global AI-based forecast product distributed through the same THREDDS +ecosystem used by other RocketPy forecast inputs. It is useful when you want a +global forecast alternative to traditional physics-only models. + +RocketPy supports the latest AIGFS global forecast through THREDDS. + +.. jupyter-execute:: + + env_aigfs = Environment(date=tomorrow) + env_aigfs.set_atmospheric_model(type="forecast", file="AIGFS") + env_aigfs.plots.atmospheric_model() + +.. note:: + + AIGFS is currently available as a global 0.25 degree forecast product on + UCAR THREDDS. + + North American Mesoscale Forecast System (NAM) ---------------------------------------------- +NAM is a regional forecast model focused on North America. It is best suited +for launches inside its coverage area when you want finer regional detail than +global models typically provide. + You can also request the latest forecasts from NAM. Since this is a regional model for North America, you need to specify latitude and longitude points within North America. @@ -78,6 +107,10 @@ We will use **SpacePort America** for this, represented by coordinates Rapid Refresh (RAP) ------------------- +RAP is a short-range, high-frequency regional model for North America. It is +especially useful for near-term operations, where fast update cycles are more +important than long forecast horizon. + The Rapid Refresh (RAP) model is another regional model for North America. It is similar to NAM, but with a higher resolution and a shorter forecast range. The same coordinates for SpacePort America will be used. @@ -111,6 +144,17 @@ The same coordinates for SpacePort America will be used. High Resolution Window (HIRESW) ------------------------------- +HIRESW is a convection-allowing, high-resolution regional system designed to +resolve local weather structure better than coarser grids. It is most useful +for short-range, local analysis where small-scale wind and weather features +matter. + +The High Resolution Window (HIRESW) model is a sophisticated weather forecasting +system that operates at a high spatial resolution of approximately 3 km. +It utilizes two main dynamical cores: the Advanced Research WRF (WRF-ARW) and +the Finite Volume Cubed Sphere (FV3), each designed to enhance the accuracy of +weather predictions. + .. danger:: **HIRESW shortcut unavailable**: ``file="HIRESW"`` is currently disabled in @@ -121,6 +165,33 @@ you can still load it explicitly by passing the path/URL in ``file`` and an appropriate mapping in ``dictionary``. +High-Resolution Rapid Refresh (HRRR) +------------------------------------ + +HRRR is a high-resolution, short-range forecast model for North America with +hourly updates. It is generally best for day-of-launch weather assessment and +rapidly changing local conditions. + +RocketPy supports HRRR through a dedicated THREDDS shortcut. +Like NAM and RAP, HRRR is a regional model over North America. + +If you have a HIRESW-compatible dataset from another provider (or a local copy), +you can still load it explicitly by passing the path/URL in ``file`` and an +appropriate mapping in ``dictionary``. + + env_hrrr = Environment( + date=now_plus_twelve, + latitude=32.988528, + longitude=-106.975056, + ) + env_hrrr.set_atmospheric_model(type="forecast", file="HRRR") + env_hrrr.plots.atmospheric_model() + +.. note:: + + HRRR is a high-resolution regional model with approximately 2.5 km grid + spacing over CONUS. Availability depends on upstream THREDDS data services. + Using Windy Atmosphere ---------------------- @@ -154,6 +225,10 @@ to EuRoC's launch area in Portugal. ECMWF ^^^^^ +ECMWF (HRES) is a global, high-skill forecast model known for strong +medium-range performance. It is often a good choice for mission planning when +you need reliable synoptic-scale forecasts several days ahead. + We can use the ``ECMWF`` model from Windy.com. .. jupyter-execute:: @@ -173,6 +248,10 @@ We can use the ``ECMWF`` model from Windy.com. GFS ^^^ +Windy's GFS option provides NOAA's global model through Windy's interface. It +is a practical baseline for global coverage and for comparing against other +models when assessing forecast uncertainty. + The ``GFS`` model is also available on Windy.com. This is the same model as described in the :ref:`global-forecast-system` section. @@ -186,6 +265,10 @@ described in the :ref:`global-forecast-system` section. ICON ^^^^ +ICON is DWD's global weather model, available in Windy for broad-scale +forecasting. It is useful as an independent global model source to cross-check +wind and temperature trends against GFS or ECMWF. + The ICON model is a global weather forecasting model already available on Windy.com. .. jupyter-execute:: @@ -203,6 +286,10 @@ The ICON model is a global weather forecasting model already available on Windy. ICON-EU ^^^^^^^ +ICON-EU is the regional European configuration of ICON, with higher spatial +detail over Europe than ICON-Global. It is best for European launch sites when +regional structure is important. + The ICON-EU model is a regional weather forecasting model available on Windy.com. .. code-block:: python @@ -228,4 +315,4 @@ Also, the servers may be down or may face high traffic. .. seealso:: To browse available NCEP model collections on UCAR THREDDS, visit - `THREDDS NCEP Catalog `_. \ No newline at end of file + `THREDDS NCEP Catalog `_. diff --git a/docs/user/environment/1-atm-models/soundings.rst b/docs/user/environment/1-atm-models/soundings.rst index 279750df5..4cf82543f 100644 --- a/docs/user/environment/1-atm-models/soundings.rst +++ b/docs/user/environment/1-atm-models/soundings.rst @@ -59,7 +59,7 @@ Integrated Global Radiosonde Archive (IGRA). These options can be retrieved as a text file in GSD format. However, RocketPy no longer provides a dedicated ``set_atmospheric_model`` type for -NOAA RUC Soundings. +NOAA RUC Soundings, since NOAA has discontinued the OPENDAP service. .. note:: diff --git a/docs/user/environment/3-further/other_apis.rst b/docs/user/environment/3-further/other_apis.rst index 01d4b9a30..37a9a0949 100644 --- a/docs/user/environment/3-further/other_apis.rst +++ b/docs/user/environment/3-further/other_apis.rst @@ -89,8 +89,10 @@ Instead of a custom dictionary, you can pass a built-in mapping name in the - ``"ECMWF_v0"`` - ``"NOAA"`` - ``"GFS"`` +- ``"AIGFS"`` - ``"NAM"`` - ``"RAP"`` +- ``"HRRR"`` - ``"HIRESW"`` (mapping available; latest-model shortcut currently disabled) - ``"GEFS"`` (mapping available; latest-model shortcut currently disabled) - ``"MERRA2"`` @@ -116,10 +118,7 @@ legacy aliases: - ``"NAM_LEGACY"`` - ``"NOAA_LEGACY"`` - ``"RAP_LEGACY"`` -- ``"CMC_LEGACY"`` - ``"GEFS_LEGACY"`` -- ``"HIRESW_LEGACY"`` -- ``"MERRA2_LEGACY"`` Legacy aliases primarily cover older variable naming patterns such as ``lev``, ``tmpprs``, ``hgtprs``, ``ugrdprs`` and ``vgrdprs``. diff --git a/docs/user/rocket/rocket_axes.rst b/docs/user/rocket/rocket_axes.rst index 56fe17cad..aa99cb231 100644 --- a/docs/user/rocket/rocket_axes.rst +++ b/docs/user/rocket/rocket_axes.rst @@ -61,6 +61,8 @@ The following figure shows the two possibilities for the user input coordinate s the same as the **Body Axes Coordinate System**. The origin of the coordinate \ system may still be different. +.. _angular_position: + Angular Position Inputs ~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/rocketpy/__init__.py b/rocketpy/__init__.py index 852a16aef..129317bea 100644 --- a/rocketpy/__init__.py +++ b/rocketpy/__init__.py @@ -29,8 +29,11 @@ AeroSurface, AirBrakes, Components, + EllipticalFin, EllipticalFins, + Fin, Fins, + FreeFormFin, FreeFormFins, GenericSurface, LinearGenericSurface, @@ -40,6 +43,7 @@ RailButtons, Rocket, Tail, + TrapezoidalFin, TrapezoidalFins, ) from .sensitivity import SensitivityModel diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py index 8e379800c..3b585fd44 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -11,10 +11,12 @@ import pytz from rocketpy.environment.fetchers import ( + fetch_aigfs_file_return_dataset, fetch_atmospheric_data_from_windy, fetch_gefs_ensemble, fetch_gfs_file_return_dataset, fetch_hiresw_file_return_dataset, + fetch_hrrr_file_return_dataset, fetch_nam_file_return_dataset, fetch_open_elevation, fetch_rap_file_return_dataset, @@ -369,9 +371,11 @@ def __initialize_constants(self): self.__weather_model_map = WeatherModelMapping() self.__atm_type_file_to_function_map = { "forecast": { + "AIGFS": fetch_aigfs_file_return_dataset, "GFS": fetch_gfs_file_return_dataset, "NAM": fetch_nam_file_return_dataset, "RAP": fetch_rap_file_return_dataset, + "HRRR": fetch_hrrr_file_return_dataset, "HIRESW": fetch_hiresw_file_return_dataset, }, "ensemble": { @@ -665,9 +669,11 @@ def __resolve_dictionary_for_dataset(self, dictionary, dataset): def __validate_dictionary(self, file, dictionary): # removed CMC until it is fixed. available_models = [ + "AIGFS", "GFS", "NAM", "RAP", + "HRRR", "HIRESW", "GEFS", "ERA5", @@ -1114,6 +1120,7 @@ def set_atmospheric_model( # pylint: disable=too-many-statements temperature=None, wind_u=0, wind_v=0, + pressure_conversion_factor="Pa", ): """Define the atmospheric model for this Environment. @@ -1132,8 +1139,9 @@ def set_atmospheric_model( # pylint: disable=too-many-statements - ``"windy"``: one of ``"ECMWF"``, ``"GFS"``, ``"ICON"`` or ``"ICONEU"``. - ``"forecast"``: local path, OPeNDAP URL, open - ``netCDF4.Dataset``, or one of ``"GFS"``, ``"NAM"`` or ``"RAP"`` - for the latest available forecast. + ``netCDF4.Dataset``, or one of ``"AIGFS"``, ``"GFS"``, + ``"NAM"``, ``"RAP"``, ``"HRRR"`` or ``"HIRESW"`` for the + latest available forecast. - ``"reanalysis"``: local path, OPeNDAP URL, or open ``netCDF4.Dataset``. - ``"ensemble"``: local path, OPeNDAP URL, open @@ -1143,8 +1151,9 @@ def set_atmospheric_model( # pylint: disable=too-many-statements Variable-name mapping for ``"forecast"``, ``"reanalysis"`` and ``"ensemble"``. It may be a custom dictionary or a built-in mapping name (for example: ``"ECMWF"``, ``"ECMWF_v0"``, - ``"NOAA"``, ``"GFS"``, ``"NAM"``, ``"RAP"``, ``"HIRESW"``, - ``"GEFS"``, ``"MERRA2"`` or ``"CMC"``). + ``"NOAA"``, ``"AIGFS"``, ``"GFS"``, ``"NAM"``, ``"RAP"``, + ``"HRRR"``, ``"HIRESW"``, ``"GEFS"``, ``"MERRA2"`` or + ``"CMC"``). If ``dictionary`` is omitted and ``file`` is one of RocketPy's latest-model shortcuts, the matching built-in mapping is selected @@ -1208,6 +1217,12 @@ def set_atmospheric_model( # pylint: disable=too-many-statements m/s). Finally, a callable or function is also accepted. The function should take one argument, the height above sea level in meters and return a corresponding wind-v in m/s. + pressure_conversion_factor : string, int, float + This defines the pressure conversion factor to Pa when type is + ``forecast`` or ``reanalysis``. The pressure unit from the data may + not be in Pascal, so the correction is necessary. Valid strings are + ("mbar", "hPa", "Pa"), or a strictly positive number if using a + custom pressure unit. Returns ------- @@ -1257,6 +1272,28 @@ def set_atmospheric_model( # pylint: disable=too-many-statements case "windy": self.process_windy_atmosphere(file) case "forecast" | "reanalysis" | "ensemble": + conversion_factor = 1 + if not isinstance(pressure_conversion_factor, (float, int, str)): + raise ValueError( + "Argument 'pressure_conversion_factor' must be numeric or a standard pressure unit ('mbar', 'hPa', 'Pa')!" + ) + if isinstance(pressure_conversion_factor, (float, int)): + if pressure_conversion_factor <= 0: + raise ValueError( + "Argument 'pressure_conversion_factor' must be strictly positive!" + ) + else: + conversion_factor = pressure_conversion_factor + if isinstance(pressure_conversion_factor, str): + if pressure_conversion_factor.lower() in ("mbar", "hpa"): + conversion_factor = 100 + elif pressure_conversion_factor.lower() == "pa": + conversion_factor = 1 + else: + raise ValueError( + "Argument 'pressure_conversion_factor' unit must be a standard pressure unit ('mbar', 'hPa', 'Pa')!" + ) + if isinstance(file, str): shortcut_map = self.__atm_type_file_to_function_map.get(type, {}) matching_shortcut = next( @@ -1297,9 +1334,11 @@ def set_atmospheric_model( # pylint: disable=too-many-statements dataset = fetch_function() if fetch_function is not None else file if type in ["forecast", "reanalysis"]: - self.process_forecast_reanalysis(dataset, dictionary) + self.process_forecast_reanalysis( + dataset, dictionary, conversion_factor=conversion_factor + ) else: - self.process_ensemble(dataset, dictionary) + self.process_ensemble(dataset, dictionary, conversion_factor) case _: # pragma: no cover raise ValueError(f"Unknown model type '{type}'.") @@ -1678,7 +1717,7 @@ def process_wyoming_sounding(self, file): # pylint: disable=too-many-statements # Save maximum expected height self._max_expected_height = data_array[-1, 1] - def process_forecast_reanalysis(self, file, dictionary): # pylint: disable=too-many-locals,too-many-statements + def process_forecast_reanalysis(self, file, dictionary, conversion_factor): # pylint: disable=too-many-locals,too-many-statements """Import and process atmospheric data from weather forecasts and reanalysis given as ``netCDF`` or ``OPeNDAP`` files. Sets pressure, temperature, wind-u and wind-v @@ -1730,6 +1769,9 @@ def process_forecast_reanalysis(self, file, dictionary): # pylint: disable=too- "u_wind": "ugrdprs", "v_wind": "vgrdprs", } + conversion_factor : float, int + Specifies the factor by which the pressure will be multiplied + in order to transform it to Pascal. Returns ------- @@ -1761,13 +1803,17 @@ def process_forecast_reanalysis(self, file, dictionary): # pylint: disable=too- # Some THREDDS datasets use projected x/y coordinates. if dictionary.get("projection") is not None: projection_variable = data.variables[dictionary["projection"]] - x_units = getattr(lon_array, "units", "m") - target_lon, target_lat = geodesic_to_lambert_conformal( - self.latitude, - self.longitude, - projection_variable, - x_units=x_units, - ) + if dictionary.get("projection") == "LambertConformal_Projection": + x_units = getattr(lon_array, "units", "m") + target_lon, target_lat = geodesic_to_lambert_conformal( + self.latitude, + self.longitude, + projection_variable, + x_units=x_units, + ) + else: + target_lon = self.longitude + target_lat = self.latitude else: target_lon = self.longitude target_lat = self.latitude @@ -1778,7 +1824,7 @@ def process_forecast_reanalysis(self, file, dictionary): # pylint: disable=too- _, lat_index = find_latitude_index(target_lat, lat_array) # Get pressure level data from file - levels = get_pressure_levels_from_file(data, dictionary) + levels = get_pressure_levels_from_file(data, dictionary, conversion_factor) # Get geopotential data from file try: @@ -1979,7 +2025,7 @@ def process_forecast_reanalysis(self, file, dictionary): # pylint: disable=too- # Close weather data data.close() - def process_ensemble(self, file, dictionary): # pylint: disable=too-many-locals,too-many-statements + def process_ensemble(self, file, dictionary, conversion_factor): # pylint: disable=too-many-locals,too-many-statements """Import and process atmospheric data from weather ensembles given as ``netCDF`` or ``OPeNDAP`` files. Sets pressure, temperature, wind-u and wind-v profiles and surface elevation obtained from a weather @@ -2030,6 +2076,9 @@ def process_ensemble(self, file, dictionary): # pylint: disable=too-many-locals "u_wind": "ugrdprs", "v_wind": "vgrdprs", } + conversion_factor : float, int + Specifies the factor by which the pressure will be multiplied + in order to transform it to Pascal. See also -------- @@ -2065,13 +2114,17 @@ class for some dictionary examples. # coordinate system before locating the nearest grid cell. if dictionary.get("projection") is not None: projection_variable = data.variables[dictionary["projection"]] - x_units = getattr(lon_array, "units", "m") - target_lon, target_lat = geodesic_to_lambert_conformal( - self.latitude, - self.longitude, - projection_variable, - x_units=x_units, - ) + if dictionary.get("projection") == "LambertConformal_Projection": + x_units = getattr(lon_array, "units", "m") + target_lon, target_lat = geodesic_to_lambert_conformal( + self.latitude, + self.longitude, + projection_variable, + x_units=x_units, + ) + else: + target_lon = self.longitude + target_lat = self.latitude else: target_lon = self.longitude target_lat = self.latitude @@ -2090,7 +2143,7 @@ class for some dictionary examples. num_members = 1 # Get pressure level data from file - levels = get_pressure_levels_from_file(data, dictionary) + levels = get_pressure_levels_from_file(data, dictionary, conversion_factor) inverse_dictionary = {v: k for k, v in dictionary.items()} param_dictionary = { diff --git a/rocketpy/environment/fetchers.py b/rocketpy/environment/fetchers.py index de63d53ad..64d5dd479 100644 --- a/rocketpy/environment/fetchers.py +++ b/rocketpy/environment/fetchers.py @@ -195,6 +195,78 @@ def fetch_rap_file_return_dataset(max_attempts=10, base_delay=2): raise RuntimeError("Unable to load latest weather data for RAP through " + file_url) +def fetch_hrrr_file_return_dataset(max_attempts=10, base_delay=2): + """Fetches the latest HRRR (High-Resolution Rapid Refresh) dataset from + the NOAA's GrADS data server using the OpenDAP protocol. + + Parameters + ---------- + max_attempts : int, optional + The maximum number of attempts to fetch the dataset. Default is 10. + base_delay : int, optional + The base delay in seconds between attempts. Default is 2. + + Returns + ------- + netCDF4.Dataset + The HRRR dataset. + + Raises + ------ + RuntimeError + If unable to load the latest weather data for HRRR. + """ + file_url = "https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/HRRR/CONUS_2p5km/Best" + attempt_count = 0 + while attempt_count < max_attempts: + try: + return netCDF4.Dataset(file_url) + except OSError: + attempt_count += 1 + time.sleep(base_delay**attempt_count) + + raise RuntimeError( + "Unable to load latest weather data for HRRR through " + file_url + ) + + +def fetch_aigfs_file_return_dataset(max_attempts=10, base_delay=2): + """Fetches the latest AIGFS (Artificial Intelligence GFS) dataset from + the NOAA's GrADS data server using the OpenDAP protocol. + + Parameters + ---------- + max_attempts : int, optional + The maximum number of attempts to fetch the dataset. Default is 10. + base_delay : int, optional + The base delay in seconds between attempts. Default is 2. + + Returns + ------- + netCDF4.Dataset + The AIGFS dataset. + + Raises + ------ + RuntimeError + If unable to load the latest weather data for AIGFS. + """ + file_url = ( + "https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/AIGFS/Global_0p25deg/Best" + ) + attempt_count = 0 + while attempt_count < max_attempts: + try: + return netCDF4.Dataset(file_url) + except OSError: + attempt_count += 1 + time.sleep(base_delay**attempt_count) + + raise RuntimeError( + "Unable to load latest weather data for AIGFS through " + file_url + ) + + def fetch_hiresw_file_return_dataset(max_attempts=10, base_delay=2): """Fetches the latest HiResW (High-Resolution Window) dataset from the NOAA's GrADS data server using the OpenDAP protocol. diff --git a/rocketpy/environment/tools.py b/rocketpy/environment/tools.py index 4a06ef2c4..816d9dd12 100644 --- a/rocketpy/environment/tools.py +++ b/rocketpy/environment/tools.py @@ -169,7 +169,7 @@ def geodesic_to_lambert_conformal(lat, lon, projection_variable, x_units="m"): ## These functions are meant to be used with netcdf4 datasets -def get_pressure_levels_from_file(data, dictionary): +def get_pressure_levels_from_file(data, dictionary, conversion_factor): """Extracts pressure levels from a netCDF4 dataset and converts them to Pa. Parameters @@ -178,6 +178,9 @@ def get_pressure_levels_from_file(data, dictionary): The netCDF4 dataset containing the pressure level data. dictionary : dict A dictionary mapping variable names to dataset keys. + conversion_factor : float, int + Specifies the factor by which the pressure will be multiplied + in order to transform it to Pascal. Returns ------- @@ -190,8 +193,7 @@ def get_pressure_levels_from_file(data, dictionary): If the pressure levels cannot be read from the file. """ try: - # Convert mbar to Pa - levels = 100 * data.variables[dictionary["level"]][:] + levels = conversion_factor * data.variables[dictionary["level"]][:] except KeyError as e: raise ValueError( "Unable to read pressure levels from file. Check file and dictionary." diff --git a/rocketpy/environment/weather_model_mapping.py b/rocketpy/environment/weather_model_mapping.py index b054a35c4..c8617a523 100644 --- a/rocketpy/environment/weather_model_mapping.py +++ b/rocketpy/environment/weather_model_mapping.py @@ -48,11 +48,11 @@ class WeatherModelMapping: "u_wind": "ugrdprs", "v_wind": "vgrdprs", } - NAM = { + AIGFS = { "time": "time", - "latitude": "y", - "longitude": "x", - "projection": "LambertConformal_Projection", + "latitude": "lat", + "longitude": "lon", + "projection": "LatLon_Projection", "level": "isobaric", "temperature": "Temperature_isobaric", "surface_geopotential_height": None, @@ -73,6 +73,19 @@ class WeatherModelMapping: "u_wind": "ugrdprs", "v_wind": "vgrdprs", } + NAM = { + "time": "time", + "latitude": "y", + "longitude": "x", + "projection": "LambertConformal_Projection", + "level": "isobaric", + "temperature": "Temperature_isobaric", + "surface_geopotential_height": None, + "geopotential_height": "Geopotential_height_isobaric", + "geopotential": None, + "u_wind": "u-component_of_wind_isobaric", + "v_wind": "v-component_of_wind_isobaric", + } ECMWF_v0 = { "time": "time", "latitude": "latitude", @@ -148,20 +161,20 @@ class WeatherModelMapping: "u_wind": "ugrdprs", "v_wind": "vgrdprs", } - CMC = { + HRRR = { "time": "time", - "latitude": "lat", - "longitude": "lon", - "level": "lev", - "ensemble": "ens", - "temperature": "tmpprs", + "latitude": "y", + "longitude": "x", + "projection": "LambertConformal_Projection", + "level": "isobaric", + "temperature": "Temperature_isobaric", "surface_geopotential_height": None, - "geopotential_height": "hgtprs", + "geopotential_height": "Geopotential_height_isobaric", "geopotential": None, - "u_wind": "ugrdprs", - "v_wind": "vgrdprs", + "u_wind": "u-component_of_wind_isobaric", + "v_wind": "v-component_of_wind_isobaric", } - CMC_LEGACY = { + CMC = { "time": "time", "latitude": "lat", "longitude": "lon", @@ -211,17 +224,6 @@ class WeatherModelMapping: "u_wind": "ugrdprs", "v_wind": "vgrdprs", } - HIRESW_LEGACY = { - "time": "time", - "latitude": "lat", - "longitude": "lon", - "level": "lev", - "temperature": "tmpprs", - "surface_geopotential_height": "hgtsfc", - "geopotential_height": "hgtprs", - "u_wind": "ugrdprs", - "v_wind": "vgrdprs", - } MERRA2 = { "time": "time", "latitude": "lat", @@ -235,19 +237,6 @@ class WeatherModelMapping: "u_wind": "U", "v_wind": "V", } - MERRA2_LEGACY = { - "time": "time", - "latitude": "lat", - "longitude": "lon", - "level": "lev", - "temperature": "T", - "surface_geopotential_height": None, - "surface_geopotential": "PHIS", - "geopotential_height": "H", - "geopotential": None, - "u_wind": "U", - "v_wind": "V", - } def __init__(self): """Build the lookup table with default and legacy mapping aliases.""" @@ -255,6 +244,7 @@ def __init__(self): self.all_dictionaries = { "GFS": self.GFS, "GFS_LEGACY": self.GFS_LEGACY, + "AIGFS": self.AIGFS, "NAM": self.NAM, "NAM_LEGACY": self.NAM_LEGACY, "ECMWF_v0": self.ECMWF_v0, @@ -263,14 +253,12 @@ def __init__(self): "NOAA_LEGACY": self.NOAA_LEGACY, "RAP": self.RAP, "RAP_LEGACY": self.RAP_LEGACY, + "HRRR": self.HRRR, "CMC": self.CMC, - "CMC_LEGACY": self.CMC_LEGACY, "GEFS": self.GEFS, "GEFS_LEGACY": self.GEFS_LEGACY, "HIRESW": self.HIRESW, - "HIRESW_LEGACY": self.HIRESW_LEGACY, "MERRA2": self.MERRA2, - "MERRA2_LEGACY": self.MERRA2_LEGACY, } def get(self, model): diff --git a/rocketpy/plots/aero_surface_plots.py b/rocketpy/plots/aero_surface_plots.py index 397ad8f55..eb97ce19b 100644 --- a/rocketpy/plots/aero_surface_plots.py +++ b/rocketpy/plots/aero_surface_plots.py @@ -1,3 +1,5 @@ +# pylint: disable=too-many-statements + from abc import ABC, abstractmethod import matplotlib.pyplot as plt @@ -139,14 +141,97 @@ class _FinsPlots(_AeroSurfacePlots): """Abstract class that contains all fin plots. This class inherits from the _AeroSurfacePlots class.""" - @abstractmethod - def draw(self, *, filename=None): - pass + def airfoil(self, *, filename=None): + """Plots the airfoil information when the fin has an airfoil shape. If + the fin does not have an airfoil shape, this method does nothing. + + Parameters + ---------- + filename : str | None, optional + The path the plot should be saved to. By default None, in which case + the plot will be shown instead of saved. + + Returns + ------- + None + """ + + if self.aero_surface.airfoil: + print("Airfoil lift curve:") + self.aero_surface.airfoil_cl.plot_1d(force_data=True, filename=filename) + + def roll(self, *, filename=None): + """Plots the roll parameters of the fin set. + + Parameters + ---------- + filename : str | None, optional + The path the plot should be saved to. By default None, in which case + the plot will be shown instead of saved. + + Returns + ------- + None + """ + print("Roll parameters:") + self.aero_surface.roll_parameters[0](filename=filename) + self.aero_surface.roll_parameters[1](filename=filename) + + def lift(self, *, filename=None): + """Plots the lift coefficient of the aero surface as a function of Mach + and the angle of attack. A 3D plot is expected. See the rocketpy.Function + class for more information on how this plot is made. Also, this method + plots the lift coefficient considering a single fin and the lift + coefficient considering all fins. + + Parameters + ---------- + filename : str | None, optional + The path the plot should be saved to. By default None, in which case + the plot will be shown instead of saved. + + Returns + ------- + None + """ + print("Lift coefficient:") + self.aero_surface.cl(filename=filename) + self.aero_surface.clalpha_single_fin(filename=filename) + self.aero_surface.clalpha_multiple_fins(filename=filename) + + def all(self, *, filename=None): + """Plots all available fin plots. + + Parameters + ---------- + filename : str | None, optional + The path the plot should be saved to. By default None, in which case + the plot will be shown instead of saved. + + Returns + ------- + None + """ + self.draw(filename=filename) + self.airfoil(filename=filename) + self.roll(filename=filename) + self.lift(filename=filename) - def airfoil(self): + +class _FinPlots(_AeroSurfacePlots): + """Abstract class that contains all fin plots. This class inherits from the + _AeroSurfacePlots class.""" + + def airfoil(self, *, filename=None): """Plots the airfoil information when the fin has an airfoil shape. If the fin does not have an airfoil shape, this method does nothing. + Parameters + ---------- + filename : str | None, optional + The path the plot should be saved to. By default None, in which case + the plot will be shown instead of saved. + Returns ------- None @@ -154,52 +239,68 @@ def airfoil(self): if self.aero_surface.airfoil: print("Airfoil lift curve:") - self.aero_surface.airfoil_cl.plot_1d(force_data=True) + self.aero_surface.airfoil_cl.plot_1d(force_data=True, filename=filename) - def roll(self): + def roll(self, *, filename=None): """Plots the roll parameters of the fin set. + Parameters + ---------- + filename : str | None, optional + The path the plot should be saved to. By default None, in which case + the plot will be shown instead of saved. + Returns ------- None """ print("Roll parameters:") - self.aero_surface.roll_parameters[0]() - self.aero_surface.roll_parameters[1]() + self.aero_surface.roll_parameters[0](filename=filename) + self.aero_surface.roll_parameters[1](filename=filename) - def lift(self): + def lift(self, *, filename=None): """Plots the lift coefficient of the aero surface as a function of Mach and the angle of attack. A 3D plot is expected. See the rocketpy.Function class for more information on how this plot is made. Also, this method plots the lift coefficient considering a single fin and the lift coefficient considering all fins. + Parameters + ---------- + filename : str | None, optional + The path the plot should be saved to. By default None, in which case + the plot will be shown instead of saved. + Returns ------- None """ print("Lift coefficient:") - self.aero_surface.cl() - self.aero_surface.clalpha_single_fin() - self.aero_surface.clalpha_multiple_fins() + self.aero_surface.cl(filename=filename) + self.aero_surface.clalpha_single_fin(filename=filename) - def all(self): + def all(self, *, filename=None): """Plots all available fin plots. + Parameters + ---------- + filename : str | None, optional + The path the plot should be saved to. By default None, in which case + the plot will be shown instead of saved. + Returns ------- None """ - self.draw() - self.airfoil() - self.roll() - self.lift() + self.draw(filename=filename) + self.airfoil(filename=filename) + self.roll(filename=filename) + self.lift(filename=filename) class _TrapezoidalFinsPlots(_FinsPlots): """Class that contains all trapezoidal fin plots.""" - # pylint: disable=too-many-statements def draw(self, *, filename=None): """Draw the fin shape along with some important information, including the center line, the quarter line and the center of pressure position. @@ -325,10 +426,129 @@ def draw(self, *, filename=None): show_or_save_plot(filename) +class _TrapezoidalFinPlots(_FinPlots): + """Class that contains all trapezoidal fin plots.""" + + def draw(self, *, filename=None): + """Draw the fin shape along with some important information, including + the center line, the quarter line and the center of pressure position. + + Returns + ------- + None + """ + # Color cycle [#348ABD, #A60628, #7A68A6, #467821, #D55E00, #CC79A7, + # #56B4E9, #009E73, #F0E442, #0072B2] + # Fin + leading_edge = plt.Line2D( + (0, self.aero_surface.sweep_length), + (0, self.aero_surface.span), + color="#A60628", + ) + tip = plt.Line2D( + ( + self.aero_surface.sweep_length, + self.aero_surface.sweep_length + self.aero_surface.tip_chord, + ), + (self.aero_surface.span, self.aero_surface.span), + color="#A60628", + ) + back_edge = plt.Line2D( + ( + self.aero_surface.sweep_length + self.aero_surface.tip_chord, + self.aero_surface.root_chord, + ), + (self.aero_surface.span, 0), + color="#A60628", + ) + root = plt.Line2D((self.aero_surface.root_chord, 0), (0, 0), color="#A60628") + + # Center and Quarter line + center_line = plt.Line2D( + ( + self.aero_surface.root_chord / 2, + self.aero_surface.sweep_length + self.aero_surface.tip_chord / 2, + ), + (0, self.aero_surface.span), + color="#7A68A6", + alpha=0.35, + linestyle="--", + label="Center Line", + ) + quarter_line = plt.Line2D( + ( + self.aero_surface.root_chord / 4, + self.aero_surface.sweep_length + self.aero_surface.tip_chord / 4, + ), + (0, self.aero_surface.span), + color="#7A68A6", + alpha=1, + linestyle="--", + label="Quarter Line", + ) + + # Center of pressure + cp_point = [self.aero_surface.cpz, self.aero_surface.Yma] + + # Mean Aerodynamic Chord + yma_start = ( + self.aero_surface.sweep_length + * (self.aero_surface.root_chord + 2 * self.aero_surface.tip_chord) + / (3 * (self.aero_surface.root_chord + self.aero_surface.tip_chord)) + ) + yma_end = ( + 2 * self.aero_surface.root_chord**2 + + self.aero_surface.root_chord * self.aero_surface.sweep_length + + 2 * self.aero_surface.root_chord * self.aero_surface.tip_chord + + 2 * self.aero_surface.sweep_length * self.aero_surface.tip_chord + + 2 * self.aero_surface.tip_chord**2 + ) / (3 * (self.aero_surface.root_chord + self.aero_surface.tip_chord)) + yma_line = plt.Line2D( + (yma_start, yma_end), + (self.aero_surface.Yma, self.aero_surface.Yma), + color="#467821", + linestyle="--", + label="Mean Aerodynamic Chord", + ) + + # Plotting + fig = plt.figure(figsize=(7, 4)) + with plt.style.context("bmh"): + ax = fig.add_subplot(111) + + # Fin + ax.add_line(leading_edge) + ax.add_line(tip) + ax.add_line(back_edge) + ax.add_line(root) + + ax.add_line(center_line) + ax.add_line(quarter_line) + ax.add_line(yma_line) + ax.scatter(*cp_point, label="Center of Pressure", color="red", s=100, zorder=10) + ax.scatter(*cp_point, facecolors="none", edgecolors="red", s=500, zorder=10) + + # Plot settings + xlim = ( + self.aero_surface.root_chord + if self.aero_surface.sweep_length + self.aero_surface.tip_chord + < self.aero_surface.root_chord + else self.aero_surface.sweep_length + self.aero_surface.tip_chord + ) + ax.set_xlim(0, xlim * 1.1) + ax.set_ylim(0, self.aero_surface.span * 1.1) + ax.set_xlabel("Root chord (m)") + ax.set_ylabel("Span (m)") + ax.set_title("Trapezoidal Fin Cross Section") + ax.legend(bbox_to_anchor=(1.05, 1.0), loc="upper left") + + plt.tight_layout() + show_or_save_plot(filename) + + class _EllipticalFinsPlots(_FinsPlots): """Class that contains all elliptical fin plots.""" - # pylint: disable=too-many-statements def draw(self, *, filename=None): """Draw the fin shape along with some important information. These being: the center line and the center of pressure position. @@ -404,10 +624,153 @@ def draw(self, *, filename=None): show_or_save_plot(filename) +class _EllipticalFinPlots(_FinPlots): + """Class that contains all elliptical fin plots.""" + + def draw(self, *, filename=None): + """Draw the fin shape along with some important information. + These being: the center line and the center of pressure position. + + Returns + ------- + None + """ + # Ellipse + ellipse = Ellipse( + (self.aero_surface.root_chord / 2, 0), + self.aero_surface.root_chord, + self.aero_surface.span * 2, + fill=False, + edgecolor="#A60628", + linewidth=2, + ) + + # Mean Aerodynamic Chord # From Barrowman's theory + yma_length = 8 * self.aero_surface.root_chord / (3 * np.pi) + yma_start = (self.aero_surface.root_chord - yma_length) / 2 + yma_end = ( + self.aero_surface.root_chord + - (self.aero_surface.root_chord - yma_length) / 2 + ) + yma_line = plt.Line2D( + (yma_start, yma_end), + (self.aero_surface.Yma, self.aero_surface.Yma), + label="Mean Aerodynamic Chord", + color="#467821", + ) + + # Center Line + center_line = plt.Line2D( + (self.aero_surface.root_chord / 2, self.aero_surface.root_chord / 2), + (0, self.aero_surface.span), + color="#7A68A6", + alpha=0.35, + linestyle="--", + label="Center Line", + ) + + # Center of pressure + cp_point = [self.aero_surface.cpz, self.aero_surface.Yma] + + # Plotting + fig = plt.figure(figsize=(7, 4)) + with plt.style.context("bmh"): + ax = fig.add_subplot(111) + ax.add_patch(ellipse) + ax.add_line(yma_line) + ax.add_line(center_line) + ax.scatter(*cp_point, label="Center of Pressure", color="red", s=100, zorder=10) + ax.scatter(*cp_point, facecolors="none", edgecolors="red", s=500, zorder=10) + + # Plot settings + ax.set_xlim(0, self.aero_surface.root_chord) + ax.set_ylim(0, self.aero_surface.span * 1.1) + ax.set_xlabel("Root chord (m)") + ax.set_ylabel("Span (m)") + ax.set_title("Elliptical Fin Cross Section") + ax.legend(bbox_to_anchor=(1.05, 1.0), loc="upper left") + + plt.tight_layout() + show_or_save_plot(filename) + + class _FreeFormFinsPlots(_FinsPlots): """Class that contains all free form fin plots.""" - # pylint: disable=too-many-statements + def draw(self, *, filename=None): + """Draw the fin shape along with some important information. + These being: the center line and the center of pressure position. + + Parameters + ---------- + filename : str | None, optional + The path the plot should be saved to. By default None, in which case + the plot will be shown instead of saved. Supported file endings are: + eps, jpg, jpeg, pdf, pgf, png, ps, raw, rgba, svg, svgz, tif, tiff + and webp (these are the formats supported by matplotlib). + + Returns + ------- + None + """ + # Color cycle [#348ABD, #A60628, #7A68A6, #467821, #D55E00, #CC79A7, + # #56B4E9, #009E73, #F0E442, #0072B2] + + # Center of pressure + cp_point = [self.aero_surface.cpz, self.aero_surface.Yma] + + # Mean Aerodynamic Chord + yma_line = plt.Line2D( + ( + self.aero_surface.mac_lead, + self.aero_surface.mac_lead + self.aero_surface.mac_length, + ), + (self.aero_surface.Yma, self.aero_surface.Yma), + color="#467821", + linestyle="--", + label="Mean Aerodynamic Chord", + ) + + # Plotting + fig = plt.figure(figsize=(7, 4)) + with plt.style.context("bmh"): + ax = fig.add_subplot(111) + + # Fin + ax.scatter( + self.aero_surface.shape_vec[0], + self.aero_surface.shape_vec[1], + color="#A60628", + ) + ax.plot( + self.aero_surface.shape_vec[0], + self.aero_surface.shape_vec[1], + color="#A60628", + ) + # line from the last point to the first point + ax.plot( + [self.aero_surface.shape_vec[0][-1], self.aero_surface.shape_vec[0][0]], + [self.aero_surface.shape_vec[1][-1], self.aero_surface.shape_vec[1][0]], + color="#A60628", + ) + + ax.add_line(yma_line) + ax.scatter(*cp_point, label="Center of Pressure", color="red", s=100, zorder=10) + ax.scatter(*cp_point, facecolors="none", edgecolors="red", s=500, zorder=10) + + # Plot settings + ax.set_xlabel("Root chord (m)") + ax.set_ylabel("Span (m)") + ax.set_title("Free Form Fin Cross Section") + ax.legend(bbox_to_anchor=(1.05, 1.0), loc="upper left") + + plt.tight_layout() + show_or_save_plot(filename) + + +class _FreeFormFinPlots(_FinPlots): + """Class that contains all free form fin plots.""" + def draw(self, *, filename=None): """Draw the fin shape along with some important information. These being: the center line and the center of pressure position. diff --git a/rocketpy/plots/rocket_plots.py b/rocketpy/plots/rocket_plots.py index e208c775f..47da8a78b 100644 --- a/rocketpy/plots/rocket_plots.py +++ b/rocketpy/plots/rocket_plots.py @@ -1,8 +1,9 @@ import matplotlib.pyplot as plt import numpy as np +from rocketpy.mathutils.vector_matrix import Vector from rocketpy.motors import EmptyMotor, HybridMotor, LiquidMotor, SolidMotor -from rocketpy.rocket.aero_surface import Fins, NoseCone, Tail +from rocketpy.rocket.aero_surface import Fin, Fins, NoseCone, Tail from rocketpy.rocket.aero_surface.generic_surface import GenericSurface from .plot_helpers import show_or_save_plot @@ -183,7 +184,7 @@ def draw(self, vis_args=None, plane="xz", *, filename=None): and webp (these are the formats supported by matplotlib). """ - self.__validate_aerodynamic_surfaces() + self.__validate_aerodynamic_surfaces(plane) if vis_args is None: vis_args = { @@ -203,9 +204,9 @@ def draw(self, vis_args=None, plane="xz", *, filename=None): csys = self.rocket._csys reverse = csys == 1 - self.rocket.aerodynamic_surfaces.sort_by_position(reverse=reverse) + surfaces = self.rocket.aerodynamic_surfaces.sort_by_position(reverse=reverse) - drawn_surfaces = self._draw_aerodynamic_surfaces(ax, vis_args, plane) + drawn_surfaces = self._draw_aerodynamic_surfaces(ax, vis_args, plane, surfaces) last_radius, last_x = self._draw_tubes(ax, drawn_surfaces, vis_args) self._draw_motor(last_radius, last_x, ax, vis_args) self._draw_rail_buttons(ax, vis_args) @@ -221,13 +222,15 @@ def draw(self, vis_args=None, plane="xz", *, filename=None): plt.tight_layout() show_or_save_plot(filename) - def __validate_aerodynamic_surfaces(self): + def __validate_aerodynamic_surfaces(self, plane): if not self.rocket.aerodynamic_surfaces: raise ValueError( "The rocket must have at least one aerodynamic surface to be drawn." ) + if plane not in ("xz", "yz"): + raise ValueError("The plane must be 'xz' or 'yz'. The default is 'xz'.") - def _draw_aerodynamic_surfaces(self, ax, vis_args, plane): + def _draw_aerodynamic_surfaces(self, ax, vis_args, plane, surfaces): """Draws the aerodynamic surfaces and saves the position of the points of interest for the tubes.""" # List of drawn surfaces with the position of points of interest @@ -240,13 +243,17 @@ def _draw_aerodynamic_surfaces(self, ax, vis_args, plane): # diameter changes. The final point of the last surface is the final # point of the last tube - for surface, position in self.rocket.aerodynamic_surfaces: + for surface, position in surfaces: if isinstance(surface, NoseCone): self._draw_nose_cone(ax, surface, position.z, drawn_surfaces, vis_args) elif isinstance(surface, Tail): self._draw_tail(ax, surface, position.z, drawn_surfaces, vis_args) elif isinstance(surface, Fins): - self._draw_fins(ax, surface, position.z, drawn_surfaces, vis_args) + self._draw_fins( + ax, surface, position.z, drawn_surfaces, vis_args, plane + ) + elif isinstance(surface, Fin): + self._draw_fin(ax, surface, position, drawn_surfaces, vis_args, plane) elif isinstance(surface, GenericSurface): self._draw_generic_surface( ax, surface, position, drawn_surfaces, vis_args, plane @@ -309,13 +316,15 @@ def _draw_tail(self, ax, surface, position, drawn_surfaces, vis_args): # Add the tail to the list of drawn surfaces drawn_surfaces.append((surface, position, surface.bottom_radius, x_tail[-1])) - def _draw_fins(self, ax, surface, position, drawn_surfaces, vis_args): + def _draw_fins(self, ax, surface, position, drawn_surfaces, vis_args, plane): """Draws the fins and saves the position of the points of interest for the tubes.""" num_fins = surface.n x_fin = -self.rocket._csys * surface.shape_vec[0] + position y_fin = surface.shape_vec[1] + surface.rocket_radius - rotation_angles = [2 * np.pi * i / num_fins for i in range(num_fins)] + rotation_angles = np.array([2 * np.pi * i / num_fins for i in range(num_fins)]) + if plane == "xz": + rotation_angles -= np.pi / 2 for angle in rotation_angles: # Create a rotation matrix for the current angle around the x-axis @@ -327,13 +336,6 @@ def _draw_fins(self, ax, surface, position, drawn_surfaces, vis_args): # Extract x and y coordinates of the rotated points x_rotated, y_rotated = rotated_points_2d - # Project points above the XY plane back into the XY plane (set z-coordinate to 0) - x_rotated = np.where( - rotated_points_2d[1] > 0, rotated_points_2d[0], x_rotated - ) - y_rotated = np.where( - rotated_points_2d[1] > 0, rotated_points_2d[1], y_rotated - ) ax.plot( x_rotated, y_rotated, @@ -343,6 +345,56 @@ def _draw_fins(self, ax, surface, position, drawn_surfaces, vis_args): drawn_surfaces.append((surface, position, surface.rocket_radius, x_rotated[-1])) + def _draw_fin(self, ax, surface, position, drawn_surfaces, vis_args, plane): + """Draws individual fins.""" + + # Get shape vec + xs = surface.shape_vec[0] + ys = surface.shape_vec[1] + zs = np.zeros_like(xs) + + # Define shape in fin coordinate system + x_fin = -zs + y_fin = ys + z_fin = xs + points = np.column_stack((x_fin, y_fin, z_fin)) + + # Move drawing coordinates to center of fin for cant angle rotation + xd = np.array([0, 0, max(xs) / 2]) + points -= xd + + # Rotate to body coordinate system + for i, p in enumerate(points): + points[i] = surface._rotation_fin_to_body @ Vector(p) + + rotated_xd = surface._rotation_fin_to_body @ Vector(xd) + points += np.array(rotated_xd) + + # Back to the drawing system + x_fin_rotated = points[:, 0] + y_fin_rotated = points[:, 1] + z_fin_rotated = points[:, 2] + + if plane == "xz": + x_rotated = self.rocket._csys * z_fin_rotated + position.z + y_rotated = x_fin_rotated + position.x + elif plane == "yz": + x_rotated = self.rocket._csys * z_fin_rotated + position.z + y_rotated = y_fin_rotated + position.y + else: # pragma: no cover + raise ValueError("Plane must be 'xz' or 'yz'.") + + ax.plot( + x_rotated, + y_rotated, + color=vis_args["fins"], + linewidth=vis_args["line_width"], + ) + + drawn_surfaces.append( + (surface, position.z, surface.rocket_radius, x_rotated[-1]) + ) + def _draw_generic_surface( self, ax, @@ -445,7 +497,7 @@ def _draw_motor(self, last_radius, last_x, ax, vis_args): self._draw_nozzle_tube(last_radius, last_x, nozzle_position, ax, vis_args) - def _generate_motor_patches(self, total_csys, ax): # pylint: disable=unused-argument + def _generate_motor_patches(self, total_csys, ax): """Generates motor patches for drawing""" motor_patches = [] @@ -654,7 +706,7 @@ def all(self): # Rocket draw if len(self.rocket.aerodynamic_surfaces) > 0: - print("\nRocket Draw") + print("\nRocket Drawing") print("-" * 40) self.draw() diff --git a/rocketpy/prints/aero_surface_prints.py b/rocketpy/prints/aero_surface_prints.py index 4eb42b08d..cc36f1b01 100644 --- a/rocketpy/prints/aero_surface_prints.py +++ b/rocketpy/prints/aero_surface_prints.py @@ -1,6 +1,7 @@ from abc import ABC, abstractmethod +# TODO: the rocketpy/prints/aero_surface_prints.py file could be separated into different, smaller files. class _AeroSurfacePrints(ABC): def __init__(self, aero_surface): self.aero_surface = aero_surface @@ -77,10 +78,8 @@ def geometry(self): print("-------------------------------------") print(f"Number of fins: {self.aero_surface.n}") print(f"Reference rocket radius: {self.aero_surface.rocket_radius:.3f} m") - try: + if hasattr(self.aero_surface, "tip_chord"): print(f"Tip chord: {self.aero_surface.tip_chord:.3f} m") - except AttributeError: - pass # it isn't a trapezoidal fin, just don't worry about tip chord print(f"Root chord: {self.aero_surface.root_chord:.3f} m") print(f"Span: {self.aero_surface.span:.3f} m") print( @@ -156,9 +155,102 @@ def lift(self): "Lift Coefficient derivative (single fin) at Mach 0 and AoA 0: " f"{self.aero_surface.clalpha_single_fin(0):.3f}" ) + + def all(self): + """Prints all information of the fin set. + + Returns + ------- + None + """ + self.identity() + self.geometry() + self.airfoil() + self.roll() + self.lift() + + +class _FinPrints(_AeroSurfacePrints): + def geometry(self): + print("Geometric information of the fin set:") + print("-------------------------------------") + print(f"Reference rocket radius: {self.aero_surface.rocket_radius:.3f} m") + if hasattr(self.aero_surface, "tip_chord"): + print(f"Tip chord: {self.aero_surface.tip_chord:.3f} m") + print(f"Root chord: {self.aero_surface.root_chord:.3f} m") + print(f"Span: {self.aero_surface.span:.3f} m") + print( + f"Cant angle: {self.aero_surface.cant_angle:.3f} ° or " + f"{self.aero_surface.cant_angle_rad:.3f} rad" + ) + print(f"Longitudinal section area: {self.aero_surface.Af:.3f} m²") + print(f"Aspect ratio: {self.aero_surface.AR:.3f} ") + print(f"Gamma_c: {self.aero_surface.gamma_c:.3f} m") + print(f"Mean aerodynamic chord: {self.aero_surface.Yma:.3f} m\n") + + def airfoil(self): + """Prints out airfoil related information of the fin set. + + Returns + ------- + None + """ + if self.aero_surface.airfoil: + print("Airfoil information:") + print("--------------------") + print( + "Number of points defining the lift curve: " + f"{len(self.aero_surface.airfoil_cl.x_array)}" + ) + print( + "Lift coefficient derivative at Mach 0 and AoA 0: " + f"{self.aero_surface.clalpha(0):.5f} 1/rad\n" + ) + + def roll(self): + """Prints out information about roll parameters + of the fin set. + + Returns + ------- + None + """ + print("Roll information of the fin set:") + print("--------------------------------") + print( + f"Geometric constant: {self.aero_surface.roll_geometrical_constant:.3f} m" + ) + print( + "Damping interference factor: " + f"{self.aero_surface.roll_damping_interference_factor:.3f} rad" + ) + print( + "Forcing interference factor: " + f"{self.aero_surface.roll_forcing_interference_factor:.3f} rad\n" + ) + + def lift(self): + """Prints out information about lift parameters + of the fin set. + + Returns + ------- + None + """ + print("Lift information of the fin set:") + print("--------------------------------") print( - "Lift Coefficient derivative (fin set) at Mach 0 and AoA 0: " - f"{self.aero_surface.clalpha_multiple_fins(0):.3f}" + "Lift interference factor: " + f"{self.aero_surface.lift_interference_factor:.3f} m" + ) + print( + "Center of Pressure position in local coordinates: " + f"({self.aero_surface.cpx:.3f}, {self.aero_surface.cpy:.3f}, " + f"{self.aero_surface.cpz:.3f})" + ) + print( + "Lift Coefficient derivative (single fin) at Mach 0 and AoA 0: " + f"{self.aero_surface.clalpha_single_fin(0):.3f}" ) def all(self): @@ -179,14 +271,26 @@ class _TrapezoidalFinsPrints(_FinsPrints): """Class that contains all trapezoidal fins prints.""" +class _TrapezoidalFinPrints(_FinPrints): + """Class that contains all trapezoidal fin prints.""" + + class _EllipticalFinsPrints(_FinsPrints): """Class that contains all elliptical fins prints.""" +class _EllipticalFinPrints(_FinPrints): + """Class that contains all elliptical fin prints.""" + + class _FreeFormFinsPrints(_FinsPrints): """Class that contains all free form fins prints.""" +class _FreeFormFinPrints(_FinPrints): + """Class that contains all free form fins prints.""" + + class _TailPrints(_AeroSurfacePrints): """Class that contains all tail prints.""" diff --git a/rocketpy/rocket/__init__.py b/rocketpy/rocket/__init__.py index 463cbe3b3..afb7f0bb6 100644 --- a/rocketpy/rocket/__init__.py +++ b/rocketpy/rocket/__init__.py @@ -2,14 +2,18 @@ from rocketpy.rocket.aero_surface import ( AeroSurface, AirBrakes, + EllipticalFin, EllipticalFins, + Fin, Fins, + FreeFormFin, FreeFormFins, GenericSurface, LinearGenericSurface, NoseCone, RailButtons, Tail, + TrapezoidalFin, TrapezoidalFins, ) from rocketpy.rocket.components import Components diff --git a/rocketpy/rocket/aero_surface/__init__.py b/rocketpy/rocket/aero_surface/__init__.py index ad784f8d0..7634d3500 100644 --- a/rocketpy/rocket/aero_surface/__init__.py +++ b/rocketpy/rocket/aero_surface/__init__.py @@ -1,9 +1,13 @@ from rocketpy.rocket.aero_surface.aero_surface import AeroSurface from rocketpy.rocket.aero_surface.air_brakes import AirBrakes from rocketpy.rocket.aero_surface.fins import ( + EllipticalFin, EllipticalFins, + Fin, Fins, + FreeFormFin, FreeFormFins, + TrapezoidalFin, TrapezoidalFins, ) from rocketpy.rocket.aero_surface.generic_surface import GenericSurface diff --git a/rocketpy/rocket/aero_surface/aero_surface.py b/rocketpy/rocket/aero_surface/aero_surface.py index 15ca14f1d..6727476c6 100644 --- a/rocketpy/rocket/aero_surface/aero_surface.py +++ b/rocketpy/rocket/aero_surface/aero_surface.py @@ -2,6 +2,8 @@ import numpy as np +from rocketpy.mathutils.vector_matrix import Matrix + class AeroSurface(ABC): """Abstract class used to define aerodynamic surfaces.""" @@ -15,6 +17,14 @@ def __init__(self, name, reference_area, reference_length): self.cpy = 0 self.cpz = 0 + self._rotation_surface_to_body = Matrix( + [ + [-1, 0, 0], + [0, 1, 0], + [0, 0, -1], + ] + ) + @staticmethod def _beta(mach): """Defines a parameter that is often used in aerodynamic @@ -130,7 +140,7 @@ def compute_forces_and_moments( R1, R2, R3, M1, M2, M3 = 0, 0, 0, 0, 0, 0 cpz = cp[2] stream_vx, stream_vy, stream_vz = stream_velocity - if stream_vx**2 + stream_vy**2 != 0: # TODO: maybe try/except + if stream_vx**2 + stream_vy**2 != 0: # Normalize component stream velocity in body frame stream_vzn = stream_vz / stream_speed if -1 * stream_vzn < 1: diff --git a/rocketpy/rocket/aero_surface/fins/__init__.py b/rocketpy/rocket/aero_surface/fins/__init__.py index 941aa5465..dd678c625 100644 --- a/rocketpy/rocket/aero_surface/fins/__init__.py +++ b/rocketpy/rocket/aero_surface/fins/__init__.py @@ -1,4 +1,8 @@ +from rocketpy.rocket.aero_surface.fins.elliptical_fin import EllipticalFin from rocketpy.rocket.aero_surface.fins.elliptical_fins import EllipticalFins +from rocketpy.rocket.aero_surface.fins.fin import Fin from rocketpy.rocket.aero_surface.fins.fins import Fins +from rocketpy.rocket.aero_surface.fins.free_form_fin import FreeFormFin from rocketpy.rocket.aero_surface.fins.free_form_fins import FreeFormFins +from rocketpy.rocket.aero_surface.fins.trapezoidal_fin import TrapezoidalFin from rocketpy.rocket.aero_surface.fins.trapezoidal_fins import TrapezoidalFins diff --git a/rocketpy/rocket/aero_surface/fins/_base_fin.py b/rocketpy/rocket/aero_surface/fins/_base_fin.py new file mode 100644 index 000000000..f6b09f797 --- /dev/null +++ b/rocketpy/rocket/aero_surface/fins/_base_fin.py @@ -0,0 +1,344 @@ +import math +from abc import abstractmethod + +import numpy as np + +from rocketpy.mathutils.function import Function + +from ..aero_surface import AeroSurface + + +class _BaseFin(AeroSurface): + """ + Base class for fins, shared by both Fin and Fins classes. + Inherits from AeroSurface. + + Handles shared initialization logic and common properties. + """ + + def __init__( + self, name, rocket_radius, root_chord, span, airfoil=None, cant_angle=0 + ): + """ + Initialize the base fin class. + + Parameters + ---------- + name : str + Name of the fin or fin set. + rocket_radius : float + Rocket radius in meters. + root_chord : float + Root chord of the fin in meters. + span : float + Span of the fin in meters. + airfoil : tuple, optional + Tuple containing airfoil data and unit ('degrees' or 'radians'). + cant_angle : float, optional + Cant angle in degrees. + """ + self.name = name + self._rocket_radius = rocket_radius + self._root_chord = root_chord + self._span = span + self._airfoil = airfoil + self._cant_angle = cant_angle + self._cant_angle_rad = math.radians(cant_angle) + self.geometry = None + + self.reference_area = np.pi * rocket_radius**2 + + super().__init__(name, self.reference_area, self.rocket_diameter) + + def _update_reference_quantities(self): + """Update quantities that depend on rocket radius.""" + self.reference_area = np.pi * self._rocket_radius**2 + self.reference_length = self.rocket_diameter + + def _update_geometry_chain(self): + """Update geometry-dependent quantities in dependency order.""" + self.evaluate_geometrical_parameters() + self.evaluate_center_of_pressure() + self.evaluate_lift_coefficient() + self.evaluate_roll_parameters() + + @property + def rocket_radius(self): + """Rocket radius in meters. + + Returns + ------- + float + Rocket radius in meters. + """ + return self._rocket_radius + + @rocket_radius.setter + def rocket_radius(self, value): + """Set rocket radius and update dependent properties. + + Parameters + ---------- + value : float + Rocket radius in meters. + """ + self._rocket_radius = value + self._update_reference_quantities() + self._update_geometry_chain() + + @property + def rocket_diameter(self): + """Reference rocket diameter in meters.""" + return 2 * self._rocket_radius + + @rocket_diameter.setter + def rocket_diameter(self, value): + """Set reference rocket diameter in meters.""" + self.rocket_radius = value / 2 + + @property + def diameter(self): + """Backward-compatible alias for :attr:`rocket_diameter`.""" + return self.rocket_diameter + + @diameter.setter + def diameter(self, value): + """Backward-compatible alias setter for :attr:`rocket_diameter`.""" + self.rocket_diameter = value + + @property + def d(self): + """Backward-compatible alias for :attr:`rocket_diameter`.""" + return self.rocket_diameter + + @d.setter + def d(self, value): + """Backward-compatible alias setter for :attr:`rocket_diameter`.""" + self.rocket_diameter = value + + @property + def ref_area(self): + """Backward-compatible alias for :attr:`reference_area`.""" + return self.reference_area + + @ref_area.setter + def ref_area(self, value): + """Backward-compatible alias setter for :attr:`reference_area`.""" + self.reference_area = value + + @property + def root_chord(self): + """Root chord length in meters. + + Returns + ------- + float + Root chord length in meters. + """ + return self._root_chord + + @root_chord.setter + def root_chord(self, value): + """Set root chord and update dependent properties. + + Parameters + ---------- + value : float + Root chord length in meters. + """ + self._root_chord = value + self._update_geometry_chain() + self.evaluate_shape() + + @property + def span(self): + """Fin span in meters. + + Returns + ------- + float + Fin span in meters. + """ + return self._span + + @span.setter + def span(self, value): + """Set fin span and update dependent properties. + + Parameters + ---------- + value : float + Fin span in meters. + """ + self._span = value + self._update_geometry_chain() + self.evaluate_shape() + + @property + def cant_angle(self): + """Cant angle in degrees. + + Returns + ------- + float + Cant angle in degrees. + """ + return self._cant_angle + + @cant_angle.setter + def cant_angle(self, value): + """Set cant angle and update radian representation. + + Parameters + ---------- + value : float + Cant angle in degrees. + """ + self._cant_angle = value + self.cant_angle_rad = math.radians(value) + + @property + def cant_angle_rad(self): + """Cant angle in radians. + + Returns + ------- + float + Cant angle in radians. + """ + return self._cant_angle_rad + + @cant_angle_rad.setter + def cant_angle_rad(self, value): + """Set cant angle in radians and update dependent properties. + + Parameters + ---------- + value : float + Cant angle in radians. + """ + self._cant_angle_rad = value + self._update_geometry_chain() + + @property + def airfoil(self): + """Airfoil data for the fin. + + Returns + ------- + tuple or None + Tuple containing airfoil data and unit ('degrees' or 'radians'), + or None if using planar fin. + """ + return self._airfoil + + @airfoil.setter + def airfoil(self, value): + """Set airfoil data and update dependent properties. + + Parameters + ---------- + value : tuple or None + Tuple containing airfoil data and unit ('degrees' or 'radians'), + or None for planar fin. + """ + self._airfoil = value + self._update_geometry_chain() + + def info(self): + """Print fin geometry and lift information.""" + self.prints.geometry() + self.prints.lift() + + def all_info(self): + """Print all available fin information and show all fin plots.""" + self.prints.all() + self.plots.all() + + def evaluate_single_fin_lift_coefficient(self): + """Evaluate the lift coefficient derivative for a single fin. + + Computes the lift coefficient derivative (clalpha) considering the + fin's geometry, airfoil characteristics (if provided), and Mach number + effects using Prandtl-Glauert compressibility correction and + Diederich's planform correlation. + + Sets the `clalpha_single_fin` attribute as a Function of Mach number. + """ + if not self.airfoil: + # Defines clalpha2D as 2*pi for planar fins + clalpha2D_incompressible = 2 * np.pi + else: + # Defines clalpha2D as the derivative of the lift coefficient curve + # for the specific airfoil + self.airfoil_cl = Function( + self.airfoil[0], + title="Airfoil lift coefficient", + interpolation="linear", + ) + + # Differentiating at alpha = 0 to get cl_alpha + clalpha2D_incompressible = self.airfoil_cl.differentiate(x=1e-3, dx=1e-3) + + # Convert to radians if needed + if self.airfoil[1] == "degrees": + clalpha2D_incompressible *= 180 / np.pi + + # Correcting for compressible flow (apply Prandtl-Glauert correction) + clalpha2D = Function(lambda mach: clalpha2D_incompressible / self._beta(mach)) + + # Diederich's Planform Correlation Parameter + planform_correlation_parameter = ( + 2 * np.pi * self.AR / (clalpha2D * np.cos(self.gamma_c)) + ) + + # Lift coefficient derivative for a single fin + def lift_source(mach): + return ( + clalpha2D(mach) + * planform_correlation_parameter(mach) + * (self.Af / self.reference_area) + * np.cos(self.gamma_c) + ) / ( + 2 + + planform_correlation_parameter(mach) + * np.sqrt(1 + (2 / planform_correlation_parameter(mach)) ** 2) + ) + + self.clalpha_single_fin = Function( + lift_source, + "Mach", + "Lift coefficient derivative for a single fin", + ) + + @abstractmethod + def evaluate_lift_coefficient(self): + """Evaluate the lift coefficient for the fin.""" + + @abstractmethod + def evaluate_roll_parameters(self): + """Evaluate roll-related parameters for the fin.""" + + @abstractmethod + def evaluate_center_of_pressure(self): + """Evaluate the center of pressure for the fin.""" + + def evaluate_geometrical_parameters(self): + """Evaluate geometric parameters of the fin. + + This method delegates to the configured geometry strategy. + """ + geometry_data = self.geometry.evaluate_geometrical_parameters() + for key, value in geometry_data.items(): + setattr(self, key, value) + + def evaluate_shape(self): + """Evaluate the shape representation of the fin. + + This method delegates to the configured geometry strategy. + """ + self.shape_vec = self.geometry.evaluate_shape() + + @abstractmethod + def draw(self): + """Draw or render the fin.""" diff --git a/rocketpy/rocket/aero_surface/fins/_geometry.py b/rocketpy/rocket/aero_surface/fins/_geometry.py new file mode 100644 index 000000000..71b213909 --- /dev/null +++ b/rocketpy/rocket/aero_surface/fins/_geometry.py @@ -0,0 +1,564 @@ +"""Geometry strategy classes for fin aerodynamic surfaces.""" + +import warnings +from abc import ABC, abstractmethod + +import numpy as np + + +class _FinGeometry(ABC): + """Base geometry strategy for fin shapes.""" + + def __init__(self, owner): + self.owner = owner + + @abstractmethod + def evaluate_geometrical_parameters(self): + """Evaluate and return geometry-dependent aerodynamic parameters.""" + + @abstractmethod + def evaluate_shape(self): + """Evaluate and return the shape vector used by plotting and outputs.""" + + def get_data(self, include_outputs=False): + """Return geometry-specific serialization data.""" + _ = include_outputs + return {} + + +class _TrapezoidalGeometry(_FinGeometry): + """Geometry strategy for trapezoidal fins.""" + + def __init__( + self, + owner, + tip_chord, + sweep_length=None, + sweep_angle=None, + ): + super().__init__(owner) + if sweep_length is not None and sweep_angle is not None: + raise ValueError("Cannot use sweep_length and sweep_angle together") + + if sweep_angle is not None: + sweep_length = np.tan(np.radians(sweep_angle)) * owner.span + elif sweep_length is None: + sweep_length = owner.root_chord - tip_chord + + self._tip_chord = tip_chord + self._sweep_length = sweep_length + self._sweep_angle = sweep_angle + + @property + def tip_chord(self): + return self._tip_chord + + @tip_chord.setter + def tip_chord(self, value): + self._tip_chord = value + + @property + def sweep_length(self): + return self._sweep_length + + @sweep_length.setter + def sweep_length(self, value): + self._sweep_length = value + + @property + def sweep_angle(self): + return self._sweep_angle + + @sweep_angle.setter + def sweep_angle(self, value): + self._sweep_angle = value + self._sweep_length = np.tan(np.radians(value)) * self.owner.span + + def evaluate_geometrical_parameters(self): + """Calculate trapezoidal fin geometric parameters.""" + # pylint: disable=invalid-name + owner = self.owner + Yr = owner.root_chord + self.tip_chord + Af = Yr * owner.span / 2 + AR = 2 * owner.span**2 / Af + gamma_c = np.arctan( + (self.sweep_length + 0.5 * self.tip_chord - 0.5 * owner.root_chord) + / owner.span + ) + Yma = (owner.span / 3) * (owner.root_chord + 2 * self.tip_chord) / Yr + + tau = (owner.span + owner.rocket_radius) / owner.rocket_radius + lift_interference_factor = 1 + 1 / tau + lambda_ = self.tip_chord / owner.root_chord + + roll_geometrical_constant = ( + (owner.root_chord + 3 * self.tip_chord) * owner.span**3 + + 4 + * (owner.root_chord + 2 * self.tip_chord) + * owner.rocket_radius + * owner.span**2 + + 6 + * (owner.root_chord + self.tip_chord) + * owner.span + * owner.rocket_radius**2 + ) / 12 + roll_damping_interference_factor = 1 + ( + ((tau - lambda_) / tau) - ((1 - lambda_) / (tau - 1)) * np.log(tau) + ) / ( + ((tau + 1) * (tau - lambda_)) / 2 + - ((1 - lambda_) * (tau**3 - 1)) / (3 * (tau - 1)) + ) + roll_forcing_interference_factor = (1 / np.pi**2) * ( + (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) + + (np.pi * (tau**2 + 1) ** 2 / (tau**2 * (tau - 1) ** 2)) + * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + - (2 * np.pi * (tau + 1)) / (tau * (tau - 1)) + + ((tau**2 + 1) ** 2 / (tau**2 * (tau - 1) ** 2)) + * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2 + - (4 * (tau + 1)) + / (tau * (tau - 1)) + * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau)) + ) + + self.Yr = Yr + self.Af = Af + self.AR = AR + self.gamma_c = gamma_c + self.Yma = Yma + self.roll_geometrical_constant = roll_geometrical_constant + self.tau = tau + self.lift_interference_factor = lift_interference_factor + self.λ = lambda_ # pylint: disable=non-ascii-name + self.roll_damping_interference_factor = roll_damping_interference_factor + self.roll_forcing_interference_factor = roll_forcing_interference_factor + + return { + "Yr": Yr, + "Af": Af, + "AR": AR, + "gamma_c": gamma_c, + "Yma": Yma, + "roll_geometrical_constant": roll_geometrical_constant, + "tau": tau, + "lift_interference_factor": lift_interference_factor, + "λ": lambda_, # pylint: disable=non-ascii-name + "roll_damping_interference_factor": roll_damping_interference_factor, + "roll_forcing_interference_factor": roll_forcing_interference_factor, + } + + def evaluate_shape(self): + owner = self.owner + if self.sweep_length: + points = [ + (0, 0), + (self.sweep_length, owner.span), + (self.sweep_length + self.tip_chord, owner.span), + (owner.root_chord, 0), + ] + else: + points = [ + (0, 0), + (owner.root_chord - self.tip_chord, owner.span), + (owner.root_chord, owner.span), + (owner.root_chord, 0), + ] + + x_array, y_array = zip(*points) + shape_vec = [np.array(x_array), np.array(y_array)] + self.shape_vec = shape_vec + return shape_vec + + def get_data(self, include_outputs=False): + data = { + "tip_chord": self.tip_chord, + "sweep_length": self.sweep_length, + "sweep_angle": self.sweep_angle, + } + if include_outputs: + data.update( + { + "shape_vec": getattr(self, "shape_vec", None), + "Af": getattr(self, "Af", None), + "AR": getattr(self, "AR", None), + "gamma_c": getattr(self, "gamma_c", None), + "Yma": getattr(self, "Yma", None), + "roll_geometrical_constant": getattr( + self, "roll_geometrical_constant", None + ), + "tau": getattr(self, "tau", None), + "lift_interference_factor": getattr( + self, "lift_interference_factor", None + ), + "roll_damping_interference_factor": getattr( + self, "roll_damping_interference_factor", None + ), + "roll_forcing_interference_factor": getattr( + self, "roll_forcing_interference_factor", None + ), + } + ) + return data + + +class _EllipticalGeometry(_FinGeometry): + """Geometry strategy for elliptical fins.""" + + def evaluate_geometrical_parameters(self): # pylint: disable=too-many-statements + """Calculate elliptical fin geometric parameters.""" + owner = self.owner + + # pylint: disable=invalid-name + Af = (np.pi * owner.root_chord / 2 * owner.span) / 2 + gamma_c = 0 + AR = 2 * owner.span**2 / Af + Yma = owner.span / (3 * np.pi) * np.sqrt(9 * np.pi**2 - 64) + roll_geometrical_constant = ( + owner.root_chord + * owner.span + * ( + 3 * np.pi * owner.span**2 + + 32 * owner.rocket_radius * owner.span + + 12 * np.pi * owner.rocket_radius**2 + ) + / 48 + ) + + tau = (owner.span + owner.rocket_radius) / owner.rocket_radius + lift_interference_factor = 1 + 1 / tau + if owner.span > owner.rocket_radius: + roll_damping_interference_factor = 1 + ( + owner.rocket_radius**2 + * ( + 2 + * owner.rocket_radius**2 + * np.sqrt(owner.span**2 - owner.rocket_radius**2) + * np.log( + ( + 2 + * owner.span + * np.sqrt(owner.span**2 - owner.rocket_radius**2) + + 2 * owner.span**2 + ) + / owner.rocket_radius + ) + - 2 + * owner.rocket_radius**2 + * np.sqrt(owner.span**2 - owner.rocket_radius**2) + * np.log(2 * owner.span) + + 2 * owner.span**3 + - np.pi * owner.rocket_radius * owner.span**2 + - 2 * owner.rocket_radius**2 * owner.span + + np.pi * owner.rocket_radius**3 + ) + ) / ( + 2 + * owner.span**2 + * (owner.span / 3 + np.pi * owner.rocket_radius / 4) + * (owner.span**2 - owner.rocket_radius**2) + ) + elif owner.span < owner.rocket_radius: + roll_damping_interference_factor = 1 - ( + owner.rocket_radius**2 + * ( + 2 * owner.span**3 + - np.pi * owner.span**2 * owner.rocket_radius + - 2 * owner.span * owner.rocket_radius**2 + + np.pi * owner.rocket_radius**3 + + 2 + * owner.rocket_radius**2 + * np.sqrt(-(owner.span**2) + owner.rocket_radius**2) + * np.arctan( + owner.span / np.sqrt(-(owner.span**2) + owner.rocket_radius**2) + ) + - np.pi + * owner.rocket_radius**2 + * np.sqrt(-(owner.span**2) + owner.rocket_radius**2) + ) + ) / ( + 2 + * owner.span + * (-(owner.span**2) + owner.rocket_radius**2) + * (owner.span**2 / 3 + np.pi * owner.span * owner.rocket_radius / 4) + ) + else: + roll_damping_interference_factor = (28 - 3 * np.pi) / (4 + 3 * np.pi) + + roll_forcing_interference_factor = (1 / np.pi**2) * ( + (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) + + (np.pi * (tau**2 + 1) ** 2 / (tau**2 * (tau - 1) ** 2)) + * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + - (2 * np.pi * (tau + 1)) / (tau * (tau - 1)) + + ((tau**2 + 1) ** 2 / (tau**2 * (tau - 1) ** 2)) + * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2 + - (4 * (tau + 1)) + / (tau * (tau - 1)) + * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau)) + ) + + self.Af = Af + self.AR = AR + self.gamma_c = gamma_c + self.Yma = Yma + self.roll_geometrical_constant = roll_geometrical_constant + self.tau = tau + self.lift_interference_factor = lift_interference_factor + self.roll_damping_interference_factor = roll_damping_interference_factor + self.roll_forcing_interference_factor = roll_forcing_interference_factor + + return { + "Af": Af, + "AR": AR, + "gamma_c": gamma_c, + "Yma": Yma, + "roll_geometrical_constant": roll_geometrical_constant, + "tau": tau, + "lift_interference_factor": lift_interference_factor, + "roll_damping_interference_factor": roll_damping_interference_factor, + "roll_forcing_interference_factor": roll_forcing_interference_factor, + } + + def evaluate_shape(self): + owner = self.owner + angles = np.arange(0, 180, 5) + x_array = owner.root_chord / 2 + owner.root_chord / 2 * np.cos( + np.radians(angles) + ) + y_array = owner.span * np.sin(np.radians(angles)) + shape_vec = [x_array, y_array] + self.shape_vec = shape_vec + return shape_vec + + def get_data(self, include_outputs=False): + if not include_outputs: + return {} + return { + "Af": getattr(self, "Af", None), + "AR": getattr(self, "AR", None), + "gamma_c": getattr(self, "gamma_c", None), + "Yma": getattr(self, "Yma", None), + "roll_geometrical_constant": getattr( + self, "roll_geometrical_constant", None + ), + "tau": getattr(self, "tau", None), + "lift_interference_factor": getattr(self, "lift_interference_factor", None), + "roll_damping_interference_factor": getattr( + self, "roll_damping_interference_factor", None + ), + "roll_forcing_interference_factor": getattr( + self, "roll_forcing_interference_factor", None + ), + } + + +class _FreeFormGeometry(_FinGeometry): + """Geometry strategy for free-form fins.""" + + def __init__(self, owner, shape_points): + super().__init__(owner) + self.shape_points = shape_points + + @staticmethod + def infer_dimensions(shape_points): + """Infer root chord and span from free-form points.""" + down = False + for i in range(1, len(shape_points)): + if shape_points[i][1] > shape_points[i - 1][1] and down: + warnings.warn( + "Jagged fin shape detected. This may cause small " + "inaccuracies center of pressure and pitch moment " + "calculations." + ) + break + if shape_points[i][1] < shape_points[i - 1][1]: + down = True + + root_chord = abs(shape_points[0][0] - shape_points[-1][0]) + ys = [point[1] for point in shape_points] + span = max(ys) - min(ys) + return root_chord, span + + def evaluate_geometrical_parameters( + self, + ): # pylint: disable=too-many-statements,too-many-locals,invalid-name + """Calculate free-form fin geometric parameters.""" + owner = self.owner + + Af = 0 + for i in range(len(self.shape_points) - 1): + x1, y1 = self.shape_points[i] + x2, y2 = self.shape_points[i + 1] + Af += (y1 + y2) * (x1 - x2) + Af = abs(Af) / 2 + if Af < 1e-6: + raise ValueError("Fin area is too small. Check the shape_points.") + + AR = 2 * owner.span**2 / Af + tau = (owner.span + owner.rocket_radius) / owner.rocket_radius + lift_interference_factor = 1 + 1 / tau + + roll_forcing_interference_factor = (1 / np.pi**2) * ( + (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) + + (np.pi * (tau**2 + 1) ** 2 / (tau**2 * (tau - 1) ** 2)) + * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + - (2 * np.pi * (tau + 1)) / (tau * (tau - 1)) + + ((tau**2 + 1) ** 2 / (tau**2 * (tau - 1) ** 2)) + * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2 + - (4 * (tau + 1)) + / (tau * (tau - 1)) + * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau)) + ) + + points_per_line = 40 + chord_lead = np.ones(points_per_line) * np.inf + chord_trail = np.ones(points_per_line) * -np.inf + chord_length = np.zeros(points_per_line) + + for p in range(1, len(self.shape_points)): + x1, y1 = self.shape_points[p - 1] + x2, y2 = self.shape_points[p] + + prev_idx = int(y1 / owner.span * (points_per_line - 1)) + curr_idx = int(y2 / owner.span * (points_per_line - 1)) + prev_idx = np.clip(prev_idx, 0, points_per_line - 1) + curr_idx = np.clip(curr_idx, 0, points_per_line - 1) + + if prev_idx > curr_idx: + prev_idx, curr_idx = curr_idx, prev_idx + + for i in range(prev_idx, curr_idx + 1): + y = i * owner.span / (points_per_line - 1) + if y1 != y2: + x = np.clip( + (y - y2) / (y1 - y2) * x1 + (y1 - y) / (y1 - y2) * x2, + min(x1, x2), + max(x1, x2), + ) + else: + x = x1 + + chord_lead[i] = min(chord_lead[i], x) + chord_trail[i] = max(chord_trail[i], x) + + if y1 < y2: + chord_length[i] -= x + else: + chord_length[i] += x + + invalid_lead = np.isnan(chord_lead) | np.isinf(chord_lead) + invalid_trail = np.isnan(chord_trail) | np.isinf(chord_trail) + chord_lead[invalid_lead | invalid_trail] = 0 + chord_trail[invalid_lead | invalid_trail] = 0 + + chord_length[chord_length < 0] = 0 + chord_length[np.isnan(chord_length)] = 0 + max_chord = chord_trail - chord_lead + chord_length = np.minimum(chord_length, max_chord) + + radius = owner.rocket_radius + total_area = 0 + mac_length = 0 + mac_lead = 0 + mac_span = 0 + cos_gamma_sum = 0 + roll_geometrical_constant = 0 + roll_damping_numerator = 0 + roll_damping_denominator = 0 + + dy = owner.span / (points_per_line - 1) + for i in range(points_per_line): + chord = chord_trail[i] - chord_lead[i] + y = i * dy + + mac_length += chord * chord + mac_span += y * chord + mac_lead += chord_lead[i] * chord + total_area += chord + roll_geometrical_constant += chord_length[i] * (radius + y) ** 2 + roll_damping_numerator += radius**3 * chord / (radius + y) ** 2 + roll_damping_denominator += (radius + y) * chord + + if i > 0: + dx = (chord_trail[i] + chord_lead[i]) / 2 - ( + chord_trail[i - 1] + chord_lead[i - 1] + ) / 2 + cos_gamma_sum += dy / np.hypot(dx, dy) + + mac_length *= dy + mac_span *= dy + mac_lead *= dy + total_area *= dy + roll_geometrical_constant *= dy + roll_damping_numerator *= dy + roll_damping_denominator *= dy + + mac_length /= total_area + mac_span /= total_area + mac_lead /= total_area + cos_gamma = cos_gamma_sum / (points_per_line - 1) + + gamma_c = np.arccos(cos_gamma) + + self.Af = Af + self.AR = AR + self.gamma_c = gamma_c + self.Yma = mac_span + self.mac_length = mac_length + self.mac_lead = mac_lead + self.tau = tau + self.roll_geometrical_constant = roll_geometrical_constant + self.lift_interference_factor = lift_interference_factor + self.roll_forcing_interference_factor = roll_forcing_interference_factor + self.roll_damping_interference_factor = 1 + ( + roll_damping_numerator / roll_damping_denominator + ) + + return { + "Af": Af, + "AR": AR, + "gamma_c": gamma_c, + "Yma": mac_span, + "mac_length": mac_length, + "mac_lead": mac_lead, + "tau": tau, + "roll_geometrical_constant": roll_geometrical_constant, + "lift_interference_factor": lift_interference_factor, + "roll_forcing_interference_factor": roll_forcing_interference_factor, + "roll_damping_interference_factor": self.roll_damping_interference_factor, + } + + def evaluate_shape(self): + x_array, y_array = zip(*self.shape_points) + shape_vec = [np.array(x_array), np.array(y_array)] + self.shape_vec = shape_vec + return shape_vec + + def get_data(self, include_outputs=False): + data = {"shape_points": self.shape_points} + if include_outputs: + data.update( + { + "Af": getattr(self, "Af", None), + "AR": getattr(self, "AR", None), + "gamma_c": getattr(self, "gamma_c", None), + "Yma": getattr(self, "Yma", None), + "mac_length": getattr(self, "mac_length", None), + "mac_lead": getattr(self, "mac_lead", None), + "roll_geometrical_constant": getattr( + self, "roll_geometrical_constant", None + ), + "tau": getattr(self, "tau", None), + "lift_interference_factor": getattr( + self, "lift_interference_factor", None + ), + "roll_forcing_interference_factor": getattr( + self, "roll_forcing_interference_factor", None + ), + "roll_damping_interference_factor": getattr( + self, "roll_damping_interference_factor", None + ), + } + ) + return data diff --git a/rocketpy/rocket/aero_surface/fins/elliptical_fin.py b/rocketpy/rocket/aero_surface/fins/elliptical_fin.py new file mode 100644 index 000000000..f809bca29 --- /dev/null +++ b/rocketpy/rocket/aero_surface/fins/elliptical_fin.py @@ -0,0 +1,193 @@ +from rocketpy.plots.aero_surface_plots import _EllipticalFinPlots +from rocketpy.prints.aero_surface_prints import _EllipticalFinPrints +from rocketpy.rocket.aero_surface.fins._geometry import _EllipticalGeometry +from rocketpy.rocket.aero_surface.fins.fin import Fin + + +class EllipticalFin(Fin): + """Class that defines and holds information for an elliptical fin set. + + This class inherits from the Fin class. + + Note + ---- + Local coordinate system: + - Origin located at the top of the root chord. + - Z axis along the longitudinal axis of symmetry, positive downwards (top -> bottom). + - Y axis perpendicular to the Z axis, in the span direction, positive upwards. + - X axis completes the right-handed coordinate system. + + See Also + -------- + Fin + + Attributes + ---------- + EllipticalFin.rocket_radius : float + The reference rocket radius used for lift coefficient normalization, in + meters. + EllipticalFin.airfoil : tuple + Tuple of two items. First is the airfoil lift curve. + Second is the unit of the curve (radians or degrees) + EllipticalFin.cant_angle : float + Fins cant angle with respect to the rocket centerline, in degrees. + EllipticalFin.cant_angle_rad : float + Fins cant angle with respect to the rocket centerline, in radians. + EllipticalFin.root_chord : float + Fin root chord in meters. + EllipticalFin.span : float + Fin span in meters. + EllipticalFin.name : string + Name of fin set. + EllipticalFin.sweep_length : float + Fins sweep length in meters. By sweep length, understand the axial + distance between the fin root leading edge and the fin tip leading edge + measured parallel to the rocket centerline. + EllipticalFin.sweep_angle : float + Fins sweep angle with respect to the rocket centerline. Must + be given in degrees. + EllipticalFin.rocket_diameter : float + Reference diameter of the rocket, in meters. + EllipticalFin.reference_area : float + Reference area of the rocket. + EllipticalFin.Af : float + Area of the longitudinal section of each fin in the set. + EllipticalFin.AR : float + Aspect ratio of the fin. + EllipticalFin.gamma_c : float + Fin mid-chord sweep angle. + EllipticalFin.Yma : float + Span wise position of the mean aerodynamic chord. + EllipticalFin.roll_geometrical_constant : float + Geometrical constant used in roll calculations. + EllipticalFin.tau : float + Geometrical relation used to simplify lift and roll calculations. + EllipticalFin.lift_interference_factor : float + Factor of Fin-Body interference in the lift coefficient. + EllipticalFin.cp : tuple + Tuple with the x, y and z local coordinates of the fin set center of + pressure. Has units of length and is given in meters. + EllipticalFin.cpx : float + Fin set local center of pressure x coordinate. Has units of length and + is given in meters. + EllipticalFin.cpy : float + Fin set local center of pressure y coordinate. Has units of length and + is given in meters. + EllipticalFin.cpz : float + Fin set local center of pressure z coordinate. Has units of length and + is given in meters. + EllipticalFin.cl : Function + Function which defines the lift coefficient as a function of the angle + of attack and the Mach number. Takes as input the angle of attack in + radians and the Mach number. Returns the lift coefficient. + EllipticalFin.clalpha : float + Lift coefficient slope. Has units of 1/rad. + """ + + def __init__( + self, + angular_position, + root_chord, + span, + rocket_radius, + cant_angle=0, + airfoil=None, + name="Elliptical Fin", + ): + """Initialize EllipticalFin class. + + Parameters + ---------- + angular_position : float + Angular position of the fin in degrees measured as the rotation + around the symmetry axis of the rocket relative to one of the other + principal axis. See :ref:`Angular Position Inputs ` + root_chord : int, float + Fin root chord in meters. + span : int, float + Fin span in meters. + rocket_radius : int, float + Reference radius to calculate lift coefficient, in meters. + cant_angle : int, float, optional + Fins cant angle with respect to the rocket centerline. Must + be given in degrees. + sweep_length : int, float, optional + Fins sweep length in meters. By sweep length, understand the axial + distance between the fin root leading edge and the fin tip leading + edge measured parallel to the rocket centerline. If not given, the + sweep length is assumed to be equal the root chord minus the tip + chord, in which case the fin is a right trapezoid with its base + perpendicular to the rocket's axis. Cannot be used in conjunction + with sweep_angle. + sweep_angle : int, float, optional + Fins sweep angle with respect to the rocket centerline. Must + be given in degrees. If not given, the sweep angle is automatically + calculated, in which case the fin is assumed to be a right trapezoid + with its base perpendicular to the rocket's axis. + Cannot be used in conjunction with sweep_length. + airfoil : tuple, optional + Default is null, in which case fins will be treated as flat plates. + Otherwise, if tuple, fins will be considered as airfoils. The + tuple's first item specifies the airfoil's lift coefficient + by angle of attack and must be either a .csv, .txt, ndarray + or callable. The .csv and .txt files can contain a single line + header and the first column must specify the angle of attack, while + the second column must specify the lift coefficient. The + ndarray should be as [(x0, y0), (x1, y1), (x2, y2), ...] + where x0 is the angle of attack and y0 is the lift coefficient. + If callable, it should take an angle of attack as input and + return the lift coefficient at that angle of attack. + The tuple's second item is the unit of the angle of attack, + accepting either "radians" or "degrees". + name : str + Name of elliptical fin. + + Returns + ------- + None + """ + + super().__init__( + angular_position, + root_chord, + span, + rocket_radius, + cant_angle, + airfoil, + name, + ) + + self.geometry = _EllipticalGeometry(self) + self._update_geometry_chain() + self.evaluate_shape() + + self.prints = _EllipticalFinPrints(self) + self.plots = _EllipticalFinPlots(self) + + def evaluate_center_of_pressure(self): + """Calculates and returns the center of pressure of the fin in local + coordinates. The center of pressure position is saved and stored as a + tuple.""" + # Barrowman elliptical-fin center of pressure location. + cpz = 0.288 * self.root_chord + self.cpx = 0 + self.cpy = self.Yma + self.cpz = cpz + self.cp = (self.cpx, self.cpy, self.cpz) + + def to_dict(self, include_outputs=False): + data = super().to_dict(include_outputs=include_outputs) + data.update(self.geometry.get_data(include_outputs=include_outputs)) + return data + + @classmethod + def from_dict(cls, data): + return cls( + angular_position=data["angular_position"], + root_chord=data["root_chord"], + span=data["span"], + rocket_radius=data["rocket_radius"], + cant_angle=data["cant_angle"], + airfoil=data["airfoil"], + name=data["name"], + ) diff --git a/rocketpy/rocket/aero_surface/fins/elliptical_fins.py b/rocketpy/rocket/aero_surface/fins/elliptical_fins.py index 61d98bc0d..4576bd1f3 100644 --- a/rocketpy/rocket/aero_surface/fins/elliptical_fins.py +++ b/rocketpy/rocket/aero_surface/fins/elliptical_fins.py @@ -1,7 +1,6 @@ -import numpy as np - from rocketpy.plots.aero_surface_plots import _EllipticalFinsPlots from rocketpy.prints.aero_surface_prints import _EllipticalFinsPrints +from rocketpy.rocket.aero_surface.fins._geometry import _EllipticalGeometry from .fins import Fins @@ -35,9 +34,6 @@ class EllipticalFins(Fins): Second is the unit of the curve (radians or degrees) EllipticalFins.cant_angle : float Fins cant angle with respect to the rocket centerline, in degrees. - EllipticalFins.changing_attribute_dict : dict - Dictionary that stores the name and the values of the attributes that - may be changed during a simulation. Useful for control systems. EllipticalFins.cant_angle_rad : float Fins cant angle with respect to the rocket centerline, in radians. EllipticalFins.root_chord : float @@ -53,9 +49,9 @@ class EllipticalFins(Fins): EllipticalFins.sweep_angle : float Fins sweep angle with respect to the rocket centerline. Must be given in degrees. - EllipticalFins.d : float + EllipticalFins.rocket_diameter : float Reference diameter of the rocket, in meters. - EllipticalFins.ref_area : float + EllipticalFins.reference_area : float Reference area of the rocket. EllipticalFins.Af : float Area of the longitudinal section of each fin in the set. @@ -162,10 +158,9 @@ def __init__( name, ) - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() + self.geometry = _EllipticalGeometry(self) + self._update_geometry_chain() + self.evaluate_shape() self.prints = _EllipticalFinsPrints(self) self.plots = _EllipticalFinsPlots(self) @@ -179,160 +174,18 @@ def evaluate_center_of_pressure(self): ------- None """ - # Center of pressure position in local coordinates + # Barrowman elliptical-fin center of pressure location. cpz = 0.288 * self.root_chord self.cpx = 0 self.cpy = 0 self.cpz = cpz self.cp = (self.cpx, self.cpy, self.cpz) - def evaluate_geometrical_parameters(self): # pylint: disable=too-many-statements - """Calculates and saves fin set's geometrical parameters such as the - fins' area, aspect ratio and parameters for roll movement. - - Returns - ------- - None - """ - - # Compute auxiliary geometrical parameters - # pylint: disable=invalid-name - Af = (np.pi * self.root_chord / 2 * self.span) / 2 # Fin area - gamma_c = 0 # Zero for elliptical fins - AR = 2 * self.span**2 / Af # Fin aspect ratio - Yma = ( - self.span / (3 * np.pi) * np.sqrt(9 * np.pi**2 - 64) - ) # Span wise coord of mean aero chord - roll_geometrical_constant = ( - self.root_chord - * self.span - * ( - 3 * np.pi * self.span**2 - + 32 * self.rocket_radius * self.span - + 12 * np.pi * self.rocket_radius**2 - ) - / 48 - ) - - # Fin–body interference correction parameters - tau = (self.span + self.rocket_radius) / self.rocket_radius - lift_interference_factor = 1 + 1 / tau - if self.span > self.rocket_radius: - roll_damping_interference_factor = 1 + ( - (self.rocket_radius**2) - * ( - 2 - * (self.rocket_radius**2) - * np.sqrt(self.span**2 - self.rocket_radius**2) - * np.log( - ( - 2 - * self.span - * np.sqrt(self.span**2 - self.rocket_radius**2) - + 2 * self.span**2 - ) - / self.rocket_radius - ) - - 2 - * (self.rocket_radius**2) - * np.sqrt(self.span**2 - self.rocket_radius**2) - * np.log(2 * self.span) - + 2 * self.span**3 - - np.pi * self.rocket_radius * self.span**2 - - 2 * (self.rocket_radius**2) * self.span - + np.pi * self.rocket_radius**3 - ) - ) / ( - 2 - * (self.span**2) - * (self.span / 3 + np.pi * self.rocket_radius / 4) - * (self.span**2 - self.rocket_radius**2) - ) - elif self.span < self.rocket_radius: - roll_damping_interference_factor = 1 - ( - self.rocket_radius**2 - * ( - 2 * self.span**3 - - np.pi * self.span**2 * self.rocket_radius - - 2 * self.span * self.rocket_radius**2 - + np.pi * self.rocket_radius**3 - + 2 - * self.rocket_radius**2 - * np.sqrt(-(self.span**2) + self.rocket_radius**2) - * np.arctan( - (self.span) / (np.sqrt(-(self.span**2) + self.rocket_radius**2)) - ) - - np.pi - * self.rocket_radius**2 - * np.sqrt(-(self.span**2) + self.rocket_radius**2) - ) - ) / ( - 2 - * self.span - * (-(self.span**2) + self.rocket_radius**2) - * (self.span**2 / 3 + np.pi * self.span * self.rocket_radius / 4) - ) - else: - roll_damping_interference_factor = (28 - 3 * np.pi) / (4 + 3 * np.pi) - - roll_forcing_interference_factor = (1 / np.pi**2) * ( - (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) - + ((np.pi * (tau**2 + 1) ** 2) / (tau**2 * (tau - 1) ** 2)) - * np.arcsin((tau**2 - 1) / (tau**2 + 1)) - - (2 * np.pi * (tau + 1)) / (tau * (tau - 1)) - + ((tau**2 + 1) ** 2) - / (tau**2 * (tau - 1) ** 2) - * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2 - - (4 * (tau + 1)) - / (tau * (tau - 1)) - * np.arcsin((tau**2 - 1) / (tau**2 + 1)) - + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau)) - ) - - # Store values - # pylint: disable=invalid-name - self.Af = Af # Fin area - self.AR = AR # Fin aspect ratio - self.gamma_c = gamma_c # Mid chord angle - self.Yma = Yma # Span wise coord of mean aero chord - self.roll_geometrical_constant = roll_geometrical_constant - self.tau = tau - self.lift_interference_factor = lift_interference_factor - self.roll_damping_interference_factor = roll_damping_interference_factor - self.roll_forcing_interference_factor = roll_forcing_interference_factor - - self.evaluate_shape() - - def evaluate_shape(self): - angles = np.arange(0, 180, 5) - x_array = self.root_chord / 2 + self.root_chord / 2 * np.cos(np.radians(angles)) - y_array = self.span * np.sin(np.radians(angles)) - self.shape_vec = [x_array, y_array] - - def info(self): - self.prints.geometry() - self.prints.lift() - - def all_info(self): - self.prints.all() - self.plots.all() - def to_dict(self, **kwargs): data = super().to_dict(**kwargs) - if kwargs.get("include_outputs", False): - data.update( - { - "Af": self.Af, - "AR": self.AR, - "gamma_c": self.gamma_c, - "Yma": self.Yma, - "roll_geometrical_constant": self.roll_geometrical_constant, - "tau": self.tau, - "lift_interference_factor": self.lift_interference_factor, - "roll_damping_interference_factor": self.roll_damping_interference_factor, - "roll_forcing_interference_factor": self.roll_forcing_interference_factor, - } - ) + data.update( + self.geometry.get_data(include_outputs=kwargs.get("include_outputs", False)) + ) return data @classmethod diff --git a/rocketpy/rocket/aero_surface/fins/fin.py b/rocketpy/rocket/aero_surface/fins/fin.py new file mode 100644 index 000000000..5fec8a099 --- /dev/null +++ b/rocketpy/rocket/aero_surface/fins/fin.py @@ -0,0 +1,464 @@ +import math + +import numpy as np + +from rocketpy.mathutils.function import Function +from rocketpy.mathutils.vector_matrix import Matrix, Vector +from rocketpy.rocket.aero_surface.fins._base_fin import _BaseFin + + +class Fin(_BaseFin): + """Abstract class that holds common methods for the individual fin classes. + Cannot be instantiated. + + Note + ---- + Local coordinate system: + - Origin located at the top of the root chord. + - Z axis along the longitudinal axis of symmetry, positive downwards (top -> bottom). + - Y axis perpendicular to the Z axis, in the span direction, positive upwards. + - X axis completes the right-handed coordinate system. + + Attributes + ---------- + Fin.rocket_radius : float + The reference rocket radius used for lift coefficient normalization, + in meters. + Fin.airfoil : tuple + Tuple of two items. First is the airfoil lift curve. + Second is the unit of the curve (radians or degrees). + Fin.cant_angle : float + Fin cant angle with respect to the rocket centerline, in degrees. + Fin.changing_attribute_dict : dict + Dictionary that stores the name and the values of the attributes that + may be changed during a simulation. Useful for control systems. + Fin.cant_angle_rad : float + Fin cant angle with respect to the rocket centerline, in radians. + Fin.root_chord : float + Fin root chord in meters. + Fin.tip_chord : float + Fin tip chord in meters. + Fin.span : float + Fin span in meters. + Fin.name : string + Name of fin set. + Fin.sweep_length : float + Fin sweep length in meters. By sweep length, understand the axial + distance between the fin root leading edge and the fin tip leading edge + measured parallel to the rocket centerline. + Fin.sweep_angle : float + Fin sweep angle with respect to the rocket centerline. Must + be given in degrees. + Fin.rocket_diameter : float + Reference diameter of the rocket. Has units of length and is given + in meters. + Fin.reference_area : float + Reference area of the rocket. + Fin.Af : float + Area of the longitudinal section of each fin in the set. + Fin.AR : float + Aspect ratio of each fin in the set. + Fin.gamma_c : float + Fin mid-chord sweep angle. + Fin.Yma : float + Span wise position of the mean aerodynamic chord. + Fin.roll_geometrical_constant : float + Geometrical constant used in roll calculations. + Fin.tau : float + Geometrical relation used to simplify lift and roll calculations. + Fin.lift_interference_factor : float + Factor of Fin-Body interference in the lift coefficient. + Fin.cp : tuple + Tuple with the x, y and z local coordinates of the fin set center of + pressure. Has units of length and is given in meters. + Fin.cpx : float + Fin set local center of pressure x coordinate. Has units of length and + is given in meters. + Fin.cpy : float + Fin set local center of pressure y coordinate. Has units of length and + is given in meters. + Fin.cpz : float + Fin set local center of pressure z coordinate. Has units of length and + is given in meters. + Fin.cl : Function + Function which defines the lift coefficient as a function of the angle + of attack and the Mach number. Takes as input the angle of attack in + radians and the Mach number. Returns the lift coefficient. + Fin.clalpha : float + Lift coefficient slope. Has units of 1/rad. + Fin.roll_parameters : list + List containing the roll moment lift coefficient, the roll moment + damping coefficient and the cant angle in radians. + """ + + def __init__( + self, + angular_position, + root_chord, + span, + rocket_radius, + cant_angle=0, + airfoil=None, + name="Fin", + ): + """Initialize Fin class. + + Parameters + ---------- + angular_position : float + Angular position of the fin in degrees measured as the rotation + around the symmetry axis of the rocket relative to one of the other + principal axis. See :ref:`Angular Position Inputs ` + root_chord : int, float + Fin root chord in meters. + span : int, float + Fin span in meters. + rocket_radius : int, float + Reference rocket radius used for lift coefficient normalization. + cant_angle : int, float, optional + Fin cant angle with respect to the rocket centerline. Must + be given in degrees. + airfoil : tuple, optional + Default is null, in which case fins will be treated as flat plates. + Otherwise, if tuple, fins will be considered as airfoils. The + tuple's first item specifies the airfoil's lift coefficient + by angle of attack and must be either a .csv, .txt, ndarray + or callable. The .csv and .txt files can contain a single line + header and the first column must specify the angle of attack, while + the second column must specify the lift coefficient. The + ndarray should be as [(x0, y0), (x1, y1), (x2, y2), ...] + where x0 is the angle of attack and y0 is the lift coefficient. + If callable, it should take an angle of attack as input and + return the lift coefficient at that angle of attack. + The tuple's second item is the unit of the angle of attack, + accepting either "radians" or "degrees". + name : str + Name of fin. + """ + super().__init__( + name=name, + rocket_radius=rocket_radius, + root_chord=root_chord, + span=span, + airfoil=airfoil, + cant_angle=cant_angle, + ) + + # Store values + self._angular_position = angular_position + self._angular_position_rad = math.radians(angular_position) + + @property + def cant_angle(self): + return self._cant_angle + + @cant_angle.setter + def cant_angle(self, value): + self._cant_angle = value + self.cant_angle_rad = math.radians(value) + + @property + def cant_angle_rad(self): + return self._cant_angle_rad + + @cant_angle_rad.setter + def cant_angle_rad(self, value): + self._cant_angle_rad = value + self.evaluate_geometrical_parameters() + self.evaluate_center_of_pressure() + self.evaluate_lift_coefficient() + self.evaluate_roll_parameters() + self.evaluate_rotation_matrix() + + @property + def angular_position(self): + return self._angular_position + + @angular_position.setter + def angular_position(self, value): + self._angular_position = value + self.angular_position_rad = math.radians(value) + + @property + def angular_position_rad(self): + return self._angular_position_rad + + @angular_position_rad.setter + def angular_position_rad(self, value): + self._angular_position_rad = value + self.evaluate_rotation_matrix() + + def evaluate_lift_coefficient(self): + """Calculates and returns the fin set's lift coefficient. + The lift coefficient is saved and returned. This function + also calculates and saves the lift coefficient derivative + for a single fin and the lift coefficient derivative for + a number of n fins corrected for Fin-Body interference. + + Returns + ------- + None + """ + self.evaluate_single_fin_lift_coefficient() + + self.clalpha = self.clalpha_single_fin * self.lift_interference_factor + + # Cl = clalpha * alpha + self.cl = Function( + lambda alpha, mach: alpha * self.clalpha(mach), + ["Alpha (rad)", "Mach"], + "Lift coefficient", + ) + + return self.cl + + def evaluate_roll_parameters(self): + """Calculates and returns the fin set's roll coefficients. + The roll coefficients are saved in a list. + + Returns + ------- + self.roll_parameters : list + List containing the roll moment lift coefficient, the + roll moment damping coefficient and the cant angle in + radians + """ + clf_delta = ( + self.roll_forcing_interference_factor + * (self.Yma + self.rocket_radius) + * self.clalpha_single_fin + / self.reference_length + ) # Function of mach number + clf_delta.set_inputs("Mach") + clf_delta.set_outputs("Roll moment forcing coefficient derivative") + clf_delta.set_title( + "Roll moment forcing coefficient derivative vs. Mach number" + ) + cld_omega = -( + 2 + * self.roll_damping_interference_factor + * self.clalpha_single_fin + * np.cos(self.cant_angle_rad) + * self.roll_geometrical_constant + / (self.reference_area * self.reference_length**2) + ) # Function of mach number + cld_omega.set_inputs("Mach") + cld_omega.set_outputs("Roll moment damping coefficient derivative") + cld_omega.set_title( + "Roll moment damping coefficient derivative vs. Mach number" + ) + self.roll_parameters = [clf_delta, cld_omega, self.cant_angle_rad] + return self.roll_parameters + + def evaluate_rotation_matrix(self): + """Calculates and returns the rotation matrix from the rocket body frame + to the fin frame. + + Note + ---- + Local coordinate system: + + - Origin located at the leading edge of the root chord. + - Z axis along the longitudinal axis of the fin, positive + downwards (leading edge -> trailing edge). + - Y axis perpendicular to the Z axis, in the span direction, + positive upwards (root chord -> tip chord). + - X axis completes the right-handed coordinate system. + + + Returns + ------- + None + + References + ---------- + :ref:`Individual Fin Model ` + """ + phi = self.angular_position_rad + delta = self.cant_angle_rad + sin_phi = math.sin(phi) + cos_phi = math.cos(phi) + sin_delta = math.sin(delta) + cos_delta = math.cos(delta) + + # Rotation about body Z by angular position + R_phi = Matrix( + [ + [cos_phi, -sin_phi, 0], + [sin_phi, cos_phi, 0], + [0, 0, 1], + ] + ) + + # Cant rotation about body Y + R_delta = Matrix( + [ + [cos_delta, 0, -sin_delta], + [0, 1, 0], + [sin_delta, 0, cos_delta], + ] + ) + + # 180 flip about Y to align fin leading/trailing edge + R_pi = Matrix( + [ + [-1, 0, 0], + [0, 1, 0], + [0, 0, -1], + ] + ) + + # Uncanted body to fin, then apply cant + R_uncanted = R_phi @ R_pi + R_body_to_fin = R_delta @ R_uncanted + + # Store for downstream transforms + self._rotation_fin_to_body_uncanted = R_uncanted.transpose + self._rotation_body_to_fin = R_body_to_fin + self._rotation_fin_to_body = R_body_to_fin.transpose + self._rotation_surface_to_body = self._rotation_fin_to_body + + def compute_forces_and_moments( + self, + stream_velocity, + stream_speed, + stream_mach, + rho, + cp, + omega, + *args, + ): # pylint: disable=arguments-differ + """Computes the forces and moments acting on the aerodynamic surface. + + Parameters + ---------- + stream_velocity : tuple of float + The velocity of the airflow relative to the surface. + stream_speed : float + The magnitude of the airflow speed. + stream_mach : float + The Mach number of the airflow. + rho : float + Air density. + cp : Vector + Center of pressure coordinates in the body frame. + omega: tuple[float, float, float] + Tuple containing angular velocities around the x, y, z axes. + + Returns + ------- + tuple of float + The aerodynamic forces (lift, side_force, drag) and moments + (pitch, yaw, roll) in the body frame. + """ + R1, R2, R3, M1, M2, M3 = 0, 0, 0, 0, 0, 0 + + # stream velocity in fin frame + stream_velocity_f = self._rotation_body_to_fin @ stream_velocity + + attack_angle = np.arctan2(stream_velocity_f[0], stream_velocity_f[2]) + # Force in the X direction of the fin + X = ( + 0.5 + * rho + * stream_speed**2 + * self.reference_area + * self.cl.get_value_opt(attack_angle, stream_mach) + ) + # Force in body frame + R1, R2, R3 = self._rotation_fin_to_body @ Vector([X, 0, 0]) + # Moments + M1, M2, M3 = cp ^ Vector([R1, R2, R3]) + # Apply roll interference factor, disregarding lift interference factor + M3 *= self.roll_forcing_interference_factor / self.lift_interference_factor + + # Roll damping + _, cld_omega, _ = self.roll_parameters + M3_damping = ( + (1 / 2 * rho * stream_speed) + * self.reference_area + * (self.reference_length) ** 2 + * cld_omega.get_value_opt(stream_mach) + * omega[2] # omega3 + / 2 + ) + M3 += M3_damping + return R1, R2, R3, M1, M2, M3 + + def _compute_leading_edge_position(self, position, _csys): + """Computes the position of the fin leading edge in a rocket's user, + given its position in a rocket.""" + # Point from deflection from cant angle in the plane perpendicular to + # the fuselage where the fin is located in the fin frame + p = Vector( + [ + -self.root_chord / 2 * np.sin(self.cant_angle_rad), + 0, + self.root_chord / 2 * (1 - np.cos(self.cant_angle_rad)), + ] + ) + # Rotate the point to the body frame orientation + p = self._rotation_fin_to_body_uncanted @ p + + # Rotate the point to the user-defined coordinate system + p = Vector([p.x * _csys, p.y, p.z * _csys]) + + # Build the leading-edge position in the user frame as if no cant + # angle was applied. Scalars are interpreted as z coordinates only, + # while vectors/tuples/lists are interpreted as full (x, y, z) + # coordinates. + if isinstance(position, (Vector, tuple, list)): + position = Vector(position) + else: + position = Vector( + [ + -self.rocket_radius * math.sin(self.angular_position_rad) * _csys, + self.rocket_radius * math.cos(self.angular_position_rad), + position, + ] + ) + + # Translate the position of the fin leading edge to the position of the + # fin leading edge with cant angle + position += p + return position + + def to_dict(self, include_outputs=False): + data = { + "angular_position": self.angular_position, + "root_chord": self.root_chord, + "span": self.span, + "rocket_radius": self.rocket_radius, + "cant_angle": self.cant_angle, + "airfoil": self.airfoil, + "name": self.name, + } + + if include_outputs: + data.update( + { + "cp": self.cp, + "cl": self.cl, + "roll_parameters": self.roll_parameters, + "rocket_diameter": self.rocket_diameter, + "diameter": self.rocket_diameter, + "d": self.rocket_diameter, + "reference_area": self.reference_area, + "ref_area": self.reference_area, + } + ) + return data + + def draw(self, *, filename=None): + """Draw the fin shape along with some important information, including + the center line, the quarter line and the center of pressure position. + + Parameters + ---------- + filename : str | None, optional + The path the plot should be saved to. By default None, in which case + the plot will be shown instead of saved. Supported file endings are: + eps, jpg, jpeg, pdf, pgf, png, ps, raw, rgba, svg, svgz, tif, tiff + and webp (these are the formats supported by matplotlib). + """ + self.plots.draw(filename=filename) diff --git a/rocketpy/rocket/aero_surface/fins/fins.py b/rocketpy/rocket/aero_surface/fins/fins.py index fe24a84f4..b61bfcc19 100644 --- a/rocketpy/rocket/aero_surface/fins/fins.py +++ b/rocketpy/rocket/aero_surface/fins/fins.py @@ -1,11 +1,10 @@ import numpy as np from rocketpy.mathutils.function import Function +from rocketpy.rocket.aero_surface.fins._base_fin import _BaseFin -from ..aero_surface import AeroSurface - -class Fins(AeroSurface): +class Fins(_BaseFin): """Abstract class that holds common methods for the fin classes. Cannot be instantiated. @@ -49,10 +48,10 @@ class Fins(AeroSurface): Fins.sweep_angle : float Fins sweep angle with respect to the rocket centerline. Must be given in degrees. - Fins.d : float + Fins.rocket_diameter : float Reference diameter of the rocket. Has units of length and is given in meters. - Fins.ref_area : float + Fins.reference_area : float Reference area of the rocket. Fins.Af : float Area of the longitudinal section of each fin in the set. @@ -132,27 +131,18 @@ def __init__( accepting either "radians" or "degrees". name : str Name of fin set. - - Returns - ------- - None """ - # Compute auxiliary geometrical parameters - d = 2 * rocket_radius - ref_area = np.pi * rocket_radius**2 # Reference area - - super().__init__(name, ref_area, d) + super().__init__( + name=name, + rocket_radius=rocket_radius, + root_chord=root_chord, + span=span, + airfoil=airfoil, + cant_angle=-cant_angle, + ) # Store values self._n = n - self._rocket_radius = rocket_radius - self._airfoil = airfoil - self._cant_angle = cant_angle - self._root_chord = root_chord - self._span = span - self.name = name - self.d = d - self.ref_area = ref_area # Reference area @property def n(self): @@ -166,134 +156,26 @@ def n(self, value): self.evaluate_lift_coefficient() self.evaluate_roll_parameters() - @property - def root_chord(self): - return self._root_chord - - @root_chord.setter - def root_chord(self, value): - self._root_chord = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - - @property - def span(self): - return self._span - - @span.setter - def span(self, value): - self._span = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - - @property - def rocket_radius(self): - return self._rocket_radius - - @rocket_radius.setter - def rocket_radius(self, value): - self._rocket_radius = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - - @property - def cant_angle(self): - return self._cant_angle - - @cant_angle.setter - def cant_angle(self, value): - self._cant_angle = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - - @property - def airfoil(self): - return self._airfoil - - @airfoil.setter - def airfoil(self, value): - self._airfoil = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() - def evaluate_lift_coefficient(self): """Calculates and returns the fin set's lift coefficient. The lift coefficient is saved and returned. This function also calculates and saves the lift coefficient derivative for a single fin and the lift coefficient derivative for a number of n fins corrected for Fin-Body interference. - - Returns - ------- - None """ - if not self.airfoil: - # Defines clalpha2D as 2*pi for planar fins - clalpha2D_incompressible = 2 * np.pi - else: - # Defines clalpha2D as the derivative of the lift coefficient curve - # for the specific airfoil - self.airfoil_cl = Function( - self.airfoil[0], - interpolation="linear", - ) - - # Differentiating at alpha = 0 to get cl_alpha - clalpha2D_incompressible = self.airfoil_cl.differentiate_complex_step( - x=1e-3, dx=1e-3 - ) - - # Convert to radians if needed - if self.airfoil[1] == "degrees": - clalpha2D_incompressible *= 180 / np.pi - - # Correcting for compressible flow (apply Prandtl-Glauert correction) - clalpha2D = Function(lambda mach: clalpha2D_incompressible / self._beta(mach)) - - # Diederich's Planform Correlation Parameter - planform_correlation_parameter = ( - 2 * np.pi * self.AR / (clalpha2D * np.cos(self.gamma_c)) - ) - - # Lift coefficient derivative for a single fin - def lift_source(mach): - return ( - clalpha2D(mach) - * planform_correlation_parameter(mach) - * (self.Af / self.ref_area) - * np.cos(self.gamma_c) - ) / ( - 2 - + planform_correlation_parameter(mach) - * np.sqrt(1 + (2 / planform_correlation_parameter(mach)) ** 2) - ) - - self.clalpha_single_fin = Function( - lift_source, - "Mach", - "Lift coefficient derivative for a single fin", - ) + self.evaluate_single_fin_lift_coefficient() # Lift coefficient derivative for n fins corrected with Fin-Body interference self.clalpha_multiple_fins = ( - self.lift_interference_factor - * self.fin_num_correction(self.n) + self.fin_num_correction(self.n) + * self.lift_interference_factor * self.clalpha_single_fin ) # Function of mach number self.clalpha_multiple_fins.set_inputs("Mach") self.clalpha_multiple_fins.set_outputs( f"Lift coefficient derivative for {self.n:.0f} fins" ) + self.clalpha = self.clalpha_multiple_fins # Cl = clalpha * alpha @@ -316,31 +198,32 @@ def evaluate_roll_parameters(self): roll moment damping coefficient and the cant angle in radians """ - - self.cant_angle_rad = np.radians(self.cant_angle) - clf_delta = ( self.roll_forcing_interference_factor - * self.n + * self.fin_num_correction(self.n) * (self.Yma + self.rocket_radius) * self.clalpha_single_fin - / self.d + / self.reference_length ) # Function of mach number clf_delta.set_inputs("Mach") clf_delta.set_outputs("Roll moment forcing coefficient derivative") - clf_delta.set_title(None) - cld_omega = ( + clf_delta.set_title( + "Roll moment forcing coefficient derivative vs. Mach number" + ) + cld_omega = -( 2 * self.roll_damping_interference_factor * self.n * self.clalpha_single_fin * np.cos(self.cant_angle_rad) * self.roll_geometrical_constant - / (self.ref_area * self.d**2) + / (self.reference_area * self.reference_length**2) ) # Function of mach number cld_omega.set_inputs("Mach") cld_omega.set_outputs("Roll moment damping coefficient derivative") - cld_omega.set_title(None) + cld_omega.set_title( + "Roll moment damping coefficient derivative vs. Mach number" + ) self.roll_parameters = [clf_delta, cld_omega, self.cant_angle_rad] return self.roll_parameters @@ -425,7 +308,7 @@ def compute_forces_and_moments( * omega[2] / 2 ) - M3 = M3_forcing - M3_damping + M3 = M3_forcing + M3_damping return R1, R2, R3, M1, M2, M3 def to_dict(self, **kwargs): @@ -463,8 +346,11 @@ def to_dict(self, **kwargs): "cp": self.cp, "cl": cl, "roll_parameters": self.roll_parameters, - "d": self.d, - "ref_area": self.ref_area, + "rocket_diameter": self.rocket_diameter, + "diameter": self.rocket_diameter, + "d": self.rocket_diameter, + "reference_area": self.reference_area, + "ref_area": self.reference_area, } ) @@ -481,9 +367,5 @@ def draw(self, *, filename=None): the plot will be shown instead of saved. Supported file endings are: eps, jpg, jpeg, pdf, pgf, png, ps, raw, rgba, svg, svgz, tif, tiff and webp (these are the formats supported by matplotlib). - - Returns - ------- - None """ self.plots.draw(filename=filename) diff --git a/rocketpy/rocket/aero_surface/fins/free_form_fin.py b/rocketpy/rocket/aero_surface/fins/free_form_fin.py new file mode 100644 index 000000000..5099d322b --- /dev/null +++ b/rocketpy/rocket/aero_surface/fins/free_form_fin.py @@ -0,0 +1,189 @@ +from rocketpy.plots.aero_surface_plots import _FreeFormFinPlots +from rocketpy.prints.aero_surface_prints import _FreeFormFinPrints +from rocketpy.rocket.aero_surface.fins._geometry import _FreeFormGeometry +from rocketpy.rocket.aero_surface.fins.fin import Fin + + +class FreeFormFin(Fin): + """Class that defines and holds information for a free form fin set. + + This class inherits from the Fin class. + + Note + ---- + Local coordinate system: + - Origin located at the top of the root chord. + - Z axis along the longitudinal axis of symmetry, positive downwards (top -> bottom). + - Y axis perpendicular to the Z axis, in the span direction, positive upwards. + - X axis completes the right-handed coordinate system. + + See Also + -------- + Fin + + Attributes + ---------- + FreeFormFin.n : int + Number of fins in fin set. + FreeFormFin.rocket_radius : float + The reference rocket radius used for lift coefficient normalization, in + meters. + FreeFormFin.airfoil : tuple + Tuple of two items. First is the airfoil lift curve. + Second is the unit of the curve (radians or degrees). + FreeFormFin.cant_angle : float + Fins cant angle with respect to the rocket centerline, in degrees. + FreeFormFin.cant_angle_rad : float + Fins cant angle with respect to the rocket centerline, in radians. + FreeFormFin.root_chord : float + Fin root chord in meters. + FreeFormFin.span : float + Fin span in meters. + FreeFormFin.name : string + Name of fin set. + FreeFormFin.rocket_diameter : float + Reference diameter of the rocket, in meters. + FreeFormFin.reference_area : float + Reference area of the rocket, in m². + FreeFormFin.Af : float + Area of the longitudinal section of each fin in the set. + FreeFormFin.AR : float + Aspect ratio of the fin. + FreeFormFin.gamma_c : float + Fin mid-chord sweep angle. + FreeFormFin.Yma : float + Span wise position of the mean aerodynamic chord. + FreeFormFin.roll_geometrical_constant : float + Geometrical constant used in roll calculations. + FreeFormFin.tau : float + Geometrical relation used to simplify lift and roll calculations. + FreeFormFin.lift_interference_factor : float + Factor of Fin-Body interference in the lift coefficient. + FreeFormFin.cp : tuple + Tuple with the x, y and z local coordinates of the fin set center of + pressure. Has units of length and is given in meters. + FreeFormFin.cpx : float + Fin set local center of pressure x coordinate. Has units of length and + is given in meters. + FreeFormFin.cpy : float + Fin set local center of pressure y coordinate. Has units of length and + is given in meters. + FreeFormFin.cpz : float + Fin set local center of pressure z coordinate. Has units of length and + is given in meters. + FreeFormFin.cl : Function + Function which defines the lift coefficient as a function of the angle + of attack and the Mach number. Takes as input the angle of attack in + radians and the Mach number. Returns the lift coefficient. + FreeFormFin.clalpha : float + Lift coefficient slope. Has units of 1/rad. + FreeFormFin.mac_length : float + Mean aerodynamic chord length of the fin set. + FreeFormFin.mac_lead : float + Mean aerodynamic chord leading edge x coordinate. + """ + + def __init__( + self, + angular_position, + shape_points, + rocket_radius, + cant_angle=0, + airfoil=None, + name="Free Form Fin", + ): + """Initialize FreeFormFin class. + + Parameters + ---------- + angular_position : float + Angular position of the fin in degrees measured as the rotation + around the symmetry axis of the rocket relative to one of the other + principal axis. See :ref:`Angular Position Inputs ` + shape_points : list + List of tuples (x, y) containing the coordinates of the fin's + geometry defining points. The point (0, 0) is the root leading edge. + Positive x is rearwards, positive y is upwards (span direction). + The shape will be interpolated between the points, in the order + they are given. The last point connects to the first point, and + represents the trailing edge. + rocket_radius : int, float + Reference radius to calculate lift coefficient, in meters. + cant_angle : int, float, optional + Fins cant angle with respect to the rocket centerline. Must + be given in degrees. + airfoil : tuple, optional + Default is null, in which case fins will be treated as flat plates. + Otherwise, if tuple, fins will be considered as airfoils. The + tuple's first item specifies the airfoil's lift coefficient + by angle of attack and must be either a .csv, .txt, ndarray + or callable. The .csv and .txt files can contain a single line + header and the first column must specify the angle of attack, while + the second column must specify the lift coefficient. The + ndarray should be as [(x0, y0), (x1, y1), (x2, y2), ...] + where x0 is the angle of attack and y0 is the lift coefficient. + If callable, it should take an angle of attack as input and + return the lift coefficient at that angle of attack. + The tuple's second item is the unit of the angle of attack, + accepting either "radians" or "degrees". + name : str + Name of the free form fin. + + Returns + ------- + None + """ + root_chord, span = _FreeFormGeometry.infer_dimensions(shape_points) + + super().__init__( + angular_position, + root_chord, + span, + rocket_radius, + cant_angle, + airfoil, + name, + ) + + self.geometry = _FreeFormGeometry(self, shape_points) + self._update_geometry_chain() + self.evaluate_shape() + + self.prints = _FreeFormFinPrints(self) + self.plots = _FreeFormFinPlots(self) + + def evaluate_center_of_pressure(self): + """Calculates and returns the center of pressure of the fin in local + coordinates. The center of pressure position is saved and stored as a + tuple. + + Returns + ------- + None + """ + # Center of pressure position in local coordinates + cpz = self.mac_lead + 0.25 * self.mac_length + self.cpx = 0 + self.cpy = self.Yma + self.cpz = cpz + self.cp = (self.cpx, self.cpy, self.cpz) + + @property + def shape_points(self): + return self.geometry.shape_points + + def to_dict(self, include_outputs=False): + data = super().to_dict(include_outputs=include_outputs) + data.update(self.geometry.get_data(include_outputs=include_outputs)) + return data + + @classmethod + def from_dict(cls, data): + return cls( + angular_position=data["angular_position"], + shape_points=data["shape_points"], + rocket_radius=data["rocket_radius"], + cant_angle=data["cant_angle"], + airfoil=data["airfoil"], + name=data["name"], + ) diff --git a/rocketpy/rocket/aero_surface/fins/free_form_fins.py b/rocketpy/rocket/aero_surface/fins/free_form_fins.py index 72758171e..d7c7e9512 100644 --- a/rocketpy/rocket/aero_surface/fins/free_form_fins.py +++ b/rocketpy/rocket/aero_surface/fins/free_form_fins.py @@ -1,9 +1,6 @@ -import warnings - -import numpy as np - from rocketpy.plots.aero_surface_plots import _FreeFormFinsPlots from rocketpy.prints.aero_surface_prints import _FreeFormFinsPrints +from rocketpy.rocket.aero_surface.fins._geometry import _FreeFormGeometry from .fins import Fins @@ -45,9 +42,9 @@ class FreeFormFins(Fins): Fin span in meters. FreeFormFins.name : string Name of fin set. - FreeFormFins.d : float + FreeFormFins.rocket_diameter : float Reference diameter of the rocket, in meters. - FreeFormFins.ref_area : float + FreeFormFins.reference_area : float Reference area of the rocket, in m². FreeFormFins.Af : float Area of the longitudinal section of each fin in the set. @@ -135,23 +132,7 @@ def __init__( ------- None """ - self.shape_points = shape_points - - down = False - for i in range(1, len(shape_points)): - if shape_points[i][1] > shape_points[i - 1][1] and down: - warnings.warn( - "Jagged fin shape detected. This may cause small inaccuracies " - "center of pressure and pitch moment calculations." - ) - break - if shape_points[i][1] < shape_points[i - 1][1]: - down = True - i += 1 - - root_chord = abs(shape_points[0][0] - shape_points[-1][0]) - ys = [point[1] for point in shape_points] - span = max(ys) - min(ys) + root_chord, span = _FreeFormGeometry.infer_dimensions(shape_points) super().__init__( n, @@ -163,10 +144,9 @@ def __init__( name, ) - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() + self.geometry = _FreeFormGeometry(self, shape_points) + self._update_geometry_chain() + self.evaluate_shape() self.prints = _FreeFormFinsPrints(self) self.plots = _FreeFormFinsPlots(self) @@ -187,215 +167,24 @@ def evaluate_center_of_pressure(self): self.cpz = cpz self.cp = (self.cpx, self.cpy, self.cpz) - def evaluate_geometrical_parameters(self): # pylint: disable=too-many-statements - """ - Calculates and saves the fin set's geometrical parameters such as the - fin area, aspect ratio, and parameters related to roll movement. This - method uses the same calculations to those in OpenRocket for free-form - fin shapes. - - Returns - ------- - None - """ - # pylint: disable=invalid-name - # pylint: disable=too-many-locals - # Calculate the fin area (Af) using the Shoelace theorem (polygon area formula) - Af = 0 - for i in range(len(self.shape_points) - 1): - x1, y1 = self.shape_points[i] - x2, y2 = self.shape_points[i + 1] - Af += (y1 + y2) * (x1 - x2) - Af = abs(Af) / 2 - if Af < 1e-6: - raise ValueError("Fin area is too small. Check the shape_points.") - - # Calculate aspect ratio (AR) and lift interference factors - AR = 2 * self.span**2 / Af # Aspect ratio - tau = (self.span + self.rocket_radius) / self.rocket_radius - lift_interference_factor = 1 + 1 / tau - - # Calculate roll forcing interference factor using OpenRocket's approach - roll_forcing_interference_factor = (1 / np.pi**2) * ( - (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) - + ((np.pi * (tau**2 + 1) ** 2) / (tau**2 * (tau - 1) ** 2)) - * np.arcsin((tau**2 - 1) / (tau**2 + 1)) - - (2 * np.pi * (tau + 1)) / (tau * (tau - 1)) - + ((tau**2 + 1) ** 2 / (tau**2 * (tau - 1) ** 2)) - * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2 - - (4 * (tau + 1)) - / (tau * (tau - 1)) - * np.arcsin((tau**2 - 1) / (tau**2 + 1)) - + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau)) - ) - - # Define number of interpolation points along the span of the fin - points_per_line = 40 # Same as OpenRocket - - # Initialize arrays for leading/trailing edge and chord lengths - chord_lead = np.ones(points_per_line) * np.inf # Leading edge x coordinates - chord_trail = np.ones(points_per_line) * -np.inf # Trailing edge x coordinates - chord_length = np.zeros( - points_per_line - ) # Chord length for each spanwise section - - # Discretize fin shape and calculate chord length, leading, and trailing edges - for p in range(1, len(self.shape_points)): - x1, y1 = self.shape_points[p - 1] - x2, y2 = self.shape_points[p] - - # Compute corresponding points along the fin span (clamp to valid range) - prev_idx = int(y1 / self.span * (points_per_line - 1)) - curr_idx = int(y2 / self.span * (points_per_line - 1)) - prev_idx = np.clip(prev_idx, 0, points_per_line - 1) - curr_idx = np.clip(curr_idx, 0, points_per_line - 1) - - if prev_idx > curr_idx: - prev_idx, curr_idx = curr_idx, prev_idx - - # Compute intersection of fin edge with each spanwise section - for i in range(prev_idx, curr_idx + 1): - y = i * self.span / (points_per_line - 1) - if y1 != y2: - x = np.clip( - (y - y2) / (y1 - y2) * x1 + (y1 - y) / (y1 - y2) * x2, - min(x1, x2), - max(x1, x2), - ) - else: - x = x1 # Handle horizontal segments - - # Update leading and trailing edge positions - chord_lead[i] = min(chord_lead[i], x) - chord_trail[i] = max(chord_trail[i], x) - - # Update chord length - if y1 < y2: - chord_length[i] -= x - else: - chord_length[i] += x - - # Replace infinities and handle invalid values in chord_lead and chord_trail - for i in range(points_per_line): - if ( - np.isinf(chord_lead[i]) - or np.isinf(chord_trail[i]) - or np.isnan(chord_lead[i]) - or np.isnan(chord_trail[i]) - ): - chord_lead[i] = 0 - chord_trail[i] = 0 - if chord_length[i] < 0 or np.isnan(chord_length[i]): - chord_length[i] = 0 - if chord_length[i] > chord_trail[i] - chord_lead[i]: - chord_length[i] = chord_trail[i] - chord_lead[i] - - # Initialize integration variables for various aerodynamic and roll properties - radius = self.rocket_radius - total_area = 0 - mac_length = 0 # Mean aerodynamic chord length - mac_lead = 0 # Mean aerodynamic chord leading edge - mac_span = 0 # Mean aerodynamic chord spanwise position (Yma) - cos_gamma_sum = 0 # Sum of cosine of the sweep angle - roll_geometrical_constant = 0 - roll_damping_numerator = 0 - roll_damping_denominator = 0 - - # Perform integration over spanwise sections - dy = self.span / (points_per_line - 1) - for i in range(points_per_line): - chord = chord_trail[i] - chord_lead[i] - y = i * dy - - # Update integration variables - mac_length += chord * chord - mac_span += y * chord - mac_lead += chord_lead[i] * chord - total_area += chord - roll_geometrical_constant += chord_length[i] * (radius + y) ** 2 - roll_damping_numerator += radius**3 * chord / (radius + y) ** 2 - roll_damping_denominator += (radius + y) * chord - - # Update cosine of sweep angle (cos_gamma) - if i > 0: - dx = (chord_trail[i] + chord_lead[i]) / 2 - ( - chord_trail[i - 1] + chord_lead[i - 1] - ) / 2 - cos_gamma_sum += dy / np.hypot(dx, dy) - - # Finalize mean aerodynamic chord properties - mac_length *= dy - mac_span *= dy - mac_lead *= dy - total_area *= dy - roll_geometrical_constant *= dy - roll_damping_numerator *= dy - roll_damping_denominator *= dy - - mac_length /= total_area - mac_span /= total_area - mac_lead /= total_area - cos_gamma = cos_gamma_sum / (points_per_line - 1) - - # Store computed values - self.Af = Af # Fin area - self.AR = AR # Aspect ratio - self.gamma_c = np.arccos(cos_gamma) # Sweep angle - self.Yma = mac_span # Mean aerodynamic chord spanwise position - self.mac_length = mac_length - self.mac_lead = mac_lead - self.tau = tau - self.roll_geometrical_constant = roll_geometrical_constant - self.lift_interference_factor = lift_interference_factor - self.roll_forcing_interference_factor = roll_forcing_interference_factor - self.roll_damping_interference_factor = 1 + ( - roll_damping_numerator / roll_damping_denominator - ) - - # Evaluate the shape and finalize geometry - self.evaluate_shape() - - def evaluate_shape(self): - x_array, y_array = zip(*self.shape_points) - self.shape_vec = [np.array(x_array), np.array(y_array)] + @property + def shape_points(self): + return self.geometry.shape_points def to_dict(self, **kwargs): data = super().to_dict(**kwargs) - data["shape_points"] = self.shape_points - - if kwargs.get("include_outputs", False): - data.update( - { - "Af": self.Af, - "AR": self.AR, - "gamma_c": self.gamma_c, - "Yma": self.Yma, - "mac_length": self.mac_length, - "mac_lead": self.mac_lead, - "roll_geometrical_constant": self.roll_geometrical_constant, - "tau": self.tau, - "lift_interference_factor": self.lift_interference_factor, - "roll_forcing_interference_factor": self.roll_forcing_interference_factor, - "roll_damping_interference_factor": self.roll_damping_interference_factor, - } - ) + data.update( + self.geometry.get_data(include_outputs=kwargs.get("include_outputs", False)) + ) return data @classmethod def from_dict(cls, data): return cls( - data["n"], - data["shape_points"], - data["rocket_radius"], - data["cant_angle"], - data["airfoil"], - data["name"], + n=data["n"], + shape_points=data["shape_points"], + rocket_radius=data["rocket_radius"], + cant_angle=data["cant_angle"], + airfoil=data["airfoil"], + name=data["name"], ) - - def info(self): - self.prints.geometry() - self.prints.lift() - - def all_info(self): - self.prints.all() - self.plots.all() diff --git a/rocketpy/rocket/aero_surface/fins/trapezoidal_fin.py b/rocketpy/rocket/aero_surface/fins/trapezoidal_fin.py new file mode 100644 index 000000000..c58055945 --- /dev/null +++ b/rocketpy/rocket/aero_surface/fins/trapezoidal_fin.py @@ -0,0 +1,240 @@ +from rocketpy.plots.aero_surface_plots import _TrapezoidalFinPlots +from rocketpy.prints.aero_surface_prints import _TrapezoidalFinPrints +from rocketpy.rocket.aero_surface.fins._geometry import _TrapezoidalGeometry + +from .fin import Fin + + +class TrapezoidalFin(Fin): + """A class used to represent a single trapezoidal fin. + + This class inherits from the Fin class. + + Note + ---- + Local coordinate system: + - Origin located at the top of the root chord. + - Z axis along the longitudinal axis of symmetry, positive downwards (top -> bottom). + - Y axis perpendicular to the Z axis, in the span direction, positive upwards. + - X axis completes the right-handed coordinate system. + + See Also + -------- + Fin : Parent class + + Attributes + ---------- + TrapezoidalFin.angular_position : float + Angular position of the fin set with respect to the rocket centerline, + in degrees. + TrapezoidalFin.rocket_radius : float + The reference rocket radius used for lift coefficient normalization, in + meters. + TrapezoidalFin.airfoil : tuple + Tuple of two items. First is the airfoil lift curve. + Second is the unit of the curve (radians or degrees). + TrapezoidalFin.cant_angle : float + Fins cant angle with respect to the rocket centerline, in degrees. + TrapezoidalFin.cant_angle_rad : float + Fins cant angle with respect to the rocket centerline, in radians. + TrapezoidalFin.root_chord : float + Fin root chord in meters. + TrapezoidalFin.tip_chord : float + Fin tip chord in meters. + TrapezoidalFin.span : float + Fin span in meters. + TrapezoidalFin.name : string + Name of fin set. + TrapezoidalFin.sweep_length : float + Fins sweep length in meters. By sweep length, understand the axial + distance between the fin root leading edge and the fin tip leading edge + measured parallel to the rocket centerline. + TrapezoidalFin.sweep_angle : float + Fins sweep angle with respect to the rocket centerline. Must + be given in degrees. + TrapezoidalFin.rocket_diameter : float + Reference diameter of the rocket, in meters. + TrapezoidalFins.fin_area : float + Area of the longitudinal section of each fin in the set. + TrapezoidalFins.AR : float + Aspect ratio of the fin. + TrapezoidalFin.gamma_c : float + Fin mid-chord sweep angle. + TrapezoidalFin.yma : float + Span wise position of the mean aerodynamic chord. + TrapezoidalFin.roll_geometrical_constant : float + Geometrical constant used in roll calculations. + TrapezoidalFin.tau : float + Geometrical relation used to simplify lift and roll calculations. + TrapezoidalFin.lift_interference_factor : float + Factor of Fin-Body interference in the lift coefficient. + TrapezoidalFin.cp : tuple + Tuple with the x, y and z local coordinates of the fin set center of + pressure. Has units of length and is given in meters. + TrapezoidalFin.cpx : float + Fin set local center of pressure x coordinate. Has units of length and + is given in meters. + TrapezoidalFin.cpy : float + Fin set local center of pressure y coordinate. Has units of length and + is given in meters. + TrapezoidalFin.cpz : float + Fin set local center of pressure z coordinate. Has units of length and + is given in meters. + """ + + def __init__( + self, + angular_position, + root_chord, + tip_chord, + span, + rocket_radius, + cant_angle=0, + sweep_length=None, + sweep_angle=None, + airfoil=None, + name="Trapezoidal Fin", + ): + """Initializes the TrapezoidalFin class. + + Parameters + ---------- + angular_position : float + Angular position of the fin in degrees measured as the rotation + around the symmetry axis of the rocket relative to one of the other + principal axis. See :ref:`Angular Position Inputs ` + root_chord : int, float + Fin root chord in meters. + tip_chord : int, float + Fin tip chord in meters. + span : int, float + Fin span in meters. + rocket_radius : int, float + Reference radius to calculate lift coefficient, in meters. + cant_angle : int, float, optional + Fins cant angle with respect to the rocket centerline. Must + be given in degrees. + sweep_length : int, float, optional + Fins sweep length in meters. By sweep length, understand the axial + distance between the fin root leading edge and the fin tip leading + edge measured parallel to the rocket centerline. If not given, the + sweep length is assumed to be equal the root chord minus the tip + chord, in which case the fin is a right trapezoid with its base + perpendicular to the rocket's axis. Cannot be used in conjunction + with sweep_angle. + sweep_angle : int, float, optional + Fins sweep angle with respect to the rocket centerline. Must + be given in degrees. If not given, the sweep angle is automatically + calculated, in which case the fin is assumed to be a right trapezoid + with its base perpendicular to the rocket's axis. + Cannot be used in conjunction with sweep_length. + airfoil : tuple, optional + Default is null, in which case fins will be treated as flat plates. + Otherwise, if tuple, fins will be considered as airfoils. The + tuple's first item specifies the airfoil's lift coefficient + by angle of attack and must be either a .csv, .txt, ndarray + or callable. The .csv and .txt files can contain a single line + header and the first column must specify the angle of attack, while + the second column must specify the lift coefficient. The + ndarray should be as [(x0, y0), (x1, y1), (x2, y2), ...] + where x0 is the angle of attack and y0 is the lift coefficient. + If callable, it should take an angle of attack as input and + return the lift coefficient at that angle of attack. + The tuple's second item is the unit of the angle of attack, + accepting either "radians" or "degrees". + name : str + Name of the trapezoidal fin. + """ + super().__init__( + angular_position, + root_chord, + span, + rocket_radius, + cant_angle, + airfoil, + name, + ) + + self.geometry = _TrapezoidalGeometry( + self, + tip_chord=tip_chord, + sweep_length=sweep_length, + sweep_angle=sweep_angle, + ) + self._update_geometry_chain() + self.evaluate_shape() + self.evaluate_rotation_matrix() + + self.prints = _TrapezoidalFinPrints(self) + self.plots = _TrapezoidalFinPlots(self) + + @property + def tip_chord(self): + return self.geometry.tip_chord + + @tip_chord.setter + def tip_chord(self, value): + self.geometry.tip_chord = value + self._update_geometry_chain() + self.evaluate_shape() + + @property + def sweep_angle(self): + return self.geometry.sweep_angle + + @sweep_angle.setter + def sweep_angle(self, value): + self.geometry.sweep_angle = value + self._update_geometry_chain() + self.evaluate_shape() + + @property + def sweep_length(self): + return self.geometry.sweep_length + + @sweep_length.setter + def sweep_length(self, value): + self.geometry.sweep_length = value + self._update_geometry_chain() + self.evaluate_shape() + + def evaluate_center_of_pressure(self): + """Calculates and returns the center of pressure of the fin in local + coordinates. The center of pressure position is saved and stored as a + tuple. + + Returns + ------- + None + """ + # Center of pressure position in local coordinates + cpz = (self.sweep_length / 3) * ( + (self.root_chord + 2 * self.tip_chord) / (self.root_chord + self.tip_chord) + ) + (1 / 6) * ( + self.root_chord + + self.tip_chord + - self.root_chord * self.tip_chord / (self.root_chord + self.tip_chord) + ) + self.cpx = 0 + self.cpy = self.Yma + self.cpz = cpz + self.cp = (self.cpx, self.cpy, self.cpz) + + def to_dict(self, include_outputs=False): + data = super().to_dict(include_outputs=include_outputs) + data.update(self.geometry.get_data(include_outputs=include_outputs)) + return data + + @classmethod + def from_dict(cls, data): + return cls( + angular_position=data["angular_position"], + root_chord=data["root_chord"], + tip_chord=data["tip_chord"], + span=data["span"], + rocket_radius=data["rocket_radius"], + cant_angle=data["cant_angle"], + sweep_length=data.get("sweep_length"), + airfoil=data["airfoil"], + name=data["name"], + ) diff --git a/rocketpy/rocket/aero_surface/fins/trapezoidal_fins.py b/rocketpy/rocket/aero_surface/fins/trapezoidal_fins.py index c6b4ea633..2c9adea58 100644 --- a/rocketpy/rocket/aero_surface/fins/trapezoidal_fins.py +++ b/rocketpy/rocket/aero_surface/fins/trapezoidal_fins.py @@ -1,7 +1,6 @@ -import numpy as np - from rocketpy.plots.aero_surface_plots import _TrapezoidalFinsPlots from rocketpy.prints.aero_surface_prints import _TrapezoidalFinsPrints +from rocketpy.rocket.aero_surface.fins._geometry import _TrapezoidalGeometry from .fins import Fins @@ -35,9 +34,6 @@ class TrapezoidalFins(Fins): Second is the unit of the curve (radians or degrees). TrapezoidalFins.cant_angle : float Fins cant angle with respect to the rocket centerline, in degrees. - TrapezoidalFins.changing_attribute_dict : dict - Dictionary that stores the name and the values of the attributes that - may be changed during a simulation. Useful for control systems. TrapezoidalFins.cant_angle_rad : float Fins cant angle with respect to the rocket centerline, in radians. TrapezoidalFins.root_chord : float @@ -55,9 +51,9 @@ class TrapezoidalFins(Fins): TrapezoidalFins.sweep_angle : float Fins sweep angle with respect to the rocket centerline. Must be given in degrees. - TrapezoidalFins.d : float + TrapezoidalFins.rocket_diameter : float Reference diameter of the rocket, in meters. - TrapezoidalFins.ref_area : float + TrapezoidalFins.reference_area : float Reference area of the rocket, in m². TrapezoidalFins.Af : float Area of the longitudinal section of each fin in the set. @@ -169,62 +165,47 @@ def __init__( name, ) - # Check if sweep angle or sweep length is given - if sweep_length is not None and sweep_angle is not None: - raise ValueError("Cannot use sweep_length and sweep_angle together") - elif sweep_angle is not None: - sweep_length = np.tan(sweep_angle * np.pi / 180) * span - elif sweep_length is None: - sweep_length = root_chord - tip_chord - - self._tip_chord = tip_chord - self._sweep_length = sweep_length - self._sweep_angle = sweep_angle - - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() + self.geometry = _TrapezoidalGeometry( + self, + tip_chord=tip_chord, + sweep_length=sweep_length, + sweep_angle=sweep_angle, + ) + self._update_geometry_chain() + self.evaluate_shape() self.prints = _TrapezoidalFinsPrints(self) self.plots = _TrapezoidalFinsPlots(self) @property def tip_chord(self): - return self._tip_chord + return self.geometry.tip_chord @tip_chord.setter def tip_chord(self, value): - self._tip_chord = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() + self.geometry.tip_chord = value + self._update_geometry_chain() + self.evaluate_shape() @property def sweep_angle(self): - return self._sweep_angle + return self.geometry.sweep_angle @sweep_angle.setter def sweep_angle(self, value): - self._sweep_angle = value - self._sweep_length = np.tan(value * np.pi / 180) * self.span - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() + self.geometry.sweep_angle = value + self._update_geometry_chain() + self.evaluate_shape() @property def sweep_length(self): - return self._sweep_length + return self.geometry.sweep_length @sweep_length.setter def sweep_length(self, value): - self._sweep_length = value - self.evaluate_geometrical_parameters() - self.evaluate_center_of_pressure() - self.evaluate_lift_coefficient() - self.evaluate_roll_parameters() + self.geometry.sweep_length = value + self._update_geometry_chain() + self.evaluate_shape() def evaluate_center_of_pressure(self): """Calculates and returns the center of pressure of the fin set in local @@ -248,124 +229,11 @@ def evaluate_center_of_pressure(self): self.cpz = cpz self.cp = (self.cpx, self.cpy, self.cpz) - def evaluate_geometrical_parameters(self): # pylint: disable=too-many-statements - """Calculates and saves fin set's geometrical parameters such as the - fins' area, aspect ratio and parameters for roll movement. - - Returns - ------- - None - """ - # pylint: disable=invalid-name - Yr = self.root_chord + self.tip_chord - Af = Yr * self.span / 2 # Fin area - AR = 2 * self.span**2 / Af # Fin aspect ratio - gamma_c = np.arctan( - (self.sweep_length + 0.5 * self.tip_chord - 0.5 * self.root_chord) - / (self.span) - ) - Yma = ( - (self.span / 3) * (self.root_chord + 2 * self.tip_chord) / Yr - ) # Span wise coord of mean aero chord - - # Fin–body interference correction parameters - tau = (self.span + self.rocket_radius) / self.rocket_radius - lift_interference_factor = 1 + 1 / tau - lambda_ = self.tip_chord / self.root_chord - - # Parameters for Roll Moment. - # Documented at: https://docs.rocketpy.org/en/latest/technical/ - roll_geometrical_constant = ( - (self.root_chord + 3 * self.tip_chord) * self.span**3 - + 4 - * (self.root_chord + 2 * self.tip_chord) - * self.rocket_radius - * self.span**2 - + 6 * (self.root_chord + self.tip_chord) * self.span * self.rocket_radius**2 - ) / 12 - roll_damping_interference_factor = 1 + ( - ((tau - lambda_) / (tau)) - ((1 - lambda_) / (tau - 1)) * np.log(tau) - ) / ( - ((tau + 1) * (tau - lambda_)) / (2) - - ((1 - lambda_) * (tau**3 - 1)) / (3 * (tau - 1)) - ) - roll_forcing_interference_factor = (1 / np.pi**2) * ( - (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) - + ((np.pi * (tau**2 + 1) ** 2) / (tau**2 * (tau - 1) ** 2)) - * np.arcsin((tau**2 - 1) / (tau**2 + 1)) - - (2 * np.pi * (tau + 1)) / (tau * (tau - 1)) - + ((tau**2 + 1) ** 2) - / (tau**2 * (tau - 1) ** 2) - * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2 - - (4 * (tau + 1)) - / (tau * (tau - 1)) - * np.arcsin((tau**2 - 1) / (tau**2 + 1)) - + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau)) - ) - - # Store values - self.Yr = Yr - self.Af = Af # Fin area - self.AR = AR # Aspect Ratio - self.gamma_c = gamma_c # Mid chord angle - self.Yma = Yma # Span wise coord of mean aero chord - self.roll_geometrical_constant = roll_geometrical_constant - self.tau = tau - self.lift_interference_factor = lift_interference_factor - self.λ = lambda_ # pylint: disable=non-ascii-name - self.roll_damping_interference_factor = roll_damping_interference_factor - self.roll_forcing_interference_factor = roll_forcing_interference_factor - - self.evaluate_shape() - - def evaluate_shape(self): - if self.sweep_length: - points = [ - (0, 0), - (self.sweep_length, self.span), - (self.sweep_length + self.tip_chord, self.span), - (self.root_chord, 0), - ] - else: - points = [ - (0, 0), - (self.root_chord - self.tip_chord, self.span), - (self.root_chord, self.span), - (self.root_chord, 0), - ] - - x_array, y_array = zip(*points) - self.shape_vec = [np.array(x_array), np.array(y_array)] - - def info(self): - self.prints.geometry() - self.prints.lift() - - def all_info(self): - self.prints.all() - self.plots.all() - 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 kwargs.get("include_outputs", False): - data.update( - { - "shape_vec": self.shape_vec, - "Af": self.Af, - "AR": self.AR, - "gamma_c": self.gamma_c, - "Yma": self.Yma, - "roll_geometrical_constant": self.roll_geometrical_constant, - "tau": self.tau, - "lift_interference_factor": self.lift_interference_factor, - "roll_damping_interference_factor": self.roll_damping_interference_factor, - "roll_forcing_interference_factor": self.roll_forcing_interference_factor, - } - ) + data.update( + self.geometry.get_data(include_outputs=kwargs.get("include_outputs", False)) + ) return data @classmethod diff --git a/rocketpy/rocket/aero_surface/generic_surface.py b/rocketpy/rocket/aero_surface/generic_surface.py index 8ab438620..d1ae3a8f0 100644 --- a/rocketpy/rocket/aero_surface/generic_surface.py +++ b/rocketpy/rocket/aero_surface/generic_surface.py @@ -85,6 +85,8 @@ def __init__( self.cpz = center_of_pressure[2] self.name = name + self._rotation_surface_to_body = Matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) + default_coefficients = self._get_default_coefficients() self._check_coefficients(coefficients, default_coefficients) coefficients = self._complete_coefficients(coefficients, default_coefficients) diff --git a/rocketpy/rocket/aero_surface/rail_buttons.py b/rocketpy/rocket/aero_surface/rail_buttons.py index 89331c99f..7d3a9bd30 100644 --- a/rocketpy/rocket/aero_surface/rail_buttons.py +++ b/rocketpy/rocket/aero_surface/rail_buttons.py @@ -17,6 +17,7 @@ class RailButtons(AeroSurface): Angular position of the rail buttons in degrees measured as the rotation around the symmetry axis of the rocket relative to one of the other principal axis. + See :ref:`Angular Position Inputs ` RailButtons.angular_position_rad : float Angular position of the rail buttons in radians. RailButtons.button_height : float, optional diff --git a/rocketpy/rocket/components.py b/rocketpy/rocket/components.py index 9ce8595b3..57e4d12f8 100644 --- a/rocketpy/rocket/components.py +++ b/rocketpy/rocket/components.py @@ -1,4 +1,5 @@ from collections import namedtuple +from copy import deepcopy class Components: @@ -186,7 +187,7 @@ def clear(self): self._components.clear() def sort_by_position(self, reverse=False): - """Sort the list of components by z axis position. + """Returns a new Components object sorted components by z axis position. Parameters ---------- @@ -196,9 +197,12 @@ def sort_by_position(self, reverse=False): Returns ------- - None + Components + A new Components object sorted by component position. """ - self._components.sort(key=lambda x: x.position.z, reverse=reverse) + components = deepcopy(self) + components._components.sort(key=lambda x: x.position.z, reverse=reverse) + return components def to_dict(self, **kwargs): # pylint: disable=unused-argument return { diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index e3692d2e8..84655b6d7 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -21,7 +21,10 @@ Tail, TrapezoidalFins, ) +from rocketpy.rocket.aero_surface.fins.elliptical_fin import EllipticalFin +from rocketpy.rocket.aero_surface.fins.free_form_fin import FreeFormFin from rocketpy.rocket.aero_surface.fins.free_form_fins import FreeFormFins +from rocketpy.rocket.aero_surface.fins.trapezoidal_fin import TrapezoidalFin from rocketpy.rocket.aero_surface.generic_surface import GenericSurface from rocketpy.rocket.components import Components from rocketpy.rocket.parachute import Parachute @@ -661,14 +664,20 @@ def __evaluate_single_surface_cp_to_cdm(self, surface, position): """Calculates the relative position of each aerodynamic surface center of pressure to the rocket's center of dry mass in Body Axes Coordinate System.""" - pos = Vector( + # position of the surfaces coordinate system origin in body frame + pos_origin = Vector( [ - (position.x - self.cm_eccentricity_x) * self._csys - surface.cpx, - (position.y - self.cm_eccentricity_y) - surface.cpy, - (position.z - self.center_of_dry_mass_position) * self._csys - - surface.cpz, + (position.x - self.cm_eccentricity_x) * self._csys, + (position.y - self.cm_eccentricity_y), + (position.z - self.center_of_dry_mass_position) * self._csys, ] ) + # position of the center of pressure in body frame + pos = ( + surface._rotation_surface_to_body + @ Vector([surface.cpx, surface.cpy, surface.cpz]) + + pos_origin + ) # TODO: this should be recomputed whenever cant angle changes for fin self.surfaces_cp_to_cdm[surface] = pos def evaluate_stability_margin(self): @@ -1066,11 +1075,20 @@ def __add_single_surface(self, surface, position): """Adds a single aerodynamic surface to the rocket. Makes checks for rail buttons case, and position type. """ - position = ( - Vector([0, 0, position]) - if not isinstance(position, (Vector, tuple, list)) - else Vector(position) - ) + if isinstance(surface, (TrapezoidalFin, EllipticalFin, FreeFormFin)): + # TODO: the leading edge position should be recomputed whenever cant + # angle of the fin changes, but currently it is only computed at the + # moment the fin is added to the rocket. Detecting when the cant + # angle changes is hard, because it is a parameter of the fin, while + # the leading edge position is only defined on the rocket + position = surface._compute_leading_edge_position(position, self._csys) + else: + position = ( + Vector([0, 0, position]) + if not isinstance(position, (Vector, tuple, list)) + else Vector(position) + ) + if isinstance(surface, RailButtons): self.rail_buttons = Components() self.rail_buttons.add(surface, position) @@ -1085,17 +1103,23 @@ def add_surfaces(self, surfaces, positions): Parameters ---------- - surfaces : list, AeroSurface, NoseCone, TrapezoidalFins, EllipticalFins, Tail, RailButtons + surfaces : list[AeroSurface], AeroSurface Aerodynamic surface to be added to the rocket. Can be a list of AeroSurface if more than one surface is to be added. - positions : int, float, list, tuple, Vector - Position, in m, of the aerodynamic surface's center of pressure - relative to the user defined rocket coordinate system. - If a list is passed, it will correspond to the position of each item - in the surfaces list. - For NoseCone type, position is relative to the nose cone tip. - For Fins type, position is relative to the point belonging to - the root chord which is highest in the rocket coordinate system. + positions : int, float, tuple, list, Vector + Position(s) of the aerodynamic surface's reference point. Can be: + + - a single number (int or float) giving the z-coordinate along + the rocket axis. + - a sequence of three numbers (x, y, z) representing the full + position in the user-defined coordinate system. + + If passing multiple surfaces, provide a list of positions matching + each surface in order. + For NoseCone type, position is the tip coordinate along the axis. + For Fins type, position refers to the z-coordinate of the root + chord leading-edge point closest to the nose cone, before any + cant-angle offset is considered. For Tail type, position is relative to the point belonging to the tail which is highest in the rocket coordinate system. For RailButtons type, position is relative to the lower rail button. @@ -1108,10 +1132,18 @@ def add_surfaces(self, surfaces, positions): ------- None """ - try: + if isinstance(surfaces, Iterable): + if isinstance(positions, Iterable): + if len(surfaces) != len(positions): + raise ValueError( + "The number of surfaces and positions must be the same." + ) + else: + positions = [positions] * len(surfaces) + for surface, position in zip(surfaces, positions): self.__add_single_surface(surface, position) - except TypeError: + else: self.__add_single_surface(surfaces, positions) self.evaluate_center_of_pressure() @@ -1285,10 +1317,10 @@ def add_trapezoidal_fins( tip_chord : int, float Fin tip chord in meters. position : int, float - Fin set position relative to the rocket's coordinate system. - By fin set position, understand the point belonging to the root - chord which is highest in the rocket coordinate system (i.e. - the point closest to the nose cone tip). + Fin set position in the z coordinate of the user defined rocket + coordinate system. By fin set position, understand the point + belonging to the root chord which is highest in the rocket + coordinate system (i.e. the point closest to the nose cone tip). See Also -------- @@ -1334,6 +1366,12 @@ def add_trapezoidal_fins( fin_set : TrapezoidalFins Fin set object created. """ + if n <= 2: + raise ValueError( + "Number of fins must be greater than 2. " + "For 1 or 2 fins, create a FreeFormFin object " + "and add it to the rocket using the add_surfaces method." + ) # Modify radius if not given, use rocket radius, otherwise use given. radius = radius if radius is not None else self.radius @@ -1381,10 +1419,10 @@ def add_elliptical_fins( span : int, float Fin span in meters. position : int, float - Fin set position relative to the rocket's coordinate system. By fin - set position, understand the point belonging to the root chord which - is highest in the rocket coordinate system (i.e. the point - closest to the nose cone tip). + Fin set position in the z coordinate of the user defined rocket + coordinate system. By fin set position, understand the point + belonging to the root chord which is highest in the rocket + coordinate system (i.e. the point closest to the nose cone tip). See Also -------- @@ -1420,6 +1458,13 @@ def add_elliptical_fins( fin_set : EllipticalFins Fin set object created. """ + if n <= 2: + raise ValueError( + "Number of fins must be greater than 2. " + "For 1 or 2 fins, create a FreeFormFin object " + "and add it to the rocket using the add_surfaces method." + ) + radius = radius if radius is not None else self.radius fin_set = EllipticalFins(n, root_chord, span, radius, cant_angle, airfoil, name) self.add_surfaces(fin_set, position) @@ -1451,10 +1496,10 @@ def add_free_form_fins( The shape will be interpolated between the points, in the order they are given. The last point connects to the first point. position : int, float - Fin set position relative to the rocket's coordinate system. - By fin set position, understand the point belonging to the root - chord which is highest in the rocket coordinate system (i.e. - the point closest to the nose cone tip). + Fin set position in the z coordinate of the user defined rocket + coordinate system. By fin set position, understand the point + belonging to the root chord which is highest in the rocket + coordinate system (i.e. the point closest to the nose cone tip). See Also -------- @@ -1486,6 +1531,12 @@ def add_free_form_fins( fin_set : FreeFormFins Fin set object created. """ + if n <= 2: + raise ValueError( + "Number of fins must be greater than 2. " + "For 1 or 2 fins, create a FreeFormFin object " + "and add it to the rocket using the add_surfaces method." + ) # Modify radius if not given, use rocket radius, otherwise use given. radius = radius if radius is not None else self.radius @@ -1630,8 +1681,7 @@ def add_sensor(self, sensor, position): must be in the format (x, y, z) where x, y, and z are defined in the rocket's user defined coordinate system. If a single value is passed, it is assumed to be along the z-axis (centerline) of the - rocket's user defined coordinate system and angular_position and - radius must be given. + rocket's user defined coordinate system. Returns ------- @@ -1825,7 +1875,7 @@ def set_rail_buttons( as the rotation around the symmetry axis of the rocket relative to one of the other principal axis. Default value is 45 degrees, generally used in rockets with - 4 fins. + 4 fins. See :ref:`Angular Position Inputs ` radius : int, float, optional Fuselage radius where the rail buttons are located. diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 1443d1d80..2293d9706 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -698,15 +698,6 @@ def __simulate(self, verbose): self.__process_sensors_and_controllers_at_current_node(node, phase) - for controller in node._controllers: - controller( - self.t, - self.y_sol, - self.solution, - self.sensors, - self.env, - ) - for parachute in node.parachutes: # Calculate and save pressure signal ( @@ -853,7 +844,7 @@ def __process_sensors_and_controllers_at_current_node(self, node, phase): phase : FlightPhase The current flight phase. """ - if self.sensors: + if node._component_sensors: u_dot = phase.derivative(self.t, self.y_sol) self.__measure_sensors(node._component_sensors, u_dot) diff --git a/rocketpy/simulation/monte_carlo.py b/rocketpy/simulation/monte_carlo.py index e10789a7d..42a566b7b 100644 --- a/rocketpy/simulation/monte_carlo.py +++ b/rocketpy/simulation/monte_carlo.py @@ -525,6 +525,73 @@ def estimate_confidence_interval( return res.confidence_interval + def simulate_convergence( + self, + target_attribute="apogee_time", + target_confidence=0.95, + tolerance=0.5, + max_simulations=1000, + batch_size=50, + parallel=False, + n_workers=None, + ): + """Run Monte Carlo simulations in batches until the confidence interval + width converges within the specified tolerance or the maximum number of + simulations is reached. + + Parameters + ---------- + target_attribute : str + The target attribute to track its convergence (e.g., "apogee", "apogee_time", etc.). + target_confidence : float, optional + The confidence level for the interval (between 0 and 1). Default is 0.95. + tolerance : float, optional + The desired width of the confidence interval in seconds, meters, or other units. Default is 0.5. + max_simulations : int, optional + The maximum number of simulations to run to avoid infinite loops. Default is 1000. + batch_size : int, optional + The number of simulations to run in each batch. Default is 50. + parallel : bool, optional + Whether to run simulations in parallel. Default is False. + n_workers : int, optional + The number of worker processes to use if running in parallel. Default is None. + + Returns + ------- + confidence_interval_history : list of float + History of confidence interval widths, one value per batch of simulations. + The last element corresponds to the width when the simulation stopped for + either meeting the tolerance or reaching the maximum number of simulations. + """ + + self.import_outputs(self.filename.with_suffix(".outputs.txt")) + confidence_interval_history = [] + + while self.num_of_loaded_sims < max_simulations: + total_sims = min(self.num_of_loaded_sims + batch_size, max_simulations) + + self.simulate( + number_of_simulations=total_sims, + append=True, + include_function_data=False, + parallel=parallel, + n_workers=n_workers, + ) + + self.import_outputs(self.filename.with_suffix(".outputs.txt")) + + ci = self.estimate_confidence_interval( + attribute=target_attribute, + confidence_level=target_confidence, + ) + + confidence_interval_history.append(float(ci.high - ci.low)) + + if float(ci.high - ci.low) <= tolerance: + break + + return confidence_interval_history + def __evaluate_flight_inputs(self, sim_idx): """Evaluates the inputs of a single flight simulation. diff --git a/tests/acceptance/test_bella_lui_rocket.py b/tests/acceptance/test_bella_lui_rocket.py index bcfe325bc..cdccb67d8 100644 --- a/tests/acceptance/test_bella_lui_rocket.py +++ b/tests/acceptance/test_bella_lui_rocket.py @@ -67,6 +67,7 @@ def test_bella_lui_rocket_data_asserts_acceptance(): type="Reanalysis", file="data/weather/bella_lui_weather_data_ERA5.nc", dictionary="ECMWF", + pressure_conversion_factor="hPa", ) env.max_expected_height = 2000 diff --git a/tests/acceptance/test_ndrt_2020_rocket.py b/tests/acceptance/test_ndrt_2020_rocket.py index 7874ee164..04dae8725 100644 --- a/tests/acceptance/test_ndrt_2020_rocket.py +++ b/tests/acceptance/test_ndrt_2020_rocket.py @@ -83,6 +83,7 @@ def test_ndrt_2020_rocket_data_asserts_acceptance(env_file): type="Reanalysis", file=env_file, dictionary="ECMWF", + pressure_conversion_factor="hPa", ) env.max_expected_height = 2000 diff --git a/tests/fixtures/environment/environment_fixtures.py b/tests/fixtures/environment/environment_fixtures.py index 818f434c7..5a3f3cd64 100644 --- a/tests/fixtures/environment/environment_fixtures.py +++ b/tests/fixtures/environment/environment_fixtures.py @@ -111,6 +111,7 @@ def environment_spaceport_america_2023(): type="Reanalysis", file="data/weather/spaceport_america_pressure_levels_2023_hourly.nc", dictionary="ECMWF", + pressure_conversion_factor="hPa", ) env.max_expected_height = 6000 diff --git a/tests/fixtures/surfaces/surface_fixtures.py b/tests/fixtures/surfaces/surface_fixtures.py index c48627478..3819595dc 100644 --- a/tests/fixtures/surfaces/surface_fixtures.py +++ b/tests/fixtures/surfaces/surface_fixtures.py @@ -1,11 +1,14 @@ import pytest from rocketpy.rocket.aero_surface import ( + EllipticalFin, EllipticalFins, + FreeFormFin, FreeFormFins, NoseCone, RailButtons, Tail, + TrapezoidalFin, TrapezoidalFins, ) @@ -69,6 +72,29 @@ def calisto_trapezoidal_fins(): ) +@pytest.fixture +def calisto_trapezoidal_fin(): + """A single trapezoidal fin based on Calisto dimensions. + + Returns + ------- + rocketpy.TrapezoidalFin + A single trapezoidal fin. + """ + return TrapezoidalFin( + angular_position=0, + span=0.100, + root_chord=0.120, + tip_chord=0.040, + rocket_radius=0.0635, + name="calisto_trapezoidal_fin", + cant_angle=0, + sweep_length=None, + sweep_angle=None, + airfoil=None, + ) + + @pytest.fixture def calisto_trapezoidal_fins_custom_sweep_length(): """The trapezoidal fins of the Calisto rocket with @@ -130,6 +156,23 @@ def calisto_free_form_fins(): ) +@pytest.fixture +def calisto_free_form_fin(): + """A single free-form fin based on Calisto-like dimensions. + + Returns + ------- + rocketpy.FreeFormFin + A single free-form fin. + """ + return FreeFormFin( + angular_position=0, + shape_points=[(0, 0), (0.08, 0.1), (0.12, 0.1), (0.12, 0)], + rocket_radius=0.0635, + name="calisto_free_form_fin", + ) + + @pytest.fixture def calisto_rail_buttons(): """The rail buttons of the Calisto rocket. @@ -157,3 +200,23 @@ def elliptical_fin_set(): airfoil=None, name="Test Elliptical Fins", ) + + +@pytest.fixture +def calisto_elliptical_fin(): + """A single elliptical fin based on Calisto-like dimensions. + + Returns + ------- + rocketpy.EllipticalFin + A single elliptical fin. + """ + return EllipticalFin( + angular_position=0, + span=0.100, + root_chord=0.120, + rocket_radius=0.0635, + cant_angle=0, + airfoil=None, + name="calisto_elliptical_fin", + ) diff --git a/tests/integration/environment/test_environment.py b/tests/integration/environment/test_environment.py index 5802650dc..271f7574e 100644 --- a/tests/integration/environment/test_environment.py +++ b/tests/integration/environment/test_environment.py @@ -1,5 +1,5 @@ import time -from datetime import date, datetime, timezone +from datetime import date, datetime, timedelta, timezone from unittest.mock import patch import numpy as np @@ -45,6 +45,7 @@ def test_era5_atmosphere(mock_show, example_spaceport_env): # pylint: disable=u type="Reanalysis", file="data/weather/SpaceportAmerica_2018_ERA-5.nc", dictionary="ECMWF", + pressure_conversion_factor="hPa", ) assert example_spaceport_env.all_info() is None @@ -146,6 +147,7 @@ def test_wind_plots_wrapping_direction(mock_show, example_plain_env): # pylint: assert example_plain_env.plots.atmospheric_model() is None +@pytest.mark.slow @pytest.mark.parametrize( "model_name", [ @@ -195,6 +197,42 @@ def test_gfs_atmosphere(mock_show, example_spaceport_env): # pylint: disable=un assert example_spaceport_env.all_info() is None +@pytest.mark.slow +@patch("matplotlib.pyplot.show") +def test_aigfs_atmosphere(mock_show, example_spaceport_env): # pylint: disable=unused-argument + """Tests the Forecast model with the AIGFS file. + + Parameters + ---------- + mock_show : mock + Mock object to replace matplotlib.pyplot.show() method. + example_spaceport_env : rocketpy.Environment + Example environment object to be tested. + """ + example_spaceport_env.set_atmospheric_model(type="Forecast", file="AIGFS") + assert example_spaceport_env.all_info() is None + + +@pytest.mark.slow +@patch("matplotlib.pyplot.show") +def test_hrrr_atmosphere(mock_show, example_spaceport_env): # pylint: disable=unused-argument + """Tests the Forecast model with the HRRR file. + + Parameters + ---------- + mock_show : mock + Mock object to replace matplotlib.pyplot.show() method. + example_spaceport_env : rocketpy.Environment + Example environment object to be tested. + """ + # Sometimes the HRRR latest-model can fail due to not having at least 24 + # hours in the future in the forecast, so we try with 12 hours in the future + # only. + example_spaceport_env.set_date(datetime.now() + timedelta(hours=12)) + example_spaceport_env.set_atmospheric_model(type="Forecast", file="HRRR") + assert example_spaceport_env.all_info() is None + + @pytest.mark.slow @patch("matplotlib.pyplot.show") def test_nam_atmosphere(mock_show, example_spaceport_env): # pylint: disable=unused-argument diff --git a/tests/integration/simulation/test_monte_carlo.py b/tests/integration/simulation/test_monte_carlo.py index 4b1b82392..98af2431d 100644 --- a/tests/integration/simulation/test_monte_carlo.py +++ b/tests/integration/simulation/test_monte_carlo.py @@ -236,3 +236,30 @@ def invalid_data_collector(flight): monte_carlo_calisto.simulate(number_of_simulations=10, append=False) finally: _post_test_file_cleanup() + + +@pytest.mark.slow +def test_monte_carlo_simulate_convergence(monte_carlo_calisto): + """Tests the simulate_convergence method of the MonteCarlo class. + + Parameters + ---------- + monte_carlo_calisto : MonteCarlo + The MonteCarlo object, this is a pytest fixture. + """ + try: + ci_history = monte_carlo_calisto.simulate_convergence( + target_attribute="apogee", + target_confidence=0.95, + tolerance=5.0, + max_simulations=20, + batch_size=5, + parallel=False, + ) + + assert isinstance(ci_history, list) + assert all(isinstance(width, float) for width in ci_history) + assert len(ci_history) >= 1 + assert monte_carlo_calisto.num_of_loaded_sims <= 20 + finally: + _post_test_file_cleanup() diff --git a/tests/integration/test_encoding.py b/tests/integration/test_encoding.py index 38d5aef41..cd2de7ee1 100644 --- a/tests/integration/test_encoding.py +++ b/tests/integration/test_encoding.py @@ -273,7 +273,7 @@ def test_trapezoidal_fins_encoder(fin_name, request): 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: + 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) diff --git a/tests/integration/test_rocket.py b/tests/integration/test_rocket.py index c47096617..e36dd933e 100644 --- a/tests/integration/test_rocket.py +++ b/tests/integration/test_rocket.py @@ -19,7 +19,7 @@ def test_airfoil( calisto.add_surfaces(calisto_tail, -1.313) test_rocket.add_trapezoidal_fins( - 2, + 3, span=0.100, root_chord=0.120, tip_chord=0.040, @@ -28,7 +28,7 @@ def test_airfoil( name="NACA0012", ) test_rocket.add_trapezoidal_fins( - 2, + 3, span=0.100, root_chord=0.120, tip_chord=0.040, diff --git a/tests/unit/environment/test_environment.py b/tests/unit/environment/test_environment.py index beb6d5ac6..222eb9a2d 100644 --- a/tests/unit/environment/test_environment.py +++ b/tests/unit/environment/test_environment.py @@ -336,8 +336,8 @@ def test_resolve_dictionary_keeps_compatible_mapping(example_plain_env): assert resolved is gfs_mapping -def test_resolve_dictionary_falls_back_to_legacy_mapping(example_plain_env): - """Fallback to a compatible built-in mapping for legacy NOMADS-style files.""" +def test_resolve_dictionary_falls_back_to_first_compatible_mapping(example_plain_env): + """Fallback to the first compatible built-in mapping for legacy-style files.""" thredds_gfs_mapping = example_plain_env._Environment__weather_model_map.get("GFS") dataset = _DummyDataset( [ @@ -356,7 +356,6 @@ def test_resolve_dictionary_falls_back_to_legacy_mapping(example_plain_env): thredds_gfs_mapping, dataset ) - # Explicit legacy mappings should be preferred over unrelated model mappings. assert resolved == example_plain_env._Environment__weather_model_map.get( "GFS_LEGACY" ) @@ -578,9 +577,10 @@ def test_set_atmospheric_model_normalizes_shortcut_case_for_forecast(example_pla called_arguments = {} - def fake_process_forecast_reanalysis(dataset, dictionary): + def fake_process_forecast_reanalysis(dataset, dictionary, conversion_factor): called_arguments["dataset"] = dataset called_arguments["dictionary"] = dictionary + called_arguments["conversion_factr"] = conversion_factor environment.process_forecast_reanalysis = fake_process_forecast_reanalysis @@ -602,3 +602,52 @@ def test_set_atmospheric_model_raises_for_unknown_model_type(example_plain_env): # Act / Assert with pytest.raises(ValueError, match="Unknown model type"): environment.set_atmospheric_model(type="unknown_type") + + +@pytest.mark.parametrize("shortcut_name", ["AIGFS", "HRRR"]) +def test_forecast_shortcut_and_dictionary_are_case_insensitive( + monkeypatch, shortcut_name +): + """Ensure forecast shortcuts and built-in dictionaries ignore input casing.""" + # Arrange + env = Environment(date=(2026, 3, 17, 12), latitude=32.99, longitude=-106.97) + + sentinel_dataset = object() + env._Environment__atm_type_file_to_function_map["forecast"][shortcut_name] = ( + lambda: sentinel_dataset + ) + + captured = {} + + def fake_process_forecast_reanalysis(file, dictionary, conversion_factor): + captured["file"] = file + captured["dictionary"] = dictionary + captured["conversion_factor"] = conversion_factor + + monkeypatch.setattr( + env, "process_forecast_reanalysis", fake_process_forecast_reanalysis + ) + monkeypatch.setattr(env, "calculate_density_profile", lambda: None) + monkeypatch.setattr(env, "calculate_speed_of_sound_profile", lambda: None) + monkeypatch.setattr(env, "calculate_dynamic_viscosity", lambda: None) + + # Act + env.set_atmospheric_model( + type="forecast", + file=shortcut_name.lower(), + dictionary=shortcut_name.lower(), + ) + + # Assert + expected_dictionary = env._Environment__weather_model_map.get(shortcut_name) + assert captured["file"] is sentinel_dataset + assert captured["dictionary"] == expected_dictionary + assert env.atmospheric_model_file == shortcut_name + assert env.atmospheric_model_dict == expected_dictionary + + +def test_weather_model_mapping_get_is_case_insensitive(): + """Ensure built-in mapping names are resolved regardless of casing.""" + mapping = WeatherModelMapping() + assert mapping.get("aigfs") == mapping.get("AIGFS") + assert mapping.get("ecmwf_v0") == mapping.get("ECMWF_v0") diff --git a/tests/unit/rocket/aero_surface/test_fin_geometry.py b/tests/unit/rocket/aero_surface/test_fin_geometry.py new file mode 100644 index 000000000..9531a5c98 --- /dev/null +++ b/tests/unit/rocket/aero_surface/test_fin_geometry.py @@ -0,0 +1,317 @@ +"""Unit tests for fin geometry strategy classes.""" + +import numpy as np +import pytest + +from rocketpy.rocket.aero_surface.fins._geometry import ( + _EllipticalGeometry, + _FreeFormGeometry, + _TrapezoidalGeometry, +) + + +def test_trapezoidal_geometry_evaluate_geometrical_parameters( + calisto_trapezoidal_fin, +): + """Ensure trapezoidal geometry populates the derived fin parameters.""" + # Arrange + geometry = calisto_trapezoidal_fin.geometry + + # Act + geometry.evaluate_geometrical_parameters() + + # Assert + owner = calisto_trapezoidal_fin + expected_area = (owner.root_chord + owner.tip_chord) * owner.span / 2 + expected_aspect_ratio = 2 * owner.span**2 / expected_area + expected_gamma_c = np.arctan( + (geometry.sweep_length + 0.5 * owner.tip_chord - 0.5 * owner.root_chord) + / owner.span + ) + expected_mid_aerodynamic_span = ( + owner.span + / 3 + * (owner.root_chord + 2 * owner.tip_chord) + / (owner.root_chord + owner.tip_chord) + ) + tau = (owner.span + owner.rocket_radius) / owner.rocket_radius + lambda_ = owner.tip_chord / owner.root_chord + expected_roll_constant = ( + (owner.root_chord + 3 * owner.tip_chord) * owner.span**3 + + 4 + * (owner.root_chord + 2 * owner.tip_chord) + * owner.rocket_radius + * owner.span**2 + + 6 * (owner.root_chord + owner.tip_chord) * owner.span * owner.rocket_radius**2 + ) / 12 + expected_lift_interference_factor = 1 + 1 / tau + expected_roll_damping_factor = 1 + ( + ((tau - lambda_) / tau) - ((1 - lambda_) / (tau - 1)) * np.log(tau) + ) / ( + ((tau + 1) * (tau - lambda_)) / 2 + - ((1 - lambda_) * (tau**3 - 1)) / (3 * (tau - 1)) + ) + expected_roll_forcing_factor = (1 / np.pi**2) * ( + (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) + + (np.pi * (tau**2 + 1) ** 2 / (tau**2 * (tau - 1) ** 2)) + * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + - (2 * np.pi * (tau + 1)) / (tau * (tau - 1)) + + ((tau**2 + 1) ** 2 / (tau**2 * (tau - 1) ** 2)) + * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2 + - (4 * (tau + 1)) / (tau * (tau - 1)) * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau)) + ) + + assert isinstance(geometry, _TrapezoidalGeometry) + assert owner.Af == pytest.approx(expected_area) + assert owner.AR == pytest.approx(expected_aspect_ratio) + assert owner.gamma_c == pytest.approx(expected_gamma_c) + assert owner.Yma == pytest.approx(expected_mid_aerodynamic_span) + assert owner.roll_geometrical_constant == pytest.approx(expected_roll_constant) + assert owner.tau == pytest.approx(tau) + assert owner.lift_interference_factor == pytest.approx( + expected_lift_interference_factor + ) + assert owner.roll_damping_interference_factor == pytest.approx( + expected_roll_damping_factor + ) + assert owner.roll_forcing_interference_factor == pytest.approx( + expected_roll_forcing_factor + ) + + +def test_trapezoidal_geometry_evaluate_shape_sets_expected_points( + calisto_trapezoidal_fin, +): + """Ensure trapezoidal geometry shape points match the configured sweep.""" + # Arrange + geometry = calisto_trapezoidal_fin.geometry + + # Act + geometry.evaluate_shape() + + # Assert + np.testing.assert_allclose( + calisto_trapezoidal_fin.shape_vec[0], + np.array([0.0, 0.08, 0.12, 0.12]), + ) + np.testing.assert_allclose( + calisto_trapezoidal_fin.shape_vec[1], + np.array([0.0, 0.1, 0.1, 0.0]), + ) + + +def test_trapezoidal_geometry_get_data_returns_inputs_and_outputs( + calisto_trapezoidal_fin, +): + """Ensure trapezoidal geometry serialization includes optional outputs.""" + # Arrange + geometry = calisto_trapezoidal_fin.geometry + + # Act + geometry.evaluate_geometrical_parameters() + data_without_outputs = geometry.get_data() + data_with_outputs = geometry.get_data(include_outputs=True) + + # Assert + assert data_without_outputs["tip_chord"] == pytest.approx( + calisto_trapezoidal_fin.tip_chord + ) + assert data_without_outputs["sweep_length"] == pytest.approx( + calisto_trapezoidal_fin.sweep_length + ) + assert data_without_outputs["sweep_angle"] is None + assert set(data_with_outputs) >= { + "tip_chord", + "sweep_length", + "sweep_angle", + "shape_vec", + "Af", + "AR", + "gamma_c", + "Yma", + "roll_geometrical_constant", + "tau", + "lift_interference_factor", + "roll_damping_interference_factor", + "roll_forcing_interference_factor", + } + np.testing.assert_allclose( + data_with_outputs["shape_vec"][0], calisto_trapezoidal_fin.shape_vec[0] + ) + + +def test_elliptical_geometry_evaluate_geometrical_parameters( + calisto_elliptical_fin, +): + """Ensure elliptical geometry populates the derived fin parameters.""" + # Arrange + geometry = calisto_elliptical_fin.geometry + + # Act + geometry.evaluate_geometrical_parameters() + + # Assert + owner = calisto_elliptical_fin + expected_area = np.pi * owner.root_chord / 2 * owner.span / 2 + expected_aspect_ratio = 2 * owner.span**2 / expected_area + expected_mid_aerodynamic_span = ( + owner.span / (3 * np.pi) * np.sqrt(9 * np.pi**2 - 64) + ) + expected_roll_constant = ( + owner.root_chord + * owner.span + * ( + 3 * np.pi * owner.span**2 + + 32 * owner.rocket_radius * owner.span + + 12 * np.pi * owner.rocket_radius**2 + ) + / 48 + ) + tau = (owner.span + owner.rocket_radius) / owner.rocket_radius + expected_lift_interference_factor = 1 + 1 / tau + + assert isinstance(geometry, _EllipticalGeometry) + assert owner.Af == pytest.approx(expected_area) + assert owner.AR == pytest.approx(expected_aspect_ratio) + assert owner.gamma_c == pytest.approx(0) + assert owner.Yma == pytest.approx(expected_mid_aerodynamic_span) + assert owner.roll_geometrical_constant == pytest.approx(expected_roll_constant) + assert owner.tau == pytest.approx(tau) + assert owner.lift_interference_factor == pytest.approx( + expected_lift_interference_factor + ) + assert owner.roll_damping_interference_factor > 0 + assert owner.roll_forcing_interference_factor > 0 + + +def test_elliptical_geometry_evaluate_shape_sets_expected_points( + calisto_elliptical_fin, +): + """Ensure elliptical geometry evaluates the expected semi-ellipse shape.""" + # Arrange + geometry = calisto_elliptical_fin.geometry + + # Act + geometry.evaluate_shape() + + # Assert + angles = np.arange(0, 180, 5) + expected_x = calisto_elliptical_fin.root_chord / 2 + ( + calisto_elliptical_fin.root_chord / 2 * np.cos(np.radians(angles)) + ) + expected_y = calisto_elliptical_fin.span * np.sin(np.radians(angles)) + np.testing.assert_allclose(calisto_elliptical_fin.shape_vec[0], expected_x) + np.testing.assert_allclose(calisto_elliptical_fin.shape_vec[1], expected_y) + + +def test_elliptical_geometry_get_data_returns_expected_outputs( + calisto_elliptical_fin, +): + """Ensure elliptical geometry serialization includes optional outputs.""" + # Arrange + geometry = calisto_elliptical_fin.geometry + + # Act + geometry.evaluate_geometrical_parameters() + data_without_outputs = geometry.get_data() + data_with_outputs = geometry.get_data(include_outputs=True) + + # Assert + assert data_without_outputs == {} + assert data_with_outputs["Af"] == pytest.approx(calisto_elliptical_fin.Af) + assert data_with_outputs["AR"] == pytest.approx(calisto_elliptical_fin.AR) + assert data_with_outputs["gamma_c"] == pytest.approx(calisto_elliptical_fin.gamma_c) + assert data_with_outputs["Yma"] == pytest.approx(calisto_elliptical_fin.Yma) + assert data_with_outputs["roll_geometrical_constant"] == pytest.approx( + calisto_elliptical_fin.roll_geometrical_constant + ) + assert data_with_outputs["tau"] == pytest.approx(calisto_elliptical_fin.tau) + + +def test_free_form_geometry_infer_dimensions_warns_on_jagged_shape(): + """Ensure jagged free-form fins emit a warning while dimensions are inferred.""" + # Arrange + shape_points = [(0, 0), (0.05, 0.1), (0.06, 0.05), (0.09, 0.07), (0.12, 0)] + + # Act + with pytest.warns(UserWarning, match="Jagged fin shape detected"): + root_chord, span = _FreeFormGeometry.infer_dimensions(shape_points) + + # Assert + assert root_chord == pytest.approx(0.12) + assert span == pytest.approx(0.1) + + +def test_free_form_geometry_evaluate_geometrical_parameters(calisto_free_form_fin): + """Ensure free-form geometry populates the derived fin parameters.""" + # Arrange + geometry = calisto_free_form_fin.geometry + + # Act + geometry.evaluate_geometrical_parameters() + + # Assert + owner = calisto_free_form_fin + assert isinstance(geometry, _FreeFormGeometry) + assert owner.Af == pytest.approx(0.008) + assert owner.AR == pytest.approx(2.5) + assert owner.gamma_c > 0 + assert owner.Yma > 0 + assert owner.mac_length > 0 + assert owner.mac_lead >= 0 + assert owner.roll_geometrical_constant > 0 + assert owner.tau > 0 + assert owner.lift_interference_factor > 1 + assert owner.roll_damping_interference_factor > 1 + assert owner.roll_forcing_interference_factor > 0 + + +def test_free_form_geometry_evaluate_shape_sets_expected_points( + calisto_free_form_fin, +): + """Ensure free-form geometry exposes the configured shape points.""" + # Arrange + geometry = calisto_free_form_fin.geometry + + # Act + geometry.evaluate_shape() + + # Assert + np.testing.assert_allclose( + calisto_free_form_fin.shape_vec[0], + np.array([0.0, 0.08, 0.12, 0.12]), + ) + np.testing.assert_allclose( + calisto_free_form_fin.shape_vec[1], + np.array([0.0, 0.1, 0.1, 0.0]), + ) + + +def test_free_form_geometry_get_data_returns_expected_outputs( + calisto_free_form_fin, +): + """Ensure free-form geometry serialization includes optional outputs.""" + # Arrange + geometry = calisto_free_form_fin.geometry + + # Act + geometry.evaluate_geometrical_parameters() + data_without_outputs = geometry.get_data() + data_with_outputs = geometry.get_data(include_outputs=True) + + # Assert + assert data_without_outputs == { + "shape_points": [(0, 0), (0.08, 0.1), (0.12, 0.1), (0.12, 0)], + } + assert data_with_outputs["shape_points"] == calisto_free_form_fin.shape_points + assert data_with_outputs["Af"] == pytest.approx(calisto_free_form_fin.Af) + assert data_with_outputs["AR"] == pytest.approx(calisto_free_form_fin.AR) + assert data_with_outputs["gamma_c"] == pytest.approx(calisto_free_form_fin.gamma_c) + assert data_with_outputs["Yma"] == pytest.approx(calisto_free_form_fin.Yma) + assert data_with_outputs["mac_length"] == pytest.approx( + calisto_free_form_fin.mac_length + ) + assert data_with_outputs["mac_lead"] == pytest.approx( + calisto_free_form_fin.mac_lead + ) diff --git a/tests/unit/rocket/aero_surface/test_individual_fins.py b/tests/unit/rocket/aero_surface/test_individual_fins.py new file mode 100644 index 000000000..d232e0772 --- /dev/null +++ b/tests/unit/rocket/aero_surface/test_individual_fins.py @@ -0,0 +1,424 @@ +"""Unit tests for individual fin classes.""" + +from unittest.mock import patch + +import numpy as np +import pytest + +from rocketpy import ( + EllipticalFin, + FreeFormFin, + Rocket, + TrapezoidalFin, + TrapezoidalFins, +) +from rocketpy.mathutils.vector_matrix import Vector + + +@pytest.mark.parametrize( + "fixture_name,expected_class", + [ + ("calisto_trapezoidal_fin", TrapezoidalFin), + ("calisto_elliptical_fin", EllipticalFin), + ("calisto_free_form_fin", FreeFormFin), + ], +) +def test_individual_fin_info_returns_none(request, fixture_name, expected_class): + """Ensure info() executes for all individual fin classes.""" + # Arrange + fin = request.getfixturevalue(fixture_name) + + # Act + result = fin.info() + + # Assert + assert isinstance(fin, expected_class) + assert result is None + + +@patch("matplotlib.pyplot.show") +@pytest.mark.parametrize( + "fixture_name", + [ + "calisto_trapezoidal_fin", + "calisto_elliptical_fin", + "calisto_free_form_fin", + ], +) +def test_individual_fin_draw_returns_none(mock_show, request, fixture_name): # pylint: disable=unused-argument + """Ensure draw() executes for all individual fin classes.""" + # Arrange + fin = request.getfixturevalue(fixture_name) + + # Act + result = fin.draw(filename=None) + + # Assert + assert result is None + + +@pytest.mark.parametrize( + "fixture_name", + [ + "calisto_trapezoidal_fin", + "calisto_elliptical_fin", + "calisto_free_form_fin", + ], +) +def test_individual_fin_angular_position_updates_radians(request, fixture_name): + """Ensure angular_position setter updates angular_position_rad.""" + # Arrange + fin = request.getfixturevalue(fixture_name) + + # Act + fin.angular_position = 45 + + # Assert + assert fin.angular_position == 45 + np.testing.assert_allclose(fin.angular_position_rad, np.pi / 4) + + +def test_trapezoidal_fin_setters_update_geometry(calisto_trapezoidal_fin): + """Ensure trapezoidal fin geometry setters update exposed values.""" + # Arrange + fin = calisto_trapezoidal_fin + + # Act + fin.tip_chord = 0.05 + fin.sweep_angle = 12.0 + fin.sweep_length = 0.03 + + # Assert + np.testing.assert_allclose(fin.tip_chord, 0.05) + np.testing.assert_allclose(fin.sweep_angle, 12.0) + np.testing.assert_allclose(fin.sweep_length, 0.03) + + +def test_individual_fin_rocket_diameter_aliases_are_kept_in_sync( + calisto_trapezoidal_fin, +): + """Ensure rocket_diameter is canonical and old aliases remain compatible.""" + # Arrange + fin = calisto_trapezoidal_fin + + # Act + fin.rocket_diameter = 0.15 + + # Assert + np.testing.assert_allclose(fin.rocket_diameter, 0.15) + np.testing.assert_allclose(fin.diameter, 0.15) + np.testing.assert_allclose(fin.d, 0.15) + np.testing.assert_allclose(fin.rocket_radius, 0.075) + np.testing.assert_allclose(fin.reference_length, 0.15) + + # Act + fin.d = 0.20 + + # Assert + np.testing.assert_allclose(fin.rocket_diameter, 0.20) + np.testing.assert_allclose(fin.diameter, 0.20) + np.testing.assert_allclose(fin.d, 0.20) + np.testing.assert_allclose(fin.rocket_radius, 0.10) + np.testing.assert_allclose(fin.reference_length, 0.20) + + +def test_individual_fin_reference_area_and_ref_area_alias_are_kept_in_sync( + calisto_trapezoidal_fin, +): + """Ensure reference_area is canonical and ref_area remains compatible.""" + # Arrange + fin = calisto_trapezoidal_fin + + # Act + fin.reference_area = 0.123 + + # Assert + np.testing.assert_allclose(fin.reference_area, 0.123) + np.testing.assert_allclose(fin.ref_area, 0.123) + + # Act + fin.ref_area = 0.456 + + # Assert + np.testing.assert_allclose(fin.reference_area, 0.456) + np.testing.assert_allclose(fin.ref_area, 0.456) + + +def test_individual_fin_to_dict_include_outputs_exposes_diameter_aliases( + calisto_trapezoidal_fin, +): + """Ensure output serialization exposes canonical and alias diameter keys.""" + # Arrange + fin = calisto_trapezoidal_fin + + # Act + data = fin.to_dict(include_outputs=True) + + # Assert + np.testing.assert_allclose(data["rocket_diameter"], fin.rocket_diameter) + np.testing.assert_allclose(data["diameter"], fin.rocket_diameter) + np.testing.assert_allclose(data["d"], fin.rocket_diameter) + np.testing.assert_allclose(data["reference_area"], fin.reference_area) + np.testing.assert_allclose(data["ref_area"], fin.reference_area) + + +def test_trapezoidal_fin_rejects_inconsistent_sweep_inputs(): + """Ensure trapezoidal fin rejects sweep_length with sweep_angle together.""" + # Arrange / Act / Assert + with pytest.raises( + ValueError, match="Cannot use sweep_length and sweep_angle together" + ): + TrapezoidalFin( + angular_position=0, + root_chord=0.12, + tip_chord=0.04, + span=0.1, + rocket_radius=0.0635, + sweep_length=0.02, + sweep_angle=10.0, + ) + + +def test_free_form_fin_shape_points_property(calisto_free_form_fin): + """Ensure free-form fin exposes the original shape points.""" + # Arrange + fin = calisto_free_form_fin + + # Act + shape_points = fin.shape_points + + # Assert + assert shape_points == [(0, 0), (0.08, 0.1), (0.12, 0.1), (0.12, 0)] + + +@pytest.mark.parametrize( + "fixture_name,required_keys", + [ + ( + "calisto_trapezoidal_fin", + { + "angular_position", + "root_chord", + "span", + "rocket_radius", + "cant_angle", + "airfoil", + "name", + "tip_chord", + "sweep_length", + "sweep_angle", + }, + ), + ( + "calisto_elliptical_fin", + { + "angular_position", + "root_chord", + "span", + "rocket_radius", + "cant_angle", + "airfoil", + "name", + }, + ), + ( + "calisto_free_form_fin", + { + "angular_position", + "rocket_radius", + "cant_angle", + "airfoil", + "name", + "shape_points", + }, + ), + ], +) +def test_individual_fin_to_dict_contains_expected_keys( + request, fixture_name, required_keys +): + """Ensure to_dict for each individual fin exposes expected input keys.""" + # Arrange + fin = request.getfixturevalue(fixture_name) + + # Act + data = fin.to_dict() + + # Assert + assert required_keys.issubset(data.keys()) + + +@pytest.mark.parametrize( + "fixture_name,fin_class,comparisons", + [ + ( + "calisto_trapezoidal_fin", + TrapezoidalFin, + ["angular_position", "root_chord", "tip_chord", "span", "rocket_radius"], + ), + ( + "calisto_elliptical_fin", + EllipticalFin, + ["angular_position", "root_chord", "span", "rocket_radius"], + ), + ( + "calisto_free_form_fin", + FreeFormFin, + ["angular_position", "rocket_radius"], + ), + ], +) +def test_individual_fin_from_dict_roundtrip( + request, fixture_name, fin_class, comparisons +): + """Ensure each individual fin can be reconstructed with from_dict.""" + # Arrange + fin = request.getfixturevalue(fixture_name) + data = fin.to_dict() + + # Act + reconstructed = fin_class.from_dict(data) + + # Assert + assert isinstance(reconstructed, fin_class) + for field in comparisons: + np.testing.assert_allclose(getattr(reconstructed, field), getattr(fin, field)) + + if fin_class is FreeFormFin: + assert reconstructed.shape_points == fin.shape_points + + +def test_trapezoidal_fin_from_dict_roundtrip_preserves_sweep_length(): + """Ensure TrapezoidalFin round-trip preserves non-default sweep geometry.""" + # Arrange + original = TrapezoidalFin( + angular_position=0, + root_chord=0.12, + tip_chord=0.04, + span=0.1, + rocket_radius=0.0635, + cant_angle=0, + sweep_angle=15.0, + name="roundtrip_trapezoidal_fin", + ) + data = original.to_dict() + + # Act + reconstructed = TrapezoidalFin.from_dict(data) + + # Assert + np.testing.assert_allclose(reconstructed.sweep_length, original.sweep_length) + + +def test_calisto_finset_vs_four_individual_fins_close(): + """Ensure a 4-fin set and 4 individual fins produce close aerodynamics. + + Notes + ----- + A fin set model includes finite-set lift correction for the number of fins. + For 4 fins, this correction is equivalent to scaling the sum of 4 + individual-fin lift derivatives by 1/2. + """ + # Arrange + finset_rocket = Rocket( + radius=0.0635, + 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", + ) + finset_rocket.add_surfaces( + TrapezoidalFins( + n=4, + span=0.100, + root_chord=0.120, + tip_chord=0.040, + rocket_radius=0.0635, + name="calisto_trapezoidal_fins", + cant_angle=0, + sweep_length=None, + sweep_angle=None, + airfoil=None, + ), + -1.168, + ) + + individual_fins_rocket = Rocket( + radius=0.0635, + 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", + ) + + individual_fins = [ + TrapezoidalFin( + angular_position=angle, + root_chord=0.120, + tip_chord=0.040, + span=0.100, + rocket_radius=0.0635, + name=f"calisto_trapezoidal_fin_{i}", + cant_angle=0, + sweep_length=None, + sweep_angle=None, + airfoil=None, + ) + for i, angle in enumerate((0, 90, 180, 270), start=1) + ] + individual_fins_rocket.add_surfaces(individual_fins, [-1.168] * 4) + + mach_grid = np.linspace(0, 2, 21) + + # Act + cp_finset = finset_rocket.cp_position(mach_grid) + cp_individual = individual_fins_rocket.cp_position(mach_grid) + clalpha_finset = finset_rocket.total_lift_coeff_der(mach_grid) + clalpha_individual = individual_fins_rocket.total_lift_coeff_der(mach_grid) + lift_correction = TrapezoidalFins.fin_num_correction(4) / 4 + clalpha_individual_corrected = np.array(clalpha_individual) * lift_correction + + # Assert + np.testing.assert_allclose(cp_individual, cp_finset, rtol=1e-6, atol=1e-6) + np.testing.assert_allclose(clalpha_individual_corrected, clalpha_finset) + + +@pytest.mark.parametrize( + "position_input", + [ + (0.02, -0.01, -1.2), + Vector([0.02, -0.01, -1.2]), + ], +) +def test_add_individual_fin_accepts_full_3d_position(position_input): + """Ensure individual fins accept full (x, y, z) position inputs.""" + # Arrange + rocket = Rocket( + radius=0.0635, + 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", + ) + fin = TrapezoidalFin( + angular_position=30, + root_chord=0.120, + tip_chord=0.040, + span=0.100, + rocket_radius=0.0635, + cant_angle=0, + name="position_test_fin", + ) + + # Act + rocket.add_surfaces(fin, position_input) + stored_position = rocket.aerodynamic_surfaces[0].position + + # Assert + assert stored_position == Vector([0.02, -0.01, -1.2]) diff --git a/tests/unit/test_tools.py b/tests/unit/test_tools.py index fcf67ad37..a1e96eb9e 100644 --- a/tests/unit/test_tools.py +++ b/tests/unit/test_tools.py @@ -1,6 +1,7 @@ import numpy as np import pytest +from rocketpy import Environment from rocketpy.tools import ( calculate_cubic_hermite_coefficients, euler313_to_quaternions, @@ -100,3 +101,39 @@ def test_tuple_handler(input_value, expected_output): def test_tuple_handler_exceptions(input_value, expected_exception): with pytest.raises(expected_exception): tuple_handler(input_value) + + +@pytest.mark.parametrize("pressure_conversion_factor", ["hPa", "mbar", "Pa", 100]) +def test_valid_pressure_conversion_factor(pressure_conversion_factor): + env = Environment( + gravity=9.81, + latitude=47.213476, + longitude=9.003336, + date=(2020, 2, 22, 13), + elevation=407, + ) + env.set_atmospheric_model( + type="Reanalysis", + file="data/weather/bella_lui_weather_data_ERA5.nc", + dictionary="ECMWF", + pressure_conversion_factor=pressure_conversion_factor, + ) + + +@pytest.mark.parametrize("pressure_conversion_factor", [-1, "mPa"]) +def test_invalid_pressure_conversion_factor(pressure_conversion_factor): + env = Environment( + gravity=9.81, + latitude=47.213476, + longitude=9.003336, + date=(2020, 2, 22, 13), + elevation=407, + ) + + with pytest.raises(ValueError): + env.set_atmospheric_model( + type="Reanalysis", + file="data/weather/bella_lui_weather_data_ERA5.nc", + dictionary="ECMWF", + pressure_conversion_factor=pressure_conversion_factor, + )