diff --git a/docs/notebooks/deployable_payload_example.ipynb b/docs/notebooks/deployable_payload_example.ipynb
new file mode 100644
index 000000000..50f15525f
--- /dev/null
+++ b/docs/notebooks/deployable_payload_example.ipynb
@@ -0,0 +1,698 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Deployable Payload Flight Simulation Example"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Today we try to demonstrate how to use RocketPy to simulate a flight of a rocket\n",
+ "that presents a deployable payload."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "To run this notebook, we will need:\n",
+ "\n",
+ "* RocketPy\n",
+ "* netCDF4 (to get weather forecasts)\n",
+ "* Data files (we will clone RocketPy's repository for these)\n",
+ "\n",
+ "Therefore, let's run the following lines of code:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "!pip install rocketpy netCDF4\n",
+ "!git clone https://github.com/RocketPy-Team/RocketPy.git"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import os\n",
+ "\n",
+ "os.chdir(\"RocketPy/docs/notebooks\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now we can start!\n",
+ "\n",
+ "Here we go through a simplified rocket trajectory simulation to get you started. Let's start by importing the rocketpy module."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 31,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from rocketpy import Environment, SolidMotor, Rocket, Flight, Function, utilities"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "If you are using a version of Jupyter Notebooks, it is recommended to run the following lines of code to make plots which will be shown later interactive and/or higher quality."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 32,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Using Google Colab? Uncomment the following line:\n",
+ "%config InlineBackend.figure_formats = ['svg']\n",
+ "%matplotlib inline\n",
+ "\n",
+ "# Using Jupyter Notebook/Lab or VSCode? Uncomment the following line:\n",
+ "# %matplotlib widget"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Setting Up a Simulation"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Creating an Environment for Spaceport America"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 33,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "Env = Environment(\n",
+ " railLength=5.2, latitude=32.990254, longitude=-106.974998, elevation=1400\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "To get weather data from the GFS forecast, available online, we run the following lines.\n",
+ "See [Environment Class Usage](environment_class_usage.ipynb) for more information on how to use the Environment class."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 34,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import datetime\n",
+ "\n",
+ "tomorrow = datetime.date.today() + datetime.timedelta(days=1)\n",
+ "\n",
+ "Env.setDate((tomorrow.year, tomorrow.month, tomorrow.day, 12)) # Hour given in UTC time\n",
+ "\n",
+ "Env.setAtmosphericModel(type=\"Forecast\", file=\"GFS\")\n",
+ "Env.maxExpectedHeight = 8000"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 35,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Launch Site Details\n",
+ "\n",
+ "Launch Rail Length: 5.2 m\n",
+ "Launch Date: 2022-10-06 12:00:00 UTC\n",
+ "Launch Site Latitude: 32.99025°\n",
+ "Launch Site Longitude: -106.97500°\n",
+ "Reference Datum: SIRGAS2000\n",
+ "Launch Site UTM coordinates: 315468.64 W 3651938.65 N\n",
+ "Launch Site UTM zone: 13S\n",
+ "Launch Site Surface Elevation: 1471.5 m\n",
+ "\n",
+ "\n",
+ "Atmospheric Model Details\n",
+ "\n",
+ "Atmospheric Model Type: Forecast\n",
+ "Forecast Maximum Height: 8.000 km\n",
+ "Forecast Time Period: From 2022-10-05 00:00:00 to 2022-10-21 00:00:00 UTC\n",
+ "Forecast Hour Interval: 3 hrs\n",
+ "Forecast Latitude Range: From -90.0 ° To 90.0 °\n",
+ "Forecast Longitude Range: From 0.0 ° To 359.75 °\n",
+ "\n",
+ "\n",
+ "Surface Atmospheric Conditions\n",
+ "\n",
+ "Surface Wind Speed: 3.53 m/s\n",
+ "Surface Wind Direction: 9.63°\n",
+ "Surface Wind Heading: 189.63°\n",
+ "Surface Pressure: 858.08 hPa\n",
+ "Surface Temperature: 287.52 K\n",
+ "Surface Air Density: 1.040 kg/m³\n",
+ "Surface Speed of Sound: 339.92 m/s\n",
+ "\n",
+ "\n",
+ "Atmospheric Model Plots\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": "\n\n\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "Env.info()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Creating a Motor\n",
+ "\n",
+ "A solid rocket motor is used in this case. See [Solid Motor Class Usage](solid_motor_class_usage.ipynb) for more information on how to use the Motor class."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 36,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "Pro75M1670 = SolidMotor(\n",
+ " thrustSource=\"../../data/motors/Cesaroni_M1670.eng\",\n",
+ " burnOut=3.9,\n",
+ " grainNumber=5,\n",
+ " grainSeparation=5 / 1000,\n",
+ " grainDensity=1815,\n",
+ " grainOuterRadius=33 / 1000,\n",
+ " grainInitialInnerRadius=15 / 1000,\n",
+ " grainInitialHeight=120 / 1000,\n",
+ " nozzleRadius=33 / 1000,\n",
+ " throatRadius=11 / 1000,\n",
+ " interpolationMethod=\"linear\",\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 37,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ "Motor Details\n",
+ "Total Burning Time: 3.9 s\n",
+ "Total Propellant Mass: 2.956 kg\n",
+ "Propellant Exhaust Velocity: 2038.745 m/s\n",
+ "Average Thrust: 1545.218 N\n",
+ "Maximum Thrust: 2200.0 N at 0.15 s after ignition.\n",
+ "Total Impulse: 6026.350 Ns\n",
+ "\n",
+ "Plots\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": "\n\n\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "Pro75M1670.info()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Simulating the First Flight Stage"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Let's start to simulate our rocket's flight. We will use the Environment and Motor objects we created before.\n",
+ "\n",
+ "We will assume that the payload is be ejected at apogee, however, this can be modified if needed."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We start by defining the value of each relevant mass, ensuring the are correct before continuing."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 48,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Rocket dry mass: 16.24 kg (with Payload)\n",
+ "Propellant Mass: 2.956 kg\n",
+ "Payload Mass: 4.5 kg\n",
+ "Fully loaded Rocket Mass: 19.2 kg\n"
+ ]
+ }
+ ],
+ "source": [
+ "# 16.241 is the mass of the rocket including the payload but without the propellant\n",
+ "PayloadMass = 4.5 # in kg\n",
+ "RocketMass = 16.241 - PayloadMass # in kg\n",
+ "\n",
+ "print(\"Rocket dry mass: {:.4} kg (with Payload)\".format(RocketMass + PayloadMass))\n",
+ "print(\"Propellant Mass: {:.4} kg\".format(Pro75M1670.mass(0)))\n",
+ "print(\"Payload Mass: {:.4} kg\".format(PayloadMass))\n",
+ "print(\n",
+ " \"Fully loaded Rocket Mass: {:.4} kg\".format(\n",
+ " RocketMass + Pro75M1670.mass(0) + PayloadMass\n",
+ " )\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 49,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "Rocket1 = Rocket(\n",
+ " motor=Pro75M1670,\n",
+ " radius=127 / 2000,\n",
+ " mass=RocketMass + PayloadMass,\n",
+ " inertiaI=6.60,\n",
+ " inertiaZ=0.0351,\n",
+ " distanceRocketNozzle=-1.255,\n",
+ " distanceRocketPropellant=-0.85704,\n",
+ " powerOffDrag=\"../../data/calisto/powerOffDragCurve.csv\",\n",
+ " powerOnDrag=\"../../data/calisto/powerOnDragCurve.csv\",\n",
+ ")\n",
+ "\n",
+ "Rocket1.setRailButtons([0.2, -0.5])\n",
+ "\n",
+ "NoseCone_Rocket1 = Rocket1.addNose(\n",
+ " length=0.55829, kind=\"vonKarman\", distanceToCM=0.71971\n",
+ ")\n",
+ "\n",
+ "FinSet_Rocket1 = Rocket1.addFins(\n",
+ " 4, span=0.100, rootChord=0.120, tipChord=0.040, distanceToCM=-1.04956\n",
+ ")\n",
+ "\n",
+ "Tail_Rocket1 = Rocket1.addTail(\n",
+ " topRadius=0.0635, bottomRadius=0.0435, length=0.060, distanceToCM=-1.194656\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 50,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Inertia Details\n",
+ "Rocket Dry Mass: 16.241 kg (No Propellant)\n",
+ "Rocket Total Mass: 19.196911961392022 kg (With Propellant)\n",
+ "\n",
+ "Geometrical Parameters\n",
+ "Rocket Radius: 0.0635 m\n",
+ "\n",
+ "Aerodynamics Stability\n",
+ "Initial Static Margin: 2.051 c\n",
+ "Final Static Margin: 3.090 c\n",
+ "\n",
+ "Aerodynamics Plots\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": "\n\n\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "Rocket1.info()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 51,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "RocketFlight1 = Flight(\n",
+ " rocket=Rocket1, environment=Env, inclination=85, heading=25, terminateOnApogee=True\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Start the Second Flight Stage"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now we will simulate the second flight stage, which is the landing phase of our Rocket.\n",
+ "Here we will consider that the payload was ejected at the apogee of the first stage.\n",
+ "Therefore we should be careful with the value of its mass."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 52,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "Rocket2 = Rocket(\n",
+ " motor=Pro75M1670, # This motor will not be used\n",
+ " radius=127 / 2000,\n",
+ " mass=RocketMass,\n",
+ " inertiaI=6.60,\n",
+ " inertiaZ=0.0351,\n",
+ " distanceRocketNozzle=-1.255,\n",
+ " distanceRocketPropellant=-0.85704,\n",
+ " powerOffDrag=1,\n",
+ " powerOnDrag=1,\n",
+ ")\n",
+ "\n",
+ "\n",
+ "def drogueTrigger(p, y):\n",
+ " # p = pressure\n",
+ " # y = [x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3]\n",
+ " # activate drogue when vz < 0 m/s.\n",
+ " return True if y[5] < 0 else False\n",
+ "\n",
+ "\n",
+ "def mainTrigger(p, y):\n",
+ " # p = pressure\n",
+ " # y = [x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3]\n",
+ " # activate main when vz < 0 m/s and z < 800 + 1400 m (+1400 due to surface elevation).\n",
+ " return True if y[5] < 0 and y[2] < 800 + 1400 else False\n",
+ "\n",
+ "\n",
+ "# Define Parachutes for the rocket\n",
+ "Main_Rocket2 = Rocket2.addParachute(\n",
+ " \"Main\",\n",
+ " CdS=7.2,\n",
+ " trigger=mainTrigger,\n",
+ " samplingRate=105,\n",
+ " lag=1.5,\n",
+ " noise=(0, 8.3, 0.5),\n",
+ ")\n",
+ "\n",
+ "Drogue_Rocket2 = Rocket2.addParachute(\n",
+ " \"Drogue\",\n",
+ " CdS=0.72,\n",
+ " trigger=drogueTrigger,\n",
+ " samplingRate=105,\n",
+ " lag=1.5,\n",
+ " noise=(0, 8.3, 0.5),\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 53,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Inertia Details\n",
+ "Rocket Dry Mass: 11.741 kg (No Propellant)\n",
+ "Rocket Total Mass: 14.696911961392022 kg (With Propellant)\n",
+ "\n",
+ "Geometrical Parameters\n",
+ "Rocket Radius: 0.0635 m\n",
+ "\n",
+ "Aerodynamics Stability\n",
+ "Initial Static Margin: -1.357 c\n",
+ "Final Static Margin: -0.000 c\n",
+ "\n",
+ "Main Parachute\n",
+ "CdS Coefficient: 7.2 m2\n",
+ "\n",
+ "Drogue Parachute\n",
+ "CdS Coefficient: 0.72 m2\n",
+ "\n",
+ "Aerodynamics Plots\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": "\n\n\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "Rocket2.info()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The magic line `initialSolution=RocketFlight1` will make the simulation start from the end of the first stage."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 54,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "RocketFlight2 = Flight(\n",
+ " rocket=Rocket2,\n",
+ " environment=Env,\n",
+ " inclination=0,\n",
+ " heading=0,\n",
+ " maxTime=600,\n",
+ " initialSolution=RocketFlight1,\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Simulating the 3 Flight Stage - Payload Flight"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Here we will simulate the payload flight, which is the third flight stage of our Rocket.\n",
+ "The Payload will be ejected at the apogee of the first stage.\n",
+ "Here, it will be modelled as a \"dummy\" rocket, which does not have any aerodynamic surfaces to stabilize it, neither a motor which ignites.\n",
+ "It does, however, have parachutes."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 55,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Define the \"Payload Rocket\"\n",
+ "\n",
+ "PayloadRocket = Rocket(\n",
+ " motor=Pro75M1670, # This motor will not be used\n",
+ " radius=127 / 2000,\n",
+ " mass=PayloadMass,\n",
+ " inertiaI=6.60,\n",
+ " inertiaZ=0.0351,\n",
+ " distanceRocketNozzle=-1.255,\n",
+ " distanceRocketPropellant=-0.85704,\n",
+ " powerOffDrag=0.5,\n",
+ " powerOnDrag=0.5,\n",
+ ")\n",
+ "\n",
+ "\n",
+ "def drogueTrigger(p, y):\n",
+ " # p = pressure\n",
+ " # y = [x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3]\n",
+ " # activate drogue when vz < 0 m/s.\n",
+ " return True if y[5] < 0 else False\n",
+ "\n",
+ "\n",
+ "def mainTrigger(p, y):\n",
+ " # p = pressure\n",
+ " # y = [x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3]\n",
+ " # activate main when vz < 0 m/s and z < 800 + 1400 m (+1400 due to surface elevation).\n",
+ " return True if y[5] < 0 and y[2] < 800 + 1400 else False\n",
+ "\n",
+ "\n",
+ "PayloadDrogue = PayloadRocket.addParachute(\n",
+ " \"Drogue\",\n",
+ " CdS=0.35,\n",
+ " trigger=drogueTrigger,\n",
+ " samplingRate=105,\n",
+ " lag=1.5,\n",
+ " noise=(0, 8.3, 0.5),\n",
+ ")\n",
+ "\n",
+ "PayloadMain = PayloadRocket.addParachute(\n",
+ " \"Main\",\n",
+ " CdS=4.0,\n",
+ " trigger=mainTrigger,\n",
+ " samplingRate=105,\n",
+ " lag=1.5,\n",
+ " noise=(0, 8.3, 0.5),\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The magic line `initialSolution=RocketFlight1` will make the simulation start from the end of the first stage."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 56,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "PayloadFlight = Flight(\n",
+ " rocket=PayloadRocket,\n",
+ " environment=Env,\n",
+ " inclination=0,\n",
+ " heading=0,\n",
+ " maxTime=600,\n",
+ " initialSolution=RocketFlight1,\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Plotting Everything together"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We will invoke a method from RocketPy's utilities class in order to visualize \n",
+ "the trajectory."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 57,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": "\n\n\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "utilities.compareFlightTrajectories(\n",
+ " flight_list=[RocketFlight1, RocketFlight2, PayloadFlight],\n",
+ " names=[\"Rocket - 1st Stage\", \"Rocket - 2nd Stage\", \"Payload Flight\"],\n",
+ ")"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3.10.6 64-bit",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.6"
+ },
+ "orig_nbformat": 4,
+ "vscode": {
+ "interpreter": {
+ "hash": "26de051ba29f2982a8de78e945f0abaf191376122a1563185a90213a26c5da77"
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/docs/requirements.txt b/docs/requirements.txt
index 2caba5aec..64fbb52dc 100644
--- a/docs/requirements.txt
+++ b/docs/requirements.txt
@@ -1,3 +1,3 @@
nbsphinx>=0.8.0
-pydata-sphinx-theme==0.6.3
+pydata-sphinx-theme==0.10.1
m2r2>=0.2.1
diff --git a/docs/user/index.rst b/docs/user/index.rst
index a6e1c88d6..c0274e623 100644
--- a/docs/user/index.rst
+++ b/docs/user/index.rst
@@ -11,7 +11,9 @@ Welcome to RocketPy's user documentation!
../notebooks/environment_class_usage.ipynb
../notebooks/environment_analysis_class_usage.ipynb
../notebooks/environment_analysis_EuroC_example.ipynb
+ ../notebooks/solid_motor_class_usage.ipynb
../notebooks/dispersion_analysis/dispersion_analysis.ipynb
../notebooks/utilities_usage.ipynb
+ ../notebooks/deployable_payload_example.ipynb
../matlab/matlab.rst
diff --git a/rocketpy/Flight.py b/rocketpy/Flight.py
index 8be8a2342..ad83a67c2 100644
--- a/rocketpy/Flight.py
+++ b/rocketpy/Flight.py
@@ -8,8 +8,6 @@
import math
import time
-import warnings
-from functools import cached_property
import matplotlib.pyplot as plt
import numpy as np
@@ -86,6 +84,7 @@ class Flight:
Scipy LSODA integration scheme.
State Space Vector Definition:
+ (Only available after Flight.postProcess has been called.)
Flight.x : Function
Rocket's X coordinate (positive east) as a function of time.
Flight.y : Function
@@ -202,6 +201,7 @@ class Flight:
Stores and manages flight phases.
Solution Secondary Attributes:
+ (Only available after Flight.postProcess has been called.)
Atmospheric:
Flight.windVelocityX : Function
Wind velocity X (East) experienced by the rocket as a
@@ -539,15 +539,18 @@ def __init__(
heading : int, float, optional
Heading angle relative to north given in degrees.
Default is 90, which points in the x direction.
- initialSolution : array, optional
+ initialSolution : array, Flight, optional
Initial solution array to be used. Format is
- initialSolution = [self.tInitial,
- xInit, yInit, zInit,
- vxInit, vyInit, vzInit,
- e0Init, e1Init, e2Init, e3Init,
- w1Init, w2Init, w3Init].
- If None, the initial solution will start with all null values,
- except for the euler parameters which will be calculated based
+ initialSolution = [
+ self.tInitial,
+ xInit, yInit, zInit,
+ vxInit, vyInit, vzInit,
+ e0Init, e1Init, e2Init, e3Init,
+ w1Init, w2Init, w3Init
+ ].
+ If a Flight object is used, the last state vector will be used as
+ initial solution. If None, the initial solution will start with
+ all null values, except for the euler parameters which will be calculated based
on given values of inclination and heading. Default is None.
terminateOnApogee : boolean, optional
Whether to terminate simulation when rocket reaches apogee.
@@ -605,82 +608,23 @@ def __init__(
self.terminateOnApogee = terminateOnApogee
# Modifying Rail Length for a better out of rail condition
- upperRButton = max(self.rocket.railButtons[0])
- lowerRButton = min(self.rocket.railButtons[0])
- nozzle = self.rocket.distanceRocketNozzle
- self.effective1RL = self.env.rL - abs(nozzle - upperRButton)
- self.effective2RL = self.env.rL - abs(nozzle - lowerRButton)
+ if self.rocket.railButtons is not None:
+ upperRButton = max(self.rocket.railButtons[0])
+ lowerRButton = min(self.rocket.railButtons[0])
+ nozzle = self.rocket.distanceRocketNozzle
+ self.effective1RL = self.env.rL - abs(nozzle - upperRButton)
+ self.effective2RL = self.env.rL - abs(nozzle - lowerRButton)
+ self.effective1RL = self.effective2RL = self.env.rL
# Flight initialization
self.__init_post_process_variables()
- # Initialize solution monitors
- self.outOfRailTime = 0
- self.outOfRailState = np.array([0])
- self.outOfRailVelocity = 0
- self.apogeeState = np.array([0])
- self.apogeeTime = 0
- self.apogeeX = 0
- self.apogeeY = 0
- self.apogee = 0
- self.xImpact = 0
- self.yImpact = 0
- self.impactVelocity = 0
- self.impactState = np.array([0])
- self.parachuteEvents = []
- self.postProcessed = False
- self._drift = Function(0)
- self._bearing = Function(0)
- self._latitude = Function(0)
- self._longitude = Function(0)
- # Initialize solver monitors
- self.functionEvaluations = []
- self.functionEvaluationsPerTimeStep = []
- self.timeSteps = []
- # Initialize solution state
- self.solution = []
- if self.initialSolution is None:
- # Initialize time and state variables
- self.tInitial = 0
- xInit, yInit, zInit = 0, 0, self.env.elevation
- vxInit, vyInit, vzInit = 0, 0, 0
- w1Init, w2Init, w3Init = 0, 0, 0
- # Initialize attitude
- psiInit = -heading * (np.pi / 180) # Precession / Heading Angle
- thetaInit = (inclination - 90) * (np.pi / 180) # Nutation Angle
- e0Init = np.cos(psiInit / 2) * np.cos(thetaInit / 2)
- e1Init = np.cos(psiInit / 2) * np.sin(thetaInit / 2)
- e2Init = np.sin(psiInit / 2) * np.sin(thetaInit / 2)
- e3Init = np.sin(psiInit / 2) * np.cos(thetaInit / 2)
- # Store initial conditions
- self.initialSolution = [
- self.tInitial,
- xInit,
- yInit,
- zInit,
- vxInit,
- vyInit,
- vzInit,
- e0Init,
- e1Init,
- e2Init,
- e3Init,
- w1Init,
- w2Init,
- w3Init,
- ]
- # Set initial derivative for rail phase
- self.initialDerivative = self.uDotRail1
- else:
- # Initial solution given, ignore rail phase
- # TODO: Check if rocket is actually out of rail. Otherwise, start at rail
- self.outOfRailState = self.initialSolution[1:]
- self.outOfRailTime = self.initialSolution[0]
- self.initialDerivative = self.uDot
+ self.__init_solution_monitors()
+ self.__init_flight_state()
self.tInitial = self.initialSolution[0]
self.solution.append(self.initialSolution)
self.t = self.solution[-1][0]
- self.ySol = self.solution[-1][1:]
+ self.y = self.solution[-1][1:]
# Calculate normal and lateral surface wind
windU = self.env.windVelocityX(self.env.elevation)
@@ -716,7 +660,7 @@ def __init__(
phase.solver = integrate.LSODA(
phase.derivative,
t0=phase.t,
- y0=self.ySol,
+ y0=self.y,
t_bound=phase.timeBound,
min_step=self.minTimeStep,
max_step=self.maxTimeStep,
@@ -725,7 +669,7 @@ def __init__(
)
# print('\n\tSolver Initialization Details')
# print('\tInitial Time: ', phase.t)
- # print('\tInitial State: ', self.ySol)
+ # print('\tInitial State: ', self.y)
# print('\tTime Bound: ', phase.timeBound)
# print('\tMin Step: ', self.minTimeStep)
# print('\tMax Step: ', self.maxTimeStep)
@@ -771,13 +715,13 @@ def __init__(
for parachute in node.parachutes:
# Calculate and save pressure signal
- pressure = self.env.pressure.getValueOpt(self.ySol[2])
+ pressure = self.env.pressure.getValueOpt(self.y[2])
parachute.cleanPressureSignal.append([node.t, pressure])
# Calculate and save noise
noise = parachute.noiseFunction()
parachute.noiseSignal.append([node.t, noise])
parachute.noisyPressureSignal.append([node.t, pressure + noise])
- if parachute.trigger(pressure + noise, self.ySol):
+ if parachute.trigger(pressure + noise, self.y):
# print('\nEVENT DETECTED')
# print('Parachute Triggered')
# print('Name: ', parachute.name, ' | Lag: ', parachute.lag)
@@ -811,6 +755,10 @@ def __init__(
while phase.solver.status == "running":
# Step
phase.solver.step()
+ # Normalize quaternions
+ phase.solver.y[6:10] = phase.solver.y[6:10] / sum(
+ phase.solver.y[6:10] ** 2
+ )
# Save step result
self.solution += [[phase.solver.t, *phase.solver.y]]
# Step step metrics
@@ -821,7 +769,7 @@ def __init__(
self.timeSteps.append(phase.solver.step_size)
# Update time and state
self.t = phase.solver.t
- self.ySol = phase.solver.y
+ self.y = phase.solver.y
if verbose:
print(
"Current Simulation Time: {:3.4f} s".format(self.t),
@@ -835,9 +783,9 @@ def __init__(
# Check for first out of rail event
if len(self.outOfRailState) == 1 and (
- self.ySol[0] ** 2
- + self.ySol[1] ** 2
- + (self.ySol[2] - self.env.elevation) ** 2
+ self.y[0] ** 2
+ + self.y[1] ** 2
+ + (self.y[2] - self.env.elevation) ** 2
>= self.effective1RL**2
):
# Rocket is out of rail
@@ -900,12 +848,12 @@ def __init__(
# Determine final state when upper button is going out of rail
self.t = valid_t_root[0] + self.solution[-2][0]
interpolator = phase.solver.dense_output()
- self.ySol = interpolator(self.t)
- self.solution[-1] = [self.t, *self.ySol]
+ self.y = interpolator(self.t)
+ self.solution[-1] = [self.t, *self.y]
self.outOfRailTime = self.t
- self.outOfRailState = self.ySol
+ self.outOfRailState = self.y
self.outOfRailVelocity = (
- self.ySol[3] ** 2 + self.ySol[4] ** 2 + self.ySol[5] ** 2
+ self.y[3] ** 2 + self.y[4] ** 2 + self.y[5] ** 2
) ** (0.5)
# Create new flight phase
self.flightPhases.addPhase(
@@ -917,7 +865,7 @@ def __init__(
phase.solver.status = "finished"
# Check for apogee event
- if len(self.apogeeState) == 1 and self.ySol[5] < 0:
+ if len(self.apogeeState) == 1 and self.y[5] < 0:
# print('\nPASSIVE EVENT DETECTED')
# print('Rocket Has Reached Apogee!')
# Apogee reported
@@ -950,7 +898,7 @@ def __init__(
phase.timeNodes.addNode(self.t, [], [])
phase.solver.status = "finished"
# Check for impact event
- if self.ySol[2] < self.env.elevation:
+ if self.y[2] < self.env.elevation:
# print('\nPASSIVE EVENT DETECTED')
# print('Rocket Has Reached Ground!')
# Impact reported
@@ -995,11 +943,11 @@ def __init__(
# Determine impact state at t_root
self.t = valid_t_root[0] + self.solution[-2][0]
interpolator = phase.solver.dense_output()
- self.ySol = interpolator(self.t)
+ self.y = interpolator(self.t)
# Roll back solution
- self.solution[-1] = [self.t, *self.ySol]
+ self.solution[-1] = [self.t, *self.y]
# Save impact state
- self.impactState = self.ySol
+ self.impactState = self.y
self.xImpact = self.impactState[0]
self.yImpact = self.impactState[1]
self.zImpact = self.impactState[2]
@@ -1091,7 +1039,7 @@ def __init__(
)
# Rollback history
self.t = overshootableNode.t
- self.ySol = overshootableNode.y
+ self.y = overshootableNode.y
self.solution[-1] = [
overshootableNode.t,
*overshootableNode.y,
@@ -1114,85 +1062,164 @@ def __init_post_process_variables(self):
"""Initialize post-process variables."""
# Initialize all variables created during Flight.postProcess()
# Important to do so that MATLAB® can access them
- self._windVelocityX = Function(0)
- self._windVelocityY = Function(0)
- self._density = Function(0)
- self._pressure = Function(0)
- self._dynamicViscosity = Function(0)
- self._speedOfSound = Function(0)
- self._ax = Function(0)
- self._ay = Function(0)
- self._az = Function(0)
- self._alpha1 = Function(0)
- self._alpha2 = Function(0)
- self._alpha3 = Function(0)
- self._speed = Function(0)
- self._maxSpeed = 0
- self._maxSpeedTime = 0
- self._horizontalSpeed = Function(0)
- self._Acceleration = Function(0)
- self._maxAcceleration = 0
- self._maxAccelerationTime = 0
- self._pathAngle = Function(0)
- self._attitudeVectorX = Function(0)
- self._attitudeVectorY = Function(0)
- self._attitudeVectorZ = Function(0)
- self._attitudeAngle = Function(0)
- self._lateralAttitudeAngle = Function(0)
- self._phi = Function(0)
- self._theta = Function(0)
- self._psi = Function(0)
- self._R1 = Function(0)
- self._R2 = Function(0)
- self._R3 = Function(0)
- self._M1 = Function(0)
- self._M2 = Function(0)
- self._M3 = Function(0)
- self._aerodynamicLift = Function(0)
- self._aerodynamicDrag = Function(0)
- self._aerodynamicBendingMoment = Function(0)
- self._aerodynamicSpinMoment = Function(0)
- self._railButton1NormalForce = Function(0)
- self._maxRailButton1NormalForce = 0
- self._railButton1ShearForce = Function(0)
- self._maxRailButton1ShearForce = 0
- self._railButton2NormalForce = Function(0)
- self._maxRailButton2NormalForce = 0
- self._railButton2ShearForce = Function(0)
- self._maxRailButton2ShearForce = 0
- self._rotationalEnergy = Function(0)
- self._translationalEnergy = Function(0)
- self._kineticEnergy = Function(0)
- self._potentialEnergy = Function(0)
- self._totalEnergy = Function(0)
- self._thrustPower = Function(0)
- self._dragPower = Function(0)
- self._attitudeFrequencyResponse = Function(0)
- self._omega1FrequencyResponse = Function(0)
- self._omega2FrequencyResponse = Function(0)
- self._omega3FrequencyResponse = Function(0)
- self._streamVelocityX = Function(0)
- self._streamVelocityY = Function(0)
- self._streamVelocityZ = Function(0)
- self._freestreamSpeed = Function(0)
- self._apogeeFreestreamSpeed = 0
- self._MachNumber = Function(0)
- self._maxMachNumber = 0
- self._maxMachNumberTime = 0
- self._ReynoldsNumber = Function(0)
- self._maxReynoldsNumber = 0
- self._maxReynoldsNumberTime = 0
- self._dynamicPressure = Function(0)
- self._maxDynamicPressure = 0
- self._maxDynamicPressureTime = 0
- self._totalPressure = Function(0)
- self._maxTotalPressure = 0
- self._maxTotalPressureTime = 0
- self._angleOfAttack = Function(0)
+ self.windVelocityX = Function(0)
+ self.windVelocityY = Function(0)
+ self.density = Function(0)
+ self.pressure = Function(0)
+ self.dynamicViscosity = Function(0)
+ self.speedOfSound = Function(0)
+ self.ax = Function(0)
+ self.ay = Function(0)
+ self.az = Function(0)
+ self.alpha1 = Function(0)
+ self.alpha2 = Function(0)
+ self.alpha3 = Function(0)
+ self.speed = Function(0)
+ self.maxSpeed = 0
+ self.maxSpeedTime = 0
+ self.horizontalSpeed = Function(0)
+ self.Acceleration = Function(0)
+ self.maxAcceleration = 0
+ self.maxAccelerationTime = 0
+ self.pathAngle = Function(0)
+ self.attitudeVectorX = Function(0)
+ self.attitudeVectorY = Function(0)
+ self.attitudeVectorZ = Function(0)
+ self.attitudeAngle = Function(0)
+ self.lateralAttitudeAngle = Function(0)
+ self.phi = Function(0)
+ self.theta = Function(0)
+ self.psi = Function(0)
+ self.R1 = Function(0)
+ self.R2 = Function(0)
+ self.R3 = Function(0)
+ self.M1 = Function(0)
+ self.M2 = Function(0)
+ self.M3 = Function(0)
+ self.aerodynamicLift = Function(0)
+ self.aerodynamicDrag = Function(0)
+ self.aerodynamicBendingMoment = Function(0)
+ self.aerodynamicSpinMoment = Function(0)
+ self.railButton1NormalForce = Function(0)
+ self.maxRailButton1NormalForce = 0
+ self.railButton1ShearForce = Function(0)
+ self.maxRailButton1ShearForce = 0
+ self.railButton2NormalForce = Function(0)
+ self.maxRailButton2NormalForce = 0
+ self.railButton2ShearForce = Function(0)
+ self.maxRailButton2ShearForce = 0
+ self.rotationalEnergy = Function(0)
+ self.translationalEnergy = Function(0)
+ self.kineticEnergy = Function(0)
+ self.potentialEnergy = Function(0)
+ self.totalEnergy = Function(0)
+ self.thrustPower = Function(0)
+ self.dragPower = Function(0)
+ self.attitudeFrequencyResponse = Function(0)
+ self.omega1FrequencyResponse = Function(0)
+ self.omega2FrequencyResponse = Function(0)
+ self.omega3FrequencyResponse = Function(0)
+ self.streamVelocityX = Function(0)
+ self.streamVelocityY = Function(0)
+ self.streamVelocityZ = Function(0)
+ self.freestreamSpeed = Function(0)
+ self.apogeeFreestreamSpeed = 0
+ self.MachNumber = Function(0)
+ self.maxMachNumber = 0
+ self.maxMachNumberTime = 0
+ self.ReynoldsNumber = Function(0)
+ self.maxReynoldsNumber = 0
+ self.maxReynoldsNumberTime = 0
+ self.dynamicPressure = Function(0)
+ self.maxDynamicPressure = 0
+ self.maxDynamicPressureTime = 0
+ self.totalPressure = Function(0)
+ self.maxTotalPressure = 0
+ self.maxTotalPressureTime = 0
+ self.angleOfAttack = Function(0)
self.flutterMachNumber = Function(0)
self.difference = Function(0)
self.safetyFactor = Function(0)
+ def __init_solution_monitors(self):
+ """Initialize all variables related to solution monitoring."""
+ # Initialize trajectory monitors
+ self.outOfRailTime = 0
+ self.outOfRailState = np.array([0])
+ self.outOfRailVelocity = 0
+ self.apogeeState = np.array([0])
+ self.apogeeTime = 0
+ self.apogeeX = 0
+ self.apogeeY = 0
+ self.apogee = 0
+ self.xImpact = 0
+ self.yImpact = 0
+ self.impactVelocity = 0
+ self.impactState = np.array([0])
+ self.parachuteEvents = []
+ self.postProcessed = False
+ self.latitude = 0 # Function(0)
+ self.longitude = 0 # Function(0)
+ # Initialize solver monitors
+ self.functionEvaluations = []
+ self.functionEvaluationsPerTimeStep = []
+ self.timeSteps = []
+ # Initialize solution state
+ self.solution = []
+
+ # TODO: Need unit tests
+ def __init_flight_state(self):
+ """Initialize flight state."""
+ if self.initialSolution is None:
+ # Initialize time and state variables based on heading and inclination
+ self.tInitial = 0
+ xInit, yInit, zInit = 0, 0, self.env.elevation
+ vxInit, vyInit, vzInit = 0, 0, 0
+ w1Init, w2Init, w3Init = 0, 0, 0
+ # Initialize attitude
+ psiInit = -self.heading * (np.pi / 180) # Precession / self. Angle
+ thetaInit = (self.inclination - 90) * (np.pi / 180) # Nutation Angle
+ e0Init = np.cos(psiInit / 2) * np.cos(thetaInit / 2)
+ e1Init = np.cos(psiInit / 2) * np.sin(thetaInit / 2)
+ e2Init = np.sin(psiInit / 2) * np.sin(thetaInit / 2)
+ e3Init = np.sin(psiInit / 2) * np.cos(thetaInit / 2)
+ # Store initial conditions
+ self.initialSolution = [
+ self.tInitial,
+ xInit,
+ yInit,
+ zInit,
+ vxInit,
+ vyInit,
+ vzInit,
+ e0Init,
+ e1Init,
+ e2Init,
+ e3Init,
+ w1Init,
+ w2Init,
+ w3Init,
+ ]
+ # Set initial derivative for rail phase
+ self.initialDerivative = self.uDotRail1
+ # TODO: Need unit tests
+ elif isinstance(self.initialSolution, Flight):
+ # TODO: Add optional argument to specify wether to merge/concatenate flight or not
+ # Initialize time and state variables based on last solution of
+ # previous flight
+ self.initialSolution = self.initialSolution.solution[-1]
+ # Set unused monitors
+ self.outOfRailState = self.initialSolution[1:]
+ self.outOfRailTime = self.initialSolution[0]
+ # Set initial derivative for 6-DOF flight phase
+ self.initialDerivative = self.uDot
+ else:
+ # Initial solution given, ignore rail phase
+ # TODO: Check if rocket is actually out of rail. Otherwise, start at rail
+ self.outOfRailState = self.initialSolution[1:]
+ self.outOfRailTime = self.initialSolution[0]
+ self.initialDerivative = self.uDot
+
def uDotRail1(self, t, u, postProcessing=False):
"""Calculates derivative of u state vector with respect to time
when rocket is flying in 1 DOF motion in the rail.
@@ -1511,21 +1538,19 @@ def uDot(self, t, u, postProcessing=False):
if postProcessing:
# 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])
+ self.R1.append([t, R1])
+ self.R2.append([t, R2])
+ self.R3.append([t, R3])
+ self.M1.append([t, M1])
+ self.M2.append([t, M2])
+ self.M3.append([t, M3])
# Atmospheric Conditions
- self.windVelocityX_list.append([t, self.env.windVelocityX.getValueOpt(z)])
- self.windVelocityY_list.append([t, self.env.windVelocityY.getValueOpt(z)])
- self.density_list.append([t, self.env.density.getValueOpt(z)])
- self.dynamicViscosity_list.append(
- [t, self.env.dynamicViscosity.getValueOpt(z)]
- )
- self.pressure_list.append([t, self.env.pressure.getValueOpt(z)])
- self.speedOfSound_list.append([t, self.env.speedOfSound.getValueOpt(z)])
+ self.windVelocityX.append([t, self.env.windVelocityX(z)])
+ self.windVelocityY.append([t, self.env.windVelocityY(z)])
+ self.density.append([t, self.env.density(z)])
+ self.dynamicViscosity.append([t, self.env.dynamicViscosity(z)])
+ self.pressure.append([t, self.env.pressure(z)])
+ self.speedOfSound.append([t, self.env.speedOfSound(z)])
return uDot
@@ -1589,713 +1614,206 @@ def uDotParachute(self, t, u, postProcessing=False):
if postProcessing:
# 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])
+ self.R1.append([t, Dx])
+ self.R2.append([t, Dy])
+ self.R3.append([t, Dz])
+ self.M1.append([t, 0])
+ self.M2.append([t, 0])
+ self.M3.append([t, 0])
# Atmospheric Conditions
- self.windVelocityX_list.append([t, self.env.windVelocityX(z)])
- self.windVelocityY_list.append([t, self.env.windVelocityY(z)])
- self.density_list.append([t, self.env.density(z)])
- self.dynamicViscosity_list.append([t, self.env.dynamicViscosity(z)])
- self.pressure_list.append([t, self.env.pressure(z)])
- self.speedOfSound_list.append([t, self.env.speedOfSound(z)])
+ self.windVelocityX.append([t, self.env.windVelocityX(z)])
+ self.windVelocityY.append([t, self.env.windVelocityY(z)])
+ self.density.append([t, self.env.density(z)])
+ self.dynamicViscosity.append([t, self.env.dynamicViscosity(z)])
+ self.pressure.append([t, self.env.pressure(z)])
+ self.speedOfSound.append([t, self.env.speedOfSound(z)])
return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0]
- # Process first type of outputs - state vector
- # Transform solution array into Functions
- @cached_property
- def x(self, interpolation="spline", extrapolation="natural"):
- """Rocket x position as a function of time.
+ def postProcess(self, interpolation="spline", extrapolation="natural"):
+ """Post-process all Flight information produced during
+ simulation. Includes the calculation of maximum values,
+ calculation of secondary values such as energy and conversion
+ of lists to Function objects to facilitate plotting.
Parameters
----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps
- interpolation, 'constant', which returns the value of the function
- at the edge of the interval, and 'zero', which returns zero for all
- points outside of source range. Default is 'natural'.
+ interpolation: string
+ Interpolation method to be used in the Function objects.
+ extrapolation: string
+ Extrapolation method to be used in the Function objects.
- Returns
- -------
- x: Function
- Rocket x position as a function of time.
+ Return
+ ------
+ None
"""
- return Function(
- np.array(self.solution)[:, [0, 1]],
- "Time (s)",
- "X (m)",
- interpolation,
- extrapolation,
+ # Process first type of outputs - state vector
+ # Transform solution array into Functions
+ sol = np.array(self.solution)
+ self.x = Function(
+ sol[:, [0, 1]], "Time (s)", "X (m)", interpolation, extrapolation
)
-
- @cached_property
- def y(self, interpolation="spline", extrapolation="natural"):
- """ocket y position as a function of time.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps
- interpolation, 'constant', which returns the value of the function
- at the edge of the interval, and 'zero', which returns zero for all
- points outside of source range. Default is 'natural'.
-
- Returns
- -------
- y: Function
- Rocket y position as a function of time.
- """
- return Function(
- np.array(self.solution)[:, [0, 2]],
- "Time (s)",
- "Y (m)",
- interpolation,
- extrapolation,
+ self.y = Function(
+ sol[:, [0, 2]], "Time (s)", "Y (m)", interpolation, extrapolation
)
-
- @cached_property
- def z(self, interpolation="spline", extrapolation="natural"):
- """Rocket z position as a function of time.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps
- interpolation, 'constant', which returns the value of the function
- at the edge of the interval, and 'zero', which returns zero for all
- points outside of source range. Default is 'natural'.
-
- Returns
- -------
- z: Function
- Rocket z position as a function of time.
- """
- return Function(
- np.array(self.solution)[:, [0, 3]],
- "Time (s)",
- "Z (m)",
- interpolation,
- extrapolation,
+ self.z = Function(
+ sol[:, [0, 3]], "Time (s)", "Z (m)", interpolation, extrapolation
)
-
- @cached_property
- def vx(self, interpolation="spline", extrapolation="natural"):
- """Rocket x velocity as a function of time.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps
- interpolation, 'constant', which returns the value of the function
- at the edge of the interval, and 'zero', which returns zero for all
- points outside of source range. Default is 'natural'.
-
- Returns
- -------
- vx: Function
- Rocket x velocity as a function of time.
- """
-
- return Function(
- np.array(self.solution)[:, [0, 4]],
- "Time (s)",
- "Vx (m/s)",
- interpolation,
- extrapolation,
+ self.vx = Function(
+ sol[:, [0, 4]], "Time (s)", "Vx (m/s)", interpolation, extrapolation
)
-
- @cached_property
- def vy(self, interpolation="spline", extrapolation="natural"):
- """Rocket y velocity as a function of time.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps
- interpolation, 'constant', which returns the value of the function
- at the edge of the interval, and 'zero', which returns zero for all
- points outside of source range. Default is 'natural'.
-
- Returns
- -------
- vy: Function
- Rocket y velocity as a function of time.
- """
- return Function(
- np.array(self.solution)[:, [0, 5]],
- "Time (s)",
- "Vy (m/s)",
- interpolation,
- extrapolation,
+ self.vy = Function(
+ sol[:, [0, 5]], "Time (s)", "Vy (m/s)", interpolation, extrapolation
)
-
- @cached_property
- def vz(self, interpolation="spline", extrapolation="natural"):
- """Rocket z velocity as a function of time.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps
- interpolation, 'constant', which returns the value of the function
- at the edge of the interval, and 'zero', which returns zero for all
- points outside of source range. Default is 'natural'.
-
- Returns
- -------
- vz: Function
- Rocket z velocity as a function of time.
- """
- return Function(
- np.array(self.solution)[:, [0, 6]],
- "Time (s)",
- "Vz (m/s)",
- interpolation,
- extrapolation,
- )
-
- @cached_property
- def e0(self, interpolation="spline", extrapolation="natural"):
- return Function(
- np.array(self.solution)[:, [0, 7]],
- "Time (s)",
- "e0",
- interpolation,
- extrapolation,
+ self.vz = Function(
+ sol[:, [0, 6]], "Time (s)", "Vz (m/s)", interpolation, extrapolation
)
-
- @cached_property
- def e1(self, interpolation="spline", extrapolation="natural"):
- return Function(
- np.array(self.solution)[:, [0, 8]],
- "Time (s)",
- "e1",
- interpolation,
- extrapolation,
- )
-
- @cached_property
- def e2(self, interpolation="spline", extrapolation="natural"):
- return Function(
- np.array(self.solution)[:, [0, 9]],
- "Time (s)",
- "e2",
- interpolation,
- extrapolation,
+ self.e0 = Function(
+ sol[:, [0, 7]], "Time (s)", "e0", interpolation, extrapolation
)
-
- @cached_property
- def e3(self, interpolation="spline", extrapolation="natural"):
- return Function(
- np.array(self.solution)[:, [0, 10]],
- "Time (s)",
- "e3",
- interpolation,
- extrapolation,
+ self.e1 = Function(
+ sol[:, [0, 8]], "Time (s)", "e1", interpolation, extrapolation
)
-
- @cached_property
- def w1(self, interpolation="spline", extrapolation="natural"):
- return Function(
- np.array(self.solution)[:, [0, 11]],
- "Time (s)",
- "ω1 (rad/s)",
- interpolation,
- extrapolation,
+ self.e2 = Function(
+ sol[:, [0, 9]], "Time (s)", "e2", interpolation, extrapolation
)
-
- @cached_property
- def w2(self, interpolation="spline", extrapolation="natural"):
- return Function(
- np.array(self.solution)[:, [0, 12]],
- "Time (s)",
- "ω2 (rad/s)",
- interpolation,
- extrapolation,
+ self.e3 = Function(
+ sol[:, [0, 10]], "Time (s)", "e3", interpolation, extrapolation
)
-
- @cached_property
- def w3(self, interpolation="spline", extrapolation="natural"):
- return Function(
- np.array(self.solution)[:, [0, 13]],
- "Time (s)",
- "ω3 (rad/s)",
- interpolation,
- extrapolation,
+ self.w1 = Function(
+ sol[:, [0, 11]], "Time (s)", "ω1 (rad/s)", interpolation, extrapolation
)
-
- # Process second type of outputs - accelerations components
- @cached_property
- def ax(self, interpolation="spline", extrapolation="natural"):
- ax = self.retrieve_acceleration_arrays[0]
- # Convert accelerations to functions
- ax = Function(ax, "Time (s)", "Ax (m/s2)", interpolation, extrapolation)
- return ax
-
- @cached_property
- def ay(self, interpolation="spline", extrapolation="natural"):
- ay = self.retrieve_acceleration_arrays[1]
- # Convert accelerations to functions
- ay = Function(ay, "Time (s)", "Ay (m/s2)", interpolation, extrapolation)
- return ay
-
- @cached_property
- def az(self, interpolation="spline", extrapolation="natural"):
- az = self.retrieve_acceleration_arrays[2]
- # Convert accelerations to functions
- az = Function(az, "Time (s)", "Az (m/s2)", interpolation, extrapolation)
- return az
-
- @cached_property
- def alpha1(self, interpolation="spline", extrapolation="natural"):
- alpha1 = self.retrieve_acceleration_arrays[3]
- # Convert accelerations to functions
- alpha1 = Function(
- alpha1, "Time (s)", "α1 (rad/s2)", interpolation, extrapolation
+ self.w2 = Function(
+ sol[:, [0, 12]], "Time (s)", "ω2 (rad/s)", interpolation, extrapolation
)
- return alpha1
-
- @cached_property
- def alpha2(self, interpolation="spline", extrapolation="natural"):
- alpha2 = self.retrieve_acceleration_arrays[4]
- # Convert accelerations to functions
- alpha2 = Function(
- alpha2, "Time (s)", "α2 (rad/s2)", interpolation, extrapolation
+ self.w3 = Function(
+ sol[:, [0, 13]], "Time (s)", "ω3 (rad/s)", interpolation, extrapolation
)
- return alpha2
- @cached_property
- def alpha3(self, interpolation="spline", extrapolation="natural"):
- alpha3 = self.retrieve_acceleration_arrays[5]
+ # Process second type of outputs - accelerations
+ # Initialize acceleration arrays
+ self.ax, self.ay, self.az = [], [], []
+ self.alpha1, self.alpha2, self.alpha3 = [], [], []
+ # Go through each time step and calculate accelerations
+ # Get flight phases
+ for phase_index, phase in self.timeIterator(self.flightPhases):
+ initTime = phase.t
+ finalTime = self.flightPhases[phase_index + 1].t
+ currentDerivative = 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 initTime < step[0] <= finalTime or (
+ initTime == self.tInitial and step[0] == self.tInitial
+ ):
+ # Get derivatives
+ uDot = currentDerivative(step[0], step[1:])
+ # Get accelerations
+ ax, ay, az = uDot[3:6]
+ alpha1, alpha2, alpha3 = uDot[10:]
+ # Save accelerations
+ self.ax.append([step[0], ax])
+ self.ay.append([step[0], ay])
+ self.az.append([step[0], az])
+ self.alpha1.append([step[0], alpha1])
+ self.alpha2.append([step[0], alpha2])
+ self.alpha3.append([step[0], alpha3])
# Convert accelerations to functions
- alpha3 = Function(
- alpha3, "Time (s)", "α3 (rad/s2)", interpolation, extrapolation
- )
- return alpha3
-
- # Process third type of outputs - Temporary values
- @cached_property
- def R1(self, interpolation="spline", extrapolation="natural"):
- """Aerodynamic force along the first axis that is perpendicular to the
- rocket's axis of symmetry.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps interpolation, 'constant',
- which returns the value of the function at the edge of the interval,
- and 'zero', which returns zero for all points outside of source
- range. Default is 'natural'.
-
- Returns
- -------
- R1: Function
- Aero force along the first axis that is perpendicular to the
- rocket's axis of symmetry.
- """
- R1 = self.retrieve_temporary_values_arrays[0]
- R1 = Function(R1, "Time (s)", "R1 (m)", interpolation, extrapolation)
-
- return R1
-
- @cached_property
- def R2(self, interpolation="spline", extrapolation="natural"):
- """Aerodynamic force along the second axis that is perpendicular to the
- rocket's axis of symmetry.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps interpolation, 'constant',
- which returns the value of the function at the edge of the interval,
- and 'zero', which returns zero for all points outside of source
- range. Default is 'natural'.
-
- Returns
- -------
- R2: Function
- Aero force along the second axis that is perpendicular to the
- rocket's axis of symmetry.
- """
- R2 = self.retrieve_temporary_values_arrays[1]
- R2 = Function(R2, "Time (s)", "R2 (m)", interpolation, extrapolation)
-
- return R2
-
- @cached_property
- def R3(self, interpolation="spline", extrapolation="natural"):
- """Aerodynamic force along the rocket's axis of symmetry.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps interpolation, 'constant',
- which returns the value of the function at the edge of the interval,
- and 'zero', which returns zero for all points outside of source
- range. Default is 'natural'.
-
- Returns
- -------
- R3: Function
- Aerodynamic force along the rocket's axis of symmetry.
- """
- R3 = self.retrieve_temporary_values_arrays[2]
- R3 = Function(R3, "Time (s)", "R3 (m)", interpolation, extrapolation)
-
- return R3
-
- @cached_property
- def M1(self, interpolation="spline", extrapolation="natural"):
- """Aerodynamic bending moment in the same direction as the axis that is
- perpendicular to the rocket's axis of symmetry.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps interpolation, 'constant',
- which returns the value of the function at the edge of the interval,
- and 'zero', which returns zero for all points outside of source
- range. Default is 'natural'.
-
- Returns
- -------
- M1: Function
- Aero moment along the first axis that is perpendicular to the
- rocket's axis of symmetry.
- """
- M1 = self.retrieve_temporary_values_arrays[3]
- M1 = Function(M1, "Time (s)", "M1 (Nm)", interpolation, extrapolation)
-
- return M1
-
- @cached_property
- def M2(self, interpolation="spline", extrapolation="natural"):
- """Aerodynamic moment in the same direction as the second axis that is
- perpendicular to the rocket's axis of symmetry.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps interpolation, 'constant',
- which returns the value of the function at the edge of the interval,
- and 'zero', which returns zero for all points outside of source
- range. Default is 'natural'.
-
- Returns
- -------
- M2: Function
- Aero moment along the second axis that is perpendicular to the
- rocket's axis of symmetry.
- """
- M2 = self.retrieve_temporary_values_arrays[4]
- M2 = Function(M2, "Time (s)", "M2 (Nm)", interpolation, extrapolation)
-
- return M2
-
- @cached_property
- def M3(self, interpolation="spline", extrapolation="natural"):
- """Aerodynamic bending in the rocket's axis of symmetry direction.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps interpolation, 'constant',
- which returns the value of the function at the edge of the interval,
- and 'zero', which returns zero for all points outside of source
- range. Default is 'natural'.
-
- Returns
- -------
- M3: Function
- Aero moment in the same direction as the rocket's axis of symmetry.
- """
- M3 = self.retrieve_temporary_values_arrays[5]
- M3 = Function(M3, "Time (s)", "M3 (Nm)", interpolation, extrapolation)
-
- return M3
-
- @cached_property
- def pressure(self, interpolation="spline", extrapolation="natural"):
- """Air pressure at each time step.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps interpolation, 'constant',
- which returns the value of the function at the edge of the interval,
- and 'zero', which returns zero for all points outside of source
- range. Default is 'natural'.
-
- Returns
- -------
- pressure: Function
- Air pressure at each time step.
- """
- pressure = self.retrieve_temporary_values_arrays[6]
- pressure = Function(
- pressure, "Time (s)", "Pressure (Pa)", interpolation, extrapolation
- )
-
- return pressure
-
- @cached_property
- def density(self, interpolation="spline", extrapolation="natural"):
- """Air density at each time step.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps interpolation, 'constant',
- which returns the value of the function at the edge of the interval,
- and 'zero', which returns zero for all points outside of source
- range. Default is 'natural'.
-
- Returns
- -------
- density: Function
- Air density at each time step.
- """
- density = self.retrieve_temporary_values_arrays[7]
- density = Function(
- density, "Time (s)", "Density (kg/m³)", interpolation, extrapolation
- )
-
- return density
-
- @cached_property
- def dynamicViscosity(self, interpolation="spline", extrapolation="natural"):
- """Air dynamic viscosity at each time step.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps interpolation, 'constant',
- which returns the value of the function at the edge of the interval,
- and 'zero', which returns zero for all points outside of source
- range. Default is 'natural'.
-
- Returns
- -------
- dynamicViscosity: Function
- Air dynamic viscosity at each time step.
- """
- dynamicViscosity = self.retrieve_temporary_values_arrays[8]
- dynamicViscosity = Function(
- dynamicViscosity,
- "Time (s)",
- "Dynamic Viscosity (Pa s)",
- interpolation,
- extrapolation,
- )
-
- return dynamicViscosity
-
- @cached_property
- def speedOfSound(self, interpolation="spline", extrapolation="natural"):
- """Speed of sound at each time step.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps interpolation, 'constant',
- which returns the value of the function at the edge of the interval,
- and 'zero', which returns zero for all points outside of source
- range. Default is 'natural'.
-
- Returns
- -------
- speedOfSound: Function
- Speed of sound at each time step.
- """
- speedOfSound = self.retrieve_temporary_values_arrays[9]
- speedOfSound = Function(
- speedOfSound,
- "Time (s)",
- "Speed of Sound (m/s)",
- interpolation,
- extrapolation,
- )
-
- return speedOfSound
-
- @cached_property
- def windVelocityX(self, interpolation="spline", extrapolation="natural"):
- """Wind velocity in the X direction at each time step.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps interpolation, 'constant',
- which returns the value of the function at the edge of the interval,
- and 'zero', which returns zero for all points outside of source
- range. Default is 'natural'.
-
- Returns
- -------
- windVelocityX: Function
- Wind velocity in the X direction at each time step.
- """
- windVelocityX = self.retrieve_temporary_values_arrays[10]
- windVelocityX = Function(
- windVelocityX,
+ self.ax = Function(self.ax, "Time (s)", "Ax (m/s2)", interpolation)
+ self.ay = Function(self.ay, "Time (s)", "Ay (m/s2)", interpolation)
+ self.az = Function(self.az, "Time (s)", "Az (m/s2)", interpolation)
+ self.alpha1 = Function(self.alpha1, "Time (s)", "α1 (rad/s2)", interpolation)
+ self.alpha2 = Function(self.alpha2, "Time (s)", "α2 (rad/s2)", interpolation)
+ self.alpha3 = Function(self.alpha3, "Time (s)", "α3 (rad/s2)", interpolation)
+
+ # Process third type of outputs - temporary values calculated during integration
+ # Initialize force and atmospheric arrays
+ self.R1, self.R2, self.R3, self.M1, self.M2, self.M3 = [], [], [], [], [], []
+ self.pressure, self.density, self.dynamicViscosity, self.speedOfSound = (
+ [],
+ [],
+ [],
+ [],
+ )
+ self.windVelocityX, self.windVelocityY = [], []
+ # Go through each time step and calculate forces and atmospheric values
+ # Get flight phases
+ for phase_index, phase in self.timeIterator(self.flightPhases):
+ initTime = phase.t
+ finalTime = self.flightPhases[phase_index + 1].t
+ currentDerivative = 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 initTime < step[0] <= finalTime or (
+ initTime == self.tInitial and step[0] == self.tInitial
+ ):
+ # Call derivatives in post processing mode
+ uDot = currentDerivative(step[0], step[1:], postProcessing=True)
+ # Convert forces and atmospheric arrays to functions
+ self.R1 = Function(self.R1, "Time (s)", "R1 (N)", interpolation)
+ self.R2 = Function(self.R2, "Time (s)", "R2 (N)", interpolation)
+ self.R3 = Function(self.R3, "Time (s)", "R3 (N)", interpolation)
+ self.M1 = Function(self.M1, "Time (s)", "M1 (Nm)", interpolation)
+ self.M2 = Function(self.M2, "Time (s)", "M2 (Nm)", interpolation)
+ self.M3 = Function(self.M3, "Time (s)", "M3 (Nm)", interpolation)
+ self.windVelocityX = Function(
+ self.windVelocityX,
"Time (s)",
"Wind Velocity X (East) (m/s)",
interpolation,
- extrapolation,
)
-
- return windVelocityX
-
- @cached_property
- def windVelocityY(self, interpolation="spline", extrapolation="natural"):
- """Wind velocity in the Y direction at each time step.
-
- Parameters
- ----------
- extrapolation : str, optional
- Function extrapolation mode. Options are 'linear', 'polynomial',
- 'akima' and 'spline'. Default is 'spline'.
- interpolation : str, optional
- Function extrapolation mode. Options are 'natural', which keeps interpolation, 'constant',
- which returns the value of the function at the edge of the interval,
- and 'zero', which returns zero for all points outside of source
- range. Default is 'natural'.
-
- Returns
- -------
- windVelocityY: Function
- Wind velocity in the Y direction at each time step.
- """
- windVelocityY = self.retrieve_temporary_values_arrays[10]
- windVelocityY = Function(
- windVelocityY,
+ self.windVelocityY = Function(
+ self.windVelocityY,
"Time (s)",
"Wind Velocity Y (North) (m/s)",
interpolation,
- extrapolation,
+ )
+ self.density = Function(
+ self.density, "Time (s)", "Density (kg/m³)", interpolation
+ )
+ self.pressure = Function(
+ self.pressure, "Time (s)", "Pressure (Pa)", interpolation
+ )
+ self.dynamicViscosity = Function(
+ self.dynamicViscosity, "Time (s)", "Dynamic Viscosity (Pa s)", interpolation
+ )
+ self.speedOfSound = Function(
+ self.speedOfSound, "Time (s)", "Speed of Sound (m/s)", interpolation
)
- return windVelocityY
-
- # Process fourth type of output - values calculated from previous outputs
-
- # Kinematics functions and values
- # Velocity Magnitude
- @cached_property
- def speed(self):
- speed = (self.vx**2 + self.vy**2 + self.vz**2) ** 0.5
- speed.setOutputs("Speed - Velocity Magnitude (m/s)")
- return speed
+ # Process fourth type of output - values calculated from previous outputs
- @cached_property
- def maxSpeedTime(self):
+ # Kinematics functions and values
+ # Velocity Magnitude
+ self.speed = (self.vx**2 + self.vy**2 + self.vz**2) ** 0.5
+ self.speed.setOutputs("Speed - Velocity Magnitude (m/s)")
maxSpeedTimeIndex = np.argmax(self.speed[:, 1])
- return self.speed[maxSpeedTimeIndex, 0]
-
- @cached_property
- def maxSpeed(self):
- return self.speed(self.maxSpeedTime)
-
- # Accelerations
- @cached_property
- def acceleration(self):
- acceleration = (self.ax**2 + self.ay**2 + self.az**2) ** 0.5
- acceleration.setOutputs("Acceleration Magnitude (m/s²)")
- return acceleration
-
- @cached_property
- def maxAcceleration(self):
+ self.maxSpeed = self.speed[maxSpeedTimeIndex, 1]
+ self.maxSpeedTime = self.speed[maxSpeedTimeIndex, 0]
+ # Acceleration
+ self.acceleration = (self.ax**2 + self.ay**2 + self.az**2) ** 0.5
+ self.acceleration.setOutputs("Acceleration Magnitude (m/s²)")
maxAccelerationTimeIndex = np.argmax(self.acceleration[:, 1])
- return self.acceleration[maxAccelerationTimeIndex, 1]
-
- @cached_property
- def maxAccelerationTime(self):
- maxAccelerationTimeIndex = np.argmax(self.acceleration[:, 1])
- return self.acceleration[maxAccelerationTimeIndex, 0]
-
- # Path Angle
- @cached_property
- def horizontalSpeed(self):
- return (self.vx**2 + self.vy**2) ** 0.5
-
- @cached_property
- def pathAngle(self):
+ self.maxAcceleration = self.acceleration[maxAccelerationTimeIndex, 1]
+ self.maxAccelerationTime = self.acceleration[maxAccelerationTimeIndex, 0]
+ # Path Angle
+ self.horizontalSpeed = (self.vx**2 + self.vy**2) ** 0.5
pathAngle = (180 / np.pi) * np.arctan2(
self.vz[:, 1], self.horizontalSpeed[:, 1]
)
pathAngle = np.column_stack([self.vz[:, 0], pathAngle])
- return Function(pathAngle, "Time (s)", "Path Angle (°)")
-
- # Attitude Angle
- @cached_property
- def attitudeVectorX(self):
- return 2 * (self.e1 * self.e3 + self.e0 * self.e2) # a13
-
- @cached_property
- def attitudeVectorY(self):
- return 2 * (self.e2 * self.e3 - self.e0 * self.e1) # a23
-
- @cached_property
- def attitudeVectorZ(self):
- return 1 - 2 * (self.e1**2 + self.e2**2) # a33
-
- @cached_property
- def attitudeAngle(self):
+ self.pathAngle = Function(
+ pathAngle, "Time (s)", "Path Angle (°)", interpolation
+ )
+ # Attitude Angle
+ self.attitudeVectorX = 2 * (self.e1 * self.e3 + self.e0 * self.e2) # a13
+ self.attitudeVectorY = 2 * (self.e2 * self.e3 - self.e0 * self.e1) # a23
+ self.attitudeVectorZ = 1 - 2 * (self.e1**2 + self.e2**2) # a33
horizontalAttitudeProj = (
self.attitudeVectorX**2 + self.attitudeVectorY**2
) ** 0.5
@@ -2303,11 +1821,10 @@ def attitudeAngle(self):
self.attitudeVectorZ[:, 1], horizontalAttitudeProj[:, 1]
)
attitudeAngle = np.column_stack([self.attitudeVectorZ[:, 0], attitudeAngle])
- return Function(attitudeAngle, "Time (s)", "Attitude Angle (°)")
-
- # Lateral Attitude Angle
- @cached_property
- def lateralAttitudeAngle(self):
+ self.attitudeAngle = Function(
+ attitudeAngle, "Time (s)", "Attitude Angle (°)", interpolation
+ )
+ # Lateral Attitude Angle
lateralVectorAngle = (np.pi / 180) * (self.heading - 90)
lateralVectorX = np.sin(lateralVectorAngle)
lateralVectorY = np.cos(lateralVectorAngle)
@@ -2331,333 +1848,175 @@ def lateralAttitudeAngle(self):
lateralAttitudeAngle = np.column_stack(
[self.attitudeVectorZ[:, 0], lateralAttitudeAngle]
)
- return Function(lateralAttitudeAngle, "Time (s)", "Lateral Attitude Angle (°)")
-
- # Euler Angles
- @cached_property
- def psi(self):
+ self.lateralAttitudeAngle = Function(
+ lateralAttitudeAngle,
+ "Time (s)",
+ "Lateral Attitude Angle (°)",
+ interpolation,
+ )
+ # Euler Angles
psi = (180 / np.pi) * (
np.arctan2(self.e3[:, 1], self.e0[:, 1])
+ np.arctan2(-self.e2[:, 1], -self.e1[:, 1])
) # Precession angle
psi = np.column_stack([self.e1[:, 0], psi]) # Precession angle
- return Function(psi, "Time (s)", "Precession Angle - ψ (°)")
+ self.psi = Function(psi, "Time (s)", "Precession Angle - ψ (°)", interpolation)
- @cached_property
- def phi(self):
phi = (180 / np.pi) * (
np.arctan2(self.e3[:, 1], self.e0[:, 1])
- np.arctan2(-self.e2[:, 1], -self.e1[:, 1])
) # Spin angle
phi = np.column_stack([self.e1[:, 0], phi]) # Spin angle
- return Function(phi, "Time (s)", "Spin Angle - φ (°)")
+ self.phi = Function(phi, "Time (s)", "Spin Angle - φ (°)", interpolation)
- @cached_property
- def theta(self):
theta = (
(180 / np.pi)
* 2
- * np.arcsin(-((self.e1[:, 1] ** 2 + self.e2[:, 1] ** 2) ** 0.5))
+ * np.arcsin(
+ np.clip(-((self.e1[:, 1] ** 2 + self.e2[:, 1] ** 2) ** 0.5), -1, 1)
+ )
) # Nutation angle
theta = np.column_stack([self.e1[:, 0], theta]) # Nutation angle
- return Function(theta, "Time (s)", "Nutation Angle - θ (°)")
-
- # Fluid Mechanics variables
- # Freestream Velocity
- @cached_property
- def streamVelocityX(self, interpolation="spline", extrapolation="natural"):
- streamVelocityX = np.column_stack(
- (self.vx[:, 0], self.windVelocityX[:, 1] - self.vx[:, 1])
- )
- streamVelocityX = Function(
- streamVelocityX,
- "Time (s)",
- "Freestream Velocity X (m/s)",
- interpolation,
- extrapolation,
- )
- return streamVelocityX
-
- @cached_property
- def streamVelocityY(self, interpolation="spline", extrapolation="natural"):
- streamVelocityY = np.column_stack(
- (self.vy[:, 0], self.windVelocityY[:, 1] - self.vy[:, 1])
- )
- streamVelocityY = Function(
- streamVelocityY,
- "Time (s)",
- "Freestream Velocity Y (m/s)",
- interpolation,
- extrapolation,
- )
- return streamVelocityY
-
- @cached_property
- def streamVelocityZ(self, interpolation="spline", extrapolation="natural"):
- streamVelocityZ = np.column_stack((self.vz[:, 0], -self.vz[:, 1]))
- streamVelocityZ = Function(
- streamVelocityZ,
- "Time (s)",
- "Freestream Velocity Z (m/s)",
- interpolation,
- extrapolation,
- )
- return streamVelocityZ
-
- @cached_property
- def freestreamSpeed(self):
- freestreamSpeed = (
- self.streamVelocityX**2
- + self.streamVelocityY**2
- + self.streamVelocityZ**2
- ) ** 0.5
- freestreamSpeed.setInputs("Time (s)")
- freestreamSpeed.setOutputs("Freestream Speed (m/s)")
- return freestreamSpeed
-
- # Apogee Freestream speed
- @cached_property
- def apogeeFreestreamSpeed(self):
- return self.freestreamSpeed(self.apogeeTime)
-
- # Mach Number
- @cached_property
- def MachNumber(self):
- MachNumber = self.freestreamSpeed / self.speedOfSound
- MachNumber.setInputs("Time (s)")
- MachNumber.setOutputs("Mach Number")
- return MachNumber
-
- @cached_property
- def maxMachNumberTime(self):
- maxMachNumberTimeIndex = np.argmax(self.MachNumber[:, 1])
- return self.MachNumber[maxMachNumberTimeIndex, 0]
-
- @cached_property
- def maxMachNumber(self):
- return self.MachNumber(self.maxMachNumberTime)
-
- # Reynolds Number
- @cached_property
- def ReynoldsNumber(self):
- ReynoldsNumber = (
- self.density * self.freestreamSpeed / self.dynamicViscosity
- ) * (2 * self.rocket.radius)
- ReynoldsNumber.setInputs("Time (s)")
- ReynoldsNumber.setOutputs("Reynolds Number")
- return ReynoldsNumber
-
- @cached_property
- def maxReynoldsNumberTime(self):
- maxReynoldsNumberTimeIndex = np.argmax(self.ReynoldsNumber[:, 1])
- return self.ReynoldsNumber[maxReynoldsNumberTimeIndex, 0]
-
- @cached_property
- def maxReynoldsNumber(self):
- return self.ReynoldsNumber(self.maxReynoldsNumberTime)
-
- # Dynamic Pressure
- @cached_property
- def dynamicPressure(self):
- dynamicPressure = 0.5 * self.density * self.freestreamSpeed**2
- dynamicPressure.setInputs("Time (s)")
- dynamicPressure.setOutputs("Dynamic Pressure (Pa)")
- return dynamicPressure
-
- @cached_property
- def maxDynamicPressureTime(self):
- maxDynamicPressureTimeIndex = np.argmax(self.dynamicPressure[:, 1])
- return self.dynamicPressure[maxDynamicPressureTimeIndex, 0]
-
- @cached_property
- def maxDynamicPressure(self):
- return self.dynamicPressure(self.maxDynamicPressureTime)
-
- # Total Pressure
- @cached_property
- def totalPressure(self):
- totalPressure = self.pressure * (1 + 0.2 * self.MachNumber**2) ** (3.5)
- totalPressure.setInputs("Time (s)")
- totalPressure.setOutputs("Total Pressure (Pa)")
- return totalPressure
-
- @cached_property
- def maxtotalPressureTime(self):
- maxtotalPressureTimeIndex = np.argmax(self.totalPressure[:, 1])
- return self.totalPressure[maxtotalPressureTimeIndex, 0]
-
- @cached_property
- def maxtotalPressure(self):
- return self.totalPressure(self.maxtotalPressureTime)
-
- # Dynamics functions and variables
-
- # Aerodynamic Lift and Drag
- @cached_property
- def aerodynamicLift(self):
- aerodynamicLift = (self.R1**2 + self.R2**2) ** 0.5
- aerodynamicLift.setInputs("Time (s)")
- aerodynamicLift.setOutputs("Aerodynamic Lift Force (N)")
- return aerodynamicLift
-
- @cached_property
- def aerodynamicDrag(self):
- aerodynamicDrag = -1 * self.R3
- aerodynamicDrag.setInputs("Time (s)")
- aerodynamicDrag.setOutputs("Aerodynamic Drag Force (N)")
- return aerodynamicDrag
-
- @cached_property
- def aerodynamicBendingMoment(self):
-
- aerodynamicBendingMoment = (self.M1**2 + self.M2**2) ** 0.5
- aerodynamicBendingMoment.setInputs("Time (s)")
- aerodynamicBendingMoment.setOutputs("Aerodynamic Bending Moment (N m)")
- return aerodynamicBendingMoment
-
- @cached_property
- def aerodynamicSpinMoment(self):
- aerodynamicSpinMoment = self.M3
- aerodynamicSpinMoment.setInputs("Time (s)")
- aerodynamicSpinMoment.setOutputs("Aerodynamic Spin Moment (N m)")
- return aerodynamicSpinMoment
-
- # Energy
- # Kinetic Energy
- @cached_property
- def rotationalEnergy(self):
+ self.theta = Function(
+ theta, "Time (s)", "Nutation Angle - θ (°)", interpolation
+ )
+
+ # Dynamics functions and variables
+ if self.rocket.railButtons is not None:
+ # Rail Button Forces
+ alpha = self.rocket.railButtons.angularPosition * (
+ np.pi / 180
+ ) # Rail buttons angular position
+ D1 = self.rocket.railButtons.distanceToCM[
+ 0
+ ] # Distance from Rail Button 1 (upper) to CM
+ D2 = self.rocket.railButtons.distanceToCM[
+ 1
+ ] # Distance from Rail Button 2 (lower) to CM
+ F11 = (self.R1 * D2 - self.M2) / (
+ D1 + D2
+ ) # Rail Button 1 force in the 1 direction
+ F12 = (self.R2 * D2 + self.M1) / (
+ D1 + D2
+ ) # Rail Button 1 force in the 2 direction
+ F21 = (self.R1 * D1 + self.M2) / (
+ D1 + D2
+ ) # Rail Button 2 force in the 1 direction
+ F22 = (self.R2 * D1 - self.M1) / (
+ D1 + D2
+ ) # Rail Button 2 force in the 2 direction
+ outOfRailTimeIndex = np.searchsorted(
+ F11[:, 0], self.outOfRailTime
+ ) # Find out of rail time index
+ # F11 = F11[:outOfRailTimeIndex + 1, :] # Limit force calculation to when rocket is in rail
+ # F12 = F12[:outOfRailTimeIndex + 1, :] # Limit force calculation to when rocket is in rail
+ # F21 = F21[:outOfRailTimeIndex + 1, :] # Limit force calculation to when rocket is in rail
+ # F22 = F22[:outOfRailTimeIndex + 1, :] # Limit force calculation to when rocket is in rail
+ self.railButton1NormalForce = F11 * np.cos(alpha) + F12 * np.sin(alpha)
+ self.railButton1NormalForce.setOutputs("Upper Rail Button Normal Force (N)")
+ self.railButton1ShearForce = F11 * -np.sin(alpha) + F12 * np.cos(alpha)
+ self.railButton1ShearForce.setOutputs("Upper Rail Button Shear Force (N)")
+ self.railButton2NormalForce = F21 * np.cos(alpha) + F22 * np.sin(alpha)
+ self.railButton2NormalForce.setOutputs("Lower Rail Button Normal Force (N)")
+ self.railButton2ShearForce = F21 * -np.sin(alpha) + F22 * np.cos(alpha)
+ self.railButton2ShearForce.setOutputs("Lower Rail Button Shear Force (N)")
+ # Rail Button Maximum Forces
+ if outOfRailTimeIndex == 0:
+ self.maxRailButton1NormalForce = 0
+ self.maxRailButton1ShearForce = 0
+ self.maxRailButton2NormalForce = 0
+ self.maxRailButton2ShearForce = 0
+ else:
+ self.maxRailButton1NormalForce = np.amax(
+ self.railButton1NormalForce[:outOfRailTimeIndex]
+ )
+ self.maxRailButton1ShearForce = np.amax(
+ self.railButton1ShearForce[:outOfRailTimeIndex]
+ )
+ self.maxRailButton2NormalForce = np.amax(
+ self.railButton2NormalForce[:outOfRailTimeIndex]
+ )
+ self.maxRailButton2ShearForce = np.amax(
+ self.railButton2ShearForce[:outOfRailTimeIndex]
+ )
+ # Aerodynamic Lift and Drag
+ self.aerodynamicLift = (self.R1**2 + self.R2**2) ** 0.5
+ self.aerodynamicLift.setOutputs("Aerodynamic Lift Force (N)")
+ self.aerodynamicDrag = -1 * self.R3
+ self.aerodynamicDrag.setOutputs("Aerodynamic Drag Force (N)")
+ self.aerodynamicBendingMoment = (self.M1**2 + self.M2**2) ** 0.5
+ self.aerodynamicBendingMoment.setOutputs("Aerodynamic Bending Moment (N m)")
+ self.aerodynamicSpinMoment = self.M3
+ self.aerodynamicSpinMoment.setOutputs("Aerodynamic Spin Moment (N m)")
+ # Energy
b = -self.rocket.distanceRocketPropellant
- mu = self.rocket.reducedMass
- Rz = self.rocket.inertiaZ
- Ri = self.rocket.inertiaI
- Tz = self.rocket.motor.inertiaZ
- Ti = self.rocket.motor.inertiaI
- I1, I2, I3 = (Ri + Ti + mu * b**2), (Ri + Ti + mu * b**2), (Rz + Tz)
- # Redefine I1, I2 and I3 grid
- grid = self.vx[:, 0]
- I1 = Function(np.column_stack([grid, I1(grid)]), "Time (s)")
- I2 = Function(np.column_stack([grid, I2(grid)]), "Time (s)")
- I3 = Function(np.column_stack([grid, I3(grid)]), "Time (s)")
- rotationalEnergy = 0.5 * (
- I1 * self.w1**2 + I2 * self.w2**2 + I3 * self.w3**2
- )
- rotationalEnergy.setInputs("Time (s)")
- rotationalEnergy.setOutputs("Rotational Kinetic Energy (J)")
- return rotationalEnergy
-
- @cached_property
- def translationalEnergy(self):
- # Redefine total mass grid
- totalMass = self.rocket.totalMass
- grid = self.vx[:, 0]
- totalMass = Function(np.column_stack([grid, totalMass(grid)]), "Time (s)")
- translationalEnergy = (
- 0.5 * totalMass * (self.vx**2 + self.vy**2 + self.vz**2)
- )
- translationalEnergy.setInputs("Time (s)")
- translationalEnergy.setOutputs("Translational Kinetic Energy (J)")
- return translationalEnergy
-
- @cached_property
- def kineticEnergy(self):
- kineticEnergy = self.rotationalEnergy + self.translationalEnergy
- kineticEnergy.setInputs("Time (s)")
- kineticEnergy.setOutputs("Kinetic Energy (J)")
- return kineticEnergy
-
- # Potential Energy
- @cached_property
- def potentialEnergy(self):
- # Redefine total mass grid
- totalMass = self.rocket.totalMass
- grid = self.vx[:, 0]
- totalMass = Function(np.column_stack([grid, totalMass(grid)]), "Time (s)")
- potentialEnergy = totalMass * self.env.g * self.z
- potentialEnergy.setInputs("Time (s)")
- return potentialEnergy
-
- # Total Mechanical Energy
- @cached_property
- def totalEnergy(self):
- totalEnergy = self.kineticEnergy + self.potentialEnergy
- totalEnergy.setInputs("Time (s)")
- totalEnergy.setOutputs("Total Mechanical Energy (J)")
- return totalEnergy
-
- # Thrust Power
- @cached_property
- def thrustPower(self):
- grid = self.vx[:, 0]
- # Redefine thrust grid
- thrust = Function(
- np.column_stack([grid, self.rocket.motor.thrust(grid)]), "Time (s)"
- )
- thrustPower = thrust * self.speed
- thrustPower.setOutputs("Thrust Power (W)")
- return thrustPower
-
- # Drag Power
- @cached_property
- def dragPower(self):
- dragPower = self.R3 * self.speed
- dragPower.setOutputs("Drag Power (W)")
- return dragPower
-
- # Angle of Attack
- @cached_property
- def angleOfAttack(self, interpolation="spline", extrapolation="natural"):
- """Angle of attack of the rocket with respect to the freestream
- velocity vector.
-
- Parameters
- ----------
- interpolation : str, optional
- Interpolation method, by default "spline"
- extrapolation : str, optional
- Extrapolation method, by default "natural"
-
- Returns
- -------
- angleOfAttack: Function
- Angle of attack of the rocket with respect to the freestream
- velocity vector.
- """
- dotProduct = [
- -self.attitudeVectorX.getValueOpt(i) * self.streamVelocityX.getValueOpt(i)
- - self.attitudeVectorY.getValueOpt(i) * self.streamVelocityY.getValueOpt(i)
- - self.attitudeVectorZ.getValueOpt(i) * self.streamVelocityZ.getValueOpt(i)
- for i in self.x[:, 0]
- ]
- # Define freestream speed list
- freestreamSpeed = [self.freestreamSpeed(i) for i in self.x[:, 0]]
- freestreamSpeed = np.nan_to_num(freestreamSpeed)
-
- # Normalize dot product
- dotProductNormalized = [
- i / j if j > 1e-6 else 0 for i, j in zip(dotProduct, freestreamSpeed)
- ]
- dotProductNormalized = np.nan_to_num(dotProductNormalized)
- dotProductNormalized = np.clip(dotProductNormalized, -1, 1)
-
- # Calculate angle of attack and convert to degrees
- angleOfAttack = np.rad2deg(np.arccos(dotProductNormalized))
- angleOfAttack = np.column_stack([self.x[:, 0], angleOfAttack])
-
- # Convert to Function
- angleOfAttack = Function(
- angleOfAttack,
- "Time (s)",
- "Angle Of Attack (°)",
- interpolation,
- extrapolation,
+ totalMass = self.rocket.totalMass
+ mu = self.rocket.reducedMass
+ Rz = self.rocket.inertiaZ
+ Ri = self.rocket.inertiaI
+ Tz = self.rocket.motor.inertiaZ
+ Ti = self.rocket.motor.inertiaI
+ I1, I2, I3 = (Ri + Ti + mu * b**2), (Ri + Ti + mu * b**2), (Rz + Tz)
+ # Redefine I1, I2 and I3 grid
+ grid = self.vx[:, 0]
+ I1 = Function(
+ source=np.column_stack([grid, I1(grid)]),
+ inputs="Time (s)",
+ outputs="Moment of Inertia (kg*m^2)",
+ interpolation=interpolation,
+ )
+ I2 = Function(
+ source=np.column_stack([grid, I2(grid)]),
+ inputs="Time (s)",
+ outputs="Moment of Inertia (kg*m^2)",
+ interpolation=interpolation,
+ )
+ I3 = Function(
+ source=np.column_stack([grid, I3(grid)]),
+ inputs="Time (s)",
+ outputs="Moment of Inertia (kg*m^2)",
+ interpolation=interpolation,
)
-
- return angleOfAttack
-
- # Stability and Control variables
- # Angular velocities frequency response - Fourier Analysis
- @cached_property
- def omega1FrequencyResponse(self):
+ # Redefine total mass grid
+ totalMass = Function(
+ source=np.column_stack([grid, totalMass(grid)]),
+ inputs="Time (s)",
+ outputs="Mass (kg)",
+ interpolation=interpolation,
+ )
+ # Redefine thrust grid
+ thrust = Function(
+ source=np.column_stack([grid, self.rocket.motor.thrust(grid)]),
+ inputs="Time (s)",
+ outputs="Thrust (N)",
+ interpolation=interpolation,
+ )
+ # Get some nicknames
+ vx, vy, vz = self.vx, self.vy, self.vz
+ w1, w2, w3 = self.w1, self.w2, self.w3
+ # Kinetic Energy
+ self.rotationalEnergy = 0.5 * (I1 * w1**2 + I2 * w2**2 + I3 * w3**2)
+ self.rotationalEnergy.setOutputs("Rotational Kinetic Energy (J)")
+ self.translationalEnergy = 0.5 * totalMass * (vx**2 + vy**2 + vz**2)
+ self.translationalEnergy.setOutputs("Translational Kinetic Energy (J)")
+ self.kineticEnergy = self.rotationalEnergy + self.translationalEnergy
+ self.kineticEnergy.setOutputs("Kinetic Energy (J)")
+ # Potential Energy
+ self.potentialEnergy = totalMass * self.env.g * self.z
+ self.potentialEnergy.setInputs("Time (s)")
+ # Total Mechanical Energy
+ self.totalEnergy = self.kineticEnergy + self.potentialEnergy
+ self.totalEnergy.setOutputs("Total Mechanical Energy (J)")
+ # Thrust Power
+ self.thrustPower = thrust * self.speed
+ self.thrustPower.setOutputs("Thrust Power (W)")
+ # Drag Power
+ self.dragPower = self.R3 * self.speed
+ self.dragPower.setOutputs("Drag Power (W)")
+
+ # Stability and Control variables
+ # Angular velocities frequency response - Fourier Analysis
+ # Omega 1 - w1
Fs = 100.0
# sampling rate
Ts = 1.0 / Fs
@@ -2673,13 +2032,13 @@ def omega1FrequencyResponse(self):
Y = np.fft.fft(y) / n # fft computing and normalization
Y = Y[range(n // 2)]
omega1FrequencyResponse = np.column_stack([frq, abs(Y)])
- omega1FrequencyResponse = Function(
- omega1FrequencyResponse, "Frequency (Hz)", "Omega 1 Angle Fourier Amplitude"
+ self.omega1FrequencyResponse = Function(
+ omega1FrequencyResponse,
+ "Frequency (Hz)",
+ "Omega 1 Angle Fourier Amplitude",
+ interpolation,
)
- return omega1FrequencyResponse
-
- @cached_property
- def omega2FrequencyResponse(self):
+ # Omega 2 - w2
Fs = 100.0
# sampling rate
Ts = 1.0 / Fs
@@ -2695,13 +2054,13 @@ def omega2FrequencyResponse(self):
Y = np.fft.fft(y) / n # fft computing and normalization
Y = Y[range(n // 2)]
omega2FrequencyResponse = np.column_stack([frq, abs(Y)])
- omega2FrequencyResponse = Function(
- omega2FrequencyResponse, "Frequency (Hz)", "Omega 2 Angle Fourier Amplitude"
+ self.omega2FrequencyResponse = Function(
+ omega2FrequencyResponse,
+ "Frequency (Hz)",
+ "Omega 2 Angle Fourier Amplitude",
+ interpolation,
)
- return omega2FrequencyResponse
-
- @cached_property
- def omega3FrequencyResponse(self):
+ # Omega 3 - w3
Fs = 100.0
# sampling rate
Ts = 1.0 / Fs
@@ -2717,13 +2076,13 @@ def omega3FrequencyResponse(self):
Y = np.fft.fft(y) / n # fft computing and normalization
Y = Y[range(n // 2)]
omega3FrequencyResponse = np.column_stack([frq, abs(Y)])
- omega3FrequencyResponse = Function(
- omega3FrequencyResponse, "Frequency (Hz)", "Omega 3 Angle Fourier Amplitude"
+ self.omega3FrequencyResponse = Function(
+ omega3FrequencyResponse,
+ "Frequency (Hz)",
+ "Omega 3 Angle Fourier Amplitude",
+ interpolation,
)
- return omega3FrequencyResponse
-
- @cached_property
- def attitudeFrequencyResponse(self):
+ # Attitude Frequency Response
Fs = 100.0
# sampling rate
Ts = 1.0 / Fs
@@ -2739,500 +2098,182 @@ def attitudeFrequencyResponse(self):
Y = np.fft.fft(y) / n # fft computing and normalization
Y = Y[range(n // 2)]
attitudeFrequencyResponse = np.column_stack([frq, abs(Y)])
- attitudeFrequencyResponse = Function(
+ self.attitudeFrequencyResponse = Function(
attitudeFrequencyResponse,
"Frequency (Hz)",
"Attitude Angle Fourier Amplitude",
+ interpolation,
)
- return attitudeFrequencyResponse
-
- # Static Margin
- @cached_property
- def staticMargin(self):
- return self.rocket.staticMargin
-
- # Rail Button Forces
- @cached_property
- def railButton1NormalForce(self):
- """Upper rail button normal force as a function of time
-
- Returns
- -------
- railButton1NormalForce: Function
- Upper rail button normal force as a function of time
- """
-
- if isinstance(self.calculate_rail_button_forces, tuple):
- F11, F12 = self.calculate_rail_button_forces[0:2]
- else:
- F11, F12 = self.calculate_rail_button_forces()[0:2]
- alpha = self.rocket.railButtons.angularPosition * (np.pi / 180)
- railButton1NormalForce = F11 * np.cos(alpha) + F12 * np.sin(alpha)
- railButton1NormalForce.setOutputs("Upper Rail Button Normal Force (N)")
-
- return railButton1NormalForce
-
- @cached_property
- def railButton1ShearForce(self):
- """Upper rail button shear force as a function of time
-
- Returns
- -------
- _type_
- _description_
- """
- if isinstance(self.calculate_rail_button_forces, tuple):
- F11, F12 = self.calculate_rail_button_forces[0:2]
- else:
- F11, F12 = self.calculate_rail_button_forces()[0:2]
- alpha = self.rocket.railButtons.angularPosition * (
- np.pi / 180
- ) # Rail buttons angular position
- railButton1ShearForce = F11 * -np.sin(alpha) + F12 * np.cos(alpha)
- railButton1ShearForce.setOutputs("Upper Rail Button Shear Force (N)")
-
- return railButton1ShearForce
-
- @cached_property
- def railButton2NormalForce(self):
- """Lower rail button normal force as a function of time
-
- Returns
- -------
- railButton2NormalForce: Function
- Lower rail button normal force as a function of time
- """
- if isinstance(self.calculate_rail_button_forces, tuple):
- F21, F22 = self.calculate_rail_button_forces[2:4]
- else:
- F21, F22 = self.calculate_rail_button_forces()[2:4]
- alpha = self.rocket.railButtons.angularPosition * (np.pi / 180)
- railButton2NormalForce = F21 * np.cos(alpha) + F22 * np.sin(alpha)
- railButton2NormalForce.setOutputs("Lower Rail Button Normal Force (N)")
-
- return railButton2NormalForce
-
- @cached_property
- def railButton2ShearForce(self):
- """Lower rail button shear force as a function of time
-
- Returns
- -------
- railButton2ShearForce: Function
- Lower rail button shear force as a function of time
- """
- if isinstance(self.calculate_rail_button_forces, tuple):
- F21, F22 = self.calculate_rail_button_forces[2:4]
- else:
- F21, F22 = self.calculate_rail_button_forces()[2:4]
- alpha = self.rocket.railButtons.angularPosition * (np.pi / 180)
- railButton2ShearForce = F21 * -np.sin(alpha) + F22 * np.cos(alpha)
- railButton2ShearForce.setOutputs("Lower Rail Button Shear Force (N)")
-
- return railButton2ShearForce
-
- @cached_property
- def maxRailButton1NormalForce(self):
- """Maximum upper rail button normal force, in Newtons
-
- Returns
- -------
- maxRailButton1NormalForce: float
- Maximum upper rail button normal force, in Newtons
- """
- if isinstance(self.calculate_rail_button_forces, tuple):
- F11 = self.calculate_rail_button_forces[0]
- else:
- F11 = self.calculate_rail_button_forces()[0]
- outOfRailTimeIndex = np.searchsorted(F11[:, 0], self.outOfRailTime)
- if outOfRailTimeIndex == 0:
- return 0
- else:
- return np.max(self.railButton1NormalForce[:outOfRailTimeIndex])
-
- @cached_property
- def maxRailButton1ShearForce(self):
- """Maximum upper rail button shear force, in Newtons
-
- Returns
- -------
- maxRailButton1ShearForce: float
- Maximum upper rail button shear force, in Newtons
- """
- if isinstance(self.calculate_rail_button_forces, tuple):
- F11 = self.calculate_rail_button_forces[0]
- else:
- F11 = self.calculate_rail_button_forces()[0]
- outOfRailTimeIndex = np.searchsorted(F11[:, 0], self.outOfRailTime)
- if outOfRailTimeIndex == 0:
- return 0
- else:
- return np.max(self.railButton1ShearForce[:outOfRailTimeIndex])
-
- @cached_property
- def maxRailButton2NormalForce(self):
- """Maximum lower rail button normal force, in Newtons
-
- Returns
- -------
- maxRailButton2NormalForce: float
- Maximum lower rail button normal force, in Newtons
- """
- if isinstance(self.calculate_rail_button_forces, tuple):
- F11 = self.calculate_rail_button_forces[0]
- else:
- F11 = self.calculate_rail_button_forces()[0]
- outOfRailTimeIndex = np.searchsorted(F11[:, 0], self.outOfRailTime)
- if outOfRailTimeIndex == 0:
- return 0
- else:
- return np.max(self.railButton2NormalForce[:outOfRailTimeIndex])
-
- @cached_property
- def maxRailButton2ShearForce(self):
- """Maximum lower rail button shear force, in Newtons
-
- Returns
- -------
- maxRailButton2ShearForce: float
- Maximum lower rail button shear force, in Newtons
- """
- if isinstance(self.calculate_rail_button_forces, tuple):
- F11 = self.calculate_rail_button_forces[0]
- else:
- F11 = self.calculate_rail_button_forces()[0]
- outOfRailTimeIndex = np.searchsorted(F11[:, 0], self.outOfRailTime)
- if outOfRailTimeIndex == 0:
- return 0
- else:
- return np.max(self.railButton2ShearForce[:outOfRailTimeIndex])
-
- @cached_property
- def drift(self):
- """Rocket horizontal distance to tha launch point, in meters, at each
- time step.
-
- Returns
- -------
- drift: Function
- Rocket horizontal distance to tha launch point, in meters, at each
- time step.
- """
- drift = np.column_stack(
- (self.x[:, 0], (self.x[:, 1] ** 2 + self.y[:, 1] ** 2) ** 0.5)
- )
- drift = Function(
- drift,
- "Time (s)",
- "Horizontal Distance to Launch Point (m)",
+ # Static Margin
+ self.staticMargin = self.rocket.staticMargin
+
+ # Fluid Mechanics variables
+ # Freestream Velocity
+ self.streamVelocityX = self.windVelocityX - self.vx
+ self.streamVelocityX.setOutputs("Freestream Velocity X (m/s)")
+ self.streamVelocityY = self.windVelocityY - self.vy
+ self.streamVelocityY.setOutputs("Freestream Velocity Y (m/s)")
+ self.streamVelocityZ = -1 * self.vz
+ self.streamVelocityZ.setOutputs("Freestream Velocity Z (m/s)")
+ self.freestreamSpeed = (
+ self.streamVelocityX**2
+ + self.streamVelocityY**2
+ + self.streamVelocityZ**2
+ ) ** 0.5
+ self.freestreamSpeed.setOutputs("Freestream Speed (m/s)")
+ # Apogee Freestream speed
+ self.apogeeFreestreamSpeed = self.freestreamSpeed(self.apogeeTime)
+ # Mach Number
+ self.MachNumber = self.freestreamSpeed / self.speedOfSound
+ self.MachNumber.setOutputs("Mach Number")
+ maxMachNumberTimeIndex = np.argmax(self.MachNumber[:, 1])
+ self.maxMachNumberTime = self.MachNumber[maxMachNumberTimeIndex, 0]
+ self.maxMachNumber = self.MachNumber[maxMachNumberTimeIndex, 1]
+ # Reynolds Number
+ self.ReynoldsNumber = (
+ self.density * self.freestreamSpeed / self.dynamicViscosity
+ ) * (2 * self.rocket.radius)
+ self.ReynoldsNumber.setOutputs("Reynolds Number")
+ maxReynoldsNumberTimeIndex = np.argmax(self.ReynoldsNumber[:, 1])
+ self.maxReynoldsNumberTime = self.ReynoldsNumber[maxReynoldsNumberTimeIndex, 0]
+ self.maxReynoldsNumber = self.ReynoldsNumber[maxReynoldsNumberTimeIndex, 1]
+ # Dynamic Pressure
+ self.dynamicPressure = 0.5 * self.density * self.freestreamSpeed**2
+ self.dynamicPressure.setOutputs("Dynamic Pressure (Pa)")
+ maxDynamicPressureTimeIndex = np.argmax(self.dynamicPressure[:, 1])
+ self.maxDynamicPressureTime = self.dynamicPressure[
+ maxDynamicPressureTimeIndex, 0
+ ]
+ self.maxDynamicPressure = self.dynamicPressure[maxDynamicPressureTimeIndex, 1]
+ # Total Pressure
+ self.totalPressure = self.pressure * (1 + 0.2 * self.MachNumber**2) ** (3.5)
+ self.totalPressure.setOutputs("Total Pressure (Pa)")
+ maxtotalPressureTimeIndex = np.argmax(self.totalPressure[:, 1])
+ self.maxtotalPressureTime = self.totalPressure[maxtotalPressureTimeIndex, 0]
+ self.maxtotalPressure = self.totalPressure[maxDynamicPressureTimeIndex, 1]
+ # Angle of Attack
+ angleOfAttack = []
+ for i in range(len(self.attitudeVectorX[:, 1])):
+ dotProduct = -(
+ self.attitudeVectorX[i, 1] * self.streamVelocityX[i, 1]
+ + self.attitudeVectorY[i, 1] * self.streamVelocityY[i, 1]
+ + self.attitudeVectorZ[i, 1] * self.streamVelocityZ[i, 1]
+ )
+ if self.freestreamSpeed[i, 1] < 1e-6:
+ angleOfAttack.append([self.freestreamSpeed[i, 0], 0])
+ else:
+ dotProductNormalized = dotProduct / self.freestreamSpeed[i, 1]
+ dotProductNormalized = (
+ 1 if dotProductNormalized > 1 else dotProductNormalized
+ )
+ dotProductNormalized = (
+ -1 if dotProductNormalized < -1 else dotProductNormalized
+ )
+ angleOfAttack.append(
+ [
+ self.freestreamSpeed[i, 0],
+ (180 / np.pi) * np.arccos(dotProductNormalized),
+ ]
+ )
+ self.angleOfAttack = Function(
+ angleOfAttack, "Time (s)", "Angle Of Attack (°)", "linear"
)
- return drift
-
- @cached_property
- def bearing(self, interpolation="spline", extrapolation="natural"):
- """Rocket bearing compass, in degrees, at each time step.
-
- Returns
- -------
- bearing: Function
- Rocket bearing, in degrees, at each time step.
- """
-
- # Get some nicknames
- t = self.x[:, 0]
- x, y = self.x[:, 1], self.y[:, 1]
+ # Converts x and y positions to lat and lon
+ # We are currently considering the earth as a sphere.
bearing = []
- for i in range(len(t)):
- # Forcing arctan2(0, 0) = self.heading
- if abs(x[i]) < 1e-6 and abs(y[i]) < 1e-6:
- bearing.append(np.deg2rad(self.heading))
- elif abs(x[i]) < 1e-6: # check if the rocket is on x axis
- if y[i] > 0:
- bearing.append(0)
+ distance = [
+ (self.x[i][1] ** 2 + self.y[i][1] ** 2) ** 0.5 for i in range(len(self.x))
+ ]
+ for i in range(len(self.x)):
+ # Check if the point is over the grid (i.e. if x*y == 0)
+ if self.x[i][1] == 0:
+ if self.y[i][1] < 0:
+ bearing.append(3.14159265359)
else:
- bearing.append(np.pi)
- elif abs(y[i]) < 1e-6: # check if the rocket is on x axis
- if x[i] > 0:
- bearing.append(np.pi / 2)
+ bearing.append(0)
+ continue
+ if self.y[i][1] == 0:
+ if self.x[i][1] < 0:
+ bearing.append(3 * 3.14159265359 / 2)
+ elif self.x[i][1] > 0:
+ bearing.append(3.14159265359 / 2)
else:
- bearing.append(3 * np.pi / 2)
- else:
- # Calculate bearing as the azimuth considering different quadrants
- if x[i] * y[i] < 0 and x[i] < 0: # Fourth quadrant
- bearing.append(-np.pi / 2 + np.arctan(abs(y[i]) / abs(x[i])))
- elif x[i] * y[i] < 0 and x[i] > 0: # Second quadrant
- bearing.append(np.pi / 2 + np.arctan(abs(x[i]) / abs(y[i])))
- elif x[i] * y[i] > 0 and x[i] < 0: # Third quadrant
- bearing.append(np.pi + np.arctan(abs(x[i]) / abs(y[i])))
- else: # First quadrant
- bearing.append(np.arctan(abs(x[i]) / abs(y[i])))
-
- bearing = np.rad2deg(bearing)
- bearing = np.column_stack((t, bearing))
- print(bearing)
- bearing = Function(
- bearing,
- "Time (s)",
- "Bearing (deg)",
- interpolation,
- extrapolation,
- )
-
- return bearing
-
- @cached_property
- def latitude(self):
- """Rocket latitude coordinate, in degrees, at each time step.
-
- Returns
- -------
- latitude: Function
- Rocket latitude coordinate, in degrees, at each time step.
- """
- lat1 = np.deg2rad(self.env.lat) # Launch lat point converted to radians
-
- # Applies the haversine equation to find final lat/lon coordinates
- latitude = np.rad2deg(
- np.arcsin(
- np.sin(lat1) * np.cos(self.drift[:, 1] / self.env.earthRadius)
- + np.cos(lat1)
- * np.sin(self.drift[:, 1] / self.env.earthRadius)
- * np.cos(np.deg2rad(self.bearing[:, 1]))
- )
- )
- latitude = np.column_stack((self.x[:, 0], latitude))
- latitude = Function(latitude, "Time (s)", "Latitude (°)", "linear")
-
- return latitude
+ continue
+ continue
+
+ # Calculate bearing as the azimuth considering different quadrants
+ if self.x[i][1] * self.y[i][1] > 0 and self.x[i][1] > 0:
+ bearing.append(math.atan(abs(self.x[i][1]) / abs(self.y[i][1])))
+ elif self.x[i][1] * self.y[i][1] < 0 and self.x[i][1] > 0:
+ bearing.append(
+ 3.14159265359 / 2 + math.atan(abs(self.y[i][1]) / abs(self.x[i][1]))
+ )
+ elif self.x[i][1] * self.y[i][1] > 0 and self.x[i][1] < 0:
+ bearing.append(
+ 3.14159265359 + math.atan(abs(self.x[i][1]) / abs(self.y[i][1]))
+ )
+ elif self.x[i][1] * self.y[i][1] < 0 and self.x[i][1] < 0:
+ bearing.append(
+ 3 * 3.14159265359 / 2
+ + math.atan(abs(self.y[i][1]) / abs(self.x[i][1]))
+ )
- @cached_property
- def longitude(self):
- """Rocket longitude coordinate, in degrees, at each time step.
+ # Store values of distance and bearing using appropriate units
+ # self.distance = distance # Must be in meters
+ # self.bearing = bearing # Must be in radians
- Returns
- -------
- longitude: Function
- Rocket longitude coordinate, in degrees, at each time step.
- """
- lat1 = np.deg2rad(self.env.lat) # Launch lat point converted to radians
- lon1 = np.deg2rad(self.env.lon) # Launch lon point converted to radians
+ lat1 = (
+ 3.14159265359 * self.env.lat / 180
+ ) # Launch lat point converted to radians
+ lon1 = (
+ 3.14159265359 * self.env.lon / 180
+ ) # Launch long point converted to radians
+ R = self.env.earthRadius
+ lat2 = []
+ lon2 = []
# Applies the haversine equation to find final lat/lon coordinates
- longitude = np.rad2deg(
- lon1
- + np.arctan2(
- np.sin(np.deg2rad(self.bearing[:, 1]))
- * np.sin(self.drift[:, 1] / self.env.earthRadius)
- * np.cos(lat1),
- np.cos(self.drift[:, 1] / self.env.earthRadius)
- - np.sin(lat1) * np.sin(np.deg2rad(self.latitude[:, 1])),
- )
- )
-
- longitude = np.column_stack((self.x[:, 0], longitude))
- longitude = Function(longitude, "Time (s)", "Longitude (°)", "linear")
-
- return 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 = [], [], []
- alpha1, alpha2, alpha3 = [], [], []
- # Go through each time step and calculate accelerations
- # Get flight phases
- for phase_index, phase in self.timeIterator(self.flightPhases):
- initTime = phase.t
- finalTime = self.flightPhases[phase_index + 1].t
- currentDerivative = 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 initTime < step[0] <= finalTime:
- # Get derivatives
- uDot = currentDerivative(step[0], step[1:])
- # Get accelerations
- ax_value, ay_value, az_value = uDot[3:6]
- alpha1_value, alpha2_value, alpha3_value = uDot[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
- - dynamicViscosity
- - speedOfSound
-
- Parameters
- ----------
- None
-
- 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 ? direction at each time step
- self.M3_list: list
- Aerodynamic bending moment in ? 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.dynamicViscosity_list: list
- Dynamic viscosity at each time step
- elf_list._speedOfSound: list
- Speed of sound at each time step
- self.windVelocityX_list: list
- Wind velocity in x direction at each time step
- self.windVelocityY_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.dynamicViscosity_list = []
- self.speedOfSound_list = []
- self.windVelocityX_list = []
- self.windVelocityY_list = []
-
- # Go through each time step and calculate forces and atmospheric values
- # Get flight phases
- for phase_index, phase in self.timeIterator(self.flightPhases):
- initTime = phase.t
- finalTime = self.flightPhases[phase_index + 1].t
- currentDerivative = 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 initTime < step[0] <= finalTime or (initTime == 0 and step[0] == 0):
- # Call derivatives in post processing mode
- uDot = currentDerivative(step[0], step[1:], postProcessing=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.dynamicViscosity_list,
- self.speedOfSound_list,
- self.windVelocityX_list,
- self.windVelocityY_list,
- ]
-
- return temporary_values
+ for i in range(len(self.x)):
+ # Please notice that for distances lower than 1 centimeter the difference on latitude or longitude too small
+ if abs(self.y[i][1]) < 1e-2:
+ lat2.append(self.env.lat)
+ else:
+ lat2.append(
+ (180 / 3.14159265359)
+ * math.asin(
+ math.sin(lat1) * math.cos(distance[i] / R)
+ + math.cos(lat1)
+ * math.sin(distance[i] / R)
+ * math.cos(bearing[i])
+ )
+ )
+ if abs(self.x[i][1]) < 1e-2:
+ lon2.append(self.env.lon)
+ else:
+ lon2.append(
+ (180 / 3.14159265359)
+ * (
+ lon1
+ + math.atan2(
+ math.sin(bearing[i])
+ * math.sin(distance[i] / R)
+ * math.cos(lat1),
+ math.cos(distance[i] / R)
+ - math.sin(lat1) * math.sin(lat2[i]),
+ )
+ )
+ )
- @cached_property
- def calculate_rail_button_forces(self):
- """Calculate the forces applied to the rail buttons while rocket is still
- on the launch rail. It will return 0 if none rail buttons are defined.
+ latitude = [[self.solution[i][0], lat2[i]] for i in range(len(self.solution))]
+ longitude = [[self.solution[i][0], lon2[i]] for i in range(len(self.solution))]
- Returns
- -------
- F11: Function
- Rail Button 1 force in the 1 direction
- F12: Function
- Rail Button 1 force in the 2 direction
- F21: Function
- Rail Button 2 force in the 1 direction
- F22: Function
- Rail Button 2 force in the 2 direction
- """
+ # Store final values of lat/lon as a function of time
+ self.latitude = Function(latitude, "Time (s)", "Latitude (°)", "linear")
+ self.longitude = Function(longitude, "Time (s)", "Longitude (°)", "linear")
- if self.rocket.railButtons is None:
- warnings.warn(
- "Trying to calculate rail button forces without rail buttons defined."
- )
- return 0, 0, 0, 0
-
- # Distance from Rail Button 1 (upper) to CM
- D1 = self.rocket.railButtons.distanceToCM[0]
- # Distance from Rail Button 2 (lower) to CM
- D2 = self.rocket.railButtons.distanceToCM[1]
- F11 = (self.R1 * D2 - self.M2) / (D1 + D2)
- F11.setOutputs("Upper button force direction 1 (m)")
- F12 = (self.R2 * D2 + self.M1) / (D1 + D2)
- F12.setOutputs("Upper button force direction 2 (m)")
- F21 = (self.R1 * D1 + self.M2) / (D1 + D2)
- F21.setOutputs("Lower button force direction 1 (m)")
- F22 = (self.R2 * D1 - self.M1) / (D1 + D2)
- F22.setOutputs("Lower button force direction 2 (m)")
-
- # F11 = F11[:outOfRailTimeIndex + 1, :] # Limit force calculation to when rocket is in rail
- # F12 = F12[:outOfRailTimeIndex + 1, :] # Limit force calculation to when rocket is in rail
- # F21 = F21[:outOfRailTimeIndex + 1, :] # Limit force calculation to when rocket is in rail
- # F22 = F22[:outOfRailTimeIndex + 1, :] # Limit force calculation to when rocket is in rail
-
- return F11, F12, F21, F22
-
- @cached_property
- def __calculate_pressure_signal(self):
- """Calculate the pressure signal from the pressure sensor.
- It creates a SignalFunction attribute in the parachute object.
- Parachute works as a subclass of Rocket class.
+ # Post process other quantities
- Returns
- -------
- None
- """
# Transform parachute sensor feed into functions
for parachute in self.rocket.parachutes:
parachute.cleanPressureSignalFunction = Function(
@@ -3251,23 +2292,6 @@ def __calculate_pressure_signal(self):
parachute.noiseSignal, "Time (s)", "Pressure Noise (Pa)", "linear"
)
- return None
-
- def postProcess(self, interpolation="spline", extrapolation="natural"):
- """Post-process all Flight information produced during
- simulation. Includes the calculation of maximum values,
- calculation of secondary values such as energy and conversion
- of lists to Function objects to facilitate plotting.
-
- Parameters
- ----------
- None
-
- Return
- ------
- None
- """
-
# Register post processing
self.postProcessed = True
@@ -3284,6 +2308,9 @@ def info(self):
------
None
"""
+ # Post-process results
+ if self.postProcessed is False:
+ self.postProcess()
# Get index of out of rail time
outOfRailTimeIndexes = np.nonzero(self.x[:, 0] == self.outOfRailTime)
@@ -3477,6 +2504,9 @@ def printInitialConditionsData(self):
------
None
"""
+ # Post-process results
+ if self.postProcessed is False:
+ self.postProcess()
print(
"Position - x: {:.2f} m | y: {:.2f} m | z: {:.2f} m".format(
@@ -3590,6 +2620,9 @@ def plot3dTrajectory(self):
------
None
"""
+ # Post-process results
+ if self.postProcessed is False:
+ self.postProcess()
# Get max and min x and y
maxZ = max(self.z[:, 1] - self.env.elevation)
@@ -3645,6 +2678,9 @@ def plotLinearKinematicsData(self):
------
None
"""
+ # Post-process results
+ if self.postProcessed is False:
+ self.postProcess()
# Velocity and acceleration plots
fig2 = plt.figure(figsize=(9, 12))
@@ -3720,6 +2756,9 @@ def plotAttitudeData(self):
------
None
"""
+ # Post-process results
+ if self.postProcessed is False:
+ self.postProcess()
# Get index of time before parachute event
if len(self.parachuteEvents) > 0:
@@ -3785,6 +2824,9 @@ def plotFlightPathAngleData(self):
------
None
"""
+ # Post-process results
+ if self.postProcessed is False:
+ self.postProcess()
# Get index of time before parachute event
if len(self.parachuteEvents) > 0:
@@ -3837,6 +2879,9 @@ def plotAngularKinematicsData(self):
------
None
"""
+ # Post-process results
+ if self.postProcessed is False:
+ self.postProcess()
# Get index of time before parachute event
if len(self.parachuteEvents) > 0:
@@ -3918,6 +2963,9 @@ def plotTrajectoryForceData(self):
------
None
"""
+ # Post-process results
+ if self.postProcessed is False:
+ self.postProcess()
# Get index of out of rail time
outOfRailTimeIndexes = np.nonzero(self.x[:, 0] == self.outOfRailTime)
@@ -3934,46 +2982,51 @@ def plotTrajectoryForceData(self):
eventTimeIndex = -1
# Rail Button Forces
- fig6 = plt.figure(figsize=(9, 6))
-
- ax1 = plt.subplot(211)
- ax1.plot(
- self.railButton1NormalForce[:outOfRailTimeIndex, 0],
- self.railButton1NormalForce[:outOfRailTimeIndex, 1],
- label="Upper Rail Button",
- )
- ax1.plot(
- self.railButton2NormalForce[:outOfRailTimeIndex, 0],
- self.railButton2NormalForce[:outOfRailTimeIndex, 1],
- label="Lower Rail Button",
- )
- ax1.set_xlim(0, self.outOfRailTime if self.outOfRailTime > 0 else self.tFinal)
- ax1.legend()
- ax1.grid(True)
- ax1.set_xlabel("Time (s)")
- ax1.set_ylabel("Normal Force (N)")
- ax1.set_title("Rail Buttons Normal Force")
-
- ax2 = plt.subplot(212)
- ax2.plot(
- self.railButton1ShearForce[:outOfRailTimeIndex, 0],
- self.railButton1ShearForce[:outOfRailTimeIndex, 1],
- label="Upper Rail Button",
- )
- ax2.plot(
- self.railButton2ShearForce[:outOfRailTimeIndex, 0],
- self.railButton2ShearForce[:outOfRailTimeIndex, 1],
- label="Lower Rail Button",
- )
- ax2.set_xlim(0, self.outOfRailTime if self.outOfRailTime > 0 else self.tFinal)
- ax2.legend()
- ax2.grid(True)
- ax2.set_xlabel("Time (s)")
- ax2.set_ylabel("Shear Force (N)")
- ax2.set_title("Rail Buttons Shear Force")
+ if self.rocket.railButtons is not None:
+ fig6 = plt.figure(figsize=(9, 6))
+
+ ax1 = plt.subplot(211)
+ ax1.plot(
+ self.railButton1NormalForce[:outOfRailTimeIndex, 0],
+ self.railButton1NormalForce[:outOfRailTimeIndex, 1],
+ label="Upper Rail Button",
+ )
+ ax1.plot(
+ self.railButton2NormalForce[:outOfRailTimeIndex, 0],
+ self.railButton2NormalForce[:outOfRailTimeIndex, 1],
+ label="Lower Rail Button",
+ )
+ ax1.set_xlim(
+ 0, self.outOfRailTime if self.outOfRailTime > 0 else self.tFinal
+ )
+ ax1.legend()
+ ax1.grid(True)
+ ax1.set_xlabel("Time (s)")
+ ax1.set_ylabel("Normal Force (N)")
+ ax1.set_title("Rail Buttons Normal Force")
+
+ ax2 = plt.subplot(212)
+ ax2.plot(
+ self.railButton1ShearForce[:outOfRailTimeIndex, 0],
+ self.railButton1ShearForce[:outOfRailTimeIndex, 1],
+ label="Upper Rail Button",
+ )
+ ax2.plot(
+ self.railButton2ShearForce[:outOfRailTimeIndex, 0],
+ self.railButton2ShearForce[:outOfRailTimeIndex, 1],
+ label="Lower Rail Button",
+ )
+ ax2.set_xlim(
+ 0, self.outOfRailTime if self.outOfRailTime > 0 else self.tFinal
+ )
+ ax2.legend()
+ ax2.grid(True)
+ ax2.set_xlabel("Time (s)")
+ ax2.set_ylabel("Shear Force (N)")
+ ax2.set_title("Rail Buttons Shear Force")
- plt.subplots_adjust(hspace=0.5)
- plt.show()
+ plt.subplots_adjust(hspace=0.5)
+ plt.show()
# Aerodynamic force and moment plots
fig7 = plt.figure(figsize=(9, 12))
@@ -4042,6 +3095,9 @@ def plotEnergyData(self):
-------
None
"""
+ # Post-process results
+ if self.postProcessed is False:
+ self.postProcess()
# Get index of out of rail time
outOfRailTimeIndexes = np.nonzero(self.x[:, 0] == self.outOfRailTime)
@@ -4137,6 +3193,9 @@ def plotFluidMechanicsData(self):
------
None
"""
+ # Post-process results
+ if self.postProcessed is False:
+ self.postProcess()
# Get index of out of rail time
outOfRailTimeIndexes = np.nonzero(self.x[:, 0] == self.outOfRailTime)
@@ -4219,6 +3278,9 @@ def calculateFinFlutterAnalysis(self, finThickness, shearModulus):
------
None
"""
+ # Post-process results
+ if self.postProcessed is False:
+ self.postProcess()
s = (self.rocket.tipChord + self.rocket.rootChord) * self.rocket.span / 2
ar = self.rocket.span * self.rocket.span / s
@@ -4364,6 +3426,9 @@ def plotStabilityAndControlData(self):
------
None
"""
+ # Post-process results
+ if self.postProcessed is False:
+ self.postProcess()
fig9 = plt.figure(figsize=(9, 6))
@@ -4435,6 +3500,9 @@ def plotPressureSignals(self):
------
None
"""
+ # Post-process results
+ if self.postProcessed is False:
+ self.postProcess()
if len(self.rocket.parachutes) == 0:
plt.figure()
@@ -4482,6 +3550,8 @@ def exportPressures(self, fileName, timeStep):
------
None
"""
+ if self.postProcessed is False:
+ self.postProcess()
timePoints = np.arange(0, self.tFinal, timeStep)
@@ -4528,6 +3598,8 @@ def exportData(self, fileName, *variables, timeStep=None):
will be exported. Otherwise, linear interpolation is carried out to
calculate values at the desired time steps. Example: 0.001.
"""
+ if self.postProcessed is False:
+ self.postProcess()
# Fast evaluation for the most basic scenario
if timeStep is None and len(variables) == 0:
@@ -4580,19 +3652,7 @@ def exportData(self, fileName, *variables, timeStep=None):
# Loop through variables, get points and names (for the header)
for variable in variables:
- if variable in self.__dict__.keys():
- variableFunction = self.__dict__[variable]
- # Deal with variables that are not Flight attributes (e.g. 'angleOfAttack')
- # Usually, these are properties of the Flight class
- else:
- try:
- obj = getattr(self.__class__, variable)
- if isinstance(obj, property) or isinstance(obj, cached_property):
- variableFunction = obj.__get__(self, self.__class__)
- except AttributeError:
- raise AttributeError(
- "Variable '{}' not found in Flight class".format(variable)
- )
+ variableFunction = self.__dict__[variable]
variablePoints = variableFunction(timePoints)
exportedMatrix += [variablePoints]
exportedHeader += ", " + variableFunction.__outputs__[0]
@@ -4648,6 +3708,8 @@ def exportKML(
None
"""
# Define time points vector
+ if self.postProcessed is False:
+ self.postProcess()
if timeStep is None:
# Get time from any Function, should all be the same
timePoints = self.z[:, 0]
@@ -4701,6 +3763,9 @@ def allInfo(self):
------
None
"""
+ # Post-process results
+ if self.postProcessed is False:
+ self.postProcess()
# Print initial conditions
print("Initial Conditions\n")
diff --git a/rocketpy/utilities.py b/rocketpy/utilities.py
index 02ee2ca51..7671b414b 100644
--- a/rocketpy/utilities.py
+++ b/rocketpy/utilities.py
@@ -3,12 +3,15 @@
__copyright__ = "Copyright 20XX, RocketPy Team"
__license__ = "MIT"
+import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
+import matplotlib.pyplot as plt
from .Environment import Environment
from .Function import Function
+# Parachutes related functions
# TODO: Needs tests
def compute_CdS_from_drop_test(
@@ -198,6 +201,155 @@ def du(z, u):
return altitudeFunction, velocityFunction, final_sol
+def compareTrajectories(
+ trajectory_list,
+ names=None,
+ legend=True,
+):
+ """Creates a trajectory plot combining the trajectories listed.
+ This function was created based two source-codes:
+ - Mateus Stano: https://github.com/RocketPy-Team/Hackathon_2020/pull/123
+ - Dyllon Preston: https://github.com/Dyllon-P/MBS-Template/blob/main/MBS.py
+
+ Parameters
+ ----------
+ trajectory_list : list, array
+ List of trajectories. Must be in the form of [trajectory_1, trajectory_2, ..., trajectory_n]
+ where each element is a list with the arrays regarding positions in x, y and z [x, y, z].
+ The trajectories must be in the same reference frame. The z coordinate must be referenced
+ to the ground or to the sea level, but it is important that all trajectories are passed
+ in the same reference.
+ names : list, optional
+ List of strings with the name of each trajectory inputted. The names must be in
+ the same order as the trajectories in trajectory_list. If no names are passed, the
+ trajectories will be named as "Trajectory 1", "Trajectory 2", ..., "Trajectory n".
+ legend : boolean, optional
+ Whether legend will or will not be plotted. Default is True
+
+ Returns
+ -------
+ None
+
+ """
+ # TODO: Allow the user to catch different planes (xy, yz, xz) from the main plot (this can be done in a separate function)
+ # TODO: Allow the user to set the colors or color style
+ # TODO: Allow the user to set the line style
+
+ # Initialize variables
+ maxX, maxY, maxZ, minX, minY, minZ, maxXY, minXY = 0, 0, 0, 0, 0, 0, 0, 0
+
+ names = (
+ [("Trajectory " + str(i + 1)) for i in range(len(trajectory_list))]
+ if names == None
+ else names
+ )
+
+ # Create the figure
+ fig1 = plt.figure(figsize=(9, 9))
+ ax1 = plt.subplot(111, projection="3d")
+
+ # Iterate through trajectories
+ for i, trajectory in enumerate(trajectory_list):
+
+ x, y, z = trajectory
+
+ # Find max/min values for each component
+ maxX = max(x) if max(x) > maxX else maxX
+ maxY = max(y) if max(y) > maxY else maxY
+ maxZ = max(z) if max(z) > maxZ else maxZ
+ minX = min(x) if min(x) < minX else minX
+ minY = min(x) if min(x) < minX else minX
+ minZ = min(z) if min(z) < minZ else minZ
+ maxXY = max(maxX, maxY) if max(maxX, maxY) > maxXY else maxXY
+ minXY = min(minX, minY) if min(minX, minY) > minXY else minXY
+
+ # Add Trajectory as a plot in main figure
+ ax1.plot(x, y, z, linewidth="2", label=names[i])
+
+ # Plot settings
+ ax1.scatter(0, 0, 0)
+ ax1.set_xlabel("X - East (m)")
+ ax1.set_ylabel("Y - North (m)")
+ ax1.set_zlabel("Z - Altitude (m)")
+ ax1.set_title("Flight Trajectories Comparison")
+ ax1.set_zlim3d([minZ, maxZ])
+ ax1.set_ylim3d([minXY, maxXY])
+ ax1.set_xlim3d([minXY, maxXY])
+ ax1.view_init(15, 45)
+ if legend:
+ plt.legend()
+ plt.show()
+
+ return None
+
+
+def compareFlightTrajectories(
+ flight_list,
+ names=None,
+ legend=True,
+):
+ """Creates a trajectory plot that is the combination of the trajectories of
+ the Flight objects passed via a Python list.
+
+ Parameters
+ ----------
+ flight_list : list, array
+ List of FLight objects. The flights must be in the same reference frame.
+ names : list, optional
+ List of strings with the name of each trajectory inputted. The names must be in
+ the same order as the trajectories in flight_list
+ legend : boolean, optional
+ Whether legend will or will not be included. Default is True
+
+ Returns
+ -------
+ None
+
+ """
+ # TODO: Allow the user to catch different planes (xy, yz, xz) from the main plot
+ # TODO: Allow the user to set the colors or color style
+ # TODO: Allow the user to set the line style
+
+ # Iterate through Flight objects and create a list of trajectories
+ trajectory_list = []
+ for flight in flight_list:
+
+ # Check post process
+ if flight.postProcessed is False:
+ flight.postProcess()
+
+ # Get trajectories
+ x = flight.x[:, 1]
+ y = flight.y[:, 1]
+ z = flight.z[:, 1] - flight.env.elevation
+ trajectory_list.append([x, y, z])
+
+ # Call compareTrajectories function to do the hard work
+ compareTrajectories(trajectory_list, names, legend)
+
+ return None
+
+
+def compareAllInfo(flight_list, names=None):
+ """Creates a plot with the altitude, velocity and acceleration of the
+ Flight objects passed via a Python list.
+
+ Parameters
+ ----------
+ flight_list : list, array
+ List of FLight objects. The flights must be in the same reference frame.
+ names : list, optional
+ List of strings with the name of each trajectory inputted. The names must be in
+ the same order as the trajectories in flight_list
+
+ Returns
+ -------
+ None
+
+ """
+ return None
+
+
def create_dispersion_dictionary(filename):
"""Creates a dictionary with the rocket data provided by a .csv file.
File should be organized in four columns: attribute_class, parameter_name,
diff --git a/setup.py b/setup.py
index b94d4e669..6237e49a4 100644
--- a/setup.py
+++ b/setup.py
@@ -17,7 +17,6 @@
"pytz",
"timezonefinder",
"simplekml",
- "jsonpickle",
],
maintainer="RocketPy Developers",
author="Giovani Hidalgo Ceotto",
diff --git a/tests/test_flight.py b/tests/test_flight.py
index 9c66a7b3a..b08763339 100644
--- a/tests/test_flight.py
+++ b/tests/test_flight.py
@@ -249,9 +249,9 @@ def mainTrigger(p, y):
10,
0.0,
0.0,
- 0.0,
- 0.0,
- 0.0,
+ 0.9990482215818579,
+ -0.043619387365336,
+ -0.0,
0.0,
0.0,
0.0,
@@ -533,7 +533,7 @@ def test_export_data():
def test_export_KML():
- "Tests weather the method Flight.exportKML is working as intended"
+ "Tests wether the method Flight.exportKML is working as intended"
test_env = Environment(
railLength=5,