ENH: nosecone draw and bluffness

rocketpy/AeroSurface.py

@@ -1,14 +1,16 @@
-__author__ = "Guilherme Fernandes Alves, Mateus Stano Junqueira"
+__author__ = "Guilherme Fernandes Alves, Mateus Stano Junqueira, Giovani Hidalgo Ceotto, Franz MasatoshiYuri, Calebe Gomes Teles"
 __copyright__ = "Copyright 20XX, RocketPy Team"
 __license__ = "MIT"
 
 from abc import ABC, abstractmethod
 
 import matplotlib.pyplot as plt
 import numpy as np
-from matplotlib.patches import Ellipse
 
 from .Function import Function
+from abc import ABC, abstractmethod, abstractproperty
+from matplotlib.patches import Ellipse
+from scipy.optimize import fsolve
 
 
 class AeroSurface(ABC):
@@ -96,7 +98,8 @@ class NoseCone(AeroSurface):
     NoseCone.rocketRadius : float
         The reference rocket radius used for lift coefficient normalization, in meters.
     NoseCone.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".
     NoseCone.name : string
         Nose cone name. Has no impact in simulation, as it is only used to
         display data in a more organized matter.
@@ -125,6 +128,7 @@ def __init__(
         length,
         kind,
         baseRadius=None,
+        bluffiness=0,
         rocketRadius=None,
         name="Nose Cone",
     ):
@@ -139,9 +143,13 @@ def __init__(
             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.
+            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
             assumed as 1.
+        bluffiness : float, optional
+            Ratio between the radius of the circle on the tip of the ogive and
+            the radius of the base of the ogive.
         rocketRadius : int, float, optional
             The reference rocket radius used for lift coefficient normalization.
             If not given, the ratio between baseRadius and rocketRadius will be
@@ -159,6 +167,7 @@ def __init__(
         self._rocketRadius = rocketRadius
         self._baseRadius = baseRadius
         self._length = length
+        self.bluffiness = bluffiness
         self.kind = kind
 
         self.evaluateGeometricalParameters()
@@ -202,15 +211,27 @@ def kind(self):
 
     @kind.setter
     def kind(self, value):
-        # Analyze type
+        # Analyzes nosecone type
+        # Sets the k for Cp calculation
+        # Sets the function which creates the respective curve
         self._kind = value
+        value = (value.replace(" ", "")).lower()
+
         if value == "conical":
             self.k = 2 / 3
-        elif value == "ogive":
-            self.k = 0.466
+            self.y_nosecone = Function(lambda x: x * self.baseRadius / self.length)
+
         elif value == "lvhaack":
             self.k = 0.563
-        elif value == "tangent":
+            theta = lambda x: np.arccos(1 - 2 * max(min(x / self.length, 1), 0))
+            self.y_nosecone = Function(
+                lambda x: self.baseRadius
+                * (theta(x) - np.sin(2 * theta(x)) / 2 + (np.sin(theta(x)) ** 3) / 3)
+                ** (0.5)
+                / (np.pi**0.5)
+            )
+
+        elif value in ["tangent", "tangentogive", "ogive"]:
             rho = (self.baseRadius**2 + self.length**2) / (2 * self.baseRadius)
             volume = np.pi * (
                 self.length * rho**2
@@ -219,10 +240,85 @@ def kind(self, value):
             )
             area = np.pi * self.baseRadius**2
             self.k = 1 - volume / (area * self.length)
+            self.y_nosecone = Function(
+                lambda x: np.sqrt(rho**2 - (min(x - self.length, 0)) ** 2)
+                + (self.baseRadius - rho)
+            )
+
         elif value == "elliptical":
             self.k = 1 / 3
+            self.y_nosecone = Function(
+                lambda x: self.baseRadius
+                * np.sqrt(1 - ((x - self.length) / self.length) ** 2)
+            )
+
+        elif value == "vonkarman":
+            self.k = 0.5
+            theta = lambda x: np.arccos(1 - 2 * max(min(x / self.length, 1), 0))
+            self.y_nosecone = Function(
+                lambda x: self.baseRadius
+                * (theta(x) - np.sin(2 * theta(x)) / 2) ** (0.5)
+                / (np.pi**0.5)
+            )
+        elif value == "parabolic":
+            self.k = 0.5
+            self.y_nosecone = Function(
+                lambda x: self.baseRadius
+                * ((2 * x / self.length - (x / self.length) ** 2) / (2 - 1))
+            )
+
         else:
-            self.k = 0.5  # Parabolic and Von Karman
+            raise ValueError(
+                f"Nose Cone kind '{self.kind}' not found, "
+                + "please use one of the following Nose Cone kinds:"
+                + '\n\t"conical"'
+                + '\n\t"ogive"'
+                + '\n\t"lvhaack"'
+                + '\n\t"tangent"'
+                + '\n\t"vonkarman"'
+                + '\n\t"elliptical"'
+                + '\n\t"parabolic"\n'
+            )
+
+        n = 127  # Points on the final curve.
+        p = 3  # Density modifier. Greater n makes more points closer to 0. n=1 -> points equally spaced.
+
+        # Finds the tangential intersection point between the circle and nosecone curve.
+        yPrimeNosecone = lambda x: self.y_nosecone.differentiate(x)
+        xIntercept = lambda x: x + self.y_nosecone(x) * yPrimeNosecone(x)
+        radius = lambda x: (self.y_nosecone(x) ** 2 + (x - xIntercept(x)) ** 2) ** 0.5
+        xInit = fsolve(
+            lambda x: radius(x) - self.bluffiness * self.baseRadius if x > 3e-7 else 0,
+            self.bluffiness * self.baseRadius,
+            xtol=1e-7,
+        )[0]
+
+        # Corrects circle radius if it's too small.
+        if xInit > 0:
+            r = self.bluffiness * self.baseRadius
+        else:
+            r = 0
+            print(
+                "ATTENTION: The chosen bluffiness ratio is insufficient for the selected nosecone category, thereby the effective bluffiness will be 0."
+            )
+
+        # Creates the circle at correct position.
+        circleCenter = xIntercept(xInit)
+        circle = lambda x: abs(r**2 - (x - circleCenter) ** 2) ** 0.5
+
+        # Function defining final shape of curve with circle o the tip.
+        finalShape = Function(lambda x: self.y_nosecone(x) if x >= xInit else circle(x))
+        finalShapeVec = np.vectorize(finalShape)
+
+        # Creates the vectors X and Y with the points of the curve.
+        self.nosecone_Xs = (self.length - (circleCenter - r)) * (
+            np.linspace(0, 1, n) ** p
+        )
+        self.nosecone_Ys = finalShapeVec(self.nosecone_Xs + (circleCenter - r))
+
+        # Evaluates final geometry parameters.
+        self.length = self.nosecone_Xs[-1]
+        self.FinenessRatio = self.length / (2 * self.baseRadius)
         self.evaluateCenterOfPressure()
 
     def evaluateGeometricalParameters(self):
@@ -291,6 +387,62 @@ def evaluateCenterOfPressure(self):
         self.cp = (self.cpx, self.cpy, self.cpz)
         return self.cp
 
+    def draw(self):
+        # Figure creation and set up
+        fig_Ogive, ax = plt.subplots()
+        ax.set_xlim(-0.05, self.length * 1.02)  # Horizontal size
+        ax.set_ylim(-self.baseRadius * 1.05, self.baseRadius * 1.05)  # Vertical size
+        ax.set_aspect("equal")  # Makes the graduation be the same on both axis
+        ax.set_facecolor("#EEEEEE")  # Background colour
+        ax.grid(True, linestyle="--", linewidth=0.5)
+
+        cp_plot = (self.cpz, 0)
+        # Plotting
+        ax.plot(
+            self.nosecone_Xs, self.nosecone_Ys, linestyle="-", color="#A60628"
+        )  # Ogive's upper side
+        ax.plot(
+            self.nosecone_Xs, -self.nosecone_Ys, linestyle="-", color="#A60628"
+        )  # Ogive's lower side
+        ax.scatter(
+            *cp_plot, label="Center Of Pressure", color="red", s=100, zorder=10
+        )  # Center of pressure inner circle
+        ax.scatter(
+            *cp_plot, facecolors="none", edgecolors="red", s=500, zorder=10
+        )  # Center of pressure outer circle
+        # Center Line
+        ax.plot(
+            [0, self.nosecone_Xs[len(self.nosecone_Xs) - 1]],
+            [0, 0],
+            linestyle="--",
+            color="#7A68A6",
+            linewidth=1.5,
+            label="Center Line",
+        )
+        # Vertical base line
+        ax.plot(
+            [
+                self.nosecone_Xs[len(self.nosecone_Xs) - 1],
+                self.nosecone_Xs[len(self.nosecone_Xs) - 1],
+            ],
+            [
+                self.nosecone_Ys[len(self.nosecone_Ys) - 1],
+                -self.nosecone_Ys[len(self.nosecone_Ys) - 1],
+            ],
+            linestyle="-",
+            color="#A60628",
+            linewidth=1.5,
+        )
+
+        # Labels and legend
+        ax.set_xlabel("Length")
+        ax.set_ylabel("Radius")
+        ax.set_title(self.kind + " Nose Cone")
+        ax.legend(bbox_to_anchor=(1, -0.2))
+        # Show Plot
+        plt.show()
+        return None
+
     def geometricalInfo(self):
         """Prints out all the geometric information of the nose cone.
 

rocketpy/Rocket.py

@@ -543,7 +543,7 @@ def addTail(
         # Return self
         return tail
 
-    def addNose(self, length, kind, position, name="Nosecone"):
+    def addNose(self, length, kind, position, bluffiness=0, name="Nose Cone"):
         """Creates a nose cone, storing its parameters as part of the
         aerodynamicSurfaces list. Its parameters are the axial position
         along the rocket and its derivative of the coefficient of lift
@@ -561,6 +561,9 @@ def addNose(self, length, kind, position, name="Nosecone"):
         position : int, float
             Nose cone tip coordinate relative to the rocket's coordinate system.
             See `Rocket.coordinateSystemOrientation` for more information.
+        bluffiness : float, optional
+            Ratio between the radius of the circle on the tip of the ogive and
+            the radius of the base of the ogive.
         name : string
             Nose cone name. Default is "Nose Cone".
 
@@ -570,7 +573,14 @@ def addNose(self, length, kind, position, name="Nosecone"):
             Nose cone object created.
         """
         # Create a nose as an object of NoseCone class
-        nose = NoseCone(length, kind, self.radius, self.radius, name)
+        nose = NoseCone(
+            length=length,
+            kind=kind,
+            baseRadius=self.radius,
+            rocketRadius=self.radius,
+            bluffiness=bluffiness,
+            name=name,
+        )
 
         # Add nose to the list of aerodynamic surfaces
         self.addSurfaces(nose, position)

tests/test_rocket.py

@@ -370,13 +370,7 @@ def test_evaluate_static_margin_assert_cp_equals_cm(kg, m, dimensionless_rocket)
 
 @pytest.mark.parametrize(
     "k, type",
-    (
-        [2 / 3, "conical"],
-        [0.466, "ogive"],
-        [0.563, "lvhaack"],
-        [0.5, "default"],
-        [0.5, "not a mapped string, to show default case"],
-    ),
+    ([2 / 3, "conical"], [0.46469957130675876, "ogive"], [0.563, "lvhaack"]),
 )
 def test_add_nose_assert_cp_cm_plus_nose(k, type, rocket, dimensionless_rocket, m):
     rocket.addNose(length=0.55829, kind=type, position=1.278)