diff --git a/rocketpy/AeroSurfaces.py b/rocketpy/AeroSurfaces.py index d8ca3435e..775175d53 100644 --- a/rocketpy/AeroSurfaces.py +++ b/rocketpy/AeroSurfaces.py @@ -11,21 +11,33 @@ class AeroSurfaces: - """Class used to hold multiple aerodynamic surfaces and their positions.""" + """Class used to hold one or more aerodynamic surfaces.""" def __init__(self): self._aeroSurfaces = [] - def append(self, aeroSurface, position): - self._aeroSurfaces.append((aeroSurface, position)) + def append(self, aeroSurface, position=None): + if position: + self._aeroSurfaces.append((aeroSurface, position)) + else: + self._aeroSurfaces.append(aeroSurface) def remove(self, aeroSurface): - for surface, position in self._aeroSurfaces: - if surface == aeroSurface: - self._aeroSurfaces.remove((aeroSurface, position)) + for surface in self._aeroSurfaces: + if isinstance(surface, tuple): + if surface[0] == aeroSurface: + self._aeroSurfaces.remove((surface[0], surface[1])) + else: + if surface == aeroSurface: + self._aeroSurfaces.remove(surface) def pop(self, index=-1): - return self._aerosurfaces.pop(index) + return self._aeroSurfaces.pop(index) + + def __repr__(self): + if len(self._aeroSurfaces) == 1: + return self._aeroSurfaces[0].__repr__() + return self._aeroSurfaces.__repr__() def __len__(self): return len(self._aeroSurfaces) @@ -33,6 +45,12 @@ def __len__(self): def __getitem__(self, index): return self._aeroSurfaces[index] + def __getattr__(self, name): + if len(self._aeroSurfaces) == 1: + attr = getattr(self._aeroSurfaces[0], name) + return attr + return getattr(self._aeroSurfaces, name) + def __iter__(self): return iter(self._aeroSurfaces) @@ -75,7 +93,12 @@ class NoseCone: """ def __init__( - self, length, kind, baseRadius=None, rocketRadius=None, name="Nose Cone" + self, + length, + kind, + baseRadius=None, + rocketRadius=None, + name="Nose Cone", ): """Initializes the nose cone. It is used to define the nose cone length, kind, center of pressure and lift coefficient curve. @@ -85,7 +108,8 @@ def __init__( length : float Nose cone length. Has units of length and must be given in meters. kind : string - Nose cone kind. Can be "conical", "ogive" or "lvhaack". + Nose cone kind. Can be "conical", "ogive", "elliptical", "tangent", + "von karman", "parabolic" or "lvhaack". baseRadius : float, optional Nose cone base radius. Has units of length and must be given in meters. If not given, the ratio between baseRadius and rocketRadius will be @@ -102,34 +126,120 @@ def __init__( ------- None """ - self.length = length + self.cpy = 0 + self.cpx = 0 + self.position = None # in relation to rocket + + self._rocketRadius = rocketRadius + self._baseRadius = baseRadius + self._length = length self.kind = kind self.name = name - self.baseRadius = baseRadius - self.rocketRadius = rocketRadius - if self.baseRadius is None or self.rocketRadius is None: - self.radiusRatio = 1 - else: - self.radiusRatio = self.baseRadius / self.rocketRadius + self.evaluateGeometricalParameters() + self.evaluateLiftCoefficient() + self.evaluateCenterOfPressure() + # # Store values + # nose = {"cp": (0, 0, self.cpz), "cl": self.cl, "name": name} + + return None + + @property + def rocketRadius(self): + return self._rocketRadius + + @rocketRadius.setter + def rocketRadius(self, value): + self._rocketRadius = value + self.evaluateGeometricalParameters() + self.evaluateLiftCoefficient() + + @property + def baseRadius(self): + return self._baseRadius + + @baseRadius.setter + def baseRadius(self, value): + self._baseRadius = value + self.evaluateGeometricalParameters() + self.evaluateLiftCoefficient() + + @property + def length(self): + return self._length + + @length.setter + def length(self, value): + self._length = value + self.evaluateCenterOfPressure() + + @property + def kind(self): + return self._kind + @kind.setter + def kind(self, value): # Analyze type - if self.kind == "conical": + self._kind = value + if value == "conical": self.k = 2 / 3 - elif self.kind == "ogive": + elif value == "ogive": self.k = 0.466 - elif self.kind == "lvhaack": + elif value == "lvhaack": self.k = 0.563 + elif value == "tangent": + rho = (self.baseRadius**2 + self.length**2) / (2 * self.baseRadius) + volume = np.pi * ( + self.length * rho**2 + - (self.length**3) / 3 + - (rho - self.baseRadius) * rho**2 * np.arcsin(self.length / rho) + ) + area = np.pi * self.baseRadius**2 + self.k = 1 - volume / (area * self.length) + elif value == "elliptical": + self.k = 1 / 3 else: - self.k = 0.5 + self.k = 0.5 # Parabolic and Von Karman + self.evaluateCenterOfPressure() - # Calculate cp position in local coordinates - # Local coordinate origin is found at the tip of the nose cone - self.cpz = self.k * length - self.cpy = 0 - self.cpx = 0 - self.cp = (self.cpx, self.cpy, self.cpz) + def __repr__(self): + rep = f"NoseCone Object -> Name: {self.name}, kind: {self.kind}" + + return rep + + def evaluateGeometricalParameters(self): + """Calculates and saves nose cone's radius ratio. + + Parameters + ---------- + None + + Returns + ------- + None + """ + if self.baseRadius is None or self.rocketRadius is None: + self.radiusRatio = 1 + else: + self.radiusRatio = self.baseRadius / self.rocketRadius + + def evaluateLiftCoefficient(self): + """Calculates and returns nose cone's lift coefficient. + The lift coefficient is saved and returned. This function + also calculates and saves its lift coefficient derivative. + + Parameters + ---------- + None + Returns + ------- + self.cl : Function + Function of the angle of attack (Alpha) and the mach number + (Mach) expressing the lift coefficient of the nose cone. The inputs + are the angle of attack (in radians) and the mach number. + The output is the lift coefficient of the nose cone. + """ # Calculate clalpha self.clalpha = 2 * self.radiusRatio**2 self.cl = Function( @@ -137,12 +247,29 @@ def __init__( ["Alpha (rad)", "Mach"], "Cl", ) - # # Store values - # nose = {"cp": (0, 0, self.cpz), "cl": self.cl, "name": name} + return self.cl - return None + def evaluateCenterOfPressure(self): + """Calculates and returns the center of pressure of the nose cone in local + coordinates. The center of pressure position is saved and stored as a tuple. + Local coordinate origin is found at the tip of the nose cone. + Parameters + ---------- + None - def geometricInfo(self): + Returns + ------- + self.cp : tuple + Tuple containing cpx, cpy, cpz. + """ + + self.cpz = self.k * self.length + self.cpy = 0 + self.cpx = 0 + self.cp = (self.cpx, self.cpy, self.cpz) + return self.cp + + def geometricalInfo(self): """Prints out all the geometric information of the nose cone. Parameters @@ -153,8 +280,10 @@ def geometricInfo(self): ------- None """ - print("Nose Cone Geometric Information of Nose: {}".format(self.name)) + print(f"\nGeometric Information of {self.name}") print("-------------------------------") + if self.position: + print(f"Position: {self.position:.3f} m") print(f"Length: {self.length:.3f} m") print(f"Kind: {self.kind}") print(f"Base Radius: {self.baseRadius:.3f} m") @@ -174,7 +303,7 @@ def aerodynamicInfo(self): ------- None """ - print(f"Nose Cone Aerodynamic Information of Nose: {self.name}") + print(f"\nAerodynamic Information of {self.name}") print("-------------------------------") print(f"Center of Pressure Position in Local Coordinates: {self.cp} m") print(f"Lift Coefficient Slope: {self.clalpha:.3f} 1/rad") @@ -194,7 +323,7 @@ def allInfo(self): ------- None """ - self.geometricInfo() + self.geometricalInfo() self.aerodynamicInfo() return None @@ -328,18 +457,91 @@ def __init__( Aref = np.pi * rocketRadius**2 # Reference area # Store values - self.n = n - self.rocketRadius = rocketRadius - self.airfoil = airfoil - self.cantAngle = cantAngle - self.rootChord = rootChord - self.span = span + self._n = n + self._rocketRadius = rocketRadius + self._airfoil = airfoil + self._cantAngle = cantAngle + self._rootChord = rootChord + self._span = span self.name = name self.d = d self.Aref = Aref # Reference area + self.position = None # in relation to rocket return None + @property + def n(self): + return self._n + + @n.setter + def n(self, value): + self._n = value + self.evaluateGeometricalParameters() + self.evaluateCenterOfPressure() + self.evaluateLiftCoefficient() + self.evaluateRollParameters() + + @property + def rootChord(self): + return self._rootChord + + @rootChord.setter + def rootChord(self, value): + self._rootChord = value + self.evaluateGeometricalParameters() + self.evaluateCenterOfPressure() + self.evaluateLiftCoefficient() + self.evaluateRollParameters() + + @property + def span(self): + return self._span + + @span.setter + def span(self, value): + self._span = value + self.evaluateGeometricalParameters() + self.evaluateCenterOfPressure() + self.evaluateLiftCoefficient() + self.evaluateRollParameters() + + @property + def rocketRadius(self): + return self._rocketRadius + + @rocketRadius.setter + def rocketRadius(self, value): + self._rocketRadius = value + self.evaluateGeometricalParameters() + self.evaluateCenterOfPressure() + self.evaluateLiftCoefficient() + self.evaluateRollParameters() + + @property + def cantAngle(self): + return self._cantAngle + + @cantAngle.setter + def cantAngle(self, value): + self._cantAngle = value + self.evaluateGeometricalParameters() + self.evaluateCenterOfPressure() + self.evaluateLiftCoefficient() + self.evaluateRollParameters() + + @property + def airfoil(self): + return self._airfoil + + @airfoil.setter + def airfoil(self, value): + self._airfoil = value + self.evaluateGeometricalParameters() + self.evaluateCenterOfPressure() + self.evaluateLiftCoefficient() + self.evaluateRollParameters() + @abstractmethod def evaluateCenterOfPressure(self): """Calculates and returns the fin set's center of pressure position in local @@ -357,6 +559,22 @@ def evaluateCenterOfPressure(self): pass + @abstractmethod + def evaluateGeometricalParameters(self): + """Calculates and saves fin set's geometrical parameters such as the + fins' area, aspect ratio and parameters for roll movement. + + Parameters + ---------- + None + + Returns + ------- + None + """ + + pass + def evaluateLiftCoefficient(self): """Calculates and returns the finset's lift coefficient. The lift coefficient is saved and returned. This function @@ -382,13 +600,13 @@ def evaluateLiftCoefficient(self): else: # Defines clalpha2D as the derivative of the # lift coefficient curve for a specific airfoil - airfoilCl = Function( + self.airfoilCl = Function( self.airfoil[0], interpolation="linear", ) # Differentiating at x = 0 to get cl_alpha - clalpha2D_Mach0 = airfoilCl.differentiate(x=1e-3, dx=1e-3) + clalpha2D_Mach0 = self.airfoilCl.differentiate(x=1e-3, dx=1e-3) # Convert to radians if needed if self.airfoil[1] == "degrees": @@ -551,18 +769,19 @@ def geometricalInfo(self): None """ - print("\n\nGeometrical Parameters\n") + print("\nGeometrical Parameters\n") if isinstance(self, TrapezoidalFins): print("Fin Type: Trapezoidal") print("Tip Chord: {:.3f} m".format(self.tipChord)) else: print("Fin Type: Elliptical") - print("Root Chord: {:.3f} m".format(self.rootChord)) print("Span: {:.3f} m".format(self.span)) + if self.position: + print("Position: {:.3f} m".format(self.position)) print("Cant Angle: {:.3f} °".format(self.cantAngle)) - print("Longitudinal Section Area: {:.3f} m".format(self.Af)) - print("Aspect Ratio: {:.3f} m".format(self.AR)) + print("Longitudinal Section Area: {:.3f} m²".format(self.Af)) + print("Aspect Ratio: {:.3f} ".format(self.AR)) print("Gamma_c: {:.3f} m".format(self.gamma_c)) print("Mean Aerodynamic Chord: {:.3f} m".format(self.Yma)) print( @@ -583,7 +802,7 @@ def aerodynamicInfo(self): ------ None """ - print("\n\nAerodynamic Information") + print("\nAerodynamic Information") print("----------------") print("Lift Interference Factor: {:.3f} m".format(self.liftInterferenceFactor)) print( @@ -591,15 +810,21 @@ def aerodynamicInfo(self): self.cpx, self.cpy, self.cpz ) ) + print() print( "Lift Coefficient derivative as a Function of Alpha and Mach for Single Fin" ) + print() self.clalphaSingleFin() + print() print( "Lift Coefficient derivative as a Function of Alpha and Mach for the Fin Set" ) + print() self.clalphaMultipleFins() + print() print("Lift Coefficient as a Function of Alpha and Mach for the Fin Set") + print() self.cl() return None @@ -633,12 +858,36 @@ def rollInfo(self): self.rollForcingInterferenceFactor ) ) - + # lacks a title for the plot self.rollParameters[0]() + # lacks a title for the plot self.rollParameters[1]() return None + def airfoilInfo(self): + """Prints out airfoil related information of the + fin set. + + Parameters + ---------- + None + + Return + ------ + None + """ + if self.airfoil is not None: + print("\n\nAerodynamic Information\n") + print( + "Airfoil's Lift Curve as a Function of Alpha ({}))".format( + self.airfoil[1] + ) + ) + self.airfoilCl.plot1D() + + return None + def allInfo(self): """Prints out all data and graphs available about the Rocket. @@ -663,6 +912,7 @@ def allInfo(self): self.geometricalInfo() self.aerodynamicInfo() self.rollInfo() + self.airfoilInfo() return None @@ -828,64 +1078,57 @@ def __init__( # Sweep length is given pass - Yr = rootChord + tipChord - Af = Yr * span / 2 # Fin area - AR = 2 * span**2 / Af # Fin aspect ratio - gamma_c = np.arctan((sweepLength + 0.5 * tipChord - 0.5 * rootChord) / (span)) - Yma = ( - (span / 3) * (rootChord + 2 * tipChord) / Yr - ) # Span wise coord of mean aero chord + self._tipChord = tipChord + self._sweepLength = sweepLength + self._sweepAngle = sweepAngle - # Fin–body interference correction parameters - tau = (span + rocketRadius) / rocketRadius - liftInterferenceFactor = 1 + 1 / tau - λ = tipChord / rootChord + self.evaluateGeometricalParameters() + self.evaluateCenterOfPressure() + self.evaluateLiftCoefficient() + self.evaluateRollParameters() - # Parameters for Roll Moment. - # Documented at: https://github.com/RocketPy-Team/RocketPy/blob/master/docs/technical/aerodynamics/Roll_Equations.pdf - rollGeometricalConstant = ( - (rootChord + 3 * tipChord) * span**3 - + 4 * (rootChord + 2 * tipChord) * rocketRadius * span**2 - + 6 * (rootChord + tipChord) * span * rocketRadius**2 - ) / 12 - rollDampingInterferenceFactor = 1 + ( - ((tau - λ) / (tau)) - ((1 - λ) / (tau - 1)) * np.log(tau) - ) / ( - ((tau + 1) * (tau - λ)) / (2) - ((1 - λ) * (tau**3 - 1)) / (3 * (tau - 1)) - ) - rollForcingInterferenceFactor = (1 / np.pi**2) * ( - (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) - + ((np.pi * (tau**2 + 1) ** 2) / (tau**2 * (tau - 1) ** 2)) - * np.arcsin((tau**2 - 1) / (tau**2 + 1)) - - (2 * np.pi * (tau + 1)) / (tau * (tau - 1)) - + ((tau**2 + 1) ** 2) - / (tau**2 * (tau - 1) ** 2) - * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2 - - (4 * (tau + 1)) - / (tau * (tau - 1)) - * np.arcsin((tau**2 - 1) / (tau**2 + 1)) - + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau)) - ) + @property + def tipChord(self): + return self._tipChord - self.tipChord = tipChord - self.sweepLength = sweepLength - self.sweepAngle = sweepAngle - self.Yr = Yr - self.Af = Af # Fin area - self.AR = AR # Aspect Ratio - self.gamma_c = gamma_c # Mid chord angle - self.Yma = Yma # Span wise coord of mean aero chord - self.rollGeometricalConstant = rollGeometricalConstant - self.tau = tau - self.liftInterferenceFactor = liftInterferenceFactor - self.λ = λ - self.rollDampingInterferenceFactor = rollDampingInterferenceFactor - self.rollForcingInterferenceFactor = rollForcingInterferenceFactor + @tipChord.setter + def tipChord(self, value): + self._tipChord = value + self.evaluateGeometricalParameters() + self.evaluateCenterOfPressure() + self.evaluateLiftCoefficient() + self.evaluateRollParameters() + @property + def sweepAngle(self): + return self._sweepAngle + + @sweepAngle.setter + def sweepAngle(self, value): + self._sweepAngle = value + self._sweepLength = np.tan(value * np.pi / 180) * self.span + self.evaluateGeometricalParameters() self.evaluateCenterOfPressure() self.evaluateLiftCoefficient() self.evaluateRollParameters() + @property + def sweepLength(self): + return self._sweepLength + + @sweepLength.setter + def sweepLength(self, value): + self._sweepLength = value + self.evaluateGeometricalParameters() + self.evaluateCenterOfPressure() + self.evaluateLiftCoefficient() + self.evaluateRollParameters() + + def __repr__(self): + rep = f"TrapezoidalFins Object -> Name: {self.name}" + + return rep + def evaluateCenterOfPressure(self): """Calculates and returns the center of pressure of the fin set in local coordinates. The center of pressure position is saved and stored as a tuple. @@ -984,7 +1227,7 @@ def draw(self): ) # Plotting - fig3 = plt.figure(figsize=(4, 4)) + fig3 = plt.figure(figsize=(7, 4)) with plt.style.context("bmh"): ax1 = fig3.add_subplot(111) @@ -1015,10 +1258,82 @@ def draw(self): ax1.set_title("Trapezoidal Fin") ax1.legend(bbox_to_anchor=(1.05, 1.0), loc="upper left") + plt.tight_layout() plt.show() return None + def evaluateGeometricalParameters(self): + """Calculates and saves fin set's geometrical parameters such as the + fins' area, aspect ratio and parameters for roll movement. + + Parameters + ---------- + None + + Returns + ------- + None + """ + + Yr = self.rootChord + self.tipChord + Af = Yr * self.span / 2 # Fin area + AR = 2 * self.span**2 / Af # Fin aspect ratio + gamma_c = np.arctan( + (self.sweepLength + 0.5 * self.tipChord - 0.5 * self.rootChord) + / (self.span) + ) + Yma = ( + (self.span / 3) * (self.rootChord + 2 * self.tipChord) / Yr + ) # Span wise coord of mean aero chord + + # Fin–body interference correction parameters + tau = (self.span + self.rocketRadius) / self.rocketRadius + liftInterferenceFactor = 1 + 1 / tau + λ = self.tipChord / self.rootChord + + # Parameters for Roll Moment. + # Documented at: https://github.com/RocketPy-Team/RocketPy/blob/master/docs/technical/aerodynamics/Roll_Equations.pdf + rollGeometricalConstant = ( + (self.rootChord + 3 * self.tipChord) * self.span**3 + + 4 + * (self.rootChord + 2 * self.tipChord) + * self.rocketRadius + * self.span**2 + + 6 * (self.rootChord + self.tipChord) * self.span * self.rocketRadius**2 + ) / 12 + rollDampingInterferenceFactor = 1 + ( + ((tau - λ) / (tau)) - ((1 - λ) / (tau - 1)) * np.log(tau) + ) / ( + ((tau + 1) * (tau - λ)) / (2) - ((1 - λ) * (tau**3 - 1)) / (3 * (tau - 1)) + ) + rollForcingInterferenceFactor = (1 / np.pi**2) * ( + (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) + + ((np.pi * (tau**2 + 1) ** 2) / (tau**2 * (tau - 1) ** 2)) + * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + - (2 * np.pi * (tau + 1)) / (tau * (tau - 1)) + + ((tau**2 + 1) ** 2) + / (tau**2 * (tau - 1) ** 2) + * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2 + - (4 * (tau + 1)) + / (tau * (tau - 1)) + * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau)) + ) + + # Store values + self.Yr = Yr + self.Af = Af # Fin area + self.AR = AR # Aspect Ratio + self.gamma_c = gamma_c # Mid chord angle + self.Yma = Yma # Span wise coord of mean aero chord + self.rollGeometricalConstant = rollGeometricalConstant + self.tau = tau + self.liftInterferenceFactor = liftInterferenceFactor + self.λ = λ + self.rollDampingInterferenceFactor = rollDampingInterferenceFactor + self.rollForcingInterferenceFactor = rollForcingInterferenceFactor + class EllipticalFins(Fins): """Class that defines and holds information for an elliptical fin set. @@ -1163,83 +1478,18 @@ def __init__( name, ) - # Compute auxiliary geometrical parameters - Af = (np.pi * rootChord / 2 * span) / 2 # Fin area - gamma_c = 0 # Zero for elliptical fins - AR = 2 * span**2 / Af # Fin aspect ratio - Yma = ( - span / (3 * np.pi) * np.sqrt(9 * np.pi**2 - 64) - ) # Span wise coord of mean aero chord - rollGeometricalConstant = ( - rootChord - * span - * ( - 3 * np.pi * span**2 - + 32 * rocketRadius * span - + 12 * np.pi * rocketRadius**2 - ) - / 48 - ) - - # Fin–body interference correction parameters - tau = (span + rocketRadius) / rocketRadius - liftInterferenceFactor = 1 + 1 / tau - rollDampingInterferenceFactor = 1 + ( - (rocketRadius**2) - * ( - 2 - * (rocketRadius**2) - * np.sqrt(span**2 - rocketRadius**2) - * np.log( - (2 * span * np.sqrt(span**2 - rocketRadius**2) + 2 * span**2) - / rocketRadius - ) - - 2 - * (rocketRadius**2) - * np.sqrt(span**2 - rocketRadius**2) - * np.log(2 * span) - + 2 * span**3 - - np.pi * rocketRadius * span**2 - - 2 * (rocketRadius**2) * span - + np.pi * rocketRadius**3 - ) - ) / ( - 2 - * (span**2) - * (span / 3 + np.pi * rocketRadius / 4) - * (span**2 - rocketRadius**2) - ) - rollForcingInterferenceFactor = (1 / np.pi**2) * ( - (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) - + ((np.pi * (tau**2 + 1) ** 2) / (tau**2 * (tau - 1) ** 2)) - * np.arcsin((tau**2 - 1) / (tau**2 + 1)) - - (2 * np.pi * (tau + 1)) / (tau * (tau - 1)) - + ((tau**2 + 1) ** 2) - / (tau**2 * (tau - 1) ** 2) - * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2 - - (4 * (tau + 1)) - / (tau * (tau - 1)) - * np.arcsin((tau**2 - 1) / (tau**2 + 1)) - + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau)) - ) - - # Store values - self.Af = Af # Fin area - self.AR = AR # Fin aspect ratio - self.gamma_c = gamma_c # Mid chord angle - self.Yma = Yma # Span wise coord of mean aero chord - self.rollGeometricalConstant = rollGeometricalConstant - self.tau = tau - self.liftInterferenceFactor = liftInterferenceFactor - self.rollDampingInterferenceFactor = rollDampingInterferenceFactor - self.rollForcingInterferenceFactor = rollForcingInterferenceFactor - + self.evaluateGeometricalParameters() self.evaluateCenterOfPressure() self.evaluateLiftCoefficient() self.evaluateRollParameters() return None + def __repr__(self): + rep = f"EllipticalFins Object -> Name: {self.name}" + + return rep + def evaluateCenterOfPressure(self): """Calculates and returns the center of pressure of the fin set in local coordinates. The center of pressure position is saved and stored as a tuple. @@ -1259,7 +1509,7 @@ def evaluateCenterOfPressure(self): self.cpy = 0 self.cpz = cpz self.cp = (self.cpx, self.cpy, self.cpz) - return self + return self.cp def draw(self): """Draw the fin shape along with some important information. @@ -1309,7 +1559,7 @@ def draw(self): cp_point = [self.cpz, self.Yma] # Plotting - fig3 = plt.figure(figsize=(4, 4)) + fig3 = plt.figure(figsize=(7, 4)) with plt.style.context("bmh"): ax1 = fig3.add_subplot(111) ax1.add_patch(el) @@ -1328,10 +1578,98 @@ def draw(self): ax1.set_title("Elliptical Fin") ax1.legend(bbox_to_anchor=(1.05, 1.0), loc="upper left") + plt.tight_layout() plt.show() return None + def evaluateGeometricalParameters(self): + """Calculates and saves fin set's geometrical parameters such as the + fins' area, aspect ratio and parameters for roll movement. + + Parameters + ---------- + None + + Returns + ------- + None + """ + + # Compute auxiliary geometrical parameters + Af = (np.pi * self.rootChord / 2 * self.span) / 2 # Fin area + gamma_c = 0 # Zero for elliptical fins + AR = 2 * self.span**2 / Af # Fin aspect ratio + Yma = ( + self.span / (3 * np.pi) * np.sqrt(9 * np.pi**2 - 64) + ) # Span wise coord of mean aero chord + rollGeometricalConstant = ( + self.rootChord + * self.span + * ( + 3 * np.pi * self.span**2 + + 32 * self.rocketRadius * self.span + + 12 * np.pi * self.rocketRadius**2 + ) + / 48 + ) + + # Fin–body interference correction parameters + tau = (self.span + self.rocketRadius) / self.rocketRadius + liftInterferenceFactor = 1 + 1 / tau + rollDampingInterferenceFactor = 1 + ( + (self.rocketRadius**2) + * ( + 2 + * (self.rocketRadius**2) + * np.sqrt(self.span**2 - self.rocketRadius**2) + * np.log( + ( + 2 * self.span * np.sqrt(self.span**2 - self.rocketRadius**2) + + 2 * self.span**2 + ) + / self.rocketRadius + ) + - 2 + * (self.rocketRadius**2) + * np.sqrt(self.span**2 - self.rocketRadius**2) + * np.log(2 * self.span) + + 2 * self.span**3 + - np.pi * self.rocketRadius * self.span**2 + - 2 * (self.rocketRadius**2) * self.span + + np.pi * self.rocketRadius**3 + ) + ) / ( + 2 + * (self.span**2) + * (self.span / 3 + np.pi * self.rocketRadius / 4) + * (self.span**2 - self.rocketRadius**2) + ) + rollForcingInterferenceFactor = (1 / np.pi**2) * ( + (np.pi**2 / 4) * ((tau + 1) ** 2 / tau**2) + + ((np.pi * (tau**2 + 1) ** 2) / (tau**2 * (tau - 1) ** 2)) + * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + - (2 * np.pi * (tau + 1)) / (tau * (tau - 1)) + + ((tau**2 + 1) ** 2) + / (tau**2 * (tau - 1) ** 2) + * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2 + - (4 * (tau + 1)) + / (tau * (tau - 1)) + * np.arcsin((tau**2 - 1) / (tau**2 + 1)) + + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau)) + ) + + # Store values + self.Af = Af # Fin area + self.AR = AR # Fin aspect ratio + self.gamma_c = gamma_c # Mid chord angle + self.Yma = Yma # Span wise coord of mean aero chord + self.rollGeometricalConstant = rollGeometricalConstant + self.tau = tau + self.liftInterferenceFactor = liftInterferenceFactor + self.rollDampingInterferenceFactor = rollDampingInterferenceFactor + self.rollForcingInterferenceFactor = rollForcingInterferenceFactor + class Tail: """Class that defines a tail. Currently only accepts conical tails. @@ -1404,14 +1742,77 @@ def __init__(self, topRadius, bottomRadius, length, rocketRadius, name="Tail"): """ # Store arguments as attributes - self.topRadius = topRadius - self.bottomRadius = bottomRadius - self.length = length + self._topRadius = topRadius + self._bottomRadius = bottomRadius + self._length = length + self._rocketRadius = rocketRadius self.name = name - self.rocketRadius = rocketRadius + self.position = None # in relation to rocket + + # Calculate geometrical parameters + self.evaluateGeometricalParameters() + self.evaluateLiftCoefficient() + self.evaluateCenterOfPressure() + + return None + + @property + def topRadius(self): + return self._topRadius + + @topRadius.setter + def topRadius(self, value): + self._topRadius = value + self.evaluateGeometricalParameters() + self.evaluateLiftCoefficient() + self.evaluateCenterOfPressure() - # Calculate ratio between top and bottom radius + @property + def bottomRadius(self): + return self._bottomRadius + + @bottomRadius.setter + def bottomRadius(self, value): + self._bottomRadius = value + self.evaluateGeometricalParameters() + self.evaluateLiftCoefficient() + self.evaluateCenterOfPressure() + + @property + def length(self): + return self._length + + @length.setter + def length(self, value): + self._length = value + self.evaluateGeometricalParameters() + self.evaluateCenterOfPressure() + @property + def rocketRadius(self): + return self._rocketRadius + + @rocketRadius.setter + def rocketRadius(self, value): + self._rocketRadius = value + self.evaluateLiftCoefficient() + + def __repr__(self): + rep = f"Tail Object -> Name: {self.name}" + + return rep + + def evaluateGeometricalParameters(self): + """Calculates and saves tail's slant length and surface area. + + Parameters + ---------- + None + + Returns + ------- + None + """ # Calculate tail slant length self.slantLength = np.sqrt( (self.length) ** 2 + (self.topRadius - self.bottomRadius) ** 2 @@ -1421,31 +1822,61 @@ def __init__(self, topRadius, bottomRadius, length, rocketRadius, name="Tail"): np.pi * self.slantLength * (self.topRadius + self.bottomRadius) ) - # Calculate cp position in local coordinates - r = topRadius / bottomRadius - cpz = (length / 3) * (1 + (1 - r) / (1 - r**2)) + def evaluateLiftCoefficient(self): + """Calculates and returns tail's lift coefficient. + The lift coefficient is saved and returned. This function + also calculates and saves its lift coefficient derivative. + Parameters + ---------- + None + + Returns + ------- + self.cl : Function + Function of the angle of attack (Alpha) and the mach number + (Mach) expressing the lift coefficient of the tail. The inputs + are the angle of attack (in radians) and the mach number. + The output is the lift coefficient of the tail. + """ # Calculate clalpha clalpha = 2 * ( - (bottomRadius / rocketRadius) ** 2 - (topRadius / rocketRadius) ** 2 + (self.bottomRadius / self.rocketRadius) ** 2 + - (self.topRadius / self.rocketRadius) ** 2 ) cl = Function( lambda alpha, mach: clalpha * alpha, ["Alpha (rad)", "Mach"], "Cl", ) + self.cl = cl + self.clalpha = clalpha + return self.cl + + def evaluateCenterOfPressure(self): + """Calculates and returns the center of pressure of the tail in local + coordinates. The center of pressure position is saved and stored as a tuple. + Parameters + ---------- + None + + Returns + ------- + self.cp : tuple + Tuple containing cpx, cpy, cpz. + """ + # Calculate cp position in local coordinates + r = self.topRadius / self.bottomRadius + cpz = (self.length / 3) * (1 + (1 - r) / (1 - r**2)) # Store values as class attributes self.cpx = 0 self.cpy = 0 self.cpz = cpz self.cp = (self.cpx, self.cpy, self.cpz) - self.cl = cl - self.clalpha = clalpha - - return None + return self.cp - def geometricInfo(self): + def geometricalInfo(self): """Prints out all the geometric information of the tail. Returns @@ -1453,18 +1884,22 @@ def geometricInfo(self): None """ - print(f"\nTail name: {self.name}") + print(f"\nGeometric Information of {self.name}") + print("-------------------------------") + if self.position: + print(f"Tail Position: {self.position:.3f} m") print(f"Tail Top Radius: {self.topRadius:.3f} m") print(f"Tail Bottom Radius: {self.bottomRadius:.3f} m") print(f"Tail Length: {self.length:.3f} m") print(f"Reference Radius: {2*self.rocketRadius:.3f} m") print(f"Tail Slant Length: {self.slantLength:.3f} m") - print(f"Tail Surface Area: {self.surfaceArea:.6f} m^2") + print(f"Tail Surface Area: {self.surfaceArea:.6f} m²") return None def aerodynamicInfo(self): - print(f"\nTail name: {self.name}") + print(f"\nAerodynamic Information of {self.name}") + print("-------------------------------") print(f"Tail Center of Pressure Position in Local Coordinates: {self.cp} m") print(f"Tail Lift Coefficient Slope: {self.clalpha:.3f} 1/rad") print("Tail Lift Coefficient as a function of Alpha and Mach:") @@ -1479,7 +1914,7 @@ def allInfo(self): ------- None """ - self.geometricInfo() + self.geometricalInfo() self.aerodynamicInfo() return None diff --git a/rocketpy/Flight.py b/rocketpy/Flight.py index 57bf816a6..d4af252b9 100644 --- a/rocketpy/Flight.py +++ b/rocketpy/Flight.py @@ -1784,7 +1784,7 @@ def density(self): @funcify_method("Time (s)", "Dynamic Viscosity (Pa s)", "spline", "constant") def dynamicViscosity(self): """Air dynamic viscosity felt by the rocket as a rocketpy.Function of time.""" - return self.retrieve_temporary_values_arrays[7] + return self.retrieve_temporary_values_arrays[8] @funcify_method("Time (s)", "Speed of Sound (m/s)", "spline", "constant") def speedOfSound(self): diff --git a/rocketpy/Function.py b/rocketpy/Function.py index 7ef4c8c54..955e2bfb7 100644 --- a/rocketpy/Function.py +++ b/rocketpy/Function.py @@ -8,7 +8,7 @@ import matplotlib.pyplot as plt import numpy as np -from scipy import integrate, linalg +from scipy import integrate, linalg, optimize class Function: @@ -206,6 +206,13 @@ def source(x): # Do things if domDim is 1 if self.__domDim__ == 1: source = source[source[:, 0].argsort()] + + self.xArray = source[:, 0] + self.xmin, self.xmax = self.xArray[0], self.xArray[-1] + + self.yArray = source[:, 1] + self.ymin, self.ymax = self.yArray[0], self.yArray[-1] + # Finally set data source as source self.source = source # Set default interpolation for point source if it hasn't @@ -216,6 +223,15 @@ def source(x): self.setInterpolation(self.__interpolation__) # Do things if function is multivariate else: + self.xArray = source[:, 0] + self.xmin, self.xmax = self.xArray[0], self.xArray[-1] + + self.yArray = source[:, 1] + self.ymin, self.ymax = self.yArray[0], self.yArray[-1] + + self.zArray = source[:, 2] + self.zmin, self.zmax = self.zArray[0], self.zArray[-1] + # Finally set data source as source self.source = source if self.__interpolation__ is None: @@ -290,9 +306,9 @@ def setGetValueOpt(self): self : Function """ # Retrieve general info - xData = self.source[:, 0] - yData = self.source[:, 1] - xmin, xmax = xData[0], xData[-1] + xData = self.xArray + yData = self.yArray + xmin, xmax = self.xmin, self.xmax if self.__extrapolation__ == "zero": extrapolation = 0 # Extrapolation is zero elif self.__extrapolation__ == "natural": @@ -473,7 +489,7 @@ def setDiscrete( if self.__domDim__ == 1: Xs = np.linspace(lower, upper, samples) Ys = self.getValue(Xs.tolist()) if oneByOne else self.getValue(Xs) - self.source = np.concatenate(([Xs], [Ys])).transpose() + self.setSource(np.concatenate(([Xs], [Ys])).transpose()) self.setInterpolation(interpolation) self.setExtrapolation(extrapolation) elif self.__domDim__ == 2: @@ -488,7 +504,7 @@ def setDiscrete( mesh = [[Xs[i], Ys[i]] for i in range(len(Xs))] # Evaluate function at all mesh nodes and convert it to matrix Zs = np.array(self.getValue(mesh)) - self.source = np.concatenate(([Xs], [Ys], [Zs])).transpose() + self.setSource(np.concatenate(([Xs], [Ys], [Zs])).transpose()) self.__interpolation__ = "shepard" return self @@ -572,7 +588,7 @@ def setDiscreteBasedOnModel(self, modelFunction, oneByOne=True): if self.__domDim__ == 1: Xs = modelFunction.source[:, 0] Ys = self.getValue(Xs.tolist()) if oneByOne else self.getValue(Xs) - self.source = np.concatenate(([Xs], [Ys])).transpose() + self.setSource(np.concatenate(([Xs], [Ys])).transpose()) elif self.__domDim__ == 2: # Create nodes to evaluate function Xs = modelFunction.source[:, 0] @@ -582,7 +598,7 @@ def setDiscreteBasedOnModel(self, modelFunction, oneByOne=True): mesh = [[Xs[i], Ys[i]] for i in range(len(Xs))] # Evaluate function at all mesh nodes and convert it to matrix Zs = np.array(self.getValue(mesh)) - self.source = np.concatenate(([Xs], [Ys], [Zs])).transpose() + self.setSource(np.concatenate(([Xs], [Ys], [Zs])).transpose()) self.setInterpolation(modelFunction.__interpolation__) self.setExtrapolation(modelFunction.__extrapolation__) @@ -735,9 +751,9 @@ def getValue(self, *args): if isinstance(args[0], (int, float)): args = [list(args)] x = np.array(args[0]) - xData = self.source[:, 0] - yData = self.source[:, 1] - xmin, xmax = xData[0], xData[-1] + xData = self.xArray + yData = self.yArray + xmin, xmax = self.xmin, self.xmax coeffs = self.__polynomialCoefficients__ A = np.zeros((len(args[0]), coeffs.shape[0])) for i in range(coeffs.shape[0]): @@ -752,13 +768,13 @@ def getValue(self, *args): return ans if len(ans) > 1 else ans[0] # Returns value for spline, akima or linear interpolation function type elif self.__interpolation__ in ["spline", "akima", "linear"]: - if isinstance(args[0], (int, float, complex)): + if isinstance(args[0], (int, float, complex, np.integer)): args = [list(args)] x = [arg for arg in args[0]] - xData = self.source[:, 0] - yData = self.source[:, 1] + xData = self.xArray + yData = self.yArray xIntervals = np.searchsorted(xData, x) - xmin, xmax = xData[0], xData[-1] + xmin, xmax = self.xmin, self.xmax if self.__interpolation__ == "spline": coeffs = self.__splineCoefficients__ for i in range(len(x)): @@ -856,9 +872,9 @@ def getValueOpt_deprecated(self, *args): # Interpolated Function # Retrieve general info - xData = self.source[:, 0] - yData = self.source[:, 1] - xmin, xmax = xData[0], xData[-1] + xData = self.xArray + yData = self.yArray + xmin, xmax = self.xmin, self.xmax if self.__extrapolation__ == "zero": extrapolation = 0 # Extrapolation is zero elif self.__extrapolation__ == "natural": @@ -1006,8 +1022,8 @@ def getValueOpt2(self, *args): # Returns value for spline, akima or linear interpolation function type elif self.__interpolation__ in ["spline", "akima", "linear"]: x = args[0] - xData = self.source[:, 0] - yData = self.source[:, 1] + xData = self.xArray + yData = self.yArray # Hunt in intervals near the last interval which was used. xInterval = self.last_interval if xData[xInterval - 1] <= x <= xData[xInterval]: @@ -1016,7 +1032,7 @@ def getValueOpt2(self, *args): xInterval = np.searchsorted(xData, x) self.last_interval = xInterval if xInterval < len(xData) else 0 # Interval found... keep going - xmin, xmax = xData[0], xData[-1] + xmin, xmax = self.xmin, self.xmax if self.__interpolation__ == "spline": coeffs = self.__splineCoefficients__ if x == xmin or x == xmax: @@ -1100,6 +1116,17 @@ def __len__(self): """ return len(self.source) + def __bool__(self): + """Returns true if self exists. This is to avoid getting into __len__ + method in boolean statements. + + Returns + ------- + bool : bool + Always True. + """ + return True + # Define all conversion methods def toFrequencyDomain(self, lower, upper, samplingFrequency, removeDC=True): """Performs the conversion of the Function to the Frequency Domain and returns @@ -1233,6 +1260,7 @@ def plot1D( forceData=False, forcePoints=False, returnObject=False, + equalAxis=False, ): """Plot 1-Dimensional Function, from a lower limit to an upper limit, by sampling the Function several times in the interval. The title of @@ -1277,8 +1305,8 @@ def plot1D( upper = 10 if upper is None else upper else: # Determine boundaries - xData = self.source[:, 0] - xmin, xmax = xData[0], xData[-1] + xData = self.xArray + xmin, xmax = self.xmin, self.xmax lower = xmin if lower is None else lower upper = xmax if upper is None else upper # Plot data points if forceData = True @@ -1295,6 +1323,8 @@ def plot1D( # Plots function if forcePoints: plt.scatter(x, y, marker="o") + if equalAxis: + plt.axis("equal") plt.plot(x, y) # Turn on grid and set title and axis plt.grid(True) @@ -1362,8 +1392,8 @@ def plot2D( upper = 2 * [upper] if isinstance(upper, (int, float)) else upper else: # Determine boundaries - xData = self.source[:, 0] - yData = self.source[:, 1] + xData = self.xArray + yData = self.yArray xMin, xMax = xData.min(), xData.max() yMin, yMax = yData.min(), yData.max() lower = [xMin, yMin] if lower is None else lower @@ -1556,8 +1586,8 @@ def __interpolatePolynomial__(self): # Find the degree of the polynomial interpolation degree = self.source.shape[0] - 1 # Get x and y values for all supplied points. - x = self.source[:, 0] - y = self.source[:, 1] + x = self.xArray + y = self.yArray # Check if interpolation requires large numbers if np.amax(x) ** degree > 1e308: print( @@ -1576,8 +1606,8 @@ def __interpolatePolynomial__(self): def __interpolateSpline__(self): """Calculate natural spline coefficients that fit the data exactly.""" # Get x and y values for all supplied points - x = self.source[:, 0] - y = self.source[:, 1] + x = self.xArray + y = self.yArray mdim = len(x) h = [x[i + 1] - x[i] for i in range(0, mdim - 1)] # Initialize the matrix @@ -1606,8 +1636,8 @@ def __interpolateSpline__(self): def __interpolateAkima__(self): """Calculate akima spline coefficients that fit the data exactly""" # Get x and y values for all supplied points - x = self.source[:, 0] - y = self.source[:, 1] + x = self.xArray + y = self.yArray # Estimate derivatives at each point d = [0] * len(x) d[0] = (y[1] - y[0]) / (x[1] - x[0]) @@ -1677,7 +1707,7 @@ def __truediv__(self, other): and isinstance(self.source, np.ndarray) and self.__interpolation__ == other.__interpolation__ and self.__inputs__ == other.__inputs__ - and np.array_equal(self.source[:, 0], other.source[:, 0]) + and np.array_equal(self.xArray, other.xArray) ): # Operate on grid values with np.errstate(divide="ignore"): @@ -1700,8 +1730,8 @@ def __truediv__(self, other): # Check if Function object source is array or callable if isinstance(self.source, np.ndarray): # Operate on grid values - Ys = self.source[:, 1] / other - Xs = self.source[:, 0] + Ys = self.yArray / other + Xs = self.xArray source = np.concatenate(([Xs], [Ys])).transpose() # Retrieve inputs, outputs and interpolation inputs = self.__inputs__[:] @@ -1735,8 +1765,8 @@ def __rtruediv__(self, other): if isinstance(other, (float, int, complex)): if isinstance(self.source, np.ndarray): # Operate on grid values - Ys = other / self.source[:, 1] - Xs = self.source[:, 0] + Ys = other / self.yArray + Xs = self.xArray source = np.concatenate(([Xs], [Ys])).transpose() # Retrieve inputs, outputs and interpolation inputs = self.__inputs__[:] @@ -1781,11 +1811,12 @@ def __pow__(self, other): and isinstance(self.source, np.ndarray) and self.__interpolation__ == other.__interpolation__ and self.__inputs__ == other.__inputs__ - and np.array_equal(self.source[:, 0], other.source[:, 0]) + and np.any(self.xArray - other.xArray) == False + and np.array_equal(self.xArray, other.xArray) ): # Operate on grid values - Ys = self.source[:, 1] ** other.source[:, 1] - Xs = self.source[:, 0] + Ys = self.yArray**other.yArray + Xs = self.xArray source = np.concatenate(([Xs], [Ys])).transpose() # Retrieve inputs, outputs and interpolation inputs = self.__inputs__[:] @@ -1802,8 +1833,8 @@ def __pow__(self, other): # Check if Function object source is array or callable if isinstance(self.source, np.ndarray): # Operate on grid values - Ys = self.source[:, 1] ** other - Xs = self.source[:, 0] + Ys = self.yArray**other + Xs = self.xArray source = np.concatenate(([Xs], [Ys])).transpose() # Retrieve inputs, outputs and interpolation inputs = self.__inputs__[:] @@ -1837,8 +1868,8 @@ def __rpow__(self, other): if isinstance(other, (float, int, complex)): if isinstance(self.source, np.ndarray): # Operate on grid values - Ys = other ** self.source[:, 1] - Xs = self.source[:, 0] + Ys = other**self.yArray + Xs = self.xArray source = np.concatenate(([Xs], [Ys])).transpose() # Retrieve inputs, outputs and interpolation inputs = self.__inputs__[:] @@ -1883,11 +1914,11 @@ def __mul__(self, other): and isinstance(self.source, np.ndarray) and self.__interpolation__ == other.__interpolation__ and self.__inputs__ == other.__inputs__ - and np.array_equal(self.source[:, 0], other.source[:, 0]) + and np.array_equal(self.xArray, other.xArray) ): # Operate on grid values - Ys = self.source[:, 1] * other.source[:, 1] - Xs = self.source[:, 0] + Ys = self.yArray * other.yArray + Xs = self.xArray source = np.concatenate(([Xs], [Ys])).transpose() # Retrieve inputs, outputs and interpolation inputs = self.__inputs__[:] @@ -1904,8 +1935,8 @@ def __mul__(self, other): # Check if Function object source is array or callable if isinstance(self.source, np.ndarray): # Operate on grid values - Ys = self.source[:, 1] * other - Xs = self.source[:, 0] + Ys = self.yArray * other + Xs = self.xArray source = np.concatenate(([Xs], [Ys])).transpose() # Retrieve inputs, outputs and interpolation inputs = self.__inputs__[:] @@ -1939,8 +1970,8 @@ def __rmul__(self, other): if isinstance(other, (float, int, complex)): if isinstance(self.source, np.ndarray): # Operate on grid values - Ys = other * self.source[:, 1] - Xs = self.source[:, 0] + Ys = other * self.yArray + Xs = self.xArray source = np.concatenate(([Xs], [Ys])).transpose() # Retrieve inputs, outputs and interpolation inputs = self.__inputs__[:] @@ -1985,11 +2016,11 @@ def __add__(self, other): and isinstance(self.source, np.ndarray) and self.__interpolation__ == other.__interpolation__ and self.__inputs__ == other.__inputs__ - and np.array_equal(self.source[:, 0], other.source[:, 0]) + and np.array_equal(self.xArray, other.xArray) ): # Operate on grid values - Ys = self.source[:, 1] + other.source[:, 1] - Xs = self.source[:, 0] + Ys = self.yArray + other.yArray + Xs = self.xArray source = np.concatenate(([Xs], [Ys])).transpose() # Retrieve inputs, outputs and interpolation inputs = self.__inputs__[:] @@ -2006,8 +2037,8 @@ def __add__(self, other): # Check if Function object source is array or callable if isinstance(self.source, np.ndarray): # Operate on grid values - Ys = self.source[:, 1] + other - Xs = self.source[:, 0] + Ys = self.yArray + other + Xs = self.xArray source = np.concatenate(([Xs], [Ys])).transpose() # Retrieve inputs, outputs and interpolation inputs = self.__inputs__[:] @@ -2041,8 +2072,8 @@ def __radd__(self, other): if isinstance(other, (float, int, complex)): if isinstance(self.source, np.ndarray): # Operate on grid values - Ys = other + self.source[:, 1] - Xs = self.source[:, 0] + Ys = other + self.yArray + Xs = self.xArray source = np.concatenate(([Xs], [Ys])).transpose() # Retrieve inputs, outputs and interpolation inputs = self.__inputs__[:] @@ -2087,11 +2118,11 @@ def __sub__(self, other): and isinstance(self.source, np.ndarray) and self.__interpolation__ == other.__interpolation__ and self.__inputs__ == other.__inputs__ - and np.array_equal(self.source[:, 0], other.source[:, 0]) + and np.array_equal(self.xArray, other.xArray) ): # Operate on grid values - Ys = self.source[:, 1] - other.source[:, 1] - Xs = self.source[:, 0] + Ys = self.yArray - other.yArray + Xs = self.xArray source = np.concatenate(([Xs], [Ys])).transpose() # Retrieve inputs, outputs and interpolation inputs = self.__inputs__[:] @@ -2108,8 +2139,8 @@ def __sub__(self, other): # Check if Function object source is array or callable if isinstance(self.source, np.ndarray): # Operate on grid values - Ys = self.source[:, 1] - other - Xs = self.source[:, 0] + Ys = self.yArray - other + Xs = self.xArray source = np.concatenate(([Xs], [Ys])).transpose() # Retrieve inputs, outputs and interpolation inputs = self.__inputs__[:] @@ -2143,8 +2174,8 @@ def __rsub__(self, other): if isinstance(other, (float, int, complex)): if isinstance(self.source, np.ndarray): # Operate on grid values - Ys = other - self.source[:, 1] - Xs = self.source[:, 0] + Ys = other - self.yArray + Xs = self.xArray source = np.concatenate(([Xs], [Ys])).transpose() # Retrieve inputs, outputs and interpolation inputs = self.__inputs__[:] @@ -2187,17 +2218,24 @@ def integral(self, a, b, numerical=False): a, b = b, a # Different implementations depending on interpolation if self.__interpolation__ == "spline" and numerical is False: - xData = self.source[:, 0] - yData = self.source[:, 1] + xData = self.xArray + yData = self.yArray coeffs = self.__splineCoefficients__ ans = 0 # Check to see if interval starts before point data if a < xData[0]: if self.__extrapolation__ == "constant": - ans += yData[0] * (xData[0] - a) + ans += yData[0] * (min(xData[0], b) - a) elif self.__extrapolation__ == "natural": c = coeffs[:, 0] - subB = a - xData[0] # subA = 0 + subB = a - xData[0] + subA = min(b, xData[0]) - xData[0] + ans += ( + (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) @@ -2207,29 +2245,43 @@ def integral(self, a, b, numerical=False): else: # self.__extrapolation__ = 'zero' pass + # Integrate in subintervals between Xs of given data up to b - i = 0 + i = max(np.searchsorted(xData, a, side="left") - 1, 0) while i < len(xData) - 1 and xData[i] < b: - if b < xData[i + 1]: - subB = b - xData[i] # subA = 0 + if xData[i] <= a <= xData[i + 1] and xData[i] <= b <= xData[i + 1]: + subA = a - xData[i] + subB = b - xData[i] + elif xData[i] <= a <= xData[i + 1]: + subA = a - xData[i] + subB = xData[i + 1] - xData[i] + elif b <= xData[i + 1]: + subA = 0 + subB = b - xData[i] else: - subB = xData[i + 1] - xData[i] # subA = 0 + subA = 0 + subB = xData[i + 1] - xData[i] 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[0] * subB ) + ans -= ( + (c[3] * subA**4) / 4 + + (c[2] * subA**3 / 3) + + (c[1] * subA**2 / 2) + + c[0] * subA + ) i += 1 # Check to see if interval ends after point data if b > xData[-1]: if self.__extrapolation__ == "constant": - ans += yData[-1] * (b - xData[-1]) + ans += yData[-1] * (b - max(xData[-1], a)) elif self.__extrapolation__ == "natural": c = coeffs[:, -1] - subA = xData[-1] - xData[-2] + subA = max(xData[-1], a) - xData[-2] subB = b - xData[-2] ans -= ( (c[3] * subA**4) / 4 @@ -2248,8 +2300,8 @@ def integral(self, a, b, numerical=False): pass elif self.__interpolation__ == "linear" and numerical is False: # Integrate from a to b using np.trapz - xData = self.source[:, 0] - yData = self.source[:, 1] + xData = self.xArray + yData = self.yArray # Get data in interval xIntegrationData = xData[(xData >= a) & (xData <= b)] yIntegrationData = yData[(xData >= a) & (xData <= b)] @@ -2270,15 +2322,426 @@ def integral(self, a, b, numerical=False): ans = np.trapz(yIntegrationData, xIntegrationData) else: # Integrate numerically - ans, _ = integrate.quad(self, a, b, epsabs=0.1, limit=10000) + ans, _ = integrate.quad(self, a, b, epsabs=0.001, limit=10000) return integrationSign * ans - # Not implemented - def differentiate(self, x, dx=1e-6): - return (self.getValue(x + dx) - self.getValue(x - dx)) / (2 * dx) - # h = (10)**-300 - # z = x + h*1j - # return self(z).imag/h + def differentiate(self, x, dx=1e-6, order=1): + """Differentiate a Function object at a given point. + + Parameters + ---------- + x : float + Point at which to differentiate. + dx : float + Step size to use for numerical differentiation. + order : int + Order of differentiation. + + Returns + ------- + ans : float + Evaluated derivative. + """ + if order == 1: + return (self.getValue(x + dx) - self.getValue(x - dx)) / (2 * dx) + elif order == 2: + return ( + self.getValue(x + dx) - 2 * self.getValue(x) + self.getValue(x - dx) + ) / dx**2 + + def derivativeFunction(self): + """Returns a Function object which gives the derivative of the Function object. + + Returns + ------- + result : Function + A Function object which gives the derivative of self. + """ + # Check if Function object source is array + if isinstance(self.source, np.ndarray): + # Operate on grid values + Ys = np.diff(self.yArray) / np.diff(self.xArray) + Xs = self.source[:-1, 0] + np.diff(self.xArray) / 2 + source = np.column_stack((Xs, Ys)) + # Retrieve inputs, outputs and interpolation + inputs = self.__inputs__[:] + outputs = f"d({self.__outputs__[0]})/d({inputs[0]})" + else: + source = lambda x: self.differentiate(x) + inputs = self.__inputs__[:] + outputs = f"d({self.__outputs__[0]})/d({inputs[0]})" + + # Create new Function object + return Function(source, inputs, outputs, self.__interpolation__) + + def integralFunction(self, lower=None, upper=None, datapoints=100): + """Returns a Function object representing the integral of the Function object. + + Parameters + ---------- + lower : scalar, optional + The lower limit of the interval in which the function is to be + evaluated at. If the Function is given by a dataset, the default + value is the start of the dataset. + upper : scalar, optional + The upper limit of the interval in which the function is to be + evaluated at. If the Function is given by a dataset, the default + value is the end of the dataset. + datapoints : int, optional + The number of points in which the integral will be evaluated for + plotting it, which draws lines between each evaluated point. + The default value is 100. + + Returns + ------- + result : Function + The integral of the Function object. + """ + if isinstance(self.source, np.ndarray): + lower = self.source[0, 0] if lower is None else lower + upper = self.source[-1, 0] if upper is None else upper + xData = np.linspace(lower, upper, datapoints) + yData = np.zeros(datapoints) + for i in range(datapoints): + yData[i] = self.integral(lower, xData[i]) + return Function( + np.column_stack((xData, yData)), + inputs=self.__inputs__, + outputs=[o + " Integral" for o in self.__outputs__], + ) + else: + lower = 0 if lower is None else lower + return Function( + lambda x: self.integral(lower, x), + inputs=self.__inputs__, + outputs=[o + " Integral" for o in self.__outputs__], + ) + + def isBijective(self): + """Checks whether the Function is bijective. Only applicable to + Functions whose source is a list of points, raises an error otherwise. + + Returns + ------- + result : bool + True if the Function is bijective, False otherwise. + """ + if isinstance(self.source, np.ndarray): + xDataDistinct = set(self.xArray) + yDataDistinct = set(self.yArray) + distinctMap = set(zip(xDataDistinct, yDataDistinct)) + return len(distinctMap) == len(xDataDistinct) == len(yDataDistinct) + else: + raise TypeError( + "Only Functions whose source is a list of points can be " + "checked for bijectivity." + ) + + def isStrictlyBijective(self): + """Checks whether the Function is "strictly" bijective. + Only applicable to Functions whose source is a list of points, + raises an error otherwise. + + Notes + ----- + By "strictly" bijective, this implementation considers the + list-of-points-defined Function bijective between each consecutive pair + of points. Therefore, the Function may be flagged as not bijective even + if the mapping between the set of points which define the Function is + bijective. + + Returns + ------- + result : bool + True if the Function is "strictly" bijective, False otherwise. + + Examples + -------- + >>> f = Function([[0, 0], [1, 1], [2, 4]]) + >>> f.isBijective() + True + >>> f.isStrictlyBijective() + True + + >>> f = Function([[-1, 1], [0, 0], [1, 1], [2, 4]]) + >>> f.isBijective() + False + >>> f.isStrictlyBijective() + False + + A Function which is not "strictly" bijective, but is bijective, can be + constructed as x^2 defined at -1, 0 and 2. + + >>> f = Function([[-1, 1], [0, 0], [2, 4]]) + >>> f.isBijective() + True + >>> f.isStrictlyBijective() + False + """ + if isinstance(self.source, np.ndarray): + # Assuming domain is sorted, range must also be + yData = self.yArray + # Both ascending and descending order means Function is bijective + yDataDiff = np.diff(yData) + return np.all(yDataDiff >= 0) or np.all(yDataDiff <= 0) + else: + raise TypeError( + "Only Functions whose source is a list of points can be " + "checked for bijectivity." + ) + + def inverseFunction(self, approxFunc=None, tol=1e-4): + """ + Returns the inverse of the Function. The inverse function of F is a + function that undoes the operation of F. The inverse of F exists if + and only if F is bijective. Makes the domain the range and the range + the domain. + + If the Function is given by a list of points, its bijectivity is + checked and an error is raised if it is not bijective. + If the Function is given by a function, its bijection is not + checked and may lead to innacuracies outside of its bijective region. + + Parameters + ---------- + approxFunc : callable, optional + A function that approximates the inverse of the Function. This + function is used to find the starting guesses for the inverse + root finding algorithm. This is better used when the inverse + in complex but has a simple approximation or when the root + finding algorithm performs poorly due to default start point. + The default is None in which case the starting point is zero. + + tol : float, optional + The tolerance for the inverse root finding algorithm. The default + is 1e-4. + + Returns + ------- + result : Function + A Function whose domain and range have been inverted. + """ + if isinstance(self.source, np.ndarray): + if self.isStrictlyBijective(): + # Swap the columns + source = np.flip(self.source, axis=1) + else: + raise ValueError( + "Function is not bijective, so it does not have an inverse." + ) + else: + if approxFunc is not None: + source = lambda x: self.findInput(x, start=approxFunc(x), tol=tol) + else: + source = lambda x: self.findInput(x, start=0, tol=tol) + return Function( + source, + inputs=self.__outputs__, + outputs=self.__inputs__, + interpolation=self.__interpolation__, + ) + + def findInput(self, val, start, tol=1e-4): + """ + Finds the optimal input for a given output. + + Parameters + ---------- + val : int, float + The value of the output. + start : int, float + Initial guess of the output. + tol : int, float + Tolerance for termination. + + Returns + ------- + result : ndarray + The value of the input which gives the output closest to val. + """ + return optimize.root( + lambda x: self.getValue(x) - val, + start, + tol=tol, + ).x[0] + + def average(self, lower, upper): + """ + Returns the average of the function. + + Parameters + ---------- + lower : float + Lower point of the region that the average will be calculated at. + upper : float + Upper point of the region that the average will be calculated at. + + Returns + ------- + result : float + The average of the function. + """ + return self.integral(lower, upper) / (upper - lower) + + def averageFunction(self, lower=None): + """ + Returns a Function object representing the average of the Function object. + + Parameters + ---------- + lower : float + Lower limit of the new domain. Only required if the Function's source + is a callable instead of a list of points. + + Returns + ------- + result : Function + The average of the Function object. + """ + if isinstance(self.source, np.ndarray): + if lower is None: + lower = self.source[0, 0] + upper = self.source[-1, 0] + xData = np.linspace(lower, upper, 100) + yData = np.zeros(100) + yData[0] = self.source[:, 1][0] + for i in range(1, 100): + yData[i] = self.average(lower, xData[i]) + return Function( + np.concatenate(([xData], [yData])).transpose(), + inputs=self.__inputs__, + outputs=[o + " Average" for o in self.__outputs__], + ) + else: + if lower is None: + lower = 0 + return Function( + lambda x: self.average(lower, x), + inputs=self.__inputs__, + outputs=[o + " Average" for o in self.__outputs__], + ) + + def compose(self, func, extrapolate=False): + """ + Returns a Function object which is the result of inputing a function + into a function (i.e. f(g(x))). The domain will become the domain of + the input function and the range will become the range of the original + function. + + Parameters + ---------- + func : Function + The function to be inputed into the function. + + extrapolate : bool, optional + Whether or not to extrapolate the function if the input function's + range is outside of the original function's domain. The default is + False. + + Returns + ------- + result : Function + The result of inputing the function into the function. + """ + # Check if the input is a function + if not isinstance(func, Function): + raise TypeError("Input must be a Function object.") + + if isinstance(self.source, np.ndarray) and isinstance(func.source, np.ndarray): + # Perform bounds check for composition + if not extrapolate: + if func.ymin < self.xmin and func.ymax > self.xmax: + raise ValueError( + f"Input Function image {func.ymin, func.ymax} must be within " + f"the domain of the Function {self.xmin, self.xmax}." + ) + + return Function( + np.concatenate(([func.xArray], [self(func.yArray)])).T, + inputs=func.__inputs__, + outputs=self.__outputs__, + interpolation=self.__interpolation__, + extrapolation=self.__extrapolation__, + ) + else: + return Function( + lambda x: self(func(x)), + inputs=func.__inputs__, + outputs=self.__outputs__, + interpolation=self.__interpolation__, + extrapolation=self.__extrapolation__, + ) + + +class PiecewiseFunction(Function): + def __new__( + cls, + source, + inputs=["Scalar"], + outputs=["Scalar"], + interpolation="akima", + extrapolation=None, + datapoints=50, + ): + """ + Creates a piecewise function from a dictionary of functions. The keys of the dictionary + must be tuples that represent the domain of the function. The domains must be disjoint. + The piecewise function will be evaluated at datapoints points to create Function object. + + Parameters + ---------- + source: dictionary + A dictionary of Function objects, where the keys are the domains. + inputs : list + A list of strings that represent the inputs of the function. + outputs: list + A list of strings that represent the outputs of the function. + interpolation: str + The type of interpolation to use. The default value is 'akima'. + extrapolation: str + The type of extrapolation to use. The default value is None. + datapoints: int + The number of points in which the piecewise function will be + evaluated to create a base function. The default value is 100. + """ + # Check if source is a dictionary + if not isinstance(source, dict): + raise TypeError("source must be a dictionary") + # Check if all keys are tuples + for key in source.keys(): + if not isinstance(key, tuple): + raise TypeError("keys of source must be tuples") + # Check if all domains are disjoint + for key1 in source.keys(): + for key2 in source.keys(): + if key1 != key2: + if key1[0] < key2[1] and key1[1] > key2[0]: + raise ValueError("domains must be disjoint") + + # Crate Function + def calcOutput(func, inputs): + o = np.zeros(len(inputs)) + for j in range(len(inputs)): + o[j] = func.getValue(inputs[j]) + return o + + inputData = [] + outputData = [] + for key in sorted(source.keys()): + i = np.linspace(key[0], key[1], datapoints) + i = i[~np.in1d(i, inputData)] + inputData = np.concatenate((inputData, i)) + + f = Function(source[key]) + outputData = np.concatenate((outputData, calcOutput(f, i))) + + return Function( + np.concatenate(([inputData], [outputData])).T, + inputs=inputs, + outputs=outputs, + interpolation=interpolation, + extrapolation=extrapolation, + ) def funcify_method(*args, **kwargs): diff --git a/rocketpy/Rocket.py b/rocketpy/Rocket.py index 6e8386bb9..381e6321b 100644 --- a/rocketpy/Rocket.py +++ b/rocketpy/Rocket.py @@ -214,6 +214,9 @@ def __init__( # Aerodynamic data initialization self.aerodynamicSurfaces = AeroSurfaces() + self.nosecone = AeroSurfaces() + self.fins = AeroSurfaces() + self.tail = AeroSurfaces() self.cpPosition = 0 self.staticMargin = Function( lambda x: 0, inputs="Time (s)", outputs="Static Margin (c)" @@ -481,9 +484,12 @@ def addTail( # Create new tail as an object of the Tail class tail = Tail(topRadius, bottomRadius, length, radius, name) + # Saves position on object for practicality + tail.position = position - # Add tail to aerodynamic surfaces list + # Add tail to aerodynamic surfaces and tail list self.aerodynamicSurfaces.append(aeroSurface=tail, position=position) + self.tail.append(tail) # Refresh static margin calculation self.evaluateStaticMargin() @@ -519,9 +525,12 @@ def addNose(self, length, kind, position, name="Nose Cone"): """ # Create a nose as an object of NoseCone class nose = NoseCone(length, kind, self.radius, self.radius, name) + # Saves position on object for practicality + nose.position = position # Add nose to the list of aerodynamic surfaces self.aerodynamicSurfaces.append(aeroSurface=nose, position=position) + self.nosecone.append(nose) # Refresh static margin calculation self.evaluateStaticMargin() @@ -630,9 +639,12 @@ def addTrapezoidalFins( airfoil, name, ) + # Saves position on object for practicality + finSet.position = position # Add fin set to the list of aerodynamic surfaces self.aerodynamicSurfaces.append(aeroSurface=finSet, position=position) + self.fins.append(finSet) # Refresh static margin calculation self.evaluateStaticMargin() @@ -703,9 +715,12 @@ def addEllipticalFins( # Create a fin set as an object of EllipticalFins class finSet = EllipticalFins(n, rootChord, span, radius, cantAngle, airfoil, name) + # Saves position on object for practicality + finSet.position = position # Add fin set to the list of aerodynamic surfaces self.aerodynamicSurfaces.append(aeroSurface=finSet, position=position) + self.fins.append(finSet) # Refresh static margin calculation self.evaluateStaticMargin() diff --git a/tests/conftest.py b/tests/conftest.py index 650a8a78a..528e0ce2f 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -199,6 +199,38 @@ def linear_func(): ) +@pytest.fixture +def linearly_interpolated_func(): + """Create a linearly interpolated function based on a list of points. + + Returns + ------- + Function + Piece-wise linearly interpolated, with constant extrapolation + """ + return Function( + [[0, 0], [1, 7], [2, -3], [3, -1], [4, 3]], + interpolation="spline", + extrapolation="constant", + ) + + +@pytest.fixture +def spline_interpolated_func(): + """Create a spline interpolated function based on a list of points. + + Returns + ------- + Function + Spline interpolated, with natural extrapolation + """ + return Function( + [[0, 0], [1, 7], [2, -3], [3, -1], [4, 3]], + interpolation="spline", + extrapolation="natural", + ) + + @pytest.fixture def func_from_csv(): func = Function( diff --git a/tests/test_function.py b/tests/test_function.py index 5e56f12c3..3be1dc110 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -134,7 +134,9 @@ def test_extrapolation_methods(linear_func): assert np.isclose(linear_func.getValue(-1), -1, atol=1e-6) -def test_integral(linear_func): +@pytest.mark.parametrize("a", [-1, 0, 0.5, 1, 2, 2.5, 3.5, 4, 5]) +@pytest.mark.parametrize("b", [-1, 0, 0.5, 1, 2, 2.5, 3.5, 4, 5]) +def test_integral_linear_interpolation(linearly_interpolated_func, a, b): """Test the integral method of the Function class. Parameters @@ -143,4 +145,33 @@ def test_integral(linear_func): A Function object created from a list of values. """ # Test integral - assert np.isclose(linear_func.integral(0, 1), 0.5, atol=1e-6) + assert np.isclose( + linearly_interpolated_func.integral(a, b, numerical=False), + linearly_interpolated_func.integral(a, b, numerical=True), + atol=1e-3, + ) + + +@pytest.mark.parametrize("func", ["linear_func", "spline_interpolated_func"]) +@pytest.mark.parametrize("a", [-1, -0.5, 0, 0.5, 1, 2, 2.5, 3.5, 4, 5]) +@pytest.mark.parametrize("b", [-1, -0.5, 0, 0.5, 1, 2, 2.5, 3.5, 4, 5]) +def test_integral_spline_interpolation(request, func, a, b): + """Test the integral method of the Function class. + + Parameters + ---------- + spline_func : rocketpy.Function + A Function object created from a list of values. + a : float + Lower limit of the integral. + b : float + Upper limit of the integral. + """ + # Test integral + # Get the function from the fixture + func = request.getfixturevalue(func) + assert np.isclose( + func.integral(a, b, numerical=False), + func.integral(a, b, numerical=True), + atol=1e-3, + ) diff --git a/tests/test_utilities.py b/tests/test_utilities.py index 1dffd93fd..b595bc8bd 100644 --- a/tests/test_utilities.py +++ b/tests/test_utilities.py @@ -29,8 +29,8 @@ def test_fin_flutter_analysis(mock_show, flight): see_graphs=False, ) - assert abs(flutter_mach(15) - 1.085295573) < 1e-6 - assert abs(safety_factor(15) - 3.373824095) < 1e-6 + assert abs(flutter_mach(15) - 1.085295573) < 1e-3 + assert abs(safety_factor(15) - 3.373824095) < 1e-3 assert ( utilities.fin_flutter_analysis( fin_thickness=2 / 1000,