diff --git a/TODO.md b/TODO.md deleted file mode 100644 index 1ecf698a5..000000000 --- a/TODO.md +++ /dev/null @@ -1,41 +0,0 @@ -# New Features that Contributors can Implement - -| Feature Name | Description | Programing Difficulty | Theoretical Knowledge Required | Current stage | -| ----------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | --------------------- | ---------------------------------- | --------------------------------------- | -| Rail Buttons | Finish rail buttons implementation, allowing for the calculation of loads they carry during launch. | Easy | Basic Dynamics | Not started | -| Custom Aerodynamic Surfaces | Create a new method in the Rocket Class that allows for custom aerodynamic coefficients such as Drag, Lift and Moment. | Easy | Basic Dynamics and Aerodynamics | Not started | -| Custom JSON weather files | Add option to import weather conditions as JOSN files such as those used by Windy.tv | Easy | JSON | Almost finished | -| Airfoil profiles | Currently, the Rocket class supports only flat fins, although the fins in real life can have more sophisticated aerodynamic profiles. To start solving this problem, we could allow the user to change the CLalpha of the fins before simulating. | Easy | None | Finished by @brunosorban | -| Geodesic models | Make RocketPy capable of converting impact coordinates to useful lat / lon coordinates, based on different models of cartographic representation, such as WGS84 or NAD83 data. | Easy | None | Almost finished | -| Recovery Class | Create a new class regarding all parachute calculations and analysis. | Easy | Bascic OOP | Not started | -| Rocket Add Discrete Controller Method | Add a method to the Rocket Class that allows the user to set a time-discrete, closed-loop, control function to be simulated during the flight. | Easy | Basic Dynamics and Basic Control | Not started | -| Rocket Add Continuous Controller Method | Add a method to the Rocket Class that allows the user to set a time-continuous, closed-loop, control function to be simulated during the flight. | Easy | Basic Dynamics and Basic Control | Not started | -| Parachute initial force | Calculate the Parachute deployment impact forces | Medium | Rocketry Parachute literature | Not started | -| Graphical User Interface | Desktop or Web graphical interface that can be used to run RocketPy in a more user friendly way, possibly also writing the code so the user can learn how to used RocketPy as a module as well. | Medium | None | Not started | -| 3D Rocket Attitude Animation | Create a 3D animation of the rocket's attitude during flight, allowing the engineer to better visualize the behavior of the rocket. Can also me used for marketing purposes. | Medium | Quaternions | Almost finished | -| Elevation profile data | Since the Earth is not flat, we can implement some functions that allow us to more accurately determine the impact coordinates by researching the point at which the rocket crosses the Earth's relief. Good candidates for this job are SRTM data. | Medium | None | Started by @Gui-FernandesBR | -| Export pressures | Improve the exportation of pressure values to facilitate the simulation of the apoggee detection algorithm. It currently supports the "noisy" pressure output only after parachutes have been defined. | Medium | None | Not started | -| PyRX + RocketPy | Integrate the GUI telemetry software PyRX with RocketPy | Medium | None | Not started | -| Maximum distance | Calculate the maximum distance between the rocket and the Telemetry base station. The objective is to improve the choosing and testing of the antennas. | Medium | None | Not started | -| Optimize apogee detection algorithm | Optimize the detection filter coeffitients for the specific launch. | Hard | Adaptive Digital Filtering and C++ | Started by @guilhermebene and @Lucas-KB | -| Plane Motion Flight Phase | Modified dynamics equations. | Hard | Advanced Dynamics | Not started | -| Hybrid/Liquid Motor | Derive the equations of motion for non-constant propellant center of mass. | Hard | Advanced Dynamics. | Started by @lucasfourier | -| Integration with Fusion360 (or similar) | Make RocketPy able to read a 3D geometry file that contains the entire Rocket description. This could save us a few minutes, as we would no longer need to measure multiple distances in our CAD software. | Hard | Autodesk Fusion 360 API | Not started | -| Parachute oscillations | Analysis of all Parachute stability | Hard | Rocketry Parachute literature | Not started | -| Multibody Dynamics Parachute Flight Phase | Currently, descent under parachute is simulated as a 3 degree of freedom (DOF) system. However, parachutes can and should be modeled separately from the rocket, creating a 6 (DOF) system, 2 for the parachute (spherical constrain) and 3 for the rocket. This can greatly enhance the prediction of the landing point. | Hard | Dynamics | Not started | -| Deployable Payload simulation | Allow RocketPy to simulate the trajectory of objects deployed during flight | Hard | Dynamics | Not started | - -# New Usage Examples and Documentation that Contributors can Create - -| Feature Name | Description | Programing Difficulty | Theoretical Knowledge Required | Current stage | -| ---------------------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------- | --------------------- | ------------------------------------------------------------- | ------------------------- | -| Rocket Class Notebook | Notebook describing functionalities of the Rocket Class. | Easy | Basic | Not started | -| SolidMotor Class Notebook | Notebook describing functionalities of the SolidMotor Class. | Easy | Basic | Almost finished | -| Flight Class Notebook | Notebook describing functionalities of the Flight Class. | Easy | Basic | Not started | -| Dispersion Notebook | Notebook describing how to run dispersion analysis and process results. | Medium | Monte Carlo | Finished | -| Multi-Stage Rockets | Notebook showing how RocketPy can be used to simulate multistage rockets. | Medium | Basic | Not started | -| Neural Network Controls | Implementation of neural networks to control the 6 degrees of freedom of the rocket during flight. | Hard | Neural Networks, Very Basic Control Theory | Not started | -| Two-way Coupling with Fluent | Make use of discrete control method to couple an aerodynamics CFD simulation in Ansys Fluent and flight simulation using RocketPy, getting accurate aerodynamic forces. | Hard | Ansys As a Server (AAS), Ansys Fluent, Aerodynamics, Dynamics | Not started | -| Communication platforms | What if we create a server on discord or WhatsApp so everyone interested could talk and contribute more easialy? | Easy | None | Started by @giovaniceotto | -| JOSS publication | Write and publish a paper on the Journal of Open Source Software | Hard | Scientific writing | Started by @giovaniceotto | -| Github page | Create a website page using the github tools. See more in https://pages.github.com/ | Medium | None | Not Started | diff --git a/docs/technical/aerodynamics/Roll_Equations.pdf b/docs/technical/aerodynamics/Roll_Equations.pdf index 758a6d616..c6c71e646 100644 Binary files a/docs/technical/aerodynamics/Roll_Equations.pdf and b/docs/technical/aerodynamics/Roll_Equations.pdf differ diff --git a/requirements.txt b/requirements.txt index 30001b109..4a59f377b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,4 +2,5 @@ numpy>=1.0 scipy>=1.0 matplotlib>=3.0 netCDF4>=1.4 -requests \ No newline at end of file +requests +pytz diff --git a/rocketpy/Environment.py b/rocketpy/Environment.py index 675803c40..0ce542171 100644 --- a/rocketpy/Environment.py +++ b/rocketpy/Environment.py @@ -9,6 +9,7 @@ import bisect import warnings import time +import pytz from datetime import datetime, timedelta from inspect import signature, getsourcelines from collections import namedtuple @@ -25,12 +26,27 @@ try: import netCDF4 - from netCDF4 import Dataset except ImportError: + has_netCDF4 = False warnings.warn( "Unable to load netCDF4. NetCDF files and OPeNDAP will not be imported.", ImportWarning, ) +else: + has_netCDF4 = True + + +def requires_netCDF4(func): + def wrapped_func(*args, **kwargs): + if has_netCDF4: + func(*args, **kwargs) + else: + raise ImportError( + "This feature requires netCDF4 to be installed. Install it with `pip install netCDF4`" + ) + + return wrapped_func + from .Function import Function @@ -80,7 +96,11 @@ class Environment: Environment.elevation : float Launch site elevation. Environment.date : datetime - Date time of launch. + Date time of launch in UTC. + Environment.localDate : datetime + Date time of launch in the local time zone, defined by Environment.timeZone. + Environment.timeZone : string + Local time zone specification. See pytz for time zone info. Topographic information: Environment.elevLonArray: array @@ -290,6 +310,7 @@ def __init__( longitude=0, elevation=0, datum="SIRGAS2000", + timeZone="UTC", ): """Initialize Environment class, saving launch rail length, launch date, location coordinates and elevation. Note that @@ -326,7 +347,13 @@ def __init__( 'Open-Elevation' which uses the Open-Elevation API to find elevation data. For this option, latitude and longitude must also be specified. Default value is 0. - datum: + 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. + timeZone : string, optional + Name of the time zone. To see all time zones, import pytz and run Returns ------- @@ -343,12 +370,14 @@ def __init__( # Save date if date != None: - self.setDate(date) + self.setDate(date, timeZone) else: self.date = None + self.localDate = None + self.timeZone = None # Initialize constants - self.earthRadius = 6.3781 * (10 ** 6) + self.earthRadius = 6.3781 * (10**6) self.airGasConstant = 287.05287 # in J/K/Kg # Initialize atmosphere @@ -378,21 +407,30 @@ def __init__( return None - def setDate(self, date): + def setDate(self, date, timeZone="UTC"): """Set date and time of launch and update weather conditions if date dependent atmospheric model is used. Parameters ---------- - date : Date - Date object specifying launch date and time. + date : Datetime + Datetime object specifying launch date and time. + timeZone : string, optional + Name of the time zone. To see all time zones, import pytz and run + print(pytz.all_timezones). Default time zone is "UTC". Return ------ None """ - # Store date - self.date = datetime(*date) + # Store date and configure time zone + self.timeZone = timeZone + tz = pytz.timezone(self.timeZone) + localDate = datetime(*date) + if localDate.tzinfo == None: + localDate = tz.localize(localDate) + self.localDate = localDate + self.date = self.localDate.astimezone(pytz.UTC) # Update atmospheric conditions if atmosphere type is Forecast, # Reanalysis or Ensemble @@ -478,7 +516,7 @@ def setElevation(self, elevation="Open-Elevation"): response = requests.get(requestURL) results = response.json()["results"] self.elevation = results[0]["elevation"] - print("Elevation received: ", self.elevation) + print("Elevation received:", self.elevation) except: raise RuntimeError("Unabel to reach Open-Elevation API servers.") else: @@ -487,6 +525,7 @@ def setElevation(self, elevation="Open-Elevation"): " Open-Elevation API. See Environment.setLocation." ) + @requires_netCDF4 def setTopographicProfile(self, type, file, dictionary="netCDF4", crs=None): """[UNDER CONSTRUCTION] Defines the Topographic profile, importing data from previous downloaded files. Mainly data from the Shuttle Radar @@ -515,7 +554,7 @@ def setTopographicProfile(self, type, file, dictionary="netCDF4", crs=None): if type == "NASADEM_HGT": if dictionary == "netCDF4": - rootgrp = Dataset(file, "r", format="NETCDF4") + rootgrp = netCDF4.Dataset(file, "r", format="NETCDF4") self.elevLonArray = rootgrp.variables["lon"][:].tolist() self.elevLatArray = rootgrp.variables["lat"][:].tolist() self.elevArray = rootgrp.variables["NASADEM_HGT"][:].tolist() @@ -1686,6 +1725,7 @@ def processNOAARUCSounding(self, file): # Save maximum expected height self.maxExpectedHeight = pressure_array[-1, 0] + @requires_netCDF4 def processForecastReanalysis(self, file, dictionary): """Import and process atmospheric data from weather forecasts and reanalysis given as netCDF or OPeNDAP files. @@ -1944,7 +1984,7 @@ def processForecastReanalysis(self, file, dictionary): windV = ((y2 - y) / (y2 - y1)) * f_x_y1 + ((y - y1) / (y2 - y1)) * f_x_y2 # Determine wind speed, heading and direction - windSpeed = np.sqrt(windU ** 2 + windV ** 2) + windSpeed = np.sqrt(windU**2 + windV**2) windHeading = np.arctan2(windU, windV) * (180 / np.pi) % 360 windDirection = (windHeading - 180) % 360 @@ -2079,6 +2119,7 @@ def processForecastReanalysis(self, file, dictionary): return None + @requires_netCDF4 def processEnsemble(self, file, dictionary): """Import and process atmospheric data from weather ensembles given as netCDF or OPeNDAP files. @@ -2361,7 +2402,7 @@ def processEnsemble(self, file, dictionary): windV = ((y2 - y) / (y2 - y1)) * f_x_y1 + ((y - y1) / (y2 - y1)) * f_x_y2 # Determine wind speed, heading and direction - windSpeed = np.sqrt(windU ** 2 + windV ** 2) + windSpeed = np.sqrt(windU**2 + windV**2) windHeading = np.arctan2(windU, windV) * (180 / np.pi) % 360 windDirection = (windHeading - 180) % 360 @@ -2820,9 +2861,18 @@ def info(self): """ # Print launch site details print("Launch Site Details") - print("\nLaunch Rail Length: ", self.rL, " m") - if self.date != None: - print("Launch Date: ", self.date, " UTC") + print("\nLaunch Rail Length:", self.rL, " m") + time_format = "%Y-%m-%d %H:%M:%S" + if self.date != None and "UTC" not in self.timeZone: + print( + "Launch Date:", + self.date.strftime(time_format), + "UTC |", + self.localDate.strftime(time_format), + self.timeZone, + ) + elif self.date != None: + print("Launch Date:", self.date.strftime(time_format), "UTC") if self.lat != None and self.lon != None: print("Launch Site Latitude: {:.5f}°".format(self.lat)) print("Launch Site Longitude: {:.5f}°".format(self.lon)) @@ -2839,7 +2889,7 @@ def info(self): # Print atmospheric model details print("\n\nAtmospheric Model Details") modelType = self.atmosphericModelType - print("\nAtmospheric Model Type: ", modelType) + print("\nAtmospheric Model Type:", modelType) print( modelType + " Maximum Height: {:.3f} km".format(self.maxExpectedHeight / 1000) @@ -2850,7 +2900,7 @@ def info(self): endDate = self.atmosphericModelEndDate interval = self.atmosphericModelInterval print(modelType + " Time Period: From ", initDate, " to ", endDate, " UTC") - print(modelType + " Hour Interval: ", interval, " hrs") + print(modelType + " Hour Interval:", interval, " hrs") # Determine latitude and longitude range initLat = self.atmosphericModelInitLat endLat = self.atmosphericModelEndLat @@ -2859,8 +2909,8 @@ def info(self): print(modelType + " Latitude Range: From ", initLat, "° To ", endLat, "°") print(modelType + " Longitude Range: From ", initLon, "° To ", endLon, "°") if modelType == "Ensemble": - print("Number of Ensemble Members: ", self.numEnsembleMembers) - print("Selected Ensemble Member: ", self.ensembleMember, " (Starts from 0)") + print("Number of Ensemble Members:", self.numEnsembleMembers) + print("Selected Ensemble Member:", self.ensembleMember, " (Starts from 0)") # Print atmospheric conditions print("\n\nSurface Atmospheric Conditions") @@ -2947,9 +2997,18 @@ def allInfo(self): # Print launch site details print("\n\nLaunch Site Details") - print("\nLaunch Rail Length: ", self.rL, " m") - if self.date != None: - print("Launch Date: ", self.date, " UTC") + print("\nLaunch Rail Length:", self.rL, " m") + time_format = "%Y-%m-%d %H:%M:%S" + if self.date != None and "UTC" not in self.timeZone: + print( + "Launch Date:", + self.date.strftime(time_format), + "UTC |", + self.localDate.strftime(time_format), + self.timeZone, + ) + elif self.date != None: + print("Launch Date:", self.date.strftime(time_format), "UTC") if self.lat != None and self.lon != None: print("Launch Site Latitude: {:.5f}°".format(self.lat)) print("Launch Site Longitude: {:.5f}°".format(self.lon)) @@ -2958,7 +3017,7 @@ def allInfo(self): # Print atmospheric model details print("\n\nAtmospheric Model Details") modelType = self.atmosphericModelType - print("\nAtmospheric Model Type: ", modelType) + print("\nAtmospheric Model Type:", modelType) print( modelType + " Maximum Height: {:.3f} km".format(self.maxExpectedHeight / 1000) @@ -2969,7 +3028,7 @@ def allInfo(self): endDate = self.atmosphericModelEndDate interval = self.atmosphericModelInterval print(modelType + " Time Period: From ", initDate, " to ", endDate, " UTC") - print(modelType + " Hour Interval: ", interval, " hrs") + print(modelType + " Hour Interval:", interval, " hrs") # Determine latitude and longitude range initLat = self.atmosphericModelInitLat endLat = self.atmosphericModelEndLat @@ -2978,8 +3037,8 @@ def allInfo(self): print(modelType + " Latitude Range: From ", initLat, "° To ", endLat, "°") print(modelType + " Longitude Range: From ", initLon, "° To ", endLon, "°") if modelType == "Ensemble": - print("Number of Ensemble Members: ", self.numEnsembleMembers) - print("Selected Ensemble Member: ", self.ensembleMember, " (Starts from 0)") + print("Number of Ensemble Members:", self.numEnsembleMembers) + print("Selected Ensemble Member:", self.ensembleMember, " (Starts from 0)") # Print atmospheric conditions print("\n\nSurface Atmospheric Conditions") @@ -3249,7 +3308,7 @@ def geodesicToUtm(self, lat, lon, datum): # Evaluate reference parameters K0 = 1 - 1 / 2500 - e2 = 2 * flattening - flattening ** 2 + e2 = 2 * flattening - flattening**2 e2lin = e2 / (1 - e2) # Evaluate auxiliary parameters @@ -3272,9 +3331,9 @@ def geodesicToUtm(self, lat, lon, datum): # 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 + 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 + M = (61 - 58 * t + t * t + 600 * c - 330 * e2lin) * (ag**6) / 720 # Evaluate the final coordinates x = 500000 + K0 * n * (ag + J + K) @@ -3347,7 +3406,7 @@ def utmToGeodesic(self, x, y, utmZone, hemis, datum): # Calculate reference values K0 = 1 - 1 / 2500 - e2 = 2 * flattening - flattening ** 2 + e2 = 2 * flattening - flattening**2 e2lin = e2 / (1 - e2) e1 = (1 - (1 - e2) ** 0.5) / (1 + (1 - e2) ** 0.5) @@ -3371,20 +3430,20 @@ def utmToGeodesic(self, x, y, utmZone, hemis, datum): 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) + 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) + * (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) + * (d**5) / 120 ) @@ -3447,8 +3506,8 @@ def calculateEarthRadius(self, lat, datum): # 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)) ** 2 + + (np.sin(lat) * (semiMinorAxis**2)) ** 2 ) / ((np.cos(lat) * semiMajorAxis) ** 2 + (np.sin(lat) * semiMinorAxis) ** 2) ) @@ -3506,7 +3565,7 @@ def printEarthDetails(self): # print("Launch Site UTM coordinates: {:.2f} ".format(self.initialEast) # + self.initialEW + " {:.2f} ".format(self.initialNorth) + self.initialHemisphere # ) - # print("Launch Site UTM zone number: ", self.initialUtmZone) + # print("Launch Site UTM zone number:", self.initialUtmZone) # print("Launch Site Surface Elevation: {:.1f} m".format(self.elevation)) print("Earth Radius at Launch site: {:.1f} m".format(self.earthRadius)) print("Gravity acceleration at launch site: Still not implemented :(") diff --git a/rocketpy/Flight.py b/rocketpy/Flight.py index 2ac4d1a7e..3824af5ab 100644 --- a/rocketpy/Flight.py +++ b/rocketpy/Flight.py @@ -833,7 +833,7 @@ def __init__( self.y[0] ** 2 + self.y[1] ** 2 + (self.y[2] - self.env.elevation) ** 2 - >= self.effective1RL ** 2 + >= self.effective1RL**2 ): # Rocket is out of rail # Check exactly when it went out using root finding @@ -847,7 +847,7 @@ def __init__( # Get points y0 = ( sum([self.solution[-2][i] ** 2 for i in [1, 2, 3]]) - - self.effective1RL ** 2 + - self.effective1RL**2 ) yp0 = 2 * sum( [ @@ -858,7 +858,7 @@ def __init__( t1 = self.solution[-1][0] - self.solution[-2][0] y1 = ( sum([self.solution[-1][i] ** 2 for i in [1, 2, 3]]) - - self.effective1RL ** 2 + - self.effective1RL**2 ) yp1 = 2 * sum( [ @@ -873,15 +873,15 @@ def __init__( D = float(phase.solver.step_size) d = float(y0) c = float(yp0) - b = float((3 * y1 - yp1 * D - 2 * c * D - 3 * d) / (D ** 2)) - a = float(-(2 * y1 - yp1 * D - c * D - 2 * d) / (D ** 3)) + 1e-5 + b = float((3 * y1 - yp1 * D - 2 * c * D - 3 * d) / (D**2)) + a = float(-(2 * y1 - yp1 * D - c * D - 2 * d) / (D**3)) + 1e-5 # Find roots - d0 = b ** 2 - 3 * a * c - d1 = 2 * b ** 3 - 9 * a * b * c + 27 * d * a ** 2 - c1 = ((d1 + (d1 ** 2 - 4 * d0 ** 3) ** (0.5)) / 2) ** (1 / 3) + d0 = b**2 - 3 * a * c + d1 = 2 * b**3 - 9 * a * b * c + 27 * d * a**2 + c1 = ((d1 + (d1**2 - 4 * d0**3) ** (0.5)) / 2) ** (1 / 3) t_roots = [] for k in [0, 1, 2]: - c2 = c1 * (-1 / 2 + 1j * (3 ** 0.5) / 2) ** k + c2 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** k t_roots.append(-(1 / (3 * a)) * (b + c2 + d0 / c2)) # Find correct root valid_t_root = [] @@ -968,15 +968,15 @@ def __init__( D = float(phase.solver.step_size) d = float(y0) c = float(yp0) - b = float((3 * y1 - yp1 * D - 2 * c * D - 3 * d) / (D ** 2)) - a = float(-(2 * y1 - yp1 * D - c * D - 2 * d) / (D ** 3)) + b = float((3 * y1 - yp1 * D - 2 * c * D - 3 * d) / (D**2)) + a = float(-(2 * y1 - yp1 * D - c * D - 2 * d) / (D**3)) # Find roots - d0 = b ** 2 - 3 * a * c - d1 = 2 * b ** 3 - 9 * a * b * c + 27 * d * a ** 2 - c1 = ((d1 + (d1 ** 2 - 4 * d0 ** 3) ** (0.5)) / 2) ** (1 / 3) + d0 = b**2 - 3 * a * c + d1 = 2 * b**3 - 9 * a * b * c + 27 * d * a**2 + c1 = ((d1 + (d1**2 - 4 * d0**3) ** (0.5)) / 2) ** (1 / 3) t_roots = [] for k in [0, 1, 2]: - c2 = c1 * (-1 / 2 + 1j * (3 ** 0.5) / 2) ** k + c2 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** k t_roots.append(-(1 / (3 * a)) * (b + c2 + d0 / c2)) # Find correct root valid_t_root = [] @@ -1151,14 +1151,14 @@ def uDotRail1(self, t, u, postProcessing=False): # Calculate Forces Thrust = self.rocket.motor.thrust.getValueOpt(t) rho = self.env.density.getValueOpt(z) - R3 = -0.5 * rho * (freestreamSpeed ** 2) * self.rocket.area * (dragCoeff) + R3 = -0.5 * rho * (freestreamSpeed**2) * self.rocket.area * (dragCoeff) # Calculate Linear acceleration - a3 = (R3 + Thrust) / M - (e0 ** 2 - e1 ** 2 - e2 ** 2 + e3 ** 2) * self.env.g + a3 = (R3 + Thrust) / M - (e0**2 - e1**2 - e2**2 + e3**2) * self.env.g if a3 > 0: ax = 2 * (e1 * e3 + e0 * e2) * a3 ay = 2 * (e2 * e3 - e0 * e1) * a3 - az = (1 - 2 * (e1 ** 2 + e2 ** 2)) * a3 + az = (1 - 2 * (e1**2 + e2**2)) * a3 else: ax, ay, az = 0, 0, 0 @@ -1257,15 +1257,15 @@ def uDot(self, t, u, postProcessing=False): a = b * Mt / M rN = self.rocket.motor.nozzleRadius # Prepare transformation matrix - a11 = 1 - 2 * (e2 ** 2 + e3 ** 2) + a11 = 1 - 2 * (e2**2 + e3**2) a12 = 2 * (e1 * e2 - e0 * e3) a13 = 2 * (e1 * e3 + e0 * e2) a21 = 2 * (e1 * e2 + e0 * e3) - a22 = 1 - 2 * (e1 ** 2 + e3 ** 2) + a22 = 1 - 2 * (e1**2 + e3**2) a23 = 2 * (e2 * e3 - e0 * e1) a31 = 2 * (e1 * e3 - e0 * e2) a32 = 2 * (e2 * e3 + e0 * e1) - a33 = 1 - 2 * (e1 ** 2 + e2 ** 2) + a33 = 1 - 2 * (e1**2 + e2**2) # Transformation matrix: (123) -> (XYZ) K = [[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]] # Transformation matrix: (XYZ) -> (123) or K transpose @@ -1287,7 +1287,7 @@ def uDot(self, t, u, postProcessing=False): else: dragCoeff = self.rocket.powerOffDrag.getValueOpt(freestreamMach) rho = self.env.density.getValueOpt(z) - R3 = -0.5 * rho * (freestreamSpeed ** 2) * self.rocket.area * (dragCoeff) + R3 = -0.5 * rho * (freestreamSpeed**2) * self.rocket.area * (dragCoeff) # Off center moment M1 += self.rocket.cpEccentricityY * R3 M2 -= self.rocket.cpEccentricityX * R3 @@ -1297,7 +1297,7 @@ def uDot(self, t, u, postProcessing=False): vzB = a13 * vx + a23 * vy + a33 * vz # Calculate lift and moment for each component of the rocket for aerodynamicSurface in self.rocket.aerodynamicSurfaces: - compCp = aerodynamicSurface[0][2] + compCp = aerodynamicSurface["cp"][2] # Component absolute velocity in body frame compVxB = vxB + compCp * omega2 compVyB = vyB - compCp * omega1 @@ -1314,23 +1314,23 @@ def uDot(self, t, u, postProcessing=False): compStreamVyB = compWindVyB - compVyB compStreamVzB = compWindVzB - compVzB compStreamSpeed = ( - compStreamVxB ** 2 + compStreamVyB ** 2 + compStreamVzB ** 2 + compStreamVxB**2 + compStreamVyB**2 + compStreamVzB**2 ) ** 0.5 # Component attack angle and lift force compAttackAngle = 0 compLift, compLiftXB, compLiftYB = 0, 0, 0 - if compStreamVxB ** 2 + compStreamVyB ** 2 != 0: + if compStreamVxB**2 + compStreamVyB**2 != 0: # Normalize component stream velocity in body frame compStreamVzBn = compStreamVzB / compStreamSpeed if -1 * compStreamVzBn < 1: compAttackAngle = np.arccos(-compStreamVzBn) - cLift = abs(aerodynamicSurface[1](compAttackAngle)) + cLift = aerodynamicSurface["cl"](compAttackAngle, freestreamMach) # Component lift force magnitude compLift = ( - 0.5 * rho * (compStreamSpeed ** 2) * self.rocket.area * cLift + 0.5 * rho * (compStreamSpeed**2) * self.rocket.area * cLift ) # Component lift force components - liftDirNorm = (compStreamVxB ** 2 + compStreamVyB ** 2) ** 0.5 + liftDirNorm = (compStreamVxB**2 + compStreamVyB**2) ** 0.5 compLiftXB = compLift * (compStreamVxB / liftDirNorm) compLiftYB = compLift * (compStreamVyB / liftDirNorm) # Add to total lift force @@ -1340,37 +1340,50 @@ def uDot(self, t, u, postProcessing=False): M1 -= (compCp + a) * compLiftYB M2 += (compCp + a) * compLiftXB # Calculates Roll Moment - if aerodynamicSurface[-1] == "Fins": - Clfdelta, Cldomega, cantAngleRad = aerodynamicSurface[2] - if cantAngleRad != 0: - Clf = Clfdelta * cantAngleRad - Cld = Cldomega * omega3 / freestreamSpeed - M3 += Clf - Cld + if aerodynamicSurface["name"] == "Fins": + Clfdelta, Cldomega, cantAngleRad = aerodynamicSurface["roll parameters"] + M3f = ( + (1 / 2 * rho * freestreamSpeed**2) + * self.rocket.area + * 2 + * self.rocket.radius + * Clfdelta(freestreamMach) + * cantAngleRad + ) + M3d = ( + (1 / 2 * rho * freestreamSpeed) + * self.rocket.area + * (2 * self.rocket.radius) ** 2 + * Cldomega(freestreamMach) + * omega3 + / 2 + ) + M3 += M3f - M3d # Calculate derivatives # Angular acceleration alpha1 = ( M1 - ( - omega2 * omega3 * (Rz + Tz - Ri - Ti - mu * b ** 2) + omega2 * omega3 * (Rz + Tz - Ri - Ti - mu * b**2) + omega1 * ( (TiDot + MtDot * (Mr - 1) * (b / M) ** 2) - MtDot * ((rN / 2) ** 2 + (c - b * mu / Mr) ** 2) ) ) - ) / (Ri + Ti + mu * b ** 2) + ) / (Ri + Ti + mu * b**2) alpha2 = ( M2 - ( - omega1 * omega3 * (Ri + Ti + mu * b ** 2 - Rz - Tz) + omega1 * omega3 * (Ri + Ti + mu * b**2 - Rz - Tz) + omega2 * ( (TiDot + MtDot * (Mr - 1) * (b / M) ** 2) - MtDot * ((rN / 2) ** 2 + (c - b * mu / Mr) ** 2) ) ) - ) / (Ri + Ti + mu * b ** 2) - alpha3 = (M3 - omega3 * (TzDot - MtDot * (rN ** 2) / 2)) / (Rz + Tz) + ) / (Ri + Ti + mu * b**2) + alpha3 = (M3 - omega3 * (TzDot - MtDot * (rN**2) / 2)) / (Rz + Tz) # Euler parameters derivative e0Dot = 0.5 * (-omega1 * e1 - omega2 * e2 - omega3 * e3) e1Dot = 0.5 * (omega1 * e0 + omega3 * e2 - omega2 * e3) @@ -1379,7 +1392,7 @@ def uDot(self, t, u, postProcessing=False): # Linear acceleration L = [ - (R1 - b * Mt * (omega2 ** 2 + omega3 ** 2) - 2 * c * MtDot * omega2) / M, + (R1 - b * Mt * (omega2**2 + omega3**2) - 2 * c * MtDot * omega2) / M, (R2 + b * Mt * (alpha3 + omega1 * omega2) + 2 * c * MtDot * omega1) / M, (R3 - b * Mt * (alpha2 - omega1 * omega3) + Thrust) / M, ] @@ -1450,11 +1463,11 @@ def uDotParachute(self, t, u, postProcessing=False): R = 1.5 rho = self.env.density.getValueOpt(u[2]) to = 1.2 - ma = ka * rho * (4 / 3) * np.pi * R ** 3 + ma = ka * rho * (4 / 3) * np.pi * R**3 mp = self.rocket.mass eta = 1 - Rdot = (6 * R * (1 - eta) / (1.2 ** 6)) * ( - (1 - eta) * t ** 5 + eta * (to ** 3) * (t ** 2) + Rdot = (6 * R * (1 - eta) / (1.2**6)) * ( + (1 - eta) * t**5 + eta * (to**3) * (t**2) ) Rdot = 0 # Get relevant state data @@ -1470,7 +1483,7 @@ def uDotParachute(self, t, u, postProcessing=False): freestreamZ = vz # Determine drag force pseudoD = ( - -0.5 * rho * CdS * freestreamSpeed - ka * rho * 4 * np.pi * (R ** 2) * Rdot + -0.5 * rho * CdS * freestreamSpeed - ka * rho * 4 * np.pi * (R**2) * Rdot ) Dx = pseudoD * freestreamX Dy = pseudoD * freestreamY @@ -1497,7 +1510,7 @@ def uDotParachute(self, t, u, postProcessing=False): return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0] - def postProcess(self): + def postProcess(self, interpolation="spline", extrapolation="natural"): """Post-process all Flight information produced during simulation. Includes the calculation of maximum values, calculation of secondary values such as energy and conversion @@ -1515,43 +1528,43 @@ def postProcess(self): # Transform solution array into Functions sol = np.array(self.solution) self.x = Function( - sol[:, [0, 1]], "Time (s)", "X (m)", "spline", extrapolation="natural" + sol[:, [0, 1]], "Time (s)", "X (m)", interpolation, extrapolation ) self.y = Function( - sol[:, [0, 2]], "Time (s)", "Y (m)", "spline", extrapolation="natural" + sol[:, [0, 2]], "Time (s)", "Y (m)", interpolation, extrapolation ) self.z = Function( - sol[:, [0, 3]], "Time (s)", "Z (m)", "spline", extrapolation="natural" + sol[:, [0, 3]], "Time (s)", "Z (m)", interpolation, extrapolation ) self.vx = Function( - sol[:, [0, 4]], "Time (s)", "Vx (m/s)", "spline", extrapolation="natural" + sol[:, [0, 4]], "Time (s)", "Vx (m/s)", interpolation, extrapolation ) self.vy = Function( - sol[:, [0, 5]], "Time (s)", "Vy (m/s)", "spline", extrapolation="natural" + sol[:, [0, 5]], "Time (s)", "Vy (m/s)", interpolation, extrapolation ) self.vz = Function( - sol[:, [0, 6]], "Time (s)", "Vz (m/s)", "spline", extrapolation="natural" + sol[:, [0, 6]], "Time (s)", "Vz (m/s)", interpolation, extrapolation ) self.e0 = Function( - sol[:, [0, 7]], "Time (s)", "e0", "spline", extrapolation="natural" + sol[:, [0, 7]], "Time (s)", "e0", interpolation, extrapolation ) self.e1 = Function( - sol[:, [0, 8]], "Time (s)", "e1", "spline", extrapolation="natural" + sol[:, [0, 8]], "Time (s)", "e1", interpolation, extrapolation ) self.e2 = Function( - sol[:, [0, 9]], "Time (s)", "e2", "spline", extrapolation="natural" + sol[:, [0, 9]], "Time (s)", "e2", interpolation, extrapolation ) self.e3 = Function( - sol[:, [0, 10]], "Time (s)", "e3", "spline", extrapolation="natural" + sol[:, [0, 10]], "Time (s)", "e3", interpolation, extrapolation ) self.w1 = Function( - sol[:, [0, 11]], "Time (s)", "ω1 (rad/s)", "spline", extrapolation="natural" + sol[:, [0, 11]], "Time (s)", "ω1 (rad/s)", interpolation, extrapolation ) self.w2 = Function( - sol[:, [0, 12]], "Time (s)", "ω2 (rad/s)", "spline", extrapolation="natural" + sol[:, [0, 12]], "Time (s)", "ω2 (rad/s)", interpolation, extrapolation ) self.w3 = Function( - sol[:, [0, 13]], "Time (s)", "ω3 (rad/s)", "spline", extrapolation="natural" + sol[:, [0, 13]], "Time (s)", "ω3 (rad/s)", interpolation, extrapolation ) # Process second type of outputs - accelerations @@ -1583,12 +1596,12 @@ def postProcess(self): self.alpha2.append([step[0], alpha2]) self.alpha3.append([step[0], alpha3]) # Convert accelerations to functions - self.ax = Function(self.ax, "Time (s)", "Ax (m/s2)", "spline") - self.ay = Function(self.ay, "Time (s)", "Ay (m/s2)", "spline") - self.az = Function(self.az, "Time (s)", "Az (m/s2)", "spline") - self.alpha1 = Function(self.alpha1, "Time (s)", "α1 (rad/s2)", "spline") - self.alpha2 = Function(self.alpha2, "Time (s)", "α2 (rad/s2)", "spline") - self.alpha3 = Function(self.alpha3, "Time (s)", "α3 (rad/s2)", "spline") + self.ax = Function(self.ax, "Time (s)", "Ax (m/s2)", interpolation) + self.ay = Function(self.ay, "Time (s)", "Ay (m/s2)", interpolation) + self.az = Function(self.az, "Time (s)", "Az (m/s2)", interpolation) + self.alpha1 = Function(self.alpha1, "Time (s)", "α1 (rad/s2)", interpolation) + self.alpha2 = Function(self.alpha2, "Time (s)", "α2 (rad/s2)", interpolation) + self.alpha3 = Function(self.alpha3, "Time (s)", "α3 (rad/s2)", interpolation) # Process third type of outputs - temporary values calculated during integration # Initialize force and atmospheric arrays @@ -1615,44 +1628,54 @@ def postProcess(self): # Call derivatives in post processing mode uDot = currentDerivative(step[0], step[1:], postProcessing=True) # Convert forces and atmospheric arrays to functions - self.R1 = Function(self.R1, "Time (s)", "R1 (N)", "spline") - self.R2 = Function(self.R2, "Time (s)", "R2 (N)", "spline") - self.R3 = Function(self.R3, "Time (s)", "R3 (N)", "spline") - self.M1 = Function(self.M1, "Time (s)", "M1 (Nm)", "spline") - self.M2 = Function(self.M2, "Time (s)", "M2 (Nm)", "spline") - self.M3 = Function(self.M3, "Time (s)", "M3 (Nm)", "spline") + self.R1 = Function(self.R1, "Time (s)", "R1 (N)", interpolation) + self.R2 = Function(self.R2, "Time (s)", "R2 (N)", interpolation) + self.R3 = Function(self.R3, "Time (s)", "R3 (N)", interpolation) + self.M1 = Function(self.M1, "Time (s)", "M1 (Nm)", interpolation) + self.M2 = Function(self.M2, "Time (s)", "M2 (Nm)", interpolation) + self.M3 = Function(self.M3, "Time (s)", "M3 (Nm)", interpolation) self.windVelocityX = Function( - self.windVelocityX, "Time (s)", "Wind Velocity X (East) (m/s)", "spline" + self.windVelocityX, + "Time (s)", + "Wind Velocity X (East) (m/s)", + interpolation, ) self.windVelocityY = Function( - self.windVelocityY, "Time (s)", "Wind Velocity Y (North) (m/s)", "spline" + self.windVelocityY, + "Time (s)", + "Wind Velocity Y (North) (m/s)", + interpolation, + ) + self.density = Function( + self.density, "Time (s)", "Density (kg/m³)", interpolation + ) + self.pressure = Function( + self.pressure, "Time (s)", "Pressure (Pa)", interpolation ) - self.density = Function(self.density, "Time (s)", "Density (kg/m³)", "spline") - self.pressure = Function(self.pressure, "Time (s)", "Pressure (Pa)", "spline") self.dynamicViscosity = Function( - self.dynamicViscosity, "Time (s)", "Dynamic Viscosity (Pa s)", "spline" + self.dynamicViscosity, "Time (s)", "Dynamic Viscosity (Pa s)", interpolation ) self.speedOfSound = Function( - self.speedOfSound, "Time (s)", "Speed of Sound (m/s)", "spline" + self.speedOfSound, "Time (s)", "Speed of Sound (m/s)", interpolation ) # Process fourth type of output - values calculated from previous outputs # Kinematics functions and values # Velocity Magnitude - self.speed = (self.vx ** 2 + self.vy ** 2 + self.vz ** 2) ** 0.5 + self.speed = (self.vx**2 + self.vy**2 + self.vz**2) ** 0.5 self.speed.setOutputs("Speed - Velocity Magnitude (m/s)") maxSpeedTimeIndex = np.argmax(self.speed[:, 1]) self.maxSpeed = self.speed[maxSpeedTimeIndex, 1] self.maxSpeedTime = self.speed[maxSpeedTimeIndex, 0] # Acceleration - self.acceleration = (self.ax ** 2 + self.ay ** 2 + self.az ** 2) ** 0.5 + self.acceleration = (self.ax**2 + self.ay**2 + self.az**2) ** 0.5 self.acceleration.setOutputs("Acceleration Magnitude (m/s²)") maxAccelerationTimeIndex = np.argmax(self.acceleration[:, 1]) self.maxAcceleration = self.acceleration[maxAccelerationTimeIndex, 1] self.maxAccelerationTime = self.acceleration[maxAccelerationTimeIndex, 0] # Path Angle - self.horizontalSpeed = (self.vx ** 2 + self.vy ** 2) ** 0.5 + self.horizontalSpeed = (self.vx**2 + self.vy**2) ** 0.5 pathAngle = (180 / np.pi) * np.arctan2( self.vz[:, 1], self.horizontalSpeed[:, 1] ) @@ -1661,9 +1684,9 @@ def postProcess(self): # Attitude Angle self.attitudeVectorX = 2 * (self.e1 * self.e3 + self.e0 * self.e2) # a13 self.attitudeVectorY = 2 * (self.e2 * self.e3 - self.e0 * self.e1) # a23 - self.attitudeVectorZ = 1 - 2 * (self.e1 ** 2 + self.e2 ** 2) # a33 + self.attitudeVectorZ = 1 - 2 * (self.e1**2 + self.e2**2) # a33 horizontalAttitudeProj = ( - self.attitudeVectorX ** 2 + self.attitudeVectorY ** 2 + self.attitudeVectorX**2 + self.attitudeVectorY**2 ) ** 0.5 attitudeAngle = (180 / np.pi) * np.arctan2( self.attitudeVectorZ[:, 1], horizontalAttitudeProj[:, 1] @@ -1684,9 +1707,9 @@ def postProcess(self): attitudeLateralPlaneProjY = self.attitudeVectorY[:, 1] - attitudeLateralProjY attitudeLateralPlaneProjZ = self.attitudeVectorZ[:, 1] attitudeLateralPlaneProj = ( - attitudeLateralPlaneProjX ** 2 - + attitudeLateralPlaneProjY ** 2 - + attitudeLateralPlaneProjZ ** 2 + attitudeLateralPlaneProjX**2 + + attitudeLateralPlaneProjY**2 + + attitudeLateralPlaneProjZ**2 ) ** 0.5 lateralAttitudeAngle = (180 / np.pi) * np.arctan2( attitudeLateralProj, attitudeLateralPlaneProj @@ -1778,11 +1801,11 @@ def postProcess(self): self.railButton2ShearForce[:outOfRailTimeIndex] ) # Aerodynamic Lift and Drag - self.aerodynamicLift = (self.R1 ** 2 + self.R2 ** 2) ** 0.5 + self.aerodynamicLift = (self.R1**2 + self.R2**2) ** 0.5 self.aerodynamicLift.setOutputs("Aerodynamic Lift Force (N)") self.aerodynamicDrag = -1 * self.R3 self.aerodynamicDrag.setOutputs("Aerodynamic Drag Force (N)") - self.aerodynamicBendingMoment = (self.M1 ** 2 + self.M2 ** 2) ** 0.5 + self.aerodynamicBendingMoment = (self.M1**2 + self.M2**2) ** 0.5 self.aerodynamicBendingMoment.setOutputs("Aerodynamic Bending Moment (N m)") self.aerodynamicSpinMoment = self.M3 self.aerodynamicSpinMoment.setOutputs("Aerodynamic Spin Moment (N m)") @@ -1794,7 +1817,7 @@ def postProcess(self): Ri = self.rocket.inertiaI Tz = self.rocket.motor.inertiaZ Ti = self.rocket.motor.inertiaI - I1, I2, I3 = (Ri + Ti + mu * b ** 2), (Ri + Ti + mu * b ** 2), (Rz + Tz) + I1, I2, I3 = (Ri + Ti + mu * b**2), (Ri + Ti + mu * b**2), (Rz + Tz) # Redefine I1, I2 and I3 grid grid = self.vx[:, 0] I1 = Function(np.column_stack([grid, I1(grid)]), "Time (s)") @@ -1810,9 +1833,9 @@ def postProcess(self): vx, vy, vz = self.vx, self.vy, self.vz w1, w2, w3 = self.w1, self.w2, self.w3 # Kinetic Energy - self.rotationalEnergy = 0.5 * (I1 * w1 ** 2 + I2 * w2 ** 2 + I3 * w3 ** 2) + self.rotationalEnergy = 0.5 * (I1 * w1**2 + I2 * w2**2 + I3 * w3**2) self.rotationalEnergy.setOutputs("Rotational Kinetic Energy (J)") - self.translationalEnergy = 0.5 * totalMass * (vx ** 2 + vy ** 2 + vz ** 2) + self.translationalEnergy = 0.5 * totalMass * (vx**2 + vy**2 + vz**2) self.translationalEnergy.setOutputs("Translational Kinetic Energy (J)") self.kineticEnergy = self.rotationalEnergy + self.translationalEnergy self.kineticEnergy.setOutputs("Kinetic Energy (J)") @@ -1921,9 +1944,9 @@ def postProcess(self): self.streamVelocityZ = -1 * self.vz self.streamVelocityZ.setOutputs("Freestream Velocity Z (m/s)") self.freestreamSpeed = ( - self.streamVelocityX ** 2 - + self.streamVelocityY ** 2 - + self.streamVelocityZ ** 2 + self.streamVelocityX**2 + + self.streamVelocityY**2 + + self.streamVelocityZ**2 ) ** 0.5 self.freestreamSpeed.setOutputs("Freestream Speed (m/s)") # Apogee Freestream speed @@ -1943,7 +1966,7 @@ def postProcess(self): self.maxReynoldsNumberTime = self.ReynoldsNumber[maxReynoldsNumberTimeIndex, 0] self.maxReynoldsNumber = self.ReynoldsNumber[maxReynoldsNumberTimeIndex, 1] # Dynamic Pressure - self.dynamicPressure = 0.5 * self.density * self.freestreamSpeed ** 2 + self.dynamicPressure = 0.5 * self.density * self.freestreamSpeed**2 self.dynamicPressure.setOutputs("Dynamic Pressure (Pa)") maxDynamicPressureTimeIndex = np.argmax(self.dynamicPressure[:, 1]) self.maxDynamicPressureTime = self.dynamicPressure[ @@ -1951,7 +1974,7 @@ def postProcess(self): ] self.maxDynamicPressure = self.dynamicPressure[maxDynamicPressureTimeIndex, 1] # Total Pressure - self.totalPressure = self.pressure * (1 + 0.2 * self.MachNumber ** 2) ** (3.5) + self.totalPressure = self.pressure * (1 + 0.2 * self.MachNumber**2) ** (3.5) self.totalPressure.setOutputs("Total Pressure (Pa)") maxtotalPressureTimeIndex = np.argmax(self.totalPressure[:, 1]) self.maxtotalPressureTime = self.totalPressure[maxtotalPressureTimeIndex, 0] @@ -2305,7 +2328,7 @@ def calculateStallWindVelocity(self, stallAngle): wV = ( 2 * vF * math.cos(theta) / c + ( - 4 * vF * vF * math.cos(theta) * math.cos(theta) / (c ** 2) + 4 * vF * vF * math.cos(theta) * math.cos(theta) / (c**2) + 4 * 1 * vF * vF / c ) ** 0.5 @@ -2996,7 +3019,7 @@ def calculateFinFlutterAnalysis(self, finThickness, shearModulus): # Calculate the Fin Flutter Mach Number self.flutterMachNumber = ( (shearModulus * 2 * (ar + 2) * (finThickness / self.rocket.rootChord) ** 3) - / (1.337 * (ar ** 3) * (la + 1) * self.pressure) + / (1.337 * (ar**3) * (la + 1) * self.pressure) ) ** 0.5 # Calculate difference between Fin Flutter Mach Number and the Rocket Speed diff --git a/rocketpy/Function.py b/rocketpy/Function.py index acdebd479..446dca183 100644 --- a/rocketpy/Function.py +++ b/rocketpy/Function.py @@ -300,7 +300,7 @@ def getValueOpt(x): xInterval = xInterval if xInterval != 0 else 1 a = coeffs[:, xInterval - 1] x = x - xData[xInterval - 1] - y = a[3] * x ** 3 + a[2] * x ** 2 + a[1] * x + a[0] + y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] else: # Extrapolate if extrapolation == 0: # Extrapolation == zero @@ -308,7 +308,7 @@ def getValueOpt(x): elif extrapolation == 1: # Extrapolation == natural a = coeffs[:, 0] if x < xmin else coeffs[:, -1] x = x - xData[0] if x < xmin else x - xData[-2] - y = a[3] * x ** 3 + a[2] * x ** 2 + a[1] * x + a[0] + y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] else: # Extrapolation is set to constant y = yData[0] if x < xmin else yData[-1] return y @@ -352,14 +352,14 @@ def getValueOpt(x): # Interpolate xInterval = xInterval if xInterval != 0 else 1 a = coeffs[4 * xInterval - 4 : 4 * xInterval] - y = a[3] * x ** 3 + a[2] * x ** 2 + a[1] * x + a[0] + y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] else: # Extrapolate if extrapolation == 0: # Extrapolation == zero y = 0 elif extrapolation == 1: # Extrapolation == natural a = coeffs[:4] if x < xmin else coeffs[-4:] - y = a[3] * x ** 3 + a[2] * x ** 2 + a[1] * x + a[0] + y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] else: # Extrapolation is set to constant y = yData[0] if x < xmin else yData[-1] return y @@ -375,7 +375,7 @@ def getValueOpt(x): # Interpolate y = 0 for i in range(len(coeffs)): - y += coeffs[i] * (x ** i) + y += coeffs[i] * (x**i) else: # Extrapolate if extrapolation == 0: # Extrapolation == zero @@ -383,7 +383,7 @@ def getValueOpt(x): elif extrapolation == 1: # Extrapolation == natural y = 0 for i in range(len(coeffs)): - y += coeffs[i] * (x ** i) + y += coeffs[i] * (x**i) else: # Extrapolation is set to constant y = yData[0] if x < xmin else yData[-1] return y @@ -576,7 +576,7 @@ def getValue(self, *args): coeffs = self.__polynomialCoefficients__ A = np.zeros((len(args[0]), coeffs.shape[0])) for i in range(coeffs.shape[0]): - A[:, i] = x ** i + A[:, i] = x**i ans = A.dot(coeffs).tolist() for i in range(len(x)): if not (xmin <= x[i] <= xmax): @@ -713,7 +713,7 @@ def getValueOpt_deprecated(self, *args): xInterval = xInterval if xInterval != 0 else 1 a = coeffs[:, xInterval - 1] x = x - xData[xInterval - 1] - y = a[3] * x ** 3 + a[2] * x ** 2 + a[1] * x + a[0] + y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] else: # Extrapolate if extrapolation == 0: # Extrapolation == zero @@ -721,7 +721,7 @@ def getValueOpt_deprecated(self, *args): elif extrapolation == 1: # Extrapolation == natural a = coeffs[:, 0] if x < xmin else coeffs[:, -1] x = x - xData[0] if x < xmin else x - xData[-2] - y = a[3] * x ** 3 + a[2] * x ** 2 + a[1] * x + a[0] + y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] else: # Extrapolation is set to constant y = yData[0] if x < xmin else yData[-1] return y @@ -757,14 +757,14 @@ def getValueOpt_deprecated(self, *args): # Interpolate xInterval = xInterval if xInterval != 0 else 1 a = coeffs[4 * xInterval - 4 : 4 * xInterval] - y = a[3] * x ** 3 + a[2] * x ** 2 + a[1] * x + a[0] + y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] else: # Extrapolate if extrapolation == 0: # Extrapolation == zero y = 0 elif extrapolation == 1: # Extrapolation == natural a = coeffs[:4] if x < xmin else coeffs[-4:] - y = a[3] * x ** 3 + a[2] * x ** 2 + a[1] * x + a[0] + y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] else: # Extrapolation is set to constant y = yData[0] if x < xmin else yData[-1] return y @@ -777,7 +777,7 @@ def getValueOpt_deprecated(self, *args): # Interpolate y = 0 for i in range(len(coeffs)): - y += coeffs[i] * (x ** i) + y += coeffs[i] * (x**i) else: # Extrapolate if extrapolation == 0: # Extrapolation == zero @@ -785,7 +785,7 @@ def getValueOpt_deprecated(self, *args): elif extrapolation == 1: # Extrapolation == natural y = 0 for i in range(len(coeffs)): - y += coeffs[i] * (x ** i) + y += coeffs[i] * (x**i) else: # Extrapolation is set to constant y = yData[0] if x < xmin else yData[-1] return y @@ -863,7 +863,7 @@ def getValueOpt2(self, *args): else: a = coeffs[:, xInterval - 1] x = x - xData[xInterval - 1] - x = a[3] * x ** 3 + a[2] * x ** 2 + a[1] * x + a[0] + x = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] else: # Extrapolate if self.__extrapolation__ == "zero": @@ -896,10 +896,10 @@ def getValueOpt2(self, *args): x = yData[xInterval] elif xmin < x < xmax: a = coeffs[4 * xInterval - 4 : 4 * xInterval] - x = a[3] * x ** 3 + a[2] * x ** 2 + a[1] * x + a[0] + x = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] elif self.__extrapolation__ == "natural": a = coeffs[:4] if x < xmin else coeffs[-4:] - x = a[3] * x ** 3 + a[2] * x ** 2 + a[1] * x + a[0] + x = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] else: # Extrapolate if self.__extrapolation__ == "zero": @@ -1353,7 +1353,7 @@ def __interpolatePolynomial__(self): # Create coefficient matrix1 A = np.zeros((degree + 1, degree + 1)) for i in range(degree + 1): - A[:, i] = x ** i + A[:, i] = x**i # Solve the system and store the resultant coefficients self.__polynomialCoefficients__ = np.linalg.solve(A, y) @@ -1408,10 +1408,10 @@ def __interpolateAkima__(self): dl, dr = d[i], d[i + 1] A = np.array( [ - [1, xl, xl ** 2, xl ** 3], - [1, xr, xr ** 2, xr ** 3], - [0, 1, 2 * xl, 3 * xl ** 2], - [0, 1, 2 * xr, 3 * xr ** 2], + [1, xl, xl**2, xl**3], + [1, xr, xr**2, xr**3], + [0, 1, 2 * xl, 3 * xl**2], + [0, 1, 2 * xr, 3 * xr**2], ] ) Y = np.array([yl, yr, dl, dr]).T @@ -1977,9 +1977,9 @@ def integral(self, a, b, numerical=False): c = coeffs[:, 0] subB = a - xData[0] # subA = 0 ans -= ( - (c[3] * subB ** 4) / 4 - + (c[2] * subB ** 3 / 3) - + (c[1] * subB ** 2 / 2) + (c[3] * subB**4) / 4 + + (c[2] * subB**3 / 3) + + (c[1] * subB**2 / 2) + c[0] * subB ) else: @@ -1995,9 +1995,9 @@ def integral(self, a, b, numerical=False): c = coeffs[:, i] subB = xData[i + 1] - xData[i] # subA = 0 ans += ( - (c[3] * subB ** 4) / 4 - + (c[2] * subB ** 3 / 3) - + (c[1] * subB ** 2 / 2) + (c[3] * subB**4) / 4 + + (c[2] * subB**3 / 3) + + (c[1] * subB**2 / 2) + c[0] * subB ) i += 1 @@ -2010,15 +2010,15 @@ def integral(self, a, b, numerical=False): subA = xData[-1] - xData[-2] subB = b - xData[-2] ans -= ( - (c[3] * subA ** 4) / 4 - + (c[2] * subA ** 3 / 3) - + (c[1] * subA ** 2 / 2) + (c[3] * subA**4) / 4 + + (c[2] * subA**3 / 3) + + (c[1] * subA**2 / 2) + c[0] * subA ) ans += ( - (c[3] * subB ** 4) / 4 - + (c[2] * subB ** 3 / 3) - + (c[1] * subB ** 2 / 2) + (c[3] * subB**4) / 4 + + (c[2] * subB**3 / 3) + + (c[1] * subB**2 / 2) + c[0] * subB ) else: diff --git a/rocketpy/Rocket.py b/rocketpy/Rocket.py index b42ae84c2..585184201 100644 --- a/rocketpy/Rocket.py +++ b/rocketpy/Rocket.py @@ -220,7 +220,7 @@ def __init__( # Define rocket geometrical parameters in SI units self.radius = radius - self.area = np.pi * self.radius ** 2 + self.area = np.pi * self.radius**2 # Center of mass distance to points of interest self.distanceRocketNozzle = distanceRocketNozzle @@ -370,12 +370,14 @@ def evaluateStaticMargin(self): # Calculate total lift coefficient derivative and center of pressure if len(self.aerodynamicSurfaces) > 0: for aerodynamicSurface in self.aerodynamicSurfaces: - self.totalLiftCoeffDer += aerodynamicSurface[1].differentiate( - x=1e-2, dx=1e-3 - ) + self.totalLiftCoeffDer += Function( + lambda alpha: aerodynamicSurface["cl"](alpha, 0) + ).differentiate(x=1e-2, dx=1e-3) self.cpPosition += ( - aerodynamicSurface[1].differentiate(x=1e-2, dx=1e-3) - * aerodynamicSurface[0][2] + Function( + lambda alpha: aerodynamicSurface["cl"](alpha, 0) + ).differentiate(x=1e-2, dx=1e-3) + * aerodynamicSurface["cp"][2] ) self.cpPosition /= self.totalLiftCoeffDer @@ -416,8 +418,13 @@ def addTail(self, topRadius, bottomRadius, length, distanceToCM): Returns ------- - cldata : Function - Object of the Function class. Contains tail's lift data. + cl : Function + Function of the angle of attack (Alpha) and the mach number + (Mach) expressing the tail's lift coefficient. The inputs + are the angle of attack (in radians) and the mach number. + The output is the tail's lift coefficient. In the current + implementation, the tail's lift coefficient does not vary + with mach. self : Rocket Object of the Rocket class. """ @@ -429,22 +436,20 @@ def addTail(self, topRadius, bottomRadius, length, distanceToCM): # Calculate cp position relative to cm if distanceToCM < 0: - cpz = distanceToCM - (length / 3) * (1 + (1 - r) / (1 - r ** 2)) + cpz = distanceToCM - (length / 3) * (1 + (1 - r) / (1 - r**2)) else: - cpz = distanceToCM + (length / 3) * (1 + (1 - r) / (1 - r ** 2)) + cpz = distanceToCM + (length / 3) * (1 + (1 - r) / (1 - r**2)) # Calculate clalpha clalpha = -2 * (1 - r ** (-2)) * (topRadius / rref) ** 2 - cldata = Function( - lambda x: clalpha * x, - "Alpha (rad)", + cl = Function( + lambda alpha, mach: clalpha * alpha, + ["Alpha (rad)", "Mach"], "Cl", - interpolation="linear", - extrapolation="natural", ) # Store values as new aerodynamic surface - tail = [(0, 0, cpz), cldata, "Tail"] + tail = {"cp": (0, 0, cpz), "cl": cl, "name": "Tail"} self.aerodynamicSurfaces.append(tail) # Refresh static margin calculation @@ -476,8 +481,13 @@ def addNose(self, length, kind, distanceToCM): Returns ------- - cldata : Function - Object of the Function class. Contains nose's lift data. + cl : Function + Function of the angle of attack (Alpha) and the mach number + (Mach) expressing the nose cone's lift coefficient. The inputs + are the angle of attack (in radians) and the mach number. + The output is the nose cone's lift coefficient. In the current + implementation, the nose cone's lift coefficient does not vary + with mach self : Rocket Object of the Rocket class. """ @@ -499,16 +509,14 @@ def addNose(self, length, kind, distanceToCM): # Calculate clalpha clalpha = 2 - cldata = Function( - lambda x: clalpha * x, - "Alpha (rad)", + cl = Function( + lambda alpha, mach: clalpha * alpha, + ["Alpha (rad)", "Mach"], "Cl", - interpolation="linear", - extrapolation="natural", ) # Store values - nose = [(0, 0, cpz), cldata, "Nose Cone"] + nose = {"cp": (0, 0, cpz), "cl": cl, "name": "Nose Cone"} self.aerodynamicSurfaces.append(nose) # Refresh static margin calculation @@ -564,8 +572,11 @@ def addFins( Returns ------- - cldata : Function - Object of the Function class. Contains fin's lift data. + cl : Function + Function of the angle of attack (Alpha) and the mach number + (Mach) expressing the fin's lift coefficient. The inputs + are the angle of attack (in radians) and the mach number. + The output is the fin's lift coefficient. self : Rocket Object of the Rocket class. """ @@ -576,16 +587,40 @@ def addFins( Yr = rootChord + tipChord s = span Af = Yr * s / 2 # fin area - Ymac = ( + Yma = ( (s / 3) * (Cr + 2 * Ct) / Yr ) # span wise position of fin's mean aerodynamic chord - Lf = np.sqrt((rootChord / 2 - tipChord / 2) ** 2 + span ** 2) + gamac = np.arctan((Cr - Ct) / (2 * span)) + Lf = np.sqrt((rootChord / 2 - tipChord / 2) ** 2 + span**2) radius = self.radius if radius == 0 else radius d = 2 * radius + Aref = np.pi * radius**2 + AR = 2 * s**2 / Af # Barrowman's convention for fin's aspect ratio cantAngleRad = np.radians(cantAngle) - trapezoidalConstant = ((Yr) / 2) * (radius ** 2) * s - trapezoidalConstant += ((Cr + 2 * Ct) / 3) * radius * (s ** 2) - trapezoidalConstant += ((Cr + 3 * Ct) / 12) * (s ** 3) + trapezoidalConstant = ( + (Cr + 3 * Ct) * s**3 + + 4 * (Cr + 2 * Ct) * radius * s**2 + + 6 * (Cr + Ct) * s * radius**2 + ) / 12 + + # Fin–body interference correction parameters + τ = (s + radius) / radius + λ = Ct / Cr + liftInterferenceFactor = 1 + 1 / τ + rollForcingInterferenceFactor = (1 / np.pi**2) * ( + (np.pi**2 / 4) * ((τ + 1) ** 2 / τ**2) + + ((np.pi * (τ**2 + 1) ** 2) / (τ**2 * (τ - 1) ** 2)) + * np.arcsin((τ**2 - 1) / (τ**2 + 1)) + - (2 * np.pi * (τ + 1)) / (τ * (τ - 1)) + + ((τ**2 + 1) ** 2) + / (τ**2 * (τ - 1) ** 2) + * (np.arcsin((τ**2 - 1) / (τ**2 + 1))) ** 2 + - (4 * (τ + 1)) / (τ * (τ - 1)) * np.arcsin((τ**2 - 1) / (τ**2 + 1)) + + (8 / (τ - 1) ** 2) * np.log((τ**2 + 1) / (2 * τ)) + ) + rollDampingInterferenceFactor = 1 + ( + ((τ - λ) / (τ)) - ((1 - λ) / (τ - 1)) * np.log(τ) + ) / (((τ + 1) * (τ - λ)) / (2) - ((1 - λ) * (τ**3 - 1)) / (3 * (τ - 1))) # Save geometric parameters for later Fin Flutter Analysis and Roll Moment Calculation self.rootChord = Cr @@ -593,6 +628,55 @@ def addFins( self.span = s self.distanceRocketFins = distanceToCM + # Auxiliary functions + + # Defines beta parameter + def beta(mach): + """Defines a parameter that is commonly used in aerodynamic + equations. It is commonly used in the Prandtl factor which + corrects subsonic force coefficients for compressible flow. + + Parameters + ---------- + mach : int, float + Number of mach. + + Returns + ------- + beta : int, float + Value that characterizes flow speed based on the mach number. + """ + + if mach < 0.8: + return np.sqrt(1 - mach**2) + elif mach < 1.1: + return np.sqrt(1 - 0.8**2) + else: + return np.sqrt(mach**2 - 1) + + # Defines number of fins correction + def finNumCorrection(n): + """Calculates a corrector factor for the lift coefficient of multiple fins. + The specifics values are documented at: + Niskanen, S. (2013). “OpenRocket technical documentation”. In: Development + of an Open Source model rocket simulation software. + + Parameters + ---------- + n : int + Number of fins. + + Returns + ------- + Corrector factor : int + Factor that accounts for the number of fins. + """ + correctorFactor = [2.37, 2.74, 2.99, 3.24] + if n >= 5 and n <= 8: + return correctorFactor[n - 5] + else: + return n / 2 + # Calculate cp position relative to cm if distanceToCM < 0: cpz = distanceToCM - ( @@ -607,44 +691,35 @@ def addFins( # Calculate lift parameters for planar fins if not airfoil: - # Calculate clalpha - clalpha = (4 * n * (s / d) ** 2) / (1 + np.sqrt(1 + (2 * Lf / Yr) ** 2)) - clalpha *= 1 + radius / (s + radius) - # # Create a function of lift values by attack angle - cldata = Function( - lambda x: clalpha * x, "Alpha (rad)", "Cl", interpolation="linear" - ) - # Parameters for Roll Moment. Documented at: https://github.com/Projeto-Jupiter/RocketPy/blob/develop/docs/technical/aerodynamics/Roll_Equations.pdf - clfDelta = n * (Ymac + radius) * clalpha / d - cldOmega = ( - n * clalpha * np.cos(cantAngleRad) * trapezoidalConstant / (Af * d) + clalphaSingleFin = Function( + lambda mach: 2 + * np.pi + * AR + * (Af / Aref) + / (2 + np.sqrt(4 + ((beta(mach) * AR) / (np.cos(gamac))) ** 2)), ) - rollParameters = ( - [clfDelta, cldOmega, cantAngleRad] if cantAngleRad != 0 else [0, 0, 0] - ) - - # Store values - fin = [(0, 0, cpz), cldata, rollParameters, "Fins"] - self.aerodynamicSurfaces.append(fin) - # Refresh static margin calculation - self.evaluateStaticMargin() + clalphaMultipleFins = ( + liftInterferenceFactor * finNumCorrection(n) * clalphaSingleFin + ) # Function of mach number - # Return self - return self.aerodynamicSurfaces[-1] + # Calculates clalpha * alpha + cl = Function( + lambda alpha, mach: alpha * clalphaMultipleFins(mach), + ["Alpha (rad)", "Mach"], + "Cl", + ) else: def cnalfa1(cn): """Calculates the normal force coefficient derivative of a 3D airfoil for a given Cnalfa0 - Parameters ---------- cn : int Normal force coefficient derivative of a 2D airfoil. - Returns ------- Cnalfa1 : int @@ -654,7 +729,7 @@ def cnalfa1(cn): # Retrieve parameters for calculations Af = (Cr + Ct) * span / 2 # fin area - AR = 2 * (span ** 2) / Af # Aspect ratio + AR = 2 * (span**2) / Af # Aspect ratio gamac = np.arctan((Cr - Ct) / (2 * span)) # mid chord angle FD = 2 * np.pi * AR / (cn * np.cos(gamac)) @@ -663,7 +738,7 @@ def cnalfa1(cn): * FD * (Af / self.area) * np.cos(gamac) - / (2 + FD * (1 + (4 / FD ** 2)) ** 0.5) + / (2 + FD * (1 + (4 / FD**2)) ** 0.5) ) return Cnalfa1 @@ -672,7 +747,7 @@ def cnalfa1(cn): # Applies number of fins to lift coefficient data data = [[cl[0], (n / 2) * cnalfa1(cl[1])] for cl in read] - cldata = Function( + cl = Function( data, "Alpha (rad)", "Cl", @@ -681,26 +756,39 @@ def cnalfa1(cn): ) # Takes an approximation to an angular coefficient - clalpha = cldata.differentiate(x=0, dx=1e-2) - - # Parameters for Roll Moment. Documented at: https://github.com/Projeto-Jupiter/RocketPy/blob/develop/docs/technical/aerodynamics/Roll_Equations.pdf - clfDelta = n * (Ymac + radius) * clalpha / d - cldOmega = ( - n * clalpha * np.cos(cantAngleRad) * trapezoidalConstant / (Af * d) - ) - rollParameters = ( - [clfDelta, cldOmega, cantAngleRad] if cantAngleRad != 0 else [0, 0, 0] - ) + clalpha = cl.differentiate(x=0, dx=1e-2) + + # Parameters for Roll Moment. + # Documented at: https://github.com/Projeto-Jupiter/RocketPy/blob/develop/docs/technical/aerodynamics/Roll_Equations.pdf + clfDelta = ( + rollForcingInterferenceFactor * n * (Yma + radius) * clalphaSingleFin / d + ) # Function of mach number + cldOmega = ( + 2 + * rollDampingInterferenceFactor + * n + * clalphaSingleFin + * np.cos(cantAngleRad) + * trapezoidalConstant + / (Aref * d**2) + ) + # Function of mach number + rollParameters = [clfDelta, cldOmega, cantAngleRad] - # Store values - fin = [(0, 0, cpz), cldata, rollParameters, "Fins"] - self.aerodynamicSurfaces.append(fin) + # Store values + fin = { + "cp": (0, 0, cpz), + "cl": cl, + "roll parameters": rollParameters, + "name": "Fins", + } + self.aerodynamicSurfaces.append(fin) - # Refresh static margin calculation - self.evaluateStaticMargin() + # Refresh static margin calculation + self.evaluateStaticMargin() - # Return self - return self.aerodynamicSurfaces[-1] + # Return self + return self.aerodynamicSurfaces[-1] def addParachute( self, name, CdS, trigger, samplingRate=100, lag=0, noise=(0, 0, 0) @@ -986,16 +1074,18 @@ def allInfo(self): # Print rocket aerodynamics quantities print("\nAerodynamics Lift Coefficient Derivatives") for aerodynamicSurface in self.aerodynamicSurfaces: - name = aerodynamicSurface[-1] - clalpha = aerodynamicSurface[1].differentiate(x=1e-2, dx=1e-3) + name = aerodynamicSurface["name"] + clalpha = Function( + lambda alpha: aerodynamicSurface["cl"](alpha, 0), + ).differentiate(x=1e-2, dx=1e-3) print( name + " Lift Coefficient Derivative: {:.3f}".format(clalpha) + "/rad" ) print("\nAerodynamics Center of Pressure") for aerodynamicSurface in self.aerodynamicSurfaces: - name = aerodynamicSurface[-1] - cpz = aerodynamicSurface[0][2] + name = aerodynamicSurface["name"] + cpz = aerodynamicSurface["cp"][2] print(name + " Center of Pressure to CM: {:.3f}".format(cpz) + " m") print( "Distance - Center of Pressure to CM: " diff --git a/rocketpy/SolidMotor.py b/rocketpy/SolidMotor.py index 8a780719f..3fe4a5713 100644 --- a/rocketpy/SolidMotor.py +++ b/rocketpy/SolidMotor.py @@ -253,7 +253,7 @@ def __init__( self.grainInitialVolume = ( self.grainInitialHeight * np.pi - * (self.grainOuterRadius ** 2 - self.grainInitialInnerRadius ** 2) + * (self.grainOuterRadius**2 - self.grainInitialInnerRadius**2) ) self.grainInitialMass = self.grainDensity * self.grainInitialVolume self.propellantInitialMass = self.grainNumber * self.grainInitialMass @@ -423,7 +423,7 @@ def evaluateMass(self): @property def throatArea(self): - return np.pi * self.throatRadius ** 2 + return np.pi * self.throatRadius**2 def evaluateGeometry(self): """Calculates grain inner radius and grain height as a @@ -462,9 +462,9 @@ def geometryDot(y, t): grainMassDot = self.massDot(t) / self.grainNumber rI, h = y rIDot = ( - -0.5 * grainMassDot / (density * np.pi * (rO ** 2 - rI ** 2 + rI * h)) + -0.5 * grainMassDot / (density * np.pi * (rO**2 - rI**2 + rI * h)) ) - hDot = 1.0 * grainMassDot / (density * np.pi * (rO ** 2 - rI ** 2 + rI * h)) + hDot = 1.0 * grainMassDot / (density * np.pi * (rO**2 - rI**2 + rI * h)) return [rIDot, hDot] # Solve the system of differential equations @@ -511,8 +511,8 @@ def evaluateBurnArea(self): 2 * np.pi * ( - self.grainOuterRadius ** 2 - - self.grainInnerRadius ** 2 + self.grainOuterRadius**2 + - self.grainInnerRadius**2 + self.grainInnerRadius * self.grainHeight ) * self.grainNumber @@ -583,8 +583,8 @@ def evaluateInertia(self): grainMassDot = self.massDot / self.grainNumber grainNumber = self.grainNumber grainInertiaI = grainMass * ( - (1 / 4) * (self.grainOuterRadius ** 2 + self.grainInnerRadius ** 2) - + (1 / 12) * self.grainHeight ** 2 + (1 / 4) * (self.grainOuterRadius**2 + self.grainInnerRadius**2) + + (1 / 12) * self.grainHeight**2 ) # Calculate each grain's distance d to propellant center of mass @@ -593,7 +593,7 @@ def evaluateInertia(self): d = d * (self.grainInitialHeight + self.grainSeparation) # Calculate inertia for all grains - self.inertiaI = grainNumber * grainInertiaI + grainMass * np.sum(d ** 2) + self.inertiaI = grainNumber * grainInertiaI + grainMass * np.sum(d**2) self.inertiaI.setOutputs("Propellant Inertia I (kg*m2)") # Inertia I Dot @@ -601,8 +601,8 @@ def evaluateInertia(self): grainInertiaIDot = ( grainMassDot * ( - (1 / 4) * (self.grainOuterRadius ** 2 + self.grainInnerRadius ** 2) - + (1 / 12) * self.grainHeight ** 2 + (1 / 4) * (self.grainOuterRadius**2 + self.grainInnerRadius**2) + + (1 / 12) * self.grainHeight**2 ) + grainMass * ((1 / 2) * self.grainInnerRadius - (1 / 3) * self.grainHeight) @@ -611,7 +611,7 @@ def evaluateInertia(self): # Calculate inertia I dot for all grains self.inertiaIDot = grainNumber * grainInertiaIDot + grainMassDot * np.sum( - d ** 2 + d**2 ) self.inertiaIDot.setOutputs("Propellant Inertia I Dot (kg*m2/s)") @@ -619,13 +619,13 @@ def evaluateInertia(self): self.inertiaZ = ( (1 / 2.0) * self.mass - * (self.grainOuterRadius ** 2 + self.grainInnerRadius ** 2) + * (self.grainOuterRadius**2 + self.grainInnerRadius**2) ) self.inertiaZ.setOutputs("Propellant Inertia Z (kg*m2)") # Inertia Z Dot self.inertiaZDot = (1 / 2.0) * self.massDot * ( - self.grainOuterRadius ** 2 + self.grainInnerRadius ** 2 + self.grainOuterRadius**2 + self.grainInnerRadius**2 ) + self.mass * self.grainInnerRadius * self.burnRate self.inertiaZDot.setOutputs("Propellant Inertia Z Dot (kg*m2/s)") diff --git a/rocketpy/utilities.py b/rocketpy/utilities.py index 0051c0622..7bf5a7232 100644 --- a/rocketpy/utilities.py +++ b/rocketpy/utilities.py @@ -31,5 +31,5 @@ def compute_CdS_from_drop_test( """ - CdS = 2 * rocket_mass * gravity / ((terminal_velocity ** 2) * air_density) + CdS = 2 * rocket_mass * gravity / ((terminal_velocity**2) * air_density) return CdS diff --git a/setup.py b/setup.py index 7910006b5..4ecc61447 100644 --- a/setup.py +++ b/setup.py @@ -6,7 +6,13 @@ setuptools.setup( name="rocketpy", version="0.9.9", - install_requires=["numpy>=1.0", "scipy>=1.0", "matplotlib>=3.0", "requests"], + install_requires=[ + "numpy>=1.0", + "scipy>=1.0", + "matplotlib>=3.0", + "requests", + "pytz", + ], maintainer="RocketPy Developers", author="Giovani Hidalgo Ceotto", author_email="ghceotto@gmail.com", diff --git a/tests/test_environment.py b/tests/test_environment.py index 8015d3620..12753df13 100644 --- a/tests/test_environment.py +++ b/tests/test_environment.py @@ -1,4 +1,5 @@ import datetime +import pytz from unittest.mock import patch import pytest @@ -32,7 +33,7 @@ def test_env_set_date(example_env): tomorrow = datetime.date.today() + datetime.timedelta(days=1) example_env.setDate((tomorrow.year, tomorrow.month, tomorrow.day, 12)) assert example_env.date == datetime.datetime( - tomorrow.year, tomorrow.month, tomorrow.day, 12 + tomorrow.year, tomorrow.month, tomorrow.day, 12, tzinfo=pytz.utc ) diff --git a/tests/test_flight.py b/tests/test_flight.py index eea06b013..075da6579 100644 --- a/tests/test_flight.py +++ b/tests/test_flight.py @@ -1,9 +1,56 @@ import datetime from unittest.mock import patch +import numpy as np import pytest +from rocketpy import Environment, Flight, Rocket, SolidMotor, Function +from scipy import optimize -from rocketpy import Environment, SolidMotor, Rocket, Flight +# Helper functions +def setup_rocket_with_given_static_margin(rocket, static_margin): + """Takes any rocket, removes its aerodynamic surfaces and adds a set of + nose, fins and tail specially designed to have a given static margin. + The rocket is modified in place. + + Parameters + ---------- + rocket : Rocket + Rocket to be modified + static_margin : float + Static margin that the given rocket shall have + + Returns + ------- + rocket : Rocket + Rocket with the given static margin. + """ + + def compute_static_margin_error_given_distance(distanceToCM, static_margin, rocket): + rocket.aerodynamicSurfaces = [] + rocket.addNose(length=0.5, kind="vonKarman", distanceToCM=1.0) + rocket.addFins( + 4, + span=0.100, + rootChord=0.100, + tipChord=0.100, + distanceToCM=distanceToCM, + ) + rocket.addTail( + topRadius=0.0635, + bottomRadius=0.0435, + length=0.060, + distanceToCM=-1.194656, + ) + return rocket.staticMargin(0) - static_margin + + sol = optimize.root_scalar( + compute_static_margin_error_given_distance, + bracket=[-2.0, 2.0], + method="brentq", + args=(static_margin, rocket), + ) + + return rocket @patch("matplotlib.pyplot.show") @@ -207,3 +254,181 @@ def mainTrigger(p, y): ) assert test_flight.allInfo() == None + + +@pytest.mark.parametrize("wind_u, wind_v", [(0, 10), (0, -10), (10, 0), (-10, 0)]) +@pytest.mark.parametrize( + "static_margin, max_time", + [(-0.1, 2), (-0.01, 5), (0, 5), (0.01, 20), (0.1, 20), (1.0, 20)], +) +def test_stability_static_margins(wind_u, wind_v, static_margin, max_time): + """Test stability margins for a constant velocity flight, 100 m/s, wind a + lateral wind speed of 10 m/s. Rocket has infinite mass to prevent side motion. + Check if a restoring moment exists depending on static margins.""" + + # Create an environment with ZERO gravity to keep the rocket's speed constant + Env = Environment(gravity=0, railLength=0, latitude=0, longitude=0, elevation=0) + Env.setAtmosphericModel( + type="CustomAtmosphere", + wind_u=wind_u, + wind_v=wind_v, + pressure=101325, + temperature=300, + ) + # Make sure that the freestreamMach will always be 0, so that the rocket + # behaves as the STATIC (freestreamMach=0) margin predicts + Env.speedOfSound = Function(1e16) + + # Create a motor with ZERO thrust and ZERO mass to keep the rocket's speed constant + DummyMotor = SolidMotor( + thrustSource=1e-300, + burnOut=1e-10, + grainNumber=5, + grainSeparation=5 / 1000, + grainDensity=1e-300, + grainOuterRadius=33 / 1000, + grainInitialInnerRadius=15 / 1000, + grainInitialHeight=120 / 1000, + nozzleRadius=33 / 1000, + throatRadius=11 / 1000, + ) + + # Create a rocket with ZERO drag and HUGE mass to keep the rocket's speed constant + DummyRocket = Rocket( + motor=DummyMotor, + radius=127 / 2000, + mass=1e16, + inertiaI=1, + inertiaZ=0.0351, + distanceRocketNozzle=-1.255, + distanceRocketPropellant=-0.85704, + powerOffDrag=0, + powerOnDrag=0, + ) + DummyRocket.setRailButtons([0.2, -0.5]) + setup_rocket_with_given_static_margin(DummyRocket, static_margin) + + # Simulate + init_pos = [0, 0, 100] # Start at 100 m of altitude + init_vel = [0, 0, 100] # Start at 100 m/s + init_att = [1, 0, 0, 0] # Inclination of 90 deg and heading of 0 deg + init_angvel = [0, 0, 0] + initial_solution = [0] + init_pos + init_vel + init_att + init_angvel + TestFlight = Flight( + rocket=DummyRocket, + environment=Env, + initialSolution=initial_solution, + maxTime=max_time, + maxTimeStep=1e-2, + verbose=False, + ) + TestFlight.postProcess(interpolation="linear") + + # Check stability according to static margin + if wind_u == 0: + moments = TestFlight.M1.source[:, 1] + wind_sign = np.sign(wind_v) + else: # wind_v == 0 + moments = TestFlight.M2.source[:, 1] + wind_sign = -np.sign(wind_u) + + assert ( + (static_margin > 0 and np.max(moments) * np.min(moments) < 0) + or (static_margin < 0 and np.all(moments / wind_sign <= 0)) + or (static_margin == 0 and np.all(np.abs(moments) <= 1e-10)) + ) + + +@patch("matplotlib.pyplot.show") +def test_rolling_flight(mock_show): + test_env = Environment( + railLength=5, + latitude=32.990254, + longitude=-106.974998, + elevation=1400, + datum="WGS84", + ) + tomorrow = datetime.date.today() + datetime.timedelta(days=1) + test_env.setDate( + (tomorrow.year, tomorrow.month, tomorrow.day, 12) + ) # Hour given in UTC time + test_env.setAtmosphericModel(type="StandardAtmosphere") + + test_motor = SolidMotor( + thrustSource="data/motors/Cesaroni_M1670.eng", + burnOut=3.9, + grainNumber=5, + grainSeparation=5 / 1000, + grainDensity=1815, + grainOuterRadius=33 / 1000, + grainInitialInnerRadius=15 / 1000, + grainInitialHeight=120 / 1000, + nozzleRadius=33 / 1000, + throatRadius=11 / 1000, + interpolationMethod="linear", + ) + + test_rocket = Rocket( + motor=test_motor, + radius=127 / 2000, + mass=19.197 - 2.956, + inertiaI=6.60, + inertiaZ=0.0351, + distanceRocketNozzle=-1.255, + distanceRocketPropellant=-0.85704, + powerOffDrag="data/calisto/powerOffDragCurve.csv", + powerOnDrag="data/calisto/powerOnDragCurve.csv", + ) + + test_rocket.setRailButtons([0.2, -0.5]) + + NoseCone = test_rocket.addNose( + length=0.55829, kind="vonKarman", distanceToCM=0.71971 + ) + FinSet = test_rocket.addFins( + 4, + span=0.100, + rootChord=0.120, + tipChord=0.040, + distanceToCM=-1.04956, + cantAngle=0.5, + ) + Tail = test_rocket.addTail( + topRadius=0.0635, bottomRadius=0.0435, length=0.060, distanceToCM=-1.194656 + ) + + def drogueTrigger(p, y): + # p = pressure + # y = [x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3] + # activate drogue when vz < 0 m/s. + return True if y[5] < 0 else False + + def mainTrigger(p, y): + # p = pressure + # y = [x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3] + # activate main when vz < 0 m/s and z < 800 m. + return True if y[5] < 0 and y[2] < 800 else False + + Main = test_rocket.addParachute( + "Main", + CdS=10.0, + trigger=mainTrigger, + samplingRate=105, + lag=1.5, + noise=(0, 8.3, 0.5), + ) + + Drogue = test_rocket.addParachute( + "Drogue", + CdS=1.0, + trigger=drogueTrigger, + samplingRate=105, + lag=1.5, + noise=(0, 8.3, 0.5), + ) + + test_flight = Flight( + rocket=test_rocket, environment=test_env, inclination=85, heading=0 + ) + + assert test_flight.allInfo() == None diff --git a/tests/test_solidmotor.py b/tests/test_solidmotor.py index 8f5a9efc3..da8b6411d 100644 --- a/tests/test_solidmotor.py +++ b/tests/test_solidmotor.py @@ -38,7 +38,7 @@ def test_motor(mock_show): def test_initialize_motor_asserts_dynamic_values(solid_motor): grain_vol = grainInitialHeight * ( - np.pi * (grainOuterRadius ** 2 - grainInitialInnerRadius ** 2) + np.pi * (grainOuterRadius**2 - grainInitialInnerRadius**2) ) grain_mass = grain_vol * grainDensity @@ -73,7 +73,7 @@ def test_grain_geometry_progession_asserts_extreme_values(solid_motor): def test_mass_curve_asserts_extreme_values(solid_motor): grain_vol = grainInitialHeight * ( - np.pi * (grainOuterRadius ** 2 - grainInitialInnerRadius ** 2) + np.pi * (grainOuterRadius**2 - grainInitialInnerRadius**2) ) grain_mass = grain_vol * grainDensity @@ -86,8 +86,8 @@ def test_burn_area_asserts_extreme_values(solid_motor): 2 * np.pi * ( - grainOuterRadius ** 2 - - grainInitialInnerRadius ** 2 + grainOuterRadius**2 + - grainInitialInnerRadius**2 + grainInitialInnerRadius * grainInitialHeight ) * grainNumber @@ -108,20 +108,20 @@ def test_burn_area_asserts_extreme_values(solid_motor): def test_evaluate_inertia_I_asserts_extreme_values(solid_motor): grain_vol = grainInitialHeight * ( - np.pi * (grainOuterRadius ** 2 - grainInitialInnerRadius ** 2) + np.pi * (grainOuterRadius**2 - grainInitialInnerRadius**2) ) grain_mass = grain_vol * grainDensity grainInertiaI_initial = grain_mass * ( - (1 / 4) * (grainOuterRadius ** 2 + grainInitialInnerRadius ** 2) - + (1 / 12) * grainInitialHeight ** 2 + (1 / 4) * (grainOuterRadius**2 + grainInitialInnerRadius**2) + + (1 / 12) * grainInitialHeight**2 ) initialValue = (grainNumber - 1) / 2 d = np.linspace(-initialValue, initialValue, grainNumber) d = d * (grainInitialHeight + grainSeparation) - inertiaI_initial = grainNumber * grainInertiaI_initial + grain_mass * np.sum(d ** 2) + inertiaI_initial = grainNumber * grainInertiaI_initial + grain_mass * np.sum(d**2) assert np.allclose( solid_motor.inertiaI.getSource()[0][-1], inertiaI_initial, atol=0.01 @@ -131,12 +131,12 @@ def test_evaluate_inertia_I_asserts_extreme_values(solid_motor): def test_evaluate_inertia_Z_asserts_extreme_values(solid_motor): grain_vol = grainInitialHeight * ( - np.pi * (grainOuterRadius ** 2 - grainInitialInnerRadius ** 2) + np.pi * (grainOuterRadius**2 - grainInitialInnerRadius**2) ) grain_mass = grain_vol * grainDensity grainInertiaZ_initial = ( - grain_mass * (1 / 2.0) * (grainInitialInnerRadius ** 2 + grainOuterRadius ** 2) + grain_mass * (1 / 2.0) * (grainInitialInnerRadius**2 + grainOuterRadius**2) ) assert np.allclose( @@ -173,7 +173,7 @@ def tests_import_eng_asserts_read_values_correctly(solid_motor): def tests_export_eng_asserts_exported_values_correct(solid_motor): - grain_vol = 0.12 * (np.pi * (0.033 ** 2 - 0.015 ** 2)) + grain_vol = 0.12 * (np.pi * (0.033**2 - 0.015**2)) grain_mass = grain_vol * 1815 * 5 solid_motor.exportEng(fileName="tests/solid_motor.eng", motorName="test_motor")