diff --git a/rocketpy/Environment.py b/rocketpy/Environment.py index 675803c40..cd8863714 100644 --- a/rocketpy/Environment.py +++ b/rocketpy/Environment.py @@ -25,12 +25,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 @@ -348,7 +363,7 @@ def __init__( self.date = 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 @@ -487,6 +502,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 +531,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 +1702,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 +1961,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 +2096,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 +2379,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 @@ -3249,7 +3267,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 +3290,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 +3365,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 +3389,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 +3465,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) ) diff --git a/rocketpy/Flight.py b/rocketpy/Flight.py index 2ac4d1a7e..aa8290730 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 @@ -1314,12 +1314,12 @@ 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: @@ -1327,10 +1327,10 @@ def uDot(self, t, u, postProcessing=False): cLift = abs(aerodynamicSurface[1](compAttackAngle)) # 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 @@ -1351,26 +1351,26 @@ def uDot(self, t, u, postProcessing=False): 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 +1379,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 +1450,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 +1470,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 @@ -1640,19 +1640,19 @@ def postProcess(self): # 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 +1661,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 +1684,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 +1778,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 +1794,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 +1810,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 +1921,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 +1943,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 +1951,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 +2305,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 +2996,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..43f97c876 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 @@ -429,9 +429,9 @@ 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 @@ -579,13 +579,13 @@ def addFins( Ymac = ( (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) + Lf = np.sqrt((rootChord / 2 - tipChord / 2) ** 2 + span**2) radius = self.radius if radius == 0 else radius d = 2 * radius 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 = ((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 @@ -654,7 +654,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 +663,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 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