From 0c8b389d6679ecd7cc99d7af3dfe3ec325920734 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Wed, 18 Sep 2024 01:14:25 +0200
Subject: [PATCH 01/18] ENH: add free form fins file

---
 .../aero_surface/fins/free_form_fins.py       | 367 ++++++++++++++++++
 1 file changed, 367 insertions(+)
 create mode 100644 rocketpy/rocket/aero_surface/fins/free_form_fins.py

diff --git a/rocketpy/rocket/aero_surface/fins/free_form_fins.py b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
new file mode 100644
index 000000000..007bf1c5c
--- /dev/null
+++ b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
@@ -0,0 +1,367 @@
+import warnings
+import numpy as np
+
+from rocketpy.plots.aero_surface_plots import _FreeFormFinsPlots
+from rocketpy.prints.aero_surface_prints import _FreeFormFinsPrints
+
+from .fins import Fins
+
+
+class FreeFormFins(Fins):
+    """Class that defines and holds information for a free form fin set.
+
+    This class inherits from the Fins class.
+
+    Note
+    ----
+    Local coordinate system:
+        - Origin located at the top of the root chord.
+        - Z axis along the longitudinal axis of symmetry, positive downwards (top -> bottom).
+        - Y axis perpendicular to the Z axis, in the span direction, positive upwards.
+        - X axis completes the right-handed coordinate system.
+
+    See Also
+    --------
+    Fins
+
+    Attributes
+    ----------
+    FreeFormFins.n : int
+        Number of fins in fin set.
+    FreeFormFins.rocket_radius : float
+        The reference rocket radius used for lift coefficient normalization, in
+        meters.
+    FreeFormFins.airfoil : tuple
+        Tuple of two items. First is the airfoil lift curve.
+        Second is the unit of the curve (radians or degrees).
+    FreeFormFins.cant_angle : float
+        Fins cant angle with respect to the rocket centerline, in degrees.
+    FreeFormFins.changing_attribute_dict : dict
+        Dictionary that stores the name and the values of the attributes that
+        may be changed during a simulation. Useful for control systems.
+    FreeFormFins.cant_angle_rad : float
+        Fins cant angle with respect to the rocket centerline, in radians.
+    FreeFormFins.root_chord : float
+        Fin root chord in meters.
+    FreeFormFins.tip_chord : float
+        Fin tip chord in meters.
+    FreeFormFins.span : float
+        Fin span in meters.
+    FreeFormFins.name : string
+        Name of fin set.
+    FreeFormFins.sweep_length : float
+        Fins sweep length in meters. By sweep length, understand the axial
+        distance between the fin root leading edge and the fin tip leading edge
+        measured parallel to the rocket centerline.
+    FreeFormFins.sweep_angle : float
+        Fins sweep angle with respect to the rocket centerline. Must
+        be given in degrees.
+    FreeFormFins.d : float
+        Reference diameter of the rocket, in meters.
+    FreeFormFins.ref_area : float
+        Reference area of the rocket, in m².
+    FreeFormFins.Af : float
+        Area of the longitudinal section of each fin in the set.
+    FreeFormFins.AR : float
+        Aspect ratio of each fin in the set
+    FreeFormFins.gamma_c : float
+        Fin mid-chord sweep angle.
+    FreeFormFins.Yma : float
+        Span wise position of the mean aerodynamic chord.
+    FreeFormFins.roll_geometrical_constant : float
+        Geometrical constant used in roll calculations.
+    FreeFormFins.tau : float
+        Geometrical relation used to simplify lift and roll calculations.
+    FreeFormFins.lift_interference_factor : float
+        Factor of Fin-Body interference in the lift coefficient.
+    FreeFormFins.cp : tuple
+        Tuple with the x, y and z local coordinates of the fin set center of
+        pressure. Has units of length and is given in meters.
+    FreeFormFins.cpx : float
+        Fin set local center of pressure x coordinate. Has units of length and
+        is given in meters.
+    FreeFormFins.cpy : float
+        Fin set local center of pressure y coordinate. Has units of length and
+        is given in meters.
+    FreeFormFins.cpz : float
+        Fin set local center of pressure z coordinate. Has units of length and
+        is given in meters.
+    FreeFormFins.cl : Function
+        Function which defines the lift coefficient as a function of the angle
+        of attack and the Mach number. Takes as input the angle of attack in
+        radians and the Mach number. Returns the lift coefficient.
+    FreeFormFins.clalpha : float
+        Lift coefficient slope. Has units of 1/rad.
+    FreeFormFins.mac_length : float
+        Mean aerodynamic chord length of the fin set.
+    FreeFormFins.mac_lead : float
+        Mean aerodynamic chord leading edge x coordinate.
+    """
+
+    def __init__(
+        self,
+        n,
+        shape_points,
+        rocket_radius,
+        cant_angle=0,
+        airfoil=None,
+        name="Fins",
+    ):
+        """Initialize FreeFormFins class.
+
+        Parameters
+        ----------
+        n : int
+            Number of fins, from 2 to infinity.
+        shape_points : list
+            List of tuples (x, y) containing the coordinates of the fin's
+            geometry defining points. The point (0, 0) is the root leading edge.
+            Positive x is rearwards, positive y is upwards (span direction).
+        rocket_radius : int, float
+            Reference radius to calculate lift coefficient, in meters.
+        cant_angle : int, float, optional
+            Fins cant angle with respect to the rocket centerline. Must
+            be given in degrees.
+        airfoil : tuple, optional
+            Default is null, in which case fins will be treated as flat plates.
+            Otherwise, if tuple, fins will be considered as airfoils. The
+            tuple's first item specifies the airfoil's lift coefficient
+            by angle of attack and must be either a .csv, .txt, ndarray
+            or callable. The .csv and .txt files can contain a single line
+            header and the first column must specify the angle of attack, while
+            the second column must specify the lift coefficient. The
+            ndarray should be as [(x0, y0), (x1, y1), (x2, y2), ...]
+            where x0 is the angle of attack and y0 is the lift coefficient.
+            If callable, it should take an angle of attack as input and
+            return the lift coefficient at that angle of attack.
+            The tuple's second item is the unit of the angle of attack,
+            accepting either "radians" or "degrees".
+        name : str
+            Name of fin set.
+
+        Returns
+        -------
+        None
+        """
+        self.shape_points = shape_points
+
+        down = False
+        for i in range(1, len(shape_points)):
+            if shape_points[i][1] > shape_points[i - 1][1] and down:
+                warnings.warn(
+                    "Jagged fin shape detected. This may cause small inaccuracies "
+                    "center of pressure and pitch moment calculations."
+                )
+                break
+            if shape_points[i][1] < shape_points[i - 1][1]:
+                down = True
+            i += 1
+
+        root_chord = abs(shape_points[0][0] - shape_points[-1][0])
+        ys = [point[1] for point in shape_points]
+        span = max(ys) - min(ys)
+
+        super().__init__(
+            n,
+            root_chord,
+            span,
+            rocket_radius,
+            cant_angle,
+            airfoil,
+            name,
+        )
+
+        self.evaluate_geometrical_parameters()
+        self.evaluate_center_of_pressure()
+        self.evaluate_lift_coefficient()
+        self.evaluate_roll_parameters()
+
+        self.prints = _FreeFormFinsPrints(self)
+        self.plots = _FreeFormFinsPlots(self)
+
+    def evaluate_center_of_pressure(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.
+
+        Returns
+        -------
+        None
+        """
+        # Center of pressure position in local coordinates
+        cpz = self.mac_lead + 0.25 * self.mac_length
+        self.cpx = 0
+        self.cpy = 0
+        self.cpz = cpz
+        self.cp = (self.cpx, self.cpy, self.cpz)
+
+    def evaluate_geometrical_parameters(self):  # pylint: disable=too-many-statements
+        """Calculates and saves fin set's geometrical parameters such as the
+        fins' area, aspect ratio and parameters for roll movement. The
+        calculations related to the free form fins points are done in the exact
+        same way as Open Rocket.
+
+        Returns
+        -------
+        None
+        """
+        # pylint: disable=invalid-name
+        Af = 0
+        for i in range(len(self.shape_points) - 1):
+            Af += (self.shape_points[i][1] + self.shape_points[i + 1][1]) * (
+                self.shape_points[i][0] - self.shape_points[i + 1][0]
+            )
+        Af = abs(Af) / 2
+        if Af < 1e-6:
+            raise ValueError("Fin area is too small. Check the shape_points.")
+
+        AR = 2 * self.span**2 / Af  # Fin aspect ratio
+        tau = (self.span + self.rocket_radius) / self.rocket_radius
+        lift_interference_factor = 1 + 1 / tau
+        roll_forcing_interference_factor = (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))
+        )
+
+        # Discretize fin shape interpolation into points_per_line points
+        points_per_line = 40  # Arbitrary number of points per line, same as OPR
+
+        chord_lead = np.ones(points_per_line) * np.inf  # Leading edge of the chord
+        chord_trail = np.ones(points_per_line) * -np.inf  # Trailing edge of the chord
+        chord_length = np.zeros(points_per_line)  # Chord length
+
+        # Calculate chord length and leading/trailing edge
+        # If fin is jagged, calculation will be for the equivalent shape
+        for p in range(1, len(self.shape_points)):
+            # Coordinates of the two points
+            x1 = np.float64(self.shape_points[p - 1][0])
+            y1 = np.float64(self.shape_points[p - 1][1])
+            x2 = np.float64(self.shape_points[p][0])
+            y2 = np.float64(self.shape_points[p][1])
+
+            # Correction for jagged fin
+            previous_point = int(y1 * 1.0001 / self.span * (points_per_line - 1))
+            current_point = int(y2 * 1.0001 / self.span * (points_per_line - 1))
+            previous_point = max(min(previous_point, points_per_line - 1), 0)
+            current_point = max(min(current_point, points_per_line - 1), 0)
+            if previous_point > current_point:
+                previous_point, current_point = current_point, previous_point
+
+            # Calculate chord length and leading/trailing edge
+            for i in range(previous_point, current_point + 1):
+                y = i * self.span / (points_per_line - 1)
+                x = max(
+                    min(
+                        (y - y2) / (y1 - y2) * x1 + (y1 - y) / (y1 - y2) * x2,
+                        max(x1, x2),
+                    ),
+                    min(x1, x2),
+                )
+                if x < chord_lead[i]:
+                    chord_lead[i] = x
+                if x > chord_trail[i]:
+                    chord_trail[i] = x
+
+                if y1 < y2:
+                    chord_length[i] -= x
+                else:
+                    chord_length[i] += x
+
+        # Remove infinities and nans for cases where there are equal coordinates
+        for i in range(points_per_line):
+            if (
+                np.isinf(chord_lead[i])
+                or np.isinf(chord_trail[i])
+                or np.isnan(chord_lead[i])
+                or np.isnan(chord_trail[i])
+            ):
+                chord_lead[i] = 0
+                chord_trail[i] = 0
+            if chord_length[i] < 0 or np.isnan(chord_length[i]):
+                chord_length[i] = 0
+            if chord_length[i] > chord_trail[i] - chord_lead[i]:
+                chord_length[i] = chord_trail[i] - chord_lead[i]
+
+        # Initialize integration variables
+        radius = self.rocket_radius  # Rocket radius
+        area = 0  # Area of the fin, this will be different from Af if fin is jagged
+        mac_length = 0  # Mean aerodynamic chord length)
+        mac_lead = 0  # Mean aerodynamic chord leading edge x coordinate
+        mac_span = 0  # Mean aerodynamic chord span wise position (Yma)
+        cos_gamma = 0  # Cosine of the sweep angle
+        roll_geometrical_constant = 0
+        roll_damping_interference_factor_numerator = 0
+        roll_damping_interference_factor_denominator = 0
+
+        # Perform integration
+        dy = self.span / (points_per_line - 1)
+        for i in range(points_per_line):
+            length = chord_trail[i] - chord_lead[i]
+            y = i * dy
+
+            mac_length += length * length
+            mac_span += y * length
+            mac_lead += chord_lead[i] * length
+            area += length
+            roll_geometrical_constant += chord_length[i] * (radius + y) ** 2
+            roll_damping_interference_factor_numerator += (
+                radius**3 * length / (radius + y) ** 2
+            )
+            roll_damping_interference_factor_denominator += (radius + y) * length
+
+            if i > 0:
+                dx = (chord_trail[i] + chord_lead[i]) / 2 - (
+                    chord_trail[i - 1] + chord_lead[i - 1]
+                ) / 2
+                cos_gamma += dy / (dx**2 + dy**2) ** 0.5
+
+        mac_length *= dy
+        mac_span *= dy
+        mac_lead *= dy
+        area *= dy
+        roll_geometrical_constant *= dy
+        roll_damping_interference_factor_numerator *= dy
+        roll_damping_interference_factor_denominator *= dy
+        mac_length /= area
+        mac_span /= area
+        mac_lead /= area
+        cos_gamma /= points_per_line - 1
+
+        # Store values
+        self.Af = Af  # Fin area
+        self.AR = AR  # Aspect Ratio
+        self.gamma_c = np.arccos(cos_gamma)  # Mid chord angle
+        self.Yma = mac_span  # Span wise coord of mean aero chord
+        self.mac_length = mac_length
+        self.mac_lead = mac_lead
+        self.tau = tau
+        self.roll_geometrical_constant = roll_geometrical_constant
+        self.lift_interference_factor = lift_interference_factor
+        self.roll_forcing_interference_factor = roll_forcing_interference_factor
+        self.roll_damping_interference_factor = 1 + (
+            roll_damping_interference_factor_numerator
+            / roll_damping_interference_factor_denominator
+        )
+
+        self.evaluate_shape()
+
+    def evaluate_shape(self):
+        x_array, y_array = zip(*self.shape_points)
+        self.shape_vec = [np.array(x_array), np.array(y_array)]
+
+    def info(self):
+        self.prints.geometry()
+        self.prints.lift()
+
+    def all_info(self):
+        self.prints.all()
+        self.plots.all()

From fe31262023b95e8b76f6066f2ac372f79c77ce71 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Wed, 18 Sep 2024 01:15:07 +0200
Subject: [PATCH 02/18] ENH: free form fins plots and prints

---
 rocketpy/plots/aero_surface_plots.py   | 61 ++++++++++++++++++++++++++
 rocketpy/prints/aero_surface_prints.py |  4 ++
 2 files changed, 65 insertions(+)

diff --git a/rocketpy/plots/aero_surface_plots.py b/rocketpy/plots/aero_surface_plots.py
index 6632f2de6..85f27c123 100644
--- a/rocketpy/plots/aero_surface_plots.py
+++ b/rocketpy/plots/aero_surface_plots.py
@@ -380,6 +380,67 @@ def draw(self):
         plt.show()
 
 
+class _FreeFormFinsPlots(_FinsPlots):
+    """Class that contains all free form fin plots."""
+
+    # pylint: disable=too-many-statements
+    def draw(self):
+        """Draw the fin shape along with some important information, including
+        the center line, the quarter line and the center of pressure position.
+
+        Returns
+        -------
+        None
+        """
+        # Color cycle [#348ABD, #A60628, #7A68A6, #467821, #D55E00, #CC79A7,
+        # #56B4E9, #009E73, #F0E442, #0072B2]
+
+        # Center of pressure
+        cp_point = [self.aero_surface.cpz, self.aero_surface.Yma]
+
+        # Mean Aerodynamic Chord
+        yma_line = plt.Line2D(
+            (
+                self.aero_surface.mac_lead,
+                self.aero_surface.mac_lead + self.aero_surface.mac_length,
+            ),
+            (self.aero_surface.Yma, self.aero_surface.Yma),
+            color="#467821",
+            linestyle="--",
+            label="Mean Aerodynamic Chord",
+        )
+
+        # Plotting
+        fig = plt.figure(figsize=(7, 4))
+        with plt.style.context("bmh"):
+            ax = fig.add_subplot(111)
+
+        # Fin
+        ax.scatter(
+            self.aero_surface.shape_vec[0],
+            self.aero_surface.shape_vec[1],
+            color="#A60628",
+        )
+        ax.plot(
+            self.aero_surface.shape_vec[0],
+            self.aero_surface.shape_vec[1],
+            color="#A60628",
+        )
+
+        ax.add_line(yma_line)
+        ax.scatter(*cp_point, label="Center of Pressure", color="red", s=100, zorder=10)
+        ax.scatter(*cp_point, facecolors="none", edgecolors="red", s=500, zorder=10)
+
+        # Plot settings
+        ax.set_xlabel("Root chord (m)")
+        ax.set_ylabel("Span (m)")
+        ax.set_title("Trapezoidal Fin Cross Section")
+        ax.legend(bbox_to_anchor=(1.05, 1.0), loc="upper left")
+
+        plt.tight_layout()
+        plt.show()
+
+
 class _TailPlots(_AeroSurfacePlots):
     """Class that contains all tail plots."""
 
diff --git a/rocketpy/prints/aero_surface_prints.py b/rocketpy/prints/aero_surface_prints.py
index 45b597ccb..378433511 100644
--- a/rocketpy/prints/aero_surface_prints.py
+++ b/rocketpy/prints/aero_surface_prints.py
@@ -183,6 +183,10 @@ class _EllipticalFinsPrints(_FinsPrints):
     """Class that contains all elliptical fins prints."""
 
 
+class _FreeFormFinsPrints(_FinsPrints):
+    """Class that contains all free form fins prints."""
+
+
 class _TailPrints(_AeroSurfacePrints):
     """Class that contains all tail prints."""
 

From 9a02ef8b49d6611e9e1095e1caa29681fc322220 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Wed, 18 Sep 2024 01:15:19 +0200
Subject: [PATCH 03/18] ENH: add free form fins to rocket

---
 rocketpy/rocket/rocket.py | 78 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 78 insertions(+)

diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py
index 310ca7e38..a5c975d87 100644
--- a/rocketpy/rocket/rocket.py
+++ b/rocketpy/rocket/rocket.py
@@ -18,6 +18,7 @@
     Tail,
     TrapezoidalFins,
 )
+from rocketpy.rocket.aero_surface.fins.free_form_fins import FreeFormFins
 from rocketpy.rocket.aero_surface.generic_surface import GenericSurface
 from rocketpy.rocket.components import Components
 from rocketpy.rocket.parachute import Parachute
@@ -1322,6 +1323,83 @@ def add_elliptical_fins(
         self.add_surfaces(fin_set, position)
         return fin_set
 
+    def add_free_form_fins(
+        self,
+        n,
+        shape_points,
+        position,
+        cant_angle=0.0,
+        radius=None,
+        airfoil=None,
+        name="Fins",
+    ):
+        """Create a free form fin set, storing its parameters as part of the
+        aerodynamic_surfaces list. Its parameters are the axial position along
+        the rocket and its derivative of the coefficient of lift in respect to
+        angle of attack.
+
+        Parameters
+        ----------
+        n : int
+            Number of fins, from 2 to infinity.
+        shape_points : list
+            List of tuples (x, y) containing the coordinates of the fin's
+            geometry defining points. The point (0, 0) is the root leading edge.
+            Positive x is rearwards, positive y is upwards (span direction).
+        position : int, float
+            Fin set position relative to the rocket's coordinate system.
+            By fin set position, understand the point belonging to the root
+            chord which is highest in the rocket coordinate system (i.e.
+            the point closest to the nose cone tip).
+
+            See Also
+            --------
+            :ref:`positions`
+        cant_angle : int, float, optional
+            Fins cant angle with respect to the rocket centerline. Must
+            be given in degrees.
+        radius : int, float, optional
+            Reference fuselage radius where the fins are located. This is used
+            to calculate lift coefficient and to draw the rocket. If None,
+            which is default, the rocket radius will be used.
+        airfoil : tuple, optional
+            Default is null, in which case fins will be treated as flat plates.
+            Otherwise, if tuple, fins will be considered as airfoils. The
+            tuple's first item specifies the airfoil's lift coefficient
+            by angle of attack and must be either a .csv, .txt, ndarray
+            or callable. The .csv and .txt files can contain a single line
+            header and the first column must specify the angle of attack, while
+            the second column must specify the lift coefficient. The
+            ndarray should be as [(x0, y0), (x1, y1), (x2, y2), ...]
+            where x0 is the angle of attack and y0 is the lift coefficient.
+            If callable, it should take an angle of attack as input and
+            return the lift coefficient at that angle of attack.
+            The tuple's second item is the unit of the angle of attack,
+            accepting either "radians" or "degrees".
+
+        Returns
+        -------
+        fin_set : FreeFormFins
+            Fin set object created.
+        """
+
+        # Modify radius if not given, use rocket radius, otherwise use given.
+        radius = radius if radius is not None else self.radius
+
+        # Create a fin set as an object of TrapezoidalFins class
+        fin_set = FreeFormFins(
+            n,
+            shape_points,
+            radius,
+            cant_angle,
+            airfoil,
+            name,
+        )
+
+        # Add fin set to the list of aerodynamic surfaces
+        self.add_surfaces(fin_set, position)
+        return fin_set
+
     def add_parachute(
         self, name, cd_s, trigger, sampling_rate=100, lag=0, noise=(0, 0, 0)
     ):

From 3d2880ff6cd03de2b442b941aec231a47e0452eb Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Wed, 18 Sep 2024 01:15:35 +0200
Subject: [PATCH 04/18] ENH: update init files

---
 rocketpy/__init__.py                          | 1 +
 rocketpy/rocket/__init__.py                   | 1 +
 rocketpy/rocket/aero_surface/__init__.py      | 7 ++++++-
 rocketpy/rocket/aero_surface/fins/__init__.py | 1 +
 4 files changed, 9 insertions(+), 1 deletion(-)

diff --git a/rocketpy/__init__.py b/rocketpy/__init__.py
index 0ada37891..265692659 100644
--- a/rocketpy/__init__.py
+++ b/rocketpy/__init__.py
@@ -38,6 +38,7 @@
     Rocket,
     Tail,
     TrapezoidalFins,
+    FreeFormFins,
 )
 from .sensitivity import SensitivityModel
 from .sensors import Accelerometer, Barometer, GnssReceiver, Gyroscope
diff --git a/rocketpy/rocket/__init__.py b/rocketpy/rocket/__init__.py
index 7000d9108..79b381dd2 100644
--- a/rocketpy/rocket/__init__.py
+++ b/rocketpy/rocket/__init__.py
@@ -10,6 +10,7 @@
     RailButtons,
     Tail,
     TrapezoidalFins,
+    FreeFormFins,
 )
 from rocketpy.rocket.components import Components
 from rocketpy.rocket.parachute import Parachute
diff --git a/rocketpy/rocket/aero_surface/__init__.py b/rocketpy/rocket/aero_surface/__init__.py
index 7cfab65fc..c9e43e6a0 100644
--- a/rocketpy/rocket/aero_surface/__init__.py
+++ b/rocketpy/rocket/aero_surface/__init__.py
@@ -1,6 +1,11 @@
 from rocketpy.rocket.aero_surface.aero_surface import AeroSurface
 from rocketpy.rocket.aero_surface.air_brakes import AirBrakes
-from rocketpy.rocket.aero_surface.fins import EllipticalFins, Fins, TrapezoidalFins
+from rocketpy.rocket.aero_surface.fins import (
+    EllipticalFins,
+    Fins,
+    TrapezoidalFins,
+    FreeFormFins,
+)
 from rocketpy.rocket.aero_surface.generic_surface import GenericSurface
 from rocketpy.rocket.aero_surface.linear_generic_surface import LinearGenericSurface
 from rocketpy.rocket.aero_surface.nose_cone import NoseCone
diff --git a/rocketpy/rocket/aero_surface/fins/__init__.py b/rocketpy/rocket/aero_surface/fins/__init__.py
index f1efc603a..b8c9afd92 100644
--- a/rocketpy/rocket/aero_surface/fins/__init__.py
+++ b/rocketpy/rocket/aero_surface/fins/__init__.py
@@ -1,3 +1,4 @@
 from rocketpy.rocket.aero_surface.fins.elliptical_fins import EllipticalFins
 from rocketpy.rocket.aero_surface.fins.fins import Fins
 from rocketpy.rocket.aero_surface.fins.trapezoidal_fins import TrapezoidalFins
+from rocketpy.rocket.aero_surface.fins.free_form_fins import FreeFormFins

From 384531611d04075290eadaa5e44622e7bcc8de95 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Wed, 18 Sep 2024 01:16:29 +0200
Subject: [PATCH 05/18] MNT: isort

---
 rocketpy/__init__.py                                | 2 +-
 rocketpy/rocket/__init__.py                         | 2 +-
 rocketpy/rocket/aero_surface/__init__.py            | 2 +-
 rocketpy/rocket/aero_surface/fins/__init__.py       | 2 +-
 rocketpy/rocket/aero_surface/fins/free_form_fins.py | 1 +
 5 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/rocketpy/__init__.py b/rocketpy/__init__.py
index 265692659..539b8b2cb 100644
--- a/rocketpy/__init__.py
+++ b/rocketpy/__init__.py
@@ -30,6 +30,7 @@
     Components,
     EllipticalFins,
     Fins,
+    FreeFormFins,
     GenericSurface,
     LinearGenericSurface,
     NoseCone,
@@ -38,7 +39,6 @@
     Rocket,
     Tail,
     TrapezoidalFins,
-    FreeFormFins,
 )
 from .sensitivity import SensitivityModel
 from .sensors import Accelerometer, Barometer, GnssReceiver, Gyroscope
diff --git a/rocketpy/rocket/__init__.py b/rocketpy/rocket/__init__.py
index 79b381dd2..0687b5ee5 100644
--- a/rocketpy/rocket/__init__.py
+++ b/rocketpy/rocket/__init__.py
@@ -4,13 +4,13 @@
     AirBrakes,
     EllipticalFins,
     Fins,
+    FreeFormFins,
     GenericSurface,
     LinearGenericSurface,
     NoseCone,
     RailButtons,
     Tail,
     TrapezoidalFins,
-    FreeFormFins,
 )
 from rocketpy.rocket.components import Components
 from rocketpy.rocket.parachute import Parachute
diff --git a/rocketpy/rocket/aero_surface/__init__.py b/rocketpy/rocket/aero_surface/__init__.py
index c9e43e6a0..ad784f8d0 100644
--- a/rocketpy/rocket/aero_surface/__init__.py
+++ b/rocketpy/rocket/aero_surface/__init__.py
@@ -3,8 +3,8 @@
 from rocketpy.rocket.aero_surface.fins import (
     EllipticalFins,
     Fins,
-    TrapezoidalFins,
     FreeFormFins,
+    TrapezoidalFins,
 )
 from rocketpy.rocket.aero_surface.generic_surface import GenericSurface
 from rocketpy.rocket.aero_surface.linear_generic_surface import LinearGenericSurface
diff --git a/rocketpy/rocket/aero_surface/fins/__init__.py b/rocketpy/rocket/aero_surface/fins/__init__.py
index b8c9afd92..941aa5465 100644
--- a/rocketpy/rocket/aero_surface/fins/__init__.py
+++ b/rocketpy/rocket/aero_surface/fins/__init__.py
@@ -1,4 +1,4 @@
 from rocketpy.rocket.aero_surface.fins.elliptical_fins import EllipticalFins
 from rocketpy.rocket.aero_surface.fins.fins import Fins
-from rocketpy.rocket.aero_surface.fins.trapezoidal_fins import TrapezoidalFins
 from rocketpy.rocket.aero_surface.fins.free_form_fins import FreeFormFins
+from rocketpy.rocket.aero_surface.fins.trapezoidal_fins import TrapezoidalFins
diff --git a/rocketpy/rocket/aero_surface/fins/free_form_fins.py b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
index 007bf1c5c..3aeb13d7e 100644
--- a/rocketpy/rocket/aero_surface/fins/free_form_fins.py
+++ b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
@@ -1,4 +1,5 @@
 import warnings
+
 import numpy as np
 
 from rocketpy.plots.aero_surface_plots import _FreeFormFinsPlots

From 2869e91399e0d170298818252c57f46bf13a48c0 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Wed, 18 Sep 2024 01:22:57 +0200
Subject: [PATCH 06/18] MNT: improve overall evaluate_geom

---
 .../aero_surface/fins/free_form_fins.py       | 162 +++++++++---------
 1 file changed, 84 insertions(+), 78 deletions(-)

diff --git a/rocketpy/rocket/aero_surface/fins/free_form_fins.py b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
index 3aeb13d7e..0b24bf76e 100644
--- a/rocketpy/rocket/aero_surface/fins/free_form_fins.py
+++ b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
@@ -197,35 +197,37 @@ def evaluate_center_of_pressure(self):
         self.cp = (self.cpx, self.cpy, self.cpz)
 
     def evaluate_geometrical_parameters(self):  # pylint: disable=too-many-statements
-        """Calculates and saves fin set's geometrical parameters such as the
-        fins' area, aspect ratio and parameters for roll movement. The
-        calculations related to the free form fins points are done in the exact
-        same way as Open Rocket.
+        """
+        Calculates and saves the fin set's geometrical parameters such as the
+        fin area, aspect ratio, and parameters related to roll movement. This method
+        uses similar calculations to those in OpenRocket for free-form fin shapes.
 
         Returns
         -------
         None
         """
-        # pylint: disable=invalid-name
+        # Calculate the fin area (Af) using the Shoelace theorem (polygon area formula)
         Af = 0
         for i in range(len(self.shape_points) - 1):
-            Af += (self.shape_points[i][1] + self.shape_points[i + 1][1]) * (
-                self.shape_points[i][0] - self.shape_points[i + 1][0]
-            )
+            x1, y1 = self.shape_points[i]
+            x2, y2 = self.shape_points[i + 1]
+            Af += (y1 + y2) * (x1 - x2)
         Af = abs(Af) / 2
         if Af < 1e-6:
             raise ValueError("Fin area is too small. Check the shape_points.")
 
-        AR = 2 * self.span**2 / Af  # Fin aspect ratio
+        # Calculate aspect ratio (AR) and lift interference factors
+        AR = 2 * self.span**2 / Af  # Aspect ratio
         tau = (self.span + self.rocket_radius) / self.rocket_radius
         lift_interference_factor = 1 + 1 / tau
+
+        # Calculate roll forcing interference factor using OpenRocket's approach
         roll_forcing_interference_factor = (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)
+            + ((tau**2 + 1) ** 2 / (tau**2 * (tau - 1) ** 2))
             * (np.arcsin((tau**2 - 1) / (tau**2 + 1))) ** 2
             - (4 * (tau + 1))
             / (tau * (tau - 1))
@@ -233,51 +235,53 @@ def evaluate_geometrical_parameters(self):  # pylint: disable=too-many-statement
             + (8 / (tau - 1) ** 2) * np.log((tau**2 + 1) / (2 * tau))
         )
 
-        # Discretize fin shape interpolation into points_per_line points
-        points_per_line = 40  # Arbitrary number of points per line, same as OPR
+        # Define number of interpolation points along the span of the fin
+        points_per_line = 40  # Same as OpenRocket
 
-        chord_lead = np.ones(points_per_line) * np.inf  # Leading edge of the chord
-        chord_trail = np.ones(points_per_line) * -np.inf  # Trailing edge of the chord
-        chord_length = np.zeros(points_per_line)  # Chord length
+        # Initialize arrays for leading/trailing edge and chord lengths
+        chord_lead = np.ones(points_per_line) * np.inf  # Leading edge x coordinates
+        chord_trail = np.ones(points_per_line) * -np.inf  # Trailing edge x coordinates
+        chord_length = np.zeros(
+            points_per_line
+        )  # Chord length for each spanwise section
 
-        # Calculate chord length and leading/trailing edge
-        # If fin is jagged, calculation will be for the equivalent shape
+        # Discretize fin shape and calculate chord length, leading, and trailing edges
         for p in range(1, len(self.shape_points)):
-            # Coordinates of the two points
-            x1 = np.float64(self.shape_points[p - 1][0])
-            y1 = np.float64(self.shape_points[p - 1][1])
-            x2 = np.float64(self.shape_points[p][0])
-            y2 = np.float64(self.shape_points[p][1])
-
-            # Correction for jagged fin
-            previous_point = int(y1 * 1.0001 / self.span * (points_per_line - 1))
-            current_point = int(y2 * 1.0001 / self.span * (points_per_line - 1))
-            previous_point = max(min(previous_point, points_per_line - 1), 0)
-            current_point = max(min(current_point, points_per_line - 1), 0)
-            if previous_point > current_point:
-                previous_point, current_point = current_point, previous_point
-
-            # Calculate chord length and leading/trailing edge
-            for i in range(previous_point, current_point + 1):
+            x1, y1 = self.shape_points[p - 1]
+            x2, y2 = self.shape_points[p]
+
+            # Compute corresponding points along the fin span (clamp to valid range)
+            prev_idx = int(y1 * 1.0001 / self.span * (points_per_line - 1))
+            curr_idx = int(y2 * 1.0001 / self.span * (points_per_line - 1))
+            prev_idx = np.clip(prev_idx, 0, points_per_line - 1)
+            curr_idx = np.clip(curr_idx, 0, points_per_line - 1)
+
+            if prev_idx > curr_idx:
+                prev_idx, curr_idx = curr_idx, prev_idx
+
+            # Compute intersection of fin edge with each spanwise section
+            for i in range(prev_idx, curr_idx + 1):
                 y = i * self.span / (points_per_line - 1)
-                x = max(
-                    min(
+                if y1 != y2:
+                    x = np.clip(
                         (y - y2) / (y1 - y2) * x1 + (y1 - y) / (y1 - y2) * x2,
+                        min(x1, x2),
                         max(x1, x2),
-                    ),
-                    min(x1, x2),
-                )
-                if x < chord_lead[i]:
-                    chord_lead[i] = x
-                if x > chord_trail[i]:
-                    chord_trail[i] = x
+                    )
+                else:
+                    x = x1  # Handle horizontal segments
+
+                # Update leading and trailing edge positions
+                chord_lead[i] = min(chord_lead[i], x)
+                chord_trail[i] = max(chord_trail[i], x)
 
+                # Update chord length
                 if y1 < y2:
                     chord_length[i] -= x
                 else:
                     chord_length[i] += x
 
-        # Remove infinities and nans for cases where there are equal coordinates
+        # Replace infinities and handle invalid values in chord_lead and chord_trail
         for i in range(points_per_line):
             if (
                 np.isinf(chord_lead[i])
@@ -292,56 +296,58 @@ def evaluate_geometrical_parameters(self):  # pylint: disable=too-many-statement
             if chord_length[i] > chord_trail[i] - chord_lead[i]:
                 chord_length[i] = chord_trail[i] - chord_lead[i]
 
-        # Initialize integration variables
-        radius = self.rocket_radius  # Rocket radius
-        area = 0  # Area of the fin, this will be different from Af if fin is jagged
-        mac_length = 0  # Mean aerodynamic chord length)
-        mac_lead = 0  # Mean aerodynamic chord leading edge x coordinate
-        mac_span = 0  # Mean aerodynamic chord span wise position (Yma)
-        cos_gamma = 0  # Cosine of the sweep angle
+        # Initialize integration variables for various aerodynamic and roll properties
+        radius = self.rocket_radius
+        total_area = 0
+        mac_length = 0  # Mean aerodynamic chord length
+        mac_lead = 0  # Mean aerodynamic chord leading edge
+        mac_span = 0  # Mean aerodynamic chord spanwise position (Yma)
+        cos_gamma_sum = 0  # Sum of cosine of the sweep angle
         roll_geometrical_constant = 0
-        roll_damping_interference_factor_numerator = 0
-        roll_damping_interference_factor_denominator = 0
+        roll_damping_numerator = 0
+        roll_damping_denominator = 0
 
-        # Perform integration
+        # Perform integration over spanwise sections
         dy = self.span / (points_per_line - 1)
         for i in range(points_per_line):
-            length = chord_trail[i] - chord_lead[i]
+            chord = chord_trail[i] - chord_lead[i]
             y = i * dy
 
-            mac_length += length * length
-            mac_span += y * length
-            mac_lead += chord_lead[i] * length
-            area += length
+            # Update integration variables
+            mac_length += chord * chord
+            mac_span += y * chord
+            mac_lead += chord_lead[i] * chord
+            total_area += chord
             roll_geometrical_constant += chord_length[i] * (radius + y) ** 2
-            roll_damping_interference_factor_numerator += (
-                radius**3 * length / (radius + y) ** 2
-            )
-            roll_damping_interference_factor_denominator += (radius + y) * length
+            roll_damping_numerator += radius**3 * chord / (radius + y) ** 2
+            roll_damping_denominator += (radius + y) * chord
 
+            # Update cosine of sweep angle (cos_gamma)
             if i > 0:
                 dx = (chord_trail[i] + chord_lead[i]) / 2 - (
                     chord_trail[i - 1] + chord_lead[i - 1]
                 ) / 2
-                cos_gamma += dy / (dx**2 + dy**2) ** 0.5
+                cos_gamma_sum += dy / np.hypot(dx, dy)
 
+        # Finalize mean aerodynamic chord properties
         mac_length *= dy
         mac_span *= dy
         mac_lead *= dy
-        area *= dy
+        total_area *= dy
         roll_geometrical_constant *= dy
-        roll_damping_interference_factor_numerator *= dy
-        roll_damping_interference_factor_denominator *= dy
-        mac_length /= area
-        mac_span /= area
-        mac_lead /= area
-        cos_gamma /= points_per_line - 1
-
-        # Store values
+        roll_damping_numerator *= dy
+        roll_damping_denominator *= dy
+
+        mac_length /= total_area
+        mac_span /= total_area
+        mac_lead /= total_area
+        cos_gamma = cos_gamma_sum / (points_per_line - 1)
+
+        # Store computed values
         self.Af = Af  # Fin area
-        self.AR = AR  # Aspect Ratio
-        self.gamma_c = np.arccos(cos_gamma)  # Mid chord angle
-        self.Yma = mac_span  # Span wise coord of mean aero chord
+        self.AR = AR  # Aspect ratio
+        self.gamma_c = np.arccos(cos_gamma)  # Sweep angle
+        self.Yma = mac_span  # Mean aerodynamic chord spanwise position
         self.mac_length = mac_length
         self.mac_lead = mac_lead
         self.tau = tau
@@ -349,10 +355,10 @@ def evaluate_geometrical_parameters(self):  # pylint: disable=too-many-statement
         self.lift_interference_factor = lift_interference_factor
         self.roll_forcing_interference_factor = roll_forcing_interference_factor
         self.roll_damping_interference_factor = 1 + (
-            roll_damping_interference_factor_numerator
-            / roll_damping_interference_factor_denominator
+            roll_damping_numerator / roll_damping_denominator
         )
 
+        # Evaluate the shape and finalize geometry
         self.evaluate_shape()
 
     def evaluate_shape(self):

From d2926c76cf9e8d621bc5fc77288a2e06f6faa869 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Wed, 18 Sep 2024 01:35:15 +0200
Subject: [PATCH 07/18] DOC: minor improvement to docs

---
 rocketpy/rocket/aero_surface/fins/free_form_fins.py | 2 ++
 rocketpy/rocket/rocket.py                           | 2 ++
 2 files changed, 4 insertions(+)

diff --git a/rocketpy/rocket/aero_surface/fins/free_form_fins.py b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
index 0b24bf76e..3e6ac2c95 100644
--- a/rocketpy/rocket/aero_surface/fins/free_form_fins.py
+++ b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
@@ -118,6 +118,8 @@ def __init__(
             List of tuples (x, y) containing the coordinates of the fin's
             geometry defining points. The point (0, 0) is the root leading edge.
             Positive x is rearwards, positive y is upwards (span direction).
+            The shape will be interpolated between the points, in the order
+            they are given. The last point connects to the first point.
         rocket_radius : int, float
             Reference radius to calculate lift coefficient, in meters.
         cant_angle : int, float, optional
diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py
index a5c975d87..b422b2800 100644
--- a/rocketpy/rocket/rocket.py
+++ b/rocketpy/rocket/rocket.py
@@ -1346,6 +1346,8 @@ def add_free_form_fins(
             List of tuples (x, y) containing the coordinates of the fin's
             geometry defining points. The point (0, 0) is the root leading edge.
             Positive x is rearwards, positive y is upwards (span direction).
+            The shape will be interpolated between the points, in the order
+            they are given. The last point connects to the first point.
         position : int, float
             Fin set position relative to the rocket's coordinate system.
             By fin set position, understand the point belonging to the root

From d37a3ab0cb9aab00dde507cc3e1f8437a9bf6f49 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Wed, 18 Sep 2024 01:40:42 +0200
Subject: [PATCH 08/18] MNT: lint

---
 rocketpy/rocket/aero_surface/fins/free_form_fins.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/rocketpy/rocket/aero_surface/fins/free_form_fins.py b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
index 3e6ac2c95..a274495ba 100644
--- a/rocketpy/rocket/aero_surface/fins/free_form_fins.py
+++ b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
@@ -201,13 +201,16 @@ def evaluate_center_of_pressure(self):
     def evaluate_geometrical_parameters(self):  # pylint: disable=too-many-statements
         """
         Calculates and saves the fin set's geometrical parameters such as the
-        fin area, aspect ratio, and parameters related to roll movement. This method
-        uses similar calculations to those in OpenRocket for free-form fin shapes.
+        fin area, aspect ratio, and parameters related to roll movement. This
+        method uses the same calculations to those in OpenRocket for free-form
+        fin shapes.
 
         Returns
         -------
         None
         """
+        # pylint: disable=invalid-name
+        # pylint: disable=too-many-locals
         # Calculate the fin area (Af) using the Shoelace theorem (polygon area formula)
         Af = 0
         for i in range(len(self.shape_points) - 1):

From 4ef3d0b8a1afbd2b01b4729e72e8f8399ee62e07 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Thu, 19 Sep 2024 12:25:18 +0200
Subject: [PATCH 09/18] TST: add tests

---
 tests/fixtures/surfaces/surface_fixtures.py | 18 +++++++++
 tests/integration/test_rocket.py            | 45 ++++++++++++++++++++-
 2 files changed, 61 insertions(+), 2 deletions(-)

diff --git a/tests/fixtures/surfaces/surface_fixtures.py b/tests/fixtures/surfaces/surface_fixtures.py
index 10595205f..396206bd7 100644
--- a/tests/fixtures/surfaces/surface_fixtures.py
+++ b/tests/fixtures/surfaces/surface_fixtures.py
@@ -1,6 +1,7 @@
 import pytest
 
 from rocketpy import NoseCone, RailButtons, Tail, TrapezoidalFins
+from rocketpy.rocket.aero_surface.fins.free_form_fins import FreeFormFins
 
 
 @pytest.fixture
@@ -62,6 +63,23 @@ def calisto_trapezoidal_fins():
     )
 
 
+@pytest.fixture
+def calisto_free_form_fins():
+    """The free form fins of the Calisto rocket.
+
+    Returns
+    -------
+    rocketpy.FreeFormFins
+        The free form fins of the Calisto rocket.
+    """
+    return FreeFormFins(
+        n=4,
+        shape_points=[(0, 0), (0.08, 0.1), (0.12, 0.1), (0.12, 0)],
+        rocket_radius=0.0635,
+        name="calisto_free_form_fins",
+    )
+
+
 @pytest.fixture
 def calisto_rail_buttons():
     """The rail buttons of the Calisto rocket.
diff --git a/tests/integration/test_rocket.py b/tests/integration/test_rocket.py
index 76e67f93a..3ec221e7f 100644
--- a/tests/integration/test_rocket.py
+++ b/tests/integration/test_rocket.py
@@ -126,9 +126,50 @@ def test_rocket(mock_show, calisto_robust):  # pylint: disable=unused-argument
 
 @patch("matplotlib.pyplot.show")
 def test_aero_surfaces_infos(  # pylint: disable=unused-argument
-    mock_show, calisto_nose_cone, calisto_tail, calisto_trapezoidal_fins
+    mock_show,
+    calisto_nose_cone,
+    calisto_tail,
+    calisto_trapezoidal_fins,
+    calisto_free_form_fins,
 ):
     assert calisto_nose_cone.all_info() is None
     assert calisto_trapezoidal_fins.all_info() is None
     assert calisto_tail.all_info() is None
-    assert calisto_trapezoidal_fins.draw() is None
+    assert calisto_trapezoidal_fins.all_info() is None
+    assert calisto_free_form_fins.all_info() is None
+
+
+import pytest
+
+
+@pytest.mark.parametrize(
+    "attribute, tolerance",
+    [
+        ("cpz", 1e-3),
+        ("clalpha", 1e-3),
+        ("Af", 1e-3),
+        ("gamma_c", 1e-3),
+        ("Yma", 1e-3),
+        ("tau", 1e-3),
+        ("roll_geometrical_constant", 1e-3),
+        ("lift_interference_factor", 1e-3),
+        ("roll_forcing_interference_factor", 1e-3),
+        ("roll_damping_interference_factor", 1e-2),
+    ],
+)
+def test_calisto_free_form_fins_equivalence(
+    calisto_free_form_fins, calisto_trapezoidal_fins, attribute, tolerance
+):
+    """Test the equivalence of the free form fins with the same geometric
+    characteristics as the trapezoidal fins, comparing cp, cnalpha, and
+    geometrical parameters."""
+
+    # Handle the 'clalpha' method comparison differently as it's a callable
+    if attribute == "clalpha":
+        free_form_value = calisto_free_form_fins.clalpha(0)
+        trapezoidal_value = calisto_trapezoidal_fins.clalpha(0)
+    else:
+        free_form_value = getattr(calisto_free_form_fins, attribute)
+        trapezoidal_value = getattr(calisto_trapezoidal_fins, attribute)
+
+    assert abs(free_form_value - trapezoidal_value) < tolerance

From a5910f70d387a751af6cc75e1b4b1e1c39bd7345 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Fri, 20 Sep 2024 18:27:52 +0200
Subject: [PATCH 10/18] BUG: remove division by zero possibility

---
 rocketpy/rocket/aero_surface/linear_generic_surface.py | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/rocketpy/rocket/aero_surface/linear_generic_surface.py b/rocketpy/rocket/aero_surface/linear_generic_surface.py
index 1ac7aac56..fba2b1db8 100644
--- a/rocketpy/rocket/aero_surface/linear_generic_surface.py
+++ b/rocketpy/rocket/aero_surface/linear_generic_surface.py
@@ -357,11 +357,16 @@ def _compute_from_coefficients(
         # Precompute common values
         dyn_pressure_area = 0.5 * rho * stream_speed**2 * self.reference_area
         dyn_pressure_area_damping = (
-            dyn_pressure_area * self.reference_length / (2 * stream_speed)
+            0.5 * rho * stream_speed * self.reference_area * self.reference_length / 2
         )
         dyn_pressure_area_length = dyn_pressure_area * self.reference_length
         dyn_pressure_area_length_damping = (
-            dyn_pressure_area_length * self.reference_length / (2 * stream_speed)
+            0.5
+            * rho
+            * stream_speed
+            * self.reference_area
+            * self.reference_length**2
+            / 2
         )
 
         # Compute aerodynamic forces

From 49cc5e11f655a7e1c0cd1ba4feecdb615e694f45 Mon Sep 17 00:00:00 2001
From: MateusStano <69485049+MateusStano@users.noreply.github.com>
Date: Fri, 20 Sep 2024 18:42:00 +0200
Subject: [PATCH 11/18] MNT: remove unecessary comment

Co-authored-by: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com>
---
 rocketpy/rocket/rocket.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py
index b422b2800..ec3bfc0ec 100644
--- a/rocketpy/rocket/rocket.py
+++ b/rocketpy/rocket/rocket.py
@@ -1388,7 +1388,6 @@ def add_free_form_fins(
         # Modify radius if not given, use rocket radius, otherwise use given.
         radius = radius if radius is not None else self.radius
 
-        # Create a fin set as an object of TrapezoidalFins class
         fin_set = FreeFormFins(
             n,
             shape_points,

From 62b9ea9a7a2c670da8a6de1f5c292412964666c9 Mon Sep 17 00:00:00 2001
From: MateusStano <69485049+MateusStano@users.noreply.github.com>
Date: Fri, 20 Sep 2024 18:49:26 +0200
Subject: [PATCH 12/18] DOC: fix draw plot title

Co-authored-by: Lucas Prates <57069366+Lucas-Prates@users.noreply.github.com>
---
 rocketpy/plots/aero_surface_plots.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/rocketpy/plots/aero_surface_plots.py b/rocketpy/plots/aero_surface_plots.py
index 85f27c123..80b0b450e 100644
--- a/rocketpy/plots/aero_surface_plots.py
+++ b/rocketpy/plots/aero_surface_plots.py
@@ -434,7 +434,7 @@ def draw(self):
         # Plot settings
         ax.set_xlabel("Root chord (m)")
         ax.set_ylabel("Span (m)")
-        ax.set_title("Trapezoidal Fin Cross Section")
+        ax.set_title("Free Form Fin Cross Section")
         ax.legend(bbox_to_anchor=(1.05, 1.0), loc="upper left")
 
         plt.tight_layout()

From 879bba88537405f1a6f8e24f1766692891a956aa Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Fri, 20 Sep 2024 19:01:02 +0200
Subject: [PATCH 13/18] DOC: remove unused attributes

---
 rocketpy/rocket/aero_surface/fins/free_form_fins.py | 12 ------------
 1 file changed, 12 deletions(-)

diff --git a/rocketpy/rocket/aero_surface/fins/free_form_fins.py b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
index a274495ba..badbf3b37 100644
--- a/rocketpy/rocket/aero_surface/fins/free_form_fins.py
+++ b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
@@ -37,26 +37,14 @@ class FreeFormFins(Fins):
         Second is the unit of the curve (radians or degrees).
     FreeFormFins.cant_angle : float
         Fins cant angle with respect to the rocket centerline, in degrees.
-    FreeFormFins.changing_attribute_dict : dict
-        Dictionary that stores the name and the values of the attributes that
-        may be changed during a simulation. Useful for control systems.
     FreeFormFins.cant_angle_rad : float
         Fins cant angle with respect to the rocket centerline, in radians.
     FreeFormFins.root_chord : float
         Fin root chord in meters.
-    FreeFormFins.tip_chord : float
-        Fin tip chord in meters.
     FreeFormFins.span : float
         Fin span in meters.
     FreeFormFins.name : string
         Name of fin set.
-    FreeFormFins.sweep_length : float
-        Fins sweep length in meters. By sweep length, understand the axial
-        distance between the fin root leading edge and the fin tip leading edge
-        measured parallel to the rocket centerline.
-    FreeFormFins.sweep_angle : float
-        Fins sweep angle with respect to the rocket centerline. Must
-        be given in degrees.
     FreeFormFins.d : float
         Reference diameter of the rocket, in meters.
     FreeFormFins.ref_area : float

From bff7ecb61960f3f5c637e9f3987f2aafa924dc96 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Fri, 20 Sep 2024 19:01:36 +0200
Subject: [PATCH 14/18] ENH: remove unecessary multiplication from opr

---
 rocketpy/rocket/aero_surface/fins/free_form_fins.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/rocketpy/rocket/aero_surface/fins/free_form_fins.py b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
index badbf3b37..7fa409309 100644
--- a/rocketpy/rocket/aero_surface/fins/free_form_fins.py
+++ b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
@@ -244,8 +244,8 @@ def evaluate_geometrical_parameters(self):  # pylint: disable=too-many-statement
             x2, y2 = self.shape_points[p]
 
             # Compute corresponding points along the fin span (clamp to valid range)
-            prev_idx = int(y1 * 1.0001 / self.span * (points_per_line - 1))
-            curr_idx = int(y2 * 1.0001 / self.span * (points_per_line - 1))
+            prev_idx = int(y1 / self.span * (points_per_line - 1))
+            curr_idx = int(y2 / self.span * (points_per_line - 1))
             prev_idx = np.clip(prev_idx, 0, points_per_line - 1)
             curr_idx = np.clip(curr_idx, 0, points_per_line - 1)
 

From ca460b812e8d88b677c32aaf57db283abce9175c Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Fri, 20 Sep 2024 19:01:49 +0200
Subject: [PATCH 15/18] ENH: add last line for fin plot

---
 rocketpy/plots/aero_surface_plots.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/rocketpy/plots/aero_surface_plots.py b/rocketpy/plots/aero_surface_plots.py
index 85f27c123..b3f02c957 100644
--- a/rocketpy/plots/aero_surface_plots.py
+++ b/rocketpy/plots/aero_surface_plots.py
@@ -426,6 +426,12 @@ def draw(self):
             self.aero_surface.shape_vec[1],
             color="#A60628",
         )
+        # line from the last point to the first point
+        ax.plot(
+            [self.aero_surface.shape_vec[0][-1], self.aero_surface.shape_vec[0][0]],
+            [self.aero_surface.shape_vec[1][-1], self.aero_surface.shape_vec[1][0]],
+            color="#A60628",
+        )
 
         ax.add_line(yma_line)
         ax.scatter(*cp_point, label="Center of Pressure", color="red", s=100, zorder=10)

From 378a4f9c8887da99d95f12c8c0a9486b843f7899 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Sat, 21 Sep 2024 11:38:52 +0200
Subject: [PATCH 16/18] MNT: flake

---
 tests/integration/test_rocket.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/tests/integration/test_rocket.py b/tests/integration/test_rocket.py
index 3ec221e7f..15b213189 100644
--- a/tests/integration/test_rocket.py
+++ b/tests/integration/test_rocket.py
@@ -139,9 +139,6 @@ def test_aero_surfaces_infos(  # pylint: disable=unused-argument
     assert calisto_free_form_fins.all_info() is None
 
 
-import pytest
-
-
 @pytest.mark.parametrize(
     "attribute, tolerance",
     [

From 502379edeacb0075775227ee2e9e516bd09840e5 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Sat, 21 Sep 2024 11:42:57 +0200
Subject: [PATCH 17/18] MNT: isort

---
 tests/integration/test_rocket.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/integration/test_rocket.py b/tests/integration/test_rocket.py
index 15b213189..e51fe925e 100644
--- a/tests/integration/test_rocket.py
+++ b/tests/integration/test_rocket.py
@@ -1,6 +1,7 @@
 from unittest.mock import patch
 
 import numpy as np
+import pytest
 
 
 @patch("matplotlib.pyplot.show")

From 49abe78fdbc80cf1a15937f366a2332a9eae4d89 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Sat, 21 Sep 2024 12:29:42 +0200
Subject: [PATCH 18/18] DOC: improve shape point docs

---
 rocketpy/rocket/aero_surface/fins/free_form_fins.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/rocketpy/rocket/aero_surface/fins/free_form_fins.py b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
index 7fa409309..1db2ef18d 100644
--- a/rocketpy/rocket/aero_surface/fins/free_form_fins.py
+++ b/rocketpy/rocket/aero_surface/fins/free_form_fins.py
@@ -107,7 +107,8 @@ def __init__(
             geometry defining points. The point (0, 0) is the root leading edge.
             Positive x is rearwards, positive y is upwards (span direction).
             The shape will be interpolated between the points, in the order
-            they are given. The last point connects to the first point.
+            they are given. The last point connects to the first point, and
+            represents the trailing edge.
         rocket_radius : int, float
             Reference radius to calculate lift coefficient, in meters.
         cant_angle : int, float, optional