diff --git a/.flake8 b/.flake8 deleted file mode 100644 index 01efeb095..000000000 --- a/.flake8 +++ /dev/null @@ -1,5 +0,0 @@ -[flake8] -max-line-length = 88 -max-module-lines= 3000 -exclude = .git,__pycache__ -ignore = E203, W503 diff --git a/.github/workflows/test_pytest.yaml b/.github/workflows/test_pytest.yaml index ffe524718..c3e06afe9 100644 --- a/.github/workflows/test_pytest.yaml +++ b/.github/workflows/test_pytest.yaml @@ -6,6 +6,8 @@ on: paths: - "**.py" - ".github/**" + - "pyproject.toml" + - "requirements*" defaults: run: @@ -16,8 +18,11 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - os: [ubuntu-latest, macos-latest, windows-latest] + os: [ubuntu-latest, windows-latest] python-version: [3.8, 3.12] + include: + - os: macos-latest + python-version: 3.12 env: OS: ${{ matrix.os }} PYTHON: ${{ matrix.python-version }} diff --git a/CHANGELOG.md b/CHANGELOG.md index 5a00ea820..869b49aa8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,17 +32,36 @@ and this project adheres to [Semantic Versioning](http://semver.org/). ### Added +- ENH: Exponential backoff decorator (fix #449) [#588](https://github.com/RocketPy-Team/RocketPy/pull/588) +- ENH: Function Validation Rework & Swap `np.searchsorted` to `bisect_left` [#582](https://github.com/RocketPy-Team/RocketPy/pull/582) +- ENH: Add new stability margin properties to Flight class [#572](https://github.com/RocketPy-Team/RocketPy/pull/572) - ENH: adds `Function.remove_outliers` method [#554](https://github.com/RocketPy-Team/RocketPy/pull/554) ### Changed +- BLD: Change setup.py to pyproject.toml [#589](https://github.com/RocketPy-Team/RocketPy/pull/589) +- DEP: delete deprecated rocketpy.tools.cached_property [#587](https://github.com/RocketPy-Team/RocketPy/pull/587) +- ENH: Flight simulation speed up [#581] (https://github.com/RocketPy-Team/RocketPy/pull/581) +- MNT: Modularize Rocket Draw [#580](https://github.com/RocketPy-Team/RocketPy/pull/580) - ENH: Initial solution from flights gets previous results [#568](https://github.com/RocketPy-Team/RocketPy/pull/568) +- DOC: Improvements of Environment docstring phrasing [#565](https://github.com/RocketPy-Team/RocketPy/pull/565) +- MNT: Refactor flight prints module [#579](https://github.com/RocketPy-Team/RocketPy/pull/579) - DOC: Convert CompareFlights example notebooks to .rst files [#576](https://github.com/RocketPy-Team/RocketPy/pull/576) +- MNT: Refactor inertia calculations using parallel axis theorem [#573] (https://github.com/RocketPy-Team/RocketPy/pull/573) - ENH: Optional argument to show the plot in Function.compare_plots [#563](https://github.com/RocketPy-Team/RocketPy/pull/563) ### Fixed + +- BUG: Optional Dependencies Naming in pyproject.toml. [#592](https://github.com/RocketPy-Team/RocketPy/pull/592) +- BUG: Swap rocket.total_mass.differentiate for motor.total_mass_flow rate [#585](https://github.com/RocketPy-Team/RocketPy/pull/585) - BUG: export_eng 'Motor' method would not work for liquid motors. [#559](https://github.com/RocketPy-Team/RocketPy/pull/559) +## [v1.2.2] - 2024-03-22 + +You can install this version by running `pip install rocketpy==1.2.2` + +- BUG: wrong rocket mass in parachute u dot method [#569](https://github.com/RocketPy-Team/RocketPy/pull/569) + ## [v1.2.1] - 2024-02-22 You can install this version by running `pip install rocketpy==1.2.1` diff --git a/docs/conf.py b/docs/conf.py index f2397dfa8..75f7a18a7 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -24,7 +24,7 @@ author = "RocketPy Team" # The full version, including alpha/beta/rc tags -release = "1.2.1" +release = "1.2.2" # -- General configuration --------------------------------------------------- diff --git a/docs/user/installation.rst b/docs/user/installation.rst index 16dec6886..8f4e0e183 100644 --- a/docs/user/installation.rst +++ b/docs/user/installation.rst @@ -19,7 +19,7 @@ If you want to choose a specific version to guarantee compatibility, you may ins .. code-block:: shell - pip install rocketpy==1.2.1 + pip install rocketpy==1.2.2 Optional Installation Method: ``conda`` @@ -49,7 +49,7 @@ Once you are done downloading/cloning RocketPy's repository, you can install it .. code-block:: shell - python setup.py install + python -m pip install . Development version @@ -76,7 +76,7 @@ Alternatively, you can clone RocketPy's repository, check out the branch named ` .. code-block:: shell - python setup.py install + python -m pip install -e . Requirements diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 000000000..350f7ae37 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,74 @@ +[project] +name = "rocketpy" +version = "1.2.2" +description="Advanced 6-DOF trajectory simulation for High-Power Rocketry." +dynamic = ["dependencies"] +readme = "README.md" +requires-python = ">=3.8" +authors = [ + {name = "Giovani Hidalgo Ceotto", email = "ghceotto@gmail.com"} +] +classifiers = [ + "Programming Language :: Python :: 3", + "License :: OSI Approved :: MIT License", + "Operating System :: OS Independent", + "Topic :: Rocket Flight Simulation :: Libraries :: Python Modules", +] + +[project.urls] +homepage = "https://rocketpy.org/" +documentation = "https://docs.rocketpy.org/" +repository = "https://github.com/RocketPy-Team/RocketPy" + +[build-system] +requires = ["setuptools"] +build-backend = "setuptools.build_meta" + +[tool.setuptools] +py-modules = ['rocketpy'] + +[tool.setuptools.dynamic] +dependencies = { file = ["requirements.txt"] } + +[project.optional-dependencies] +tests = [ + "pytest", + "pytest-coverage", + "black[jupyter]", + "flake8-black", + "flake8-pyproject", + "pandas", + "numericalunits==1.25", + "pylint", + "isort" +] + +env-analysis = [ + "windrose>=1.6.8", + "timezonefinder", + "jsonpickle", + "ipython", + "ipywidgets>=7.6.3" +] + +all = ["rocketpy[env-analysis]"] + +[tool.black] +line-length = 88 +include = '\.py$|\.ipynb$' +skip-string-normalization = true + +[tool.coverage.report] +# Regexes for lines to exclude from consideration +exclude_also = [ + #Don't complain about exceptions or warnings not being covered by tests + "warnings.warn*" +] + +[tool.flake8] +max-line-length = 88 +max-module-lines = 3000 +ignore = ['E203', 'W503'] +exclude = [ + '.git,__pycache__', +] diff --git a/requirements-tests.txt b/requirements-tests.txt index 42f5a8292..4fc2bc811 100644 --- a/requirements-tests.txt +++ b/requirements-tests.txt @@ -2,7 +2,8 @@ pytest pytest-coverage black[jupyter] flake8-black +flake8-pyproject pandas numericalunits==1.25 pylint -isort \ No newline at end of file +isort diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py index 546d446e8..1f50364bf 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -3,7 +3,7 @@ import re import warnings from collections import namedtuple -from datetime import datetime, timedelta +from datetime import datetime, timedelta, timezone import numpy as np import numpy.ma as ma @@ -13,6 +13,7 @@ from ..mathutils.function import Function, funcify_method from ..plots.environment_plots import _EnvironmentPlots from ..prints.environment_prints import _EnvironmentPrints +from ..tools import exponential_backoff try: import netCDF4 @@ -57,8 +58,7 @@ class Environment: Environment.datum : string The desired reference ellipsoid model, the following options are available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default - is "SIRGAS2000", then this model will be used if the user make some - typing mistake + is "SIRGAS2000". Environment.initial_east : float Launch site East UTM coordinate Environment.initial_north : float @@ -74,7 +74,7 @@ class Environment: Launch site E/W hemisphere Environment.elevation : float Launch site elevation. - Environment.date : datetime + Environment.datetime_date : datetime Date time of launch in UTC. Environment.local_date : datetime Date time of launch in the local time zone, defined by @@ -276,49 +276,70 @@ def __init__( timezone="UTC", max_expected_height=80000.0, ): - """Initialize Environment class, saving launch rail length, - launch date, location coordinates and elevation. Note that - by default the standard atmosphere is loaded until another + """Initializes the Environment class, capturing essential parameters of + the launch site, including the launch date, geographical coordinates, + and elevation. This class is designed to calculate crucial variables + for the Flight simulation, such as atmospheric air pressure, density, + and gravitational acceleration. + + Note that the default atmospheric model is the International Standard + Atmosphere as defined by ISO 2533 unless specified otherwise in + :meth:`Environment.set_atmospheric_model`. Parameters ---------- gravity : int, float, callable, string, array, optional Surface gravitational acceleration. Positive values point the - acceleration down. If None, the Somigliana formula is used to - date : array, optional - Array of length 4, stating (year, month, day, hour (UTC)) - of rocket launch. Must be given if a Forecast, Reanalysis + acceleration down. If None, the Somigliana formula is used. + See :meth:`Environment.set_gravity_model` for more information. + date : list or tuple, optional + List or tuple of length 4, stating (year, month, day, hour) in the + time zone of the parameter ``timezone``. + Alternatively, can be a ``datetime`` object specifying launch + date and time. The dates are stored as follows: + + - :attr:`Environment.local_date`: Local time of launch in + the time zone specified by the parameter ``timezone``. + + - :attr:`Environment.datetime_date`: UTC time of launch. + + Must be given if a Forecast, Reanalysis or Ensemble, will be set as an atmospheric model. + Default is None. + See :meth:`Environment.set_date` for more information. latitude : float, optional Latitude in degrees (ranging from -90 to 90) of rocket launch location. Must be given if a Forecast, Reanalysis or Ensemble will be used as an atmospheric model or if - Open-Elevation will be used to compute elevation. + Open-Elevation will be used to compute elevation. Positive + values correspond to the North. Default value is 0, which + corresponds to the equator. longitude : float, optional - Longitude in degrees (ranging from -180 to 360) of rocket + Longitude in degrees (ranging from -180 to 180) of rocket launch location. Must be given if a Forecast, Reanalysis or Ensemble will be used as an atmospheric model or if - Open-Elevation will be used to compute elevation. + Open-Elevation will be used to compute elevation. Positive + values correspond to the East. Default value is 0, which + corresponds to the Greenwich Meridian. elevation : float, optional Elevation of launch site measured as height above sea level in meters. Alternatively, can be set as 'Open-Elevation' which uses the Open-Elevation API to find elevation data. For this option, latitude and longitude must also be specified. Default value is 0. - datum : string + datum : string, optional The desired reference ellipsoidal model, the following options are available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default - is "SIRGAS2000", then this model will be used if the user make some - typing mistake. + is "SIRGAS2000". timezone : string, optional Name of the time zone. To see all time zones, import pytz and run - print(pytz.all_timezones). Default time zone is "UTC". + ``print(pytz.all_timezones)``. Default time zone is "UTC". max_expected_height : float, optional Maximum altitude in meters to keep weather data. The altitude must be above sea level (ASL). Especially useful for visualization. Can be altered as desired by doing `max_expected_height = number`. Depending on the atmospheric model, this value may be automatically - mofified. + modified. Returns ------- @@ -396,15 +417,57 @@ def set_date(self, date, timezone="UTC"): Parameters ---------- - date : Datetime - Datetime object specifying launch date and time. + date : list, tuple, datetime + List or tuple of length 4, stating (year, month, day, hour) in the + time zone of the parameter ``timezone``. See Notes for more + information. + Alternatively, can be a ``datetime`` object specifying launch + date and time. timezone : string, optional Name of the time zone. To see all time zones, import pytz and run - print(pytz.all_timezones). Default time zone is "UTC". + ``print(pytz.all_timezones)``. Default time zone is "UTC". Returns ------- None + + Notes + ----- + - If the ``date`` is given as a list or tuple, it should be in the same + time zone as specified by the ``timezone`` parameter. This local + time will be available in the attribute :attr:`Environment.local_date` + while the UTC time will be available in the attribute + :attr:`Environment.datetime_date`. + + - If the ``date`` is given as a ``datetime`` object without a time zone, + it will be assumed to be in the same time zone as specified by the + ``timezone`` parameter. However, if the ``datetime`` object has a time + zone specified in its ``tzinfo`` attribute, the ``timezone`` + parameter will be ignored. + + Examples + -------- + + Let's set the launch date as an list: + + >>> date = [2000, 1, 1, 13] # January 1st, 2000 at 13:00 UTC+1 + >>> env = Environment() + >>> env.set_date(date, timezone="Europe/Rome") + >>> print(env.datetime_date) # Get UTC time + 2000-01-01 12:00:00+00:00 + >>> print(env.local_date) + 2000-01-01 13:00:00+01:00 + + Now let's set the launch date as a ``datetime`` object: + + >>> from datetime import datetime + >>> date = datetime(2000, 1, 1, 13, 0, 0) + >>> env = Environment() + >>> env.set_date(date, timezone="Europe/Rome") + >>> print(env.datetime_date) # Get UTC time + 2000-01-01 12:00:00+00:00 + >>> print(env.local_date) + 2000-01-01 13:00:00+01:00 """ # Store date and configure time zone self.timezone = timezone @@ -458,23 +521,66 @@ def set_location(self, latitude, longitude): self.atmospheric_model_file, self.atmospheric_model_dict ) - # Return None - - def set_gravity_model(self, gravity): - """Sets the gravity model to be used in the simulation based on the - given user input to the gravity parameter. + def set_gravity_model(self, gravity=None): + """Defines the gravity model based on the given user input to the + gravity parameter. The gravity model is responsible for computing the + gravity acceleration at a given height above sea level in meters. Parameters ---------- - gravity : None or Function source - If None, the Somigliana formula is used to compute the gravity - acceleration. Otherwise, the user can provide a Function object - representing the gravity model. + gravity : int, float, callable, string, list, optional + The gravitational acceleration in m/s² to be used in the + simulation, this value is positive when pointing downwards. + The input type can be one of the following: + + - ``int`` or ``float``: The gravity acceleration is set as a\ + constant function with respect to height; + + - ``callable``: This callable should receive the height above\ + sea level in meters and return the gravity acceleration; + + - ``list``: The datapoints should be structured as\ + ``[(h_i,g_i), ...]`` where ``h_i`` is the height above sea\ + level in meters and ``g_i`` is the gravity acceleration in m/s²; + + - ``string``: The string should correspond to a path to a CSV file\ + containing the gravity acceleration data; + + - ``None``: The Somigliana formula is used to compute the gravity\ + acceleration. + + This parameter is used as a :class:`Function` object source, check\ + out the available input types for a more detailed explanation. Returns ------- Function Function object representing the gravity model. + + Notes + ----- + This method **does not** set the gravity acceleration, it only returns + a :class:`Function` object representing the gravity model. + + Examples + -------- + Let's prepare a `Environment` object with a constant gravity + acceleration: + + >>> g_0 = 9.80665 + >>> env_cte_g = Environment(gravity=g_0) + >>> env_cte_g.gravity([0, 100, 1000]) + [9.80665, 9.80665, 9.80665] + + It's also possible to variate the gravity acceleration by defining + its function of height: + + >>> R_t = 6371000 + >>> g_func = lambda h : g_0 * (R_t / (R_t + h))**2 + >>> env_var_g = Environment(gravity=g_func) + >>> g = env_var_g.gravity(1000) + >>> print(f"{g:.6f}") + 9.803572 """ if gravity is None: return self.somigliana_gravity.set_discrete( @@ -500,7 +606,7 @@ def max_expected_height(self, value): @funcify_method("height (m)", "gravity (m/s²)") def somigliana_gravity(self, height): - """Computes the gravity acceleration with the Somigliana formula. + """Computes the gravity acceleration with the Somigliana formula [1]_. An height correction is applied to the normal gravity that is accurate for heights used in aviation. The formula is based on the WGS84 ellipsoid, but is accurate for other reference ellipsoids. @@ -514,6 +620,10 @@ def somigliana_gravity(self, height): ------- Function Function object representing the gravity model. + + References + ---------- + .. [1] https://en.wikipedia.org/wiki/Theoretical_gravity#Somigliana_equation """ a = 6378137.0 # semi_major_axis f = 1 / 298.257223563 # flattening_factor @@ -571,18 +681,9 @@ def set_elevation(self, elevation="Open-Elevation"): # self.elevation = elev - elif self.latitude != None and self.longitude != None: - try: - print("Fetching elevation from open-elevation.com...") - request_url = "https://api.open-elevation.com/api/v1/lookup?locations={:f},{:f}".format( - self.latitude, self.longitude - ) - response = requests.get(request_url) - results = response.json()["results"] - self.elevation = results[0]["elevation"] - print("Elevation received:", self.elevation) - except: - raise RuntimeError("Unable to reach Open-Elevation API servers.") + elif self.latitude is not None and self.longitude is not None: + self.elevation = self.__fetch_open_elevation() + print("Elevation received: ", self.elevation) else: raise ValueError( "Latitude and longitude must be set to use" @@ -1194,26 +1295,8 @@ def set_atmospheric_model( "v_wind": "vgrdprs", } # Attempt to get latest forecast - time_attempt = datetime.utcnow() - success = False - attempt_count = 0 - while not success and attempt_count < 10: - time_attempt -= timedelta(hours=6 * attempt_count) - file = "https://nomads.ncep.noaa.gov/dods/gens_bc/gens{:04d}{:02d}{:02d}/gep_all_{:02d}z".format( - time_attempt.year, - time_attempt.month, - time_attempt.day, - 6 * (time_attempt.hour // 6), - ) - try: - self.process_ensemble(file, dictionary) - success = True - except OSError: - attempt_count += 1 - if not success: - raise RuntimeError( - "Unable to load latest weather data for GEFS through " + file - ) + self.__fetch_gefs_ensemble(dictionary) + elif file == "CMC": # Define dictionary dictionary = { @@ -1229,27 +1312,7 @@ def set_atmospheric_model( "u_wind": "ugrdprs", "v_wind": "vgrdprs", } - # Attempt to get latest forecast - time_attempt = datetime.utcnow() - success = False - attempt_count = 0 - while not success and attempt_count < 10: - time_attempt -= timedelta(hours=12 * attempt_count) - file = "https://nomads.ncep.noaa.gov/dods/cmcens/cmcens{:04d}{:02d}{:02d}/cmcens_all_{:02d}z".format( - time_attempt.year, - time_attempt.month, - time_attempt.day, - 12 * (time_attempt.hour // 12), - ) - try: - self.process_ensemble(file, dictionary) - success = True - except OSError: - attempt_count += 1 - if not success: - raise RuntimeError( - "Unable to load latest weather data for CMC through " + file - ) + self.__fetch_cmc_ensemble(dictionary) # Process other forecasts or reanalysis else: # Check if default dictionary was requested @@ -1488,24 +1551,26 @@ def process_custom_atmosphere( # Check maximum height of custom wind input if not callable(self.wind_velocity_x.source): max_expected_height = max(self.wind_velocity_x[-1, 0], max_expected_height) - if not callable(self.wind_velocity_y.source): - max_expected_height = max(self.wind_velocity_y[-1, 0], max_expected_height) - # Compute wind profile direction and heading - wind_heading = ( - lambda h: np.arctan2(self.wind_velocity_x(h), self.wind_velocity_y(h)) - * (180 / np.pi) - % 360 - ) + def wind_heading_func(h): + return ( + np.arctan2( + self.wind_velocity_x.get_value_opt(h), + self.wind_velocity_y.get_value_opt(h), + ) + * (180 / np.pi) + % 360 + ) + self.wind_heading = Function( - wind_heading, + wind_heading_func, inputs="Height Above Sea Level (m)", outputs="Wind Heading (Deg True)", interpolation="linear", ) def wind_direction(h): - return (wind_heading(h) - 180) % 360 + return (wind_heading_func(h) - 180) % 360 self.wind_direction = Function( wind_direction, @@ -1515,7 +1580,10 @@ def wind_direction(h): ) def wind_speed(h): - return np.sqrt(self.wind_velocity_x(h) ** 2 + self.wind_velocity_y(h) ** 2) + return np.sqrt( + self.wind_velocity_x.get_value_opt(h) ** 2 + + self.wind_velocity_y.get_value_opt(h) ** 2 + ) self.wind_speed = Function( wind_speed, @@ -1541,20 +1609,7 @@ def process_windy_atmosphere(self, model="ECMWF"): model. """ - # Process the model string - model = model.lower() - if model[-1] == "u": # case iconEu - model = "".join([model[:4], model[4].upper(), model[4 + 1 :]]) - # Load data from Windy.com: json file - url = f"https://node.windy.com/forecast/meteogram/{model}/{self.latitude}/{self.longitude}/?step=undefined" - try: - response = requests.get(url).json() - except: - if model == "iconEu": - raise ValueError( - "Could not get a valid response for Icon-EU from Windy. Check if the latitude and longitude coordinates set are inside Europe.", - ) - raise + response = self.__fetch_atmospheric_data_from_windy(model) # Determine time index from model time_array = np.array(response["data"]["hours"]) @@ -1715,18 +1770,7 @@ def process_wyoming_sounding(self, file): None """ # Request Wyoming Sounding from file url - response = requests.get(file) - if response.status_code != 200: - raise ImportError("Unable to load " + file + ".") - if len(re.findall("Can't get .+ Observations at", response.text)): - raise ValueError( - re.findall("Can't get .+ Observations at .+", response.text)[0] - + " Check station number and date." - ) - if response.text == "Invalid OUTPUT: specified\n": - raise ValueError( - "Invalid OUTPUT: specified. Make sure the output is Text: List." - ) + response = self.__fetch_wyoming_sounding(file) # Process Wyoming Sounding by finding data table and station info response_split_text = re.split("(<.{0,1}PRE>)", response.text) @@ -1852,9 +1896,7 @@ def process_noaaruc_sounding(self, file): None """ # Request NOAA Ruc Sounding from file url - response = requests.get(file) - if response.status_code != 200 or len(response.text) < 10: - raise ImportError("Unable to load " + file + ".") + response = self.__fetch_noaaruc_sounding(file) # Split response into lines lines = response.text.split("\n") @@ -3040,6 +3082,24 @@ def calculate_density_profile(self): Returns ------- None + + Examples + -------- + Creating an Environment object and calculating the density + at Sea Level: + + >>> env = Environment() + >>> env.calculate_density_profile() + >>> env.density(0) + 1.225000018124288 + + Creating an Environment object and calculating the density + at 1000m above Sea Level: + + >>> env = Environment() + >>> env.calculate_density_profile() + >>> env.density(1000) + 1.1116193933422585 """ # Retrieve pressure P, gas constant R and temperature T P = self.pressure @@ -3142,22 +3202,26 @@ def add_wind_gust(self, wind_gust_x, wind_gust_y): # Reset wind heading and velocity magnitude self.wind_heading = Function( lambda h: (180 / np.pi) - * np.arctan2(self.wind_velocity_x(h), self.wind_velocity_y(h)) + * np.arctan2( + self.wind_velocity_x.get_value_opt(h), + self.wind_velocity_y.get_value_opt(h), + ) % 360, "Height (m)", "Wind Heading (degrees)", extrapolation="constant", ) self.wind_speed = Function( - lambda h: (self.wind_velocity_x(h) ** 2 + self.wind_velocity_y(h) ** 2) + lambda h: ( + self.wind_velocity_x.get_value_opt(h) ** 2 + + self.wind_velocity_y.get_value_opt(h) ** 2 + ) ** 0.5, "Height (m)", "Wind Speed (m/s)", extrapolation="constant", ) - return None - def info(self): """Prints most important data and graphs available about the Environment. @@ -3429,6 +3493,110 @@ def set_earth_geometry(self, datum): f"The reference system {datum} for Earth geometry " "is not recognized." ) + # Auxiliary functions - Fetching Data from 3rd party APIs + + @exponential_backoff(max_attempts=3, base_delay=1, max_delay=60) + def __fetch_open_elevation(self): + print("Fetching elevation from open-elevation.com...") + request_url = ( + "https://api.open-elevation.com/api/v1/lookup?locations" + f"={self.latitude},{self.longitude}" + ) + try: + response = requests.get(request_url) + except Exception as e: + raise RuntimeError("Unable to reach Open-Elevation API servers.") + results = response.json()["results"] + return results[0]["elevation"] + + @exponential_backoff(max_attempts=5, base_delay=2, max_delay=60) + def __fetch_atmospheric_data_from_windy(self, model): + model = model.lower() + if model[-1] == "u": # case iconEu + model = "".join([model[:4], model[4].upper(), model[4 + 1 :]]) + url = ( + f"https://node.windy.com/forecast/meteogram/{model}/" + f"{self.latitude}/{self.longitude}/?step=undefined" + ) + try: + response = requests.get(url).json() + except Exception as e: + if model == "iconEu": + raise ValueError( + "Could not get a valid response for Icon-EU from Windy. " + "Check if the coordinates are set inside Europe." + ) + return response + + @exponential_backoff(max_attempts=5, base_delay=2, max_delay=60) + def __fetch_wyoming_sounding(self, file): + response = requests.get(file) + if response.status_code != 200: + raise ImportError(f"Unable to load {file}.") + if len(re.findall("Can't get .+ Observations at", response.text)): + raise ValueError( + re.findall("Can't get .+ Observations at .+", response.text)[0] + + " Check station number and date." + ) + if response.text == "Invalid OUTPUT: specified\n": + raise ValueError( + "Invalid OUTPUT: specified. Make sure the output is Text: List." + ) + return response + + @exponential_backoff(max_attempts=5, base_delay=2, max_delay=60) + def __fetch_noaaruc_sounding(self, file): + response = requests.get(file) + if response.status_code != 200 or len(response.text) < 10: + raise ImportError("Unable to load " + file + ".") + + @exponential_backoff(max_attempts=5, base_delay=2, max_delay=60) + def __fetch_gefs_ensemble(self, dictionary): + time_attempt = datetime.now(tz=timezone.utc) + success = False + attempt_count = 0 + while not success and attempt_count < 10: + time_attempt -= timedelta(hours=6 * attempt_count) + file = ( + f"https://nomads.ncep.noaa.gov/dods/gens_bc/gens" + f"{time_attempt.year:04d}{time_attempt.month:02d}" + f"{time_attempt.day:02d}/" + f"gep_all_{6 * (time_attempt.hour // 6):02d}z" + ) + try: + self.process_ensemble(file, dictionary) + success = True + except OSError: + attempt_count += 1 + if not success: + raise RuntimeError( + "Unable to load latest weather data for GEFS through " + file + ) + + @exponential_backoff(max_attempts=5, base_delay=2, max_delay=60) + def __fetch_cmc_ensemble(self, dictionary): + # Attempt to get latest forecast + time_attempt = datetime.now(tz=timezone.utc) + success = False + attempt_count = 0 + while not success and attempt_count < 10: + time_attempt -= timedelta(hours=12 * attempt_count) + file = ( + f"https://nomads.ncep.noaa.gov/dods/cmcens/" + f"cmcens{time_attempt.year:04d}{time_attempt.month:02d}" + f"{time_attempt.day:02d}/" + f"cmcens_all_{12 * (time_attempt.hour // 12):02d}z" + ) + try: + self.process_ensemble(file, dictionary) + success = True + except OSError: + attempt_count += 1 + if not success: + raise RuntimeError( + "Unable to load latest weather data for CMC through " + file + ) + # Auxiliary functions - Geodesic Coordinates @staticmethod diff --git a/rocketpy/environment/environment_analysis.py b/rocketpy/environment/environment_analysis.py index c15b32551..da6fde364 100644 --- a/rocketpy/environment/environment_analysis.py +++ b/rocketpy/environment/environment_analysis.py @@ -3,6 +3,7 @@ import datetime import json from collections import defaultdict +from functools import cached_property import netCDF4 import numpy as np @@ -22,11 +23,6 @@ from ..units import convert_units from .environment import Environment -try: - from functools import cached_property -except ImportError: - from ..tools import cached_property - # TODO: the average_wind_speed_profile_by_hour and similar methods could be more abstract than currently are @@ -441,7 +437,7 @@ def __localize_input_dates(self): def __find_preferred_timezone(self): if self.preferred_timezone is None: - # Use local timezone based on lat lon pair + # Use local time zone based on lat lon pair try: timezonefinder = import_optional_dependency("timezonefinder") tf = timezonefinder.TimezoneFinder() diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index cefed044d..6dc1c764b 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -5,8 +5,10 @@ """ import warnings +from bisect import bisect_left from collections.abc import Iterable from copy import deepcopy +from functools import cached_property from inspect import signature from pathlib import Path @@ -14,12 +16,15 @@ import numpy as np from scipy import integrate, linalg, optimize -try: - from functools import cached_property -except ImportError: - from ..tools import cached_property - NUMERICAL_TYPES = (float, int, complex, np.ndarray, np.integer, np.floating) +INTERPOLATION_TYPES = { + "linear": 0, + "polynomial": 1, + "akima": 2, + "spline": 3, + "shepard": 4, +} +EXTRAPOLATION_TYPES = {"zero": 0, "natural": 1, "constant": 2} class Function: @@ -110,28 +115,20 @@ def __init__( (II) Fields in CSV files may be enclosed in double quotes. If fields are not quoted, double quotes should not appear inside them. """ - if inputs is None: - inputs = ["Scalar"] - if outputs is None: - outputs = ["Scalar"] - - inputs, outputs, interpolation, extrapolation = self._check_user_input( - source, inputs, outputs, interpolation, extrapolation - ) - - # initialize variables to avoid errors when being called by other methods - self.get_value_opt = None - self.__polynomial_coefficients__ = None - self.__akima_coefficients__ = None - self.__spline_coefficients__ = None - - # store variables - self.set_inputs(inputs) - self.set_outputs(outputs) + # initialize parameters + self.source = source + self.__inputs__ = inputs + self.__outputs__ = outputs self.__interpolation__ = interpolation self.__extrapolation__ = extrapolation - self.set_source(source) - self.set_title(title) + self.title = title + self.__img_dim__ = 1 # always 1, here for backwards compatibility + + # args must be passed from self. + self.set_source(self.source) + self.set_inputs(self.__inputs__) + self.set_outputs(self.__outputs__) + self.set_title(self.title) # Define all set methods def set_inputs(self, inputs): @@ -146,8 +143,7 @@ def set_inputs(self, inputs): ------- self : Function """ - self.__inputs__ = [inputs] if isinstance(inputs, str) else list(inputs) - self.__dom_dim__ = len(self.__inputs__) + self.__inputs__ = self.__validate_inputs(inputs) return self def set_outputs(self, outputs): @@ -162,8 +158,7 @@ def set_outputs(self, outputs): ------- self : Function """ - self.__outputs__ = [outputs] if isinstance(outputs, str) else list(outputs) - self.__img_dim__ = len(self.__outputs__) + self.__outputs__ = self.__validate_outputs(outputs) return self def set_source(self, source): @@ -212,101 +207,45 @@ def set_source(self, source): self : Function Returns the Function instance. """ - *_, interpolation, extrapolation = self._check_user_input( - source, - self.__inputs__, - self.__outputs__, - self.__interpolation__, - self.__extrapolation__, - ) - # If the source is a Function - if isinstance(source, Function): - source = source.get_source() - # Import CSV if source is a string or Path and convert values to ndarray - if isinstance(source, (str, Path)): - with open(source, "r") as file: - try: - source = np.loadtxt(file, delimiter=",", dtype=float) - except ValueError: - # If an error occurs, headers are present - source = np.loadtxt(source, delimiter=",", dtype=float, skiprows=1) - except Exception as e: - raise ValueError( - "The source file is not a valid csv or txt file." - ) from e - - # Convert to ndarray if source is a list - if isinstance(source, (list, tuple)): - source = np.array(source, dtype=np.float64) - # Convert number source into vectorized lambda function - if isinstance(source, (int, float)): - temp = 1 * source - - def source_function(_): - return temp - - source = source_function + source = self.__validate_source(source) # Handle callable source or number source if callable(source): - # Set source - self.source = source - # Set get_value_opt self.get_value_opt = source + self.__interpolation__ = None + self.__extrapolation__ = None + # Set arguments name and domain dimensions parameters = signature(source).parameters self.__dom_dim__ = len(parameters) - if self.__inputs__ == ["Scalar"]: + if self.__inputs__ is None: self.__inputs__ = list(parameters) - # Set interpolation and extrapolation - self.__interpolation__ = None - self.__extrapolation__ = None + # Handle ndarray source else: - # Check to see if dimensions match incoming data set - new_total_dim = len(source[0, :]) - old_total_dim = self.__dom_dim__ + self.__img_dim__ - - # If they don't, update default values or throw error - if new_total_dim != old_total_dim: - # Update dimensions and inputs - self.__dom_dim__ = new_total_dim - 1 - self.__inputs__ = self.__dom_dim__ * self.__inputs__ + # Evaluate dimension + self.__dom_dim__ = source.shape[1] - 1 - # Do things if domDim is 1 + # set x and y. If Function is 2D, also set z if self.__dom_dim__ == 1: source = source[source[:, 0].argsort()] - self.x_array = source[:, 0] self.x_initial, self.x_final = self.x_array[0], self.x_array[-1] - self.y_array = source[:, 1] self.y_initial, self.y_final = self.y_array[0], self.y_array[-1] - - # Finally set data source as source - self.source = source - # Do things if function is multivariate - else: + self.get_value_opt = self.__get_value_opt_1d + elif self.__dom_dim__ > 1: self.x_array = source[:, 0] self.x_initial, self.x_final = self.x_array[0], self.x_array[-1] - self.y_array = source[:, 1] self.y_initial, self.y_final = self.y_array[0], self.y_array[-1] - self.z_array = source[:, 2] self.z_initial, self.z_final = self.z_array[0], self.z_array[-1] + self.get_value_opt = self.__get_value_opt_nd - # Finally set data source as source - self.source = source - # Update extrapolation method - if self.__extrapolation__ is None: - self.set_extrapolation(extrapolation) - # Set default interpolation for point source if it hasn't - if self.__interpolation__ is None: - self.set_interpolation(interpolation) - else: - # Updates interpolation coefficients - self.set_interpolation(self.__interpolation__) + self.source = source + self.set_interpolation(self.__interpolation__) + self.set_extrapolation(self.__extrapolation__) return self @cached_property @@ -346,21 +285,27 @@ def set_interpolation(self, method="spline"): ------- self : Function """ - # Set interpolation method - self.__interpolation__ = method + if not callable(self.source): + self.__interpolation__ = self.__validate_interpolation(method) + self.__update_interpolation_coefficients(self.__interpolation__) + self.__set_interpolation_func() + return self + + def __update_interpolation_coefficients(self, method): + """Update interpolation coefficients for the given method.""" # Spline, akima and polynomial need data processing # Shepard, and linear do not - if method == "spline": - self.__interpolate_spline__() - elif method == "polynomial": + if method == "polynomial": self.__interpolate_polynomial__() + self._coeffs = self.__polynomial_coefficients__ elif method == "akima": self.__interpolate_akima__() - - # Set get_value_opt - self.set_get_value_opt() - - return self + self._coeffs = self.__akima_coefficients__ + elif method == "spline" or method is None: + self.__interpolate_spline__() + self._coeffs = self.__spline_coefficients__ + else: + self._coeffs = [] def set_extrapolation(self, method="constant"): """Set extrapolation behavior of data set. @@ -379,137 +324,163 @@ def set_extrapolation(self, method="constant"): self : Function The Function object. """ - self.__extrapolation__ = method + if not callable(self.source): + self.__extrapolation__ = self.__validate_extrapolation(method) + self.__set_extrapolation_func() return self + def __set_interpolation_func(self): + """Defines interpolation function used by the Function. Each + interpolation method has its own function with exception of shepard, + which has its interpolation/extrapolation function defined in + ``Function.__interpolate_shepard__``. The function is stored in + the attribute _interpolation_func.""" + interpolation = INTERPOLATION_TYPES[self.__interpolation__] + if interpolation == 0: # linear + + def linear_interpolation(x, x_min, x_max, x_data, y_data, coeffs): + x_interval = bisect_left(x_data, x) + x_left = x_data[x_interval - 1] + y_left = y_data[x_interval - 1] + dx = float(x_data[x_interval] - x_left) + dy = float(y_data[x_interval] - y_left) + return (x - x_left) * (dy / dx) + y_left + + self._interpolation_func = linear_interpolation + + elif interpolation == 1: # polynomial + + def polynomial_interpolation(x, x_min, x_max, x_data, y_data, coeffs): + return np.sum(coeffs * x ** np.arange(len(coeffs))) + + self._interpolation_func = polynomial_interpolation + + elif interpolation == 2: # akima + + def akima_interpolation(x, x_min, x_max, x_data, y_data, coeffs): + x_interval = bisect_left(x_data, x) + x_interval = x_interval if x_interval != 0 else 1 + a = coeffs[4 * x_interval - 4 : 4 * x_interval] + return a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] + + self._interpolation_func = akima_interpolation + + elif interpolation == 3: # spline + + def spline_interpolation(x, x_min, x_max, x_data, y_data, coeffs): + x_interval = bisect_left(x_data, x) + x_interval = max(x_interval, 1) + a = coeffs[:, x_interval - 1] + x = x - x_data[x_interval - 1] + return a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] + + self._interpolation_func = spline_interpolation + + elif interpolation == 4: # shepard does not use interpolation function + self._interpolation_func = None + + def __set_extrapolation_func(self): + """Defines extrapolation function used by the Function. Each + extrapolation method has its own function. The function is stored in + the attribute _extrapolation_func.""" + interpolation = INTERPOLATION_TYPES[self.__interpolation__] + extrapolation = EXTRAPOLATION_TYPES[self.__extrapolation__] + + if interpolation == 4: # shepard does not use extrapolation function + self._extrapolation_func = None + + elif extrapolation == 0: # zero + + def zero_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): + return 0 + + self._extrapolation_func = zero_extrapolation + elif extrapolation == 1: # natural + if interpolation == 0: # linear + + def natural_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): + x_interval = 1 if x < x_min else -1 + x_left = x_data[x_interval - 1] + y_left = y_data[x_interval - 1] + dx = float(x_data[x_interval] - x_left) + dy = float(y_data[x_interval] - y_left) + return (x - x_left) * (dy / dx) + y_left + + elif interpolation == 1: # polynomial + + def natural_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): + return np.sum(coeffs * x ** np.arange(len(coeffs))) + + elif interpolation == 2: # akima + + def natural_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): + a = coeffs[:4] if x < x_min else coeffs[-4:] + return a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] + + elif interpolation == 3: # spline + + def natural_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): + if x < x_min: + a = coeffs[:, 0] + x = x - x_data[0] + else: + a = coeffs[:, -1] + x = x - x_data[-2] + return a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] + + self._extrapolation_func = natural_extrapolation + elif extrapolation == 2: # constant + + def constant_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): + return y_data[0] if x < x_min else y_data[-1] + + self._extrapolation_func = constant_extrapolation + def set_get_value_opt(self): - """Crates a method that evaluates interpolations rather quickly - when compared to other options available, such as just calling - the object instance or calling ``Function.get_value`` directly. See - ``Function.get_value_opt`` for documentation. + """Defines a method that evaluates interpolations. Returns ------- self : Function """ + if callable(self.source): + self.get_value_opt = self.source + elif self.__dom_dim__ == 1: + self.get_value_opt = self.__get_value_opt_1d + elif self.__dom_dim__ > 1: + self.get_value_opt = self.__get_value_opt_nd + return self + + def __get_value_opt_1d(self, x): + """Evaluate the Function at a single point x. This method is used + when the Function is 1-D. + + Parameters + ---------- + x : scalar + Value where the Function is to be evaluated. + + Returns + ------- + y : scalar + Value of the Function at the specified point. + """ # Retrieve general info x_data = self.x_array y_data = self.y_array x_min, x_max = self.x_initial, self.x_final - if self.__extrapolation__ == "zero": - extrapolation = 0 # Extrapolation is zero - elif self.__extrapolation__ == "natural": - extrapolation = 1 # Extrapolation is natural - elif self.__extrapolation__ == "constant": - extrapolation = 2 # Extrapolation is constant + coeffs = self._coeffs + if x_min <= x <= x_max: + y = self._interpolation_func(x, x_min, x_max, x_data, y_data, coeffs) else: - raise ValueError(f"Invalid extrapolation type {self.__extrapolation__}") + y = self._extrapolation_func(x, x_min, x_max, x_data, y_data, coeffs) + return y - # Crete method to interpolate this info for each interpolation type - if self.__interpolation__ == "spline": - coeffs = self.__spline_coefficients__ - - def get_value_opt(x): - x_interval = np.searchsorted(x_data, x) - # Interval found... interpolate... or extrapolate - if x_min <= x <= x_max: - # Interpolate - x_interval = x_interval if x_interval != 0 else 1 - a = coeffs[:, x_interval - 1] - x = x - x_data[x_interval - 1] - y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - else: - # Extrapolate - if extrapolation == 0: # Extrapolation == zero - y = 0 - elif extrapolation == 1: # Extrapolation == natural - a = coeffs[:, 0] if x < x_min else coeffs[:, -1] - x = x - x_data[0] if x < x_min else x - x_data[-2] - y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - else: # Extrapolation is set to constant - y = y_data[0] if x < x_min else y_data[-1] - return y - - elif self.__interpolation__ == "linear": - - def get_value_opt(x): - x_interval = np.searchsorted(x_data, x) - # Interval found... interpolate... or extrapolate - if x_min <= x <= x_max: - # Interpolate - dx = float(x_data[x_interval] - x_data[x_interval - 1]) - dy = float(y_data[x_interval] - y_data[x_interval - 1]) - y = (x - x_data[x_interval - 1]) * (dy / dx) + y_data[ - x_interval - 1 - ] - else: - # Extrapolate - if extrapolation == 0: # Extrapolation == zero - y = 0 - elif extrapolation == 1: # Extrapolation == natural - x_interval = 1 if x < x_min else -1 - dx = float(x_data[x_interval] - x_data[x_interval - 1]) - dy = float(y_data[x_interval] - y_data[x_interval - 1]) - y = (x - x_data[x_interval - 1]) * (dy / dx) + y_data[ - x_interval - 1 - ] - else: # Extrapolation is set to constant - y = y_data[0] if x < x_min else y_data[-1] - return y - - elif self.__interpolation__ == "akima": - coeffs = np.array(self.__akima_coefficients__) - - def get_value_opt(x): - x_interval = np.searchsorted(x_data, x) - # Interval found... interpolate... or extrapolate - if x_min <= x <= x_max: - # Interpolate - x_interval = x_interval if x_interval != 0 else 1 - a = coeffs[4 * x_interval - 4 : 4 * x_interval] - y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - else: - # Extrapolate - if extrapolation == 0: # Extrapolation == zero - y = 0 - elif extrapolation == 1: # Extrapolation == natural - a = coeffs[:4] if x < x_min else coeffs[-4:] - y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - else: # Extrapolation is set to constant - y = y_data[0] if x < x_min else y_data[-1] - return y - - elif self.__interpolation__ == "polynomial": - coeffs = self.__polynomial_coefficients__ - - def get_value_opt(x): - # Interpolate... or extrapolate - if x_min <= x <= x_max: - # Interpolate - y = 0 - for i, coef in enumerate(coeffs): - y += coef * (x**i) - else: - # Extrapolate - if extrapolation == 0: # Extrapolation == zero - y = 0 - elif extrapolation == 1: # Extrapolation == natural - y = 0 - for i, coef in enumerate(coeffs): - y += coef * (x**i) - else: # Extrapolation is set to constant - y = y_data[0] if x < x_min else y_data[-1] - return y - - elif self.__interpolation__ == "shepard": - # change the function's name to avoid mypy's error - def get_value_opt_multiple(*args): - return self.__interpolate_shepard__(args) - - get_value_opt = get_value_opt_multiple - - self.get_value_opt = get_value_opt - return self + def __get_value_opt_nd(self, *args): + """Evaluate the Function at a single point (x, y, z). This method is + used when the Function is N-D.""" + # always use shepard for N-D functions + return self.__interpolate_shepard__(args) def set_discrete( self, @@ -904,107 +875,20 @@ def get_value(self, *args): if all(isinstance(arg, Iterable) for arg in args): return [self.source(*arg) for arg in zip(*args)] - # Returns value for shepard interpolation - elif self.__interpolation__ == "shepard": - return self.__interpolate_shepard__(args) + elif self.__dom_dim__ > 1: # deals with nd functions and shepard interp + return self.get_value_opt(*args) - # Returns value for polynomial interpolation function type - elif self.__interpolation__ == "polynomial": - if isinstance(args[0], (int, float)): - args = [list(args)] - x = np.array(args[0]) - x_data = self.x_array - y_data = self.y_array - x_min, x_max = self.x_initial, self.x_final - coeffs = self.__polynomial_coefficients__ - matrix = np.zeros((len(args[0]), coeffs.shape[0])) - for i in range(coeffs.shape[0]): - matrix[:, i] = x**i - ans = matrix.dot(coeffs).tolist() - for i, _ in enumerate(x): - if not x_min <= x[i] <= x_max: - if self.__extrapolation__ == "constant": - ans[i] = y_data[0] if x[i] < x_min else y_data[-1] - elif self.__extrapolation__ == "zero": - ans[i] = 0 - return ans if len(ans) > 1 else ans[0] - # Returns value for spline, akima or linear interpolation function type - elif self.__interpolation__ in ["spline", "akima", "linear"]: + # Returns value for other interpolation type + else: # interpolation is "polynomial", "spline", "akima" or "linear" if isinstance(args[0], (int, float, complex, np.integer)): args = [list(args)] - x = list(args[0]) - x_data = self.x_array - y_data = self.y_array - x_intervals = np.searchsorted(x_data, x) - x_min, x_max = self.x_initial, self.x_final - if self.__interpolation__ == "spline": - coeffs = self.__spline_coefficients__ - for i, _ in enumerate(x): - if x[i] == x_min or x[i] == x_max: - x[i] = y_data[x_intervals[i]] - elif x_min < x[i] < x_max or (self.__extrapolation__ == "natural"): - if not x_min < x[i] < x_max: - a = coeffs[:, 0] if x[i] < x_min else coeffs[:, -1] - x[i] = ( - x[i] - x_data[0] if x[i] < x_min else x[i] - x_data[-2] - ) - else: - a = coeffs[:, x_intervals[i] - 1] - x[i] = x[i] - x_data[x_intervals[i] - 1] - x[i] = a[3] * x[i] ** 3 + a[2] * x[i] ** 2 + a[1] * x[i] + a[0] - else: - # Extrapolate - if self.__extrapolation__ == "zero": - x[i] = 0 - else: # Extrapolation is set to constant - x[i] = y_data[0] if x[i] < x_min else y_data[-1] - elif self.__interpolation__ == "linear": - for i, _ in enumerate(x): - # Interval found... interpolate... or extrapolate - inter = x_intervals[i] - if x_min <= x[i] <= x_max: - # Interpolate - dx = float(x_data[inter] - x_data[inter - 1]) - dy = float(y_data[inter] - y_data[inter - 1]) - x[i] = (x[i] - x_data[inter - 1]) * (dy / dx) + y_data[ - inter - 1 - ] - else: - # Extrapolate - if self.__extrapolation__ == "zero": # Extrapolation == zero - x[i] = 0 - elif ( - self.__extrapolation__ == "natural" - ): # Extrapolation == natural - inter = 1 if x[i] < x_min else -1 - dx = float(x_data[inter] - x_data[inter - 1]) - dy = float(y_data[inter] - y_data[inter - 1]) - x[i] = (x[i] - x_data[inter - 1]) * (dy / dx) + y_data[ - inter - 1 - ] - else: # Extrapolation is set to constant - x[i] = y_data[0] if x[i] < x_min else y_data[-1] - else: - coeffs = self.__akima_coefficients__ - for i, _ in enumerate(x): - if x[i] == x_min or x[i] == x_max: - x[i] = y_data[x_intervals[i]] - elif x_min < x[i] < x_max or (self.__extrapolation__ == "natural"): - if not x_min < x[i] < x_max: - a = coeffs[:4] if x[i] < x_min else coeffs[-4:] - else: - a = coeffs[4 * x_intervals[i] - 4 : 4 * x_intervals[i]] - x[i] = a[3] * x[i] ** 3 + a[2] * x[i] ** 2 + a[1] * x[i] + a[0] - else: - # Extrapolate - if self.__extrapolation__ == "zero": - x[i] = 0 - else: # Extrapolation is set to constant - x[i] = y_data[0] if x[i] < x_min else y_data[-1] - if isinstance(args[0], np.ndarray): - return np.array(x) - else: - return x if len(x) > 1 else x[0] + + x = list(args[0]) + x = list(map(self.get_value_opt, x)) + if isinstance(args[0], np.ndarray): + return np.array(x) + else: + return x if len(x) > 1 else x[0] def __getitem__(self, args): """Returns item of the Function source. If the source is not an array, @@ -1988,7 +1872,7 @@ def __add__(self, other): # Create new Function object return Function(source, inputs, outputs, interpolation, extrapolation) else: - return Function(lambda x: (self.get_value(x) + other(x))) + return Function(lambda x: (self.get_value_opt(x) + other(x))) # If other is Float except... except AttributeError: if isinstance(other, NUMERICAL_TYPES): @@ -2009,10 +1893,10 @@ def __add__(self, other): source, inputs, outputs, interpolation, extrapolation ) else: - return Function(lambda x: (self.get_value(x) + other)) + return Function(lambda x: (self.get_value_opt(x) + other)) # Or if it is just a callable elif callable(other): - return Function(lambda x: (self.get_value(x) + other(x))) + return Function(lambda x: (self.get_value_opt(x) + other(x))) def __radd__(self, other): """Sums 'other' and a Function object and returns a new Function @@ -2055,7 +1939,7 @@ def __sub__(self, other): try: return self + (-other) except TypeError: - return Function(lambda x: (self.get_value(x) - other(x))) + return Function(lambda x: (self.get_value_opt(x) - other(x))) def __rsub__(self, other): """Subtracts a Function object from 'other' and returns a new Function @@ -2095,54 +1979,40 @@ def __mul__(self, other): result : Function A Function object which gives the result of self(x)*other(x). """ - # If other is Function try... - try: - # Check if Function objects source is array or callable - # Check if Function objects have the same domain discretization - if ( - isinstance(other.source, np.ndarray) - and isinstance(self.source, np.ndarray) - and self.__dom_dim__ == other.__dom_dim__ - and np.array_equal(self.x_array, other.x_array) - ): - # Operate on grid values - ys = self.y_array * other.y_array - xs = self.x_array - source = np.concatenate(([xs], [ys])).transpose() - # Retrieve inputs, outputs and interpolation - inputs = self.__inputs__[:] - outputs = self.__outputs__[0] + "*" + other.__outputs__[0] - outputs = "(" + outputs + ")" - interpolation = self.__interpolation__ - extrapolation = self.__extrapolation__ - # Create new Function object - return Function(source, inputs, outputs, interpolation, extrapolation) - else: - return Function(lambda x: (self.get_value(x) * other(x))) - # If other is Float except... - except AttributeError: - if isinstance(other, NUMERICAL_TYPES): - # Check if Function object source is array or callable - if isinstance(self.source, np.ndarray): - # Operate on grid values - ys = self.y_array * other - xs = self.x_array - source = np.concatenate(([xs], [ys])).transpose() - # Retrieve inputs, outputs and interpolation - inputs = self.__inputs__[:] - outputs = self.__outputs__[0] + "*" + str(other) - outputs = "(" + outputs + ")" - interpolation = self.__interpolation__ - extrapolation = self.__extrapolation__ - # Create new Function object - return Function( - source, inputs, outputs, interpolation, extrapolation - ) - else: - return Function(lambda x: (self.get_value(x) * other)) - # Or if it is just a callable - elif callable(other): - return Function(lambda x: (self.get_value(x) * other(x))) + self_source_is_array = isinstance(self.source, np.ndarray) + other_source_is_array = ( + isinstance(other.source, np.ndarray) + if isinstance(other, Function) + else False + ) + inputs = self.__inputs__[:] + interp = self.__interpolation__ + extrap = self.__extrapolation__ + + if ( + self_source_is_array + and other_source_is_array + and np.array_equal(self.x_array, other.x_array) + ): + source = np.column_stack((self.x_array, self.y_array * other.y_array)) + outputs = f"({self.__outputs__[0]}*{other.__outputs__[0]})" + return Function(source, inputs, outputs, interp, extrap) + elif isinstance(other, NUMERICAL_TYPES): + if not self_source_is_array: + return Function(lambda x: (self.get_value_opt(x) * other), inputs) + source = np.column_stack((self.x_array, np.multiply(self.y_array, other))) + outputs = f"({self.__outputs__[0]}*{other})" + return Function( + source, + inputs, + outputs, + interp, + extrap, + ) + elif callable(other): + return Function(lambda x: (self.get_value_opt(x) * other(x)), inputs) + else: + raise TypeError("Unsupported type for multiplication") def __rmul__(self, other): """Multiplies 'other' by a Function object and returns a new Function @@ -2335,10 +2205,10 @@ def __pow__(self, other): source, inputs, outputs, interpolation, extrapolation ) else: - return Function(lambda x: (self.get_value(x) ** other)) + return Function(lambda x: (self.get_value_opt(x) ** other)) # Or if it is just a callable elif callable(other): - return Function(lambda x: (self.get_value(x) ** other(x))) + return Function(lambda x: (self.get_value_opt(x) ** other(x))) def __rpow__(self, other): """Raises 'other' to the power of a Function object and returns @@ -2371,10 +2241,10 @@ def __rpow__(self, other): # Create new Function object return Function(source, inputs, outputs, interpolation, extrapolation) else: - return Function(lambda x: (other ** self.get_value(x))) + return Function(lambda x: (other ** self.get_value_opt(x))) # Or if it is just a callable elif callable(other): - return Function(lambda x: (other(x) ** self.get_value(x))) + return Function(lambda x: (other(x) ** self.get_value_opt(x))) def __matmul__(self, other): """Operator @ as an alias for composition. Therefore, this @@ -2549,10 +2419,12 @@ def differentiate(self, x, dx=1e-6, order=1): Evaluated derivative. """ if order == 1: - return (self.get_value(x + dx) - self.get_value(x - dx)) / (2 * dx) + return (self.get_value_opt(x + dx) - self.get_value_opt(x - dx)) / (2 * dx) elif order == 2: return ( - self.get_value(x + dx) - 2 * self.get_value(x) + self.get_value(x - dx) + self.get_value_opt(x + dx) + - 2 * self.get_value_opt(x) + + self.get_value_opt(x - dx) ) / dx**2 def identity_function(self): @@ -2996,189 +2868,188 @@ def savetxt( file.write(header_line + newline) np.savetxt(file, data_points, fmt=fmt, delimiter=delimiter, newline=newline) - @staticmethod - def _check_user_input( - source, - inputs=None, - outputs=None, - interpolation=None, - extrapolation=None, - ): - """ - Validates and processes the user input parameters for creating or - modifying a Function object. This function ensures the inputs, outputs, - interpolation, and extrapolation parameters are compatible with the - given source. It converts the source to a numpy array if necessary, sets - default values and raises warnings or errors for incompatible or - ill-defined parameters. + # Input validators + def __validate_source(self, source): + """Used to validate the source parameter for creating a Function object. Parameters ---------- - source : list, np.ndarray, or callable - The source data or Function object. If a list or ndarray, it should - contain numeric data. If a Function, its inputs and outputs are - checked against the provided inputs and outputs. - inputs : list of str or None - The names of the input variables. If None, defaults are generated - based on the dimensionality of the source. - outputs : str or list of str - The name(s) of the output variable(s). If a list is provided, it - must have a single element. - interpolation : str or None - The method of interpolation to be used. For multidimensional sources - it defaults to 'shepard' if not provided. - extrapolation : str or None - The method of extrapolation to be used. For multidimensional sources - it defaults to 'natural' if not provided. + source : np.ndarray, callable, str, Path, Function, list + The source data of the Function object. This can be a numpy array, + a callable function, a string or Path object to a csv or txt file, + a Function object, or a list of numbers. Returns ------- - tuple - A tuple containing the processed inputs, outputs, interpolation, and - extrapolation parameters. + np.ndarray, callable + The validated source parameter. Raises ------ ValueError - If the dimensionality of the source does not match the combined - dimensions of inputs and outputs. If the outputs list has more than - one element. - TypeError - If the source is not a list, np.ndarray, Function object, str or - Path. - Warning - If inputs or outputs do not match for a Function source, or if - defaults are used for inputs, interpolation,and extrapolation for a - multidimensional source. + If the source is not a valid type or if the source is not a 2D array + or a callable function. + """ + if isinstance(source, Function): + return source.get_source() - Examples - -------- - >>> from rocketpy import Function - >>> source = np.array([(1, 1), (2, 4), (3, 9)]) - >>> inputs = "x" - >>> outputs = ["y"] - >>> interpolation = 'linear' - >>> extrapolation = 'zero' - >>> inputs, outputs, interpolation, extrapolation = Function._check_user_input( - ... source, inputs, outputs, interpolation, extrapolation - ... ) - >>> inputs - ['x'] - >>> outputs - ['y'] - >>> interpolation - 'linear' - >>> extrapolation - 'zero' - """ - # check output type and dimensions - if isinstance(outputs, str): - outputs = [outputs] - if isinstance(inputs, str): - inputs = [inputs] + if isinstance(source, (str, Path)): + # Read csv or txt files and create a numpy array + try: + source = np.loadtxt(source, delimiter=",", dtype=np.float64) + except ValueError: + with open(source, "r") as file: + header, *data = file.read().splitlines() + + header = [label.strip("'").strip('"') for label in header.split(",")] + source = np.loadtxt(data, delimiter=",", dtype=np.float64) + + if len(source[0]) == len(header): + if self.__inputs__ is None: + self.__inputs__ = header[:-1] + if self.__outputs__ is None: + self.__outputs__ = [header[-1]] + except Exception as e: + raise ValueError( + "Could not read the csv or txt file to create Function source." + ) from e + + if isinstance(source, list) or isinstance(source, np.ndarray): + # Triggers an error if source is not a list of numbers + source = np.array(source, dtype=np.float64) - if len(outputs) > 1: + # Checks if 2D array + if len(source.shape) != 2: + raise ValueError( + "Source must be a 2D array in the form [[x1, x2 ..., xn, y], ...]." + ) + return source + + if isinstance(source, (int, float)): + # Convert number source into vectorized lambda function + temp = 1 * source + + def source_function(_): + return temp + + return source_function + + # If source is a callable function + return source + + def __validate_inputs(self, inputs): + """Used to validate the inputs parameter for creating a Function object. + It sets a default value if it is not provided. + + Parameters + ---------- + inputs : list of str, None + The name(s) of the input variable(s). If None, defaults to "Scalar". + + Returns + ------- + list + The validated inputs parameter. + """ + if self.__dom_dim__ == 1: + if inputs is None: + return ["Scalar"] + if isinstance(inputs, str): + return [inputs] + if isinstance(inputs, (list, tuple)): + if len(inputs) == 1: + return inputs raise ValueError( - "Output must either be a string or have dimension 1, " - + f"it currently has dimension ({len(outputs)})." + "Inputs must be a string or a list of strings with " + "the length of the domain dimension." + ) + if self.__dom_dim__ > 1: + if inputs is None: + return [f"Input {i+1}" for i in range(self.__dom_dim__)] + if isinstance(inputs, list): + if len(inputs) == self.__dom_dim__ and all( + isinstance(i, str) for i in inputs + ): + return inputs + raise ValueError( + "Inputs must be a list of strings with " + "the length of the domain dimension." ) - # check source for data type - # if list or ndarray, check for dimensions, interpolation and extrapolation - if isinstance(source, Function): - source = source.get_source() - if isinstance(source, (list, np.ndarray, str, Path)): - # Deal with csv or txt - if isinstance(source, (str, Path)): - # Convert to numpy array - try: - source = np.loadtxt(source, delimiter=",", dtype=float) - except ValueError: - with open(source, "r") as file: - header, *data = file.read().splitlines() - - header = [ - label.strip("'").strip('"') for label in header.split(",") - ] - source = np.loadtxt(data, delimiter=",", dtype=float) - - if len(source[0]) == len(header): - if inputs == ["Scalar"]: - inputs = header[:-1] - if outputs == ["Scalar"]: - outputs = [header[-1]] - except Exception as e: - raise ValueError( - "The source file is not a valid csv or txt file." - ) from e - - else: - # this will also trigger an error if the source is not a list of - # numbers or if the array is not homogeneous - source = np.array(source, dtype=np.float64) - - # check dimensions - source_dim = source.shape[1] - - # check interpolation and extrapolation - - ## single dimension - if source_dim == 2: - # possible interpolation values: linear, polynomial, akima and spline - if interpolation is None: - interpolation = "spline" - elif interpolation.lower() not in [ - "spline", - "linear", - "polynomial", - "akima", - ]: - warnings.warn( - "Interpolation method for single dimensional functions was " - + f"set to 'spline', the {interpolation} method is not supported." - ) - interpolation = "spline" - - # possible extrapolation values: constant, natural, zero - if extrapolation is None: - extrapolation = "constant" - elif extrapolation.lower() not in ["constant", "natural", "zero"]: - warnings.warn( - "Extrapolation method for single dimensional functions was " - + f"set to 'constant', the {extrapolation} method is not supported." - ) - extrapolation = "constant" - - ## multiple dimensions - elif source_dim > 2: - # check for inputs and outputs - if inputs == ["Scalar"]: - inputs = [f"Input {i+1}" for i in range(source_dim - 1)] - - if interpolation not in [None, "shepard"]: - warnings.warn( - ( - "Interpolation method for multidimensional functions was" - "set to 'shepard', other methods are not supported yet." - ), - ) - interpolation = "shepard" + def __validate_outputs(self, outputs): + """Used to validate the outputs parameter for creating a Function object. + It sets a default value if it is not provided. - if extrapolation not in [None, "natural"]: - warnings.warn( - "Extrapolation method for multidimensional functions was set" - "to 'natural', other methods are not supported yet." - ) - extrapolation = "natural" + Parameters + ---------- + outputs : str, list of str, None + The name of the output variables. If None, defaults to "Scalar". - # check input dimensions - in_out_dim = len(inputs) + len(outputs) - if source_dim != in_out_dim: + Returns + ------- + list + The validated outputs parameter. + """ + if outputs is None: + return ["Scalar"] + if isinstance(outputs, str): + return [outputs] + if isinstance(outputs, (list, tuple)): + if len(outputs) > 1: raise ValueError( - f"Source dimension ({source_dim}) does not match input " - + f"and output dimension ({in_out_dim})." + "Output must either be a string or a list of strings with " + + f"one item. It currently has dimension ({len(outputs)})." + ) + return outputs + + def __validate_interpolation(self, interpolation): + if self.__dom_dim__ == 1: + # possible interpolation values: linear, polynomial, akima and spline + if interpolation is None: + interpolation = "spline" + elif interpolation.lower() not in [ + "spline", + "linear", + "polynomial", + "akima", + ]: + warnings.warn( + "Interpolation method set to 'spline' because the " + f"{interpolation} method is not supported." + ) + interpolation = "spline" + ## multiple dimensions + elif self.__dom_dim__ > 1: + if interpolation not in [None, "shepard"]: + warnings.warn( + ( + "Interpolation method set to 'shepard'. Only 'shepard' " + "interpolation is supported for multiple dimensions." + ), + ) + interpolation = "shepard" + return interpolation + + def __validate_extrapolation(self, extrapolation): + if self.__dom_dim__ == 1: + if extrapolation is None: + extrapolation = "constant" + elif extrapolation.lower() not in ["constant", "natural", "zero"]: + warnings.warn( + "Extrapolation method set to 'constant' because the " + f"{extrapolation} method is not supported." + ) + extrapolation = "constant" + + ## multiple dimensions + elif self.__dom_dim__ > 1: + if extrapolation not in [None, "natural"]: + warnings.warn( + "Extrapolation method set to 'natural'. Other methods " + "are not supported yet." ) - return inputs, outputs, interpolation, extrapolation + extrapolation = "natural" + return extrapolation class PiecewiseFunction(Function): @@ -3262,7 +3133,7 @@ def calc_output(func, inputs): """ output = np.zeros(len(inputs)) for j, value in enumerate(inputs): - output[j] = func.get_value(value) + output[j] = func.get_value_opt(value) return output input_data = [] diff --git a/rocketpy/mathutils/vector_matrix.py b/rocketpy/mathutils/vector_matrix.py index 50d827659..332e1b680 100644 --- a/rocketpy/mathutils/vector_matrix.py +++ b/rocketpy/mathutils/vector_matrix.py @@ -1,8 +1,7 @@ from cmath import isclose +from functools import cached_property from itertools import product -from rocketpy.tools import cached_property - class Vector: """Pure python basic R3 vector class designed for simple operations. diff --git a/rocketpy/motors/hybrid_motor.py b/rocketpy/motors/hybrid_motor.py index 4b28a96d2..557333fe7 100644 --- a/rocketpy/motors/hybrid_motor.py +++ b/rocketpy/motors/hybrid_motor.py @@ -1,3 +1,7 @@ +from functools import cached_property + +from rocketpy.tools import parallel_axis_theorem_from_com + from ..mathutils.function import Function, funcify_method, reset_funcified_methods from ..plots.hybrid_motor_plots import _HybridMotorPlots from ..prints.hybrid_motor_prints import _HybridMotorPrints @@ -5,11 +9,6 @@ from .motor import Motor from .solid_motor import SolidMotor -try: - from functools import cached_property -except ImportError: - from ..tools import cached_property - class HybridMotor(Motor): """Class to specify characteristics and useful operations for Hybrid @@ -455,22 +454,21 @@ def propellant_I_11(self): ---------- .. [1] https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ - solid_correction = ( - self.solid.propellant_mass - * (self.solid.center_of_propellant_mass - self.center_of_propellant_mass) - ** 2 - ) - liquid_correction = ( - self.liquid.propellant_mass - * (self.liquid.center_of_propellant_mass - self.center_of_propellant_mass) - ** 2 - ) - I_11 = ( - self.solid.propellant_I_11 - + solid_correction - + self.liquid.propellant_I_11 - + liquid_correction + solid_mass = self.solid.propellant_mass + liquid_mass = self.liquid.propellant_mass + + cm = self.center_of_propellant_mass + solid_cm_to_cm = self.solid.center_of_propellant_mass - cm + liquid_cm_to_cm = self.liquid.center_of_propellant_mass - cm + + solid_prop_inertia = self.solid.propellant_I_11 + liquid_prop_inertia = self.liquid.propellant_I_11 + + I_11 = parallel_axis_theorem_from_com( + solid_prop_inertia, solid_mass, solid_cm_to_cm + ) + parallel_axis_theorem_from_com( + liquid_prop_inertia, liquid_mass, liquid_cm_to_cm ) return I_11 diff --git a/rocketpy/motors/liquid_motor.py b/rocketpy/motors/liquid_motor.py index ac525b379..7314e11ba 100644 --- a/rocketpy/motors/liquid_motor.py +++ b/rocketpy/motors/liquid_motor.py @@ -1,4 +1,4 @@ -import warnings +from functools import cached_property import numpy as np @@ -7,16 +7,12 @@ funcify_method, reset_funcified_methods, ) +from rocketpy.tools import parallel_axis_theorem_from_com from ..plots.liquid_motor_plots import _LiquidMotorPlots from ..prints.liquid_motor_prints import _LiquidMotorPrints from .motor import Motor -try: - from functools import cached_property -except ImportError: - from ..tools import cached_property - class LiquidMotor(Motor): """Class to specify characteristics and useful operations for Liquid @@ -388,10 +384,9 @@ def propellant_I_11(self): for positioned_tank in self.positioned_tanks: tank = positioned_tank.get("tank") tank_position = positioned_tank.get("position") - I_11 += ( - tank.inertia - + tank.fluid_mass - * (tank_position + tank.center_of_mass - center_of_mass) ** 2 + distance = tank_position + tank.center_of_mass - center_of_mass + I_11 += parallel_axis_theorem_from_com( + tank.inertia, tank.fluid_mass, distance ) return I_11 diff --git a/rocketpy/motors/motor.py b/rocketpy/motors/motor.py index 6c2242a9f..9429da88e 100644 --- a/rocketpy/motors/motor.py +++ b/rocketpy/motors/motor.py @@ -1,18 +1,14 @@ import re import warnings from abc import ABC, abstractmethod +from functools import cached_property import numpy as np from ..mathutils.function import Function, funcify_method from ..plots.motor_plots import _MotorPlots from ..prints.motor_prints import _MotorPrints -from ..tools import tuple_handler - -try: - from functools import cached_property -except ImportError: - from ..tools import cached_property +from ..tools import parallel_axis_theorem_from_com, tuple_handler class Motor(ABC): @@ -513,25 +509,19 @@ def I_11(self): ---------- .. [1] https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ - # Propellant inertia tensor 11 component wrt propellant center of mass - propellant_I_11 = self.propellant_I_11 - # Dry inertia tensor 11 component wrt dry center of mass + prop_I_11 = self.propellant_I_11 dry_I_11 = self.dry_I_11 - # Steiner theorem the get inertia wrt motor center of mass - propellant_I_11 += ( - self.propellant_mass - * (self.center_of_propellant_mass - self.center_of_mass) ** 2 - ) + prop_to_cm = self.center_of_propellant_mass - self.center_of_mass + dry_to_cm = self.center_of_dry_mass_position - self.center_of_mass - dry_I_11 += ( - self.dry_mass - * (self.center_of_dry_mass_position - self.center_of_mass) ** 2 + prop_I_11 = parallel_axis_theorem_from_com( + prop_I_11, self.propellant_mass, prop_to_cm ) + dry_I_11 = parallel_axis_theorem_from_com(dry_I_11, self.dry_mass, dry_to_cm) - # Sum of inertia components - return propellant_I_11 + dry_I_11 + return prop_I_11 + dry_I_11 @funcify_method("Time (s)", "Inertia I_22 (kg m²)") def I_22(self): diff --git a/rocketpy/motors/solid_motor.py b/rocketpy/motors/solid_motor.py index 8b1c2362e..db3527a95 100644 --- a/rocketpy/motors/solid_motor.py +++ b/rocketpy/motors/solid_motor.py @@ -1,3 +1,5 @@ +from functools import cached_property + import numpy as np from scipy import integrate @@ -6,11 +8,6 @@ from ..prints.solid_motor_prints import _SolidMotorPrints from .motor import Motor -try: - from functools import cached_property -except ImportError: - from ..tools import cached_property - class SolidMotor(Motor): """Class to specify characteristics and useful operations for solid motors. diff --git a/rocketpy/motors/tank_geometry.py b/rocketpy/motors/tank_geometry.py index f1940cbea..2eb7bd27e 100644 --- a/rocketpy/motors/tank_geometry.py +++ b/rocketpy/motors/tank_geometry.py @@ -11,10 +11,7 @@ cache = lru_cache(maxsize=None) -try: - from functools import cached_property -except ImportError: - from ..tools import cached_property +from functools import cached_property class TankGeometry: diff --git a/rocketpy/plots/environment_analysis_plots.py b/rocketpy/plots/environment_analysis_plots.py index 9fbf68d95..26727aba9 100644 --- a/rocketpy/plots/environment_analysis_plots.py +++ b/rocketpy/plots/environment_analysis_plots.py @@ -1454,7 +1454,7 @@ def wind_speed_profile_grid(self, clear_range_limits=False): if self.env_analysis.forecast: forecast = self.env_analysis.forecast y = self.env_analysis.average_wind_speed_profile_by_hour[hour][1] - x = forecast[hour].wind_speed.get_value(y) * convert_units( + x = forecast[hour].wind_speed.get_value_opt(y) * convert_units( 1, "m/s", self.env_analysis.unit_system["wind_speed"] ) ax.plot(x, y, "b--") diff --git a/rocketpy/plots/flight_plots.py b/rocketpy/plots/flight_plots.py index 4ae5141c9..21266a1f3 100644 --- a/rocketpy/plots/flight_plots.py +++ b/rocketpy/plots/flight_plots.py @@ -1,11 +1,8 @@ +from functools import cached_property + import matplotlib.pyplot as plt import numpy as np -try: - from functools import cached_property -except ImportError: - from ..tools import cached_property - class _FlightPlots: """Class that holds plot methods for Flight class. diff --git a/rocketpy/plots/rocket_plots.py b/rocketpy/plots/rocket_plots.py index 7074704e6..2be7a4a73 100644 --- a/rocketpy/plots/rocket_plots.py +++ b/rocketpy/plots/rocket_plots.py @@ -189,7 +189,6 @@ def draw(self, vis_args=None): A full list of color names can be found at: https://matplotlib.org/stable/gallery/color/named_colors """ - # TODO: we need to modularize this function, it is too big if vis_args is None: vis_args = { "background": "#EEEEEE", @@ -212,10 +211,27 @@ def draw(self, vis_args=None): reverse = csys == 1 self.rocket.aerodynamic_surfaces.sort_by_position(reverse=reverse) + drawn_surfaces = self._draw_aerodynamic_surfaces(ax, vis_args) + 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) + self._draw_center_of_mass_and_pressure(ax) + + plt.title("Rocket Representation") + plt.xlim() + plt.ylim([-self.rocket.radius * 4, self.rocket.radius * 6]) + plt.xlabel("Position (m)") + plt.ylabel("Radius (m)") + plt.legend(bbox_to_anchor=(1.05, 1), loc="upper left") + plt.tight_layout() + plt.show() + + def _draw_aerodynamic_surfaces(self, ax, vis_args): + """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 # and the radius of the rocket at that point drawn_surfaces = [] - # Idea is to get the shape of each aerodynamic surface in their own # coordinate system and then plot them in the rocket coordinate system # using the position of each surface @@ -225,113 +241,105 @@ def draw(self, vis_args=None): for surface, position in self.rocket.aerodynamic_surfaces: if isinstance(surface, NoseCone): - x_nosecone = -csys * surface.shape_vec[0] + position - y_nosecone = surface.shape_vec[1] - - ax.plot( - x_nosecone, - y_nosecone, - color=vis_args["nose"], - linewidth=vis_args["line_width"], - ) - ax.plot( - x_nosecone, - -y_nosecone, - color=vis_args["nose"], - linewidth=vis_args["line_width"], - ) - # close the nosecone - ax.plot( - [x_nosecone[-1], x_nosecone[-1]], - [y_nosecone[-1], -y_nosecone[-1]], - color=vis_args["nose"], - linewidth=vis_args["line_width"], - ) - - # Add the nosecone to the list of drawn surfaces - drawn_surfaces.append( - (surface, x_nosecone[-1], surface.rocket_radius, x_nosecone[-1]) - ) - + self._draw_nose_cone(ax, surface, position, drawn_surfaces, vis_args) elif isinstance(surface, Tail): - x_tail = -csys * surface.shape_vec[0] + position - y_tail = surface.shape_vec[1] - - ax.plot( - x_tail, - y_tail, - color=vis_args["tail"], - linewidth=vis_args["line_width"], - ) - ax.plot( - x_tail, - -y_tail, - color=vis_args["tail"], - linewidth=vis_args["line_width"], - ) - # close above and below the tail - ax.plot( - [x_tail[-1], x_tail[-1]], - [y_tail[-1], -y_tail[-1]], - color=vis_args["tail"], - linewidth=vis_args["line_width"], - ) - ax.plot( - [x_tail[0], x_tail[0]], - [y_tail[0], -y_tail[0]], - color=vis_args["tail"], - linewidth=vis_args["line_width"], - ) + self._draw_tail(ax, surface, position, drawn_surfaces, vis_args) + elif isinstance(surface, Fins): + self._draw_fins(ax, surface, position, drawn_surfaces, vis_args) + return drawn_surfaces + + def _draw_nose_cone(self, ax, surface, position, drawn_surfaces, vis_args): + """Draws the nosecone and saves the position of the points of interest + for the tubes.""" + x_nosecone = -self.rocket._csys * surface.shape_vec[0] + position + y_nosecone = surface.shape_vec[1] + ax.plot( + x_nosecone, + y_nosecone, + color=vis_args["nose"], + linewidth=vis_args["line_width"], + ) + ax.plot( + x_nosecone, + -y_nosecone, + color=vis_args["nose"], + linewidth=vis_args["line_width"], + ) + # close the nosecone + ax.plot( + [x_nosecone[-1], x_nosecone[-1]], + [y_nosecone[-1], -y_nosecone[-1]], + color=vis_args["nose"], + linewidth=vis_args["line_width"], + ) + # Add the nosecone to the list of drawn surfaces + drawn_surfaces.append( + (surface, x_nosecone[-1], surface.rocket_radius, x_nosecone[-1]) + ) - # Add the tail to the list of drawn surfaces - drawn_surfaces.append( - (surface, position, surface.bottom_radius, x_tail[-1]) - ) + def _draw_tail(self, ax, surface, position, drawn_surfaces, vis_args): + """Draws the tail and saves the position of the points of interest + for the tubes.""" + x_tail = -self.rocket._csys * surface.shape_vec[0] + position + y_tail = surface.shape_vec[1] + ax.plot( + x_tail, y_tail, color=vis_args["tail"], linewidth=vis_args["line_width"] + ) + ax.plot( + x_tail, -y_tail, color=vis_args["tail"], linewidth=vis_args["line_width"] + ) + # close above and below the tail + ax.plot( + [x_tail[-1], x_tail[-1]], + [y_tail[-1], -y_tail[-1]], + color=vis_args["tail"], + linewidth=vis_args["line_width"], + ) + ax.plot( + [x_tail[0], x_tail[0]], + [y_tail[0], -y_tail[0]], + color=vis_args["tail"], + linewidth=vis_args["line_width"], + ) + # 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): + """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)] + + for angle in rotation_angles: + # Create a rotation matrix for the current angle around the x-axis + rotation_matrix = np.array([[1, 0], [0, np.cos(angle)]]) + + # Apply the rotation to the original fin points + rotated_points_2d = np.dot(rotation_matrix, np.vstack((x_fin, y_fin))) + + # 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, + color=vis_args["fins"], + linewidth=vis_args["line_width"], + ) - # Draw fins - elif isinstance(surface, Fins): - num_fins = surface.n - x_fin = -csys * surface.shape_vec[0] + position - y_fin = surface.shape_vec[1] + surface.rocket_radius - - # Calculate the rotation angles for the other two fins (symmetrically) - rotation_angles = [2 * np.pi * i / num_fins for i in range(num_fins)] - - # Apply rotation transformations to get points for the other fins in 2D space - for angle in rotation_angles: - # Create a rotation matrix for the current angle around the x-axis - rotation_matrix = np.array([[1, 0], [0, np.cos(angle)]]) - - # Apply the rotation to the original fin points - rotated_points_2d = np.dot( - rotation_matrix, np.vstack((x_fin, y_fin)) - ) - - # 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 - ) - - # Plot the fins - ax.plot( - x_rotated, - y_rotated, - color=vis_args["fins"], - linewidth=vis_args["line_width"], - ) - - # Add the fin to the list of drawn surfaces - drawn_surfaces.append( - (surface, position, surface.rocket_radius, x_rotated[-1]) - ) + drawn_surfaces.append((surface, position, surface.rocket_radius, x_rotated[-1])) - # Draw tubes + def _draw_tubes(self, ax, drawn_surfaces, vis_args): + """Draws the tubes between the aerodynamic surfaces.""" for i, d_surface in enumerate(drawn_surfaces): # Draw the tubes, from the end of the first surface to the beginning # of the next surface, with the radius of the rocket at that point @@ -368,18 +376,41 @@ def draw(self, vis_args=None): color=vis_args["body"], linewidth=vis_args["line_width"], ) + return radius, last_x - # Draw motor + def _draw_motor(self, last_radius, last_x, ax, vis_args): + """Draws the motor from motor patches""" total_csys = self.rocket._csys * self.rocket.motor._csys nozzle_position = ( self.rocket.motor_position + self.rocket.motor.nozzle_position * total_csys ) - # List of motor patches + # Get motor patches translated to the correct position + motor_patches = self._generate_motor_patches(total_csys, ax, vis_args) + + # Draw patches + if not isinstance(self.rocket.motor, EmptyMotor): + # Add nozzle last so it is in front of the other patches + nozzle = self.rocket.motor.plots._generate_nozzle( + translate=(nozzle_position, 0), csys=self.rocket._csys + ) + motor_patches += [nozzle] + + outline = self.rocket.motor.plots._generate_motor_region( + list_of_patches=motor_patches + ) + # add outline first so it is behind the other patches + ax.add_patch(outline) + for patch in motor_patches: + ax.add_patch(patch) + + self._draw_nozzle_tube(last_radius, last_x, nozzle_position, ax, vis_args) + + def _generate_motor_patches(self, total_csys, ax, vis_args): + """Generates motor patches for drawing""" motor_patches = [] - # Get motor patches translated to the correct position - if isinstance(self.rocket.motor, (SolidMotor)): + if isinstance(self.rocket.motor, SolidMotor): grains_cm_position = ( self.rocket.motor_position + self.rocket.motor.grains_center_of_mass_position * total_csys @@ -452,27 +483,17 @@ def draw(self, vis_args=None): ) motor_patches += [tank] - # add nozzle last so it is in front of the other patches - if not isinstance(self.rocket.motor, EmptyMotor): - nozzle = self.rocket.motor.plots._generate_nozzle( - translate=(nozzle_position, 0), csys=self.rocket._csys - ) - motor_patches += [nozzle] - outline = self.rocket.motor.plots._generate_motor_region( - list_of_patches=motor_patches - ) - # add outline first so it is behind the other patches - ax.add_patch(outline) - for patch in motor_patches: - ax.add_patch(patch) + return motor_patches + def _draw_nozzle_tube(self, last_radius, last_x, nozzle_position, ax, vis_args): + """Draws the tube from the last surface to the nozzle position.""" # Check if nozzle is beyond the last surface, if so draw a tube # to it, with the radius of the last surface if self.rocket._csys == 1: if nozzle_position < last_x: x_tube = [last_x, nozzle_position] - y_tube = [radius, radius] - y_tube_negated = [-radius, -radius] + y_tube = [last_radius, last_radius] + y_tube_negated = [-last_radius, -last_radius] ax.plot( x_tube, @@ -489,8 +510,8 @@ def draw(self, vis_args=None): else: # if self.rocket._csys == -1: if nozzle_position > last_x: x_tube = [last_x, nozzle_position] - y_tube = [radius, radius] - y_tube_negated = [-radius, -radius] + y_tube = [last_radius, last_radius] + y_tube_negated = [-last_radius, -last_radius] ax.plot( x_tube, @@ -505,11 +526,12 @@ def draw(self, vis_args=None): linewidth=vis_args["line_width"], ) - # Draw rail buttons + def _draw_rail_buttons(self, ax, vis_args): + """Draws the rail buttons of the rocket.""" try: buttons, pos = self.rocket.rail_buttons[0] lower = pos - upper = pos + buttons.buttons_distance * csys + upper = pos + buttons.buttons_distance * self.rocket._csys ax.scatter( lower, -self.rocket.radius, marker="s", color=vis_args["buttons"], s=15 ) @@ -519,6 +541,8 @@ def draw(self, vis_args=None): except IndexError: pass + def _draw_center_of_mass_and_pressure(self, ax): + """Draws the center of mass and center of pressure of the rocket.""" # Draw center of mass and center of pressure cm = self.rocket.center_of_mass(0) ax.scatter(cm, 0, color="#1565c0", label="Center of Mass", s=10) @@ -528,18 +552,6 @@ def draw(self, vis_args=None): cp, 0, label="Static Center of Pressure", color="red", s=10, zorder=10 ) - # Set plot attributes - plt.title("Rocket Representation") - plt.xlim() - plt.ylim([-self.rocket.radius * 4, self.rocket.radius * 6]) - plt.xlabel("Position (m)") - plt.ylabel("Radius (m)") - plt.legend(bbox_to_anchor=(1.05, 1), loc="upper left") - plt.tight_layout() - plt.show() - - return None - def all(self): """Prints out all graphs available about the Rocket. It simply calls all the other plotter methods in this class. diff --git a/rocketpy/prints/flight_prints.py b/rocketpy/prints/flight_prints.py index 5b1f9b914..06db442f8 100644 --- a/rocketpy/prints/flight_prints.py +++ b/rocketpy/prints/flight_prints.py @@ -1,3 +1,17 @@ +"""rocketpy/prints/flight_prints.py + +This module contains the _FlightPrints class, which is responsible for printing +flight information in a user-friendly manner. + +Notes +----- +- This module does not have any external dependencies (avoid importing libraries). +- We assume that all flight information is valid, no validation checks is run. +- Avoid calculating values here, only print the values from the Flight class. +- The _FlightPrints is a private class, it is subjected to change without notice. +""" + + class _FlightPrints: """Class that holds prints methods for Flight class. @@ -5,7 +19,6 @@ class _FlightPrints: ---------- _FlightPrints.flight : Flight Flight object that will be used for the prints. - """ def __init__( @@ -24,235 +37,191 @@ def __init__( None """ self.flight = flight - return None def initial_conditions(self): - """Prints all initial conditions data available about the Flight. + """Prints initial conditions data available about the flight, including + position, velocity, attitude, euler angles, angular velocity, and + stability margin. Returns ------- None """ - print("\nInitial Conditions\n") - # Post-process results - if self.flight.post_processed is False: - self.flight.post_process() + t_init = self.flight.time[0] + + print(f"Initial time: {t_init:.3f} s") print( - "Position - x: {:.2f} m | y: {:.2f} m | z: {:.2f} m".format( - self.flight.x(0), self.flight.y(0), self.flight.z(0) - ) + f"Position - x: {self.flight.x(t_init):.2f} m | " + f"y: {self.flight.y(t_init):.2f} m | " + f"z: {self.flight.z(t_init):.2f} m" ) print( - "Velocity - Vx: {:.2f} m/s | Vy: {:.2f} m/s | Vz: {:.2f} m/s".format( - self.flight.vx(0), self.flight.vy(0), self.flight.vz(0) - ) + f"Velocity - Vx: {self.flight.vx(t_init):.2f} m/s | " + f"Vy: {self.flight.vy(t_init):.2f} m/s | " + f"Vz: {self.flight.vz(t_init):.2f} m/s" ) print( - "Attitude - e0: {:.3f} | e1: {:.3f} | e2: {:.3f} | e3: {:.3f}".format( - self.flight.e0(0), - self.flight.e1(0), - self.flight.e2(0), - self.flight.e3(0), - ) + f"Attitude (quaternions) - e0: {self.flight.e0(t_init):.3f} | " + f"e1: {self.flight.e1(t_init):.3f} | " + f"e2: {self.flight.e2(t_init):.3f} | " + f"e3: {self.flight.e3(t_init):.3f}" ) print( - "Euler Angles - Spin φ : {:.2f}° | Nutation θ: {:.2f}° | Precession ψ: {:.2f}°".format( - self.flight.phi(0), self.flight.theta(0), self.flight.psi(0) - ) + f"Euler Angles - Spin φ : {self.flight.phi(t_init):.2f}° | " + f"Nutation θ: {self.flight.theta(t_init):.2f}° | " + f"Precession ψ: {self.flight.psi(t_init):.2f}°" ) print( - "Angular Velocity - ω1: {:.2f} rad/s | ω2: {:.2f} rad/s| ω3: {:.2f} rad/s".format( - self.flight.w1(0), self.flight.w2(0), self.flight.w3(0) - ) + f"Angular Velocity - ω1: {self.flight.w1(t_init):.2f} rad/s | " + f"ω2: {self.flight.w2(t_init):.2f} rad/s | " + f"ω3: {self.flight.w3(t_init):.2f} rad/s" ) - - return None + print(f"Initial Stability Margin: {self.flight.initial_stability_margin:.3f} c") def numerical_integration_settings(self): - """Prints out the Numerical Integration settings available about the - flight. + """Prints out the numerical integration settings available about the + flight, this includes the maximum allowed flight time, maximum allowed + time step, and other settings. Returns ------- None """ - print("\nNumerical Integration Settings\n") - print("Maximum Allowed Flight Time: {:f} s".format(self.flight.max_time)) - print("Maximum Allowed Time Step: {:f} s".format(self.flight.max_time_step)) - print("Minimum Allowed Time Step: {:e} s".format(self.flight.min_time_step)) - print("Relative Error Tolerance: ", self.flight.rtol) - print("Absolute Error Tolerance: ", self.flight.atol) - print("Allow Event Overshoot: ", self.flight.time_overshoot) - print("Terminate Simulation on Apogee: ", self.flight.terminate_on_apogee) - print("Number of Time Steps Used: ", len(self.flight.time_steps)) + print(f"Maximum Allowed Flight Time: {self.flight.max_time:.2f} s") + print(f"Maximum Allowed Time Step: {self.flight.max_time_step:.2f} s") + print(f"Minimum Allowed Time Step: {self.flight.min_time_step:.2e} s") + print(f"Relative Error Tolerance: {self.flight.rtol}") + print(f"Absolute Error Tolerance: {self.flight.atol}") + print(f"Allow Event Overshoot: {self.flight.time_overshoot}") + print(f"Terminate Simulation on Apogee: {self.flight.terminate_on_apogee}") + print(f"Number of Time Steps Used: {len(self.flight.time_steps)}") print( - "Number of Derivative Functions Evaluation: ", - sum(self.flight.function_evaluations_per_time_step), + "Number of Derivative Functions Evaluation: " + f"{sum(self.flight.function_evaluations_per_time_step)}" ) + avg_func_evals_per_step = sum( + self.flight.function_evaluations_per_time_step + ) / len(self.flight.time_steps) print( - "Average Function Evaluations per Time Step: {:3f}".format( - sum(self.flight.function_evaluations_per_time_step) - / len(self.flight.time_steps) - ) + f"Average Function Evaluations per Time Step: {avg_func_evals_per_step:.3f}" ) - return None - def surface_wind_conditions(self): - """Prints out the Surface Wind Conditions available about the flight. + """Prints out the Surface Wind Conditions for the flight. Returns ------- None """ - if self.flight.post_processed is False: - self.flight.post_process() print("\nSurface Wind Conditions\n") - print( - "Frontal Surface Wind Speed: {:.2f} m/s".format( - self.flight.frontal_surface_wind - ) - ) - print( - "Lateral Surface Wind Speed: {:.2f} m/s".format( - self.flight.lateral_surface_wind - ) - ) - - return None + print(f"Frontal Surface Wind Speed: {self.flight.frontal_surface_wind:.2f} m/s") + print(f"Lateral Surface Wind Speed: {self.flight.lateral_surface_wind:.2f} m/s") def launch_rail_conditions(self): - """Prints out the Launch Rail Conditions available about the flight. + """Prints out the Launch Rail Conditions available about the flight, + including the length, inclination, and heading of the launch rail. Returns ------- None """ - print("\nLaunch Rail\n") - print("Launch Rail Length:", self.flight.rail_length, " m") - print("Launch Rail Inclination: {:.2f}°".format(self.flight.inclination)) - print("Launch Rail Heading: {:.2f}°".format(self.flight.heading)) - return None + print(f"Launch Rail Length: {self.flight.rail_length} m") + print(f"Launch Rail Inclination: {self.flight.inclination:.2f}°") + print(f"Launch Rail Heading: {self.flight.heading:.2f}°") def out_of_rail_conditions(self): - """Prints out the Out of Rail Conditions available about the flight. + """Prints out the Out of Rail Conditions available about the flight, + including the time, velocity, stability margin, angle of attack, thrust + to weight ratio, and Reynolds number. Returns ------- None """ - if self.flight.post_processed is False: - self.flight.post_process() print("\nRail Departure State\n") - print("Rail Departure Time: {:.3f} s".format(self.flight.out_of_rail_time)) - print( - "Rail Departure Velocity: {:.3f} m/s".format( - self.flight.out_of_rail_velocity - ) - ) + print(f"Rail Departure Time: {self.flight.out_of_rail_time:.3f} s") + print(f"Rail Departure Velocity: {self.flight.out_of_rail_velocity:.3f} m/s") print( - "Rail Departure Stability Margin: {:.3f} c".format( - self.flight.stability_margin(self.flight.out_of_rail_time) - ) + "Rail Departure Stability Margin: " + f"{self.flight.out_of_rail_stability_margin:.3f} c" ) print( - "Rail Departure Angle of Attack: {:.3f}°".format( - self.flight.angle_of_attack(self.flight.out_of_rail_time) - ) + "Rail Departure Angle of Attack: " + f"{self.flight.angle_of_attack(self.flight.out_of_rail_time):.3f}°" ) print( - "Rail Departure Thrust-Weight Ratio: {:.3f}".format( - self.flight.rocket.thrust_to_weight(self.flight.out_of_rail_time) - ) + "Rail Departure Thrust-Weight Ratio: " + f"{self.flight.rocket.thrust_to_weight(self.flight.out_of_rail_time):.3f}" ) print( - "Rail Departure Reynolds Number: {:.3e}".format( - self.flight.reynolds_number(self.flight.out_of_rail_time) - ) + "Rail Departure Reynolds Number: " + f"{self.flight.reynolds_number(self.flight.out_of_rail_time):.3e}" ) - return None - def burn_out_conditions(self): - """Prints out the Burn Out Conditions available about the flight. + """Prints out the Burn Out Conditions available about the flight, + including the burn out time, altitude, speed, freestream speed, Mach + number, and kinetic energy. Returns ------- None """ print("\nBurn out State\n") - print("Burn out time: {:.3f} s".format(self.flight.rocket.motor.burn_out_time)) + print(f"Burn out time: {self.flight.rocket.motor.burn_out_time:.3f} s") print( - "Altitude at burn out: {:.3f} m (AGL)".format( - self.flight.z(self.flight.rocket.motor.burn_out_time) - - self.flight.env.elevation - ) + "Altitude at burn out: " + f"{self.flight.z(self.flight.rocket.motor.burn_out_time):.3f} m (ASL) | " + f"{self.flight.altitude(self.flight.rocket.motor.burn_out_time):.3f} " + "m (AGL)" ) print( - "Rocket velocity at burn out: {:.3f} m/s".format( - self.flight.speed(self.flight.rocket.motor.burn_out_time) - ) + "Rocket speed at burn out: " + f"{self.flight.speed(self.flight.rocket.motor.burn_out_time):.3f} m/s" ) - print( - "Freestream velocity at burn out: {:.3f} m/s".format( - ( - self.flight.stream_velocity_x( - self.flight.rocket.motor.burn_out_time - ) - ** 2 - + self.flight.stream_velocity_y( - self.flight.rocket.motor.burn_out_time - ) - ** 2 - + self.flight.stream_velocity_z( - self.flight.rocket.motor.burn_out_time - ) - ** 2 - ) - ** 0.5 - ) + + stream_velocity = self.flight.free_stream_speed( + self.flight.rocket.motor.burn_out_time ) + print(f"Freestream velocity at burn out: {stream_velocity:.3f} m/s") + print( - "Mach Number at burn out: {:.3f}".format( - self.flight.mach_number(self.flight.rocket.motor.burn_out_time) - ) + "Mach Number at burn out: " + f"{self.flight.mach_number(self.flight.rocket.motor.burn_out_time):.3f}" ) print( - "Kinetic energy at burn out: {:.3e} J".format( - self.flight.kinetic_energy(self.flight.rocket.motor.burn_out_time) - ) + "Kinetic energy at burn out: " + f"{self.flight.kinetic_energy(self.flight.rocket.motor.burn_out_time):.3e} " + "J" ) - return None - def apogee_conditions(self): - """Prints out the Apogee Conditions available about the flight. + """Prints out the Apogee Conditions available about the flight, + including the apogee time, altitude, freestream speed, latitude, and + longitude. Returns ------- None """ - if self.flight.post_processed is False: - self.flight.post_process() print("\nApogee State\n") + print(f"Apogee Time: {self.flight.apogee_time:.3f} s") print( - "Apogee Altitude: {:.3f} m (ASL) | {:.3f} m (AGL)".format( - self.flight.apogee, self.flight.apogee - self.flight.env.elevation - ) + f"Apogee Altitude: {self.flight.apogee:.3f} m (ASL) | " + f"{self.flight.altitude(self.flight.apogee_time):.3f} m (AGL)" ) - print("Apogee Time: {:.3f} s".format(self.flight.apogee_time)) + print(f"Apogee Freestream Speed: {self.flight.apogee_freestream_speed:.3f} m/s") + print(f"Apogee X position: {self.flight.x(self.flight.apogee_time):.3f} m") + print(f"Apogee Y position: {self.flight.y(self.flight.apogee_time):.3f} m") + print(f"Apogee latitude: {self.flight.latitude(self.flight.apogee_time):.7f}°") print( - "Apogee Freestream Speed: {:.3f} m/s".format( - self.flight.apogee_freestream_speed - ) + f"Apogee longitude: {self.flight.longitude(self.flight.apogee_time):.7f}°" ) - return None - def events_registered(self): """Prints out the Events Registered available about the flight. @@ -260,8 +229,6 @@ def events_registered(self): ------- None """ - if self.flight.post_processed is False: - self.flight.post_process() print("\nParachute Events\n") if len(self.flight.parachute_events) == 0: print("No Parachute Events Were Triggered.") @@ -269,24 +236,17 @@ def events_registered(self): trigger_time = event[0] parachute = event[1] open_time = trigger_time + parachute.lag - velocity = self.flight.free_stream_speed(open_time) + speed = self.flight.free_stream_speed(open_time) altitude = self.flight.z(open_time) name = parachute.name.title() - print(name + " Ejection Triggered at: {:.3f} s".format(trigger_time)) - print(name + " Parachute Inflated at: {:.3f} s".format(open_time)) - print( - name - + " Parachute Inflated with Freestream Speed of: {:.3f} m/s".format( - velocity - ) - ) + print(f"Parachute: {name}") + print(f"\tEjection time: {trigger_time:.3f} s") + print(f"\tInflation time: {open_time:.3f} s") + print(f"\tFreestream speed at inflation: {speed:.3f} m/s") print( - name - + " Parachute Inflated at Height of: {:.3f} m (AGL)".format( - altitude - self.flight.env.elevation - ) + f"\tAltitude at inflation: {altitude:.3f} m (ASL) | " + f"{self.flight.altitude(trigger_time):.3f} m (AGL)" ) - return None def impact_conditions(self): """Prints out the Impact Conditions available about the flight. @@ -295,27 +255,36 @@ def impact_conditions(self): ------- None """ - if self.flight.post_processed is False: - self.flight.post_process() if len(self.flight.impact_state) != 0: print("\nImpact Conditions\n") - print("X Impact: {:.3f} m".format(self.flight.x_impact)) - print("Y Impact: {:.3f} m".format(self.flight.y_impact)) - print("Latitude: {:.7f}°".format(self.flight.latitude(self.flight.t_final))) + print(f"Time of impact: {self.flight.t_final:.3f} s") + print(f"X impact: {self.flight.x_impact:.3f} m") + print(f"Y impact: {self.flight.y_impact:.3f} m") print( - "Longitude: {:.7f}°".format(self.flight.longitude(self.flight.t_final)) + f"Altitude impact: {self.flight.z(self.flight.t_final):.3f} m (ASL) | " + f"{self.flight.altitude(self.flight.t_final):.3f} m (AGL) " + ) + print(f"Latitude: {self.flight.latitude(self.flight.t_final):.7f}°") + print(f"Longitude: {self.flight.longitude(self.flight.t_final):.7f}°") + print(f"Vertical velocity at impact: {self.flight.impact_velocity:.3f} m/s") + num_parachute_events = sum( + 1 + for event in self.flight.parachute_events + if event[0] < self.flight.t_final + ) + print( + f"Number of parachutes triggered until impact: {num_parachute_events}" ) - print("Time of Impact: {:.3f} s".format(self.flight.t_final)) - print("Velocity at Impact: {:.3f} m/s".format(self.flight.impact_velocity)) elif self.flight.terminate_on_apogee is False: print("End of Simulation") - t_final = self.flight.solution[-1][0] - print("Time: {:.3f} s".format(t_final)) - print("Altitude: {:.3f} m".format(self.flight.solution[-1][3])) - print("Latitude: {:.3f}°".format(self.flight.latitude(t_final))) - print("Longitude: {:.3f}°".format(self.flight.longitude(t_final))) - - return None + t_final = self.flight.time[-1] + print(f"Time: {t_final:.3f} s") + print( + f"Altitude: {self.flight.z(t_final)} m (ASL) | " + f"{self.flight.altitude(t_final):.3f} m (AGL)" + ) + print(f"Latitude: {self.flight.latitude(t_final):.7f}°") + print(f"Longitude: {self.flight.longitude(t_final):.7f}°") def maximum_values(self): """Prints out the Maximum Values available about the flight. @@ -326,53 +295,44 @@ def maximum_values(self): """ print("\nMaximum Values\n") print( - "Maximum Speed: {:.3f} m/s at {:.2f} s".format( - self.flight.max_speed, self.flight.max_speed_time - ) + f"Maximum Speed: {self.flight.max_speed:.3f} m/s " + f"at {self.flight.max_speed_time:.2f} s" ) print( - "Maximum Mach Number: {:.3f} Mach at {:.2f} s".format( - self.flight.max_mach_number, self.flight.max_mach_number_time - ) + f"Maximum Mach Number: {self.flight.max_mach_number:.3f} Mach " + f"at {self.flight.max_mach_number_time:.2f} s" ) print( - "Maximum Reynolds Number: {:.3e} at {:.2f} s".format( - self.flight.max_reynolds_number, self.flight.max_reynolds_number_time - ) + f"Maximum Reynolds Number: {self.flight.max_reynolds_number:.3e} " + f"at {self.flight.max_reynolds_number_time:.2f} s" ) print( - "Maximum Dynamic Pressure: {:.3e} Pa at {:.2f} s".format( - self.flight.max_dynamic_pressure, self.flight.max_dynamic_pressure_time - ) + f"Maximum Dynamic Pressure: {self.flight.max_dynamic_pressure:.3e} Pa " + f"at {self.flight.max_dynamic_pressure_time:.2f} s" ) print( - "Maximum Acceleration During Motor Burn: {:.3f} m/s² at {:.2f} s".format( - self.flight.max_acceleration_power_on, - self.flight.max_acceleration_power_on_time, - ) + "Maximum Acceleration During Motor Burn: " + f"{self.flight.max_acceleration_power_on:.3f} m/s² " + f"at {self.flight.max_acceleration_power_on_time:.2f} s" ) print( - "Maximum Gs During Motor Burn: {:.3f} g at {:.2f} s".format( - self.flight.max_acceleration_power_on / self.flight.env.standard_g, - self.flight.max_acceleration_power_on_time, - ) + "Maximum Gs During Motor Burn: " + f"{self.flight.max_acceleration_power_on / self.flight.env.standard_g:.3f} " + f"g at {self.flight.max_acceleration_power_on_time:.2f} s" ) print( - "Maximum Acceleration After Motor Burn: {:.3f} m/s² at {:.2f} s".format( - self.flight.max_acceleration_power_off, - self.flight.max_acceleration_power_off_time, - ) + "Maximum Acceleration After Motor Burn: " + f"{self.flight.max_acceleration_power_off:.3f} m/s² " + f"at {self.flight.max_acceleration_power_off_time:.2f} s" ) print( - "Maximum Gs After Motor Burn: {:.3f} g at {:.2f} s".format( - self.flight.max_acceleration_power_off / self.flight.env.standard_g, - self.flight.max_acceleration_power_off_time, - ) + "Maximum Gs After Motor Burn: " + f"{self.flight.max_acceleration_power_off / self.flight.env.standard_g:.3f}" + f" Gs at {self.flight.max_acceleration_power_off_time:.2f} s" ) print( - "Maximum Stability Margin: {:.3f} c at {:.2f} s".format( - self.flight.max_stability_margin, self.flight.max_stability_margin_time - ) + f"Maximum Stability Margin: {self.flight.max_stability_margin:.3f} c " + f"at {self.flight.max_stability_margin_time:.2f} s" ) if ( @@ -382,93 +342,92 @@ def maximum_values(self): pass else: print( - "Maximum Upper Rail Button Normal Force: {:.3f} N".format( - self.flight.max_rail_button1_normal_force - ) + "Maximum Upper Rail Button Normal Force: " + f"{self.flight.max_rail_button1_normal_force:.3f} N" ) print( - "Maximum Upper Rail Button Shear Force: {:.3f} N".format( - self.flight.max_rail_button1_shear_force - ) + "Maximum Upper Rail Button Shear Force: " + f"{self.flight.max_rail_button1_shear_force:.3f} N" ) print( - "Maximum Lower Rail Button Normal Force: {:.3f} N".format( - self.flight.max_rail_button2_normal_force - ) + "Maximum Lower Rail Button Normal Force: " + f"{self.flight.max_rail_button2_normal_force:.3f} N" ) print( - "Maximum Lower Rail Button Shear Force: {:.3f} N".format( - self.flight.max_rail_button2_shear_force - ) + "Maximum Lower Rail Button Shear Force: " + f"{self.flight.max_rail_button2_shear_force:.3f} N" ) - return None def stability_margin(self): - """Prints out the maximum and minimum stability margin available - about the flight.""" + """Prints out the stability margins of the flight at different times. + + This method prints the following: Initial Stability Margin, Out of Rail + Stability Margin, Maximum Stability Margin, and Minimum Stability Margin + + Each stability margin is printed along with the time it occurred. + + Notes + ----- + The stability margin is typically measured in calibers (c), where 1 + caliber is the diameter of the rocket. + """ print("\nStability Margin\n") print( - "Maximum Stability Margin: {:.3f} c at {:.2f} s".format( - self.flight.max_stability_margin, self.flight.max_stability_margin_time - ) + f"Initial Stability Margin: {self.flight.initial_stability_margin:.3f} c " + f"at {self.flight.time[0]:.2f} s" ) print( - "Minimum Stability Margin: {:.3f} c at {:.2f} s".format( - self.flight.min_stability_margin, self.flight.min_stability_margin_time - ) + "Out of Rail Stability Margin: " + f"{self.flight.out_of_rail_stability_margin:.3f} c " + f"at {self.flight.out_of_rail_time:.2f} s" + ) + print( + f"Maximum Stability Margin: {self.flight.max_stability_margin:.3f} c " + f"at {self.flight.max_stability_margin_time:.2f} s" + ) + print( + f"Minimum Stability Margin: {self.flight.min_stability_margin:.3f} c " + f"at {self.flight.min_stability_margin_time:.2f} s" ) - return None def all(self): - """Prints out all data available about the Flight. + """Prints out all data available about the Flight. This method invokes + all other print methods in the class. Returns ------- None """ - # Print initial conditions self.initial_conditions() print() - # Print surface wind conditions self.surface_wind_conditions() print() - # Print launch rail orientation self.launch_rail_conditions() print() - # Print out of rail conditions self.out_of_rail_conditions() print() - # Print burn out conditions self.burn_out_conditions() print() - # Print apogee conditions self.apogee_conditions() print() - # Print events registered self.events_registered() print() - # Print impact conditions self.impact_conditions() print() - # Print stability margin self.stability_margin() print() - # Print maximum values self.maximum_values() print() - # Print Numerical Integration Information self.numerical_integration_settings() print() - - return None diff --git a/rocketpy/rocket/aero_surface.py b/rocketpy/rocket/aero_surface.py index c5d154f3e..d41240ac9 100644 --- a/rocketpy/rocket/aero_surface.py +++ b/rocketpy/rocket/aero_surface.py @@ -1,6 +1,5 @@ import warnings from abc import ABC, abstractmethod -from functools import cached_property import numpy as np from scipy.optimize import fsolve diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index ffcf3b1c8..c65f1eca1 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -18,6 +18,7 @@ ) from rocketpy.rocket.components import Components from rocketpy.rocket.parachute import Parachute +from rocketpy.tools import parallel_axis_theorem_from_com class Rocket: @@ -61,6 +62,9 @@ class Rocket: pointing from the rocket's nose cone to the rocket's tail. Rocket.mass : float Rocket's mass without motor and propellant, measured in kg. + Rocket.dry_mass : float + Rocket's mass without propellant, measured in kg. It does include the + motor mass. Rocket.center_of_mass : Function Position of the rocket's center of mass, including propellant, relative to the user defined rocket reference system. @@ -77,6 +81,9 @@ class Rocket: Function of time expressing the total mass of the rocket, defined as the sum of the propellant mass and the rocket mass without propellant. + Rocket.total_mass_flow_rate : Function + Time derivative of rocket's total mass in kg/s as a function + of time as obtained by the thrust source of the added motor. Rocket.thrust_to_weight : Function Function of time expressing the motor thrust force divided by rocket weight. The gravitational acceleration is assumed as 9.80665 m/s^2. @@ -394,22 +401,18 @@ def evaluate_total_mass(self): def evaluate_dry_mass(self): """Calculates and returns the rocket's dry mass. The dry mass is defined as the sum of the motor's dry mass and the - rocket mass without motor. The function returns an object - of the Function class and is defined as a function of time. + rocket mass without motor. Returns ------- - self.total_mass : Function - Function of time expressing the total mass of the rocket, - defined as the sum of the propellant mass and the rocket - mass without propellant. + self.dry_mass : float + Rocket's dry mass (Rocket + Motor) (kg) """ # Make sure there is a motor associated with the rocket if self.motor is None: print("Please associate this rocket with a motor!") return False - # Calculate total dry mass: motor (without propellant) + rocket mass self.dry_mass = self.mass + self.motor.dry_mass return self.dry_mass @@ -543,7 +546,11 @@ def evaluate_stability_margin(self): """ self.stability_margin.set_source( lambda mach, time: ( - (self.center_of_mass(time) - self.cp_position(mach)) / (2 * self.radius) + ( + self.center_of_mass.get_value_opt(time) + - self.cp_position.get_value_opt(mach) + ) + / (2 * self.radius) ) * self._csys ) @@ -561,7 +568,10 @@ def evaluate_static_margin(self): """ # Calculate static margin self.static_margin.set_source( - lambda time: (self.center_of_mass(time) - self.cp_position(0)) + lambda time: ( + self.center_of_mass.get_value_opt(time) + - self.cp_position.get_value_opt(0) + ) / (2 * self.radius) ) # Change sign if coordinate system is upside down @@ -620,6 +630,10 @@ def evaluate_dry_inertias(self): ---------- .. [1] https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ + # Get masses + motor_dry_mass = self.motor.dry_mass + mass = self.mass + # Compute axes distances noMCM_to_CDM = ( self.center_of_mass_without_motor - self.center_of_dry_mass_position @@ -629,18 +643,18 @@ def evaluate_dry_inertias(self): ) # Compute dry inertias - self.dry_I_11 = ( - self.I_11_without_motor - + self.mass * noMCM_to_CDM**2 - + self.motor.dry_I_11 - + self.motor.dry_mass * motorCDM_to_CDM**2 + self.dry_I_11 = parallel_axis_theorem_from_com( + self.I_11_without_motor, mass, noMCM_to_CDM + ) + parallel_axis_theorem_from_com( + self.motor.dry_I_11, motor_dry_mass, motorCDM_to_CDM ) - self.dry_I_22 = ( - self.I_22_without_motor - + self.mass * noMCM_to_CDM**2 - + self.motor.dry_I_22 - + self.motor.dry_mass * motorCDM_to_CDM**2 + + self.dry_I_22 = parallel_axis_theorem_from_com( + self.I_22_without_motor, mass, noMCM_to_CDM + ) + parallel_axis_theorem_from_com( + self.motor.dry_I_22, motor_dry_mass, motorCDM_to_CDM ) + self.dry_I_33 = self.I_33_without_motor + self.motor.dry_I_33 self.dry_I_12 = self.I_12_without_motor + self.motor.dry_I_12 self.dry_I_13 = self.I_13_without_motor + self.motor.dry_I_13 @@ -697,18 +711,14 @@ def evaluate_inertias(self): CM_to_CPM = self.center_of_mass - self.center_of_propellant_position # Compute inertias - self.I_11 = ( - self.dry_I_11 - + self.motor.I_11 - + dry_mass * CM_to_CDM**2 - + prop_mass * CM_to_CPM**2 - ) - self.I_22 = ( - self.dry_I_22 - + self.motor.I_22 - + dry_mass * CM_to_CDM**2 - + prop_mass * CM_to_CPM**2 - ) + self.I_11 = parallel_axis_theorem_from_com( + self.dry_I_11, dry_mass, CM_to_CDM + ) + parallel_axis_theorem_from_com(self.motor.I_11, prop_mass, CM_to_CPM) + + self.I_22 = parallel_axis_theorem_from_com( + self.dry_I_22, dry_mass, CM_to_CDM + ) + parallel_axis_theorem_from_com(self.motor.I_22, prop_mass, CM_to_CPM) + self.I_33 = self.dry_I_33 + self.motor.I_33 self.I_12 = self.dry_I_12 + self.motor.I_12 self.I_13 = self.dry_I_13 + self.motor.I_13 @@ -764,6 +774,7 @@ def add_motor(self, motor, position): self.motor.center_of_dry_mass_position * _ + self.motor_position ) self.nozzle_position = self.motor.nozzle_position * _ + self.motor_position + self.total_mass_flow_rate = self.motor.total_mass_flow_rate self.evaluate_dry_mass() self.evaluate_total_mass() self.evaluate_center_of_dry_mass() diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 23c16fcf1..c03e1e41a 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -12,7 +12,10 @@ from ..plots.flight_plots import _FlightPlots from ..prints.flight_prints import _FlightPrints from ..tools import ( + calculate_cubic_hermite_coefficients, find_closest, + find_root_linear_interpolation, + find_roots_cubic_function, quaternions_to_nutation, quaternions_to_precession, quaternions_to_spin, @@ -182,13 +185,13 @@ class Flight: Flight.function_evaluations : array List that stores number of derivative function evaluations during numerical integration in cumulative manner. - Flight.function_evaluations_per_time_step : array + Flight.function_evaluations_per_time_step : list List that stores number of derivative function evaluations per time step during numerical integration. Flight.time_steps : array List of time steps taking during numerical integration in seconds. - Flight.FlightPhases : Flight.FlightPhases + Flight.flight_phases : Flight.FlightPhases Stores and manages flight phases. Flight.wind_velocity_x : Function Wind velocity X (East) experienced by the rocket as a @@ -426,6 +429,10 @@ class Flight: Rocket's static margin during flight in calibers. Flight.stability_margin : Function Rocket's stability margin during flight, in calibers. + Flight.initial_stability_margin : float + Rocket's initial stability margin in calibers. + Flight.out_of_rail_stability_margin : float + Rocket's stability margin in calibers when it leaves the rail. Flight.stream_velocity_x : Function Freestream velocity x (East) component, in m/s, as a function of time. Can be called or accessed as array. @@ -576,13 +583,7 @@ def __init__( ------- None """ - # Fetch helper classes and functions - FlightPhases = self.FlightPhases - TimeNodes = self.TimeNodes - time_iterator = self.time_iterator - - # Save rocket, parachutes, environment, maximum simulation time - # and termination events + # Save arguments self.env = environment self.rocket = rocket self.rail_length = rail_length @@ -604,38 +605,45 @@ def __init__( self.equations_of_motion = equations_of_motion # Flight initialization - self.__init_post_process_variables() self.__init_solution_monitors() self.__init_equations_of_motion() - - # Initialize prints and plots objects - self.prints = _FlightPrints(self) - self.plots = _FlightPlots(self) - - # Initialize solver monitors self.__init_solver_monitors() # Create known flight phases - self.FlightPhases = FlightPhases() - self.FlightPhases.add_phase( + self.flight_phases = self.FlightPhases() + self.flight_phases.add_phase( self.t_initial, self.initial_derivative, clear=False ) - self.FlightPhases.add_phase(self.max_time) + self.flight_phases.add_phase(self.max_time) # Simulate flight - for phase_index, phase in time_iterator(self.FlightPhases): - # print('\nCurrent Flight Phase List') - # print(self.FlightPhases) - # print('\n\tCurrent Flight Phase') - # print('\tIndex: ', phase_index, ' | Phase: ', phase) + self.__simulate(verbose) + + # Initialize prints and plots objects + self.prints = _FlightPrints(self) + self.plots = _FlightPlots(self) + + def __repr__(self): + return ( + f"" + ) + + def __simulate(self, verbose): + """Simulate the flight trajectory.""" + for phase_index, phase in self.time_iterator(self.flight_phases): # Determine maximum time for this flight phase - phase.time_bound = self.FlightPhases[phase_index + 1].t + phase.time_bound = self.flight_phases[phase_index + 1].t # Evaluate callbacks for callback in phase.callbacks: callback(self) - # Create solver for this flight phase + # Create solver for this flight phase # TODO: allow different integrators self.function_evaluations.append(0) phase.solver = integrate.LSODA( phase.derivative, @@ -647,46 +655,36 @@ def __init__( rtol=self.rtol, atol=self.atol, ) - # print('\n\tSolver Initialization Details') - # print('\t_initial Time: ', phase.t) - # print('\t_initial State: ', self.y_sol) - # print('\tTime Bound: ', phase.time_bound) - # print('\tMin Step: ', self.min_time_step) - # print('\tMax Step: ', self.max_time_step) - # print('\tTolerances: ', self.rtol, self.atol) # Initialize phase time nodes - phase.TimeNodes = TimeNodes() - # Add first time node to permanent list - phase.TimeNodes.add_node(phase.t, [], []) + phase.time_nodes = self.TimeNodes() + # Add first time node to the time_nodes list + phase.time_nodes.add_node(phase.t, [], []) # Add non-overshootable parachute time nodes if self.time_overshoot is False: - phase.TimeNodes.add_parachutes( + # TODO: move parachutes to controllers + phase.time_nodes.add_parachutes( self.parachutes, phase.t, phase.time_bound ) - phase.TimeNodes.add_controllers( + phase.time_nodes.add_controllers( self._controllers, phase.t, phase.time_bound ) - # Add lst time node to permanent list - phase.TimeNodes.add_node(phase.time_bound, [], []) - # Sort time nodes - phase.TimeNodes.sort() - # Merge equal time nodes - phase.TimeNodes.merge() + # Add last time node to the time_nodes list + phase.time_nodes.add_node(phase.time_bound, [], []) + # Organize time nodes with sort() and merge() + phase.time_nodes.sort() + phase.time_nodes.merge() # Clear triggers from first time node if necessary if phase.clear: - phase.TimeNodes[0].parachutes = [] - phase.TimeNodes[0].callbacks = [] - - # print('\n\tPhase Time Nodes') - # print('\tTime Nodes Length: ', str(len(phase.TimeNodes)), ' | Time Nodes Preview: ', phase.TimeNodes[0:3]) + phase.time_nodes[0].parachutes = [] + phase.time_nodes[0].callbacks = [] # Iterate through time nodes - for node_index, node in time_iterator(phase.TimeNodes): - # print('\n\t\tCurrent Time Node') - # print('\t\tIndex: ', node_index, ' | Time Node: ', node) + for node_index, node in self.time_iterator(phase.time_nodes): # Determine time bound for this time node - node.time_bound = phase.TimeNodes[node_index + 1].t + node.time_bound = phase.time_nodes[node_index + 1].t + # NOTE: Setting the time bound and status for the phase solver, + # and updating its internal state for the next integration step. phase.solver.t_bound = node.time_bound phase.solver._lsoda_solver._integrator.rwork[0] = phase.solver.t_bound phase.solver._lsoda_solver._integrator.call_args[4] = ( @@ -695,6 +693,7 @@ def __init__( phase.solver.status = "running" # Feed required parachute and discrete controller triggers + # TODO: parachutes should be moved to controllers for callback in node.callbacks: callback(self) @@ -703,28 +702,21 @@ def __init__( for parachute in node.parachutes: # Calculate and save pressure signal - pressure = self.env.pressure.get_value_opt(self.y_sol[2]) - parachute.clean_pressure_signal.append([node.t, pressure]) - # Calculate and save noise - noise = parachute.noise_function() - parachute.noise_signal.append([node.t, noise]) - parachute.noisy_pressure_signal.append([node.t, pressure + noise]) - # Gets height above ground level considering noise - hAGL = ( - self.env.barometric_height(pressure + noise) - - self.env.elevation + noisy_pressure, height_above_ground_level = ( + self.__calculate_and_save_pressure_signals( + parachute, node.t, self.y_sol[2] + ) ) - if parachute.triggerfunc(pressure + noise, hAGL, self.y_sol): - # print('\nEVENT DETECTED') - # print('Parachute Triggered') - # print('Name: ', parachute.name, ' | Lag: ', parachute.lag) + if parachute.triggerfunc( + noisy_pressure, height_above_ground_level, self.y_sol + ): # Remove parachute from flight parachutes self.parachutes.remove(parachute) # Create flight phase for time after detection and before inflation # Must only be created if parachute has any lag i = 1 if parachute.lag != 0: - self.FlightPhases.add_phase( + self.flight_phases.add_phase( node.t, phase.derivative, clear=True, @@ -737,7 +729,7 @@ def __init__( self, "parachute_cd_s", parachute_cd_s ) ] - self.FlightPhases.add_phase( + self.flight_phases.add_phase( node.t + parachute.lag, self.u_dot_parachute, callbacks, @@ -745,37 +737,24 @@ def __init__( index=phase_index + i, ) # Prepare to leave loops and start new flight phase - phase.TimeNodes.flush_after(node_index) - phase.TimeNodes.add_node(self.t, [], []) + phase.time_nodes.flush_after(node_index) + phase.time_nodes.add_node(self.t, [], []) phase.solver.status = "finished" # Save parachute event self.parachute_events.append([self.t, parachute]) # Step through simulation while phase.solver.status == "running": - # Step + # Execute solver step, log solution and function evaluations phase.solver.step() - # Save step result self.solution += [[phase.solver.t, *phase.solver.y]] - # Step step metrics - self.function_evaluations_per_time_step.append( - phase.solver.nfev - self.function_evaluations[-1] - ) self.function_evaluations.append(phase.solver.nfev) - self.time_steps.append(phase.solver.step_size) + # Update time and state self.t = phase.solver.t self.y_sol = phase.solver.y if verbose: - print( - "Current Simulation Time: {:3.4f} s".format(self.t), - end="\r", - ) - # print('\n\t\t\tCurrent Step Details') - # print('\t\t\tIState: ', phase.solver._lsoda_solver._integrator.istate) - # print('\t\t\tTime: ', phase.solver.t) - # print('\t\t\tAltitude: ', phase.solver.y[2]) - # print('\t\t\tEvals: ', self.function_evaluations_per_time_step[-1]) + print(f"Current Simulation Time: {self.t:3.4f} s", end="\r") # Check for first out of rail event if len(self.out_of_rail_state) == 1 and ( @@ -784,12 +763,7 @@ def __init__( + (self.y_sol[2] - self.env.elevation) ** 2 >= self.effective_1rl**2 ): - # Rocket is out of rail # Check exactly when it went out using root finding - # States before and after - # t0 -> 0 - # print('\nEVENT DETECTED') - # print('Rocket is Out of Rail!') # Disconsider elevation self.solution[-2][3] -= self.env.elevation self.solution[-1][3] -= self.env.elevation @@ -819,24 +793,23 @@ def __init__( self.solution[-2][3] += self.env.elevation self.solution[-1][3] += self.env.elevation # Cubic Hermite interpolation (ax**3 + bx**2 + cx + d) - D = float(phase.solver.step_size) - d = float(y0) - c = float(yp0) - b = float((3 * y1 - yp1 * D - 2 * c * D - 3 * d) / (D**2)) - a = float(-(2 * y1 - yp1 * D - c * D - 2 * d) / (D**3)) + 1e-5 + a, b, c, d = calculate_cubic_hermite_coefficients( + 0, + float(phase.solver.step_size), + y0, + yp0, + y1, + yp1, + ) + a += 1e-5 # TODO: why?? # Find roots - d0 = b**2 - 3 * a * c - d1 = 2 * b**3 - 9 * a * b * c + 27 * d * a**2 - c1 = ((d1 + (d1**2 - 4 * d0**3) ** (0.5)) / 2) ** (1 / 3) - t_roots = [] - for k in [0, 1, 2]: - c2 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** k - t_roots.append(-(1 / (3 * a)) * (b + c2 + d0 / c2)) + t_roots = find_roots_cubic_function(a, b, c, d) # Find correct root - valid_t_root = [] - for t_root in t_roots: - if 0 < t_root.real < t1 and abs(t_root.imag) < 0.001: - valid_t_root.append(t_root.real) + valid_t_root = [ + t_root.real + for t_root in t_roots + if 0 < t_root.real < t1 and abs(t_root.imag) < 0.001 + ] if len(valid_t_root) > 1: raise ValueError( "Multiple roots found when solving for rail exit time." @@ -853,31 +826,23 @@ def __init__( self.out_of_rail_time = self.t self.out_of_rail_time_index = len(self.solution) - 1 self.out_of_rail_state = self.y_sol - self.out_of_rail_velocity = ( - self.y_sol[3] ** 2 + self.y_sol[4] ** 2 + self.y_sol[5] ** 2 - ) ** (0.5) # Create new flight phase - self.FlightPhases.add_phase( + self.flight_phases.add_phase( self.t, self.u_dot_generalized, index=phase_index + 1, ) # Prepare to leave loops and start new flight phase - phase.TimeNodes.flush_after(node_index) - phase.TimeNodes.add_node(self.t, [], []) + phase.time_nodes.flush_after(node_index) + phase.time_nodes.add_node(self.t, [], []) phase.solver.status = "finished" # Check for apogee event if len(self.apogee_state) == 1 and self.y_sol[5] < 0: - # print('\nPASSIVE EVENT DETECTED') - # print('Rocket Has Reached Apogee!') - # Apogee reported # Assume linear vz(t) to detect when vz = 0 - vz0 = self.solution[-2][6] - t0 = self.solution[-2][0] - vz1 = self.solution[-1][6] - t1 = self.solution[-1][0] - t_root = -(t1 - t0) * vz0 / (vz1 - vz0) + t0 + t0, vz0 = self.solution[-2][0], self.solution[-2][6] + t1, vz1 = self.solution[-1][0], self.solution[-1][6] + t_root = find_root_linear_interpolation(t0, t1, vz0, vz1, 0) # Fetch state at t_root interpolator = phase.solver.dense_output() self.apogee_state = interpolator(t_root) @@ -886,88 +851,66 @@ def __init__( self.apogee_x = self.apogee_state[0] self.apogee_y = self.apogee_state[1] self.apogee = self.apogee_state[2] + if self.terminate_on_apogee: - # print('Terminate on Apogee Activated!') - self.t = t_root - self.t_final = self.t - self.state = self.apogee_state + self.t = self.t_final = t_root # Roll back solution - self.solution[-1] = [self.t, *self.state] + self.solution[-1] = [self.t, *self.apogee_state] # Set last flight phase - self.FlightPhases.flush_after(phase_index) - self.FlightPhases.add_phase(self.t) + self.flight_phases.flush_after(phase_index) + self.flight_phases.add_phase(self.t) # Prepare to leave loops and start new flight phase - phase.TimeNodes.flush_after(node_index) - phase.TimeNodes.add_node(self.t, [], []) + phase.time_nodes.flush_after(node_index) + phase.time_nodes.add_node(self.t, [], []) phase.solver.status = "finished" # Check for impact event if self.y_sol[2] < self.env.elevation: - # print('\nPASSIVE EVENT DETECTED') - # print('Rocket Has Reached Ground!') - # Impact reported - # Check exactly when it went out using root finding - # States before and after - # t0 -> 0 - # Disconsider elevation - self.solution[-2][3] -= self.env.elevation - self.solution[-1][3] -= self.env.elevation - # Get points - y0 = self.solution[-2][3] - yp0 = self.solution[-2][6] - t1 = self.solution[-1][0] - self.solution[-2][0] - y1 = self.solution[-1][3] - yp1 = self.solution[-1][6] - # Put elevation back - self.solution[-2][3] += self.env.elevation - self.solution[-1][3] += self.env.elevation + # Check exactly when it happened using root finding # Cubic Hermite interpolation (ax**3 + bx**2 + cx + d) - D = float(phase.solver.step_size) - d = float(y0) - c = float(yp0) - b = float((3 * y1 - yp1 * D - 2 * c * D - 3 * d) / (D**2)) - a = float(-(2 * y1 - yp1 * D - c * D - 2 * d) / (D**3)) + a, b, c, d = calculate_cubic_hermite_coefficients( + x0=0, # t0 + x1=float(phase.solver.step_size), # t1 - t0 + y0=float(self.solution[-2][3] - self.env.elevation), # z0 + yp0=float(self.solution[-2][6]), # vz0 + y1=float(self.solution[-1][3] - self.env.elevation), # z1 + yp1=float(self.solution[-1][6]), # vz1 + ) # Find roots - d0 = b**2 - 3 * a * c - d1 = 2 * b**3 - 9 * a * b * c + 27 * d * a**2 - c1 = ((d1 + (d1**2 - 4 * d0**3) ** (0.5)) / 2) ** (1 / 3) - t_roots = [] - for k in [0, 1, 2]: - c2 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** k - t_roots.append(-(1 / (3 * a)) * (b + c2 + d0 / c2)) + t_roots = find_roots_cubic_function(a, b, c, d) # Find correct root - valid_t_root = [] - for t_root in t_roots: - if 0 < t_root.real < t1 and abs(t_root.imag) < 0.001: - valid_t_root.append(t_root.real) + t1 = self.solution[-1][0] - self.solution[-2][0] + valid_t_root = [ + t_root.real + for t_root in t_roots + if abs(t_root.imag) < 0.001 and 0 < t_root.real < t1 + ] if len(valid_t_root) > 1: raise ValueError( "Multiple roots found when solving for impact time." ) # Determine impact state at t_root - self.t = valid_t_root[0] + self.solution[-2][0] + self.t = self.t_final = valid_t_root[0] + self.solution[-2][0] interpolator = phase.solver.dense_output() - self.y_sol = interpolator(self.t) + self.y_sol = self.impact_state = interpolator(self.t) # Roll back solution self.solution[-1] = [self.t, *self.y_sol] # Save impact state - self.impact_state = self.y_sol self.x_impact = self.impact_state[0] self.y_impact = self.impact_state[1] self.z_impact = self.impact_state[2] self.impact_velocity = self.impact_state[5] - self.t_final = self.t # Set last flight phase - self.FlightPhases.flush_after(phase_index) - self.FlightPhases.add_phase(self.t) + self.flight_phases.flush_after(phase_index) + self.flight_phases.add_phase(self.t) # Prepare to leave loops and start new flight phase - phase.TimeNodes.flush_after(node_index) - phase.TimeNodes.add_node(self.t, [], []) + phase.time_nodes.flush_after(node_index) + phase.time_nodes.add_node(self.t, [], []) phase.solver.status = "finished" # List and feed overshootable time nodes if self.time_overshoot: # Initialize phase overshootable time nodes - overshootable_nodes = TimeNodes() + overshootable_nodes = self.TimeNodes() # Add overshootable parachute time nodes overshootable_nodes.add_parachutes( self.parachutes, self.solution[-2][0], self.t @@ -975,65 +918,46 @@ def __init__( # Add last time node (always skipped) overshootable_nodes.add_node(self.t, [], []) if len(overshootable_nodes) > 1: - # Sort overshootable time nodes + # Sort and merge equal overshootable time nodes overshootable_nodes.sort() - # Merge equal overshootable time nodes overshootable_nodes.merge() # Clear if necessary if overshootable_nodes[0].t == phase.t and phase.clear: overshootable_nodes[0].parachutes = [] overshootable_nodes[0].callbacks = [] - # print('\n\t\t\t\tOvershootable Time Nodes') - # print('\t\t\t\tInterval: ', self.solution[-2][0], self.t) - # print('\t\t\t\tOvershootable Nodes Length: ', str(len(overshootable_nodes)), ' | Overshootable Nodes: ', overshootable_nodes) # Feed overshootable time nodes trigger interpolator = phase.solver.dense_output() for ( overshootable_index, overshootable_node, - ) in time_iterator(overshootable_nodes): - # print('\n\t\t\t\tCurrent Overshootable Node') - # print('\t\t\t\tIndex: ', overshootable_index, ' | Overshootable Node: ', overshootable_node) + ) in self.time_iterator(overshootable_nodes): # Calculate state at node time - overshootable_node.y = interpolator( + overshootable_node.y_sol = interpolator( overshootable_node.t ) - # Calculate and save pressure signal - pressure = self.env.pressure.get_value_opt( - overshootable_node.y[2] - ) for parachute in overshootable_node.parachutes: - # Save pressure signal - parachute.clean_pressure_signal.append( - [overshootable_node.t, pressure] - ) - # Calculate and save noise - noise = parachute.noise_function() - parachute.noise_signal.append( - [overshootable_node.t, noise] - ) - parachute.noisy_pressure_signal.append( - [overshootable_node.t, pressure + noise] - ) - # Gets height above ground level considering noise - hAGL = ( - self.env.barometric_height(pressure + noise) - - self.env.elevation + # Calculate and save pressure signal + noisy_pressure, height_above_ground_level = ( + self.__calculate_and_save_pressure_signals( + parachute, + overshootable_node.t, + overshootable_node.y_sol[2], + ) ) + # Check for parachute trigger if parachute.triggerfunc( - pressure + noise, hAGL, overshootable_node.y + noisy_pressure, + height_above_ground_level, + overshootable_node.y_sol, ): - # print('\nEVENT DETECTED') - # print('Parachute Triggered') - # print('Name: ', parachute.name, ' | Lag: ', parachute.lag) # Remove parachute from flight parachutes self.parachutes.remove(parachute) # Create flight phase for time after detection and before inflation # Must only be created if parachute has any lag i = 1 if parachute.lag != 0: - self.FlightPhases.add_phase( + self.flight_phases.add_phase( overshootable_node.t, phase.derivative, clear=True, @@ -1046,7 +970,7 @@ def __init__( self, "parachute_cd_s", parachute_cd_s ) ] - self.FlightPhases.add_phase( + self.flight_phases.add_phase( overshootable_node.t + parachute.lag, self.u_dot_parachute, callbacks, @@ -1055,17 +979,17 @@ def __init__( ) # Rollback history self.t = overshootable_node.t - self.y_sol = overshootable_node.y + self.y_sol = overshootable_node.y_sol self.solution[-1] = [ overshootable_node.t, - *overshootable_node.y, + *overshootable_node.y_sol, ] # Prepare to leave loops and start new flight phase overshootable_nodes.flush_after( overshootable_index ) - phase.TimeNodes.flush_after(node_index) - phase.TimeNodes.add_node(self.t, [], []) + phase.time_nodes.flush_after(node_index) + phase.time_nodes.add_node(self.t, [], []) phase.solver.status = "finished" # Save parachute event self.parachute_events.append( @@ -1073,30 +997,52 @@ def __init__( ) self.t_final = self.t - self._calculate_pressure_signal() + self.__transform_pressure_signals_lists_to_functions() if verbose: - print("Simulation Completed at Time: {:3.4f} s".format(self.t)) + print(f"\n>>> Simulation Completed at Time: {self.t:3.4f} s") - def __init_post_process_variables(self): - """Initialize post-process variables.""" - # Initialize all variables calculated after initialization. - # Important to do so that MATLAB® can access them - self._drift = Function(0) - self._bearing = Function(0) - self._latitude = Function(0) - self._longitude = Function(0) + def __calculate_and_save_pressure_signals(self, parachute, t, z): + """Gets noise and pressure signals and saves them in the parachute + object given the current time and altitude. + + Parameters + ---------- + parachute : Parachute + The parachute object to calculate signals for. + t : float + The current time in seconds. + z : float + The altitude above sea level in meters. + + Returns + ------- + tuple[float, float] + The noisy pressure and height above ground level. + """ + # Calculate pressure and noise + pressure = self.env.pressure.get_value_opt(z) + noise = parachute.noise_function() + noisy_pressure = pressure + noise + + # Stores in the parachute object + parachute.clean_pressure_signal.append([t, pressure]) + parachute.noise_signal.append([t, noise]) + + # Gets height above ground level considering noise + height_above_ground_level = ( + self.env.barometric_height.get_value_opt(noisy_pressure) + - self.env.elevation + ) + + return noisy_pressure, height_above_ground_level def __init_solution_monitors(self): # Initialize solution monitors self.out_of_rail_time = 0 self.out_of_rail_time_index = 0 self.out_of_rail_state = np.array([0]) - self.out_of_rail_velocity = 0 self.apogee_state = np.array([0]) self.apogee_time = 0 - self.apogee_x = 0 - self.apogee_y = 0 - self.apogee = 0 self.x_impact = 0 self.y_impact = 0 self.impact_velocity = 0 @@ -1104,8 +1050,6 @@ def __init_solution_monitors(self): self.parachute_events = [] self.post_processed = False - return None - def __init_flight_state(self): """Initialize flight state variables.""" if self.initial_solution is None: @@ -1162,8 +1106,6 @@ def __init_flight_state(self): def __init_solver_monitors(self): # Initialize solver monitors self.function_evaluations = [] - self.function_evaluations_per_time_step = [] - self.time_steps = [] # Initialize solution state self.solution = [] self.__init_flight_state() @@ -1176,11 +1118,7 @@ def __init_solver_monitors(self): def __init_equations_of_motion(self): """Initialize equations of motion.""" if self.equations_of_motion == "solid_propulsion": - self.u_dot_generalized = self.u_dot - - def __init_equations_of_motion(self): - """Initialize equations of motion.""" - if self.equations_of_motion == "solid_propulsion": + # NOTE: The u_dot is faster, but only works for solid propulsion self.u_dot_generalized = self.u_dot @cached_property @@ -1216,25 +1154,34 @@ def effective_2rl(self): @cached_property def frontal_surface_wind(self): - # Surface wind magnitude in the frontal direction at the rail's elevation - wind_u = self.env.wind_velocity_x(self.env.elevation) - wind_v = self.env.wind_velocity_y(self.env.elevation) + """Frontal wind velocity at the surface level. The frontal wind is + defined as the wind blowing in the direction of the rocket's heading. + + Returns + ------- + float + Wind velocity in the frontal direction at the surface level. + """ + wind_u = self.env.wind_velocity_x.get_value_opt(self.env.elevation) + wind_v = self.env.wind_velocity_y.get_value_opt(self.env.elevation) heading_rad = self.heading * np.pi / 180 - frontal_surface_wind = wind_u * np.sin(heading_rad) + wind_v * np.cos( - heading_rad - ) - return frontal_surface_wind + return wind_u * np.sin(heading_rad) + wind_v * np.cos(heading_rad) @cached_property def lateral_surface_wind(self): - # Surface wind magnitude in the lateral direction at the rail's elevation - wind_u = self.env.wind_velocity_x(self.env.elevation) - wind_v = self.env.wind_velocity_y(self.env.elevation) + """Lateral wind velocity at the surface level. The lateral wind is + defined as the wind blowing perpendicular to the rocket's heading. + + Returns + ------- + float + Wind velocity in the lateral direction at the surface level. + """ + wind_u = self.env.wind_velocity_x.get_value_opt(self.env.elevation) + wind_v = self.env.wind_velocity_y.get_value_opt(self.env.elevation) heading_rad = self.heading * np.pi / 180 - lateral_surface_wind = -wind_u * np.cos(heading_rad) + wind_v * np.sin( - heading_rad - ) - return lateral_surface_wind + + return -wind_u * np.cos(heading_rad) + wind_v * np.sin(heading_rad) def udot_rail1(self, t, u, post_processing=False): """Calculates derivative of u state vector with respect to time @@ -1258,11 +1205,6 @@ def udot_rail1(self, t, u, post_processing=False): e0dot, e1dot, e2dot, e3dot, alpha1, alpha2, alpha3]. """ - # Check if post processing mode is on - if post_processing: - # Use u_dot post processing code - return self.u_dot_generalized(t, u, True) - # Retrieve integration data x, y, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u @@ -1285,7 +1227,9 @@ def udot_rail1(self, t, u, post_processing=False): R3 = -0.5 * rho * (free_stream_speed**2) * self.rocket.area * (drag_coeff) # Calculate Linear acceleration - a3 = (R3 + thrust) / M - (e0**2 - e1**2 - e2**2 + e3**2) * self.env.gravity(z) + a3 = (R3 + thrust) / M - ( + e0**2 - e1**2 - e2**2 + e3**2 + ) * self.env.gravity.get_value_opt(z) if a3 > 0: ax = 2 * (e1 * e3 + e0 * e2) * a3 ay = 2 * (e2 * e3 - e0 * e1) * a3 @@ -1293,6 +1237,15 @@ def udot_rail1(self, t, u, post_processing=False): else: ax, ay, az = 0, 0, 0 + if post_processing: + self.u_dot_generalized(t, u, post_processing=True) + self.__post_processed_variables[-1][1] = ax + self.__post_processed_variables[-1][2] = ay + self.__post_processed_variables[-1][3] = az + self.__post_processed_variables[-1][4] = 0 + self.__post_processed_variables[-1][5] = 0 + self.__post_processed_variables[-1][6] = 0 + return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0] def udot_rail2(self, t, u, post_processing=False): @@ -1345,8 +1298,7 @@ def u_dot(self, t, u, post_processing=False): # Retrieve integration data x, y, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u # Determine lift force and moment - R1, R2 = 0, 0 - M1, M2, M3 = 0, 0, 0 + R1, R2, M1, M2, M3 = 0, 0, 0, 0, 0 # Determine current behavior if t < self.rocket.motor.burn_out_time: # Motor burning @@ -1366,7 +1318,6 @@ def u_dot(self, t, u, post_processing=False): M2 -= self.rocket.thrust_eccentricity_y * thrust else: # Motor stopped - # Retrieve important motor quantities # Inertias Tz, Ti, Tzdot, Tidot = 0, 0, 0, 0 # Mass @@ -1386,7 +1337,7 @@ def u_dot(self, t, u, post_processing=False): # b = -self.rocket.distance_rocket_propellant b = ( -( - self.rocket.center_of_propellant_position(0) + self.rocket.center_of_propellant_position.get_value_opt(0) - self.rocket.center_of_dry_mass_position ) * self.rocket._csys @@ -1410,8 +1361,6 @@ def u_dot(self, t, u, post_processing=False): a33 = 1 - 2 * (e1**2 + e2**2) # Transformation matrix: (123) -> (XYZ) K = [[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]] - # Transformation matrix: (XYZ) -> (123) or K transpose - Kt = [[a11, a21, a31], [a12, a22, a32], [a13, a23, a33]] # Calculate Forces and Moments # Get freestream speed @@ -1432,7 +1381,7 @@ def u_dot(self, t, u, post_processing=False): R3 = -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff for air_brakes in self.rocket.air_brakes: if air_brakes.deployment_level > 0: - air_brakes_cd = air_brakes.drag_coefficient( + air_brakes_cd = air_brakes.drag_coefficient.get_value_opt( air_brakes.deployment_level, free_stream_mach ) air_brakes_force = ( @@ -1446,7 +1395,6 @@ def u_dot(self, t, u, post_processing=False): R3 = air_brakes_force # Substitutes rocket drag coefficient else: R3 += air_brakes_force - # R3 += self.__computeDragForce(z, Vector(vx, vy, vz)) # Off center moment M1 += self.rocket.cp_eccentricity_y * R3 M2 -= self.rocket.cp_eccentricity_x * R3 @@ -1487,7 +1435,9 @@ def u_dot(self, t, u, post_processing=False): comp_stream_vz_bn = comp_stream_vz_b / comp_stream_speed if -1 * comp_stream_vz_bn < 1: comp_attack_angle = np.arccos(-comp_stream_vz_bn) - c_lift = aero_surface.cl(comp_attack_angle, free_stream_mach) + c_lift = aero_surface.cl.get_value_opt( + comp_attack_angle, free_stream_mach + ) # component lift force magnitude comp_lift = ( 0.5 * rho * (comp_stream_speed**2) * reference_area * c_lift @@ -1510,14 +1460,14 @@ def u_dot(self, t, u, post_processing=False): * reference_area * 2 * surface_radius - * clf_delta(free_stream_mach) + * clf_delta.get_value_opt(free_stream_mach) * cant_angle_rad ) M3d = ( (1 / 2 * rho * free_stream_speed) * reference_area * (2 * surface_radius) ** 2 - * cld_omega(free_stream_mach) + * cld_omega.get_value_opt(free_stream_mach) * omega3 / 2 ) @@ -1562,7 +1512,7 @@ def u_dot(self, t, u, post_processing=False): (R3 - b * Mt * (alpha2 - omega1 * omega3) + thrust) / M, ] ax, ay, az = np.dot(K, L) - az -= self.env.gravity(z) # Include gravity + az -= self.env.gravity.get_value_opt(z) # Include gravity # Create u_dot u_dot = [ @@ -1582,27 +1532,8 @@ def u_dot(self, t, u, post_processing=False): ] if post_processing: - # Dynamics variables - self.R1_list.append([t, R1]) - self.R2_list.append([t, R2]) - self.R3_list.append([t, R3]) - self.M1_list.append([t, M1]) - self.M2_list.append([t, M2]) - self.M3_list.append([t, M3]) - # Atmospheric Conditions - self.wind_velocity_x_list.append( - [t, self.env.wind_velocity_x.get_value_opt(z)] - ) - self.wind_velocity_y_list.append( - [t, self.env.wind_velocity_y.get_value_opt(z)] - ) - self.density_list.append([t, self.env.density.get_value_opt(z)]) - self.dynamic_viscosity_list.append( - [t, self.env.dynamic_viscosity.get_value_opt(z)] - ) - self.pressure_list.append([t, self.env.pressure.get_value_opt(z)]) - self.speed_of_sound_list.append( - [t, self.env.speed_of_sound.get_value_opt(z)] + self.__post_processed_variables.append( + [t, ax, ay, az, alpha1, alpha2, alpha3, R1, R2, R3, M1, M2, M3] ) return u_dot @@ -1642,8 +1573,8 @@ def u_dot_generalized(self, t, u, post_processing=False): # Retrieve necessary quantities rho = self.env.density.get_value_opt(z) total_mass = self.rocket.total_mass.get_value_opt(t) - total_mass_dot = self.rocket.total_mass.differentiate(t) - total_mass_ddot = self.rocket.total_mass.differentiate(t, order=2) + total_mass_dot = self.rocket.total_mass_flow_rate.get_value_opt(t) + total_mass_ddot = self.rocket.total_mass_flow_rate.differentiate(t) ## CM position vector and time derivatives relative to CDM in body frame r_CM_z = ( -1 @@ -1666,11 +1597,8 @@ def u_dot_generalized(self, t, u, post_processing=False): * self.rocket._csys ) S_noz_33 = 0.5 * self.rocket.motor.nozzle_radius**2 - S_noz_11 = 0.5 * S_noz_33 + 0.25 * r_NOZ**2 - S_noz_22 = S_noz_11 - S_noz_12 = 0 - S_noz_13 = 0 - S_noz_23 = 0 + S_noz_11 = S_noz_22 = 0.5 * S_noz_33 + 0.25 * r_NOZ**2 + S_noz_12, S_noz_13, S_noz_23 = 0, 0, 0 S_nozzle = Matrix( [ [S_noz_11, S_noz_12, S_noz_13], @@ -1723,7 +1651,9 @@ def u_dot_generalized(self, t, u, post_processing=False): wind_velocity_y = self.env.wind_velocity_y.get_value_opt(z) wind_velocity = Vector([wind_velocity_x, wind_velocity_y, 0]) free_stream_speed = abs((wind_velocity - Vector(v))) - free_stream_mach = free_stream_speed / self.env.speed_of_sound.get_value_opt(z) + speed_of_sound = self.env.speed_of_sound.get_value_opt(z) + free_stream_mach = free_stream_speed / speed_of_sound + if t < self.rocket.motor.burn_out_time: drag_coeff = self.rocket.power_on_drag.get_value_opt(free_stream_mach) else: @@ -1731,7 +1661,7 @@ def u_dot_generalized(self, t, u, post_processing=False): R3 += -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff for air_brakes in self.rocket.air_brakes: if air_brakes.deployment_level > 0: - air_brakes_cd = air_brakes.drag_coefficient( + air_brakes_cd = air_brakes.drag_coefficient.get_value_opt( air_brakes.deployment_level, free_stream_mach ) air_brakes_force = ( @@ -1770,18 +1700,18 @@ def u_dot_generalized(self, t, u, post_processing=False): comp_stream_velocity = comp_wind_vb - comp_vb comp_stream_vx_b, comp_stream_vy_b, comp_stream_vz_b = comp_stream_velocity comp_stream_speed = abs(comp_stream_velocity) - comp_stream_mach = ( - comp_stream_speed / self.env.speed_of_sound.get_value_opt(z) - ) + comp_stream_mach = comp_stream_speed / speed_of_sound # Component attack angle and lift force comp_attack_angle = 0 comp_lift, comp_lift_xb, comp_lift_yb = 0, 0, 0 - if comp_stream_vx_b**2 + comp_stream_vy_b**2 != 0: + if comp_stream_vx_b**2 + comp_stream_vy_b**2 != 0: # TODO: maybe try/except # Normalize component stream velocity in body frame comp_stream_vz_bn = comp_stream_vz_b / comp_stream_speed if -1 * comp_stream_vz_bn < 1: comp_attack_angle = np.arccos(-comp_stream_vz_bn) - c_lift = aero_surface.cl(comp_attack_angle, comp_stream_mach) + c_lift = aero_surface.cl.get_value_opt( + comp_attack_angle, comp_stream_mach + ) # Component lift force magnitude comp_lift = ( 0.5 * rho * (comp_stream_speed**2) * reference_area * c_lift @@ -1804,31 +1734,29 @@ def u_dot_generalized(self, t, u, post_processing=False): * reference_area * 2 * surface_radius - * clf_delta(comp_stream_mach) + * clf_delta.get_value_opt(comp_stream_mach) * cant_angle_rad ) M3d = ( (1 / 2 * rho * comp_stream_speed) * reference_area * (2 * surface_radius) ** 2 - * cld_omega(comp_stream_mach) + * cld_omega.get_value_opt(comp_stream_mach) * omega3 / 2 ) M3 += M3f - M3d except AttributeError: pass - weightB = Kt @ Vector([0, 0, -total_mass * self.env.gravity(z)]) + weightB = Kt @ Vector([0, 0, -total_mass * self.env.gravity.get_value_opt(z)]) T00 = total_mass * r_CM - T03 = ( - 2 * total_mass_dot * (Vector([0, 0, r_NOZ]) - r_CM) - - 2 * total_mass * r_CM_dot - ) + r_NOZ_vector = Vector([0, 0, r_NOZ]) + T03 = 2 * total_mass_dot * (r_NOZ_vector - r_CM) - 2 * total_mass * r_CM_dot T04 = ( - self.rocket.motor.thrust(t) * Vector([0, 0, 1]) + self.rocket.motor.thrust.get_value_opt(t) * Vector([0, 0, 1]) - total_mass * r_CM_ddot - 2 * total_mass_dot * r_CM_dot - + total_mass_ddot * (Vector([0, 0, r_NOZ]) - r_CM) + + total_mass_ddot * (r_NOZ_vector - r_CM) ) T05 = total_mass_dot * S_nozzle - I_dot @@ -1857,27 +1785,8 @@ def u_dot_generalized(self, t, u, post_processing=False): u_dot = [*r_dot, *v_dot, *e_dot, *w_dot] if post_processing: - # Dynamics variables - self.R1_list.append([t, R1]) - self.R2_list.append([t, R2]) - self.R3_list.append([t, R3]) - self.M1_list.append([t, M1]) - self.M2_list.append([t, M2]) - self.M3_list.append([t, M3]) - # Atmospheric Conditions - self.wind_velocity_x_list.append( - [t, self.env.wind_velocity_x.get_value_opt(z)] - ) - self.wind_velocity_y_list.append( - [t, self.env.wind_velocity_y.get_value_opt(z)] - ) - self.density_list.append([t, self.env.density.get_value_opt(z)]) - self.dynamic_viscosity_list.append( - [t, self.env.dynamic_viscosity.get_value_opt(z)] - ) - self.pressure_list.append([t, self.env.pressure.get_value_opt(z)]) - self.speed_of_sound_list.append( - [t, self.env.speed_of_sound.get_value_opt(z)] + self.__post_processed_variables.append( + [t, *v_dot, *w_dot, R1, R2, R3, M1, M2, M3] ) return u_dot @@ -1912,7 +1821,7 @@ def u_dot_parachute(self, t, u, post_processing=False): rho = self.env.density.get_value_opt(u[2]) to = 1.2 ma = ka * rho * (4 / 3) * np.pi * R**3 - mp = self.rocket.mass + mp = self.rocket.dry_mass eta = 1 Rdot = (6 * R * (1 - eta) / (1.2**6)) * ( (1 - eta) * t**5 + eta * (to**3) * (t**2) @@ -1941,20 +1850,9 @@ def u_dot_parachute(self, t, u, post_processing=False): az = (Dz - 9.8 * mp) / (mp + ma) if post_processing: - # Dynamics variables - self.R1_list.append([t, Dx]) - self.R2_list.append([t, Dy]) - self.R3_list.append([t, Dz]) - self.M1_list.append([t, 0]) - self.M2_list.append([t, 0]) - self.M3_list.append([t, 0]) - # Atmospheric Conditions - self.wind_velocity_x_list.append([t, self.env.wind_velocity_x(z)]) - self.wind_velocity_y_list.append([t, self.env.wind_velocity_y(z)]) - self.density_list.append([t, self.env.density(z)]) - self.dynamic_viscosity_list.append([t, self.env.dynamic_viscosity(z)]) - self.pressure_list.append([t, self.env.pressure(z)]) - self.speed_of_sound_list.append([t, self.env.speed_of_sound(z)]) + self.__post_processed_variables.append( + [t, ax, ay, az, 0, 0, 0, Dx, Dy, Dz, 0, 0, 0] + ) return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0] @@ -1963,11 +1861,29 @@ def solution_array(self): """Returns solution array of the rocket flight.""" return np.array(self.solution) + @property + def function_evaluations_per_time_step(self): + """Get the number of function evaluations per time step. This method + calculates the difference between consecutive function evaluations + during numerical integration and returns it as a list. + + Returns + ------- + list + The list of differences in function evaluations per time step. + """ + return np.diff(self.function_evaluations).tolist() + @cached_property def time(self): """Returns time array from solution.""" return self.solution_array[:, 0] + @cached_property + def time_steps(self): + """Returns time step array.""" + return np.diff(self.time) + def get_solution_at_time(self, t, atol=1e-3): """Returns the solution state vector at a given time. If the time is not found in the solution, the closest time is used and a warning is @@ -2075,51 +1991,51 @@ def w3(self): @funcify_method("Time (s)", "Ax (m/s²)", "spline", "zero") def ax(self): """Rocket x acceleration as a Function of time.""" - return self.retrieve_acceleration_arrays[0] + return self.__evaluate_post_process[:, [0, 1]] @funcify_method("Time (s)", "Ay (m/s²)", "spline", "zero") def ay(self): """Rocket y acceleration as a Function of time.""" - return self.retrieve_acceleration_arrays[1] + return self.__evaluate_post_process[:, [0, 2]] @funcify_method("Time (s)", "Az (m/s²)", "spline", "zero") def az(self): """Rocket z acceleration as a Function of time.""" - return self.retrieve_acceleration_arrays[2] + return self.__evaluate_post_process[:, [0, 3]] @funcify_method("Time (s)", "α1 (rad/s²)", "spline", "zero") def alpha1(self): """Rocket angular acceleration α1 as a Function of time.""" - return self.retrieve_acceleration_arrays[3] + return self.__evaluate_post_process[:, [0, 4]] @funcify_method("Time (s)", "α2 (rad/s²)", "spline", "zero") def alpha2(self): """Rocket angular acceleration α2 as a Function of time.""" - return self.retrieve_acceleration_arrays[4] + return self.__evaluate_post_process[:, [0, 5]] @funcify_method("Time (s)", "α3 (rad/s²)", "spline", "zero") def alpha3(self): """Rocket angular acceleration α3 as a Function of time.""" - return self.retrieve_acceleration_arrays[5] + return self.__evaluate_post_process[:, [0, 6]] # Process third type of outputs - Temporary values @funcify_method("Time (s)", "R1 (N)", "spline", "zero") def R1(self): """Aerodynamic force along the first axis that is perpendicular to the rocket's axis of symmetry as a Function of time.""" - return self.retrieve_temporary_values_arrays[0] + return self.__evaluate_post_process[:, [0, 7]] @funcify_method("Time (s)", "R2 (N)", "spline", "zero") def R2(self): """Aerodynamic force along the second axis that is perpendicular to the rocket's axis of symmetry as a Function of time.""" - return self.retrieve_temporary_values_arrays[1] + return self.__evaluate_post_process[:, [0, 8]] @funcify_method("Time (s)", "R3 (N)", "spline", "zero") def R3(self): """Aerodynamic force along the rocket's axis of symmetry as a Function of time.""" - return self.retrieve_temporary_values_arrays[2] + return self.__evaluate_post_process[:, [0, 9]] @funcify_method("Time (s)", "M1 (Nm)", "spline", "zero") def M1(self): @@ -2127,54 +2043,51 @@ def M1(self): perpendicular to the rocket's axis of symmetry as a Function of time. """ - return self.retrieve_temporary_values_arrays[3] + return self.__evaluate_post_process[:, [0, 10]] @funcify_method("Time (s)", "M2 (Nm)", "spline", "zero") def M2(self): """Aerodynamic bending moment in the same direction as the axis that is perpendicular to the rocket's axis of symmetry as a Function of time.""" - return self.retrieve_temporary_values_arrays[4] + return self.__evaluate_post_process[:, [0, 11]] @funcify_method("Time (s)", "M3 (Nm)", "spline", "zero") def M3(self): """Aerodynamic bending moment in the same direction as the rocket's axis of symmetry as a Function of time.""" - return self.retrieve_temporary_values_arrays[5] + return self.__evaluate_post_process[:, [0, 12]] @funcify_method("Time (s)", "Pressure (Pa)", "spline", "constant") def pressure(self): """Air pressure felt by the rocket as a Function of time.""" - return self.retrieve_temporary_values_arrays[6] + return [(t, self.env.pressure.get_value_opt(z)) for t, z in self.z] @funcify_method("Time (s)", "Density (kg/m³)", "spline", "constant") def density(self): """Air density felt by the rocket as a Function of time.""" - return self.retrieve_temporary_values_arrays[7] + return [(t, self.env.density.get_value_opt(z)) for t, z in self.z] @funcify_method("Time (s)", "Dynamic Viscosity (Pa s)", "spline", "constant") def dynamic_viscosity(self): """Air dynamic viscosity felt by the rocket as a Function of time.""" - return self.retrieve_temporary_values_arrays[8] + return [(t, self.env.dynamic_viscosity.get_value_opt(z)) for t, z in self.z] @funcify_method("Time (s)", "Speed of Sound (m/s)", "spline", "constant") def speed_of_sound(self): - """Speed of sound in the air felt by the rocket as a Function - of time.""" - return self.retrieve_temporary_values_arrays[9] + """Speed of sound in the air felt by the rocket as a Function of time.""" + return [(t, self.env.speed_of_sound.get_value_opt(z)) for t, z in self.z] @funcify_method("Time (s)", "Wind Velocity X (East) (m/s)", "spline", "constant") def wind_velocity_x(self): - """Wind velocity in the X direction (east) as a Function of - time.""" - return self.retrieve_temporary_values_arrays[10] + """Wind velocity in the X direction (east) as a Function of time.""" + return [(t, self.env.wind_velocity_x.get_value_opt(z)) for t, z in self.z] @funcify_method("Time (s)", "Wind Velocity Y (North) (m/s)", "spline", "constant") def wind_velocity_y(self): - """Wind velocity in the y direction (north) as a Function of - time.""" - return self.retrieve_temporary_values_arrays[11] + """Wind velocity in the y direction (north) as a Function of time.""" + return [(t, self.env.wind_velocity_y.get_value_opt(z)) for t, z in self.z] # Process fourth type of output - values calculated from previous outputs @@ -2185,6 +2098,11 @@ def speed(self): """Rocket speed, or velocity magnitude, as a Function of time.""" return (self.vx**2 + self.vy**2 + self.vz**2) ** 0.5 + @property + def out_of_rail_velocity(self): + """Velocity at which the rocket leaves the launch rail.""" + return self.speed.get_value_opt(self.out_of_rail_time) + @cached_property def max_speed_time(self): """Time at which the rocket reaches its maximum speed.""" @@ -2194,7 +2112,7 @@ def max_speed_time(self): @cached_property def max_speed(self): """Maximum speed reached by the rocket.""" - return self.speed(self.max_speed_time) + return self.speed.get_value_opt(self.max_speed_time) # Accelerations @funcify_method("Time (s)", "acceleration Magnitude (m/s²)") @@ -2220,7 +2138,7 @@ def max_acceleration_power_on_time(self): @cached_property def max_acceleration_power_on(self): """Maximum acceleration reached by the rocket during motor burn.""" - return self.acceleration(self.max_acceleration_power_on_time) + return self.acceleration.get_value_opt(self.max_acceleration_power_on_time) @cached_property def max_acceleration_power_off_time(self): @@ -2237,7 +2155,7 @@ def max_acceleration_power_off_time(self): @cached_property def max_acceleration_power_off(self): """Maximum acceleration reached by the rocket after motor burn.""" - return self.acceleration(self.max_acceleration_power_off_time) + return self.acceleration.get_value_opt(self.max_acceleration_power_off_time) @cached_property def max_acceleration_time(self): @@ -2290,13 +2208,13 @@ def attitude_angle(self): attitude_angle = (180 / np.pi) * np.arctan2( self.attitude_vector_z[:, 1], horizontal_attitude_proj[:, 1] ) - attitude_angle = np.column_stack([self.time, attitude_angle]) - return attitude_angle + return np.column_stack([self.time, attitude_angle]) # Lateral Attitude Angle @funcify_method("Time (s)", "Lateral Attitude Angle (°)") def lateral_attitude_angle(self): """Rocket lateral attitude angle as a Function of time.""" + # TODO: complex method, it should be defined elsewhere. lateral_vector_angle = (np.pi / 180) * (self.heading - 90) lateral_vector_x = np.sin(lateral_vector_angle) lateral_vector_y = np.cos(lateral_vector_angle) @@ -2352,24 +2270,17 @@ def theta(self): @funcify_method("Time (s)", "Freestream Velocity X (m/s)", "spline", "constant") def stream_velocity_x(self): """Freestream velocity X component as a Function of time.""" - stream_velocity_x = np.column_stack( - (self.time, self.wind_velocity_x[:, 1] - self.vx[:, 1]) - ) - return stream_velocity_x + return np.column_stack((self.time, self.wind_velocity_x[:, 1] - self.vx[:, 1])) @funcify_method("Time (s)", "Freestream Velocity Y (m/s)", "spline", "constant") def stream_velocity_y(self): """Freestream velocity Y component as a Function of time.""" - stream_velocity_y = np.column_stack( - (self.time, self.wind_velocity_y[:, 1] - self.vy[:, 1]) - ) - return stream_velocity_y + return np.column_stack((self.time, self.wind_velocity_y[:, 1] - self.vy[:, 1])) @funcify_method("Time (s)", "Freestream Velocity Z (m/s)", "spline", "constant") def stream_velocity_z(self): """Freestream velocity Z component as a Function of time.""" - stream_velocity_z = np.column_stack((self.time, -self.vz[:, 1])) - return stream_velocity_z + return np.column_stack((self.time, -self.vz[:, 1])) @funcify_method("Time (s)", "Freestream Speed (m/s)", "spline", "constant") def free_stream_speed(self): @@ -2385,7 +2296,7 @@ def free_stream_speed(self): @cached_property def apogee_freestream_speed(self): """Freestream speed at apogee in m/s.""" - return self.free_stream_speed(self.apogee_time) + return self.free_stream_speed.get_value_opt(self.apogee_time) # Mach Number @funcify_method("Time (s)", "Mach Number", "spline", "zero") @@ -2402,7 +2313,7 @@ def max_mach_number_time(self): @cached_property def max_mach_number(self): """Maximum Mach number.""" - return self.mach_number(self.max_mach_number_time) + return self.mach_number.get_value_opt(self.max_mach_number_time) # Stability Margin @cached_property @@ -2414,7 +2325,7 @@ def max_stability_margin_time(self): @cached_property def max_stability_margin(self): """Maximum stability margin.""" - return self.stability_margin(self.max_stability_margin_time) + return self.stability_margin.get_value_opt(self.max_stability_margin_time) @cached_property def min_stability_margin_time(self): @@ -2425,7 +2336,27 @@ def min_stability_margin_time(self): @cached_property def min_stability_margin(self): """Minimum stability margin.""" - return self.stability_margin(self.min_stability_margin_time) + return self.stability_margin.get_value_opt(self.min_stability_margin_time) + + @property + def initial_stability_margin(self): + """Stability margin at time 0. + + Returns + ------- + float + """ + return self.stability_margin.get_value_opt(self.time[0]) + + @property + def out_of_rail_stability_margin(self): + """Stability margin at the time the rocket leaves the rail. + + Returns + ------- + float + """ + return self.stability_margin.get_value_opt(self.out_of_rail_time) # Reynolds Number @funcify_method("Time (s)", "Reynolds Number", "spline", "zero") @@ -2444,7 +2375,7 @@ def max_reynolds_number_time(self): @cached_property def max_reynolds_number(self): """Maximum Reynolds number.""" - return self.reynolds_number(self.max_reynolds_number_time) + return self.reynolds_number.get_value_opt(self.max_reynolds_number_time) # Dynamic Pressure @funcify_method("Time (s)", "Dynamic Pressure (Pa)", "spline", "zero") @@ -2461,7 +2392,7 @@ def max_dynamic_pressure_time(self): @cached_property def max_dynamic_pressure(self): """Maximum dynamic pressure.""" - return self.dynamic_pressure(self.max_dynamic_pressure_time) + return self.dynamic_pressure.get_value_opt(self.max_dynamic_pressure_time) # Total Pressure @funcify_method("Time (s)", "Total Pressure (Pa)", "spline", "zero") @@ -2477,7 +2408,7 @@ def max_total_pressure_time(self): @cached_property def max_total_pressure(self): """Maximum total pressure.""" - return self.total_pressure(self.max_total_pressure_time) + return self.total_pressure.get_value_opt(self.max_total_pressure_time) # Dynamics functions and variables @@ -2520,7 +2451,7 @@ def translational_energy(self): # Redefine total_mass time grid to allow for efficient Function algebra total_mass = deepcopy(self.rocket.total_mass) total_mass.set_discrete_based_on_model(self.vz) - translational_energy = 0.5 * total_mass * (self.vx**2 + self.vy**2 + self.vz**2) + translational_energy = 0.5 * total_mass * (self.speed**2) return translational_energy @funcify_method("Time (s)", "Kinetic Energy (J)", "spline", "zero") @@ -2534,16 +2465,15 @@ def potential_energy(self): """Potential energy as a Function of time in relation to sea level.""" # Constants - GM = 3.986004418e14 + GM = 3.986004418e14 # TODO: this constant should come from Environment. # Redefine total_mass time grid to allow for efficient Function algebra total_mass = deepcopy(self.rocket.total_mass) total_mass.set_discrete_based_on_model(self.z) - potential_energy = ( + return ( GM * total_mass * (1 / (self.z + self.env.earth_radius) - 1 / self.env.earth_radius) ) - return potential_energy # Total Mechanical Energy @funcify_method("Time (s)", "Mechanical Energy (J)", "spline", "constant") @@ -2583,7 +2513,7 @@ def angle_of_attack(self): for i in self.time ] # Define freestream speed list - free_stream_speed = [self.free_stream_speed(i) for i in self.time] + free_stream_speed = [self.free_stream_speed.get_value_opt(i) for i in self.time] free_stream_speed = np.nan_to_num(free_stream_speed) # Normalize dot product @@ -2595,9 +2525,8 @@ def angle_of_attack(self): # Calculate angle of attack and convert to degrees angle_of_attack = np.rad2deg(np.arccos(dot_product_normalized)) - angle_of_attack = np.column_stack([self.time, angle_of_attack]) - return angle_of_attack + return np.column_stack([self.time, angle_of_attack]) # Frequency response and stability variables @funcify_method("Frequency (Hz)", "ω1 Fourier Amplitude", "spline", "zero") @@ -2745,6 +2674,7 @@ def latitude(self): ) return np.column_stack((self.time, latitude)) + # TODO: haversine should be defined in tools.py so we just invoke it in here. @funcify_method("Time (s)", "Longitude (°)", "linear", "constant") def longitude(self): """Rocket longitude coordinate, in degrees, as a Function of @@ -2767,142 +2697,6 @@ def longitude(self): return np.column_stack((self.time, longitude)) - @cached_property - def retrieve_acceleration_arrays(self): - """Retrieve acceleration arrays from the integration scheme - - Parameters - ---------- - - Returns - ------- - ax: list - acceleration in x direction - ay: list - acceleration in y direction - az: list - acceleration in z direction - alpha1: list - angular acceleration in x direction - alpha2: list - angular acceleration in y direction - alpha3: list - angular acceleration in z direction - """ - # Initialize acceleration arrays - ax, ay, az = [[0, 0]], [[0, 0]], [[0, 0]] - alpha1, alpha2, alpha3 = [[0, 0]], [[0, 0]], [[0, 0]] - # Go through each time step and calculate accelerations - # Get flight phases - for phase_index, phase in self.time_iterator(self.FlightPhases): - init_time = phase.t - final_time = self.FlightPhases[phase_index + 1].t - current_derivative = phase.derivative - # Call callback functions - for callback in phase.callbacks: - callback(self) - # Loop through time steps in flight phase - for step in self.solution: # Can be optimized - if init_time < step[0] <= final_time: - # Get derivatives - u_dot = current_derivative(step[0], step[1:]) - # Get accelerations - ax_value, ay_value, az_value = u_dot[3:6] - alpha1_value, alpha2_value, alpha3_value = u_dot[10:] - # Save accelerations - ax.append([step[0], ax_value]) - ay.append([step[0], ay_value]) - az.append([step[0], az_value]) - alpha1.append([step[0], alpha1_value]) - alpha2.append([step[0], alpha2_value]) - alpha3.append([step[0], alpha3_value]) - - return ax, ay, az, alpha1, alpha2, alpha3 - - @cached_property - def retrieve_temporary_values_arrays(self): - """Retrieve temporary values arrays from the integration scheme. - Currently, the following temporary values are retrieved: ``R1`` , ``R2`` - ``R3`` , ``M1`` , ``M2`` , ``M3`` , ``pressure`` , ``density`` , - ``dynamic_viscosity`` , ``speed_of_sound`` . - - Returns - ------- - self.R1_list: list - R1 values. - self.R2_list: list - R2 values. - self.R3_list: list - R3 values are the aerodynamic force values in the rocket's axis - direction. - self.M1_list: list - M1 values. - self.M2_list: list - Aerodynamic bending moment in e2 direction at each time step. - self.M3_list: list - Aerodynamic bending moment in e3 direction at each time step. - self.pressure_list: list - Air pressure at each time step. - self.density_list: list - Air density at each time step. - self.dynamic_viscosity_list: list - Dynamic viscosity at each time step. - self.speed_of_sound_list: list - Speed of sound at each time step. - self.wind_velocity_x_list: list - Wind velocity in x direction at each time step. - self.wind_velocity_y_list: list - Wind velocity in y direction at each time step. - """ - - # Initialize force and atmospheric arrays - self.R1_list = [] - self.R2_list = [] - self.R3_list = [] - self.M1_list = [] - self.M2_list = [] - self.M3_list = [] - self.pressure_list = [] - self.density_list = [] - self.dynamic_viscosity_list = [] - self.speed_of_sound_list = [] - self.wind_velocity_x_list = [] - self.wind_velocity_y_list = [] - - # Go through each time step and calculate forces and atmospheric values - # Get flight phases - for phase_index, phase in self.time_iterator(self.FlightPhases): - init_time = phase.t - final_time = self.FlightPhases[phase_index + 1].t - current_derivative = phase.derivative - # Call callback functions - for callback in phase.callbacks: - callback(self) - # Loop through time steps in flight phase - for step in self.solution: # Can be optimized - if init_time < step[0] <= final_time or ( - init_time == self.t_initial and step[0] == self.t_initial - ): - # Call derivatives in post processing mode - u_dot = current_derivative(step[0], step[1:], post_processing=True) - - temporary_values = [ - self.R1_list, - self.R2_list, - self.R3_list, - self.M1_list, - self.M2_list, - self.M3_list, - self.pressure_list, - self.density_list, - self.dynamic_viscosity_list, - self.speed_of_sound_list, - self.wind_velocity_x_list, - self.wind_velocity_y_list, - ] - - return temporary_values - def get_controller_observed_variables(self): """Retrieve the observed variables related to air brakes from the controllers. If there is only one set of observed variables, it is @@ -2918,7 +2712,7 @@ def get_controller_observed_variables(self): ) @cached_property - def __calculate_rail_button_forces(self): + def __calculate_rail_button_forces(self): # TODO: complex method. """Calculate the forces applied to the rail buttons while rocket is still on the launch rail. It will return 0 if no rail buttons are defined. @@ -3005,7 +2799,7 @@ def __calculate_rail_button_forces(self): rail_button2_shear_force, ) - def _calculate_pressure_signal(self): + def __transform_pressure_signals_lists_to_functions(self): """Calculate the pressure signal from the pressure sensor. It creates a signal_function attribute in the parachute object. Parachute works as a subclass of Rocket class. @@ -3022,17 +2816,40 @@ def _calculate_pressure_signal(self): "Pressure - Without Noise (Pa)", "linear", ) - parachute.noisy_pressure_signal_function = Function( - parachute.noisy_pressure_signal, - "Time (s)", - "Pressure - With Noise (Pa)", - "linear", - ) parachute.noise_signal_function = Function( parachute.noise_signal, "Time (s)", "Pressure Noise (Pa)", "linear" ) + parachute.noisy_pressure_signal_function = ( + parachute.clean_pressure_signal_function + + parachute.noise_signal_function + ) - return None + @cached_property + def __evaluate_post_process(self): + """Evaluate all post-processing variables by running the simulation + again but with the post-processing flag set to True. + + Returns + ------- + np.array + An array containing all post-processed variables evaluated at each + time step. Each element of the array is a list containing: + [t, ax, ay, az, alpha1, alpha2, alpha3, R1, R2, R3, M1, M2, M3] + """ + self.__post_processed_variables = [] + for phase_index, phase in self.time_iterator(self.flight_phases): + init_time = phase.t + final_time = self.flight_phases[phase_index + 1].t + current_derivative = phase.derivative + for callback in phase.callbacks: + callback(self) + for step in self.solution: + if init_time < step[0] <= final_time or ( + init_time == self.t_initial and step[0] == self.t_initial + ): + current_derivative(step[0], step[1:], post_processing=True) + + return np.array(self.__post_processed_variables) def post_process(self, interpolation="spline", extrapolation="natural"): """This method is **deprecated** and is only kept here for backwards @@ -3048,12 +2865,10 @@ def post_process(self, interpolation="spline", extrapolation="natural"): ------- None """ - # Register post processing + # TODO: add a deprecation warning maybe? self.post_processed = True - return None - - def calculate_stall_wind_velocity(self, stall_angle): + def calculate_stall_wind_velocity(self, stall_angle): # TODO: move to utilities """Function to calculate the maximum wind velocity before the angle of attack exceeds a desired angle, at the instant of departing rail launch. Can be helpful if you know the exact stall angle of all aerodynamics @@ -3064,15 +2879,15 @@ def calculate_stall_wind_velocity(self, stall_angle): stall_angle : float Angle, in degrees, for which you would like to know the maximum wind speed before the angle of attack exceeds it + Return ------ None """ v_f = self.out_of_rail_velocity - # Convert angle to radians - theta = self.inclination * 3.14159265359 / 180 - stall_angle = stall_angle * 3.14159265359 / 180 + theta = np.radians(self.inclination) + stall_angle = np.radians(stall_angle) c = (math.cos(stall_angle) ** 2 - math.cos(theta) ** 2) / math.sin( stall_angle @@ -3086,16 +2901,13 @@ def calculate_stall_wind_velocity(self, stall_angle): ** 0.5 ) / 2 - # Convert stall_angle to degrees - stall_angle = stall_angle * 180 / np.pi + stall_angle = np.degrees(stall_angle) print( "Maximum wind velocity at Rail Departure time before angle" + f" of attack exceeds {stall_angle:.3f}°: {w_v:.3f} m/s" ) - return None - - def export_pressures(self, file_name, time_step): + def export_pressures(self, file_name, time_step): # TODO: move out """Exports the pressure experienced by the rocket during the flight to an external file, the '.csv' format is recommended, as the columns will be separated by commas. It can handle flights with or without @@ -3127,12 +2939,12 @@ def export_pressures(self, file_name, time_step): if len(self.rocket.parachutes) == 0: print("No parachutes in the rocket, saving static pressure.") for t in time_points: - file.write(f"{t:f}, {self.pressure(t):.5f}\n") + file.write(f"{t:f}, {self.pressure.get_value_opt(t):.5f}\n") else: for parachute in self.rocket.parachutes: for t in time_points: - p_cl = parachute.clean_pressure_signal_function(t) - p_ns = parachute.noisy_pressure_signal_function(t) + p_cl = parachute.clean_pressure_signal_function.get_value_opt(t) + p_ns = parachute.noisy_pressure_signal_function.get_value_opt(t) file.write(f"{t:f}, {p_cl:.5f}, {p_ns:.5f}\n") # We need to save only 1 parachute data break @@ -3160,6 +2972,7 @@ class attributes which are instances of the Function class. Usage will be exported. Otherwise, linear interpolation is carried out to calculate values at the desired time steps. Example: 0.001. """ + # TODO: we should move this method to outside of class. # Fast evaluation for the most basic scenario if time_step is None and len(variables) == 0: @@ -3218,13 +3031,13 @@ class attributes which are instances of the Function class. Usage try: obj = getattr(self.__class__, variable) variable_function = obj.__get__(self, self.__class__) - except AttributeError: + except AttributeError as exc: raise AttributeError( - "Variable '{}' not found in Flight class".format(variable) - ) + f"Variable '{variable}' not found in Flight class" + ) from exc variable_points = variable_function(time_points) exported_matrix += [variable_points] - exported_header += ", " + variable_function.__outputs__[0] + exported_header += f", {variable_function.__outputs__[0]}" exported_matrix = np.array(exported_matrix).T # Fix matrix orientation @@ -3237,9 +3050,7 @@ class attributes which are instances of the Function class. Usage encoding="utf-8", ) - return - - def export_kml( + def export_kml( # TODO: should be moved out of this class. self, file_name="trajectory.kml", time_step=None, @@ -3287,27 +3098,33 @@ def export_kml( # Open kml file with simplekml library kml = simplekml.Kml(open=1) trajectory = kml.newlinestring(name="Rocket Trajectory - Powered by RocketPy") - coords = [] + if altitude_mode == "relativetoground": # In this mode the elevation data will be the Above Ground Level # elevation. Only works properly if the ground level is similar to # a plane, i.e. it might not work well if the terrain has mountains - for t in time_points: - coords.append( - ( - self.longitude(t), - self.latitude(t), - self.altitude(t), - ) + coords = [ + ( + self.longitude.get_value_opt(t), + self.latitude.get_value_opt(t), + self.altitude.get_value_opt(t), ) + for t in time_points + ] trajectory.coords = coords trajectory.altitudemode = simplekml.AltitudeMode.relativetoground else: # altitude_mode == 'absolute' # In this case the elevation data will be the Above Sea Level elevation # Ensure you use the correct value on self.env.elevation, otherwise # the trajectory path can be offset from ground - for t in time_points: - coords.append((self.longitude(t), self.latitude(t), self.z(t))) + coords = [ + ( + self.longitude.get_value_opt(t), + self.latitude.get_value_opt(t), + self.z.get_value_opt(t), + ) + for t in time_points + ] trajectory.coords = coords trajectory.altitudemode = simplekml.AltitudeMode.absolute # Modify style of trajectory linestring @@ -3319,29 +3136,13 @@ def export_kml( kml.save(file_name) print("File ", file_name, " saved with success!") - return None - def info(self): - """Prints out a summary of the data available about the Flight. - - Returns - ------- - None - """ + """Prints out a summary of the data available about the Flight.""" self.prints.all() - return None def all_info(self): - """Prints out all data and graphs available about the Flight. - - Returns - ------- - None - """ - - # Print a summary of data about the flight + """Prints out all data and graphs available about the Flight.""" self.info() - self.plots.all() return None @@ -3381,7 +3182,7 @@ def display_warning(self, *messages): """A simple function to print a warning message.""" print("WARNING:", *messages) - def add(self, flight_phase, index=None): + def add(self, flight_phase, index=None): # TODO: quite complex method """Add a flight phase to the list. It will be inserted in the correct position, according to its initial time. If no index is provided, it will be appended to the end of the list. If by any @@ -3538,6 +3339,8 @@ class FlightPhase: A flag indicating whether to clear the solution after the phase. """ + # TODO: add a "name" optional argument to the FlightPhase. Really helps. + def __init__(self, t, derivative=None, callbacks=None, clear=True): self.t = t self.derivative = derivative @@ -3545,17 +3348,19 @@ def __init__(self, t, derivative=None, callbacks=None, clear=True): self.clear = clear def __repr__(self): - if self.derivative is None: - return "{Initial Time: " + str(self.t) + " | Derivative: None}" + name = "None" if self.derivative is None else self.derivative.__name__ return ( - "{Initial Time: " - + str(self.t) - + " | Derivative: " - + self.derivative.__name__ - + "}" + f"" ) class TimeNodes: + """TimeNodes is a class that stores all the time nodes of a simulation. + It is meant to work like a python list, but it has some additional + methods that are useful for the simulation. Them items stored in are + TimeNodes object are instances of the TimeNode class. + """ + def __init__(self, init_list=[]): self.list = init_list[:] @@ -3575,20 +3380,19 @@ def add_node(self, t, parachutes, callbacks): self.list.append(self.TimeNode(t, parachutes, callbacks)) def add_parachutes(self, parachutes, t_init, t_end): - # Iterate over parachutes for parachute in parachutes: # Calculate start of sampling time nodes - pcDt = 1 / parachute.sampling_rate + sampling_interval = 1 / parachute.sampling_rate parachute_node_list = [ - self.TimeNode(i * pcDt, [parachute], []) + self.TimeNode(i * sampling_interval, [parachute], []) for i in range( - math.ceil(t_init / pcDt), math.floor(t_end / pcDt) + 1 + math.ceil(t_init / sampling_interval), + math.floor(t_end / sampling_interval) + 1, ) ] self.list += parachute_node_list def add_controllers(self, controllers, t_init, t_end): - # Iterate over controllers for controller in controllers: # Calculate start of sampling time nodes controller_time_step = 1 / controller.sampling_rate @@ -3602,29 +3406,50 @@ def add_controllers(self, controllers, t_init, t_end): self.list += controller_node_list def sort(self): - self.list.sort(key=(lambda node: node.t)) + self.list.sort() def merge(self): - # Initialize temporary list - self.tmp_list = [self.list[0]] - self.copy_list = self.list[1:] - # Iterate through all other time nodes - for node in self.copy_list: - # If there is already another node with similar time: merge - if abs(node.t - self.tmp_list[-1].t) < 1e-7: - self.tmp_list[-1].parachutes += node.parachutes - self.tmp_list[-1].callbacks += node.callbacks - # Add new node to tmp list if there is none with the same time - else: - self.tmp_list.append(node) - # Save tmp list to permanent - self.list = self.tmp_list + """Merge all the time nodes that have the same time. This is made to + avoid multiple evaluations of the same time node. This method does + not guarantee the order of the nodes in the list, so it is + recommended to sort the list before or after using this method. + """ + tmp_dict = {} + for node in self.list: + time = round(node.t, 7) + try: + # Try to access the node and merge if it exists + tmp_dict[time].parachutes += node.parachutes + tmp_dict[time].callbacks += node.callbacks + except KeyError: + # If the node does not exist, add it to the dictionary + tmp_dict[time] = node + self.list = list(tmp_dict.values()) def flush_after(self, index): del self.list[index + 1 :] class TimeNode: + """TimeNode is a class that represents a time node in the time + nodes list. It stores the time, the parachutes and the controllers + that are active at that time. This class is supposed to work + exclusively within the TimeNodes class. + """ + def __init__(self, t, parachutes, controllers): + """Create a TimeNode object. + + Parameters + ---------- + t : float + Initial time of the time node. + parachutes : list[Parachute] + List containing all the parachutes that should be evaluated + at this time node. + controllers : list[_Controller] + List containing all the controllers that should be evaluated + at this time node. + """ self.t = t self.parachutes = parachutes self.callbacks = [] @@ -3632,9 +3457,27 @@ def __init__(self, t, parachutes, controllers): def __repr__(self): return ( - "{Initial Time: " - + str(self.t) - + " | Parachutes: " - + str(len(self.parachutes)) - + "}" + f"" ) + + def __lt__(self, other): + """Allows the comparison of two TimeNode objects based on their + initial time. This is particularly useful for sorting a list of + TimeNode objects. + + Parameters + ---------- + other : TimeNode + Another TimeNode object to compare with. + + Returns + ------- + bool + True if the initial time of the current TimeNode is less + than the initial time of the other TimeNode, False + otherwise. + """ + return self.t < other.t diff --git a/rocketpy/tools.py b/rocketpy/tools.py index bcff91fb8..86ad7f17e 100644 --- a/rocketpy/tools.py +++ b/rocketpy/tools.py @@ -1,6 +1,16 @@ +"""The module rocketpy.tools contains a set of functions that are used +throughout the rocketpy package. These functions are not specific to any +particular class or module, and are used to perform general tasks that are +required by multiple classes or modules. These functions can be modified or +expanded to suit the needs of other modules and may present breaking changes +between minor versions if necessary, although this will be always avoided. +""" + +import functools import importlib import importlib.metadata import re +import time from bisect import bisect_left import numpy as np @@ -8,36 +18,10 @@ from cftime import num2pydate from packaging import version as packaging_version -_NOT_FOUND = object() - # Mapping of module name and the name of the package that should be installed INSTALL_MAPPING = {"IPython": "ipython"} -class cached_property: - def __init__(self, func): - self.func = func - self.attrname = None - self.__doc__ = func.__doc__ - - def __set_name__(self, owner, name): - self.attrname = name - - def __get__(self, instance, owner=None): - if instance is None: - return self - if self.attrname is None: - raise TypeError( - "Cannot use cached_property instance without calling __set_name__ on it." - ) - cache = instance.__dict__ - val = cache.get(self.attrname, _NOT_FOUND) - if val is _NOT_FOUND: - val = self.func(instance) - cache[self.attrname] = val - return val - - def tuple_handler(value): """Transforms the input value into a tuple that represents a range. If the input is an input or float, @@ -66,6 +50,139 @@ def tuple_handler(value): raise ValueError("value must be a list or tuple of length 1 or 2.") +def calculate_cubic_hermite_coefficients(x0, x1, y0, yp0, y1, yp1): + """Calculate the coefficients of a cubic Hermite interpolation function. + The function is defined as ax**3 + bx**2 + cx + d. + + Parameters + ---------- + x0 : float + Position of the first point. + x1 : float + Position of the second point. + y0 : float + Value of the function evaluated at the first point. + yp0 : float + Value of the derivative of the function evaluated at the first + point. + y1 : float + Value of the function evaluated at the second point. + yp1 : float + Value of the derivative of the function evaluated at the second + point. + + Returns + ------- + tuple[float, float, float, float] + The coefficients of the cubic Hermite interpolation function. + """ + dx = x1 - x0 + d = float(y0) + c = float(yp0) + b = float((3 * y1 - yp1 * dx - 2 * c * dx - 3 * d) / (dx**2)) + a = float(-(2 * y1 - yp1 * dx - c * dx - 2 * d) / (dx**3)) + return a, b, c, d + + +def find_roots_cubic_function(a, b, c, d): + """Calculate the roots of a cubic function using Cardano's method. + + This method applies Cardano's method to find the roots of a cubic + function of the form ax^3 + bx^2 + cx + d. The roots may be complex + numbers. + + Parameters + ---------- + a : float + Coefficient of the cubic term (x^3). + b : float + Coefficient of the quadratic term (x^2). + c : float + Coefficient of the linear term (x). + d : float + Constant term. + + Returns + ------- + tuple[complex, complex, complex] + A tuple containing the real and complex roots of the cubic function. + Note that the roots may be complex numbers. The roots are ordered + in the tuple as x1, x2, x3. + + References + ---------- + - Cardano's method: https://en.wikipedia.org/wiki/Cubic_function#Cardano's_method + + Examples + -------- + >>> from rocketpy.tools import find_roots_cubic_function + + First we define the coefficients of the function ax**3 + bx**2 + cx + d + >>> a, b, c, d = 1, -3, -1, 3 + >>> x1, x2, x3 = find_roots_cubic_function(a, b, c, d) + >>> x1, x2, x3 + ((-1+0j), (3+7.401486830834377e-17j), (1-1.4802973661668753e-16j)) + + To get the real part of the roots, use the real attribute of the complex + number. + >>> x1.real, x2.real, x3.real + (-1.0, 3.0, 1.0) + """ + delta_0 = b**2 - 3 * a * c + delta_1 = 2 * b**3 - 9 * a * b * c + 27 * d * a**2 + c1 = ((delta_1 + (delta_1**2 - 4 * delta_0**3) ** (0.5)) / 2) ** (1 / 3) + + c2_0 = c1 + x1 = -(1 / (3 * a)) * (b + c2_0 + delta_0 / c2_0) + + c2_1 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** 1 + x2 = -(1 / (3 * a)) * (b + c2_1 + delta_0 / c2_1) + + c2_2 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** 2 + x3 = -(1 / (3 * a)) * (b + c2_2 + delta_0 / c2_2) + + return x1, x2, x3 + + +def find_root_linear_interpolation(x0, x1, y0, y1, y): + """Calculate the root of a linear interpolation function. + + This method calculates the root of a linear interpolation function + given two points (x0, y0) and (x1, y1) and a value y. The function + is defined as y = m*x + c. + + Parameters + ---------- + x0 : float + Position of the first point. + x1 : float + Position of the second point. + y0 : float + Value of the function evaluated at the first point. + y1 : float + Value of the function evaluated at the second point. + y : float + Value of the function at the desired point. + + Returns + ------- + float + The root of the linear interpolation function. This represents the + value of x at which the function evaluates to y. + + Examples + -------- + >>> from rocketpy.tools import find_root_linear_interpolation + >>> x0, x1, y0, y1, y = 0, 1, 0, 1, 0.5 + >>> x = find_root_linear_interpolation(x0, x1, y0, y1, y) + >>> x + 0.5 + """ + m = (y1 - y0) / (x1 - x0) + c = y0 - m * x0 + return (y - c) / m + + def bilinear_interpolation(x, y, x1, x2, y1, y2, z11, z12, z21, z22): """Bilinear interpolation. It considers the values of the four points around the point to be interpolated and returns the interpolated value. @@ -153,7 +270,7 @@ def time_num_to_date_string(time_num, units, timezone, calendar="gregorian"): """Convert time number (usually hours before a certain date) into two strings: one for the date (example: 2022.04.31) and one for the hour (example: 14). See cftime.num2date for details on units and calendar. - Automatically converts time number from UTC to local timezone based on + Automatically converts time number from UTC to local time zone based on lat, lon coordinates. This function was created originally for the EnvironmentAnalysis class. @@ -382,6 +499,52 @@ def check_requirement_version(module_name, version): return True +def exponential_backoff(max_attempts, base_delay=1, max_delay=60): + def decorator(func): + @functools.wraps(func) + def wrapper(*args, **kwargs): + delay = base_delay + for i in range(max_attempts): + try: + return func(*args, **kwargs) + except Exception as e: + if i == max_attempts - 1: + raise e from None + delay = min(delay * 2, max_delay) + time.sleep(delay) + + return wrapper + + return decorator + + +def parallel_axis_theorem_from_com(com_inertia_moment, mass, distance): + """Calculates the moment of inertia of a object relative to a new axis using + the parallel axis theorem. The new axis is parallel to and at a distance + 'distance' from the original axis, which *must* passes through the object's + center of mass. + + Parameters + ---------- + com_inertia_moment : float + Moment of inertia relative to the center of mass of the object. + mass : float + Mass of the object. + distance : float + Perpendicular distance between the original and new axis. + + Returns + ------- + float + Moment of inertia relative to the new axis. + + Reference + --------- + https://en.wikipedia.org/wiki/Parallel_axis_theorem + """ + return com_inertia_moment + mass * distance**2 + + # Flight diff --git a/setup.py b/setup.py deleted file mode 100644 index 73d989de4..000000000 --- a/setup.py +++ /dev/null @@ -1,46 +0,0 @@ -import setuptools - -with open("README.md", "r") as fh: - long_description = fh.read() - -necessary_require = [ - "numpy>=1.13", - "scipy>=1.0", - "matplotlib>=3.0", - "netCDF4>=1.6.4", - "requests", - "pytz", - "simplekml", -] - -env_analysis_require = [ - "timezonefinder", - "windrose>=1.6.8", - "IPython", - "ipywidgets>=7.6.3", - "jsonpickle", -] - -setuptools.setup( - name="rocketpy", - version="1.2.1", - install_requires=necessary_require, - extras_require={ - "env_analysis": env_analysis_require, - "all": necessary_require + env_analysis_require, - }, - maintainer="RocketPy Developers", - author="Giovani Hidalgo Ceotto", - author_email="ghceotto@gmail.com", - description="Advanced 6-DOF trajectory simulation for High-Power Rocketry.", - long_description=long_description, - long_description_content_type="text/markdown", - url="https://github.com/RocketPy-Team/RocketPy", - packages=setuptools.find_packages(), - classifiers=[ - "Programming Language :: Python :: 3", - "License :: OSI Approved :: MIT License", - "Operating System :: OS Independent", - ], - python_requires=">=3.8", -) diff --git a/tests/acceptance/test_bella_lui_rocket.py b/tests/acceptance/test_bella_lui_rocket.py index 69669c66f..0041074a8 100644 --- a/tests/acceptance/test_bella_lui_rocket.py +++ b/tests/acceptance/test_bella_lui_rocket.py @@ -228,10 +228,14 @@ def drogue_trigger(p, h, y): apogee_time_measured = time_kalt[np.argmax(altitude_kalt)] apogee_time_simulated = test_flight.apogee_time - assert ( - abs(max(altitude_kalt) - test_flight.apogee + test_flight.env.elevation) - / max(altitude_kalt) - < 0.015 + apogee_error_threshold = 0.015 + apogee_error = abs( + max(altitude_kalt) - test_flight.apogee + test_flight.env.elevation + ) / max(altitude_kalt) + assert apogee_error < apogee_error_threshold, ( + f"Apogee altitude error exceeded the threshold. " + f"Expected the error to be less than {apogee_error_threshold * 100}%, " + f"but got an error of {apogee_error * 100:.1f}%." ) assert abs(max(velocity_rcp) - max(vert_vel_kalt)) / max(vert_vel_kalt) < 0.06 assert ( @@ -242,3 +246,10 @@ def drogue_trigger(p, h, y): assert ( abs(apogee_time_measured - apogee_time_simulated) / apogee_time_simulated < 0.02 ) + # Guarantee the impact velocity is within 30% of the real data. + # Use the last 5 real points to avoid outliers + assert ( + abs(test_flight.impact_velocity - np.mean(vert_vel_kalt[-5:])) + / abs(test_flight.impact_velocity) + < 0.30 + ) diff --git a/tests/fixtures/environment/environment_fixtures.py b/tests/fixtures/environment/environment_fixtures.py index 8949f9973..851be3203 100644 --- a/tests/fixtures/environment/environment_fixtures.py +++ b/tests/fixtures/environment/environment_fixtures.py @@ -1,6 +1,7 @@ from datetime import datetime, timedelta import pytest + from rocketpy import Environment, EnvironmentAnalysis @@ -54,8 +55,8 @@ def env_analysis(): EnvironmentAnalysis """ env_analysis = EnvironmentAnalysis( - start_date=datetime.datetime(2019, 10, 23), - end_date=datetime.datetime(2021, 10, 23), + start_date=datetime(2019, 10, 23), + end_date=datetime(2021, 10, 23), latitude=39.3897, longitude=-8.28896388889, start_hour=6, diff --git a/tests/test_environment.py b/tests/test_environment.py index 5fa0e2c45..7349d512b 100644 --- a/tests/test_environment.py +++ b/tests/test_environment.py @@ -1,5 +1,4 @@ import datetime -import time from unittest.mock import patch import pytest @@ -64,13 +63,8 @@ def test_wyoming_sounding_atmosphere(mock_show, example_plain_env): # "file" option, instead of receiving the URL as a string. URL = "http://weather.uwyo.edu/cgi-bin/sounding?region=samer&TYPE=TEXT%3ALIST&YEAR=2019&MONTH=02&FROM=0500&TO=0512&STNM=83779" # give it at least 5 times to try to download the file - for i in range(5): - try: - example_plain_env.set_atmospheric_model(type="wyoming_sounding", file=URL) - break - except: - time.sleep(1) # wait 1 second before trying again - pass + example_plain_env.set_atmospheric_model(type="wyoming_sounding", file=URL) + assert example_plain_env.all_info() == None assert abs(example_plain_env.pressure(0) - 93600.0) < 1e-8 assert ( diff --git a/tests/test_flight.py b/tests/test_flight.py index 506bdc561..d4ed1e114 100644 --- a/tests/test_flight.py +++ b/tests/test_flight.py @@ -604,12 +604,12 @@ def test_max_values(flight_calisto_robust): calculated by hand, it was just copied from the test results. This is because the expected values are not easy to calculate by hand, and the results are not expected to change. If the results change, the test will - fail, and the expected values must be updated. If if want to update the - values, always double check if the results are really correct. Acceptable - reasons for changes in the results are: 1) changes in the code that - improve the accuracy of the simulation, 2) a bug was found and fixed. Keep - in mind that other tests may be more accurate than this one, for example, - the acceptance tests, which are based on the results of real flights. + fail, and the expected values must be updated. If the values are updated, + always double check if the results are really correct. Acceptable reasons + for changes in the results are: 1) changes in the code that improve the + accuracy of the simulation, 2) a bug was found and fixed. Keep in mind that + other tests may be more accurate than this one, for example, the acceptance + tests, which are based on the results of real flights. Parameters ---------- @@ -618,11 +618,11 @@ def test_max_values(flight_calisto_robust): regarding this pytest fixture. """ test = flight_calisto_robust - atol = 1e-2 - assert pytest.approx(105.2774, abs=atol) == test.max_acceleration_power_on - assert pytest.approx(105.2774, abs=atol) == test.max_acceleration - assert pytest.approx(0.85999, abs=atol) == test.max_mach_number - assert pytest.approx(285.90240, abs=atol) == test.max_speed + rtol = 5e-3 + assert pytest.approx(105.1599, rel=rtol) == test.max_acceleration_power_on + assert pytest.approx(105.1599, rel=rtol) == test.max_acceleration + assert pytest.approx(0.85999, rel=rtol) == test.max_mach_number + assert pytest.approx(285.94948, rel=rtol) == test.max_speed def test_rail_buttons_forces(flight_calisto_custom_wind): @@ -691,7 +691,7 @@ def test_accelerations(flight_calisto_custom_wind, flight_time, expected_values) ("t_initial", (0, 0, 0)), ("out_of_rail_time", (0, 2.248727, 25.703072)), ("apogee_time", (-13.209436, 16.05115, -0.000257)), - ("t_final", (5, 2, -5.334289)), + ("t_final", (5, 2, -5.65998)), ], ) def test_velocities(flight_calisto_custom_wind, flight_time, expected_values): @@ -729,7 +729,7 @@ def test_velocities(flight_calisto_custom_wind, flight_time, expected_values): ("t_initial", (1.6542528, 0.65918, -0.067107)), ("out_of_rail_time", (5.05334, 2.01364, -1.7541)), ("apogee_time", (2.35291, -1.8275, -0.87851)), - ("t_final", (0, 0, 141.42421)), + ("t_final", (0, 0, 159.2212)), ], ) def test_aerodynamic_forces(flight_calisto_custom_wind, flight_time, expected_values): diff --git a/tests/test_function.py b/tests/test_function.py index 2ce94f691..6f4122e47 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -102,7 +102,8 @@ def test_setters(func_from_csv, func_2d_from_csv): func_2d_from_csv.set_interpolation("shepard") assert func_2d_from_csv.get_interpolation_method() == "shepard" func_2d_from_csv.set_extrapolation("zero") - assert func_2d_from_csv.get_extrapolation_method() == "zero" + # 2d functions do not support zero extrapolation, must change to natural + assert func_2d_from_csv.get_extrapolation_method() == "natural" @patch("matplotlib.pyplot.show") @@ -181,7 +182,32 @@ def test_extrapolation_methods(linear_func): assert linear_func.get_extrapolation_method() == "constant" assert np.isclose(linear_func.get_value(-1), 0, atol=1e-6) - # Test natural + # Test natural for linear interpolation + linear_func.set_interpolation("linear") + assert isinstance(linear_func.set_extrapolation("natural"), Function) + linear_func.set_extrapolation("natural") + assert isinstance(linear_func.get_extrapolation_method(), str) + assert linear_func.get_extrapolation_method() == "natural" + assert np.isclose(linear_func.get_value(-1), -1, atol=1e-6) + + # Test natural for spline interpolation + linear_func.set_interpolation("spline") + assert isinstance(linear_func.set_extrapolation("natural"), Function) + linear_func.set_extrapolation("natural") + assert isinstance(linear_func.get_extrapolation_method(), str) + assert linear_func.get_extrapolation_method() == "natural" + assert np.isclose(linear_func.get_value(-1), -1, atol=1e-6) + + # Test natural for akima interpolation + linear_func.set_interpolation("akima") + assert isinstance(linear_func.set_extrapolation("natural"), Function) + linear_func.set_extrapolation("natural") + assert isinstance(linear_func.get_extrapolation_method(), str) + assert linear_func.get_extrapolation_method() == "natural" + assert np.isclose(linear_func.get_value(-1), -1, atol=1e-6) + + # Test natural for polynomial interpolation + linear_func.set_interpolation("polynomial") assert isinstance(linear_func.set_extrapolation("natural"), Function) linear_func.set_extrapolation("natural") assert isinstance(linear_func.get_extrapolation_method(), str) diff --git a/tests/unit/test_environment.py b/tests/unit/test_environment.py index ac25533eb..8d676f426 100644 --- a/tests/unit/test_environment.py +++ b/tests/unit/test_environment.py @@ -1,10 +1,10 @@ +import json import os import numpy as np import numpy.ma as ma import pytest import pytz -import json from rocketpy import Environment diff --git a/tests/unit/test_flight.py b/tests/unit/test_flight.py index 87249cf44..e6ab6b8b8 100644 --- a/tests/unit/test_flight.py +++ b/tests/unit/test_flight.py @@ -258,3 +258,31 @@ def test_get_controller_observed_variables(flight_calisto_air_brakes): obs_vars = flight_calisto_air_brakes.get_controller_observed_variables() assert isinstance(obs_vars, list) assert len(obs_vars) == 0 + + +def test_initial_stability_margin(flight_calisto_custom_wind): + """Test the initial_stability_margin method of the Flight class. + + Parameters + ---------- + flight_calisto_custom_wind : rocketpy.Flight + """ + res = flight_calisto_custom_wind.initial_stability_margin + assert isinstance(res, float) + assert res == flight_calisto_custom_wind.stability_margin(0) + assert np.isclose(res, 2.05, atol=0.1) + + +def test_out_of_rail_stability_margin(flight_calisto_custom_wind): + """Test the out_of_rail_stability_margin method of the Flight class. + + Parameters + ---------- + flight_calisto_custom_wind : rocketpy.Flight + """ + res = flight_calisto_custom_wind.out_of_rail_stability_margin + assert isinstance(res, float) + assert res == flight_calisto_custom_wind.stability_margin( + flight_calisto_custom_wind.out_of_rail_time + ) + assert np.isclose(res, 2.14, atol=0.1) diff --git a/tests/unit/test_flight_time_nodes.py b/tests/unit/test_flight_time_nodes.py new file mode 100644 index 000000000..10f6b6c30 --- /dev/null +++ b/tests/unit/test_flight_time_nodes.py @@ -0,0 +1,103 @@ +"""Module to test everything related to the TimeNodes class and it's subclass +TimeNode. +""" + +import pytest + +from rocketpy.rocket import Parachute, _Controller + + +def test_time_nodes_init(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + assert len(time_nodes) == 0 + + +def test_time_nodes_getitem(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + time_nodes.add_node(1.0, [], []) + assert isinstance(time_nodes[0], flight_calisto.TimeNodes.TimeNode) + assert time_nodes[0].t == 1.0 + + +def test_time_nodes_len(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + assert len(time_nodes) == 0 + + +def test_time_nodes_add(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + example_node = flight_calisto.TimeNodes.TimeNode(1.0, [], []) + time_nodes.add(example_node) + assert len(time_nodes) == 1 + assert isinstance(time_nodes[0], flight_calisto.TimeNodes.TimeNode) + assert time_nodes[0].t == 1.0 + + +def test_time_nodes_add_node(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + time_nodes.add_node(2.0, [], []) + assert len(time_nodes) == 1 + assert time_nodes[0].t == 2.0 + assert len(time_nodes[0].parachutes) == 0 + assert len(time_nodes[0].callbacks) == 0 + + +# def test_time_nodes_add_parachutes( +# flight_calisto, calisto_drogue_chute, calisto_main_chute +# ): # TODO: implement this test + + +# def test_time_nodes_add_controllers(flight_calisto): +# TODO: implement this test + + +def test_time_nodes_sort(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + time_nodes.add_node(3.0, [], []) + time_nodes.add_node(1.0, [], []) + time_nodes.add_node(2.0, [], []) + time_nodes.sort() + assert len(time_nodes) == 3 + assert time_nodes[0].t == 1.0 + assert time_nodes[1].t == 2.0 + assert time_nodes[2].t == 3.0 + + +def test_time_nodes_merge(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + time_nodes.add_node(1.0, [], []) + time_nodes.add_node(1.0, [], []) + time_nodes.add_node(2.0, [], []) + time_nodes.merge() + assert len(time_nodes) == 2 + assert time_nodes[0].t == 1.0 + assert len(time_nodes[0].parachutes) == 0 + assert len(time_nodes[0].callbacks) == 0 + assert time_nodes[1].t == 2.0 + assert len(time_nodes[1].parachutes) == 0 + assert len(time_nodes[1].callbacks) == 0 + + +def test_time_nodes_flush_after(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + time_nodes.add_node(1.0, [], []) + time_nodes.add_node(2.0, [], []) + time_nodes.add_node(3.0, [], []) + time_nodes.flush_after(1) + assert len(time_nodes) == 2 + assert time_nodes[0].t == 1.0 + assert time_nodes[1].t == 2.0 + + +def test_time_node_init(flight_calisto): + node = flight_calisto.TimeNodes.TimeNode(1.0, [], []) + assert node.t == 1.0 + assert len(node.parachutes) == 0 + assert len(node.callbacks) == 0 + + +def test_time_node_lt(flight_calisto): + node1 = flight_calisto.TimeNodes.TimeNode(1.0, [], []) + node2 = flight_calisto.TimeNodes.TimeNode(2.0, [], []) + assert node1 < node2 + assert not node2 < node1 diff --git a/tests/unit/test_function.py b/tests/unit/test_function.py index 8bcefb818..17da59498 100644 --- a/tests/unit/test_function.py +++ b/tests/unit/test_function.py @@ -306,3 +306,31 @@ def test_remove_outliers_iqr(x, y, expected_x, expected_y): assert filtered_func.__interpolation__ == func.__interpolation__ assert filtered_func.__extrapolation__ == func.__extrapolation__ assert filtered_func.title == func.title + + +def test_set_get_value_opt(): + """Test the set_value_opt and get_value_opt methods of the Function class.""" + func = Function(lambda x: x**2) + func.source = np.array([[1, 1], [2, 4], [3, 9], [4, 16], [5, 25]]) + func.x_array = np.array([1, 2, 3, 4, 5]) + func.y_array = np.array([1, 4, 9, 16, 25]) + func.x_initial = 1 + func.x_final = 5 + func.set_interpolation("linear") + func.set_get_value_opt() + assert func.get_value_opt(2.5) == 6.5 + + +def test_get_image_dim(linear_func): + """Test the get_img_dim method of the Function class.""" + assert linear_func.get_image_dim() == 1 + + +def test_get_domain_dim(linear_func): + """Test the get_domain_dim method of the Function class.""" + assert linear_func.get_domain_dim() == 1 + + +def test_bool(linear_func): + """Test the __bool__ method of the Function class.""" + assert bool(linear_func) == True diff --git a/tests/unit/test_tools.py b/tests/unit/test_tools.py new file mode 100644 index 000000000..75a526aac --- /dev/null +++ b/tests/unit/test_tools.py @@ -0,0 +1,46 @@ +import numpy as np + +from rocketpy.tools import ( + calculate_cubic_hermite_coefficients, + find_roots_cubic_function, +) + + +def test_calculate_cubic_hermite_coefficients(): + """Test the calculate_cubic_hermite_coefficients method of the Function class.""" + # Function: f(x) = x**3 + 2x**2 -1 ; derivative: f'(x) = 3x**2 + 4x + x = np.array([-3, -2, -1, 0, 1]) + y = np.array([-10, -1, 0, -1, 2]) + + # Selects two points as x0 and x1 + x0, x1 = 0, 1 + y0, y1 = -1, 2 + yp0, yp1 = 0, 7 + + a, b, c, d = calculate_cubic_hermite_coefficients(x0, x1, y0, yp0, y1, yp1) + + assert np.isclose(a, 1) + assert np.isclose(b, 2) + assert np.isclose(c, 0) + assert np.isclose(d, -1) + assert np.allclose( + a * x**3 + b * x**2 + c * x + d, + y, + ) + + +def test_cardanos_root_finding(): + """Tests the find_roots_cubic_function method of the Function class.""" + # Function: f(x) = x**3 + 2x**2 -1 + # roots: (-1 - 5**0.5) / 2; -1; (-1 + 5**0.5) / 2 + + roots = list(find_roots_cubic_function(a=1, b=2, c=0, d=-1)) + roots.sort(key=lambda x: x.real) + + assert np.isclose(roots[0].real, (-1 - 5**0.5) / 2) + assert np.isclose(roots[1].real, -1) + assert np.isclose(roots[2].real, (-1 + 5**0.5) / 2) + + assert np.isclose(roots[0].imag, 0) + assert np.isclose(roots[1].imag, 0) + assert np.isclose(roots[2].imag, 0) diff --git a/tests/unit/test_utilities.py b/tests/unit/test_utilities.py index 68dbc3d7f..33942b445 100644 --- a/tests/unit/test_utilities.py +++ b/tests/unit/test_utilities.py @@ -153,7 +153,7 @@ def test_fin_flutter_analysis(flight_calisto_custom_wind): assert np.isclose(flutter_mach(np.inf), 1.0048188594647927, atol=5e-3) assert np.isclose(safety_factor(0), 64.78797, atol=5e-3) assert np.isclose(safety_factor(10), 2.1948620401502072, atol=5e-3) - assert np.isclose(safety_factor(np.inf), 65.40588722032527, atol=5e-3) + assert np.isclose(safety_factor(np.inf), 61.64222220697017, atol=5e-3) def test_flutter_prints(flight_calisto_custom_wind):