diff --git a/TODO.md b/TODO.md new file mode 100644 index 000000000..1ecf698a5 --- /dev/null +++ b/TODO.md @@ -0,0 +1,41 @@ +# 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 c6c71e646..758a6d616 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 4a59f377b..30001b109 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,5 +2,4 @@ numpy>=1.0 scipy>=1.0 matplotlib>=3.0 netCDF4>=1.4 -requests -pytz +requests \ No newline at end of file diff --git a/rocketpy/Environment.py b/rocketpy/Environment.py index 0ce542171..62eee271f 100644 --- a/rocketpy/Environment.py +++ b/rocketpy/Environment.py @@ -9,7 +9,6 @@ import bisect import warnings import time -import pytz from datetime import datetime, timedelta from inspect import signature, getsourcelines from collections import namedtuple @@ -26,27 +25,12 @@ 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 @@ -96,11 +80,7 @@ class Environment: Environment.elevation : float Launch site elevation. Environment.date : datetime - 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. + Date time of launch. Topographic information: Environment.elevLonArray: array @@ -310,7 +290,6 @@ 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 @@ -347,13 +326,7 @@ 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 : 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 + datum: Returns ------- @@ -370,11 +343,9 @@ def __init__( # Save date if date != None: - self.setDate(date, timeZone) + self.setDate(date) else: self.date = None - self.localDate = None - self.timeZone = None # Initialize constants self.earthRadius = 6.3781 * (10**6) @@ -407,30 +378,21 @@ def __init__( return None - def setDate(self, date, timeZone="UTC"): + def setDate(self, date): """Set date and time of launch and update weather conditions if date dependent atmospheric model is used. Parameters ---------- - 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". + date : Date + Date object specifying launch date and time. Return ------ None """ - # 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) + # Store date + self.date = datetime(*date) # Update atmospheric conditions if atmosphere type is Forecast, # Reanalysis or Ensemble @@ -516,7 +478,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: @@ -525,7 +487,6 @@ 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 @@ -554,7 +515,7 @@ def setTopographicProfile(self, type, file, dictionary="netCDF4", crs=None): if type == "NASADEM_HGT": if dictionary == "netCDF4": - rootgrp = netCDF4.Dataset(file, "r", format="NETCDF4") + rootgrp = Dataset(file, "r", format="NETCDF4") self.elevLonArray = rootgrp.variables["lon"][:].tolist() self.elevLatArray = rootgrp.variables["lat"][:].tolist() self.elevArray = rootgrp.variables["NASADEM_HGT"][:].tolist() @@ -1725,7 +1686,6 @@ 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. @@ -2119,7 +2079,6 @@ 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. @@ -2861,18 +2820,9 @@ def info(self): """ # Print launch site details print("Launch Site Details") - 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") + print("\nLaunch Rail Length: ", self.rL, " m") + if self.date != None: + print("Launch Date: ", self.date, " 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)) @@ -2889,7 +2839,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) @@ -2900,7 +2850,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 @@ -2909,8 +2859,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") @@ -2997,18 +2947,9 @@ def allInfo(self): # Print launch site details print("\n\nLaunch Site Details") - 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") + print("\nLaunch Rail Length: ", self.rL, " m") + if self.date != None: + print("Launch Date: ", self.date, " 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)) @@ -3017,7 +2958,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) @@ -3028,7 +2969,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 @@ -3037,8 +2978,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") @@ -3565,7 +3506,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 3824af5ab..aa8290730 100644 --- a/rocketpy/Flight.py +++ b/rocketpy/Flight.py @@ -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["cp"][2] + compCp = aerodynamicSurface[0][2] # Component absolute velocity in body frame compVxB = vxB + compCp * omega2 compVyB = vyB - compCp * omega1 @@ -1324,7 +1324,7 @@ def uDot(self, t, u, postProcessing=False): compStreamVzBn = compStreamVzB / compStreamSpeed if -1 * compStreamVzBn < 1: compAttackAngle = np.arccos(-compStreamVzBn) - cLift = aerodynamicSurface["cl"](compAttackAngle, freestreamMach) + cLift = abs(aerodynamicSurface[1](compAttackAngle)) # Component lift force magnitude compLift = ( 0.5 * rho * (compStreamSpeed**2) * self.rocket.area * cLift @@ -1340,25 +1340,12 @@ def uDot(self, t, u, postProcessing=False): M1 -= (compCp + a) * compLiftYB M2 += (compCp + a) * compLiftXB # Calculates Roll Moment - 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 + if aerodynamicSurface[-1] == "Fins": + Clfdelta, Cldomega, cantAngleRad = aerodynamicSurface[2] + if cantAngleRad != 0: + Clf = Clfdelta * cantAngleRad + Cld = Cldomega * omega3 / freestreamSpeed + M3 += Clf - Cld # Calculate derivatives # Angular acceleration alpha1 = ( @@ -1510,7 +1497,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, interpolation="spline", extrapolation="natural"): + def postProcess(self): """Post-process all Flight information produced during simulation. Includes the calculation of maximum values, calculation of secondary values such as energy and conversion @@ -1528,43 +1515,43 @@ def postProcess(self, interpolation="spline", extrapolation="natural"): # Transform solution array into Functions sol = np.array(self.solution) self.x = Function( - sol[:, [0, 1]], "Time (s)", "X (m)", interpolation, extrapolation + sol[:, [0, 1]], "Time (s)", "X (m)", "spline", extrapolation="natural" ) self.y = Function( - sol[:, [0, 2]], "Time (s)", "Y (m)", interpolation, extrapolation + sol[:, [0, 2]], "Time (s)", "Y (m)", "spline", extrapolation="natural" ) self.z = Function( - sol[:, [0, 3]], "Time (s)", "Z (m)", interpolation, extrapolation + sol[:, [0, 3]], "Time (s)", "Z (m)", "spline", extrapolation="natural" ) self.vx = Function( - sol[:, [0, 4]], "Time (s)", "Vx (m/s)", interpolation, extrapolation + sol[:, [0, 4]], "Time (s)", "Vx (m/s)", "spline", extrapolation="natural" ) self.vy = Function( - sol[:, [0, 5]], "Time (s)", "Vy (m/s)", interpolation, extrapolation + sol[:, [0, 5]], "Time (s)", "Vy (m/s)", "spline", extrapolation="natural" ) self.vz = Function( - sol[:, [0, 6]], "Time (s)", "Vz (m/s)", interpolation, extrapolation + sol[:, [0, 6]], "Time (s)", "Vz (m/s)", "spline", extrapolation="natural" ) self.e0 = Function( - sol[:, [0, 7]], "Time (s)", "e0", interpolation, extrapolation + sol[:, [0, 7]], "Time (s)", "e0", "spline", extrapolation="natural" ) self.e1 = Function( - sol[:, [0, 8]], "Time (s)", "e1", interpolation, extrapolation + sol[:, [0, 8]], "Time (s)", "e1", "spline", extrapolation="natural" ) self.e2 = Function( - sol[:, [0, 9]], "Time (s)", "e2", interpolation, extrapolation + sol[:, [0, 9]], "Time (s)", "e2", "spline", extrapolation="natural" ) self.e3 = Function( - sol[:, [0, 10]], "Time (s)", "e3", interpolation, extrapolation + sol[:, [0, 10]], "Time (s)", "e3", "spline", extrapolation="natural" ) self.w1 = Function( - sol[:, [0, 11]], "Time (s)", "ω1 (rad/s)", interpolation, extrapolation + sol[:, [0, 11]], "Time (s)", "ω1 (rad/s)", "spline", extrapolation="natural" ) self.w2 = Function( - sol[:, [0, 12]], "Time (s)", "ω2 (rad/s)", interpolation, extrapolation + sol[:, [0, 12]], "Time (s)", "ω2 (rad/s)", "spline", extrapolation="natural" ) self.w3 = Function( - sol[:, [0, 13]], "Time (s)", "ω3 (rad/s)", interpolation, extrapolation + sol[:, [0, 13]], "Time (s)", "ω3 (rad/s)", "spline", extrapolation="natural" ) # Process second type of outputs - accelerations @@ -1596,12 +1583,12 @@ def postProcess(self, interpolation="spline", extrapolation="natural"): 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)", 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) + 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") # Process third type of outputs - temporary values calculated during integration # Initialize force and atmospheric arrays @@ -1628,35 +1615,25 @@ def postProcess(self, interpolation="spline", extrapolation="natural"): # Call derivatives in post processing mode uDot = currentDerivative(step[0], step[1:], postProcessing=True) # Convert forces and atmospheric arrays to functions - self.R1 = Function(self.R1, "Time (s)", "R1 (N)", interpolation) - self.R2 = Function(self.R2, "Time (s)", "R2 (N)", interpolation) - self.R3 = Function(self.R3, "Time (s)", "R3 (N)", interpolation) - self.M1 = Function(self.M1, "Time (s)", "M1 (Nm)", interpolation) - self.M2 = Function(self.M2, "Time (s)", "M2 (Nm)", interpolation) - self.M3 = Function(self.M3, "Time (s)", "M3 (Nm)", interpolation) + self.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.windVelocityX = Function( - self.windVelocityX, - "Time (s)", - "Wind Velocity X (East) (m/s)", - interpolation, + self.windVelocityX, "Time (s)", "Wind Velocity X (East) (m/s)", "spline" ) self.windVelocityY = Function( - 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.windVelocityY, "Time (s)", "Wind Velocity Y (North) (m/s)", "spline" ) + 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)", interpolation + self.dynamicViscosity, "Time (s)", "Dynamic Viscosity (Pa s)", "spline" ) self.speedOfSound = Function( - self.speedOfSound, "Time (s)", "Speed of Sound (m/s)", interpolation + self.speedOfSound, "Time (s)", "Speed of Sound (m/s)", "spline" ) # Process fourth type of output - values calculated from previous outputs diff --git a/rocketpy/Rocket.py b/rocketpy/Rocket.py index 585184201..43f97c876 100644 --- a/rocketpy/Rocket.py +++ b/rocketpy/Rocket.py @@ -370,14 +370,12 @@ def evaluateStaticMargin(self): # Calculate total lift coefficient derivative and center of pressure if len(self.aerodynamicSurfaces) > 0: for aerodynamicSurface in self.aerodynamicSurfaces: - self.totalLiftCoeffDer += Function( - lambda alpha: aerodynamicSurface["cl"](alpha, 0) - ).differentiate(x=1e-2, dx=1e-3) + self.totalLiftCoeffDer += aerodynamicSurface[1].differentiate( + x=1e-2, dx=1e-3 + ) self.cpPosition += ( - Function( - lambda alpha: aerodynamicSurface["cl"](alpha, 0) - ).differentiate(x=1e-2, dx=1e-3) - * aerodynamicSurface["cp"][2] + aerodynamicSurface[1].differentiate(x=1e-2, dx=1e-3) + * aerodynamicSurface[0][2] ) self.cpPosition /= self.totalLiftCoeffDer @@ -418,13 +416,8 @@ def addTail(self, topRadius, bottomRadius, length, distanceToCM): Returns ------- - 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. + cldata : Function + Object of the Function class. Contains tail's lift data. self : Rocket Object of the Rocket class. """ @@ -442,14 +435,16 @@ def addTail(self, topRadius, bottomRadius, length, distanceToCM): # Calculate clalpha clalpha = -2 * (1 - r ** (-2)) * (topRadius / rref) ** 2 - cl = Function( - lambda alpha, mach: clalpha * alpha, - ["Alpha (rad)", "Mach"], + cldata = Function( + lambda x: clalpha * x, + "Alpha (rad)", "Cl", + interpolation="linear", + extrapolation="natural", ) # Store values as new aerodynamic surface - tail = {"cp": (0, 0, cpz), "cl": cl, "name": "Tail"} + tail = [(0, 0, cpz), cldata, "Tail"] self.aerodynamicSurfaces.append(tail) # Refresh static margin calculation @@ -481,13 +476,8 @@ def addNose(self, length, kind, distanceToCM): Returns ------- - 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 + cldata : Function + Object of the Function class. Contains nose's lift data. self : Rocket Object of the Rocket class. """ @@ -509,14 +499,16 @@ def addNose(self, length, kind, distanceToCM): # Calculate clalpha clalpha = 2 - cl = Function( - lambda alpha, mach: clalpha * alpha, - ["Alpha (rad)", "Mach"], + cldata = Function( + lambda x: clalpha * x, + "Alpha (rad)", "Cl", + interpolation="linear", + extrapolation="natural", ) # Store values - nose = {"cp": (0, 0, cpz), "cl": cl, "name": "Nose Cone"} + nose = [(0, 0, cpz), cldata, "Nose Cone"] self.aerodynamicSurfaces.append(nose) # Refresh static margin calculation @@ -572,11 +564,8 @@ def addFins( Returns ------- - 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. + cldata : Function + Object of the Function class. Contains fin's lift data. self : Rocket Object of the Rocket class. """ @@ -587,40 +576,16 @@ def addFins( Yr = rootChord + tipChord s = span Af = Yr * s / 2 # fin area - Yma = ( + Ymac = ( (s / 3) * (Cr + 2 * Ct) / Yr ) # span wise position of fin's mean aerodynamic chord - 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 = ( - (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))) + trapezoidalConstant = ((Yr) / 2) * (radius**2) * s + trapezoidalConstant += ((Cr + 2 * Ct) / 3) * radius * (s**2) + trapezoidalConstant += ((Cr + 3 * Ct) / 12) * (s**3) # Save geometric parameters for later Fin Flutter Analysis and Roll Moment Calculation self.rootChord = Cr @@ -628,55 +593,6 @@ 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 - ( @@ -691,35 +607,44 @@ def finNumCorrection(n): # 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) - clalphaSingleFin = Function( - lambda mach: 2 - * np.pi - * AR - * (Af / Aref) - / (2 + np.sqrt(4 + ((beta(mach) * AR) / (np.cos(gamac))) ** 2)), + # # 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) + ) + rollParameters = ( + [clfDelta, cldOmega, cantAngleRad] if cantAngleRad != 0 else [0, 0, 0] ) - clalphaMultipleFins = ( - liftInterferenceFactor * finNumCorrection(n) * clalphaSingleFin - ) # Function of mach number + # Store values + fin = [(0, 0, cpz), cldata, rollParameters, "Fins"] + self.aerodynamicSurfaces.append(fin) - # Calculates clalpha * alpha - cl = Function( - lambda alpha, mach: alpha * clalphaMultipleFins(mach), - ["Alpha (rad)", "Mach"], - "Cl", - ) + # Refresh static margin calculation + self.evaluateStaticMargin() + + # Return self + return self.aerodynamicSurfaces[-1] 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 @@ -747,7 +672,7 @@ def cnalfa1(cn): # Applies number of fins to lift coefficient data data = [[cl[0], (n / 2) * cnalfa1(cl[1])] for cl in read] - cl = Function( + cldata = Function( data, "Alpha (rad)", "Cl", @@ -756,39 +681,26 @@ def cnalfa1(cn): ) # Takes an approximation to an angular coefficient - 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] + clalpha = cldata.differentiate(x=0, dx=1e-2) - # Store values - fin = { - "cp": (0, 0, cpz), - "cl": cl, - "roll parameters": rollParameters, - "name": "Fins", - } - self.aerodynamicSurfaces.append(fin) + # 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] + ) - # Refresh static margin calculation - self.evaluateStaticMargin() + # Store values + fin = [(0, 0, cpz), cldata, rollParameters, "Fins"] + self.aerodynamicSurfaces.append(fin) - # Return self - return self.aerodynamicSurfaces[-1] + # Refresh static margin calculation + self.evaluateStaticMargin() + + # Return self + return self.aerodynamicSurfaces[-1] def addParachute( self, name, CdS, trigger, samplingRate=100, lag=0, noise=(0, 0, 0) @@ -1074,18 +986,16 @@ def allInfo(self): # Print rocket aerodynamics quantities print("\nAerodynamics Lift Coefficient Derivatives") for aerodynamicSurface in self.aerodynamicSurfaces: - name = aerodynamicSurface["name"] - clalpha = Function( - lambda alpha: aerodynamicSurface["cl"](alpha, 0), - ).differentiate(x=1e-2, dx=1e-3) + name = aerodynamicSurface[-1] + clalpha = aerodynamicSurface[1].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["name"] - cpz = aerodynamicSurface["cp"][2] + name = aerodynamicSurface[-1] + cpz = aerodynamicSurface[0][2] print(name + " Center of Pressure to CM: {:.3f}".format(cpz) + " m") print( "Distance - Center of Pressure to CM: " diff --git a/setup.py b/setup.py index 4ecc61447..7910006b5 100644 --- a/setup.py +++ b/setup.py @@ -6,13 +6,7 @@ setuptools.setup( name="rocketpy", version="0.9.9", - install_requires=[ - "numpy>=1.0", - "scipy>=1.0", - "matplotlib>=3.0", - "requests", - "pytz", - ], + install_requires=["numpy>=1.0", "scipy>=1.0", "matplotlib>=3.0", "requests"], 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 12753df13..8015d3620 100644 --- a/tests/test_environment.py +++ b/tests/test_environment.py @@ -1,5 +1,4 @@ import datetime -import pytz from unittest.mock import patch import pytest @@ -33,7 +32,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, tzinfo=pytz.utc + tomorrow.year, tomorrow.month, tomorrow.day, 12 ) diff --git a/tests/test_flight.py b/tests/test_flight.py index 075da6579..eea06b013 100644 --- a/tests/test_flight.py +++ b/tests/test_flight.py @@ -1,56 +1,9 @@ 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 -# 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 +from rocketpy import Environment, SolidMotor, Rocket, Flight @patch("matplotlib.pyplot.show") @@ -254,181 +207,3 @@ 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