diff --git a/docs/notebooks/dispersion_analysis/dispersion_class_usage.ipynb b/docs/notebooks/dispersion_analysis/dispersion_class_usage.ipynb new file mode 100644 index 000000000..41278e098 --- /dev/null +++ b/docs/notebooks/dispersion_analysis/dispersion_class_usage.ipynb @@ -0,0 +1,649 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Monte Carlo Dispersion Analysis with the Dispersion Class" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally the Monte Carlo Simulations can be performed using a dedicated class called Dispersion. This class is a wrapper for the Monte Carlo Simulations, and it is the recommended way to perform the simulations. Say goodbye to the long and tedious process of creating the Monte Carlo Simulations throughout jupyter notebooks!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, let's import the necessary libraries, including the newest Dispersion class!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from rocketpy import Environment, Rocket, SolidMotor, Flight, Dispersion\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If you are using Jupyter Notebooks, it is recommended to run the following line to make matplotlib plots which will be shown later interactive and higher quality." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Dispersion class allows us to perform Monte Carlo Simulations in a very simple way.\n", + "We just need to create an instance of the class, and then call the method run() to perform the simulations.\n", + "The class has a lot of capabilities, but we will only use a few of them in this example.\n", + "We encourage you to check the documentation of the class to learn more about the Dispersion.\n", + "\n", + "Also, you can check RocketPy's main reference for a better conceptual understanding \n", + "of the Monte Carlo Simulations: [RocketPy: Six Degree-of-Freedom Rocket Trajectory Simulator](https://doi.org/10.1061/(ASCE)AS.1943-5525.0001331)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will describe two different options of usage:\n", + "- Using a Flight object as input, speeding up the process of creating the simulations. (currently not completely described in this notebook)\n", + "- Using a dictionary as input, including mean values and uncertainties for each parameter \n", + "\n", + "You need only one of them to get started. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1st Option -> Use your Flight object as input for the Dispersion class" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Creating an Environment for 'Ponte de Sôr', Portugal" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Env = Environment(\n", + " railLength=5.2, latitude=39.389700, longitude=-8.288964, elevation=113\n", + ")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To get weather data from the GFS forecast, available online, we run the following lines.\n", + "\n", + "First, we set tomorrow's date." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "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" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then, we tell Env to use a GFS forecast to get the atmospheric conditions for flight.\n", + "\n", + "Don't mind the warning, it just means that not all variables, such as wind speed or atmospheric temperature, are available at all altitudes given by the forecast." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Env.setAtmosphericModel(type=\"Forecast\", file=\"GFS\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Env.info()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Creating a Motor for the Rocket" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We define a motor for the rocket, using the data from the manufacturer, and following\n", + "the [RocketPy's documentation](https://docs.rocketpy.org/en/latest/user/index.html)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "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", + ")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Pro75M1670.info()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Creating a Rocket" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Calisto = Rocket(\n", + " motor=Pro75M1670,\n", + " radius=127 / 2000,\n", + " mass=19.197 - 2.956,\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", + "Calisto.setRailButtons([0.2, -0.5])\n", + "\n", + "NoseCone = Calisto.addNose(length=0.55829, kind=\"vonKarman\", distanceToCM=0.71971)\n", + "\n", + "FinSet = Calisto.addFins(\n", + " 4, span=0.100, rootChord=0.120, tipChord=0.040, distanceToCM=-1.04956\n", + ")\n", + "\n", + "Tail = Calisto.addTail(\n", + " topRadius=0.0635, bottomRadius=0.0435, length=0.060, distanceToCM=-1.194656\n", + ")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Additionally, we set parachutes for our Rocket, as well as the trigger functions for the deployment of such parachutes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "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 < 500 + 100 m (+100 due to surface elevation).\n", + " return True if y[5] < 0 and y[2] < 500 + 100 else False\n", + "\n", + "\n", + "Main = Calisto.addParachute(\n", + " \"Main\",\n", + " CdS=10.0,\n", + " trigger=mainTrigger,\n", + " samplingRate=105,\n", + " lag=1.5,\n", + " noise=(0, 8.3, 0.5),\n", + ")\n", + "\n", + "Drogue = Calisto.addParachute(\n", + " \"Drogue\",\n", + " CdS=1.0,\n", + " trigger=drogueTrigger,\n", + " samplingRate=105,\n", + " lag=1.5,\n", + " noise=(0, 8.3, 0.5),\n", + ")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Calisto.info()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Simulate single flight" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "TestFlight = Flight(\n", + " rocket=Calisto,\n", + " environment=Env,\n", + " inclination=84,\n", + " heading=133,\n", + ")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And we can visualize the flight trajectory:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "TestFlight.plot3dTrajectory()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Starting the Monte Carlo Simulations" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, let's invoke the Dispersion class, we only need a filename to initialize it.\n", + "The filename will be used either to save the results of the simulations or to load them\n", + "from a previous ran simulation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "TestDispersion = Dispersion(filename=\"dispersion_analysis_outputs/disp_class_example\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then, we can run the simulations using the method Dispersion.run_dispersion().\n", + "But before that, we need to set some simple parameters for the simulations.\n", + "We will set them by using a dictionary, which is one of the simplest way to do it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#TODO: explain that this is std only\n", + "disp_dictionary = {\n", + " # Solid Motor Parameters\n", + " \"burnOutTime\": 0.2,\n", + " \"totalImpulse\": 0.033 * Pro75M1670.totalImpulse,\n", + " # Rocket Parameters\n", + " \"mass\": 0.100,\n", + " \"radius\": 0.001,\n", + " \"powerOffDrag\": 0.033, # Multiplier\n", + " \"powerOnDrag\": 0.033, # Multiplier\n", + " \"parachute_Main_CdS\": 1,\n", + " \"parachute_Drogue_CdS\": 0.1,\n", + " \"parachute_Main_lag\": 0.1,\n", + " \"parachute_Drogue_lag\": 0.1,\n", + " # Flight Parameters\n", + " \"inclination\": 1,\n", + " \"heading\": 2,\n", + "}\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, let's iterate over the simulations and export the data from each flight simulation!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "TestDispersion.run_dispersion(\n", + " number_of_simulations=50,\n", + " dispersion_dictionary=disp_dictionary,\n", + " flight=TestFlight,\n", + " append=False,\n", + ")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Visualizing the results" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we finally have the results of our Monte Carlo simulations loaded!\n", + "Let's play with them." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, we can print numerical information regarding the results of the simulations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "TestDispersion.import_results()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "TestDispersion.print_results()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Also, we can visualize histograms of such results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "TestDispersion.allInfo()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Export to kml so it can be visualized in Google Earth" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "TestDispersion.exportEllipsesToKML(\n", + " filename=\"dispersion_analysis_outputs/disp_class_example.kml\",\n", + " origin_lat=Env.latitude,\n", + " origin_lon=Env.longitude,\n", + " type=\"impact\",\n", + ")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2nd Option -> Running by using only a dictionary of parameters" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This second option allow us to perform the Monte Carlo Simulations without the need of a Flight object. This is useful when we want to perform the simulations for a rocket that we don't have a Flight object for, or when we want to perform the simulations for a rocket that we have a Flight object for, but we want to change some parameters of the simulations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "TestDispersion2 = Dispersion(filename=\"dispersion_analysis_outputs/disp_class_example2\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "aerodynamic_surfaces = Calisto.aerodynamicSurfaces\n", + "nose, fins, tail = aerodynamic_surfaces\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dispersion_dictionary2 = {\n", + " # Environment Parameters\n", + " \"railLength\": (Env.railLength, 0.001),\n", + " \"date\": [Env.date],\n", + " \"datum\": [\"WSG84\"],\n", + " \"elevation\": (Env.elevation, 10),\n", + " \"gravity\": (Env.gravity, 0),\n", + " \"latitude\": (Env.latitude, 0),\n", + " \"longitude\": (Env.longitude, 0),\n", + " \"timeZone\": [str(Env.timeZone)],\n", + " # Solid Motor Parameters\n", + " \"burnOutTime\": (Pro75M1670.burnOutTime, 0.2),\n", + " \"grainDensity\": (Pro75M1670.grainDensity, 0.1 * Pro75M1670.grainDensity),\n", + " \"grainInitialHeight\": (Pro75M1670.grainInitialHeight, 0.001),\n", + " \"grainInitialInnerRadius\": (Pro75M1670.grainInitialInnerRadius, 0.001),\n", + " \"grainNumber\": [Pro75M1670.grainNumber],\n", + " \"grainOuterRadius\": (Pro75M1670.grainOuterRadius, 0.001),\n", + " \"grainSeparation\": (Pro75M1670.grainSeparation, 0.001),\n", + " \"nozzleRadius\": (Pro75M1670.nozzleRadius, 0.001),\n", + " \"throatRadius\": (Pro75M1670.throatRadius, 0.001),\n", + " \"thrustSource\": [Pro75M1670.thrustSource],\n", + " \"totalImpulse\": (Pro75M1670.totalImpulse, 0.033 * Pro75M1670.totalImpulse),\n", + " # Rocket Parameters\n", + " \"mass\": (Calisto.mass, 0.100),\n", + " \"radius\": (Calisto.radius, 0.001),\n", + " \"distanceRocketNozzle\": (Calisto.distanceRocketNozzle, 0.010),\n", + " \"distanceRocketPropellant\": (Calisto.distanceRocketPropellant, 0.010),\n", + " \"inertiaI\": (Calisto.inertiaI, Calisto.inertiaI * 0.1),\n", + " \"inertiaZ\": (Calisto.inertiaZ, Calisto.inertiaZ * 0.1),\n", + " \"powerOffDrag\": (1, 0.033),\n", + " \"powerOnDrag\": (1, 0.033),\n", + " \"nose_name_kind\": [nose.kind],\n", + " \"nose_name_length\": (nose.length, 0.001),\n", + " \"nose_name_distanceToCM\": (nose.distanceToCM, 0.010),\n", + " \"finSet_name_n\": [fins.n],\n", + " \"finSet_name_rootChord\": (fins.rootChord, 0.001),\n", + " \"finSet_name_tipChord\": (fins.tipChord, 0.001),\n", + " \"finSet_name_span\": (fins.span, 0.001),\n", + " \"finSet_name_distanceToCM\": (fins.distanceToCM, 0.010),\n", + " \"finSet_name_airfoil\": [fins.airfoil],\n", + " \"tail_name_topRadius\": (tail.topRadius, 0.001),\n", + " \"tail_name_bottomRadius\": (tail.bottomRadius, 0.001),\n", + " \"tail_name_length\": (tail.length, 0.001),\n", + " \"tail_name_distanceToCM\": (tail.distanceToCM, 0.010),\n", + " \"parachute_Main_CdS\": (10,2),\n", + " \"parachute_Main_trigger\": mainTrigger,\n", + " \"parachute_Drogue_CdS\": (1,0.3),\n", + " \"parachute_Drogue_trigger\": drogueTrigger,\n", + " # Flight Parameters\n", + " \"inclination\": [85],\n", + " \"heading\": [90],\n", + "}\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "TestDispersion2.run_dispersion(\n", + " number_of_simulations=50,\n", + " dispersion_dictionary=dispersion_dictionary2,\n", + ")\n", + "# TODO: NEeds to be tested, apparently is not capturing the parachutes correctly\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "TestDispersion2.import_results()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And finally, we can export the ellipses of the results to a .kml file so it can be opened on Google Earth" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "TestDispersion2.print_results()\n" + ] + } + ], + "metadata": { + "hide_input": false, + "kernelspec": { + "display_name": "Python 3.10.5 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.5 (tags/v3.10.5:f377153, Jun 6 2022, 16:14:13) [MSC v.1929 64 bit (AMD64)]" + }, + "vscode": { + "interpreter": { + "hash": "26de051ba29f2982a8de78e945f0abaf191376122a1563185a90213a26c5da77" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/user/index.rst b/docs/user/index.rst index a6e1c88d6..65a4d68da 100644 --- a/docs/user/index.rst +++ b/docs/user/index.rst @@ -12,6 +12,7 @@ Welcome to RocketPy's user documentation! ../notebooks/environment_analysis_class_usage.ipynb ../notebooks/environment_analysis_EuroC_example.ipynb ../notebooks/dispersion_analysis/dispersion_analysis.ipynb + ../notebooks/dispersion_analysis/dispersion_class_usage.ipynb ../notebooks/utilities_usage.ipynb ../matlab/matlab.rst diff --git a/getting_started Dispersion.ipynb b/getting_started Dispersion.ipynb deleted file mode 100644 index 9e329f7db..000000000 --- a/getting_started Dispersion.ipynb +++ /dev/null @@ -1,672 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Getting Started" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "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": 2, - "metadata": {}, - "outputs": [], - "source": [ - "%load_ext autoreload\n", - "%autoreload 2" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "from rocketpy import *" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "If you are using Jupyter Notebooks, it is recommended to run the following line to make matplotlib plots which will be shown later interactive and higher quality." - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "%matplotlib widget" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "from datetime import datetime" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Running with dictionary only" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "# Set up parachutes. This rocket, named Valetudo, only has a drogue chute.\n", - "def drogueTrigger(p, y):\n", - " # Check if rocket is going down, i.e. if it has passed the apogee\n", - " vertical_velocity = y[5]\n", - " # Return true to activate parachute once the vertical velocity is negative\n", - " return True if vertical_velocity < 0 else False" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "outputs": [], - "source": [ - "analysis_parameters = {\n", - " # Mass Details\n", - " \"mass\": (\n", - " 8.257,\n", - " 0.001,\n", - " ), # Rocket's dry mass (kg) and its uncertainty (standard deviation)\n", - " # Propulsion Details - run help(SolidMotor) for more information\n", - "\n", - " # \"thrust\": [\"docs/notebooks/dispersion_analysis/dispersion_analysis_inputs/thrustCurve.csv\"],\n", - " # \"totalImpulse\": (1415.15, 35.3), # Motor total impulse (N*s)\n", - " # \"burnOutTime\": (5.274, 1), # Motor burn out time (s)\n", - " # \"nozzleRadius\": (21.642 / 1000, 0.5 / 1000), # Motor's nozzle radius (m)\n", - " # \"throatRadius\": (8 / 1000, 0.5 / 1000), # Motor's nozzle throat radius (m)\n", - " # \"grainNumber\":[6],\n", - " # \"grainSeparation\": (\n", - " # 6 / 1000,\n", - " # 1 / 1000,\n", - " # ), # Motor's grain separation (axial distance between two grains) (m)\n", - " # \"grainDensity\": (1707, 50), # Motor's grain density (kg/m^3)\n", - " # \"grainOuterRadius\": (21.4 / 1000, 0.375 / 1000), # Motor's grain outer radius (m)\n", - " # \"grainInitialInnerRadius\": (\n", - " # 9.65 / 1000,\n", - " # 0.375 / 1000,\n", - " # ), # Motor's grain inner radius (m)\n", - " # \"grainInitialHeight\": (120 / 1000, 1 / 1000), # Motor's grain height (m)\n", - " # \"interpolationMethod\":[\"linear\"],\n", - " # Aerodynamic Details - run help(Rocket) for more information\n", - " \"inertiaI\": (\n", - " 3.675,\n", - " 0.03675,\n", - " ), # Rocket's inertia moment perpendicular to its axis (kg*m^2)\n", - " \"inertiaZ\": (\n", - " 0.007,\n", - " 0.00007,\n", - " ), # Rocket's inertia moment relative to its axis (kg*m^2)\n", - " \"radius\": (40.45 / 1000, 0.001), # Rocket's radius (kg*m^2)\n", - " \"distanceRocketNozzle\": (\n", - " -1.024,\n", - " 0.001,\n", - " ), # Distance between rocket's center of dry mass and nozzle exit plane (m) (negative)\n", - " \"distanceRocketPropellant\": (\n", - " -0.571,\n", - " 0.001,\n", - " ), # Distance between rocket's center of dry mass and and center of propellant mass (m) (negative)\n", - " \"powerOffDrag\": (\n", - " 0.9081 / 1.05,\n", - " 0.033,\n", - " ), # Multiplier for rocket's drag curve. Usually has a mean value of 1 and a uncertainty of 5% to 10%\n", - " \"powerOnDrag\": (\n", - " 0.9081 / 1.05,\n", - " 0.033,\n", - " ), # Multiplier for rocket's drag curve. Usually has a mean value of 1 and a uncertainty of 5% to 10%\n", - " \"noseKind\": [\"Von Karman\"],\n", - " \"noseLength\": (0.274, 0.001), # Rocket's nose cone length (m)\n", - " \"noseDistanceToCM\": (\n", - " 1.134,\n", - " 0.001,\n", - " ), # Axial distance between rocket's center of dry mass and nearest point in its nose cone (m)\n", - " \"numberOfFins\": [4],\n", - " \"span\": (0.077, 0.0005), # Fin span (m)\n", - " \"rootChord\": (0.058, 0.0005), # Fin root chord (m)\n", - " \"tipChord\": (0.018, 0.0005), # Fin tip chord (m)\n", - " \"distanceToCM\": (\n", - " -0.906,\n", - " 0.001,\n", - " ), # Axial distance between rocket's center of dry mass and nearest point in its fin (m)\n", - " \"airfoil\": [None],\n", - "\n", - " \"noTail\": [True],\n", - "\n", - " \"positionFirstRailButton\": [0.0224],\n", - " \"positionSecondRailButton\": [-0.93],\n", - " \"railButtonAngularPosition\": [30],\n", - " \n", - " # Launch and Environment Details - run help(Environment) and help(Flight) for more information\n", - " \"inclination\": (\n", - " 84.7,\n", - " 1,\n", - " ), # Launch rail inclination angle relative to the horizontal plane (degrees)\n", - " \"heading\": (53, 2), # Launch rail heading relative to north (degrees)\n", - " \"railLength\": (5.7, 0.0005), # Launch rail length (m)\n", - " \"ensembleMember\": list(range(10)), # Members of the ensemble forecast to be used\n", - " # Parachute Details - run help(Rocket) for more information\n", - " \"parachuteNames\":[\"Drogue\"],\n", - " \"CdS\": [(\n", - " 0.349 * 1.3,\n", - " 0.07,\n", - " )], # Drag coefficient times reference area for the drogue chute (m^2)\n", - " \"trigger\":[[drogueTrigger]],\n", - " \"lag\": [(\n", - " 173,\n", - " 0.5,\n", - " )], # Time delay between parachute ejection signal is detected and parachute is inflated (s)\n", - " \"samplingRate\":[[105]],\n", - " \"noise_mean\": [[0]],\n", - " \"noise_sd\": [[8.3]],\n", - " \"noise_corr\": [[0.5]], \n", - "}" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "outputs": [], - "source": [ - "# Define basic Environment object\n", - "Env = Environment(\n", - " railLength=5.7, date=(2019, 8, 10, 21), latitude=-23.363611, longitude=-48.011389\n", - ")\n", - "Env.setElevation(668)\n", - "Env.maxExpectedHeight = 1500\n", - "Env.setAtmosphericModel(\n", - " type=\"Ensemble\",\n", - " file=\"docs/notebooks/dispersion_analysis/dispersion_analysis_inputs/LASC2019_reanalysis.nc\",\n", - " dictionary=\"ECMWF\",\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": 13, - "metadata": {}, - "outputs": [], - "source": [ - "motorDisp = SolidMotor(\n", - " thrustSource=\"docs/notebooks/dispersion_analysis/dispersion_analysis_inputs/thrustCurve.csv\",\n", - " burnOut=5.274,\n", - " grainNumber=6,\n", - " grainDensity=1707,\n", - " grainOuterRadius=21.4 / 1000,\n", - " grainInitialInnerRadius=9.65 / 1000,\n", - " grainInitialHeight=120 / 1000,\n", - " grainSeparation=6 / 1000,\n", - " nozzleRadius=21.642 / 1000,\n", - " throatRadius=8 / 1000,\n", - " reshapeThrustCurve=False,\n", - " interpolationMethod=\"linear\",\n", - " )" - ] - }, - { - "cell_type": "code", - "execution_count": 14, - "metadata": {}, - "outputs": [ - { - "ename": "AttributeError", - "evalue": "'NoneType' object has no attribute 'latitude'", - "output_type": "error", - "traceback": [ - "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[1;31mAttributeError\u001b[0m Traceback (most recent call last)", - "\u001b[1;32mc:\\Mateus\\GitHub\\RocketPy\\getting_started Dispersion.ipynb Célula: 13\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 2\u001b[0m modified_dispersion_dict \u001b[39m=\u001b[39m {i: j \u001b[39mfor\u001b[39;00m i, j \u001b[39min\u001b[39;00m analysis_parameters\u001b[39m.\u001b[39mitems()}\n\u001b[0;32m 3\u001b[0m Disp\u001b[39m.\u001b[39mmotor\u001b[39m=\u001b[39mmotorDisp\n\u001b[1;32m----> 4\u001b[0m Disp\u001b[39m.\u001b[39;49mprocessDispersionDict(modified_dispersion_dict)\n", - "File \u001b[1;32mc:\\Mateus\\GitHub\\RocketPy\\rocketpy\\Dispersion.py:294\u001b[0m, in \u001b[0;36mDispersion.processDispersionDict\u001b[1;34m(self, dispersionDict)\u001b[0m\n\u001b[0;32m 291\u001b[0m missing_input \u001b[39m=\u001b[39m \u001b[39mstr\u001b[39m(missing_input)\n\u001b[0;32m 292\u001b[0m \u001b[39m# Add to the dict\u001b[39;00m\n\u001b[0;32m 293\u001b[0m dispersionDict[missing_input] \u001b[39m=\u001b[39m [\n\u001b[1;32m--> 294\u001b[0m \u001b[39mgetattr\u001b[39;49m(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49menvironment, missing_input)\n\u001b[0;32m 295\u001b[0m ]\n\u001b[0;32m 297\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mall\u001b[39m(\n\u001b[0;32m 298\u001b[0m rocket_input \u001b[39min\u001b[39;00m dispersionDict \u001b[39mfor\u001b[39;00m rocket_input \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mrocket_inputs\n\u001b[0;32m 299\u001b[0m ):\n\u001b[0;32m 300\u001b[0m \u001b[39m# Iterate through missing inputs\u001b[39;00m\n\u001b[0;32m 301\u001b[0m \u001b[39mfor\u001b[39;00m missing_input \u001b[39min\u001b[39;00m \u001b[39mset\u001b[39m(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mrocket_inputs) \u001b[39m-\u001b[39m dispersionDict\u001b[39m.\u001b[39mkeys():\n", - "\u001b[1;31mAttributeError\u001b[0m: 'NoneType' object has no attribute 'latitude'" - ] - } - ], - "source": [ - "Disp = Dispersion('test')\n", - "modified_dispersion_dict = {i: j for i, j in analysis_parameters.items()}\n", - "Disp.motor=motorDisp\n", - "Disp.processDispersionDict(modified_dispersion_dict)" - ] - }, - { - "cell_type": "code", - "execution_count": 16, - "metadata": {}, - "outputs": [ - { - "ename": "AttributeError", - "evalue": "'SolidMotor' object has no attribute 'thrustSource'", - "output_type": "error", - "traceback": [ - "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[1;31mAttributeError\u001b[0m Traceback (most recent call last)", - "\u001b[1;32mc:\\Mateus\\GitHub\\RocketPy\\getting_started Dispersion.ipynb Célula: 14\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 1\u001b[0m Disp \u001b[39m=\u001b[39m Dispersion(\u001b[39m'\u001b[39m\u001b[39mtest\u001b[39m\u001b[39m'\u001b[39m)\n\u001b[1;32m----> 2\u001b[0m Disp\u001b[39m.\u001b[39;49mrunDispersion(\u001b[39m100\u001b[39;49m,analysis_parameters,Env,motor\u001b[39m=\u001b[39;49mmotorDisp)\n", - "File \u001b[1;32mc:\\Mateus\\GitHub\\RocketPy\\rocketpy\\Dispersion.py:474\u001b[0m, in \u001b[0;36mDispersion.runDispersion\u001b[1;34m(self, number_of_simulations, dispersionDict, environment, flight, motor, rocket, distributionType, image, realLandingPoint)\u001b[0m\n\u001b[0;32m 471\u001b[0m \u001b[39m# Creates copy of dispersionDict that will be altered\u001b[39;00m\n\u001b[0;32m 472\u001b[0m modified_dispersion_dict \u001b[39m=\u001b[39m {i: j \u001b[39mfor\u001b[39;00m i, j \u001b[39min\u001b[39;00m dispersionDict\u001b[39m.\u001b[39mitems()}\n\u001b[1;32m--> 474\u001b[0m analysis_parameters \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mprocessDispersionDict(modified_dispersion_dict)\n\u001b[0;32m 476\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mdistribuitionFunc \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39msetDistributionFunc(distributionType)\n\u001b[0;32m 477\u001b[0m \u001b[39m# Basic analysis info\u001b[39;00m\n\u001b[0;32m 478\u001b[0m \n\u001b[0;32m 479\u001b[0m \u001b[39m# Create data files for inputs, outputs and error logging\u001b[39;00m\n", - "File \u001b[1;32mc:\\Mateus\\GitHub\\RocketPy\\rocketpy\\Dispersion.py:304\u001b[0m, in \u001b[0;36mDispersion.processDispersionDict\u001b[1;34m(self, dispersionDict)\u001b[0m\n\u001b[0;32m 302\u001b[0m missing_input \u001b[39m=\u001b[39m \u001b[39mstr\u001b[39m(missing_input)\n\u001b[0;32m 303\u001b[0m \u001b[39m# Add to the dict\u001b[39;00m\n\u001b[1;32m--> 304\u001b[0m dispersionDict[missing_input] \u001b[39m=\u001b[39m [\u001b[39mgetattr\u001b[39;49m(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mrocket, missing_input)]\n\u001b[0;32m 306\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mall\u001b[39m(\n\u001b[0;32m 307\u001b[0m motor_input \u001b[39min\u001b[39;00m dispersionDict \u001b[39mfor\u001b[39;00m motor_input \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39msolidmotor_inputs\n\u001b[0;32m 308\u001b[0m ):\n\u001b[0;32m 309\u001b[0m \u001b[39m# Iterate through missing inputs\u001b[39;00m\n\u001b[0;32m 310\u001b[0m \u001b[39mfor\u001b[39;00m missing_input \u001b[39min\u001b[39;00m \u001b[39mset\u001b[39m(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39msolidmotor_inputs) \u001b[39m-\u001b[39m dispersionDict\u001b[39m.\u001b[39mkeys():\n", - "\u001b[1;31mAttributeError\u001b[0m: 'SolidMotor' object has no attribute 'thrustSource'" - ] - } - ], - "source": [ - "Disp = Dispersion('test')\n", - "Disp.runDispersion(100,analysis_parameters,Env,motor=motorDisp)" - ] - }, - { - "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": null, - "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", - "\n", - "First, we set tomorrow's date." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "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" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Then, we tell Env to use a GFS forecast to get the atmospheric conditions for flight.\n", - "\n", - "Don't mind the warning, it just means that not all variables, such as wind speed or atmospheric temperature, are available at all altitudes given by the forecast." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "Env.setAtmosphericModel(type=\"Forecast\", file=\"GFS\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can see what the weather will look like by calling the info method!" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "Env.info()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Creating a Motor\n", - "\n", - "A solid rocket motor is used in this case. To create a motor, the SolidMotor class is used and the required arguments are given.\n", - "\n", - "The SolidMotor class requires the user to have a thrust curve ready. This can come either from a .eng file for a commercial motor, such as below, or a .csv file from a static test measurement.\n", - "\n", - "Besides the thrust curve, other parameters such as grain properties and nozzle dimensions must also be given." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "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": "markdown", - "metadata": {}, - "source": [ - "To see what our thrust curve looks like, along with other import properties, we invoke the info method yet again. You may try the allInfo method if you want more information all at once!" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "Pro75M1670.info()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Creating a Rocket" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "A rocket is composed of several components. Namely, we must have a motor (good thing we have the Pro75M1670 ready), a couple of aerodynamic surfaces (nose cone, fins and tail) and parachutes (if we are not launching a missile).\n", - "\n", - "Let's start by initializing our rocket, named Calisto, supplying it with the Pro75M1670 engine, entering its inertia properties, some dimensions and also its drag curves." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "Calisto = Rocket(\n", - " motor=Pro75M1670,\n", - " radius=127 / 2000,\n", - " mass=19.197 - 2.956,\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", - "Calisto.setRailButtons([0.2, -0.5])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Adding Aerodynamic Surfaces" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now we define the aerodynamic surfaces. They are really straight forward." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "NoseCone = Calisto.addNose(length=0.55829, kind=\"vonKarman\", distanceToCM=0.71971)\n", - "\n", - "FinSet = Calisto.addFins(\n", - " 4, span=0.100, rootChord=0.120, tipChord=0.040, distanceToCM=-1.04956\n", - ")\n", - "\n", - "Tail = Calisto.addTail(\n", - " topRadius=0.0635, bottomRadius=0.0435, length=0.060, distanceToCM=-1.194656\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Adding Parachutes" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, we have parachutes! Calisto will have two parachutes, Drogue and Main.\n", - "\n", - "Both parachutes are activated by some special algorithm, which is usually really complex and a trade secret. Most algorithms are based on pressure sampling only, while some also use acceleration info.\n", - "\n", - "RocketPy allows you to define a trigger function which will decide when to activate the ejection event for each parachute. This trigger function is supplied with pressure measurement at a predefined sampling rate. This pressure signal is usually noisy, so artificial noise parameters can be given. Call help(Rocket.addParachute) for more details. Furthermore, the trigger function also receives the complete state vector of the rocket, allowing us to use velocity, acceleration or even attitude to decide when the parachute event should be triggered.\n", - "\n", - "Here, we define our trigger functions rather simply using Python. However, you can call the exact code which will fly inside your rocket as well." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "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", - "Main = Calisto.addParachute(\n", - " \"Main\",\n", - " CdS=10.0,\n", - " trigger=mainTrigger,\n", - " samplingRate=105,\n", - " lag=1.5,\n", - " noise=(0, 8.3, 0.5),\n", - ")\n", - "\n", - "Drogue = Calisto.addParachute(\n", - " \"Drogue\",\n", - " CdS=1.0,\n", - " trigger=drogueTrigger,\n", - " samplingRate=105,\n", - " lag=1.5,\n", - " noise=(0, 8.3, 0.5),\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Just be careful if you run this last cell multiple times! If you do so, your rocket will end up with lots of parachutes which activate together, which may cause problems during the flight simulation. We advise you to re-run all cells which define our rocket before running this, preventing unwanted old parachutes. Alternatively, you can run the following lines to remove parachutes.\n", - "\n", - "```python\n", - "Calisto.parachutes.remove(Drogue)\n", - "Calisto.parachutes.remove(Main)\n", - "```" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Simulating a Flight\n", - "\n", - "Simulating a flight trajectory is as simple as initializing a Flight class object givin the rocket and environnement set up above as inputs. The launch rail inclination and heading are also given here." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "TestFlight = Flight(rocket=Calisto, environment=Env, inclination=85, heading=0)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Analyzing the Results\n", - "\n", - "RocketPy gives you many plots, thats for sure! They are divided into sections to keep them organized. Alternatively, see the Flight class documentation to see how to get plots for specific variables only, instead of all of them at once." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "TestFlight.allInfo()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Export Flight Trajectory to a .kml file so it can be opened on Google Earth" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "TestDispersion = Dispersion('teste')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "dispersionDict={\n", - " 'rocketMass':(19.197 - 2.956, 0.2),\n", - " 'radius': 0.1* 127 / 2000,\n", - " 'burnOut':(3.9, 0.5), \n", - " 'parachuteNames': [\"Main\",\"Drogue\"],\n", - " 'CdS': [(10,2),(1)],\n", - " 'trigger':[[mainTrigger],[drogueTrigger]]\n", - " } " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "TestDispersion.runDispersion(number_of_simulations=10,dispersionDict=dispersionDict, \n", - " flight=TestFlight, image='dispersion_analysis_inputs/Valetudo_basemap_final.jpg')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "TestDispersion.info()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "hide_input": false, - "kernelspec": { - "display_name": "Python 3.10.5 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.5" - }, - "vscode": { - "interpreter": { - "hash": "26de051ba29f2982a8de78e945f0abaf191376122a1563185a90213a26c5da77" - } - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/rocketpy/AeroSurfaces.py b/rocketpy/AeroSurfaces.py index fcd3064a9..dc4bf98de 100644 --- a/rocketpy/AeroSurfaces.py +++ b/rocketpy/AeroSurfaces.py @@ -215,6 +215,7 @@ class Fins(ABC): Geometrical relation used to simplify lift and roll calculations Fins.liftInterferenceFactor : float Factor of Fin-Body interference in the lift coefficient + """ def __init__( @@ -230,14 +231,17 @@ def __init__( ): """Initialize Fins class. + Parameters ---------- n : int Number of fins, from 2 to infinity. + rootChord : int, float Fin root chord in meters. span : int, float Fin span in meters. + distanceToCM : int, float Fin set position relative to rocket unloaded center of mass, considering positive direction from center of mass to @@ -248,6 +252,7 @@ def __init__( cantAngle : int, float, optional Fins cant angle with respect to the rocket centerline. Must be given in degrees. + airfoil : tuple, optional Default is null, in which case fins will be treated as flat plates. Otherwise, if tuple, fins will be considered as airfoils. The @@ -277,10 +282,12 @@ def __init__( # Store values self.n = n self.rocketRadius = rocketRadius + self.airfoil = airfoil self.distanceToCM = distanceToCM self.cantAngle = cantAngle self.rootChord = rootChord + self.span = span self.name = name self.d = d @@ -328,6 +335,7 @@ def evaluateLiftCoefficient(self): if not self.airfoil: # Defines clalpha2D as 2*pi for planar fins clalpha2D = Function(lambda mach: 2 * np.pi / self.__beta(mach)) + else: # Defines clalpha2D as the derivative of the # lift coefficient curve for a specific airfoil @@ -792,6 +800,7 @@ def __init__( + 4 * (rootChord + 2 * tipChord) * rocketRadius * span**2 + 6 * (rootChord + tipChord) * span * rocketRadius**2 ) / 12 + rollDampingInterferenceFactor = 1 + ( ((tau - λ) / (tau)) - ((1 - λ) / (tau - 1)) * np.log(tau) ) / ( @@ -855,10 +864,12 @@ def evaluateCenterOfPressure(self): - self.rootChord * self.tipChord / (self.rootChord + self.tipChord) ) ) + self.cpx = 0 self.cpy = 0 self.cpz = cpz self.cp = (self.cpx, self.cpy, self.cpz) + return self.cp def draw(self): @@ -1024,6 +1035,7 @@ class EllipticalFins(Fins): Geometrical relation used to simplify lift and roll calculations Fins.liftInterferenceFactor : float Factor of Fin-Body interference in the lift coefficient + """ def __init__( @@ -1039,9 +1051,11 @@ def __init__( ): """Initialize EllipticalFins class. + Parameters ---------- n : int + Number of fins, from 2 to infinity. rootChord : int, float Fin root chord in meters. @@ -1090,9 +1104,11 @@ def __init__( name : str Name of fin set. + Returns ------- None + """ super().__init__( @@ -1152,6 +1168,7 @@ def __init__( * (span / 3 + np.pi * rocketRadius / 4) * (span**2 - rocketRadius**2) ) + rollForcingInterferenceFactor = (1 / np.pi**2) * ( (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) + ((np.pi * (tau**2 + 1) ** 2) / (tau**2 * (tau - 1) ** 2)) @@ -1278,6 +1295,7 @@ def draw(self): class Tail: + """Class that defines a tail for the rocket. Currently only accepts conical tails. @@ -1321,6 +1339,7 @@ class Tail: Tail.surfaceArea : float Surface area of the tail. Has the unit of meters squared. + """ def __init__( diff --git a/rocketpy/Dispersion.py b/rocketpy/Dispersion.py index 204849ff9..c94efcfdf 100644 --- a/rocketpy/Dispersion.py +++ b/rocketpy/Dispersion.py @@ -1,13 +1,13 @@ # -*- coding: utf-8 -*- -__author__ = "Mateus Stano Junqueira, Sofia Lopes Suesdek Rocha, Abdulklech Sorban" -__copyright__ = "Copyright 20XX, Projeto Jupiter" +__author__ = "Mateus Stano Junqueira, Sofia Lopes Suesdek Rocha, Guilherme Fernandes Alves, Bruno Abdulklech Sorban" +__copyright__ = "Copyright 20XX, RocketPy Team" __license__ = "MIT" import math import traceback -import warnings +import types from time import process_time, time import matplotlib.pyplot as plt @@ -18,12 +18,13 @@ from matplotlib.patches import Ellipse from numpy.random import * +from .AeroSurfaces import EllipticalFins, NoseCone, Tail, TrapezoidalFins from .Environment import Environment from .Flight import Flight from .Function import Function -from .Motor import HybridMotor, SolidMotor +from .Motor import SolidMotor from .Rocket import Rocket -from .utilities import invertedHaversine +from .supplement import invertedHaversine class Dispersion: @@ -31,31 +32,49 @@ class Dispersion: """Monte Carlo analysis to predict probability distributions of the rocket's landing point, apogee and other relevant information. - Attributes + Parameters ---------- - Parameters: - Dispersion.filename: string - When running a new simulation, this attribute represents the initial - part of the export filenames (e.g. 'filename.disp_outputs.txt'). - When analyzing the results of a previous simulation, this attribute - shall be the filename containing the outputs of a dispersion calculation. - Dispersion.image: string - Launch site PNG file to be plotted along with the dispersion ellipses. - Attribute needed to run a new simulation. - Dispersion.realLandingPoint: tuple - Rocket's experimental landing point relative to launch point. - Dispersion.N: integer - Number of simulations in an output file. - Other classes: - Dispersion.environment: Environment - Launch environment. - Attribute needed to run a new simulation, when Dispersion.flight remains unchanged. - Dispersion.motor: Motor - Rocket's motor. - Attribute needed to run a new simulation, when Dispersion.flight remains unchanged. - Dispersion.rocket: Rocket - Rocket with nominal values. - Attribute needed to run a new simulation, when Dispersion.flight remains unchanged. + filename : string + The name of the file containing the data to be used in the analysis. + + Attributes + ---------- # TODO: finish documentation + Dispersion.filename : string + Directory and name of dispersion files. When running a new simulation, + this parameter represents the initial part of the export filenames + (e.g. 'filename.disp_outputs.txt'). When analyzing the results of a + previous simulation, this parameter shall be the .txt filename containing + the outputs of a previous ran dispersion analysis. + Dispersion.inputs_dict : dict + Contains information regarding the input arguments of the + classes. Its keys refers to each of the classes that must be defined during + the simulation. Its values are dictionaries where the keys are the input + arguments of each class and the values are either the string "required" + (meaning it is not an optional argument) or the default value if that argument + is optional. + Dispersion.dispersion_results : dict + Holds dispersion results. + Dispersion.dispersion_dictionary : dict + Contains inputs to run dispersion + Dispersion.nose_names = [] + Dispersion.finSet_names = [] + Dispersion.tail_names = [] + Dispersion.parachute_names = [] + Dispersion.distributionFunc = None + Dispersion.distribution_type = None + Dispersion.environment = None + Dispersion.flight = None + Dispersion.motor = None + Dispersion.rocket = None + Dispersion.rocket_dispersion = None + Dispersion.number_of_simulations = 0 + Dispersion.num_of_loaded_sims = 0 + Dispersion.start_time = 0 + + Dispersion.num_of_loaded_sims : int + The number of simulations loaded from the file. + Dispersion.num_of_sims : int + The number of simulations to be performed. """ def __init__( @@ -70,44 +89,9 @@ def __init__( When running a new simulation, this parameter represents the initial part of the export filenames (e.g. 'filename.disp_outputs.txt'). When analyzing the results of a previous simulation, this parameter - shall be the .txt filename containing the outputs of a dispersion calculation. - number_of_simulations: integer, needed when running a new simulation - Number of simulations desired, must be greater than zero. - Default is zero. - flight: Flight - Original rocket's flight with nominal values. - Parameter needed to run a new simulation, when environment, - motor and rocket remain unchanged. - Default is None. - image: string, needed when running a new simulation - Launch site PNG file to be plotted along with the dispersion ellipses. - dispersionDict: dictionary, optional - Contains the information of which environment, motor, rocket and flight variables - will vary according to its standard deviation. - Format {'parameter0': (nominal value, standard deviation), 'parameter1': - (nominal value, standard deviation), ...} - (e.g. {'rocketMass':(20, 0.2), - 'burnOut': (3.9, 0.3), 'railLength': (5.2, 0.05)}) - Default is {}. - environment: Environment - Launch environment. - Parameter needed to run a new simulation, when Dispersion.flight remains unchanged. - Default is None. - motor: Motor, optional - Rocket's motor. - Parameter needed to run a new simulation, when Dispersion.flight remains unchanged. - Default is None. - rocket: Rocket, optional - Rocket with nominal values. - Parameter needed to run a new simulation, when Dispersion.flight remains unchanged. - Default is None. - distributionType: string, optional - Determines which type of distribution will be applied to variable parameters and - its respective standard deviation. - Default is 'normal' - realLandingPoint: tuple, optional - Rocket's experimental landing point relative to launch point. - Format (horizontal distance, vertical distance) + shall be the .txt filename containing the outputs of a previous ran + dispersion analysis. + Returns ------- None @@ -115,309 +99,907 @@ def __init__( # Save and initialize parameters self.filename = filename - self.number_of_simulations = 0 - self.flight = None - self.dispersionDict = {} + + # Initialize variables to be used in the analysis in case of missing inputs + self.inputs_dict = { + "environment": { + "railLength": "required", + "gravity": 9.80665, + "date": None, + "latitude": 0, + "longitude": 0, + "elevation": 0, + "datum": "WGS84", + "timeZone": "UTC", + }, + "solidmotor": { + "thrust": "required", + "burnOutTime": "required", + "totalImpulse": 0, + "grainNumber": "required", + "grainDensity": "required", + "grainOuterRadius": "required", + "grainInitialInnerRadius": "required", + "grainInitialHeight": "required", + "grainSeparation": 0, + "nozzleRadius": 0.0335, + "throatRadius": 0.0114, + }, + "rocket": { + "mass": "required", + "inertiaI": "required", + "inertiaZ": "required", + "radius": "required", + "distanceRocketNozzle": "required", + "distanceRocketPropellant": "required", + "powerOffDrag": "required", + "powerOnDrag": "required", + }, + "nose": { + "length": "required", + "kind": "Von Karman", + "distanceToCM": "required", + "name": "Nose Cone", + }, + "fins": { + "n": "required", + "rootChord": "required", + "tipChord": "required", + "span": "required", + "distanceToCM": "required", + "cantAngle": 0, + "radius": None, + "airfoil": None, + }, + "tail": { + "topRadius": "required", + "bottomRadius": "required", + "length": "required", + "distanceToCM": "required", + }, + "railbuttons": { + "positionFirstRailButton": "required", + "positionSecondRailButton": "required", + "railButtonAngularPosition": 45, + }, + "parachute": { + "CdS": "required", + "trigger": "required", + "samplingRate": 100, + "lag": 0, + "noise": (0, 0, 0), + # "noiseStd": 0, + # "noiseCorr": 0, + }, + "flight": { + "inclination": 80, + "heading": 90, + "initialSolution": None, + "terminateOnApogee": False, + "maxTime": 600, + "maxTimeStep": np.inf, + "minTimeStep": 0, + "rtol": 1e-6, + "atol": 6 * [1e-3] + 4 * [1e-6] + 3 * [1e-3], + "timeOvershoot": True, + "verbose": False, + }, + } + # TODO: add all exportable attributes to this list + self.exportable_list = [ + "apogee", + "apogeeTime", + "apogeeX", + "apogeeY", + "executionTime", + "finalStaticMargin", + "frontalSurfaceWind", + "impactVelocity", + "initialStaticMargin", + "lateralSurfaceWind", + "maxAcceleration", + "maxAccelerationTime", + "maxSpeed", + "maxSpeedTime", + "numberOfEvents", + "outOfRailStaticMargin", + "outOfRailTime", + "outOfRailVelocity", + "tFinal", + "xImpact", + "yImpact", + ] + # Initialize variables so they can be accessed by MATLAB + self.dispersion_results = {} + self.dispersion_dictionary = {} + self.nose_names = [] + self.finSet_names = [] + self.tail_names = [] + self.parachute_names = [] + self.distributionFunc = None + self.distribution_type = None self.environment = None + self.flight = None self.motor = None self.rocket = None - self.distributionType = "normal" - self.image = None - self.realLandingPoint = None - self.parachuteTriggers = [] + self.rocket_dispersion = None + self.number_of_simulations = 0 + self.num_of_loaded_sims = 0 + self.start_time = 0 - def classCheck(self): - rocketAttributes = [] - rocketInputs = [] + return None + + def __set_distribution_function(self, distribution_type): + """Sets the distribution function to be used in the analysis. + + Parameters + ---------- + distribution_type : string + The type of distribution to be used in the analysis. It can be + 'uniform', 'normal', 'lognormal', etc. - def setDistributionFunc(self, distributionType): - if distributionType == "normal" or distributionType == None: + Returns + ------- + np.random distribution function + The distribution function to be used in the analysis. + """ + if distribution_type == "normal" or distribution_type == None: return normal - elif distributionType == "beta": + elif distribution_type == "beta": return beta - elif distributionType == "binomial": + elif distribution_type == "binomial": return binomial - elif distributionType == "chisquare": + elif distribution_type == "chisquare": return chisquare - elif distributionType == "dirichlet": + elif distribution_type == "dirichlet": return dirichlet - elif distributionType == "exponential": + elif distribution_type == "exponential": return exponential - elif distributionType == "f": + elif distribution_type == "f": return f - elif distributionType == "gamma": + elif distribution_type == "gamma": return gamma - elif distributionType == "geometric": + elif distribution_type == "geometric": return geometric - elif distributionType == "gumbel": + elif distribution_type == "gumbel": return gumbel - elif distributionType == "hypergeometric": + elif distribution_type == "hypergeometric": return hypergeometric - elif distributionType == "laplace": + elif distribution_type == "laplace": return laplace - elif distributionType == "logistic": + elif distribution_type == "logistic": return logistic - elif distributionType == "lognormal": + elif distribution_type == "lognormal": return lognormal - elif distributionType == "logseries": + elif distribution_type == "logseries": return logseries - elif distributionType == "multinomial": + elif distribution_type == "multinomial": return multinomial - elif distributionType == "multivariate_normal": + elif distribution_type == "multivariate_normal": return multivariate_normal - elif distributionType == "negative_binomial": + elif distribution_type == "negative_binomial": return negative_binomial - elif distributionType == "noncentral_chisquare": + elif distribution_type == "noncentral_chisquare": return noncentral_chisquare - elif distributionType == "noncentral_f": + elif distribution_type == "noncentral_f": return noncentral_f - elif distributionType == "pareto": + elif distribution_type == "pareto": return pareto - elif distributionType == "poisson": + elif distribution_type == "poisson": return poisson - elif distributionType == "power": + elif distribution_type == "power": return power - elif distributionType == "rayleigh": + elif distribution_type == "rayleigh": return rayleigh - elif distributionType == "standard_cauchy": + elif distribution_type == "standard_cauchy": return standard_cauchy - elif distributionType == "standard_exponential": + elif distribution_type == "standard_exponential": return standard_exponential - elif distributionType == "standard_gamma": + elif distribution_type == "standard_gamma": return standard_gamma - elif distributionType == "standard_normal": + elif distribution_type == "standard_normal": return standard_normal - elif distributionType == "standard_t": + elif distribution_type == "standard_t": return standard_t - elif distributionType == "triangular": + elif distribution_type == "triangular": return triangular - elif distributionType == "uneliform": + elif distribution_type == "uneliform": return uniform - elif distributionType == "vonmises": + elif distribution_type == "vonmises": return vonmises - elif distributionType == "wald": + elif distribution_type == "wald": return wald - elif distributionType == "weibull": + elif distribution_type == "weibull": return weibull - elif distributionType == "zipf": + elif distribution_type == "zipf": return zipf else: - warnings.warn("Distribution type not supported") - - def processDispersionDict(self, dispersionDict): - # Get parachutes names - if "parachuteNames" in dispersionDict: # TODO: use only dispersionDict - for i, name in enumerate(dispersionDict["parachuteNames"]): - if "CdS" in dispersionDict: - dispersionDict["parachute_" + name + "_CdS"] = dispersionDict[ - "CdS" - ][i] - if "trigger" in dispersionDict: - dispersionDict["parachute_" + name + "_trigger"] = dispersionDict[ - "trigger" - ][i] - if "samplingRate" in dispersionDict: - dispersionDict[ - "parachute_" + name + "_samplingRate" - ] = dispersionDict["samplingRate"][i] - if "lag" in dispersionDict: - dispersionDict["parachute_" + name + "_lag"] = dispersionDict[ - "lag" - ][i] - if "noise_mean" in dispersionDict: - dispersionDict[ - "parachute_" + name + "_noise_mean" - ] = dispersionDict["noise_mean"][i] - if "noise_sd" in dispersionDict: - dispersionDict["parachute_" + name + "_noise_std"] = dispersionDict[ - "noise_sd" - ][i] - if "noise_corr" in dispersionDict: - dispersionDict[ - "parachute_" + name + "_noise_corr" - ] = dispersionDict["noise_corr"][i] - dispersionDict.pop("CdS", None) - dispersionDict.pop("trigger", None) - dispersionDict.pop("samplingRate", None) - dispersionDict.pop("lag", None) - dispersionDict.pop("noise_mean", None) - dispersionDict.pop("noise_sd", None) - dispersionDict.pop("noise_corr", None) - self.parachute_names = dispersionDict.pop("parachuteNames", None) - - for parameter_key, parameter_value in dispersionDict.items(): - if isinstance(parameter_value, (tuple, list)): - continue - else: # if parameter_value is only the standard deviation - if "parachute" in parameter_key: - _, parachute_name, parameter = parameter_key.split("_") - dispersionDict[parameter_key] = ( - getattr( - self.rocket.parachutes[ - self.parachute_names.index(parachute_name) - ], - parameter, - ), - parameter_value, - ) - else: - if parameter_key in self.environment_inputs.keys(): - try: - dispersionDict[parameter_key] = ( - getattr(self.environment, parameter_key), - parameter_value, - ) - except Exception as E: - print("Error:") - print( - "Check if parameter was inputed correctly in dispersioDict." - + " Dictionary values must be either tuple or lists." - + " If single value, the correponding Class must " - + "must be inputed in Dispersion.runDispersion method.\n" - ) - print(traceback.format_exc()) - elif parameter_key in self.solidmotor_inputs.keys(): - try: - dispersionDict[parameter_key] = ( - getattr(self.motor, parameter_key), - parameter_value, - ) - except Exception as E: - print("Error:") - print( - "Check if parameter was inputed correctly in dispersioDict." - + " Dictionary values must be either tuple or lists." - + " If single value, the correponding Class must " - + "must be inputed in Dispersion.runDispersion method.\n" - ) - print(traceback.format_exc()) - elif parameter_key in self.rocket_inputs.keys(): - try: - dispersionDict[parameter_key] = ( - getattr(self.rocket, parameter_key), - parameter_value, - ) - except Exception as E: - print("Error:") - print( - "Check if parameter was inputed correctly in dispersioDict." - + " Dictionary values must be either tuple or lists." - + " If single value, the correponding Class must " - + "must be inputed in Dispersion.runDispersion method.\n" - ) - print(traceback.format_exc()) - elif parameter_key in self.flight_inputs.keys(): - try: - dispersionDict[parameter_key] = ( - getattr(self.flight, parameter_key), - parameter_value, - ) - except Exception as E: - print("Error:") - print( - "Check if parameter was inputed correctly in dispersioDict." - + " Dictionary values must be either tuple or lists." - + " If single value, the correponding Class must " - + "must be inputed in Dispersion.runDispersion method.\n" - ) - print(traceback.format_exc()) + raise ValueError( + "Distribution type not recognized. Please use a valid distribution type." + ) + + def __process_dispersion_dict(self, dictionary): + """Read the inputted dispersion dictionary from the run_dispersion method + and return a dictionary with the processed parameters, being ready to be + used in the dispersion simulation. + + Parameters + ---------- + dictionary : dict + Dictionary containing the parameters to be varied in the dispersion + simulation. The keys of the dictionary are the names of the parameters + to be varied, and the values can be either tuple or list. If the value + is a single value, the corresponding class of the parameter need to + be passed on the run_dispersion method. + + Returns + ------- + dictionary: dict + The modified dictionary with the processed parameters. + """ + + # Now we prepare all the parachute data + dictionary = self.__process_parachute_from_dict(dictionary) # Check remaining class inputs + # Environment + dictionary = self.__process_environment_from_dict(dictionary) + + # Motor + dictionary = self.__process_motor_from_dict(dictionary) + # Rocket + dictionary = self.__process_rocket_from_dict(dictionary) + + # Rail button + dictionary = self.__process_rail_buttons_from_dict(dictionary) + + # Aerodynamic Surfaces + dictionary = self.__process_aerodynamic_surfaces_from_dict(dictionary) + + # Flight + dictionary = self.__process_flight_from_dict(dictionary) + + # Finally check the inputted data + self.__check_inputted_values_from_dict(dictionary) + + return dictionary + + def __process_flight_from_dict(self, dictionary): + """Check if all the relevant inputs for the Flight class are present in + the dispersion dictionary, input the missing ones and return the modified + dictionary. + + Parameters + ---------- + dictionary : dict + Dictionary containing the parameters to be varied in the dispersion + simulation. The keys of the dictionary are the names of the parameters + to be varied, and the values can be either tuple or list. If the value + is a single value, the corresponding class of the parameter need to + be passed on the run_dispersion method. + + Returns + ------- + dictionary: dict + Modified dictionary with the processed flight parameters. + """ + # First check if all the inputs for the flight class are present in the + # dictionary, if not, input the missing ones if not all( - environment_input in dispersionDict - for environment_input in self.environment_inputs.keys() + flight_input in dictionary + for flight_input in self.inputs_dict["flight"].keys() ): # Iterate through missing inputs for missing_input in ( - set(self.environment_inputs.keys()) - dispersionDict.keys() + set(self.inputs_dict["flight"].keys()) - dictionary.keys() ): missing_input = str(missing_input) # Add to the dict try: - dispersionDict[missing_input] = [ - getattr(self.environment, missing_input) - ] - except: - # class was not inputed + # First try to catch value from the Flight object if passed + dictionary[missing_input] = [getattr(self.flight, missing_input)] + except AttributeError: + # Flight class was not inputted + # check if missing parameter is required + if self.inputs_dict["flight"][missing_input] == "required": + raise ValueError( + "The input {} is required for the Flight class.".format( + missing_input + ) + ) + else: # if not required, uses default value + dictionary[missing_input] = [ + self.inputs_dict["flight"][missing_input] + ] + + return dictionary + + def __process_rocket_from_dict(self, dictionary): + """Check if all the relevant inputs for the Rocket class are present in + the dispersion dictionary, input the missing ones and return the modified + dictionary. + + Parameters + ---------- + dictionary : dict + Dictionary containing the parameters to be varied in the dispersion + simulation. The keys of the dictionary are the names of the parameters + to be varied, and the values can be either tuple or list. If the value + is a single value, the corresponding class of the parameter need to + be passed on the run_dispersion method. + + Returns + ------- + dictionary: dict + Modified dictionary with the processed rocket parameters. + """ + + # Checks if there are any missing rocket inputs in dictionary + if not all( + rocket_input in dictionary + for rocket_input in self.inputs_dict["rocket"].keys() + ): + # Iterate through missing inputs + for missing_input in ( + set(self.inputs_dict["rocket"].keys()) - dictionary.keys() + ): + missing_input = str(missing_input) + # Add to the dict + try: + dictionary[missing_input] = [getattr(self.rocket, missing_input)] + except AttributeError: + # class was not inputted # checks if missing parameter is required - if self.environment_inputs[missing_input] == "required": - warnings.warn(f'Missing "{missing_input}" in dispersionDict') - else: # if not uses default value - dispersionDict[missing_input] = [ - self.environment_inputs[missing_input] + if self.inputs_dict["rocket"][missing_input] == "required": + raise ValueError( + "The input {} is required for the Rocket class.".format( + missing_input + ) + ) + else: # if not, uses default value + dictionary[missing_input] = [ + self.inputs_dict["rocket"][missing_input] ] + + return dictionary + + def __process_rail_buttons_from_dict(self, dictionary): + """Check if all the relevant inputs for the RailButtons class are present + in the dispersion dictionary, input the missing ones and return the + modified dictionary. + + Parameters + ---------- + dictionary : dict + Dictionary containing the parameters to be varied in the dispersion + simulation. The keys of the dictionary are the names of the parameters + to be varied, and the values can be either tuple or list. If the value + is a single value, the corresponding class of the parameter need to + be passed on the run_dispersion method. + + Returns + ------- + dictionary: dict + Modified dictionary with the processed rail buttons parameters. + """ + if not all( - motor_input in dispersionDict - for motor_input in self.solidmotor_inputs.keys() + rail_buttons_input in dictionary + for rail_buttons_input in self.inputs_dict["railbuttons"].keys() ): # Iterate through missing inputs for missing_input in ( - set(self.solidmotor_inputs.keys()) - dispersionDict.keys() + set(self.inputs_dict["railbuttons"].keys()) - dictionary.keys() ): missing_input = str(missing_input) # Add to the dict try: - dispersionDict[missing_input] = [getattr(self.motor, missing_input)] - except: - # class was not inputed + dictionary[missing_input] = [getattr(self.rocket, missing_input)] + except AttributeError: + # class was not inputted # checks if missing parameter is required - if self.solidmotor_inputs[missing_input] == "required": - warnings.warn(f'Missing "{missing_input}" in dispersionDict') - else: # if not uses default value - dispersionDict[missing_input] = [ - self.solidmotor_inputs[missing_input] + if self.inputs_dict["railbuttons"][missing_input] == "required": + raise ValueError( + "The input {} is required for the RailButtons class.".format( + missing_input + ) + ) + else: + # if not, uses default value + dictionary[missing_input] = [ + self.inputs_dict["railbuttons"][missing_input] ] + return dictionary + + def __process_aerodynamic_surfaces_from_dict(self, dictionary): + """Check if all the relevant inputs for the AerodynamicSurfaces class + are present in the dispersion dictionary, input the missing ones and + return the modified dictionary. + + Parameters + ---------- + dictionary : dict + Dictionary containing the parameters to be varied in the dispersion + simulation. The keys of the dictionary are the names of the parameters + to be varied, and the values can be either tuple or list. If the value + is a single value, the corresponding class of the parameter need to + be passed on the run_dispersion method. + """ + + # Check the number of fin sets, noses, and tails + self.nose_names, self.finSet_names, self.tail_names = [], [], [] + + # Get names from the input dictionary + for var in dictionary.keys(): + if "nose" in var: + self.nose_names.append(var.split("_")[1]) + elif "finSet" in var: + self.finSet_names.append(var.split("_")[1]) + elif "tail" in var: + self.tail_names.append(var.split("_")[1]) + # Get names from the rocket object + for surface in self.rocket.aerodynamicSurfaces: + if isinstance(surface, NoseCone): + self.nose_names.append(surface.name) + elif isinstance(surface, (TrapezoidalFins, EllipticalFins)): + self.finSet_names.append(surface.name) + elif isinstance(surface, Tail): + self.tail_names.append(surface.name) + # Remove duplicates + self.nose_names = list(set(self.nose_names)) + self.finSet_names = list(set(self.finSet_names)) + self.tail_names = list(set(self.tail_names)) + + # Check if there are enough arguments for each kind of aero surface + # Iterate through nose names + for name in self.nose_names: + # Iterate through aerodynamic surface available at rocket object + for surface in self.rocket.aerodynamicSurfaces: + if surface.name == name and isinstance(surface, NoseCone): + # in case we find the corresponding nose, check if all the + # inputs are present in the dictionary + for input in self.inputs_dict["nose"].keys(): + if f"nose_{name}_{input}" not in dictionary: + # Try to get the value from the rocket object + try: + dictionary[f"nose_{name}_{input}"] = [ + getattr(surface, input) + ] + except AttributeError: + # If not possible, check if the input is required + if self.inputs_dict["nose"][input] == "required": + raise ValueError( + "The input {} is required for the NoseCone class.".format( + input + ) + ) + + else: + # If not required, use default value + dictionary[f"nose_{name}_{input}"] = [ + self.inputs_dict["nose"][input] + ] + + # Iterate through fin sets names + for name in self.finSet_names: + # Iterate through aerodynamic surface available at rocket object + for surface in self.rocket.aerodynamicSurfaces: + if surface.name == name and isinstance( + surface, (TrapezoidalFins, EllipticalFins) + ): + # in case we find the corresponding fin set, check if all the + # inputs are present in the dictionary + for input in self.inputs_dict["fins"].keys(): + if f"finSet_{name}_{input}" not in dictionary: + # Try to get the value from the rocket object + try: + dictionary[f"finSet_{name}_{input}"] = [ + getattr(surface, input) + ] + except AttributeError: + # If not possible, check if the input is required + if self.inputs_dict["fins"][input] == "required": + raise ValueError( + "The input {} is required for the Fins class.".format( + input + ) + ) + else: + # If not required, use default value + dictionary[f"finSet_{name}_{input}"] = [ + self.inputs_dict["fins"][input] + ] + + # Iterate through tail names + for name in self.tail_names: + # Iterate through aerodynamic surface available at rocket object + for surface in self.rocket.aerodynamicSurfaces: + if surface.name == name and isinstance(surface, Tail): + # in case we find the corresponding tail, check if all the + # inputs are present in the dictionary + for input in self.inputs_dict["tail"].keys(): + if f"tail_{name}_{input}" not in dictionary: + # Try to get the value from the rocket object + try: + dictionary[f"tail_{name}_{input}"] = [ + getattr(surface, input) + ] + except AttributeError: + # If not possible, check if the input is required + if self.inputs_dict["tail"][input] == "required": + raise ValueError( + "The input {} is required for the Tail class.".format( + input + ) + ) + else: + # If not required, use default value + dictionary[f"tail_{name}_{input}"] = [ + self.inputs_dict["tail"][input] + ] + + return dictionary + + def __process_motor_from_dict(self, dictionary): + """Check if all the relevant inputs for the Motor class are present in + the dispersion dictionary, input the missing ones and return the modified + dictionary. + + Parameters + ---------- + dictionary : dict + Dictionary containing the parameters to be varied in the dispersion + simulation. The keys of the dictionary are the names of the parameters + to be varied, and the values can be either tuple or list. If the value + is a single value, the corresponding class of the parameter need to + be passed on the run_dispersion method. + + Returns + ------- + dictionary: dict + Modified dictionary with the processed rocket parameters. + """ + # TODO: Add more options of motor (i.e. Liquid and Hybrids) + if not all( - rocket_input in dispersionDict for rocket_input in self.rocket_inputs.keys() + motor_input in dictionary + for motor_input in self.inputs_dict["solidmotor"].keys() ): # Iterate through missing inputs - for missing_input in set(self.rocket_inputs.keys()) - dispersionDict.keys(): + for missing_input in ( + set(self.inputs_dict["solidmotor"].keys()) - dictionary.keys() + ): missing_input = str(missing_input) # Add to the dict try: - dispersionDict[missing_input] = [ - getattr(self.rocket, missing_input) - ] - except: - # class was not inputed + dictionary[missing_input] = [getattr(self.motor, missing_input)] + except AttributeError: + # class was not inputted # checks if missing parameter is required - if self.rocket_inputs[missing_input] == "required": - warnings.warn(f'Missing "{missing_input}" in dispersionDict') + if self.inputs_dict["solidmotor"][missing_input] == "required": + raise ValueError( + "The input {} is required for the SolidMotor class.".format( + missing_input + ) + ) else: # if not uses default value - dispersionDict[missing_input] = [ - self.rocket_inputs[missing_input] + dictionary[missing_input] = [ + self.inputs_dict["solidmotor"][missing_input] ] + return dictionary + + def __process_environment_from_dict(self, dictionary): + """Check if all the relevant inputs for the Environment class are present in + the dispersion dictionary, input the missing ones and return the modified + dictionary. + Parameters + ---------- + dictionary : dict + Dictionary containing the parameters to be varied in the dispersion + simulation. The keys of the dictionary are the names of the parameters + to be varied, and the values can be either tuple or list. If the value + is a single value, the corresponding class of the parameter need to + be passed on the run_dispersion method. + + Returns + ------- + dictionary: dict + Modified dictionary with the processed environment parameters. + """ + # Check if there is any missing input for the environment if not all( - flight_input in dispersionDict for flight_input in self.flight_inputs.keys() + environment_input in dictionary + for environment_input in self.inputs_dict["environment"].keys() ): # Iterate through missing inputs - for missing_input in set(self.flight_inputs.keys()) - dispersionDict.keys(): + for missing_input in ( + set(self.inputs_dict["environment"].keys()) - dictionary.keys() + ): missing_input = str(missing_input) # Add to the dict try: - dispersionDict[missing_input] = [ - getattr(self.flight, missing_input) + # First try to catch value from the Environment object if passed + dictionary[missing_input] = [ + getattr(self.environment, missing_input) ] - except: - # class was not inputed + except AttributeError: + # class was not inputted # checks if missing parameter is required - if self.flight_inputs[missing_input] == "required": - warnings.warn(f'Missing "{missing_input}" in dispersionDict') - else: # if not uses default value - dispersionDict[missing_input] = [ - self.flight_inputs[missing_input] + if self.inputs_dict["environment"][missing_input] == "required": + raise ValueError( + "The input {} is required for the Environment class.".format( + missing_input + ) + ) + else: # if not required, use default value + dictionary[missing_input] = [ + self.inputs_dict["environment"][missing_input] ] + return dictionary + + def __process_parachute_from_dict(self, dictionary): + """Check if all the relevant inputs for the Parachute class are present in + the dispersion dictionary, input the missing ones and return the modified + dictionary. + + Parameters + ---------- + dictionary : dict + Dictionary containing the parameters to be varied in the dispersion + simulation. The keys of the dictionary are the names of the parameters + to be varied, and the values can be either tuple or list. If the value + is a single value, the corresponding class of the parameter need to + be passed on the run_dispersion method. + + Returns + ------- + dictionary: dict + Modified dictionary with the processed parachute parameters. + """ + # Get the number and names of parachutes + self.parachute_names = [] + for key in dictionary.keys(): + if "parachute_" in key: + self.parachute_names.append(key.split("_")[1]) + # Remove duplicates + self.parachute_names = list(set(self.parachute_names)) + + # Check if there is enough arguments for defining each parachute + for name in self.parachute_names: + for parachute_input in self.inputs_dict["parachute"].keys(): + if ( + "parachute_{}_{}".format(name, parachute_input) + not in dictionary.keys() + ): + try: # Try to get the value from the Parachute object + if len(self.rocket.parachutes) > 0: + for chute in self.rocket.parachutes: + if getattr(chute, "name") == name: + dictionary[ + "parachute_{}_{}".format(name, parachute_input) + ] = [getattr(chute, parachute_input)] + else: + raise Exception + except Exception: # Class not passed + if self.inputs_dict["parachute"][parachute_input] == "required": + raise ValueError( + "The input {} is required for the Parachute class.".format( + parachute_input + ) + ) + else: + dictionary[ + "parachute_{}_{}".format(name, parachute_input) + ] = [ + self.inputs_dict["parachute"][parachute_input], + ] - return dispersionDict + return dictionary + + def __check_inputted_values_from_dict(self, dictionary): + """Check if the inputted values are valid. If not, raise an error. + + Parameters + ---------- + dictionary : dict + Dictionary containing the parameters to be varied in the dispersion + simulation. The keys of the dictionary are the names of the parameters + to be varied, and the values can be either tuple or list. If the value + is a single value, the corresponding class of the parameter need to + be passed on the run_dispersion method. + + Returns + ------- + dictionary: dict + The modified dictionary with the processed parameters. + """ + for parameter_key, parameter_value in dictionary.items(): + if isinstance(parameter_value, (tuple, list)): + # Everything is right with the data, we have mean and stdev + continue + + # In this case the parameter_value is only the std. dev. + ## First solve the parachute values + if "parachute" in parameter_key: + _, parachute_name, parameter = parameter_key.split("_") + # try: + if isinstance(parameter_value, types.FunctionType): + # Deal with trigger functions + dictionary[parameter_key] = parameter_value + else: + dictionary[parameter_key] = ( + getattr( + self.rocket.get_parachute_by_name(parachute_name), parameter + ), + parameter_value, + ) + + ## Second corrections - Environment + elif parameter_key in self.inputs_dict["environment"].keys(): + try: + dictionary[parameter_key] = ( + getattr(self.environment, parameter_key), + parameter_value, + ) + except AttributeError: + raise AttributeError( + f"Please check if the parameter {parameter_key} was inputted" + "correctly in dispersion_dictionary." + " Dictionary values must be either tuple or lists." + " If single value, the corresponding Class must" + " be inputted in the run_dispersion method." + ) + + ## Third corrections - SolidMotor + elif parameter_key in self.inputs_dict["solidmotor"].keys(): + try: + dictionary[parameter_key] = ( + getattr(self.motor, parameter_key), + parameter_value, + ) + except AttributeError: + raise AttributeError( + f"Please check if the parameter {parameter_key} was inputted" + "correctly in dispersion_dictionary." + " Dictionary values must be either tuple or lists." + " If single value, the corresponding Class must" + " be inputted in the run_dispersion method." + ) + + # Fourth correction - Rocket + elif parameter_key in self.inputs_dict["rocket"].keys(): + try: + dictionary[parameter_key] = ( + getattr(self.rocket, parameter_key), + parameter_value, + ) + except AttributeError: + raise AttributeError( + f"Please check if the parameter {parameter_key} was inputted" + "correctly in dispersion_dictionary." + " Dictionary values must be either tuple or lists." + " If single value, the corresponding Class must" + " be inputted in the run_dispersion method." + ) + + # Fifth correction - Flight + elif parameter_key in self.inputs_dict["flight"].keys(): + try: + dictionary[parameter_key] = ( + getattr(self.flight, parameter_key), + parameter_value, + ) + except AttributeError: + raise AttributeError( + f"Please check if the parameter {parameter_key} was inputted" + "correctly in dispersion_dictionary." + " Dictionary values must be either tuple or lists." + " If single value, the corresponding Class must" + " be inputted in the run_dispersion method." + ) + + # The analysis parameter dictionary is ready! Now we have mean and stdev + # for all parameters + + return dictionary + + def __create_initial_objects(self): + """Create rocketpy objects (Environment, Motor, Rocket, Flight) in case + that they were not created yet. + + Returns + ------- + None + """ + if self.environment is None: + try: + self.environment = Environment( + railLength=self.dispersion_dictionary["railLength"][0] + ) + except: + raise TypeError( + "Cannot define basic Environment. Missing railLength value in dictionary" + ) + if self.motor is None: + try: + self.motor = SolidMotor( + thrustSource=self.dispersion_dictionary["thrustSource"][0], + burnOut=self.dispersion_dictionary["burnOutTime"][0], + grainNumber=self.dispersion_dictionary["grainNumber"][0], + grainDensity=self.dispersion_dictionary["grainDensity"][0], + grainOuterRadius=self.dispersion_dictionary["grainOuterRadius"][0], + grainInitialInnerRadius=self.dispersion_dictionary[ + "grainInitialInnerRadius" + ][0], + grainInitialHeight=self.dispersion_dictionary["grainInitialHeight"][ + 0 + ], + ) + except: + raise TypeError( + "Cannot define basic SolidMotor. Missing required parameters in dictionary" + ) + if self.rocket is None: + try: + self.rocket = Rocket( + motor=self.motor, + mass=self.dispersion_dictionary["mass"][0], + radius=self.dispersion_dictionary["radius"][0], + inertiaI=self.dispersion_dictionary["inertiaI"][0], + inertiaZ=self.dispersion_dictionary["inertiaZ"][0], + distanceRocketPropellant=self.dispersion_dictionary[ + "distanceRocketPropellant" + ][0], + distanceRocketNozzle=self.dispersion_dictionary[ + "distanceRocketNozzle" + ][0], + powerOffDrag=self.dispersion_dictionary["powerOffDrag"][0], + powerOnDrag=self.dispersion_dictionary["powerOnDrag"][0], + ) + self.rocket.setRailButtons(distanceToCM=[0.2, -0.5]) + except: + raise TypeError( + "Cannot define basic Rocket and add rail buttons. Missing required parameters in dictionary" + ) + if self.flight is None: + try: + self.flight = Flight( + rocket=self.rocket, + environment=self.environment, + inclination=self.dispersion_dictionary["inclination"][0], + heading=self.dispersion_dictionary["heading"][0], + ) + except: + raise TypeError( + "Cannot define basic Flight. Missing required parameters in dictionary" + ) + return None - def yield_flight_setting( - self, distributionFunc, analysis_parameters, number_of_simulations + def __yield_flight_setting( + self, distribution_func, analysis_parameters, number_of_simulations ): + """Yields a flight setting for the simulation + + Parameters + ---------- + distribution_func : np.random distribution function + The function that will be used to generate the random values. + analysis_parameters : dict + The dictionary with the parameters to be analyzed. This includes the + mean and standard deviation of the parameters. + number_of_simulations : int + Number of simulations desired, must be non negative. + This is needed when running a new simulation. Default is zero. + + Yields + ------ + flight_setting: dict + A dictionary with the flight setting for one simulation. - """Yields a flight setting for the simulation""" + """ i = 0 while i < number_of_simulations: @@ -425,7 +1007,12 @@ def yield_flight_setting( flight_setting = {} for parameter_key, parameter_value in analysis_parameters.items(): if type(parameter_value) is tuple: - flight_setting[parameter_key] = distributionFunc(*parameter_value) + flight_setting[parameter_key] = distribution_func(*parameter_value) + elif isinstance(parameter_value, Function): + flight_setting[parameter_key] = distribution_func(*parameter_value) + elif isinstance(parameter_value, types.FunctionType): + # Deal with parachute triggers functions + flight_setting[parameter_key] = parameter_value else: # shuffles list and gets first item shuffle(parameter_value) @@ -436,87 +1023,90 @@ def yield_flight_setting( # Yield a flight setting yield flight_setting - # TODO: Rework post process Flight method making it possible (and optmized) to - # chose what is going to be exported - def export_flight_data( + def __check_export_list(self, export_list): + """Check if export list is valid or if it is None. In case it is + None, export all possible attributes. + + Parameters + ---------- + export_list : list + List of strings with the names of the attributes to be exported + + Returns + ------- + export_list + """ + + if export_list: + for attr in export_list: + if not isinstance(attr, str): + raise TypeError("Variables must be strings.") + + # Checks if attribute is not valid + if attr not in self.export_list: + raise ValueError( + "Attribute can not be exported. Check export_list." + ) + else: + export_list = self.exportable_list + + return export_list + + def __export_flight_data( self, flight_setting, - flight_data, + flight, exec_time, dispersion_input_file, dispersion_output_file, ): + """Saves flight results in a .txt - """Saves flight results in a .txt""" - - # Generate flight results - flight_result = { - "outOfRailTime": flight_data.outOfRailTime, - "outOfRailVelocity": flight_data.outOfRailVelocity, - "apogeeTime": flight_data.apogeeTime, - "apogeeAltitude": flight_data.apogee - flight_data.env.elevation, - "apogeeX": flight_data.apogeeX, - "apogeeY": flight_data.apogeeY, - "impactTime": flight_data.tFinal, - "impactX": flight_data.xImpact, - "impactY": flight_data.yImpact, - "impactVelocity": flight_data.impactVelocity, - "initialStaticMargin": flight_data.rocket.staticMargin(0), - "outOfRailStaticMargin": flight_data.rocket.staticMargin( - flight_data.outOfRailTime - ), - "finalStaticMargin": flight_data.rocket.staticMargin( - flight_data.rocket.motor.burnOutTime - ), - "numberOfEvents": len(flight_data.parachuteEvents), - "drogueTriggerTime": [], - "drogueInflatedTime": [], - "drogueInflatedVelocity": [], - "executionTime": exec_time, - "lateralWind": flight_data.lateralSurfaceWind, - "frontalWind": flight_data.frontalSurfaceWind, - } + Parameters + ---------- + flight_setting : dict + The flight setting used in the simulation. + flight : Flight + The flight object. + exec_time : float + The execution time of the simulation. + dispersion_input_file : str + The name of the file containing all the inputs for the simulation. + dispersion_output_file : str + The name of the file containing all the outputs for the simulation. - # Calculate maximum reached velocity - sol = np.array(flight_data.solution) - flight_data.vx = Function( - sol[:, [0, 4]], - "Time (s)", - "Vx (m/s)", - "linear", - extrapolation="natural", - ) - flight_data.vy = Function( - sol[:, [0, 5]], - "Time (s)", - "Vy (m/s)", - "linear", - extrapolation="natural", - ) - flight_data.vz = Function( - sol[:, [0, 6]], - "Time (s)", - "Vz (m/s)", - "linear", - extrapolation="natural", - ) - flight_data.v = ( - flight_data.vx**2 + flight_data.vy**2 + flight_data.vz**2 - ) ** 0.5 - flight_data.maxVel = np.amax(flight_data.v.source[:, 1]) - flight_result["maxVelocity"] = flight_data.maxVel + Returns + ------- + None + """ + # TODO: This method is called at every loop of the dispersion + # so all the for loops are slowing down de dispersion + # find a more efficient way to save attributes + + # First, capture the flight data that are saved in the flight object + attributes_list = list(set(dir(flight)).intersection(self.export_list)) + flight_result = {} + for var in attributes_list: + flight_result[str(var)] = getattr(flight, var) + + # Second, capture data that needs to be calculated + for var in list(set(self.export_list) - set(attributes_list)): + if var == "executionTime": + flight_result[str(var)] = exec_time + elif var == "numberOfEvents": + flight_result[str(var)] = len(flight.parachuteEvents) + else: + raise ValueError(f"Variable {var} could not be found.") # Take care of parachute results - for trigger_time, parachute in flight_data.parachuteEvents: + for trigger_time, parachute in flight.parachuteEvents: flight_result[parachute.name + "_triggerTime"] = trigger_time flight_result[parachute.name + "_inflatedTime"] = ( trigger_time + parachute.lag ) - flight_result[parachute.name + "_inflatedVelocity"] = flight_data.v( + flight_result[parachute.name + "_inflatedVelocity"] = flight.speed( trigger_time + parachute.lag ) - else: - flight_result["parachuteInfo"] = "No Parachute Events" # Write flight setting and results to file flight_setting.pop("thrust", None) @@ -525,687 +1115,470 @@ def export_flight_data( return None - def export_flight_error(self, flight_setting, dispersion_error_file): + def __export_flight_data_error(setting, flight_setting, dispersion_error_file): + """Saves flight error in a .txt - """Saves flight error in a .txt""" + Parameters + ---------- + setting : dict + The flight setting used in the simulation. + dispersion_error_file : str + The name of the file containing all the errors for the simulation. + + Returns + ------- + None + """ dispersion_error_file.write(str(flight_setting) + "\n") return None - def runDispersion( + def run_dispersion( self, number_of_simulations, - dispersionDict, - environment, + dispersion_dictionary, + environment=None, flight=None, motor=None, rocket=None, - distributionType="normal", - image=None, - realLandingPoint=None, + distribution_type="normal", + export_list=None, + append=False, ): + """Runs the dispersion simulation and saves all data. For the simulation to be run + all classes must be defined. This can happen either trough the dispersion_dictionary + or by inputing objects - """Runs the given number of simulations and saves the data""" - - self.number_of_simulations = number_of_simulations - self.dispersionDict = dispersionDict - self.environment = environment - self.flight = flight - if flight: - self.motor = flight.rocket.motor if not motor else motor - self.rocket = flight.rocket if not rocket else rocket - self.motor = motor if motor else self.motor - self.rocket = rocket if rocket else self.rocket - self.distributionType = distributionType - self.image = image - self.realLandingPoint = realLandingPoint - - # Creates copy of dispersionDict that will be altered - modified_dispersion_dict = {i: j for i, j in dispersionDict.items()} - - analysis_parameters = self.processDispersionDict(modified_dispersion_dict) - - self.distribuitionFunc = self.setDistributionFunc(distributionType) - # Basic analysis info - - # Create data files for inputs, outputs and error logging - dispersion_error_file = open(str(self.filename) + ".disp_errors.txt", "w") - dispersion_input_file = open(str(self.filename) + ".disp_inputs.txt", "w") - dispersion_output_file = open(str(self.filename) + ".disp_outputs.txt", "w") - - # # Initialize Environment - # customAtmosphere = False - # if not self.environment: - # self.environment = Environment( - # railLength=0, - # ) - # if "envAtmosphericType" in dispersionDict: - # if dispersionDict["envAtmosphericType"] == "CustomAtmosphere": - # customAtmosphere = True - # self.environment.setDate(datetime(*dispersionDict["date"][0])) - # self.environment.setAtmosphericModel( - # type=dispersionDict["envAtmosphericType"], - # file=dispersionDict["envAtmosphericFile"] - # if "envAtmosphericFile" in dispersionDict - # else None, - # dictionary=dispersionDict["envAtmosphericDictionary"] - # if "envAtmosphericDictionary" in dispersionDict - # else None, - # ) - - # Initialize counter and timer - i = 0 - - initial_wall_time = time() - initial_cpu_time = process_time() - - # Iterate over flight settings - out = display("Starting", display_id=True) - for setting in self.yield_flight_setting( - self.distribuitionFunc, analysis_parameters, self.number_of_simulations - ): - start_time = process_time() - i += 1 - - # Creates an of environment - envDispersion = self.environment - - # Apply environment parameters variations on each iteration if possible - envDispersion.railLength = setting["railLength"] - envDispersion.gravity = setting["gravity"] - envDispersion.date = setting["date"] - envDispersion.latitude = setting["latitude"] - envDispersion.longitude = setting["longitude"] - envDispersion.elevation = setting["elevation"] - envDispersion.selectEnsembleMember(setting["ensembleMember"]) - - # Creates copy of motor - motorDispersion = self.motor - - # Apply motor parameters variations on each iteration if possible - # TODO: add hybrid motor option - motorDispersion = SolidMotor( - thrustSource=setting["thrust"], - burnOut=setting["burnOutTime"], - grainNumber=setting["grainNumber"], - grainDensity=setting["grainDensity"], - grainOuterRadius=setting["grainOuterRadius"], - grainInitialInnerRadius=setting["grainInitialInnerRadius"], - grainInitialHeight=setting["grainInitialHeight"], - grainSeparation=setting["grainSeparation"], - nozzleRadius=setting["nozzleRadius"], - throatRadius=setting["throatRadius"], - reshapeThrustCurve=(setting["burnOutTime"], setting["totalImpulse"]), - ) - - # Creates copy of rocket - rocketDispersion = self.rocket - - # Apply rocket parameters variations on each iteration if possible - rocketDispersion = Rocket( - motor=motorDispersion, - mass=setting["mass"], - inertiaI=setting["inertiaI"], - inertiaZ=setting["inertiaZ"], - radius=setting["radius"], - distanceRocketNozzle=setting["distanceRocketNozzle"], - distanceRocketPropellant=setting["distanceRocketPropellant"], - powerOffDrag=setting["powerOffDrag"], - powerOnDrag=setting["powerOnDrag"], - ) - - # Add rocket nose, fins and tail - rocketDispersion.addNose( - length=setting["noseLength"], - kind=setting["noseKind"], - distanceToCM=setting["noseDistanceToCM"], - ) - rocketDispersion.addFins( - n=setting["numberOfFins"], - rootChord=setting["rootChord"], - tipChord=setting["tipChord"], - span=setting["span"], - distanceToCM=setting["distanceToCM"], - radius=setting["radius"], - airfoil=setting["airfoil"], - ) - if not "noTail" in setting: - rocketDispersion.addTail( - topRadius=setting["topRadius"], - bottomRadius=setting["bottomRadius"], - length=setting["length"], - distanceToCM=setting["distanceToCM"], - ) - - # Add parachutes - for num, name in enumerate(self.parachute_names): - rocketDispersion.addParachute( - name=name, - CdS=setting["parachute_" + name + "_CdS"], - trigger=setting["parachute_" + name + "_trigger"], - samplingRate=setting["parachute_" + name + "_samplingRate"], - lag=setting["parachute_" + name + "_lag"], - noise=( - setting["parachute_" + name + "_noise_mean"], - setting["parachute_" + name + "_noise_std"], - setting["parachute_" + name + "_noise_corr"], - ), - ) - - rocketDispersion.setRailButtons( - distanceToCM=[ - setting["positionFirstRailButton"], - setting["positionSecondRailButton"], - ], - angularPosition=setting["railButtonAngularPosition"], - ) - - # Run trajectory simulation - try: - TestFlight = Flight( - rocket=rocketDispersion, - environment=envDispersion, - inclination=setting["inclination"], - heading=setting["heading"], - # initialSolution=setting["initialSolution"] if "initialSolution" in setting else self.flight.initialSolution, - terminateOnApogee=setting["terminateOnApogee"], - maxTime=setting["maxTime"], - maxTimeStep=setting["maxTimeStep"], - minTimeStep=setting["minTimeStep"], - rtol=setting["rtol"], - atol=setting["atol"], - timeOvershoot=setting["timeOvershoot"], - verbose=setting["verbose"], - ) - - self.export_flight_data( - setting, - TestFlight, - process_time() - start_time, - dispersion_input_file, - dispersion_output_file, - ) - except Exception as E: - print(E) - print(traceback.format_exc()) - self.export_flight_error(setting, dispersion_error_file) - - # Register time - out.update( - f"Curent iteration: {i:06d} | Average Time per Iteration: {(process_time() - initial_cpu_time)/i:2.6f} s | Estimated time left: {int((number_of_simulations - i)*((process_time() - initial_cpu_time)/i))} s" - ) - - # Done - - ## Print and save total time - final_string = f"Completed {i} iterations successfully. Total CPU time: {process_time() - initial_cpu_time} s. Total wall time: {time() - initial_wall_time} s" - out.update(final_string) - dispersion_input_file.write(final_string + "\n") - dispersion_output_file.write(final_string + "\n") - dispersion_error_file.write(final_string + "\n") - - ## Close files - dispersion_input_file.close() - dispersion_output_file.close() - dispersion_error_file.close() - - return None - - def importResults(self, dispersion_output_file): - - """Import dispersion results from .txt file""" - - # Initialize variable to store all results - dispersion_general_results = [] - - dispersion_results = { - "outOfRailTime": [], - "outOfRailVelocity": [], - "apogeeTime": [], - "apogeeAltitude": [], - "apogeeX": [], - "apogeeY": [], - "impactTime": [], - "impactX": [], - "impactY": [], - "impactVelocity": [], - "initialStaticMargin": [], - "outOfRailStaticMargin": [], - "finalStaticMargin": [], - "numberOfEvents": [], - "maxVelocity": [], - "drogueTriggerTime": [], - "drogueInflatedTime": [], - "drogueInflatedVelocity": [], - "executionTime": [], - "railDepartureAngleOfAttack": [], - "lateralWind": [], - "frontalWind": [], - } - - # Get all dispersion results - # Get file - dispersion_output_file = open(dispersion_output_file, "r+") - - # Read each line of the file and convert to dict - for line in dispersion_output_file: - # Skip comments lines - if line[0] != "{": - continue - # Eval results and store them - flight_result = eval(line) - dispersion_general_results.append(flight_result) - for parameter_key, parameter_value in flight_result.items(): - dispersion_results[parameter_key].append(parameter_value) - - # Close data file - dispersion_output_file.close() - - # Number of flights simulated - self.N = len(dispersion_general_results) - - return dispersion_results - - def meanOutOfRailTime(self, dispersion_results): - print( - f'Out of Rail Time - Mean Value: {np.mean(dispersion_results["outOfRailTime"]):0.3f} s' - ) - print( - f'Out of Rail Time - Standard Deviation: {np.std(dispersion_results["outOfRailTime"]):0.3f} s' - ) - - return None - - def plotOutOfRailTime(self, dispersion_results): - - self.meanOutOfRailTime(dispersion_results) - - plt.figure() - plt.hist(dispersion_results["outOfRailTime"], bins=int(self.N**0.5)) - plt.title("Out of Rail Time") - plt.xlabel("Time (s)") - plt.ylabel("Number of Occurences") - plt.show() - - return None - - def meanOutOfRailVelocity(self, dispersion_results): - print( - f'Out of Rail Velocity - Mean Value: {np.mean(dispersion_results["outOfRailVelocity"]):0.3f} m/s' - ) - print( - f'Out of Rail Velocity - Standard Deviation: {np.std(dispersion_results["outOfRailVelocity"]):0.3f} m/s' - ) - - return None - - def plotOutOfRailVelocity(self, dispersion_results): - - self.meanOutOfRailVelocity(dispersion_results) - - plt.figure() - plt.hist(dispersion_results["outOfRailVelocity"], bins=int(self.N**0.5)) - plt.title("Out of Rail Velocity") - plt.xlabel("Velocity (m/s)") - plt.ylabel("Number of Occurences") - plt.show() - - return None - - def meanApogeeTime(self, dispersion_results): - print( - f'Impact Time - Mean Value: {np.mean(dispersion_results["impactTime"]):0.3f} s' - ) - print( - f'Impact Time - Standard Deviation: {np.std(dispersion_results["impactTime"]):0.3f} s' - ) - - return None - - def plotApogeeTime(self, dispersion_results): - - self.meanApogeeTime(dispersion_results) - - plt.figure() - plt.hist(dispersion_results["impactTime"], bins=int(self.N**0.5)) - plt.title("Impact Time") - plt.xlabel("Time (s)") - plt.ylabel("Number of Occurences") - plt.show() - - return None - - def meanApogeeAltitude(self, dispersion_results): - print( - f'Apogee Altitude - Mean Value: {np.mean(dispersion_results["apogeeAltitude"]):0.3f} m' - ) - print( - f'Apogee Altitude - Standard Deviation: {np.std(dispersion_results["apogeeAltitude"]):0.3f} m' - ) - - return None - - def plotApogeeAltitude(self, dispersion_results): - - self.meanApogeeAltitude(dispersion_results) - - plt.figure() - plt.hist(dispersion_results["apogeeAltitude"], bins=int(self.N**0.5)) - plt.title("Apogee Altitude") - plt.xlabel("Altitude (m)") - plt.ylabel("Number of Occurences") - plt.show() - - return None - - def meanApogeeXPosition(self, dispersion_results): - print( - f'Apogee X Position - Mean Value: {np.mean(dispersion_results["apogeeX"]):0.3f} m' - ) - print( - f'Apogee X Position - Standard Deviation: {np.std(dispersion_results["apogeeX"]):0.3f} m' - ) - - return None - - def plotApogeeXPosition(self, dispersion_results): - - self.meanApogeeAltitude(dispersion_results) - - plt.figure() - plt.hist(dispersion_results["apogeeX"], bins=int(self.N**0.5)) - plt.title("Apogee X Position") - plt.xlabel("Apogee X Position (m)") - plt.ylabel("Number of Occurences") - plt.show() - - return None - - def meanApogeeYPosition(self, dispersion_results): - print( - f'Apogee Y Position - Mean Value: {np.mean(dispersion_results["apogeeY"]):0.3f} m' - ) - print( - f'Apogee Y Position - Standard Deviation: {np.std(dispersion_results["apogeeY"]):0.3f} m' - ) - - return None - - def plotApogeeYPosition(self, dispersion_results): - - self.meanApogeeAltitude(dispersion_results) - - plt.figure() - plt.hist(dispersion_results["apogeeY"], bins=int(self.N**0.5)) - plt.title("Apogee Y Position") - plt.xlabel("Apogee Y Position (m)") - plt.ylabel("Number of Occurences") - plt.show() - - return None - - def meanImpactTime(self, dispersion_results): - print( - f'Impact Time - Mean Value: {np.mean(dispersion_results["impactTime"]):0.3f} s' - ) - print( - f'Impact Time - Standard Deviation: {np.std(dispersion_results["impactTime"]):0.3f} s' - ) - - return None - - def plotImpactTime(self, dispersion_results): - - self.meanImpactTime(dispersion_results) - - plt.figure() - plt.hist(dispersion_results["impactTime"], bins=int(self.N**0.5)) - plt.title("Impact Time") - plt.xlabel("Time (s)") - plt.ylabel("Number of Occurences") - plt.show() - - return None - - def meanImpactXPosition(self, dispersion_results): - print( - f'Impact X Position - Mean Value: {np.mean(dispersion_results["impactX"]):0.3f} m' - ) - print( - f'Impact X Position - Standard Deviation: {np.std(dispersion_results["impactX"]):0.3f} m' - ) - - return None - - def plotImpactXPosition(self, dispersion_results): - - self.meanImpactXPosition(dispersion_results) - - plt.figure() - plt.hist(dispersion_results["impactX"], bins=int(self.N**0.5)) - plt.title("Impact X Position") - plt.xlabel("Impact X Position (m)") - plt.ylabel("Number of Occurences") - plt.show() - - return None - - def meanImpactYPosition(self, dispersion_results): - print( - f'Impact Y Position - Mean Value: {np.mean(dispersion_results["impactY"]):0.3f} m' - ) - print( - f'Impact Y Position - Standard Deviation: {np.std(dispersion_results["impactY"]):0.3f} m' - ) - - return None + Parameters + ---------- + number_of_simulations : int + Number of simulations to be run, must be non negative. + dispersion_dictionary : dict + The dictionary with the parameters to be analyzed. The keys must be the + names of the attributes that will be used in the dispersion simulation. + The values can either be a tuple, containing the nominal values of that + parameter and its standard deviation, a list, containing the possible + values to be randomly chosen in each simulation, or a single value (int + or float), being the standard deviation of that parameter. See example + for further explanations. + environment : Environment, optional + Environment object that will be used in the simulations. Default is None. + If none, environment must be defined via passing its attributes in the + dispersion_dictionary. Arguments related to environment will only vary + according to the distribution method if the standard deviation for the + desired attributes are on the dispersion_dictionary. + flight : Flight, optional + Flight object that will be used in the simulations. Default is None. + If none, Flight must be defined via passing its attributes in the + dispersion_dictionary. Arguments related to Flight will only vary + according to the distribution method if the standard deviation for the + desired attributes are on the dispersion_dictionary. + motor : Motor, optional + Motor object that will be used in the simulations. Default is None. + If none, Motor must be defined via passing its attributes in the + dispersion_dictionary. Arguments related to Motor will only vary + according to the distribution method if the standard deviation for the + desired attributes are on the dispersion_dictionary. + rocket : Rocket, optional + Rocket object that will be used in the simulations. Default is None. + If none, Rocket must be defined via passing its attributes in the + dispersion_dictionary. Arguments related to Rocket will only vary + according to the distribution method if the standard deviation for the + desired attributes are on the dispersion_dictionary. + distribution_type : str, optional + The probability distribution function to be used in the analysis, + by default "normal". Options are any numpy.ramdom distributions + export_list : list, optional + A list containing the name of the attributes to be saved on the dispersion + outputs file. See Examples for all possible attribues + append : bool, optional + If True, the results will be appended to the existing files. If False, + the files will be overwritten. By default False. - def plotImpactYPosition(self, dispersion_results): + Returns + ------- + None - self.meanImpactYPosition(dispersion_results) + Examples + -------- - plt.figure() - plt.hist(dispersion_results["impactY"], bins=int(self.N**0.5)) - plt.title("Impact Y Position") - plt.xlabel("Impact Y Position (m)") - plt.ylabel("Number of Occurences") - plt.show() + TODO: add list of all possible attributes in dispersion_dictionaries and + all possible attributes in export_list - return None + """ - def meanImpactVelocity(self, dispersion_results): - print( - f'Impact Velocity - Mean Value: {np.mean(dispersion_results["impactVelocity"]):0.3f} m/s' - ) - print( - f'Impact Velocity - Standard Deviation: {np.std(dispersion_results["impactVelocity"]):0.3f} m/s' - ) + # Saving the arguments as attributes + self.number_of_simulations = number_of_simulations + self.dispersion_dictionary = dispersion_dictionary + if flight: # In case a flight object is passed + self.environment = flight.env + self.motor = flight.rocket.motor + self.rocket = flight.rocket + self.flight = flight + if rocket: + self.rocket = rocket + self.motor = rocket.motor + if motor: + self.motor = motor + if environment: + self.environment = environment + + self.distribution_type = distribution_type + + # Check if there's enough object to start a flight: + ## Raise an error in case of any troubles + self.__create_initial_objects() + + # Creates copy of dispersion_dictionary that will be altered + modified_dispersion_dict = {i: j for i, j in dispersion_dictionary.items()} + + analysis_parameters = self.__process_dispersion_dict(modified_dispersion_dict) + + # TODO: This should be more flexible, allow different distributions for different parameters + self.distributionFunc = self.__set_distribution_function(self.distribution_type) - return None + # Create data files for inputs, outputs and error logging + open_mode = "a" if append else "w" + dispersion_error_file = open(f"{self.filename}.disp_errors.txt", open_mode) + dispersion_input_file = open(f"{self.filename}.disp_inputs.txt", open_mode) + dispersion_output_file = open(f"{self.filename}.disp_outputs.txt", open_mode) - def plotImpactVelocity(self, dispersion_results): + # Checks export_list + self.export_list = self.__check_export_list(export_list) - self.meanImpactVelocity(dispersion_results) + # Creates a copy of the environment + env_dispersion = self.environment - plt.figure() - plt.hist(dispersion_results["impactVelocity"], bins=int(self.N**0.5)) - plt.title("Impact Velocity") - plt.xlim(-35, 0) - plt.xlabel("Velocity (m/s)") - plt.ylabel("Number of Occurences") - plt.show() + # Creates copy of motor + motor_dispersion = self.motor - return None + # Creates copy of rocket + rocket_dispersion = self.rocket - def meanStaticMargin(self, dispersion_results): - print( - f'Initial Static Margin - Mean Value: {np.mean(dispersion_results["initialStaticMargin"]):0.3f} c' - ) - print( - f'Initial Static Margin - Standard Deviation: {np.std(dispersion_results["initialStaticMargin"]):0.3f} c' - ) + # Initialize counter and timer + i = 0 + initial_wall_time = time() + initial_cpu_time = process_time() - print( - f'Out of Rail Static Margin - Mean Value: {np.mean(dispersion_results["outOfRailStaticMargin"]):0.3f} c' - ) - print( - f'Out of Rail Static Margin - Standard Deviation: {np.std(dispersion_results["outOfRailStaticMargin"]):0.3f} c' - ) + # Iterate over flight settings, start the flight simulations + out = display("Starting", display_id=True) + for setting in self.__yield_flight_setting( + self.distributionFunc, analysis_parameters, self.number_of_simulations + ): + self.start_time = process_time() + i += 1 - print( - f'Final Static Margin - Mean Value: {np.mean(dispersion_results["finalStaticMargin"]):0.3f} c' - ) - print( - f'Final Static Margin - Standard Deviation: {np.std(dispersion_results["finalStaticMargin"]):0.3f} c' - ) + # Apply environment parameters variations on each iteration if possible + env_dispersion.railLength = setting["railLength"] + env_dispersion.gravity = setting["gravity"] + env_dispersion.date = setting["date"] + env_dispersion.latitude = setting["latitude"] + env_dispersion.longitude = setting["longitude"] + env_dispersion.elevation = setting["elevation"] + if env_dispersion.atmosphericModelType in ["Ensemble", "Reanalysis"]: + env_dispersion.selectEnsembleMember(setting["ensembleMember"]) - return None + # Apply motor parameters variations on each iteration if possible + # TODO: add hybrid and liquid motor option + motor_dispersion = SolidMotor( + thrustSource=setting["thrust"], + burnOut=setting["burnOutTime"], + grainNumber=setting["grainNumber"], + grainDensity=setting["grainDensity"], + grainOuterRadius=setting["grainOuterRadius"], + grainInitialInnerRadius=setting["grainInitialInnerRadius"], + grainInitialHeight=setting["grainInitialHeight"], + grainSeparation=setting["grainSeparation"], + nozzleRadius=setting["nozzleRadius"], + throatRadius=setting["throatRadius"], + reshapeThrustCurve=(setting["burnOutTime"], setting["totalImpulse"]), + ) - def plotStaticMargin(self, dispersion_results): + # Apply rocket parameters variations on each iteration if possible + rocket_dispersion = Rocket( + motor=motor_dispersion, + mass=setting["mass"], + inertiaI=setting["inertiaI"], + inertiaZ=setting["inertiaZ"], + radius=setting["radius"], + distanceRocketNozzle=setting["distanceRocketNozzle"], + distanceRocketPropellant=setting["distanceRocketPropellant"], + powerOffDrag=setting["powerOffDrag"], + powerOnDrag=setting["powerOnDrag"], + ) - self.meanStaticMargin(dispersion_results) + # Add rocket nose, fins and tail + # Nose + for nose in self.nose_names: + rocket_dispersion.addNose( + length=setting[f"nose_{nose}_length"], + kind=setting[f"nose_{nose}_kind"], + distanceToCM=setting[f"nose_{nose}_distanceToCM"], + name=nose, + ) - plt.figure() - plt.hist( - dispersion_results["initialStaticMargin"], - label="Initial", - bins=int(self.N**0.5), - ) - plt.hist( - dispersion_results["outOfRailStaticMargin"], - label="Out of Rail", - bins=int(self.N**0.5), - ) - plt.hist( - dispersion_results["finalStaticMargin"], - label="Final", - bins=int(self.N**0.5), - ) - plt.legend() - plt.title("Static Margin") - plt.xlabel("Static Margin (c)") - plt.ylabel("Number of Occurences") - plt.show() + # Fins + for finSet in self.finSet_names: + # TODO: Allow elliptical fins + rocket_dispersion.addTrapezoidalFins( + n=setting[f"finSet_{finSet}_n"], + rootChord=setting[f"finSet_{finSet}_rootChord"], + tipChord=setting[f"finSet_{finSet}_tipChord"], + span=setting[f"finSet_{finSet}_span"], + distanceToCM=setting[f"finSet_{finSet}_distanceToCM"], + airfoil=setting[f"finSet_{finSet}_airfoil"], + name=finSet, + ) - return None + # Tail + for tail in self.tail_names: + rocket_dispersion.addTail( + topRadius=setting[f"tail_{tail}_topRadius"], + bottomRadius=setting[f"tail_{tail}_bottomRadius"], + length=setting[f"tail_{tail}_length"], + distanceToCM=setting[f"tail_{tail}_distanceToCM"], + radius=None, + name="Tail", + ) - def meanMaximumVelocity(self, dispersion_results): - print( - f'Maximum Velocity - Mean Value: {np.mean(dispersion_results["maxVelocity"]):0.3f} m/s' - ) - print( - f'Maximum Velocity - Standard Deviation: {np.std(dispersion_results["maxVelocity"]):0.3f} m/s' - ) + # Add parachutes + rocket_dispersion.parachutes = [] # Remove existing parachutes + for name in self.parachute_names: + rocket_dispersion.addParachute( + name=name, + CdS=setting["parachute_" + name + "_CdS"], + trigger=setting["parachute_" + name + "_trigger"], + samplingRate=setting["parachute_" + name + "_samplingRate"], + lag=setting["parachute_" + name + "_lag"], + noise=setting["parachute_" + name + "_noise"], + ) - return None + rocket_dispersion.setRailButtons( + distanceToCM=[ + setting["positionFirstRailButton"], + setting["positionSecondRailButton"], + ], + angularPosition=setting["railButtonAngularPosition"], + ) - def plotMaximumVelocity(self, dispersion_results): + # Run trajectory simulation + try: + # TODO: Add initialSolution flight option + dispersion_flight = Flight( + rocket=rocket_dispersion, + environment=env_dispersion, + inclination=setting["inclination"], + heading=setting["heading"], + terminateOnApogee=setting["terminateOnApogee"], + maxTime=setting["maxTime"], + maxTimeStep=setting["maxTimeStep"], + minTimeStep=setting["minTimeStep"], + rtol=setting["rtol"], + atol=setting["atol"], + timeOvershoot=setting["timeOvershoot"], + verbose=setting["verbose"], + ) - self.meanMaximumVelocity(dispersion_results) + self.__export_flight_data( + flight_setting=setting, + flight=dispersion_flight, + exec_time=process_time() - self.start_time, + dispersion_input_file=dispersion_input_file, + dispersion_output_file=dispersion_output_file, + ) + except Exception as E: + print(E) + print(traceback.format_exc()) + self.__export_flight_data_error(setting, dispersion_error_file) - plt.figure() - plt.hist(dispersion_results["maxVelocity"], bins=int(self.N**0.5)) - plt.title("Maximum Velocity") - plt.xlabel("Velocity (m/s)") - plt.ylabel("Number of Occurences") - plt.show() + # Register time + out.update( + f"Current iteration: {i:06d} | Average Time per Iteration: " + f"{(process_time() - initial_cpu_time)/i:2.6f} s | Estimated time" + f" left: {int((number_of_simulations - i)*((process_time() - initial_cpu_time)/i))} s" + ) - return None + # Clean the house once all the simulations were already done - def meanNumberOfParachuteEvents(self, dispersion_results): - print( - f'Number of Parachute Events - Mean Value: {np.mean(dispersion_results["numberOfEvents"]):0.3f} s' - ) - print( - f'Number of Parachute Events - Standard Deviation: {np.std(dispersion_results["numberOfEvents"]):0.3f} s' + ## Print and save total time + final_string = ( + f"Completed {i} iterations successfully. Total CPU time: " + f"{process_time() - initial_cpu_time} s. Total wall time: " + f"{time() - initial_wall_time} s" ) + out.update(final_string) + dispersion_input_file.write(final_string + "\n") + dispersion_output_file.write(final_string + "\n") + dispersion_error_file.write(final_string + "\n") + + ## Close files + dispersion_input_file.close() + dispersion_output_file.close() + dispersion_error_file.close() return None - def plotNumberOfParachuteEvents(self, dispersion_results): + def import_results(self, variables=None): + """Import dispersion results from .txt file and save it into a dictionary. - self.meanNumberOfParachuteEvents(dispersion_results) + Parameters + ---------- + variables : list of str, optional + List of variables to be imported. If None, all variables will be imported. - plt.figure() - plt.hist(dispersion_results["numberOfEvents"]) - plt.title("Parachute Events") - plt.xlabel("Number of Parachute Events") - plt.ylabel("Number of Occurences") - plt.show() + Returns + ------- + None + """ + # Initialize variable to store all results + dispersion_results = {} - return None + # Get all dispersion results + # Open the file + file = open(self.filename.split(".")[0] + ".disp_outputs.txt", "r+") - def meanDrogueTriggerTime(self, dispersion_results): - print( - f'Drogue Trigger Time - Mean Value: {np.mean(dispersion_results["drogueTriggerTime"]):0.3f} s' - ) + # Read each line of the file and convert to dict + for line in file: + # Skip comments lines + if line[0] != "{": + continue + # Evaluate results and store them + flight_result = eval(line) + # Append to the list + for parameter_key, parameter_value in flight_result.items(): + if parameter_key not in dispersion_results.keys(): + # Create a new list to store the parameter + dispersion_results[parameter_key] = [parameter_value] + else: + # Append the parameter value to the list + dispersion_results[parameter_key].append(parameter_value) + + # Close data file + file.close() + + # Calculate the number of flights simulated + len_dict = {key: len(value) for key, value in dispersion_results.items()} + if min(len_dict.values()) - max(len_dict.values()) > 1: + print( + "Warning: The number of simulations imported from the file is not " + "the same for all parameters. The number of simulations will be " + "set to the minimum number of simulations found." + ) + self.num_of_loaded_sims = min(len_dict.values()) + + # Print the number of flights simulated print( - f'Drogue Trigger Time - Standard Deviation: {np.std(dispersion_results["drogueTriggerTime"]):0.3f} s' + f"A total of {self.num_of_loaded_sims} simulations were loaded from" + f" the following file: {self.filename.split('.')[0] + '.disp_outputs.txt'}" ) - return None + # Save the results as an attribute of the class + self.dispersion_results = dispersion_results - def plotDrogueTriggerTime(self, dispersion_results): + # Process the results and save them as attributes of the class + self.__process_results(variables=variables) - self.meanDrogueTriggerTime(dispersion_results) + return None - plt.figure() - plt.hist(dispersion_results["drogueTriggerTime"], bins=int(self.N**0.5)) - plt.title("Drogue Trigger Time") - plt.xlabel("Time (s)") - plt.ylabel("Number of Occurences") - plt.show() + # Start the processing analysis - return None + def __process_results(self, variables=None): + """Save the mean and standard deviation of each parameter available + in the results dictionary. Create class attributes for each parameter. - def meanDrogueFullyInflatedTime(self, dispersion_results): - print( - f'Drogue Fully Inflated Time - Mean Value: {np.mean(dispersion_results["drogueInflatedTime"]):0.3f} s' - ) - print( - f'Drogue Fully Inflated Time - Standard Deviation: {np.std(dispersion_results["drogueInflatedTime"]):0.3f} s' - ) + Parameters + ---------- + variables : list, optional + List of variables to be processed. If None, all variables will be + processed. The default is None. Example: ['outOfRailTime', 'apogeeTime'] + Returns + ------- + None + """ + if isinstance(variables, list): + for result in variables: + mean = np.mean(self.dispersion_results[result]) + stdev = np.std(self.dispersion_results[result]) + setattr(self, str(result), (mean, stdev)) + else: + for result in self.dispersion_results.keys(): + mean = np.mean(self.dispersion_results[result]) + stdev = np.std(self.dispersion_results[result]) + setattr(self, str(result), (mean, stdev)) return None - def plotDrogueFullyInflatedTime(self, dispersion_results): + # TODO: print as a table instead of prints + def print_results(self, variables=None): + """Print the mean and standard deviation of each parameter in the results + dictionary or of the variables passed as argument. - self.meanDrogueFullyInflatedTime(dispersion_results) + Parameters + ---------- + variables : list, optional + List of variables to be processed. If None, all variables will be + processed. The default is None. Example: ['outOfRailTime', 'apogee'] - plt.figure() - plt.hist(dispersion_results["drogueInflatedTime"], bins=int(self.N**0.5)) - plt.title("Drogue Fully Inflated Time") - plt.xlabel("Time (s)") - plt.ylabel("Number of Occurences") - plt.show() + Returns + ------- + None - return None + Raises + ------ + TypeError + If the variable passed as argument is not a string. + """ + # Check if the variables argument is a list, if not, use all variables + if not isinstance(variables, list): + variables = self.dispersion_results.keys() - def meanDrogueFullyVelocity(self, dispersion_results): - print( - f'Drogue Parachute Fully Inflated Velocity - Mean Value: {np.mean(dispersion_results["drogueInflatedVelocity"]):0.3f} m/s' - ) - print( - f'Drogue Parachute Fully Inflated Velocity - Standard Deviation: {np.std(dispersion_results["drogueInflatedVelocity"]):0.3f} m/s' - ) + # Check if the variables are strings + if not all(isinstance(var, str) for var in variables): + raise TypeError("The list of variables must be a list of strings.") + + for var in variables: + tp = getattr(self, var) # Get the tuple with the mean and stdev + print("{}: \u03BC = {:.3f}, \u03C3 = {:.3f}".format(var, tp[0], tp[1])) return None - def plotDrogueFullyVelocity(self, dispersion_results): + def plot_results(self, variables=None): + """Plot the results of the dispersion analysis. - self.meanDrogueFullyVelocity(dispersion_results) + Parameters + ---------- + variables : list, optional + List of variables to be plotted. If None, all variables will be + plotted. The default is None. Example: ['outOfRailTime', 'apogee'] - plt.figure() - plt.hist(dispersion_results["drogueInflatedVelocity"], bins=int(self.N**0.5)) - plt.title("Drogue Parachute Fully Inflated Velocity") - plt.xlabel("Velocity m/s)") - plt.ylabel("Number of Occurences") - plt.show() + Returns + ------- + None + """ + # Check if the variables argument is a list, if not, use all variables + if not isinstance(variables, list): + variables = self.dispersion_results.keys() + + # Check if the variables are strings + if not all(isinstance(var, str) for var in variables): + raise TypeError("The list of variables must be a list of strings.") + + for var in variables: + plt.figure() + plt.hist( + self.dispersion_results[var], + ) + plt.title("Histogram of " + var) + # plt.xlabel("Time (s)") + plt.ylabel("Number of Occurrences") + plt.show() return None - def createEllipses(self, dispersion_results): + # TODO: Create evolution plots to analyze convergence + + def __createEllipses(self, dispersion_results): """A function to create apogee and impact ellipses from the dispersion results. @@ -1213,13 +1586,38 @@ def createEllipses(self, dispersion_results): ---------- dispersion_results : dict A dictionary containing the results of the dispersion analysis. + + Returns + ------- + apogee_ellipse : Ellipse + An ellipse object representing the apogee ellipse. + impact_ellipse : Ellipse + An ellipse object representing the impact ellipse. + apogeeX : np.array + An array containing the x coordinates of the apogee ellipse. + apogeeY : np.array + An array containing the y coordinates of the apogee ellipse. + impactX : np.array + An array containing the x coordinates of the impact ellipse. + impactY : np.array + An array containing the y coordinates of the impact ellipse. """ # Retrieve dispersion data por apogee and impact XY position - apogeeX = np.array(dispersion_results["apogeeX"]) - apogeeY = np.array(dispersion_results["apogeeY"]) - impactX = np.array(dispersion_results["impactX"]) - impactY = np.array(dispersion_results["impactY"]) + try: + apogeeX = np.array(dispersion_results["apogeeX"]) + apogeeY = np.array(dispersion_results["apogeeY"]) + except KeyError: + print("No apogee data found.") + apogeeX = np.array([]) + apogeeY = np.array([]) + try: + impactX = np.array(dispersion_results["xImpact"]) + impactY = np.array(dispersion_results["yImpact"]) + except KeyError: + print("No impact data found.") + impactX = np.array([]) + impactY = np.array([]) # Define function to calculate eigen values def eigsorted(cov): @@ -1264,13 +1662,13 @@ def eigsorted(cov): ) apogeeEll.set_facecolor((0, 1, 0, 0.2)) apogee_ellipses.append(apogeeEll) - return impact_ellipses, apogee_ellipses + return impact_ellipses, apogee_ellipses, apogeeX, apogeeY, impactX, impactY def plotEllipses( self, dispersion_results, image=None, - realLandingPoint=None, + actual_landing_point=None, perimeterSize=3000, xlim=(-3000, 3000), ylim=(-3000, 3000), @@ -1283,22 +1681,35 @@ def plotEllipses( ---------- dispersion_results : dict A dictionary containing the results of the dispersion analysis - image : str + image : str, optional The path to the image to be used as the background - realLandingPoint : tuple, optional - A tuple containing the real landing point of the rocket, by default None + actual_landing_point : tuple, optional + A tuple containing the actual landing point of the rocket, if known. + Useful when comparing the dispersion results with the actual landing. + Must be given in tuple format, such as (lat, lon). By default None. # TODO: Check the order + perimeterSize : int, optional + The size of the perimeter to be plotted. The default is 3000. + xlim : tuple, optional + The limits of the x axis. The default is (-3000, 3000). + ylim : tuple, optional + The limits of the y axis. The default is (-3000, 3000). + + Returns + ------- + None """ # Import background map if image is not None: img = imread(image) - # Retrieve dispersion data por apogee and impact XY position - apogeeX = np.array(dispersion_results["apogeeX"]) - apogeeY = np.array(dispersion_results["apogeeY"]) - impactX = np.array(dispersion_results["impactX"]) - impactY = np.array(dispersion_results["impactY"]) - - impact_ellipses, apogee_ellipses = self.createEllipses(dispersion_results) + ( + impact_ellipses, + apogee_ellipses, + apogeeX, + apogeeY, + impactX, + impactY, + ) = self.__createEllipses(dispersion_results) # Create plot figure plt.figure(num=None, figsize=(8, 6), dpi=150, facecolor="w", edgecolor="k") @@ -1325,10 +1736,10 @@ def plotEllipses( label="Simulated Landing Point", ) # Draw real landing point - if realLandingPoint != None: + if actual_landing_point != None: plt.scatter( - realLandingPoint[0], - realLandingPoint[1], + actual_landing_point[0], + actual_landing_point[1], s=20, marker="X", color="red", @@ -1345,6 +1756,7 @@ def plotEllipses( ax.set_xlabel("East (m)") # Add background image to plot + # TODO: In the future, integrate with other libraries to plot the map (e.g. cartopy, ee, etc.) # You can translate the basemap by changing dx and dy (in meters) dx = 0 dy = 0 @@ -1370,7 +1782,7 @@ def plotEllipses( plt.show() return None - def prepareEllipses(self, ellipses, origin_lat, origin_lon, resolution=100): + def __prepareEllipses(self, ellipses, origin_lat, origin_lon, resolution=100): """Generate a list of latitude and longitude points for each ellipse in ellipses. @@ -1384,6 +1796,12 @@ def prepareEllipses(self, ellipses, origin_lat, origin_lon, resolution=100): Longitude of the origin of the coordinate system. resolution : int, optional Number of points to generate for each ellipse, by default 100 + + Returns + ------- + list + List of lists of tuples containing the latitude and longitude of each + point in each ellipse. """ outputs = [] @@ -1427,7 +1845,6 @@ def prepareEllipses(self, ellipses, origin_lat, origin_lon, resolution=100): def exportEllipsesToKML( self, - dispersion_results, filename, origin_lat, origin_lon, @@ -1454,18 +1871,29 @@ def exportEllipsesToKML( Number of points to be used to draw the ellipse. Default is 100. color : String Color of the ellipse. Default is 'ff0000ff', which is red. + + Returns + ------- + None """ - impact_ellipses, apogee_ellipses = self.createEllipses(dispersion_results) + ( + impact_ellipses, + apogee_ellipses, + _, + _, + _, + _, + ) = self.__createEllipses(self.dispersion_results) outputs = [] if type == "all" or type == "impact": - outputs = outputs + self.prepareEllipses( + outputs = outputs + self.__prepareEllipses( impact_ellipses, origin_lat, origin_lon, resolution=resolution ) if type == "all" or type == "apogee": - outputs = outputs + self.prepareEllipses( + outputs = outputs + self.__prepareEllipses( apogee_ellipses, origin_lat, origin_lon, resolution=resolution ) @@ -1507,173 +1935,38 @@ def exportEllipsesToKML( kml.save(filename) return None - def meanLateralWindSpeed(self, dispersion_results): - print( - f'Lateral Surface Wind Speed - Mean Value: {np.mean(dispersion_results["lateralWind"]):0.3f} m/s' - ) - print( - f'Lateral Surface Wind Speed - Standard Deviation: {np.std(dispersion_results["lateralWind"]):0.3f} m/s' - ) - - def plotLateralWindSpeed(self, dispersion_results): - - self.meanLateralWindSpeed(dispersion_results) - - plt.figure() - plt.hist(dispersion_results["lateralWind"], bins=int(self.N**0.5)) - plt.title("Lateral Surface Wind Speed") - plt.xlabel("Velocity (m/s)") - plt.ylabel("Number of Occurences") - plt.show() - - def meanFrontalWindSpeed(self, dispersion_results): - print( - f'Frontal Surface Wind Speed - Mean Value: {np.mean(dispersion_results["frontalWind"]):0.3f} m/s' - ) - print( - f'Frontal Surface Wind Speed - Standard Deviation: {np.std(dispersion_results["frontalWind"]):0.3f} m/s' - ) - - def plotFrontalWindSpeed(self, dispersion_results): - - self.meanFrontalWindSpeed(dispersion_results) - - plt.figure() - plt.hist(dispersion_results["frontalWind"], bins=int(self.N**0.5)) - plt.title("Frontal Surface Wind Speed") - plt.xlabel("Velocity (m/s)") - plt.ylabel("Number of Occurences") - plt.show() - def info(self): + """Print information about the dispersion model. - dispersion_results = self.importResults(self.filename) - - self.meanApogeeAltitude(dispersion_results) - - self.meanOutOfRailVelocity(dispersion_results) - - self.meanStaticMargin(dispersion_results) - - self.meanLateralWindSpeed(dispersion_results) - - self.meanFrontalWindSpeed(dispersion_results) - - self.meanOutOfRailTime(dispersion_results) - - self.meanApogeeTime(dispersion_results) - - self.meanApogeeXPosition(dispersion_results) - - self.meanApogeeYPosition(dispersion_results) - - self.meanImpactTime(dispersion_results) - - self.meanImpactVelocity(dispersion_results) - - self.meanImpactXPosition(dispersion_results) - - self.meanImpactYPosition(dispersion_results) - - self.meanMaximumVelocity(dispersion_results) - - self.meanNumberOfParachuteEvents(dispersion_results) - - self.meanDrogueFullyInflatedTime(dispersion_results) + Returns + ------- + None + """ - self.meanDrogueFullyVelocity(dispersion_results) + print("Monte Carlo Simulation by RocketPy") + print("Data Source: ", self.filename) + print("Number of simulations: ", self.num_of_loaded_sims) + print("Results: ") + self.print_results() - self.meanDrogueTriggerTime(dispersion_results) + return None def allInfo(self): - dispersion_results = self.importResults(self.filename) - - self.plotEllipses(dispersion_results, self.image, self.realLandingPoint) - - self.plotApogeeAltitude(dispersion_results) - - self.plotOutOfRailVelocity(dispersion_results) - - self.plotStaticMargin(dispersion_results) - - self.plotLateralWindSpeed(dispersion_results) - - self.plotFrontalWindSpeed(dispersion_results) - - self.plotOutOfRailTime(dispersion_results) - - self.plotApogeeTime(dispersion_results) - - self.plotApogeeXPosition(dispersion_results) - - self.plotApogeeYPosition(dispersion_results) + """Print and plot information about the dispersion model and the results. - self.plotImpactTime(dispersion_results) - - self.plotImpactVelocity(dispersion_results) - - self.plotImpactXPosition(dispersion_results) - - self.plotImpactYPosition(dispersion_results) - - self.plotMaximumVelocity(dispersion_results) - - self.plotNumberOfParachuteEvents(dispersion_results) - - self.plotDrogueFullyInflatedTime(dispersion_results) - - self.plotDrogueFullyVelocity(dispersion_results) - - self.plotDrogueTriggerTime(dispersion_results) - - # Variables - - environment_inputs = { - "railLength": "required", - "gravity": 9.80665, - "date": None, - "latitude": 0, - "longitude": 0, - "elevation": 0, - "datum": "SIRGAS2000", - "timeZone": "UTC", - } - - solidmotor_inputs = { - "thrust": "required", - "burnOutTime": "required", - "totalImpulse": 0, - "grainNumber": "required", - "grainDensity": "required", - "grainOuterRadius": "required", - "grainInitialInnerRadius": "required", - "grainInitialHeight": "required", - "grainSeparation": 0, - "nozzleRadius": 0.0335, - "throatRadius": 0.0114, - } + Returns + ------- + None + """ + dispersion_results = self.dispersion_results - rocket_inputs = { - "mass": "required", - "inertiaI": "required", - "inertiaZ": "required", - "radius": "required", - "distanceRocketNozzle": "required", - "distanceRocketPropellant": "required", - "powerOffDrag": "required", - "powerOnDrag": "required", - } + print("Monte Carlo Simulation by RocketPy") + print("Data Source: ", self.filename) + print("Number of simulations: ", self.num_of_loaded_sims) + print("Results: ") + self.print_results() + print("Plotting results: ") + self.plotEllipses(dispersion_results=dispersion_results) + self.plot_results() - flight_inputs = { - "inclination": 80, - "heading": 90, - "initialSolution": None, - "terminateOnApogee": False, - "maxTime": 600, - "maxTimeStep": np.inf, - "minTimeStep": 0, - "rtol": 1e-6, - "atol": 6 * [1e-3] + 4 * [1e-6] + 3 * [1e-3], - "timeOvershoot": True, - "verbose": False, - } + return None diff --git a/rocketpy/Environment.py b/rocketpy/Environment.py index ccf7452ec..e3c847d8b 100644 --- a/rocketpy/Environment.py +++ b/rocketpy/Environment.py @@ -1,7 +1,5 @@ # -*- coding: utf-8 -*- -from .Function import Function - __author__ = "Giovani Hidalgo Ceotto, Guilherme Fernandes Alves, Lucas Azevedo Pezente, Oscar Mauricio Prada Ramirez, Lucas Kierulff Balabram" __copyright__ = "Copyright 20XX, RocketPy Team" __license__ = "MIT" @@ -18,6 +16,9 @@ import pytz import requests +from .Function import Function +from .supplement import geodesicToUtm, calculateEarthRadius + try: import netCDF4 except ImportError: @@ -56,7 +57,7 @@ class Environment: Value of Air's Gas Constant = 287.05287 J/K/Kg Gravity and Launch Rail Length: - Environment.rl : float + Environment.railLength : float Launch rail length in meters. Environment.gravity : float Positive value of gravitational acceleration in m/s^2. @@ -101,7 +102,7 @@ class Environment: Environment.elevArray: array Two-dimensional Array containing the elevation information Environment.topographicProfileActivated: bool - True if the user already set a topographic plofile + True if the user already set a topographic profile Atmosphere Static Conditions: Environment.maxExpectedHeight : float @@ -351,7 +352,7 @@ def __init__( None """ # Save launch rail length - self.rL = railLength + self.railLength = railLength # Save gravity value self.gravity = gravity @@ -2762,10 +2763,10 @@ def calculateSpeedOfSoundProfile(self): # Retrieve gas constant R and temperature T R = self.airGasConstant T = self.temperature - G = 1.4 # Unused variable, why? + G = 1.4 # Compute speed of sound using sqrt(gamma*R*T) - a = (1.4 * R * T) ** 0.5 + a = (G * R * T) ** 0.5 # Set new output for the calculated speed of sound a.setOutputs("Speed of Sound (m/s)") @@ -2865,7 +2866,7 @@ def info(self): """ # Print launch site details print("Launch Site Details") - print("\nLaunch Rail Length:", self.rL, " m") + print("\nLaunch Rail Length:", self.railLength, " m") time_format = "%Y-%m-%d %H:%M:%S" if self.date != None and "UTC" not in self.timeZone: print( @@ -3001,7 +3002,7 @@ def allInfo(self): # Print launch site details print("\n\nLaunch Site Details") - print("\nLaunch Rail Length:", self.rL, " m") + print("\nLaunch Rail Length:", self.railLength, " m") time_format = "%Y-%m-%d %H:%M:%S" if self.date != None and "UTC" not in self.timeZone: print( @@ -3304,7 +3305,7 @@ def allInfoReturned(self): # Dictionary creation, if not commented follows the SI info = dict( grav=self.gravity, - launch_rail_length=self.rL, + launch_rail_length=self.railLength, elevation=self.elevation, modelType=self.atmosphericModelType, modelTypeMaxExpectedHeight=self.maxExpectedHeight, @@ -3351,9 +3352,9 @@ def exportEnvironment(self, filename="environment"): # TODO: in the future, allow the user to select which format will be used (json, csv, etc.). Default must be JSON. # TODO: add self.exportEnvDictionary to the documentation - # TODO: find a way to documennt the workaround I've used on ma.getdata(self... + # TODO: find a way to document the workaround I've used on ma.getdata(self... self.exportEnvDictionary = { - "railLength": self.rL, + "railLength": self.railLength, "gravity": self.gravity, "date": [self.date.year, self.date.month, self.date.day, self.date.hour], "latitude": self.latitude, @@ -3395,336 +3396,6 @@ def exportEnvironment(self, filename="environment"): return None - # Auxiliary functions - Geodesic Coordinates - def geodesicToUtm(self, lat, lon, datum): - """Function which converts geodetic coordinates, i.e. lat/lon, to UTM - projection coordinates. Can be used only for latitudes between -80.00° - and 84.00° - - Parameters - ---------- - lat : float - The latitude coordinates of the point of analysis, must be contained - between -80.00° and 84.00° - lon : float - The longitude coordinates of the point of analysis, must be contained - between -180.00° and 180.00° - datum : string - The desired reference ellipsoide model, the following options are - available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default - is "SIRGAS2000", then this model will be used if the user make some - typing mistake - - Returns - ------- - x: float - East coordinate, always positive - y: - North coordinate, always positive - utmZone: int - The number of the UTM zone of the point of analysis, can vary between - 1 and 60 - utmLetter: string - The letter of the UTM zone of the point of analysis, can vary between - C and X, omitting the letters "I" and "O" - hemis: string - Returns "S" for southern hemisphere and "N" for Northern hemisphere - EW: string - Returns "W" for western hemisphere and "E" for eastern hemisphere - """ - - # Calculate the central meridian of UTM zone - if lon != 0: - signal = lon / abs(lon) - if signal > 0: - aux = lon - 3 - aux = aux * signal - div = aux // 6 - lon_mc = div * 6 + 3 - EW = "E" - else: - aux = lon + 3 - aux = aux * signal - div = aux // 6 - lon_mc = (div * 6 + 3) * signal - EW = "W" - else: - lon_mc = 3 - EW = "W|E" - - # Select the desired datum (i.e. the ellipsoid parameters) - if datum == "SAD69": - semiMajorAxis = 6378160.0 - flattening = 1 / 298.25 - elif datum == "WGS84": - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257223563 - elif datum == "NAD83": - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257024899 - else: - # SIRGAS2000 - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257223563 - - # Evaluate the hemisphere and determine the N coordinate at the Equator - if lat < 0: - N0 = 10000000 - hemis = "S" - else: - N0 = 0 - hemis = "N" - - # Convert the input lat and lon to radians - lat = lat * np.pi / 180 - lon = lon * np.pi / 180 - lon_mc = lon_mc * np.pi / 180 - - # Evaluate reference parameters - K0 = 1 - 1 / 2500 - e2 = 2 * flattening - flattening**2 - e2lin = e2 / (1 - e2) - - # Evaluate auxiliary parameters - A = e2 * e2 - B = A * e2 - C = np.sin(2 * lat) - D = np.sin(4 * lat) - E = np.sin(6 * lat) - F = (1 - e2 / 4 - 3 * A / 64 - 5 * B / 256) * lat - G = (3 * e2 / 8 + 3 * A / 32 + 45 * B / 1024) * C - H = (15 * A / 256 + 45 * B / 1024) * D - I = (35 * B / 3072) * E - - # Evaluate other reference parameters - n = semiMajorAxis / ((1 - e2 * (np.sin(lat) ** 2)) ** 0.5) - t = np.tan(lat) ** 2 - c = e2lin * (np.cos(lat) ** 2) - ag = (lon - lon_mc) * np.cos(lat) - m = semiMajorAxis * (F - G + H - I) - - # Evaluate new auxiliary parameters - J = (1 - t + c) * ag * ag * ag / 6 - K = (5 - 18 * t + t * t + 72 * c - 58 * e2lin) * (ag**5) / 120 - L = (5 - t + 9 * c + 4 * c * c) * ag * ag * ag * ag / 24 - M = (61 - 58 * t + t * t + 600 * c - 330 * e2lin) * (ag**6) / 720 - - # Evaluate the final coordinates - x = 500000 + K0 * n * (ag + J + K) - y = N0 + K0 * (m + n * np.tan(lat) * (ag * ag / 2 + L + M)) - - # Convert the output lat and lon to degrees - lat = lat * 180 / np.pi - lon = lon * 180 / np.pi - lon_mc = lon_mc * 180 / np.pi - - # Calculate the UTM zone number - utmZone = int((lon_mc + 183) / 6) - - # Calculate the UTM zone letter - letters = "CDEFGHJKLMNPQRSTUVWXX" - utmLetter = letters[int(80 + lat) >> 3] - - return x, y, utmZone, utmLetter, hemis, EW - - def utmToGeodesic(self, x, y, utmZone, hemis, datum): - """Function to convert UTM coordinates to geodesic coordinates - (i.e. latitude and longitude). The latitude should be between -80° - and 84° - - Parameters - ---------- - x : float - East UTM coordinate in meters - y : float - North UTM coordinate in meters - utmZone : int - The number of the UTM zone of the point of analysis, can vary between - 1 and 60 - hemis : string - Equals to "S" for southern hemisphere and "N" for Northern hemisphere - datum : string - The desired reference ellipsoide model, the following options are - available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default - is "SIRGAS2000", then this model will be used if the user make some - typing mistake - - Returns - ------- - lat: float - latitude of the analysed point - lon: float - latitude of the analysed point - """ - - if hemis == "N": - y = y + 10000000 - - # Calculate the Central Meridian from the UTM zone number - centralMeridian = utmZone * 6 - 183 # degrees - - # Select the desired datum - if datum == "SAD69": - semiMajorAxis = 6378160.0 - flattening = 1 / 298.25 - elif datum == "WGS84": - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257223563 - elif datum == "NAD83": - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257024899 - else: - # SIRGAS2000 - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257223563 - - # Calculate reference values - K0 = 1 - 1 / 2500 - e2 = 2 * flattening - flattening**2 - e2lin = e2 / (1 - e2) - e1 = (1 - (1 - e2) ** 0.5) / (1 + (1 - e2) ** 0.5) - - # Calculate auxiliary values - A = e2 * e2 - B = A * e2 - C = e1 * e1 - D = e1 * C - E = e1 * D - - m = (y - 10000000) / K0 - mi = m / (semiMajorAxis * (1 - e2 / 4 - 3 * A / 64 - 5 * B / 256)) - - # Calculate other auxiliary values - F = (3 * e1 / 2 - 27 * D / 32) * np.sin(2 * mi) - G = (21 * C / 16 - 55 * E / 32) * np.sin(4 * mi) - H = (151 * D / 96) * np.sin(6 * mi) - - lat1 = mi + F + G + H - c1 = e2lin * (np.cos(lat1) ** 2) - t1 = np.tan(lat1) ** 2 - n1 = semiMajorAxis / ((1 - e2 * (np.sin(lat1) ** 2)) ** 0.5) - quoc = (1 - e2 * np.sin(lat1) * np.sin(lat1)) ** 3 - r1 = semiMajorAxis * (1 - e2) / (quoc**0.5) - d = (x - 500000) / (n1 * K0) - - # Calculate other auxiliary values - I = (5 + 3 * t1 + 10 * c1 - 4 * c1 * c1 - 9 * e2lin) * d * d * d * d / 24 - J = ( - (61 + 90 * t1 + 298 * c1 + 45 * t1 * t1 - 252 * e2lin - 3 * c1 * c1) - * (d**6) - / 720 - ) - K = d - (1 + 2 * t1 + c1) * d * d * d / 6 - L = ( - (5 - 2 * c1 + 28 * t1 - 3 * c1 * c1 + 8 * e2lin + 24 * t1 * t1) - * (d**5) - / 120 - ) - - # Finally calculate the coordinates in lat/lot - lat = lat1 - (n1 * np.tan(lat1) / r1) * (d * d / 2 - I + J) - lon = centralMeridian * np.pi / 180 + (K + L) / np.cos(lat1) - - # Convert final lat/lon to Degrees - lat = lat * 180 / np.pi - lon = lon * 180 / np.pi - - return lat, lon - - def calculateEarthRadius(self, lat, datum): - """Simple function to calculate the Earth Radius at a specific latitude - based on ellipsoidal reference model (datum). The earth radius here is - assumed as the distance between the ellipsoid's center of gravity and a - point on ellipsoid surface at the desired - Pay attention: The ellipsoid is an approximation for the earth model and - will obviously output an estimate of the perfect distance between earth's - relief and its center of gravity. - - Parameters - ---------- - lat : float - latitude in which the Earth radius will be calculated - datum : string - The desired reference ellipsoide model, the following options are - available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default - is "SIRGAS2000", then this model will be used if the user make some - typing mistake - - Returns - ------- - float: - Earth Radius at the desired latitude in meters - """ - # Select the desired datum (i.e. the ellipsoid parameters) - if datum == "SAD69": - semiMajorAxis = 6378160.0 - flattening = 1 / 298.25 - elif datum == "WGS84": - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257223563 - elif datum == "NAD83": - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257024899 - else: - # SIRGAS2000 - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257223563 - - # Calculate the semi minor axis length - # semiMinorAxis = semiMajorAxis - semiMajorAxis*(flattening**(-1)) - semiMinorAxis = semiMajorAxis * (1 - flattening) - - # Convert latitude to radians - lat = lat * np.pi / 180 - - # Calculate the Earth Radius in meters - eRadius = np.sqrt( - ( - (np.cos(lat) * (semiMajorAxis**2)) ** 2 - + (np.sin(lat) * (semiMinorAxis**2)) ** 2 - ) - / ((np.cos(lat) * semiMajorAxis) ** 2 + (np.sin(lat) * semiMinorAxis) ** 2) - ) - - # Convert latitude to degress - lat = lat * 180 / np.pi - - return eRadius - - def decimalDegressToArcSeconds(self, angle): - """Function to convert an angle in decimal degrees to deg/min/sec. - Converts (°) to (° ' ") - - Parameters - ---------- - angle : float - The angle that you need convert to deg/min/sec. Must be given in - decimal degrees. - - Returns - ------- - deg: float - The degrees. - min: float - The arc minutes. 1 arc-minute = (1/60)*degree - sec: float - The arc Seconds. 1 arc-second = (1/3600)*degree - """ - - if angle < 0: - signal = -1 - else: - signal = 1 - - deg = (signal * angle) // 1 - min = abs(signal * angle - deg) * 60 // 1 - sec = abs((signal * angle - deg) * 60 - min) * 60 - # print("The angle {:f} is equals to {:.0f}º {:.0f}' {:.3f}'' ".format( - # angle, signal*deg, min, sec - # )) - - return deg, min, sec - def printEarthDetails(self): """[UNDER CONSTRUCTION] Function to print information about the Earth Model used in the @@ -3745,6 +3416,3 @@ def printEarthDetails(self): print("Gravity acceleration at launch site: Still not implemented :(") return None - - -from .utilities import calculateEarthRadius, geodesicToUtm diff --git a/rocketpy/EnvironmentAnalysis.py b/rocketpy/EnvironmentAnalysis.py index 8eeacbc0f..d03034587 100644 --- a/rocketpy/EnvironmentAnalysis.py +++ b/rocketpy/EnvironmentAnalysis.py @@ -3115,6 +3115,7 @@ def exportMeanProfiles(self, filename="export_env_analysis"): Exports the mean profiles of the weather data to a file in order to it be used as inputs on Environment Class by using the CustomAtmosphere model. + OBS: Not all units are allowed as inputs of Environment Class. Parameters ---------- diff --git a/rocketpy/Flight.py b/rocketpy/Flight.py index 46be6d4be..21682d520 100644 --- a/rocketpy/Flight.py +++ b/rocketpy/Flight.py @@ -17,7 +17,7 @@ from scipy import integrate from .Function import Function -from .utilities import * +from .supplement import * class Flight: @@ -610,8 +610,8 @@ def __init__( 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.env.railLength - abs(nozzle - upperRButton) + self.effective2RL = self.env.railLength - abs(nozzle - lowerRButton) # Flight initialization self.__init_post_process_variables() @@ -1419,7 +1419,6 @@ def uDot(self, t, u, postProcessing=False): if -1 * compStreamVzBn < 1: compAttackAngle = np.arccos(-compStreamVzBn) cLift = aerodynamicSurface.cl(compAttackAngle, freestreamMach) - cLift = aerodynamicSurface.cl(compAttackAngle, freestreamMach) # Component lift force magnitude compLift = ( 0.5 * rho * (compStreamSpeed**2) * referenceArea * cLift @@ -2755,6 +2754,19 @@ def attitudeFrequencyResponse(self): def staticMargin(self): return self.rocket.staticMargin + # Save important Static Margin values + @cached_property + def initialStaticMargin(self): + return self.staticMargin(0) + + @cached_property + def outOfRailStaticMargin(self): + return self.staticMargin(self.outOfRailTime) + + @cached_property + def finalStaticMargin(self): + return self.staticMargin(self.staticMargin(0)) + # Rail Button Forces @cached_property def railButton1NormalForce(self): diff --git a/rocketpy/Motor.py b/rocketpy/Motor.py index be9d62707..bf8f86ca0 100644 --- a/rocketpy/Motor.py +++ b/rocketpy/Motor.py @@ -154,6 +154,7 @@ def __init__( # Thrust parameters self.interpolate = interpolationMethod self.burnOutTime = burnOut + self.thrustSource = thrustSource # Check if thrustSource is csv, eng, function or other if isinstance(thrustSource, str): diff --git a/rocketpy/Rocket.py b/rocketpy/Rocket.py index f76fa5651..e9d4c7a7f 100644 --- a/rocketpy/Rocket.py +++ b/rocketpy/Rocket.py @@ -621,7 +621,7 @@ def addEllipticalFins( # Create a fin set as an object of EllipticalFins class finSet = EllipticalFins( - n, rootChord, span, distanceToCM, radius, cantAngle, airfoil, name + n, rootChord, span, distanceToCM, cantAngle, radius, airfoil, name ) # Add fin set to the list of aerodynamic surfaces @@ -726,8 +726,10 @@ def setRailButtons(self, distanceToCM, angularPosition=45): distanceToCM.reverse() # Save important attributes self.railButtons = self.railButtonPair(distanceToCM, angularPosition) - self.RBdistanceToCM = distanceToCM - self.angularPosition = angularPosition + # Saving in a special format just for dispersion class + self.positionFirstRailButton = distanceToCM[0] + self.positionSecondRailButton = distanceToCM[1] + self.railButtonAngularPosition = angularPosition return None def addCMEccentricity(self, x, y): @@ -815,6 +817,24 @@ def addThrustEccentricity(self, x, y): # Return self return self + def get_parachute_by_name(self, name): + """Returns parachute object by name. + + Parameters + ---------- + name : str + Name of the parachute. + + Returns + ------- + parachute : Parachute + Parachute object with the given name. + """ + for parachute in self.parachutes: + if parachute.name == name: + return parachute + raise ValueError("No parachute with name {} found.".format(name)) + def info(self): """Prints out a summary of the data and graphs available about the Rocket. diff --git a/rocketpy/__init__.py b/rocketpy/__init__.py index 681327681..8b2d94f04 100644 --- a/rocketpy/__init__.py +++ b/rocketpy/__init__.py @@ -30,3 +30,4 @@ from .Motor import HybridMotor, SolidMotor from .Rocket import Rocket from .utilities import * +from .supplement import * diff --git a/rocketpy/supplement.py b/rocketpy/supplement.py new file mode 100644 index 000000000..4901debca --- /dev/null +++ b/rocketpy/supplement.py @@ -0,0 +1,426 @@ +# -*- coding: utf-8 -*- +__author__ = "Guilherme Fernandes Alves" +__copyright__ = "Copyright 20XX, RocketPy Team" +__license__ = "MIT" + + +import numpy as np +import math + +# Geodesic calculations functions + + +def calculateEarthRadius(lat, datum="WGS84"): + """Simple function to calculate the Earth Radius at a specific latitude + based on ellipsoidal reference model (datum). The earth radius here is + assumed as the distance between the ellipsoid's center of gravity and a + point on ellipsoid surface at the desired + Pay attention: The ellipsoid is an approximation for the earth model and + will obviously output an estimate of the perfect distance between earth's + relief and its center of gravity. + + Parameters + ---------- + lat : float + latitude in which the Earth radius will be calculated + datum : string, optional + The desired reference ellipsoid model, the following options are + available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default + is "WSG84", then this model will be used if the user make some + typing mistake. + + Returns + ------- + float: + Earth Radius at the desired latitude in meters + """ + # Select the desired datum (i.e. the ellipsoid parameters) + if datum == "SAD69": + semiMajorAxis = 6378160.0 + flattening = 1 / 298.25 + elif datum == "SIRGAS2000": + semiMajorAxis = 6378137.0 + flattening = 1 / 298.257223563 + elif datum == "NAD83": + semiMajorAxis = 6378137.0 + flattening = 1 / 298.257024899 + else: + # Use WGS84 as default + semiMajorAxis = 6378137.0 + flattening = 1 / 298.257223563 + + # Calculate the semi minor axis length + semiMinorAxis = semiMajorAxis * (1 - flattening) + + # Convert latitude to radians + lat = lat * np.pi / 180 + + # Calculate the Earth Radius in meters + eRadius = np.sqrt( + ( + (np.cos(lat) * (semiMajorAxis**2)) ** 2 + + (np.sin(lat) * (semiMinorAxis**2)) ** 2 + ) + / ((np.cos(lat) * semiMajorAxis) ** 2 + (np.sin(lat) * semiMinorAxis) ** 2) + ) + + # Convert latitude back to degrees + lat = lat * 180 / np.pi + + return eRadius + + +def Haversine(lat0, lon0, lat1, lon1, eRadius=6.3781e6): + """Returns the distance between two points in meters. + The points are defined by their latitude and longitude coordinates. + + Parameters + ---------- + lat0 : float + Latitude of the first point, in degrees. + lon0 : float + Longitude of the first point, in degrees. + lat1 : float + Latitude of the second point, in degrees. + lon1 : float + Longitude of the second point, in degrees. + eRadius : float, optional + Earth's radius in meters. Default value is 6.3781e6. + + Returns + ------- + float + Distance between the two points in meters. + + """ + lat0_rad = math.radians(lat0) + lat1_rad = math.radians(lat1) + delta_lat_rad = math.radians(lat1 - lat0) + delta_lon_rad = math.radians(lon1 - lon0) + + a = ( + math.sin(delta_lat_rad / 2) ** 2 + + math.cos(lat0_rad) * math.cos(lat1_rad) * math.sin(delta_lon_rad / 2) ** 2 + ) + c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a)) + + return eRadius * c + + +def invertedHaversine(lat0, lon0, distance, bearing, eRadius=6.3781e6): + """Returns a tuple with new latitude and longitude coordinates considering + a displacement of a given distance in a given direction (bearing compass) + starting from a point defined by (lat0, lon0). This is the opposite of + Haversine function. + + Parameters + ---------- + lat0 : float + Origin latitude coordinate, in degrees. + lon0 : float + Origin longitude coordinate, in degrees. + distance : float + Distance from the origin point, in meters. + bearing : float + Azimuth (or bearing compass) from the origin point, in degrees. + eRadius : float, optional + Earth radius, in meters. Default value is 6.3781e6. + See the supplement.calculateEarthRadius() function for more accuracy. + + Returns + ------- + lat1 : float + New latitude coordinate, in degrees. + lon1 : float + New longitude coordinate, in degrees. + + """ + + # Convert coordinates to radians + lat0_rad = np.deg2rad(lat0) + lon0_rad = np.deg2rad(lon0) + + # Apply inverted Haversine formula + lat1_rad = math.asin( + math.sin(lat0_rad) * math.cos(distance / eRadius) + + math.cos(lat0_rad) * math.sin(distance / eRadius) * math.cos(bearing) + ) + + lon1_rad = lon0_rad + math.atan2( + math.sin(bearing) * math.sin(distance / eRadius) * math.cos(lat0_rad), + math.cos(distance / eRadius) - math.sin(lat0_rad) * math.sin(lat1_rad), + ) + + # Convert back to degrees and then return + lat1_deg = np.rad2deg(lat1_rad) + lon1_deg = np.rad2deg(lon1_rad) + + return lat1_deg, lon1_deg + + +def decimalDegreesToArcSeconds(angle): + """Function to convert an angle in decimal degrees to deg/min/sec. + Converts (°) to (° ' ") + + Parameters + ---------- + angle : float + The angle that you need convert to deg/min/sec. Must be given in + decimal degrees. + + Returns + ------- + deg: float + The degrees. + min: float + The arc minutes. 1 arc-minute = (1/60)*degree + sec: float + The arc Seconds. 1 arc-second = (1/3600)*degree + """ + + if angle < 0: + signal = -1 + else: + signal = 1 + + deg = (signal * angle) // 1 + min = abs(signal * angle - deg) * 60 // 1 + sec = abs((signal * angle - deg) * 60 - min) * 60 + + return deg, min, sec + + +def geodesicToUtm(lat, lon, datum="WGS84"): + """Function which converts geodesic coordinates, i.e. lat/lon, to UTM + projection coordinates. Can be used only for latitudes between -80.00° + and 84.00° + + Parameters + ---------- + lat : float + The latitude coordinates of the point of analysis, must be contained + between -80.00° and 84.00° + lon : float + The longitude coordinates of the point of analysis, must be contained + between -180.00° and 180.00° + datum : string + The desired reference ellipsoid model, the following options are + available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default + is "WGS84", then this model will be used if the user make some + typing mistake + + Returns + ------- + x: float + East coordinate in meters, always positive + y: + North coordinate in meters, always positive + utmZone: int + The number of the UTM zone of the point of analysis, can vary between + 1 and 60 + utmLetter: string + The letter of the UTM zone of the point of analysis, can vary between + C and X, omitting the letters "I" and "O" + NorthSouthHemisphere: string + Returns "S" for southern hemisphere and "N" for Northern hemisphere + EastWestHemisphere: string + Returns "W" for western hemisphere and "E" for eastern hemisphere + """ + + # Calculate the central meridian of UTM zone + if lon != 0: + # signal = 1 for positive longitude and -1 for negative longitude + signal = lon / abs(lon) + if signal > 0: + aux = lon - 3 + aux = aux * signal + div = aux // 6 + lon_mc = div * 6 + 3 + EastWestHemisphere = "E" # Eastern hemisphere + else: + aux = lon + 3 + aux = aux * signal + div = aux // 6 + lon_mc = (div * 6 + 3) * signal + EastWestHemisphere = "W" # Western hemisphere + else: + # If the longitude is zero, the central meridian receives number 3 + lon_mc = 3 + EastWestHemisphere = "W|E" + + # Select the desired datum (i.e. the ellipsoid parameters) + # TODO: Create separate function that returns the ellipsoid parameters + if datum == "SAD69": + semiMajorAxis = 6378160.0 + flattening = 1 / 298.25 + elif datum == "SIRGAS2000": + semiMajorAxis = 6378137.0 + flattening = 1 / 298.257223563 + elif datum == "NAD83": + semiMajorAxis = 6378137.0 + flattening = 1 / 298.257024899 + else: + # WGS84 + semiMajorAxis = 6378137.0 + flattening = 1 / 298.257223563 + + # Evaluate the S/N hemisphere and determine the N coordinate at the Equator + if lat < 0: + N0 = 10000000 + NorthSouthHemisphere = "S" + else: + N0 = 0 + NorthSouthHemisphere = "N" + + # Convert the input lat and lon to radians + lat = lat * np.pi / 180 + lon = lon * np.pi / 180 + lon_mc = lon_mc * np.pi / 180 + + # Evaluate reference parameters + K0 = 1 - 1 / 2500 + e2 = 2 * flattening - flattening**2 + e2lin = e2 / (1 - e2) + + # Evaluate auxiliary parameters + A = e2 * e2 + B = A * e2 + C = np.sin(2 * lat) + D = np.sin(4 * lat) + E = np.sin(6 * lat) + F = (1 - e2 / 4 - 3 * A / 64 - 5 * B / 256) * lat + G = (3 * e2 / 8 + 3 * A / 32 + 45 * B / 1024) * C + H = (15 * A / 256 + 45 * B / 1024) * D + I = (35 * B / 3072) * E + + # Evaluate other reference parameters + n = semiMajorAxis / ((1 - e2 * (np.sin(lat) ** 2)) ** 0.5) + t = np.tan(lat) ** 2 + c = e2lin * (np.cos(lat) ** 2) + ag = (lon - lon_mc) * np.cos(lat) + m = semiMajorAxis * (F - G + H - I) + + # Evaluate new auxiliary parameters + J = (1 - t + c) * ag * ag * ag / 6 + K = (5 - 18 * t + t * t + 72 * c - 58 * e2lin) * (ag**5) / 120 + L = (5 - t + 9 * c + 4 * c * c) * ag * ag * ag * ag / 24 + M = (61 - 58 * t + t * t + 600 * c - 330 * e2lin) * (ag**6) / 720 + + # Evaluate the final coordinates + x = 500000 + K0 * n * (ag + J + K) + y = N0 + K0 * (m + n * np.tan(lat) * (ag * ag / 2 + L + M)) + + # Convert the output lat and lon to degrees + lat = lat * 180 / np.pi + lon = lon * 180 / np.pi + lon_mc = lon_mc * 180 / np.pi + + # Calculate the UTM zone number + utmZone = int((lon_mc + 183) / 6) + + # Calculate the UTM zone letter + letters = "CDEFGHJKLMNPQRSTUVWXX" + utmLetter = letters[int(80 + lat) >> 3] + + return x, y, utmZone, utmLetter, NorthSouthHemisphere, EastWestHemisphere + + +def utmToGeodesic(x, y, utmZone, NorthSouthHemisphere, datum): + """Function to convert UTM coordinates to geodesic coordinates + (i.e. latitude and longitude). The latitude should be between -80° + and 84° + + Parameters + ---------- + x : float + East UTM coordinate in meters + y : float + North UTM coordinate in meters + utmZone : int + The number of the UTM zone of the point of analysis, can vary between + 1 and 60 + NorthSouthHemisphere : string + Equals to "S" for southern hemisphere and "N" for Northern hemisphere + datum : string + The desired reference ellipsoid model, the following options are + available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default + is "WGS84", then this model will be used if the user make some + typing mistake + + Returns + ------- + lat: float + latitude of the analyzed point + lon: float + latitude of the analyzed point + """ + + if NorthSouthHemisphere == "N": + y = y + 10000000 + + # Calculate the Central Meridian from the UTM zone number + centralMeridian = utmZone * 6 - 183 # degrees + + # Select the desired datum + if datum == "SAD69": + semiMajorAxis = 6378160.0 + flattening = 1 / 298.25 + elif datum == "SIRGAS2000": + semiMajorAxis = 6378137.0 + flattening = 1 / 298.257223563 + elif datum == "NAD83": + semiMajorAxis = 6378137.0 + flattening = 1 / 298.257024899 + else: + # WGS84 + semiMajorAxis = 6378137.0 + flattening = 1 / 298.257223563 + + # Calculate reference values + K0 = 1 - 1 / 2500 + e2 = 2 * flattening - flattening**2 + e2lin = e2 / (1 - e2) + e1 = (1 - (1 - e2) ** 0.5) / (1 + (1 - e2) ** 0.5) + + # Calculate auxiliary values + A = e2 * e2 + B = A * e2 + C = e1 * e1 + D = e1 * C + E = e1 * D + + m = (y - 10000000) / K0 + mi = m / (semiMajorAxis * (1 - e2 / 4 - 3 * A / 64 - 5 * B / 256)) + + # Calculate other auxiliary values + F = (3 * e1 / 2 - 27 * D / 32) * np.sin(2 * mi) + G = (21 * C / 16 - 55 * E / 32) * np.sin(4 * mi) + H = (151 * D / 96) * np.sin(6 * mi) + + lat1 = mi + F + G + H + c1 = e2lin * (np.cos(lat1) ** 2) + t1 = np.tan(lat1) ** 2 + n1 = semiMajorAxis / ((1 - e2 * (np.sin(lat1) ** 2)) ** 0.5) + qc = (1 - e2 * np.sin(lat1) * np.sin(lat1)) ** 3 + r1 = semiMajorAxis * (1 - e2) / (qc**0.5) + d = (x - 500000) / (n1 * K0) + + # Calculate other auxiliary values + I = (5 + 3 * t1 + 10 * c1 - 4 * c1 * c1 - 9 * e2lin) * d * d * d * d / 24 + J = ( + (61 + 90 * t1 + 298 * c1 + 45 * t1 * t1 - 252 * e2lin - 3 * c1 * c1) + * (d**6) + / 720 + ) + K = d - (1 + 2 * t1 + c1) * d * d * d / 6 + L = (5 - 2 * c1 + 28 * t1 - 3 * c1 * c1 + 8 * e2lin + 24 * t1 * t1) * (d**5) / 120 + + # Finally calculate the coordinates in lat/lot + lat = lat1 - (n1 * np.tan(lat1) / r1) * (d * d / 2 - I + J) + lon = centralMeridian * np.pi / 180 + (K + L) / np.cos(lat1) + + # Convert final lat/lon to Degrees + lat = lat * 180 / np.pi + lon = lon * 180 / np.pi + + return lat, lon diff --git a/rocketpy/utilities.py b/rocketpy/utilities.py index 7de584db4..2664fd0b8 100644 --- a/rocketpy/utilities.py +++ b/rocketpy/utilities.py @@ -2,11 +2,8 @@ __author__ = "Franz Masatoshi Yuri, Lucas Kierulff Balabram, Guilherme Fernandes Alves, Bruno Abdulklech Sorban" __copyright__ = "Copyright 20XX, RocketPy Team" __license__ = "MIT" -import math -import matplotlib.pyplot as plt import numpy as np -from matplotlib.patches import Ellipse from scipy.integrate import solve_ivp from .Environment import Environment @@ -241,422 +238,3 @@ def create_dispersion_dictionary(filename): except: analysis_parameters[row[0].strip()] = "" return analysis_parameters - - -# Geodesic calculations functions - - -def calculateEarthRadius(lat, datum="WGS84"): - """Simple function to calculate the Earth Radius at a specific latitude - based on ellipsoidal reference model (datum). The earth radius here is - assumed as the distance between the ellipsoid's center of gravity and a - point on ellipsoid surface at the desired - Pay attention: The ellipsoid is an approximation for the earth model and - will obviously output an estimate of the perfect distance between earth's - relief and its center of gravity. - - Parameters - ---------- - lat : float - latitude in which the Earth radius will be calculated - datum : string, optional - The desired reference ellipsoid model, the following options are - available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default - is "WSG84", then this model will be used if the user make some - typing mistake. - - Returns - ------- - float: - Earth Radius at the desired latitude in meters - """ - # Select the desired datum (i.e. the ellipsoid parameters) - if datum == "SAD69": - semiMajorAxis = 6378160.0 - flattening = 1 / 298.25 - elif datum == "SIRGAS2000": - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257223563 - elif datum == "NAD83": - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257024899 - else: - # Use WGS84 as default - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257223563 - - # Calculate the semi minor axis length - semiMinorAxis = semiMajorAxis * (1 - flattening) - - # Convert latitude to radians - lat = lat * np.pi / 180 - - # Calculate the Earth Radius in meters - eRadius = np.sqrt( - ( - (np.cos(lat) * (semiMajorAxis**2)) ** 2 - + (np.sin(lat) * (semiMinorAxis**2)) ** 2 - ) - / ((np.cos(lat) * semiMajorAxis) ** 2 + (np.sin(lat) * semiMinorAxis) ** 2) - ) - - # Convert latitude back to degrees - lat = lat * 180 / np.pi - - return eRadius - - -def Haversine(lat0, lon0, lat1, lon1, eRadius=6.3781e6): - """Returns the distance between two points in meters. - The points are defined by their latitude and longitude coordinates. - - Parameters - ---------- - lat0 : float - Latitude of the first point, in degrees. - lon0 : float - Longitude of the first point, in degrees. - lat1 : float - Latitude of the second point, in degrees. - lon1 : float - Longitude of the second point, in degrees. - eRadius : float, optional - Earth's radius in meters. Default value is 6.3781e6. - - Returns - ------- - float - Distance between the two points in meters. - - """ - lat0_rad = math.radians(lat0) - lat1_rad = math.radians(lat1) - delta_lat_rad = math.radians(lat1 - lat0) - delta_lon_rad = math.radians(lon1 - lon0) - - a = ( - math.sin(delta_lat_rad / 2) ** 2 - + math.cos(lat0_rad) * math.cos(lat1_rad) * math.sin(delta_lon_rad / 2) ** 2 - ) - c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a)) - - return eRadius * c - - -def invertedHaversine(lat0, lon0, distance, bearing, eRadius=6.3781e6): - """Returns a tuple with new latitude and longitude coordinates considering - a displacement of a given distance in a given direction (bearing compass) - starting from a point defined by (lat0, lon0). This is the opposite of - Haversine function. - - Parameters - ---------- - lat0 : float - Origin latitude coordinate, in degrees. - lon0 : float - Origin longitude coordinate, in degrees. - distance : float - Distance from the origin point, in meters. - bearing : float - Azimuth (or bearing compass) from the origin point, in degrees. - eRadius : float, optional - Earth radius, in meters. Default value is 6.3781e6. - See the utilities.calculateEarthRadius() function for more accuracy. - - Returns - ------- - lat1 : float - New latitude coordinate, in degrees. - lon1 : float - New longitude coordinate, in degrees. - - """ - - # Convert coordinates to radians - lat0_rad = np.deg2rad(lat0) - lon0_rad = np.deg2rad(lon0) - - # Apply inverted Haversine formula - lat1_rad = math.asin( - math.sin(lat0_rad) * math.cos(distance / eRadius) - + math.cos(lat0_rad) * math.sin(distance / eRadius) * math.cos(bearing) - ) - - lon1_rad = lon0_rad + math.atan2( - math.sin(bearing) * math.sin(distance / eRadius) * math.cos(lat0_rad), - math.cos(distance / eRadius) - math.sin(lat0_rad) * math.sin(lat1_rad), - ) - - # Convert back to degrees and then return - lat1_deg = np.rad2deg(lat1_rad) - lon1_deg = np.rad2deg(lon1_rad) - - return lat1_deg, lon1_deg - - -def decimalDegreesToArcSeconds(angle): - """Function to convert an angle in decimal degrees to deg/min/sec. - Converts (°) to (° ' ") - - Parameters - ---------- - angle : float - The angle that you need convert to deg/min/sec. Must be given in - decimal degrees. - - Returns - ------- - deg: float - The degrees. - min: float - The arc minutes. 1 arc-minute = (1/60)*degree - sec: float - The arc Seconds. 1 arc-second = (1/3600)*degree - """ - - if angle < 0: - signal = -1 - else: - signal = 1 - - deg = (signal * angle) // 1 - min = abs(signal * angle - deg) * 60 // 1 - sec = abs((signal * angle - deg) * 60 - min) * 60 - - return deg, min, sec - - -def geodesicToUtm(lat, lon, datum="WGS84"): - """Function which converts geodesic coordinates, i.e. lat/lon, to UTM - projection coordinates. Can be used only for latitudes between -80.00° - and 84.00° - - Parameters - ---------- - lat : float - The latitude coordinates of the point of analysis, must be contained - between -80.00° and 84.00° - lon : float - The longitude coordinates of the point of analysis, must be contained - between -180.00° and 180.00° - datum : string - The desired reference ellipsoid model, the following options are - available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default - is "WGS84", then this model will be used if the user make some - typing mistake - - Returns - ------- - x: float - East coordinate in meters, always positive - y: - North coordinate in meters, always positive - utmZone: int - The number of the UTM zone of the point of analysis, can vary between - 1 and 60 - utmLetter: string - The letter of the UTM zone of the point of analysis, can vary between - C and X, omitting the letters "I" and "O" - NorthSouthHemisphere: string - Returns "S" for southern hemisphere and "N" for Northern hemisphere - EastWestHemisphere: string - Returns "W" for western hemisphere and "E" for eastern hemisphere - """ - - # Calculate the central meridian of UTM zone - if lon != 0: - # signal = 1 for positive longitude and -1 for negative longitude - signal = lon / abs(lon) - if signal > 0: - aux = lon - 3 - aux = aux * signal - div = aux // 6 - lon_mc = div * 6 + 3 - EastWestHemisphere = "E" # Eastern hemisphere - else: - aux = lon + 3 - aux = aux * signal - div = aux // 6 - lon_mc = (div * 6 + 3) * signal - EastWestHemisphere = "W" # Western hemisphere - else: - # If the longitude is zero, the central meridian receives number 3 - lon_mc = 3 - EastWestHemisphere = "W|E" - - # Select the desired datum (i.e. the ellipsoid parameters) - # TODO: Create separate function that returns the ellipsoid parameters - if datum == "SAD69": - semiMajorAxis = 6378160.0 - flattening = 1 / 298.25 - elif datum == "SIRGAS2000": - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257223563 - elif datum == "NAD83": - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257024899 - else: - # WGS84 - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257223563 - - # Evaluate the S/N hemisphere and determine the N coordinate at the Equator - if lat < 0: - N0 = 10000000 - NorthSouthHemisphere = "S" - else: - N0 = 0 - NorthSouthHemisphere = "N" - - # Convert the input lat and lon to radians - lat = lat * np.pi / 180 - lon = lon * np.pi / 180 - lon_mc = lon_mc * np.pi / 180 - - # Evaluate reference parameters - K0 = 1 - 1 / 2500 - e2 = 2 * flattening - flattening**2 - e2lin = e2 / (1 - e2) - - # Evaluate auxiliary parameters - A = e2 * e2 - B = A * e2 - C = np.sin(2 * lat) - D = np.sin(4 * lat) - E = np.sin(6 * lat) - F = (1 - e2 / 4 - 3 * A / 64 - 5 * B / 256) * lat - G = (3 * e2 / 8 + 3 * A / 32 + 45 * B / 1024) * C - H = (15 * A / 256 + 45 * B / 1024) * D - I = (35 * B / 3072) * E - - # Evaluate other reference parameters - n = semiMajorAxis / ((1 - e2 * (np.sin(lat) ** 2)) ** 0.5) - t = np.tan(lat) ** 2 - c = e2lin * (np.cos(lat) ** 2) - ag = (lon - lon_mc) * np.cos(lat) - m = semiMajorAxis * (F - G + H - I) - - # Evaluate new auxiliary parameters - J = (1 - t + c) * ag * ag * ag / 6 - K = (5 - 18 * t + t * t + 72 * c - 58 * e2lin) * (ag**5) / 120 - L = (5 - t + 9 * c + 4 * c * c) * ag * ag * ag * ag / 24 - M = (61 - 58 * t + t * t + 600 * c - 330 * e2lin) * (ag**6) / 720 - - # Evaluate the final coordinates - x = 500000 + K0 * n * (ag + J + K) - y = N0 + K0 * (m + n * np.tan(lat) * (ag * ag / 2 + L + M)) - - # Convert the output lat and lon to degrees - lat = lat * 180 / np.pi - lon = lon * 180 / np.pi - lon_mc = lon_mc * 180 / np.pi - - # Calculate the UTM zone number - utmZone = int((lon_mc + 183) / 6) - - # Calculate the UTM zone letter - letters = "CDEFGHJKLMNPQRSTUVWXX" - utmLetter = letters[int(80 + lat) >> 3] - - return x, y, utmZone, utmLetter, NorthSouthHemisphere, EastWestHemisphere - - -def utmToGeodesic(x, y, utmZone, NorthSouthHemisphere, datum): - """Function to convert UTM coordinates to geodesic coordinates - (i.e. latitude and longitude). The latitude should be between -80° - and 84° - - Parameters - ---------- - x : float - East UTM coordinate in meters - y : float - North UTM coordinate in meters - utmZone : int - The number of the UTM zone of the point of analysis, can vary between - 1 and 60 - NorthSouthHemisphere : string - Equals to "S" for southern hemisphere and "N" for Northern hemisphere - datum : string - The desired reference ellipsoid model, the following options are - available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default - is "WGS84", then this model will be used if the user make some - typing mistake - - Returns - ------- - lat: float - latitude of the analyzed point - lon: float - latitude of the analyzed point - """ - - if NorthSouthHemisphere == "N": - y = y + 10000000 - - # Calculate the Central Meridian from the UTM zone number - centralMeridian = utmZone * 6 - 183 # degrees - - # Select the desired datum - if datum == "SAD69": - semiMajorAxis = 6378160.0 - flattening = 1 / 298.25 - elif datum == "SIRGAS2000": - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257223563 - elif datum == "NAD83": - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257024899 - else: - # WGS84 - semiMajorAxis = 6378137.0 - flattening = 1 / 298.257223563 - - # Calculate reference values - K0 = 1 - 1 / 2500 - e2 = 2 * flattening - flattening**2 - e2lin = e2 / (1 - e2) - e1 = (1 - (1 - e2) ** 0.5) / (1 + (1 - e2) ** 0.5) - - # Calculate auxiliary values - A = e2 * e2 - B = A * e2 - C = e1 * e1 - D = e1 * C - E = e1 * D - - m = (y - 10000000) / K0 - mi = m / (semiMajorAxis * (1 - e2 / 4 - 3 * A / 64 - 5 * B / 256)) - - # Calculate other auxiliary values - F = (3 * e1 / 2 - 27 * D / 32) * np.sin(2 * mi) - G = (21 * C / 16 - 55 * E / 32) * np.sin(4 * mi) - H = (151 * D / 96) * np.sin(6 * mi) - - lat1 = mi + F + G + H - c1 = e2lin * (np.cos(lat1) ** 2) - t1 = np.tan(lat1) ** 2 - n1 = semiMajorAxis / ((1 - e2 * (np.sin(lat1) ** 2)) ** 0.5) - qc = (1 - e2 * np.sin(lat1) * np.sin(lat1)) ** 3 - r1 = semiMajorAxis * (1 - e2) / (qc**0.5) - d = (x - 500000) / (n1 * K0) - - # Calculate other auxiliary values - I = (5 + 3 * t1 + 10 * c1 - 4 * c1 * c1 - 9 * e2lin) * d * d * d * d / 24 - J = ( - (61 + 90 * t1 + 298 * c1 + 45 * t1 * t1 - 252 * e2lin - 3 * c1 * c1) - * (d**6) - / 720 - ) - K = d - (1 + 2 * t1 + c1) * d * d * d / 6 - L = (5 - 2 * c1 + 28 * t1 - 3 * c1 * c1 + 8 * e2lin + 24 * t1 * t1) * (d**5) / 120 - - # Finally calculate the coordinates in lat/lot - lat = lat1 - (n1 * np.tan(lat1) / r1) * (d * d / 2 - I + J) - lon = centralMeridian * np.pi / 180 + (K + L) / np.cos(lat1) - - # Convert final lat/lon to Degrees - lat = lat * 180 / np.pi - lon = lon * 180 / np.pi - - return lat, lon diff --git a/tests/test_rocket.py b/tests/test_rocket.py index 29710c044..c7b185375 100644 --- a/tests/test_rocket.py +++ b/tests/test_rocket.py @@ -361,7 +361,6 @@ def test_add_trapezoidal_fins_sweep_angle( # Check center of pressure translate = 0.55829 + 0.71971 - cpz = FinSet.cp[2] assert translate - cpz == pytest.approx(expected_fin_cpz, 0.01) # Check lift coefficient derivative