Repository: petrobras/ross (language: python)

Issue:
Problems with 6 dof model
Hi @raphaeltimbo,

As I mentioned in Issue #1034, I decided to conduct the same analyses performed in the ROSS tutorial part 2 [🔗](https://ross.readthedocs.io/en/latest/user_guide/tutorial_part_2.html) for the 6 dof model. The objective was to identify where errors could arise within the code. As some problems were indeed found, I opened this new issue to report them.


Discussion / hints:
Upon executing part 2 of the tutorial, I encountered the same error when invoking `run_static()` and `run_ucs()`:

```python
in Rotor.__init__(self, shaft_elements, disk_elements, bearing_elements, point_mass_elements, min_w, max_w, rated_w, tag)
    209 if len(tags_list) != len(set(tags_list)):
    210     raise ValueError("Tags should be unique.")
--> 212 self.number_dof = self._check_number_dof()
    214 ####################################################
    215 # Rotor summary
    216 ####################################################
    217 columns = [
    218     "type",
    219     "n",
   (...)
    237     "tag",
    238 ]
...
    562         "The number of degrees o freedom of all elements must be the same! There are BEARING elements with discrepant DoFs."
    563     )
    565 return int(number_dof)

Exception: The number of degrees o freedom of all elements must be the same! There are BEARING elements with discrepant DoFs.
```
which can be fixed by replacing the call of `BearingElement()`class with `elm.__class__()`. After fixing that, I compared the results obtained from the 4 dof model with the 6 dof model, and observed the following disparities:

- **Running modal analysis:**
![compare_1](https://github.com/petrobras/ross/assets/82293939/5f4fd8dd-2e16-4f8c-8a92-f4b97e3aaec1)
![compare_2](https://github.com/petrobras/ross/assets/82293939/e8479de8-e2f2-4ce7-98e9-193eec8ff4f8)

- **Plotting Campbell diagram:**
![compare_3](https://github.com/petrobras/ross/assets/82293939/f053f70b-8ba6-4fcd-aed9-7a8d75a4d2a7)

- **Plotting the deflected shape for the unbalance response:**
![compare_4](https://github.com/petrobras/ross/assets/82293939/886ed5a5-761b-495a-95d9-02009c151203)

- **Plotting the deflected shape for the unbalance response:**
![compare_5](https://github.com/petrobras/ross/assets/82293939/6418f53c-9e0d-496d-afc7-82711f7cfc1e)
So, I started debugging the result classes in `results.py` to check if there was any coding that fixed the dof number to 4. Therefore, I needed to make some changes mainly in `Shape` and `ForcedResponseResults`. After that, I got the following results for the 6 dof model:
![newplot(1)](https://github.com/petrobras/ross/assets/82293939/c938f152-91d2-47f3-9eb1-d7d4b45d5e7f)
![newplot](https://github.com/petrobras/ross/assets/82293939/e54a6cca-cfb7-48cf-94c4-b08121457f73)
![newplot(2)](https://github.com/petrobras/ross/assets/82293939/bc9c2095-06c3-42f6-a3f0-7a4982e58538)
![newplot(3)](https://github.com/petrobras/ross/assets/82293939/87ef48d3-5ddd-488e-a2d9-887291b83442)
![newplot(4)](https://github.com/petrobras/ross/assets/82293939/23d902cc-4987-4148-bfd5-a141a48cbbaf)

Fix patch:
diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py
index 2f709389c..e9e162dad 100644
--- a/ross/rotor_assembly.py
+++ b/ross/rotor_assembly.py
@@ -1014,13 +1014,13 @@ def Ksdt(self):
         Examples
         --------
         >>> rotor = rotor_example_6dof()
-        >>> np.round(rotor.Ksdt()[:6, :6]*1e6)
-        array([[     0., -23002.,      0.,   -479.,      0.,      0.],
-               [     0.,      0.,      0.,      0.,      0.,      0.],
-               [     0.,      0.,      0.,      0.,      0.,      0.],
-               [     0.,      0.,      0.,      0.,      0.,      0.],
-               [     0.,    479.,      0.,    160.,      0.,      0.],
-               [     0.,      0.,      0.,      0.,      0.,      0.]])
+        >>> np.round(rotor.Ksdt()[:6, :6] * 1e3, 2)
+        array([[  0.  , -23.  ,   0.  ,   0.48,   0.  ,   0.  ],
+               [  0.  ,   0.  ,   0.  ,   0.  ,   0.  ,   0.  ],
+               [  0.  ,   0.  ,   0.  ,   0.  ,   0.  ,   0.  ],
+               [  0.  ,   0.  ,   0.  ,   0.  ,   0.  ,   0.  ],
+               [  0.  ,  -0.48,   0.  ,   0.16,   0.  ,   0.  ],
+               [  0.  ,   0.  ,   0.  ,   0.  ,   0.  ,   0.  ]])
         """
         Ksdt0 = np.zeros((self.ndof, self.ndof))
 
diff --git a/ross/shaft_element.py b/ross/shaft_element.py
index 2e14ea4b4..2de87fd9b 100644
--- a/ross/shaft_element.py
+++ b/ross/shaft_element.py
@@ -72,6 +72,12 @@ class ShaftElement(Element):
     shear_method_calc : str, optional
         Determines which shear calculation method the user will adopt
         Default is 'cowper'
+    alpha : float, optional
+        Mass proportional damping factor.
+        Default is zero.
+    beta : float, optional
+        Stiffness proportional damping factor.
+        Default is zero.
     tag : str, optional
         Element tag.
         Default is None.
@@ -149,6 +155,8 @@ def __init__(
         gyroscopic=True,
         shear_method_calc="cowper",
         tag=None,
+        alpha=0,
+        beta=0,
     ):
         if idr is None:
             idr = idl
@@ -186,8 +194,8 @@ def __init__(
         self.odr = float(odr)
         self.color = self.material.color
 
-        self.alpha = 0.0
-        self.beta = 0.0
+        self.alpha = float(alpha)
+        self.beta = float(beta)
 
         # A_l = cross section area from the left side of the element
         # A_r = cross section area from the right side of the element
@@ -199,6 +207,9 @@ def __init__(
         # Second moment of area of the cross section from the left side
         # of the element
         Ie_l = np.pi * (odl**4 - idl**4) / 64
+        Ie_r = np.pi * (odr**4 - idr**4) / 64
+        self.Ie_l = Ie_l
+        self.Ie_r = Ie_r
 
         outer = self.odl**2 + self.odl * self.odr + self.odr**2
         inner = self.idl**2 + self.idl * self.idr + self.idr**2
@@ -235,7 +246,6 @@ def __init__(
         # the middle of the element - about the neutral plane
         Ie = Ie_l * (1 + a2 * 0.5 + b2 * 0.5**2 + gama * 0.5**3 + delta * 0.5**4)
         self.Ie = Ie
-        self.Ie_l = Ie_l
 
         phi = 0
 
@@ -796,12 +806,13 @@ def C(self):
         ...                         rotary_inertia=True,
         ...                         shear_effects=True)
         >>> Timoshenko_Element.C()[:4, :4]
-        array([[0., 0., 0., 0.],
-               [0., 0., 0., 0.],
-               [0., 0., 0., 0.],
-               [0., 0., 0., 0.]])
+        array([[ 0.,  0.,  0.,  0.],
+               [ 0.,  0., -0.,  0.],
+               [ 0., -0.,  0.,  0.],
+               [ 0.,  0.,  0.,  0.]])
         """
-        C = np.zeros((8, 8))
+        # proportional damping matrix
+        C = self.alpha * self.M() + self.beta * self.K()
 
         return C
 
@@ -1074,6 +1085,8 @@ def from_table(cls, file, sheet_type="Simple", sheet_name=0):
                         rotary_inertia=parameters["rotary_inertia"][i],
                         gyroscopic=parameters["gyroscopic"][i],
                         shear_method_calc=parameters["shear_method_calc"][i],
+                        alpha=parameters["alpha"][i] if "alpha" in parameters else 0,
+                        beta=parameters["beta"][i] if "beta" in parameters else 0,
                     )
                 )
         elif sheet_type == "Simple":
@@ -1093,6 +1106,8 @@ def from_table(cls, file, sheet_type="Simple", sheet_name=0):
                         rotary_inertia=parameters["rotary_inertia"][i],
                         gyroscopic=parameters["gyroscopic"][i],
                         shear_method_calc=parameters["shear_method_calc"][i],
+                        alpha=parameters["alpha"][i] if "alpha" in parameters else 0,
+                        beta=parameters["beta"][i] if "beta" in parameters else 0,
                     )
                 )
         return list_of_shafts
@@ -1111,6 +1126,8 @@ def section(
         shear_effects=True,
         rotary_inertia=True,
         gyroscopic=True,
+        alpha=0,
+        beta=0,
     ):
         """Shaft section constructor.
 
@@ -1147,6 +1164,12 @@ def section(
         gyroscopic : bool
             Determine if gyroscopic effects are taken into account.
             Default is False.
+        alpha : float, optional
+            Mass proportional damping factor.
+            Default is zero.
+        beta : float, optional
+            Stiffness proportional damping factor.
+            Default is zero.
 
         Returns
         -------
@@ -1185,6 +1208,8 @@ def section(
                 shear_effects=shear_effects,
                 rotary_inertia=rotary_inertia,
                 gyroscopic=gyroscopic,
+                alpha=alpha,
+                beta=beta,
             )
             for i in range(ne)
         ]
@@ -1193,86 +1218,101 @@ def section(
 
 
 class ShaftElement6DoF(ShaftElement):
-    r"""A 6 Degrees of Freedom shaft element.
+    r"""A shaft element with 6 degrees of freedom.
 
-    This class will create a shaft element that takes into
-    account shear stress, rotary inertia and gyroscopic effects.
+    This class will create a shaft element that may take into
+    account shear, rotary inertia an gyroscopic effects. The object can be
+    cylindrical or conical and the formulation is based on :cite:`genta1988conical`
     The matrices will be defined considering the following local
     coordinate vector:
 
     .. math::
 
-        [u_0, v_0, w_0, \theta_0, \psi_0, \phi_0, u_1, v_1, w_1, \theta_1, \psi_1, \phi_1]^T
-    Being the following their ordering for an element:
-
-    :math:`x_0,u_0`  - horizontal translation;
-
-    :math:`y_0,v_0`  - vertical translation;
-
-    :math:`z_0,w_0`  - axial translation;
-
-    :math:`\theta_0` - rotation around horizontal, bending on the yz plane;
-
-    :math:`\psi_0`   - rotation around vertical, bending on the xz plane;
-
-    :math:`\phi_0`   - torsion around axial, z direction.
+        [x_0, y_0, z_0, \alpha_0, \beta_0, \theta_0, x_1, y_1, z_1, \alpha_1, \beta_1, \theta_1]^T
+    Where :math:`\x_0` and :math:`\x_1` correspond to horizontal translation,
+    :math:`\y_0` and :math:`\y_1` correspond to vertical translation,
+    :math:`\z_0` and :math:`\z_1` correspond to axial translation,
+    :math:`\alpha_0` and :math:`\alpha_1` correspond to bending on the yz plane,
+    :math:`\beta_0` and :math:`\beta_1` correspond to bending on the xz plane
+    and :math:`\theta_0` and :math:`\theta_1` correspond to torsion around z axis.
 
     Parameters
     ----------
     L : float, pint.Quantity
-        Element length.
+        Element length (m).
     idl : float, pint.Quantity
-        Inner diameter of the element at the left node (m).
+        Inner diameter of the element at the left position (m).
     odl : float, pint.Quantity
-        Outer diameter of the element at the left node (m).
+        Outer diameter of the element at the left position (m).
     idr : float, pint.Quantity, optional
-        Inner diameter of the element at the right node (m).
-        Default is equal to idl value for cylindrical element.
+        Inner diameter of the element at the right position (m).
+        Default is equal to idl value (cylindrical element).
     odr : float, pint.Quantity, optional
-        Outer diameter of the element at the right node (m).
-        Default is equal to odl value for cylindrical element.
+        Outer diameter of the element at the right position (m).
+        Default is equal to odl value (cylindrical element).
     material : ross.Material
         Shaft material.
-    alpha : float, optional
-        Proportional damping coefficient, associated to the element Mass matrix
-    beta : float, optional
-        Proportional damping coefficient, associated to the element Stiffness matrix
     n : int, optional
-        Element number, coincident it's first node.
+        Element number (coincident with it's first node).
         If not given, it will be set when the rotor is assembled
         according to the element's position in the list supplied to
         the rotor constructor.
     axial_force : float, optional
         Axial force (N).
-        Default is zero.
     torque : float, optional
-        Torque moment (N*m).
-        Default is zero.
+        Torque (N*m).
     shear_effects : bool, optional
-        Determine if shear effects are taken into account;
+        Determine if shear effects are taken into account.
         Default is True.
     rotary_inertia : bool, optional
-        Determine if rotary_inertia effects are taken into account;
+        Determine if rotary_inertia effects are taken into account.
         Default is True.
     gyroscopic : bool, optional
-        Determine if gyroscopic effects are taken into account;
+        Determine if gyroscopic effects are taken into account.
         Default is True.
+    shear_method_calc : str, optional
+        Determines which shear calculation method the user will adopt
+        Default is 'cowper'
+    alpha : float, optional
+        Mass proportional damping factor.
+        Default is zero.
+    beta : float, optional
+        Stiffness proportional damping factor.
+        Default is zero.
     tag : str, optional
-        Element tag;
+        Element tag.
         Default is None.
 
     Returns
     -------
-    shaft_element : rs.ShaftElement6DoF
-        A 6 degrees of freedom shaft element, with available gyroscopic, shear and rotary inertia effects.
+    shaft_element : ross.ShaftElement6DoF
+        A 6 degrees of freedom shaft element object.
 
     Attributes
     ----------
     Poisson : float
         Poisson coefficient for the element.
+    A : float
+        Element section area at half length (m**2).
+    A_l : float
+        Element section area at left end (m**2).
+    A_r : float
+        Element section area at right end (m**2).
+    beam_cg : float
+        Element center of gravity local position (m).
+    axial_cg_pos : float
+        Element center of gravity global position (m).
+        This should be used only after the rotor is built.
+        Default is None.
+    Ie : float
+        Ie is the second moment of area of the cross section about
+        the neutral plane (m**4).
+    phi : float
+        Constant that is used according to :cite:`friswell2010dynamics` to
+        consider rotary inertia and shear effects. If these are not considered
+        :math:`\phi=0`.
     kappa : float
-        Shear coefficient for the element, determined from :cite:`Hutchingson2001`
-        formulation.
+        Shear coefficient for the element.
 
     References
     ----------
@@ -1282,315 +1322,69 @@ class ShaftElement6DoF(ShaftElement):
     Examples
     --------
     >>> from ross.materials import steel
-    >>> shaft1 = ShaftElement6DoF(L=0.5, idl=0.0, odl=0.01, idr=0.0, odr=0.01,
-    ...                           material=steel, n=0, axial_force=10, torque=30)
-    >>> shaft2 = ShaftElement6DoF(L=0.5, idl=0.05, odl=0.1, idr=0.05, odr=0.15,
-    ...                           alpha=0.01, beta=100, material=steel,
-    ...                           rotary_inertia=False, shear_effects=False)
-    >>> shaft2.kappa
-    0.7099387976608923
+    >>> # Euler-Bernoulli conical element
+    >>> Euler_Bernoulli_Element = ShaftElement6DoF(
+    ...                         material=steel, L=0.5, idl=0.05, odl=0.1,
+    ...                         idr=0.05, odr=0.15,
+    ...                         rotary_inertia=False,
+    ...                         shear_effects=False)
+    >>> Euler_Bernoulli_Element.phi
+    0
+    >>> # Timoshenko cylindrical element. In this case idr and odr are omitted.
+    >>> Timoshenko_Element = ShaftElement6DoF(
+    ...                         material=steel, L=0.5, idl=0.05, odl=0.1,
+    ...                         rotary_inertia=True,
+    ...                         shear_effects=True)
+    >>> Timoshenko_Element.phi
+    0.1571268472906404
     """
 
-    @check_units
-    def __init__(
-        self,
-        L,
-        idl,
-        odl,
-        idr=None,
-        odr=None,
-        material=None,
-        n=None,
-        axial_force=0,
-        torque=0,
-        shear_effects=True,
-        rotary_inertia=True,
-        gyroscopic=True,
-        alpha=0,
-        beta=0,
-        tag=None,
-    ):
-        if idr is None:
-            idr = idl
-        if odr is None:
-            odr = odl
-
-        if material is None:
-            raise AttributeError("Material is not defined.")
-
-        if type(material) is str:
-            self.material = Material.load_material(material)
-        else:
-            self.material = material
-
-        self.shear_effects = shear_effects
-        self.rotary_inertia = rotary_inertia
-        self.gyroscopic = gyroscopic
-        self.axial_force = axial_force
-        self.torque = torque
-        self._n = n
-        self.n_l = n
-        self.n_r = None
-        if n is not None:
-            self.n_r = n + 1
-
-        self.tag = tag
-
-        self.L = float(L)
-        self.o_d = (float(odl) + float(odr)) / 2
-        self.i_d = (float(idl) + float(idr)) / 2
-        self.idl = float(idl)
-        self.odl = float(odl)
-        self.idr = float(idr)
-        self.odr = float(odr)
-        self.color = self.material.color
-
-        # A_l = cross section area from the left side of the element
-        # A_r = cross section area from the right side of the element
-        A_l = np.pi * (odl**2 - idl**2) / 4
-        A_r = np.pi * (odr**2 - idr**2) / 4
-        self.A_l = A_l
-        self.A_r = A_r
-
-        # Second moment of area of the cross section from the left side
-        # of the element
-        Ie_l = np.pi * (odl**4 - idl**4) / 64
-
-        outer = self.odl**2 + self.odl * self.odr + self.odr**2
-        inner = self.idl**2 + self.idl * self.idr + self.idr**2
-        self.volume = np.pi * (self.L / 12) * (outer - inner)
-        self.m = self.material.rho * self.volume
-
-        roj = odl / 2
-        rij = idl / 2
-        rok = odr / 2
-        rik = idr / 2
-
-        # geometrical coefficients
-        delta_ro = rok - roj
-        delta_ri = rik - rij
-        a1 = 2 * np.pi * (roj * delta_ro - rij * delta_ri) / A_l
-        a2 = np.pi * (roj**3 * delta_ro - rij**3 * delta_ri) / Ie_l
-        b1 = np.pi * (delta_ro**2 - delta_ri**2) / A_l
-        b2 = 3 * np.pi * (roj**2 * delta_ro**2 - rij**2 * delta_ri**2) / (2 * Ie_l)
-        gama = np.pi * (roj * delta_ro**3 - rij * delta_ri**3) / Ie_l
-        delta = np.pi * (delta_ro**4 - delta_ri**4) / (4 * Ie_l)
-
-        self.a1 = a1
-        self.a2 = a2
-        self.b1 = b1
-        self.b2 = b2
-        self.gama = gama
-        self.delta = delta
-
-        # the area is calculated from the cross section located in the middle
-        # of the element
-        self.A = A_l * (1 + a1 * 0.5 + b1 * 0.5**2)
-
-        # Ie is the second moment of area of the cross section - located in
-        # the middle of the element - about the neutral plane
-        Ie = Ie_l * (1 + a2 * 0.5 + b2 * 0.5**2 + gama * 0.5**3 + delta * 0.5**4)
-        self.Ie = Ie
-        self.Ie_l = Ie_l
-
-        # geometric center
-        c1 = roj**2 + 2 * roj * rok + 3 * rok**2 - rij**2 - 2 * rij * rik - 3 * rik**2
-        c2 = (roj**2 + roj * rok + rok**2) - (rij**2 + rij * rik + rik**2)
-        self.beam_cg = L * c1 / (4 * c2)
-        self.axial_cg_pos = None
-
-        # Slenderness ratio of beam elements (G*A*L**2) / (E*I)
-        sld = (self.material.G_s * self.A * self.L**2) / (self.material.E * Ie)
-        self.slenderness_ratio = sld
-
-        # Moment of inertia
-        # fmt: off
-        self.Im = (
-            (np.pi * L * (self.m / self.volume) / 10) *
-            ((roj ** 4 + roj ** 3 * rok + roj ** 2 * rok ** 2 + roj * rok ** 3 + rok ** 4) -
-             (rij ** 4 + rij ** 3 * rik + rij ** 2 * rik ** 2 + rij * rik ** 3 + rik ** 4))
-        )
-        # fmt: on
-
-        self.alpha = float(alpha)
-        self.beta = float(beta)
-
-        # Timoshenko kappa factor determination, based on the diameters relation
-        if self.__is_circular():
-            kappa = (6 * (1 + self.material.Poisson) ** 2) / (
-                7 + 12 * self.material.Poisson + 4 * self.material.Poisson**2
-            )
-        elif self.__is_thickwall():
-            a = (self.idl + self.idr) / 2
-            b = (self.odl + self.odr) / 2
-            v = self.material.Poisson
-            kappa = (6 * (a**2 + b**2) ** 2 * (1 + v) ** 2) / (
-                7 * a**4
-                + 34 * a**2 * b**2
-                + 7 * b**4
-                + v * (12 * a**4 + 48 * a**2 * b**2 + 12 * b**4)
-                + v**2 * (4 * a**4 + 16 * a**2 * b**2 + 4 * b**4)
-            )
-        else:
-            kappa = (1 + self.material.Poisson) / (2 + self.material.Poisson)
-
-        self.kappa = kappa
-
-        self.dof_global_index = None
-
-    def __is_circular(self):
-        return self.idl == 0 and self.idr == 0
-
-    def __is_thickwall(self):
-        p = (((self.odl + self.odr) / 2) - ((self.idl + self.idr) / 2)) / (
-            (self.odl + self.odr) / 2
-        )
-        return p >= 0.2
+    def dof_mapping(self):
+        """Degrees of freedom mapping.
 
-    def __repr__(self):
-        """Return a string representation of a shaft element.
+        Returns a dictionary with a mapping between degree of freedom and its index.
 
         Returns
         -------
-        A string representation of a 6 DoF shaft object.
+        dof_mapping : dict
+            A dictionary containing the degrees of freedom and their indexes
+            with the following ordering:
+            x_0 - horizontal translation
+            y_0 - vertical translation
+            z_0 - axial translation
+            alpha_0 - rotation around horizontal
+            beta_0  - rotation around vertical
+            theta_0 - torsion around axial
 
         Examples
         --------
-        >>> from ross.materials import steel
-        >>> shaft1 = ShaftElement6DoF(
-        ...        L=0.25, idl=0, idr=0, odl=0.05, odr=0.08,
-        ...        material=steel, rotary_inertia=True, shear_effects=True
-        ... )
-        >>> shaft1 # doctest: +ELLIPSIS
-        ShaftElement6DoF(L=0.25, idl=0.0...
+        >>> sh = ShaftElement6DoF(L=0.5, idl=0.05, odl=0.1, material=steel,
+        ...                       rotary_inertia=True, shear_effects=True)
+        >>> sh.dof_mapping()["x_0"]
+        0
         """
-        return (
-            f"{self.__class__.__name__}"
-            f"(L={self.L:{0}.{5}}, idl={self.idl:{0}.{5}}, "
-            f"idr={self.idr:{0}.{5}}, odl={self.odl:{0}.{5}},  "
-            f"odr={self.odr:{0}.{5}}, material={self.material.name!r}, "
-            f"alpha={self.alpha:{0}.{5}}, beta={self.beta:{0}.{5}},  "
-            f"n={self.n})"
+        return dict(
+            x_0=0,
+            y_0=1,
+            z_0=2,
+            alpha_0=3,
+            beta_0=4,
+            theta_0=5,
+            x_1=6,
+            y_1=7,
+            z_1=8,
+            alpha_1=9,
+            beta_1=10,
+            theta_1=11,
         )
 
-    def __str__(self):
-        """Convert object into string.
+    def M(self):
+        """Mass matrix for an instance of a 6 DoF shaft element.
 
         Returns
         -------
-        The object's parameters translated to strings
-
-        Example
-        -------
-        >>> print(ShaftElement6DoF(L=0.25, idl=0, odl=0.05, odr=0.08, material=steel))
-        Element Number:             None
-        Element Lenght   (m):       0.25
-        Left Int. Diam.  (m):        0.0
-        Left Out. Diam.  (m):       0.05
-        Right Int. Diam. (m):        0.0
-        Right Out. Diam. (m):       0.08
-        Alpha damp. factor:          0.0
-        Beta damp. factor:           0.0
-        -----------------------------------
-        Steel
-        -----------------------------------
-        Density         (kg/m**3): 7810.0
-        Young`s modulus (N/m**2):  2.11e+11
-        Shear modulus   (N/m**2):  8.12e+10
-        Poisson coefficient     :  0.29926108
-        """
-        return (
-            f"Element Number:             {self.n}"
-            f"\nElement Lenght   (m): {self.L:{10}.{5}}"
-            f"\nLeft Int. Diam.  (m): {self.idl:{10}.{5}}"
-            f"\nLeft Out. Diam.  (m): {self.odl:{10}.{5}}"
-            f"\nRight Int. Diam. (m): {self.idr:{10}.{5}}"
-            f"\nRight Out. Diam. (m): {self.odr:{10}.{5}}"
-            f"\nAlpha damp. factor:   {self.alpha:{10}.{5}}"
-            f"\nBeta damp. factor:    {self.beta:{10}.{5}}"
-            f'\n{35*"-"}'
-            f"\n{self.material}"
-        )
-
-    def save(self, file):
-        signature = inspect.signature(self.__init__)
-        args_list = list(signature.parameters)
-        args = {arg: getattr(self, arg) for arg in args_list}
-
-        # add material characteristics so that the shaft element can be reconstructed
-        # even if the material is not in the available_materials file.
-        args["material"] = {
-            "name": self.material.name,
-            "rho": self.material.rho,
-            "E": self.material.E,
-            "G_s": self.material.G_s,
-            "color": self.material.color,
-        }
-
-        try:
-            data = toml.load(file)
-        except FileNotFoundError:
-            data = {}
-
-        data[f"{self.__class__.__name__}_{self.tag}"] = args
-        with open(file, "w") as f:
-            toml.dump(data, f)
-
-    @classmethod
-    def read_toml_data(cls, data):
-        data["material"] = Material(**data["material"])
-        return cls(**data)
-
-    def dof_mapping(self):
-        """Degrees of freedom mapping.
-
-        Returns a dictionary with a mapping between degree of freedom and its index.
-
-        Returns
-        -------
-        dof_mapping : dict
-            A dictionary containing the degrees of freedom and their indexes.
-
-        Examples
-        --------
-        The numbering of the degrees of freedom for each node.
-
-        Being the following their ordering for a node:
-
-        x_0 - horizontal translation
-        y_0 - vertical translation
-        z_0 - axial translation
-        alpha_0 - rotation around horizontal
-        beta_0  - rotation around vertical
-        theta_0 - torsion around axial
-
-        >>> sh = ShaftElement6DoF(L=0.5, idl=0.05, odl=0.1, material=steel,
-        ...                       rotary_inertia=True, shear_effects=True)
-        >>> sh.dof_mapping()["x_0"]
-        0
-        """
-        return dict(
-            x_0=0,
-            y_0=1,
-            z_0=2,
-            alpha_0=3,
-            beta_0=4,
-            theta_0=5,
-            x_1=6,
-            y_1=7,
-            z_1=8,
-            alpha_1=9,
-            beta_1=10,
-            theta_1=11,
-        )
-
-    def M(self):
-        """Mass matrix for an instance of a 6 DoF shaft element.
-
-        Returns
-        -------
-        M : np.ndarray
-            Mass matrix for the 6 DoF shaft element.
+        M : np.ndarray
+            Mass matrix for the 6 DoF shaft element.
 
         Examples
         --------
@@ -1598,59 +1392,143 @@ def M(self):
         >>> Timoshenko_Element.M().shape
         (12, 12)
         """
-        # temporary material and geometrical constants
+        phi = self.phi
         L = self.L
-        tempS = np.pi * (
-            ((self.odr / 2) ** 2 + (self.odl / 2) ** 2) / 2
-            - ((self.idr / 2) ** 2 + (self.idl / 2) ** 2) / 2
+        a1 = self.a1
+        a2 = self.a2
+        b1 = self.b1
+        b2 = self.b2
+        delta = self.delta
+        gama = self.gama
+
+        m1 = (
+            (468 + 882 * phi + 420 * phi**2)
+            + a1 * (108 + 210 * phi + 105 * phi**2)
+            + b1 * (38 + 78 * phi + 42 * phi**2)
         )
-        tempI = (
-            np.pi
-            / 4
-            * (
-                ((self.odr / 2) ** 4 + (self.odl / 2) ** 4) / 2
-                - ((self.idr / 2) ** 4 + (self.idl / 2) ** 4) / 2
-            )
+        m2 = (
+            (66 + 115.5 * phi + 52.5 * phi**2)
+            + a1 * (21 + 40.5 * phi + 21 * phi**2)
+            + b1 * (8.5 + 18 * phi + 10.5 * phi**2)
+        ) * L
+        m3 = (
+            (162 + 378 * phi + 210 * phi**2)
+            + a1 * (81 + 189 * phi + 105 * phi**2)
+            + b1 * (46 + 111 * phi + 63 * phi**2)
+        )
+        m4 = (
+            (39 + 94.5 * phi + 52.5 * phi**2)
+            + a1 * (18 + 40.5 * phi + 21 * phi**2)
+            + b1 * (9.5 + 21 * phi + 10.5 * phi**2)
+        ) * L
+        m5 = (
+            (12 + 21 * phi + 10.5 * phi**2)
+            + a1 * (4.5 + 9 * phi + 5.25 * phi**2)
+            + b1 * (2 + 4.5 * phi + 3 * phi**2)
+        ) * L**2
+        m6 = (
+            (39 + 94.5 * phi + 52.5 * phi**2)
+            + a1 * (21 + 54 * phi + 31.5 * phi**2)
+            + b1 * (12.5 + 34.5 * phi + 21 * phi**2)
+        ) * L
+        m7 = (
+            (9 + 21 * phi + 10.5 * phi**2)
+            + a1 * (4.5 + 10.5 * phi + 5.25 * phi**2)
+            + b1 * (2.5 + 6 * phi + 3 * phi**2)
+        ) * L**2
+        m8 = (
+            (468 + 882 * phi + 420 * phi**2)
+            + a1 * (360 + 672 * phi + 315 * phi**2)
+            + b1 * (290 + 540 * phi + 252 * phi**2)
         )
+        m9 = (
+            (66 + 115.5 * phi + 52.5 * phi**2)
+            + a1 * (45 + 75 * phi + 31.5 * phi**2)
+            + b1 * (32.5 + 52.5 * phi + 21 * phi**2)
+        ) * L
+        m10 = (
+            (12 + 21 * phi + 10.5 * phi**2)
+            + a1 * (7.5 + 12 * phi + 5.25 * phi**2)
+            + b1 * (5 + 7.5 * phi + 3 * phi**2)
+        ) * L**2
 
-        # element level matrix declaration
-        aux1 = self.material.rho * tempS * L / 420
         # fmt: off
+        M = np.array([
+            [ m1,   0, 0,   0,  m2, 0,  m3,   0, 0,   0, -m4, 0],
+            [  0,  m1, 0, -m2,   0, 0,   0,  m3, 0,  m4,   0, 0],
+            [  0,   0, 0,   0,   0, 0,   0,   0, 0,   0,   0, 0],
+            [  0, -m2, 0,  m5,   0, 0,   0, -m6, 0, -m7,   0, 0],
+            [ m2,   0, 0,   0,  m5, 0,  m6,   0, 0,   0, -m7, 0],
+            [  0,   0, 0,   0,   0, 0,   0,   0, 0,   0,   0, 0],
+            [ m3,   0, 0,   0,  m6, 0,  m8,   0, 0,   0, -m9, 0],
+            [  0,  m3, 0, -m6,   0, 0,   0,  m8, 0,  m9,   0, 0],
+            [  0,   0, 0,   0,   0, 0,   0,   0, 0,   0,   0, 0],
+            [  0,  m4, 0, -m7,   0, 0,   0,  m9, 0, m10,   0, 0],
+            [-m4,   0, 0,   0, -m7, 0, -m9,   0, 0,   0, m10, 0],
+            [  0,   0, 0,   0,   0, 0,   0,   0, 0,   0,   0, 0],
+        ]) * self.material.rho * self.A_l * L / (1260 * (1 + phi) ** 2)
+        # fmt: on
 
-        # Standard mass matrix
-        M = aux1 * np.array([
-            [  156,     0, 0,      0,  -22*L, 0,    54,     0, 0,      0,   13*L, 0],
-            [    0,   156, 0,   22*L,      0, 0,     0,    54, 0,  -13*L,      0, 0],
-            [    0,     0, 0,      0,      0, 0,     0,     0, 0,      0,      0, 0],
-            [    0,  22*L, 0, 4*L**2,      0, 0,     0,  13*L, 0,-3*L**2,      0, 0],
-            [-22*L,     0, 0,      0, 4*L**2, 0, -13*L,     0, 0,      0,-3*L**2, 0],
-            [    0,     0, 0,      0,      0, 0,     0,     0, 0,      0,      0, 0],
-            [   54,     0, 0,      0,  -13*L, 0,   156,     0, 0,      0,   22*L, 0],
-            [    0,    54, 0,   13*L,      0, 0,     0,   156, 0,  -22*L,      0, 0],
-            [    0,     0, 0,      0,      0, 0,     0,     0, 0,      0,      0, 0],
-            [    0, -13*L, 0,-3*L**2,      0, 0,     0, -22*L, 0, 4*L**2,      0, 0],
-            [ 13*L,     0, 0,      0,-3*L**2, 0,  22*L,     0, 0,      0, 4*L**2, 0],
-            [    0,     0, 0,      0,      0, 0,     0,     0, 0,      0,      0, 0],
-        ])
+        if self.rotary_inertia:
+            # fmt: off
+            m11 = 252 + 126 * a2 + 72 * b2 + 45 * gama + 30 * delta
+            m12 = (
+                21 - 105 * phi
+                + a2 * (21 - 42 * phi)
+                + b2 * (15 - 21 * phi)
+                + gama * (10.5 - 12 * phi)
+                + delta * (7.5 - 7.5 * phi)
+            ) * L
+            m13 = (
+                21 - 105 * phi
+                - 63 * a2 * phi
+                - b2 * (6 + 42 * phi)
+                - gama * (7.5 + 30 * phi)
+                - delta * (7.5 + 22.5 * phi)
+            ) * L
+            m14 = (
+                28 + 35 * phi + 70 * phi ** 2
+                + a2 * (7 - 7 * phi + 17.5 * phi ** 2)
+                + b2 * (4 - 7 * phi + 7 * phi ** 2)
+                + gama * (2.75 - 5 * phi + 3.5 * phi ** 2)
+                + delta * (2 - 3.5 * phi + 2 * phi ** 2)
+            ) * L ** 2
+            m15 = (
+                7 + 35 * phi - 35 * phi ** 2
+                + a2 * (3.5 + 17.5 * phi - 17.5 * phi ** 2)
+                + b2 * (3 + 10.5 * phi - 10.5 * phi ** 2)
+                + gama * (2.75 + 7 * phi - 7 * phi ** 2)
+                + delta * (2.5 + 5 * phi - 5 * phi ** 2)
+            ) * L ** 2
+            m16 = (
+                28 + 35 * phi + 70 * phi ** 2
+                + a2 * (21 + 42 * phi + 52.5 * phi ** 2)
+                + b2 * (18 + 42 * phi + 42 * phi ** 2)
+                + gama * (16.25 + 40 * phi + 35 * phi ** 2)
+                + delta * (15 + 37.5 * phi + 30 * phi ** 2)
+            ) * L ** 2
 
-        # Secondary inertias mass matrix
-        Ms = self.material.rho * tempI / (30 * L) * np.array([
-            [  36,   0, 0,     0,  -3*L, 0, -36,    0, 0,     0,  -3*L, 0],
-            [   0,  36, 0,   3*L,     0, 0,   0,  -36, 0,   3*L,     0, 0],
-            [   0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
-            [   0, 3*L, 0,4*L**2,     0, 0,   0, -3*L, 0, -L**2,     0, 0],
-            [-3*L,   0, 0,     0,4*L**2, 0, 3*L,    0, 0,     0, -L**2, 0],
-            [   0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
-            [ -36,   0, 0,     0,   3*L, 0,  36,    0, 0,     0,   3*L, 0],
-            [   0, -36, 0,  -3*L,     0, 0,   0,   36, 0,  -3*L,     0, 0],
-            [   0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
-            [   0, 3*L, 0, -L**2,     0, 0,   0, -3*L, 0,4*L**2,     0, 0],
-            [-3*L,   0, 0,     0, -L**2, 0, 3*L,    0, 0,     0,4*L**2, 0],
-            [   0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
-        ])
+            Mr = np.array([
+                [ m11,    0, 0,    0,  m12, 0, -m11,    0, 0,    0,  m13, 0],
+                [   0,  m11, 0, -m12,    0, 0,    0, -m11, 0, -m13,    0, 0],
+                [   0,    0, 0,    0,    0, 0,    0,    0, 0,    0,    0, 0],
+                [   0, -m12, 0,  m14,    0, 0,    0,  m12, 0, -m15,    0, 0],
+                [ m12,    0, 0,    0,  m14, 0, -m12,    0, 0,    0, -m15, 0],
+                [   0,    0, 0,    0,    0, 0,    0,    0, 0,    0,    0, 0],
+                [-m11,    0, 0,    0, -m12, 0,  m11,    0, 0,    0, -m13, 0],
+                [   0, -m11, 0,  m12,    0, 0,    0,  m11, 0,  m13,    0, 0],
+                [   0,    0, 0,    0,    0, 0,    0,    0, 0,    0,    0, 0],
+                [   0, -m13, 0, -m15,    0, 0,    0,  m13, 0,  m16,    0, 0],
+                [ m13,    0, 0,    0, -m15, 0, -m13,    0, 0,    0,  m16, 0],
+                [   0,    0, 0,    0,    0, 0,    0,    0, 0,    0,    0, 0],
+            ]) * self.material.rho * self.Ie_l / (210 * L * (1 + phi) ** 2)
+            # fmt: on
+            M += Mr
 
-        # Axial terms inertia matrix
-        Ma = self.material.rho * tempS * L / 6 * np.array([
+        # axial motion
+        Ae = (self.A_l + self.A_r) / 2
+        # fmt: off
+        Mam = np.array([
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
             [0, 0, 2, 0, 0, 0, 0, 0, 1, 0, 0, 0],
@@ -1663,10 +1541,14 @@ def M(self):
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-        ])
+        ]) * self.material.rho * Ae * L / 6
+        # fmt: on
+        M += Mam
 
-        # Torsional terms inertias matrix
-        Mr = self.material.rho * tempI * L / 6 * np.array([
+        # torsion
+        Je = 2 * (self.Ie_l + self.Ie_r) / 2
+        # fmt: off
+        Mts = np.array([
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
@@ -1679,11 +1561,9 @@ def M(self):
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
             [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 2],
-        ])
-
+        ]) * self.material.rho * Je * L / 6
         # fmt: on
-
-        M = M + Ms + Ma + Mr
+        M += Mts
 
         return M
 
@@ -1702,101 +1582,172 @@ def K(self):
         >>> Timoshenko_Element.K().shape
         (12, 12)
         """
-        # Axial force and torque applied to the element.
-        Fa = self.axial_force
-        T = self.torque
-
-        # temporary material and geometrical constants, determined as mean values
-        # from the left and right radii of the tapered shaft
         L = self.L
-        tempS = np.pi * (
-            ((self.odr / 2) ** 2 + (self.odl / 2) ** 2) / 2
-            - ((self.idr / 2) ** 2 + (self.idl / 2) ** 2) / 2
-        )
-        tempI = (
-            np.pi
-            / 4
-            * (
-                ((self.odr / 2) ** 4 + (self.odl / 2) ** 4) / 2
-                - ((self.idr / 2) ** 4 + (self.idl / 2) ** 4) / 2
-            )
-        )
-        tempJ = (
-            np.pi
-            / 2
-            * (
-                ((self.odr / 2) ** 4 + (self.odl / 2) ** 4) / 2
-                - ((self.idr / 2) ** 4 + (self.idl / 2) ** 4) / 2
-            )
-        )
+        phi = self.phi
+        a1 = self.a1
+        a2 = self.a2
+        b1 = self.b1
+        b2 = self.b2
+        delta = self.delta
+        gama = self.gama
 
-        # temporary variables
-        A = (
-            12
-            * self.material.E
-            * tempI
-            / (self.material.G_s * self.kappa * tempS * L**2)
-        )
+        # fmt: off
+        k1 = 1260 + 630 * a2 + 504 * b2 + 441 * gama + 396 * delta
+        k2 = (
+            630
+            + 210 * a2
+            + 147 * b2
+            + 126 * gama
+            + 114 * delta
+            - phi * (105 * a2 + 105 * b2 + 94.5 * gama + 84 * delta)
+        ) * L
+        k3 = (
+            630
+            + 420 * a2
+            + 357 * b2
+            + 315 * gama
+            + 282 * delta
+            + phi * (105 * a2 + 105 * b2 + 94.5 * gama + 84 * delta)
+        ) * L
+        k4 = (
+            420 + 210 * phi + 105 * phi ** 2
+            + a2 * (105 + 52.5 * phi ** 2)
+            + b2 * (56 - 35 * phi + 35 * phi ** 2)
+            + gama * (42 - 42 * phi + 26.25 * phi ** 2)
+            + delta * (36 - 42 * phi + 21 * phi ** 2)
+        ) * L ** 2
+        k5 = (
+            210 - 210 * phi - 105 * phi ** 2
+            + a2 * (105 - 105 * phi - 52.5 * phi ** 2)
+            + b2 * (91 - 70 * phi - 35 * phi ** 2)
+            + gama * (84 - 52.5 * phi - 26.25 * phi ** 2)
+            + delta * (78 - 42 * phi - 21 * phi ** 2)
+        ) * L ** 2
+        k6 = (
+            420 + 210 * phi + 105 * phi ** 2
+            + a2 * (315 + 210 * phi + 52.5 * phi ** 2)
+            + b2 * (266 + 175 * phi + 35 * phi ** 2)
+            + gama * (231 + 147 * phi + 26.25 * phi ** 2)
+            + delta * (204 + 126 * phi + 21 * phi ** 2)
+        ) * L ** 2
+        k7 = 12 + 6 * a1 + 4 * b1
+        k8 = (6 + 3 * a1 + 2 * b1) * L
+        k9 = (3 + 1.5 * a1 + b1) * L ** 2
+
+        K1 = np.array([
+            [ k1,   0,  0,   0,  k2,  0, -k1,   0,  0,   0,  k3,  0],
+            [  0,  k1,  0, -k2,   0,  0,   0, -k1,  0, -k3,   0,  0],
+            [  0,   0,  0,   0,   0,  0,   0,   0,  0,   0,   0,  0],
+            [  0, -k2,  0,  k4,   0,  0,   0,  k2,  0,  k5,   0,  0],
+            [ k2,   0,  0,   0,  k4,  0, -k2,   0,  0,   0,  k5,  0],
+            [  0,   0,  0,   0,   0,  0,   0,   0,  0,   0,   0,  0],
+            [-k1,   0,  0,   0, -k2,  0,  k1,   0,  0,   0, -k3,  0],
+            [  0, -k1,  0,  k2,   0,  0,   0,  k1,  0,  k3,   0,  0],
+            [  0,   0,  0,   0,   0,  0,   0,   0,  0,   0,   0,  0],
+            [  0, -k3,  0,  k5,   0,  0,   0,  k3,  0,  k6,   0,  0],
+            [ k3,   0,  0,   0,  k5,  0, -k3,   0,  0,   0,  k6,  0],
+            [  0,   0,  0,   0,   0,  0,   0,   0,  0,   0,   0,  0],
+        ]) * self.Ie_l / 105 
+
+        K2 = np.array([
+            [ k7,   0,  0,   0,  k8,  0, -k7,   0,  0,   0,  k8,  0],
+            [  0,  k7,  0, -k8,   0,  0,   0, -k7,  0, -k8,   0,  0],
+            [  0,   0,  0,   0,   0,  0,   0,   0,  0,   0,   0,  0],
+            [  0, -k8,  0,  k9,   0,  0,   0,  k8,  0,  k9,   0,  0],
+            [ k8,   0,  0,   0,  k9,  0, -k8,   0,  0,   0,  k9,  0],
+            [  0,   0,  0,   0,   0,  0,   0,   0,  0,   0,   0,  0],
+            [-k7,   0,  0,   0, -k8,  0,  k7,   0,  0,   0, -k8,  0],
+            [  0, -k7,  0,  k8,   0,  0,   0,  k7,  0,  k8,   0,  0],
+            [  0,   0,  0,   0,   0,  0,   0,   0,  0,   0,   0,  0],
+            [  0, -k8,  0,  k9,   0,  0,   0,  k8,  0,  k9,   0,  0],
+            [ k8,   0,  0,   0,  k9,  0, -k8,   0,  0,   0,  k9,  0],
+            [  0,   0,  0,   0,   0,  0,   0,   0,  0,   0,   0,  0],
+        ]) * self.Ie * phi * self.A_l / self.A
+        # fmt: on
 
-        # auxiliary variables
-        a1 = self.material.E * tempI / ((1 + A) * L**3)
-        a2 = self.material.G_s * tempJ / L
-        a3 = self.material.E * tempS / L
+        K = self.material.E * L ** (-3) * (1 + phi) ** (-2) * (K1 + K2)
 
+        # axial force
+        k10 = 36 + 60 * phi + 30 * phi**2
+        k11 = 3 * L
+        k12 = (4 + 5 * phi + 2.5 * phi**2) * L**2
+        k13 = (1 + 5 * phi + 2.5 * phi**2) * L**2
         # fmt: off
-        # pure stiffness matrix [Kc], added to the axial loads stiffness matrix [Ka],
-        # torsional stiffness matrix [Kr] and Timoshenko shear compensation [Ks].
-        Kc_plus = a1 * np.array([
-            [  12,   0,      0,            0,         -6*L,      0, -12,    0,      0,            0,         -6*L,      0],
-            [   0,  12,      0,          6*L,            0,      0,   0,  -12,      0,          6*L,            0,      0],
-            [   0,   0,  a3/a1,            0,            0,      0,   0,    0, -a3/a1,            0,            0,      0],
-            [   0, 6*L,      0, L**2*(A + 4),            0,      0,   0, -6*L,      0,-L**2*(A - 2),            0,      0],
-            [-6*L,   0,      0,            0, L**2*(A + 4),      0, 6*L,    0,      0,            0,-L**2*(A - 2),      0],
-            [   0,   0,      0,            0,            0,  a2/a1,   0,    0,      0,            0,            0, -a2/a1],
-            [ -12,   0,      0,            0,          6*L,      0,  12,    0,      0,            0,          6*L,      0],
-            [   0, -12,      0,         -6*L,            0,      0,   0,   12,      0,         -6*L,            0,      0],
-            [   0,   0, -a3/a1,            0,            0,      0,   0,    0,  a3/a1,            0,            0,      0],
-            [   0, 6*L,      0,-L**2*(A - 2),            0,      0,   0, -6*L,      0, L**2*(A + 4),            0,      0],
-            [-6*L,   0,      0,            0,-L**2*(A - 2),      0, 6*L,    0,      0,            0, L**2*(A + 4),      0],
-            [   0,   0,      0,            0,            0, -a2/a1,   0,    0,      0,            0,            0,  a2/a1],
-        ])
+        Kaf = np.array([
+            [ k10,    0,  0,    0,  k11,  0, -k10,    0,  0,    0,  k11,  0],
+            [   0,  k10,  0, -k11,    0,  0,    0, -k10,  0, -k11,    0,  0],
+            [   0,    0,  0,    0,    0,  0,    0,    0,  0,    0,    0,  0],
+            [   0, -k11,  0,  k12,    0,  0,    0,  k11,  0, -k13,    0,  0],
+            [ k11,    0,  0,    0,  k12,  0, -k11,    0,  0,    0, -k13,  0],
+            [   0,    0,  0,    0,    0,  0,    0,    0,  0,    0,    0,  0],
+            [-k10,    0,  0,    0, -k11,  0,  k10,    0,  0,    0, -k11,  0],
+            [   0, -k10,  0,  k11,    0,  0,    0,  k10,  0,  k11,    0,  0],
+            [   0,    0,  0,    0,    0,  0,    0,    0,  0,    0,    0,  0],
+            [   0, -k11,  0, -k13,    0,  0,    0,  k11,  0,  k12,    0,  0],
+            [ k11,    0,  0,    0, -k13,  0, -k11,    0,  0,    0,  k12,  0],
+            [   0,    0,  0,    0,    0,  0,    0,    0,  0,    0,    0,  0],
+        ]) * self.axial_force / (30 * L * (1 + phi) ** 2)
+        # fmt: on
+        K += Kaf
 
-        # stiffness matrix due to axial loading influence
-        Kf = Fa / (30 * L) * np.array([
-            [  36,   0, 0,     0,  -3*L, 0, -36,    0, 0,     0,  -3*L, 0],
-            [   0,  36, 0,   3*L,     0, 0,   0,  -36, 0,   3*L,     0, 0],
-            [   0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
-            [   0, 3*L, 0,4*L**2,     0, 0,   0, -3*L, 0, -L**2,     0, 0],
-            [-3*L,   0, 0,     0,4*L**2, 0, 3*L,    0, 0,     0, -L**2, 0],
-            [   0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
-            [ -36,   0, 0,     0,   3*L, 0,  36,    0, 0,     0,   3*L, 0],
-            [   0, -36, 0,  -3*L,     0, 0,   0,   36, 0,  -3*L,     0, 0],
-            [   0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
-            [   0, 3*L, 0, -L**2,     0, 0,   0, -3*L, 0,4*L**2,     0, 0],
-            [-3*L,   0, 0,     0, -L**2, 0, 3*L,    0, 0,     0,4*L**2, 0],
-            [   0,   0, 0,     0,     0, 0,   0,    0, 0,     0,     0, 0],
-        ])
+        # torque
+        # fmt: off
+        Ktq = np.array([
+            [ 0,  0,  0,    1,    0,  0,  0,  0,  0,   -1,    0,  0],
+            [ 0,  0,  0,    0,    1,  0,  0,  0,  0,    0,   -1,  0],
+            [ 0,  0,  0,    0,    0,  0,  0,  0,  0,    0,    0,  0],
+            [ 1,  0,  0,    0, -L/2,  0, -1,  0,  0,    0,  L/2,  0],
+            [ 0,  1,  0,  L/2,    0,  0,  0, -1,  0, -L/2,    0,  0],
+            [ 0,  0,  0,    0,    0,  0,  0,  0,  0,    0,    0,  0],
+            [ 0,  0,  0,   -1,    0,  0,  0,  0,  0,    1,    0,  0],
+            [ 0,  0,  0,    0,   -1,  0,  0,  0,  0,    0,    1,  0],
+            [ 0,  0,  0,    0,    0,  0,  0,  0,  0,    0,    0,  0],
+            [-1,  0,  0,    0, -L/2,  0,  1,  0,  0,    0,  L/2,  0],
+            [ 0, -1,  0,  L/2,    0,  0,  0,  1,  0, -L/2,    0,  0],
+            [ 0,  0,  0,    0,    0,  0,  0,  0,  0,    0,    0,  0],
+        ]) * self.torque / L
+        # fmt: on
+        K += Ktq
 
-        # stiffness matrix due to torque loading influence
-        Kt = T * np.array([
-            [   0,    0, 0, -1/L,    0, 0,    0,    0, 0,  1/L,    0, 0],
-            [   0,    0, 0,    0, -1/L, 0,    0,    0, 0,    0,  1/L, 0],
-            [   0,    0, 0,    0,    0, 0,    0,    0, 0,    0,    0, 0],
-            [-1/L,    0, 0,    0,  1/2, 0,  1/L,    0, 0,    0,  1/2, 0],
-            [   0, -1/L, 0, -1/2,    0, 0,    0,  1/L, 0, -1/2,    0, 0],
-            [   0,    0, 0,    0,    0, 0,    0,    0, 0,    0,    0, 0],
-            [   0,    0, 0,  1/L,    0, 0,    0,    0, 0, -1/L,    0, 0],
-            [   0,    0, 0,    0,  1/L, 0,    0,    0, 0,    0, -1/L, 0],
-            [   0,    0, 0,    0,    0, 0,    0,    0, 0,    0,    0, 0],
-            [ 1/L,    0, 0,    0, -1/2, 0, -1/L,    0, 0,    0, -1/2, 0],
-            [   0,  1/L, 0,  1/2,    0, 0,    0, -1/L, 0,  1/2,    0, 0],
-            [   0,    0, 0,    0,    0, 0,    0,    0, 0,    0,    0, 0],
-        ])
+        # axial motion
+        Ae = (self.A_l + self.A_r) / 2
+        # fmt: off
+        Kam = np.array([
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  1,  0,  0,  0,  0,  0, -1,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0, -1,  0,  0,  0,  0,  0,  1,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+        ]) * self.material.E * Ae / L
         # fmt: on
-        # Dynamic stiffness matrix is added independently in "def Kst"
-        # Kst = self.material.rho*tempI/(15*L)*np.array(12,12)
+        K += Kam
 
-        K = Kc_plus + Kf + Kt
+        # torsion
+        Je = 2 * (self.Ie_l + self.Ie_r) / 2
+        # fmt: off
+        Kts = np.array([
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0, -1],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [0,  0,  0,  0,  0, -1,  0,  0,  0,  0,  0,  1],
+        ]) * self.material.G_s * Je / L
+        # fmt: on
+        K += Kts
 
         return K
 
@@ -1819,34 +1770,24 @@ def Kst(self):
         >>> Timoshenko_Element.Kst().shape
         (12, 12)
         """
-        # temporary material and geometrical constants, determined as mean values
-        # from the left and right radii of the taperad shaft
         L = self.L
-        tempI = (
-            np.pi
-            / 4
-            * (
-                ((self.odr / 2) ** 4 + (self.odl / 2) ** 4) / 2
-                - ((self.idr / 2) ** 4 + (self.idl / 2) ** 4) / 2
-            )
-        )
+        Ie = (self.Ie_l + self.Ie_r) / 2
 
         # fmt: off
-        # dynamic stiffening matrix
-        Kst = self.material.rho * tempI / (15 * L) * np.array([
-            [0, -36, 0,   -3*L, 0, 0, 0,   36, 0,   -3*L, 0, 0],
-            [0,   0, 0,      0, 0, 0, 0,    0, 0,      0, 0, 0],
-            [0,   0, 0,      0, 0, 0, 0,    0, 0,      0, 0, 0],
-            [0,   0, 0,      0, 0, 0, 0,    0, 0,      0, 0, 0],
-            [0, 3*L, 0, 4*L**2, 0, 0, 0, -3*L, 0,  -L**2, 0, 0],
-            [0,   0, 0,      0, 0, 0, 0,    0, 0,      0, 0, 0],
-            [0,  36, 0,    3*L, 0, 0, 0,  -36, 0,    3*L, 0, 0],
-            [0,   0, 0,      0, 0, 0, 0,    0, 0,      0, 0, 0],
-            [0,   0, 0,      0, 0, 0, 0,    0, 0,      0, 0, 0],
-            [0,   0, 0,      0, 0, 0, 0,    0, 0,      0, 0, 0],
-            [0, 3*L, 0,  -L**2, 0, 0, 0, -3*L, 0, 4*L**2, 0, 0],
-            [0,   0, 0,      0, 0, 0, 0,    0, 0,      0, 0, 0],
-        ])
+        Kst = np.array([
+            [0,  -36,  0,    3*L,  0,  0,  0,  36,  0,    3*L,  0,  0],
+            [0,    0,  0,      0,  0,  0,  0,   0,  0,      0,  0,  0],
+            [0,    0,  0,      0,  0,  0,  0,   0,  0,      0,  0,  0],
+            [0,    0,  0,      0,  0,  0,  0,   0,  0,      0,  0,  0],
+            [0, -3*L,  0, 4*L**2,  0,  0,  0, 3*L,  0,  -L**2,  0,  0],
+            [0,    0,  0,      0,  0,  0,  0,   0,  0,      0,  0,  0],
+            [0,   36,  0,   -3*L,  0,  0,  0, -36,  0,   -3*L,  0,  0],
+            [0,    0,  0,      0,  0,  0,  0,   0,  0,      0,  0,  0],
+            [0,    0,  0,      0,  0,  0,  0,   0,  0,      0,  0,  0],
+            [0,    0,  0,      0,  0,  0,  0,   0,  0,      0,  0,  0],
+            [0, -3*L,  0,  -L**2,  0,  0,  0, 3*L,  0, 4*L**2,  0,  0],
+            [0,    0,  0,      0,  0,  0,  0,   0,  0,      0,  0,  0],
+        ]) * self.material.rho * Ie / (15 * L) 
         # fmt: on
 
         return Kst
@@ -1891,225 +1832,76 @@ def G(self):
         (12, 12)
         """
         if self.gyroscopic:
-            # temporary material and geometrical constants, determined as mean values
-            # from the left and right radii of the tapered shaft
+            phi = self.phi
             L = self.L
-            tempI = (
-                np.pi
-                / 4
-                * (
-                    ((self.odr / 2) ** 4 + (self.odl / 2) ** 4) / 2
-                    - ((self.idr / 2) ** 4 + (self.idl / 2) ** 4) / 2
-                )
-            )
-            tempS = np.pi * (
-                ((self.odr / 2) ** 2 + (self.odl / 2) ** 2) / 2
-                - ((self.idr / 2) ** 2 + (self.idl / 2) ** 2) / 2
-            )
+            a2 = self.a2
+            b2 = self.b2
+            delta = self.delta
+            gama = self.gama
 
             # fmt: off
-            # Gyroscopic effect matrix
-            phi=(
-            12
-            * self.material.E
-            * tempI
-            / (self.material.G_s * self.kappa * tempS * L ** 2)
-            )
+            g1 = 252 + 126 * a2 + 72 * b2 + 45 * gama + 30 * delta
+            g2 = (
+                21
+                - 105 * phi
+                + a2 * (21 - 42 * phi)
+                + b2 * (15 - 21 * phi)
+                + gama * (10.5 - 12 * phi)
+                + delta * (7.5 - 7.5 * phi)
+            ) * L
+            g3 = (
+                21
+                - 105 * phi
+                - 63 * a2 * phi
+                - b2 * (6 + 42 * phi)
+                - gama * (7.5 + 30 * phi)
+                - delta * (7.5 + 22.5 * phi)
+            ) * L
+            g4 = (
+                28
+                + 35 * phi
+                + 70 * phi ** 2
+                + a2 * (7 - 7 * phi + 17.5 * phi ** 2)
+                + b2 * (4 - 7 * phi + 7 * phi ** 2)
+                + gama * (2.75 - 5 * phi + 3.5 * phi ** 2)
+                + delta * (2 - 3.5 * phi + 2 * phi ** 2)
+            ) * L ** 2
+            g5 = (
+                7
+                + 35 * phi
+                - 35 * phi ** 2
+                + a2 * (3.5 + 17.5 * phi - 17.5 * phi ** 2)
+                + b2 * (3 + 10.5 * phi - 10.5 * phi ** 2)
+                + gama * (2.75 + 7 * phi - 7 * phi * 2)
+                + delta * (2.5 + 5 * phi - 5 * phi ** 2)
+            ) * L ** 2
+            g6 = (
+                28
+                + 35 * phi
+                + 70 * phi**2
+                + a2 * (21 + 42 * phi + 52.5 * phi ** 2)
+                + b2 * (18 + 42 * phi + 42 * phi ** 2)
+                + gama * (16.25 + 40 * phi + 35 * phi ** 2)
+                + delta * (15 + 37.5 * phi + 30 * phi ** 2)
+            ) * L ** 2
 
-            g1 = 36
-            g2 = (3 - 15*phi) * L
-            g3 = (4 + 5*phi + 10*phi**2) * L**2
-            g4 = (-1 - 5*phi + 5*phi**2) * L**2
-
-            G = (self.material.rho * tempI / (15*((1+phi)**2) * L)) * np.array([
-                    [   0,  g1, 0,   -g2,   0, 0,     0, -g1, 0,   -g2,   0, 0],
-                    [ -g1,   0, 0,     0, -g2, 0,    g1,   0, 0,     0, -g2, 0],
-                    [   0,   0, 0,     0,   0, 0,     0,   0, 0,     0,   0, 0],
-                    [  g2,   0, 0,     0,  g3, 0,   -g2,   0, 0,     0,  g4, 0],
-                    [   0,  g2, 0,   -g3,   0, 0,     0, -g2, 0,   -g4,   0, 0],
-                    [   0,   0, 0,     0,   0, 0,     0,   0, 0,     0,   0, 0],
-                    [   0, -g1, 0,    g2,   0, 0,     0,  g1, 0,    g2,   0, 0],
-                    [  g1,   0, 0,     0,  g2, 0,   -g1,   0, 0,     0,  g2, 0],
-                    [   0,   0, 0,     0,   0, 0,     0,   0, 0,     0,   0, 0],
-                    [  g2,   0, 0,     0,  g4, 0,   -g2,   0, 0,     0,  g3, 0],
-                    [   0,  g2, 0,   -g4,   0, 0,     0, -g2, 0,   -g3,   0, 0],
-                    [   0,   0, 0,     0,   0, 0,     0,   0, 0,     0,   0, 0],
-            ])
+            G = np.array([
+                [  0,  g1,  0, -g2,   0,  0,   0, -g1,  0, -g3,   0,  0],
+                [-g1,   0,  0,   0, -g2,  0,  g1,   0,  0,   0, -g3,  0],
+                [  0,   0,  0,   0,   0,  0,   0,   0,  0,   0,   0,  0],
+                [ g2,   0,  0,   0,  g4,  0, -g2,   0,  0,   0, -g5,  0],
+                [  0,  g2,  0, -g4,   0,  0,   0, -g2,  0,  g5,   0,  0],
+                [  0,   0,  0,   0,   0,  0,   0,   0,  0,   0,   0,  0],
+                [  0, -g1,  0,  g2,   0,  0,   0,  g1,  0,  g3,   0,  0],
+                [ g1,   0,  0,   0,  g2,  0, -g1,   0,  0,   0,  g3,  0],
+                [  0,   0,  0,   0,   0,  0,   0,   0,  0,   0,   0,  0],
+                [ g3,   0,  0,   0, -g5,  0, -g3,   0,  0,   0,  g6,  0],
+                [  0,  g3,  0,  g5,   0,  0,   0, -g3,  0, -g6,   0,  0],
+                [  0,   0,  0,   0,   0,  0,   0,   0,  0,   0,   0,  0],
+            ]) * self.material.rho * self.Ie_l * 2 / (210 * L * (1 + phi) ** 2)
             # fmt: on
+
         else:
             G = np.zeros((12, 12))
 
         return G
-
-    @classmethod
-    def from_table(cls, file, sheet_type="Simple", sheet_name=0):
-        """Instantiate one or more shafts using inputs from an Excel table.
-
-        A header with the names of the columns is required. These names should
-        match the names expected by the routine (usually the names of the
-        parameters, but also similar ones). The program will read every row
-        bellow the header until they end or it reaches a NaN.
-
-        Parameters
-        ----------
-        file : str
-            Path to the file containing the shaft parameters.
-        sheet_type : str, optional
-            Describes the kind of sheet the function should expect:
-                Simple: The input table should specify only the number of the materials
-                to be used.
-                They must be saved prior to calling the method.
-                Model: The materials parameters must be passed along with the shaft
-                parameters. Each material must have an id number and each shaft must
-                reference one of the materials ids.
-        sheet_name : int or str, optional
-            Position of the sheet in the file (starting from 0) or its name. If none is
-            passed, it is assumed to be the first sheet in the file.
-
-        Returns
-        -------
-        shaft: list
-            A list of shaft objects.
-        """
-        parameters = read_table_file(
-            file, "shaft", sheet_name=sheet_name, sheet_type=sheet_type
-        )
-        list_of_shafts = []
-        if sheet_type == "Model":
-            new_materials = {}
-            for i in range(0, len(parameters["matno"])):
-                new_material = Material(
-                    name="shaft_mat_" + str(parameters["matno"][i]),
-                    rho=parameters["rhoa"][i],
-                    E=parameters["ea"][i],
-                    G_s=parameters["ga"][i],
-                )
-                new_materials["shaft_mat_" + str(parameters["matno"][i])] = new_material
-            for i in range(0, len(parameters["L"])):
-                list_of_shafts.append(
-                    cls(
-                        L=parameters["L"][i],
-                        i_d=parameters["i_d"][i],
-                        o_d=parameters["o_d"][i],
-                        alpha=parameters["alpha"][i],
-                        beta=parameters["beta"][i],
-                        material=new_materials[parameters["material"][i]],
-                        n=parameters["n"][i],
-                        axial_force=parameters["axial_force"][i],
-                        torque=parameters["torque"][i],
-                        gyroscopic=parameters["gyroscopic"][i],
-                    )
-                )
-        elif sheet_type == "Simple":
-            for i in range(0, len(parameters["L"])):
-                list_of_shafts.append(
-                    cls(
-                        L=parameters["L"][i],
-                        i_d=parameters["i_d"][i],
-                        o_d=parameters["o_d"][i],
-                        alpha=parameters["alpha"][i],
-                        beta=parameters["beta"][i],
-                        material=parameters["material"][i],
-                        n=parameters["n"][i],
-                        axial_force=parameters["axial_force"][i],
-                        torque=parameters["torque"][i],
-                        gyroscopic=parameters["gyroscopic"][i],
-                    )
-                )
-        return list_of_shafts
-
-    @classmethod
-    def section(
-        cls,
-        L,
-        ne,
-        s_idl,
-        s_odl,
-        s_idr=None,
-        s_odr=None,
-        alpha=0,
-        beta=0,
-        material=None,
-        n=None,
-        gyroscopic=True,
-    ):
-        """Shaft section constructor.
-
-        This method will create a shaft section with length 'L' divided into
-        'ne' elements.
-
-        Parameters
-        ----------
-        odl : float, pint.Quantity
-            Outer diameter of the element at the left node.
-        idr : float, pint.Quantity, optional
-            Inner diameter of the element at the right node;
-            Default is equal to idl value for cylindrical element.
-        odr : float, pint.Quantity, optional
-            Outer diameter of the element at the right node;
-            Default is equal to odl value for cylindrical element.
-        E : float
-            Young's modulus.
-        G_s : float
-            Shear modulus.
-        alpha : float
-            Proportional damping coefficient, associated to the element Mass matrix
-        beta : float
-            Proportional damping coefficient, associated to the element Stiffness matrix
-        material : ross.material
-            Shaft material.
-        n : int, optional
-            Element number (coincident with it's first node).
-            If not given, it will be set when the rotor is assembled
-            according to the element's position in the list supplied to
-            the rotor constructor.
-        axial_force : float
-            Axial force.
-        torque : float
-            Torque.
-        gyroscopic : bool
-            Determine if gyroscopic effects are taken into account.
-            Default is False.
-
-        Returns
-        -------
-        elements : list
-            List with the 'ne' shaft elements.
-
-        Examples
-        --------
-        >>> # shaft material
-        >>> from ross.materials import steel
-        >>> # shaft inner and outer diameters
-        >>> s_idl = 0
-        >>> s_odl = 0.01585
-        >>> sec = ShaftElement.section(247.65e-3, 4, 0, 15.8e-3, material=steel)
-        >>> len(sec)
-        4
-        >>> sec[0].i_d
-        0.0
-        """
-        if s_idr is None:
-            s_idr = s_idl
-        if s_odr is None:
-            s_odr = s_odl
-
-        le = L / ne
-
-        elements = [
-            cls(
-                L=le,
-                idl=(s_idr - s_idl) * i * le / L + s_idl,
-                odl=(s_odr - s_odl) * i * le / L + s_odl,
-                idr=(s_idr - s_idl) * (i + 1) * le / L + s_idl,
-                odr=(s_odr - s_odl) * (i + 1) * le / L + s_odl,
-                alpha=alpha,
-                beta=beta,
-                material=material,
-                n=n,
-                gyroscopic=gyroscopic,
-            )
-            for i in range(ne)
-        ]
-
-        return elements


Test patch:
diff --git a/ross/tests/test_misalignment.py b/ross/tests/test_misalignment.py
index c46624706..a409d2935 100644
--- a/ross/tests/test_misalignment.py
+++ b/ross/tests/test_misalignment.py
@@ -1,82 +1,75 @@
-import os
-from pathlib import Path
-from tempfile import tempdir
-
 import numpy as np
 import pytest
-from numpy.testing import assert_allclose, assert_almost_equal
+from numpy.testing import assert_allclose
 
 import ross as rs
-from ross.defects.misalignment import MisalignmentFlex
 from ross.units import Q_
 
-steel2 = rs.Material(name="Steel", rho=7850, E=2.17e11, G_s=81.2e9)
 
-#  Rotor with 6 DoFs, with internal damping, with 10 shaft elements, 2 disks and 2 bearings.
-i_d = 0
-o_d = 0.019
-n = 33
+@pytest.fixture
+def rotor():
+    steel2 = rs.Material(name="Steel", rho=7850, E=2.17e11, Poisson=0.2992610837438423)
+
+    #  Rotor with 6 DoFs, with internal damping, with 10 shaft elements, 2 disks and 2 bearings.
+    i_d = 0
+    o_d = 0.019
 
-# fmt: off
-L = np.array(
-        [0  ,  25,  64, 104, 124, 143, 175, 207, 239, 271,
+    # fmt: off
+    L = np.array([
+        0  ,  25,  64, 104, 124, 143, 175, 207, 239, 271,
         303, 335, 345, 355, 380, 408, 436, 466, 496, 526,
         556, 586, 614, 647, 657, 667, 702, 737, 772, 807,
-        842, 862, 881, 914]
-        )/ 1000
-# fmt: on
+        842, 862, 881, 914
+    ])/ 1000
+    # fmt: on
 
-L = [L[i] - L[i - 1] for i in range(1, len(L))]
+    L = [L[i] - L[i - 1] for i in range(1, len(L))]
 
-shaft_elem = [
-    rs.ShaftElement6DoF(
-        material=steel2,
-        L=l,
-        idl=i_d,
-        odl=o_d,
-        idr=i_d,
-        odr=o_d,
-        alpha=8.0501,
-        beta=1.0e-5,
-        rotary_inertia=True,
-        shear_effects=True,
-    )
-    for l in L
-]
+    shaft_elem = [
+        rs.ShaftElement6DoF(
+            l,
+            i_d,
+            o_d,
+            material=steel2,
+            alpha=8.0501,
+            beta=1.0e-5,
+        )
+        for l in L
+    ]
 
-Id = 0.003844540885417
-Ip = 0.007513248437500
+    Id = 0.003844540885417
+    Ip = 0.007513248437500
 
-disk0 = rs.DiskElement6DoF(n=12, m=2.6375, Id=Id, Ip=Ip)
-disk1 = rs.DiskElement6DoF(n=24, m=2.6375, Id=Id, Ip=Ip)
+    disk0 = rs.DiskElement6DoF(n=12, m=2.6375, Id=Id, Ip=Ip)
+    disk1 = rs.DiskElement6DoF(n=24, m=2.6375, Id=Id, Ip=Ip)
 
-kxx1 = 4.40e5
-kyy1 = 4.6114e5
-kzz = 0
-cxx1 = 27.4
-cyy1 = 2.505
-czz = 0
-kxx2 = 9.50e5
-kyy2 = 1.09e8
-cxx2 = 50.4
-cyy2 = 100.4553
+    kxx1 = 4.40e5
+    kyy1 = 4.6114e5
+    kzz = 0
+    cxx1 = 27.4
+    cyy1 = 2.505
+    czz = 0
+    kxx2 = 9.50e5
+    kyy2 = 1.09e8
+    cxx2 = 50.4
+    cyy2 = 100.4553
 
-bearing0 = rs.BearingElement6DoF(
-    n=4, kxx=kxx1, kyy=kyy1, cxx=cxx1, cyy=cyy1, kzz=kzz, czz=czz
-)
-bearing1 = rs.BearingElement6DoF(
-    n=31, kxx=kxx2, kyy=kyy2, cxx=cxx2, cyy=cyy2, kzz=kzz, czz=czz
-)
+    bearing0 = rs.BearingElement6DoF(
+        n=4, kxx=kxx1, kyy=kyy1, kzz=kzz, cxx=cxx1, cyy=cyy1, czz=czz
+    )
+    bearing1 = rs.BearingElement6DoF(
+        n=31, kxx=kxx2, kyy=kyy2, kzz=kzz, cxx=cxx2, cyy=cyy2, czz=czz
+    )
 
-rotor = rs.Rotor(shaft_elem, [disk0, disk1], [bearing0, bearing1])
+    return rs.Rotor(shaft_elem, [disk0, disk1], [bearing0, bearing1])
 
 
 @pytest.fixture
-def mis_comb():
+def common_parameters():
     unbalance_magnitudet = np.array([5e-4, 0])
-    unbalance_phaset = np.array([-np.pi / 2, 0])
+    unbalance_phaset = Q_(np.array([-90.0, 0.0]), "degrees")
 
-    misalignment = rotor.run_misalignment(
+    return dict(
         coupling="flex",
         dt=0.1,
         tI=0,
@@ -89,42 +82,17 @@ def mis_comb():
         TD=0,
         TL=0,
         n1=0,
-        speed=125.66370614359172,
+        speed=Q_(1200, "RPM"),
         unbalance_magnitude=unbalance_magnitudet,
         unbalance_phase=unbalance_phaset,
-        mis_type="combined",
         print_progress=False,
     )
 
-    return misalignment
-
 
 @pytest.fixture
-def mis_comb_units():
-    unbalance_magnitudet = Q_(np.array([0.043398083107259365, 0]), "lb*in")
-    unbalance_phaset = Q_(np.array([-90.0, 0.0]), "degrees")
+def mis_comb(rotor, common_parameters):
 
-    misalignment = rotor.run_misalignment(
-        coupling="flex",
-        dt=0.1,
-        tI=0,
-        tF=5,
-        kd=40 * 10 ** (3),
-        ks=38 * 10 ** (3),
-        eCOUPx=2 * 10 ** (-4),
-        eCOUPy=2 * 10 ** (-4),
-        misalignment_angle=5 * np.pi / 180,
-        TD=0,
-        TL=0,
-        n1=0,
-        speed=Q_(1200, "RPM"),
-        unbalance_magnitude=unbalance_magnitudet,
-        unbalance_phase=unbalance_phaset,
-        mis_type="combined",
-        print_progress=False,
-    )
-
-    return misalignment
+    return rotor.run_misalignment(mis_type="combined", **common_parameters)
 
 
 def test_mis_comb_parameters(mis_comb):
@@ -142,718 +110,218 @@ def test_mis_comb_parameters(mis_comb):
     assert mis_comb.speed == 125.66370614359172
 
 
-def test_mis_comb_parameters_units(mis_comb_units):
-    assert mis_comb_units.dt == 0.1
-    assert mis_comb_units.tI == 0
-    assert mis_comb_units.tF == 5
-    assert mis_comb_units.kd == 40 * 10 ** (3)
-    assert mis_comb_units.ks == 38 * 10 ** (3)
-    assert mis_comb_units.eCOUPx == 2 * 10 ** (-4)
-    assert mis_comb_units.eCOUPy == 2 * 10 ** (-4)
-    assert mis_comb_units.misalignment_angle == 5 * np.pi / 180
-    assert mis_comb_units.TD == 0
-    assert mis_comb_units.TL == 0
-    assert mis_comb_units.n1 == 0
-    assert mis_comb_units.speed == 125.66370614359172
-
-
 def test_mis_comb_forces(mis_comb):
-    assert mis_comb.forces[mis_comb.n1 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [-4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_comb.forces[mis_comb.n1 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_comb.forces[mis_comb.n2 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_comb.forces[mis_comb.n2 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [-1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-
-def test_mis_comb_forces_units(mis_comb_units):
-    assert mis_comb_units.forces[mis_comb_units.n1 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [-4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,-4.40604748,-4.40604748,-4.40604748,-4.40604748,
-    -4.40604748,
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_comb_units.forces[mis_comb_units.n1 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,1.0821174,1.0821174,1.0821174,1.0821174,
-    1.0821174,
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_comb_units.forces[mis_comb_units.n2 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,4.40604748,4.40604748,4.40604748,4.40604748,
-    4.40604748,
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_comb_units.forces[mis_comb_units.n2 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [-1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,-1.0821174,-1.0821174,-1.0821174,-1.0821174,
-    -1.0821174,
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
+    # fmt: off
+    Fx_n1_comb = np.array([
+        -4.40605, -4.40605, -4.40605, -4.40605, -4.40605, -4.40605,
+        -4.40605, -4.40605, -4.40605, -4.40605, -4.40605, -4.40605,
+        -4.40605, -4.40605, -4.40605, -4.40605, -4.40605, -4.40605,
+        -4.40605, -4.40605, -4.40605, -4.40605, -4.40605, -4.40605,
+        -4.40605, -4.40605, -4.40605, -4.40605, -4.40605, -4.40605,
+        -4.40605, -4.40605, -4.40605, -4.40605, -4.40605, -4.40605,
+        -4.40605, -4.40605, -4.40605, -4.40605, -4.40605, -4.40605,
+        -4.40605, -4.40605, -4.40605, -4.40605, -4.40605, -4.40605,
+        -4.40605, -4.40605, -4.40605
+    ])
+
+    Fy_n1_comb = np.array([
+        1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212,
+        1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212,
+        1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212,
+        1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212,
+        1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212,
+        1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212,
+        1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212,
+        1.08212, 1.08212
+    ])
+
+    Fx_n2_comb = np.array([
+        4.40604748, 4.40604748, 4.40604748, 4.40604748, 4.40604748,
+        4.40604748, 4.40604748, 4.40604748, 4.40604748, 4.40604748,
+        4.40604748, 4.40604748, 4.40604748, 4.40604748, 4.40604748,
+        4.40604748, 4.40604748, 4.40604748, 4.40604748, 4.40604748,
+        4.40604748, 4.40604748, 4.40604748, 4.40604748, 4.40604748,
+        4.40604748, 4.40604748, 4.40604748, 4.40604748, 4.40604748,
+        4.40604748, 4.40604748, 4.40604748, 4.40604748, 4.40604748,
+        4.40604748, 4.40604748, 4.40604748, 4.40604748, 4.40604748,
+        4.40604748, 4.40604748, 4.40604748, 4.40604748, 4.40604748,
+        4.40604748, 4.40604748, 4.40604748, 4.40604748, 4.40604748,
+        4.40604748
+    ])
+
+    Fy_n2_comb = np.array([
+        -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
+        -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
+        -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
+        -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
+        -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
+        -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
+        -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
+        -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
+        -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
+        -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
+        -1.0821174
+    ])
+    # fmt: on
+
+    assert_allclose(mis_comb.forces[mis_comb.n1 * 6 + 0, :], Fx_n1_comb, rtol=3e-2)
+    assert_allclose(mis_comb.forces[mis_comb.n1 * 6 + 1, :], Fy_n1_comb, rtol=3e-2)
+
+    assert_allclose(mis_comb.forces[mis_comb.n2 * 6 + 0, :], Fx_n2_comb, rtol=3e-2)
+    assert_allclose(mis_comb.forces[mis_comb.n2 * 6 + 1, :], Fy_n2_comb, rtol=3e-2)
 
 
 @pytest.fixture
-def mis_parallel():
-    unbalance_magnitudet = np.array([5e-4, 0])
-    unbalance_phaset = np.array([-np.pi / 2, 0])
-
-    misalignment = rotor.run_misalignment(
-        coupling="flex",
-        dt=0.1,
-        tI=0,
-        tF=5,
-        kd=40 * 10 ** (3),
-        ks=38 * 10 ** (3),
-        eCOUPx=2 * 10 ** (-4),
-        eCOUPy=2 * 10 ** (-4),
-        misalignment_angle=5 * np.pi / 180,
-        TD=0,
-        TL=0,
-        n1=0,
-        speed=125.66370614359172,
-        unbalance_magnitude=unbalance_magnitudet,
-        unbalance_phase=unbalance_phaset,
-        mis_type="parallel",
-        print_progress=False,
-    )
-
-    return misalignment
+def mis_par(rotor, common_parameters):
+
+    return rotor.run_misalignment(mis_type="parallel", **common_parameters)
+
+
+def test_mis_par_parameters(mis_par):
+    assert mis_par.dt == 0.1
+    assert mis_par.tI == 0
+    assert mis_par.tF == 5
+    assert mis_par.kd == 40 * 10 ** (3)
+    assert mis_par.ks == 38 * 10 ** (3)
+    assert mis_par.eCOUPx == 2 * 10 ** (-4)
+    assert mis_par.eCOUPy == 2 * 10 ** (-4)
+    assert mis_par.misalignment_angle == 5 * np.pi / 180
+    assert mis_par.TD == 0
+    assert mis_par.TL == 0
+    assert mis_par.n1 == 0
+    assert mis_par.speed == 125.66370614359172
+
+
+def test_mis_par_forces(mis_par):
+    # fmt: off
+    Fx_n1_par = np.array([
+        -6.78313, -6.78313, -6.78313, -6.78313, -6.78313, -6.78313,
+        -6.78313, -6.78313, -6.78313, -6.78313, -6.78313, -6.78313,
+        -6.78313, -6.78313, -6.78313, -6.78313, -6.78313, -6.78313,
+        -6.78313, -6.78313, -6.78313, -6.78313, -6.78313, -6.78313,
+        -6.78313, -6.78313, -6.78313, -6.78313, -6.78313, -6.78313,
+        -6.78313, -6.78313, -6.78313, -6.78313, -6.78313, -6.78313,
+        -6.78313, -6.78313, -6.78313, -6.78313, -6.78313, -6.78313,
+        -6.78313, -6.78313, -6.78313, -6.78313, -6.78313, -6.78313,
+        -6.78313, -6.78313, -6.78313
+    ])
+
+    Fy_n1_par = np.array([
+        1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212,
+        1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212,
+        1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212,
+        1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212,
+        1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212,
+        1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212,
+        1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212, 1.08212,
+        1.08212, 1.08212
+    ])
+
+    Fx_n2_par = np.array([
+        6.78313, 6.78313, 6.78313, 6.78313, 6.78313, 6.78313, 6.78313,
+        6.78313, 6.78313, 6.78313, 6.78313, 6.78313, 6.78313, 6.78313,
+        6.78313, 6.78313, 6.78313, 6.78313, 6.78313, 6.78313, 6.78313,
+        6.78313, 6.78313, 6.78313, 6.78313, 6.78313, 6.78313, 6.78313,
+        6.78313, 6.78313, 6.78313, 6.78313, 6.78313, 6.78313, 6.78313,
+        6.78313, 6.78313, 6.78313, 6.78313, 6.78313, 6.78313, 6.78313,
+        6.78313, 6.78313, 6.78313, 6.78313, 6.78313, 6.78313, 6.78313,
+        6.78313, 6.78313
+    ])
+
+    Fy_n2_par = np.array([
+        -1.08212, -1.08212, -1.08212, -1.08212, -1.08212, -1.08212,
+        -1.08212, -1.08212, -1.08212, -1.08212, -1.08212, -1.08212,
+        -1.08212, -1.08212, -1.08212, -1.08212, -1.08212, -1.08212,
+        -1.08212, -1.08212, -1.08212, -1.08212, -1.08212, -1.08212,
+        -1.08212, -1.08212, -1.08212, -1.08212, -1.08212, -1.08212,
+        -1.08212, -1.08212, -1.08212, -1.08212, -1.08212, -1.08212,
+        -1.08212, -1.08212, -1.08212, -1.08212, -1.08212, -1.08212,
+        -1.08212, -1.08212, -1.08212, -1.08212, -1.08212, -1.08212,
+        -1.08212, -1.08212, -1.08212
+    ])
+    # fmt: on
+
+    assert_allclose(mis_par.forces[mis_par.n1 * 6 + 0, :], Fx_n1_par, rtol=3e-2)
+    assert_allclose(mis_par.forces[mis_par.n1 * 6 + 1, :], Fy_n1_par, rtol=3e-2)
+
+    assert_allclose(mis_par.forces[mis_par.n2 * 6 + 0, :], Fx_n2_par, rtol=3e-2)
+    assert_allclose(mis_par.forces[mis_par.n2 * 6 + 1, :], Fy_n2_par, rtol=3e-2)
 
 
 @pytest.fixture
-def mis_parallel_units():
-    unbalance_magnitudet = Q_(np.array([0.043398083107259365, 0]), "lb*in")
-    unbalance_phaset = Q_(np.array([-90.0, 0.0]), "degrees")
-
-    misalignment = rotor.run_misalignment(
-        coupling="flex",
-        dt=0.1,
-        tI=0,
-        tF=5,
-        kd=40 * 10 ** (3),
-        ks=38 * 10 ** (3),
-        eCOUPx=2 * 10 ** (-4),
-        eCOUPy=2 * 10 ** (-4),
-        misalignment_angle=5 * np.pi / 180,
-        TD=0,
-        TL=0,
-        n1=0,
-        speed=Q_(1200, "RPM"),
-        unbalance_magnitude=unbalance_magnitudet,
-        unbalance_phase=unbalance_phaset,
-        mis_type="parallel",
-        print_progress=False,
-    )
-
-    return misalignment
-
-
-def test_mis_parallel_parameters(mis_parallel):
-    assert mis_parallel.dt == 0.1
-    assert mis_parallel.tI == 0
-    assert mis_parallel.tF == 5
-    assert mis_parallel.kd == 40 * 10 ** (3)
-    assert mis_parallel.ks == 38 * 10 ** (3)
-    assert mis_parallel.eCOUPx == 2 * 10 ** (-4)
-    assert mis_parallel.eCOUPy == 2 * 10 ** (-4)
-    assert mis_parallel.misalignment_angle == 5 * np.pi / 180
-    assert mis_parallel.TD == 0
-    assert mis_parallel.TL == 0
-    assert mis_parallel.n1 == 0
-    assert mis_parallel.speed == 125.66370614359172
-
-
-def test_mis_parallel_parameters_units(mis_parallel_units):
-    assert mis_parallel_units.dt == 0.1
-    assert mis_parallel_units.tI == 0
-    assert mis_parallel_units.tF == 5
-    assert mis_parallel_units.kd == 40 * 10 ** (3)
-    assert mis_parallel_units.ks == 38 * 10 ** (3)
-    assert mis_parallel_units.eCOUPx == 2 * 10 ** (-4)
-    assert mis_parallel_units.eCOUPy == 2 * 10 ** (-4)
-    assert mis_parallel_units.misalignment_angle == 5 * np.pi / 180
-    assert mis_parallel_units.TD == 0
-    assert mis_parallel_units.TL == 0
-    assert mis_parallel_units.n1 == 0
-    assert mis_parallel_units.speed == 125.66370614359172
-
-
-def test_mis_parallel_forces(mis_parallel):
-    assert mis_parallel.forces[mis_parallel.n1 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [-6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_parallel.forces[mis_parallel.n1 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_parallel.forces[mis_parallel.n2 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_parallel.forces[mis_parallel.n2 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [-1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-
-def test_mis_parallel_forces_units(mis_parallel_units):
-    assert mis_parallel_units.forces[mis_parallel_units.n1 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [-6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529, -6.78312529, -6.78312529, -6.78312529, -6.78312529,
-    -6.78312529
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_parallel_units.forces[mis_parallel_units.n1 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174, 1.0821174,
-    1.0821174, 1.0821174, 1.0821174
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_parallel_units.forces[mis_parallel_units.n2 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529, 6.78312529, 6.78312529, 6.78312529, 6.78312529,
-    6.78312529
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_parallel_units.forces[mis_parallel_units.n2 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [-1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174, -1.0821174, -1.0821174, -1.0821174, -1.0821174,
-    -1.0821174
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
+def mis_ang(rotor, common_parameters):
+
+    return rotor.run_misalignment(mis_type="angular", **common_parameters)
+
+
+def test_mis_ang_parameters(mis_ang):
+    assert mis_ang.dt == 0.1
+    assert mis_ang.tI == 0
+    assert mis_ang.tF == 5
+    assert mis_ang.kd == 40 * 10 ** (3)
+    assert mis_ang.ks == 38 * 10 ** (3)
+    assert mis_ang.eCOUPx == 2 * 10 ** (-4)
+    assert mis_ang.eCOUPy == 2 * 10 ** (-4)
+    assert mis_ang.misalignment_angle == 5 * np.pi / 180
+    assert mis_ang.TD == 0
+    assert mis_ang.TL == 0
+    assert mis_ang.n1 == 0
+    assert mis_ang.speed == 125.66370614359172
+
+
+def test_mis_ang_forces(mis_ang):
+    # fmt: off
+    Fx_n1_ang = np.array([
+        2.37708, 2.37708, 2.37708, 2.37708, 2.37708, 2.37708, 2.37708,
+        2.37708, 2.37708, 2.37708, 2.37708, 2.37708, 2.37708, 2.37708,
+        2.37708, 2.37708, 2.37708, 2.37708, 2.37708, 2.37708, 2.37708,
+        2.37708, 2.37708, 2.37708, 2.37708, 2.37708, 2.37708, 2.37708,
+        2.37708, 2.37708, 2.37708, 2.37708, 2.37708, 2.37708, 2.37708,
+        2.37708, 2.37708, 2.37708, 2.37708, 2.37708, 2.37708, 2.37708,
+        2.37708, 2.37708, 2.37708, 2.37708, 2.37708, 2.37708, 2.37708,
+        2.37708, 2.37708
+    ])
+
+    Fy_n1_ang = np.array([
+        -0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
+         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
+         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
+         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.
+    ])
+
+    Fx_n2_ang = np.array([
+        -2.37708, -2.37708, -2.37708, -2.37708, -2.37708, -2.37708,
+        -2.37708, -2.37708, -2.37708, -2.37708, -2.37708, -2.37708,
+        -2.37708, -2.37708, -2.37708, -2.37708, -2.37708, -2.37708,
+        -2.37708, -2.37708, -2.37708, -2.37708, -2.37708, -2.37708,
+        -2.37708, -2.37708, -2.37708, -2.37708, -2.37708, -2.37708,
+        -2.37708, -2.37708, -2.37708, -2.37708, -2.37708, -2.37708,
+        -2.37708, -2.37708, -2.37708, -2.37708, -2.37708, -2.37708,
+        -2.37708, -2.37708, -2.37708, -2.37708, -2.37708, -2.37708,
+        -2.37708, -2.37708, -2.37708
+    ])
+
+    Fy_n2_ang = np.array([
+         0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
+        -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
+        -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
+        -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.
+    ])
+    # fmt: on
+
+    assert_allclose(mis_ang.forces[mis_ang.n1 * 6 + 0, :], Fx_n1_ang, rtol=3e-2)
+    assert_allclose(mis_ang.forces[mis_ang.n1 * 6 + 1, :], Fy_n1_ang, atol=1e-6)
+
+    assert_allclose(mis_ang.forces[mis_ang.n2 * 6 + 0, :], Fx_n2_ang, rtol=3e-2)
+    assert_allclose(mis_ang.forces[mis_ang.n2 * 6 + 1, :], Fy_n2_ang, atol=1e-6)
 
 
 @pytest.fixture
-def mis_angular():
+def mis_rigid(rotor):
     unbalance_magnitudet = np.array([5e-4, 0])
-    unbalance_phaset = np.array([-np.pi / 2, 0])
-
-    misalignment = rotor.run_misalignment(
-        coupling="flex",
-        dt=0.1,
-        tI=0,
-        tF=5,
-        kd=40 * 10 ** (3),
-        ks=38 * 10 ** (3),
-        eCOUPx=2 * 10 ** (-4),
-        eCOUPy=2 * 10 ** (-4),
-        misalignment_angle=5 * np.pi / 180,
-        TD=0,
-        TL=0,
-        n1=0,
-        speed=125.66370614359172,
-        unbalance_magnitude=unbalance_magnitudet,
-        unbalance_phase=unbalance_phaset,
-        mis_type="angular",
-        print_progress=False,
-    )
-
-    return misalignment
-
-
-@pytest.fixture
-def mis_angular_units():
-    unbalance_magnitudet = Q_(np.array([0.043398083107259365, 0]), "lb*in")
     unbalance_phaset = Q_(np.array([-90.0, 0.0]), "degrees")
 
-    misalignment = rotor.run_misalignment(
-        coupling="flex",
-        dt=0.1,
-        tI=0,
-        tF=5,
-        kd=40 * 10 ** (3),
-        ks=38 * 10 ** (3),
-        eCOUPx=2 * 10 ** (-4),
-        eCOUPy=2 * 10 ** (-4),
-        misalignment_angle=5 * np.pi / 180,
-        TD=0,
-        TL=0,
-        n1=0,
-        speed=Q_(1200, "RPM"),
-        unbalance_magnitude=unbalance_magnitudet,
-        unbalance_phase=unbalance_phaset,
-        mis_type="angular",
-        print_progress=False,
-    )
-
-    return misalignment
-
-
-def test_mis_angular_parameters(mis_angular):
-    assert mis_angular.dt == 0.1
-    assert mis_angular.tI == 0
-    assert mis_angular.tF == 5
-    assert mis_angular.kd == 40 * 10 ** (3)
-    assert mis_angular.ks == 38 * 10 ** (3)
-    assert mis_angular.eCOUPx == 2 * 10 ** (-4)
-    assert mis_angular.eCOUPy == 2 * 10 ** (-4)
-    assert mis_angular.misalignment_angle == 5 * np.pi / 180
-    assert mis_angular.TD == 0
-    assert mis_angular.TL == 0
-    assert mis_angular.n1 == 0
-    assert mis_angular.speed == 125.66370614359172
-
-
-def test_mis_angular_parameters_units(mis_angular_units):
-    assert mis_angular_units.dt == 0.1
-    assert mis_angular_units.tI == 0
-    assert mis_angular_units.tF == 5
-    assert mis_angular_units.kd == 40 * 10 ** (3)
-    assert mis_angular_units.ks == 38 * 10 ** (3)
-    assert mis_angular_units.eCOUPx == 2 * 10 ** (-4)
-    assert mis_angular_units.eCOUPy == 2 * 10 ** (-4)
-    assert mis_angular_units.misalignment_angle == 5 * np.pi / 180
-    assert mis_angular_units.TD == 0
-    assert mis_angular_units.TL == 0
-    assert mis_angular_units.n1 == 0
-    assert mis_angular_units.speed == 125.66370614359172
-
-
-def test_mis_angular_forces(mis_angular):
-    assert mis_angular.forces[mis_angular.n1 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_angular.forces[mis_angular.n1 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [-2.66453526e-15,  1.66147096e-11,  3.32343042e-11,  4.98774355e-11,
-    6.64779343e-11,  8.30784330e-11,  9.97633087e-11,  1.16381571e-10,
-    1.32931444e-10,  1.49549262e-10,  1.66149317e-10,  1.82851956e-10,
-    1.99536387e-10,  2.16086704e-10,  2.32772024e-10,  2.49321896e-10,
-    2.65872213e-10,  2.82556645e-10,  2.99107850e-10,  3.15792725e-10,
-    3.32375905e-10,  3.48960416e-10,  3.65779851e-10,  3.82330612e-10,
-    3.99150490e-10,  4.15429913e-10,  4.32250680e-10,  4.48800552e-10,
-    4.65620431e-10,  4.82170748e-10,  4.98721064e-10,  5.15271381e-10,
-    5.31820366e-10,  5.48641577e-10,  5.65191449e-10,  5.81740878e-10,
-    5.98291194e-10,  6.15111961e-10,  6.31661390e-10,  6.48481713e-10,
-    6.64762023e-10,  6.81716905e-10,  6.97726765e-10,  7.14276194e-10,
-    7.31366523e-10,  7.47917284e-10,  7.64467156e-10,  7.81017029e-10,
-    7.98106914e-10,  8.14657675e-10,  8.30667091e-10
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_angular.forces[mis_angular.n2 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [-2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_angular.forces[mis_angular.n2 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [ 2.66453526e-15, -1.66147096e-11, -3.32343042e-11, -4.98774355e-11,
-    -6.64779343e-11, -8.30784330e-11, -9.97633087e-11, -1.16381571e-10,
-    -1.32931444e-10, -1.49549262e-10, -1.66149317e-10, -1.82851956e-10,
-    -1.99536387e-10, -2.16086704e-10, -2.32772024e-10, -2.49321896e-10,
-    -2.65872213e-10, -2.82556645e-10, -2.99107850e-10, -3.15792725e-10,
-    -3.32375905e-10, -3.48960416e-10, -3.65779851e-10, -3.82330612e-10,
-    -3.99150490e-10, -4.15429913e-10, -4.32250680e-10, -4.48800552e-10,
-    -4.65620431e-10, -4.82170748e-10, -4.98721064e-10, -5.15271381e-10,
-    -5.31820366e-10, -5.48641577e-10, -5.65191449e-10, -5.81740878e-10,
-    -5.98291194e-10, -6.15111961e-10, -6.31661390e-10, -6.48481713e-10,
-    -6.64762023e-10, -6.81716905e-10, -6.97726765e-10, -7.14276194e-10,
-    -7.31366523e-10, -7.47917284e-10, -7.64467156e-10, -7.81017029e-10,
-    -7.98106914e-10, -8.14657675e-10, -8.30667091e-10
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-
-def test_mis_angular_forces_units(mis_angular_units):
-    assert mis_angular_units.forces[mis_angular_units.n1 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782, 2.37707782, 2.37707782, 2.37707782, 2.37707782,
-    2.37707782
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_angular_units.forces[mis_angular_units.n1 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [-2.66453526e-15,  1.66147096e-11,  3.32343042e-11,  4.98774355e-11,
-    6.64779343e-11,  8.30784330e-11,  9.97633087e-11,  1.16381571e-10,
-    1.32931444e-10,  1.49549262e-10,  1.66149317e-10,  1.82851956e-10,
-    1.99536387e-10,  2.16086704e-10,  2.32772024e-10,  2.49321896e-10,
-    2.65872213e-10,  2.82556645e-10,  2.99107850e-10,  3.15792725e-10,
-    3.32375905e-10,  3.48960416e-10,  3.65779851e-10,  3.82330612e-10,
-    3.99150490e-10,  4.15429913e-10,  4.32250680e-10,  4.48800552e-10,
-    4.65620431e-10,  4.82170748e-10,  4.98721064e-10,  5.15271381e-10,
-    5.31820366e-10,  5.48641577e-10,  5.65191449e-10,  5.81740878e-10,
-    5.98291194e-10,  6.15111961e-10,  6.31661390e-10,  6.48481713e-10,
-    6.64762023e-10,  6.81716905e-10,  6.97726765e-10,  7.14276194e-10,
-    7.31366523e-10,  7.47917284e-10,  7.64467156e-10,  7.81017029e-10,
-    7.98106914e-10,  8.14657675e-10,  8.30667091e-10
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_angular_units.forces[mis_angular_units.n2 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [-2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782, -2.37707782, -2.37707782, -2.37707782, -2.37707782,
-    -2.37707782
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-    assert mis_angular_units.forces[mis_angular_units.n2 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [ 2.66453526e-15, -1.66147096e-11, -3.32343042e-11, -4.98774355e-11,
-    -6.64779343e-11, -8.30784330e-11, -9.97633087e-11, -1.16381571e-10,
-    -1.32931444e-10, -1.49549262e-10, -1.66149317e-10, -1.82851956e-10,
-    -1.99536387e-10, -2.16086704e-10, -2.32772024e-10, -2.49321896e-10,
-    -2.65872213e-10, -2.82556645e-10, -2.99107850e-10, -3.15792725e-10,
-    -3.32375905e-10, -3.48960416e-10, -3.65779851e-10, -3.82330612e-10,
-    -3.99150490e-10, -4.15429913e-10, -4.32250680e-10, -4.48800552e-10,
-    -4.65620431e-10, -4.82170748e-10, -4.98721064e-10, -5.15271381e-10,
-    -5.31820366e-10, -5.48641577e-10, -5.65191449e-10, -5.81740878e-10,
-    -5.98291194e-10, -6.15111961e-10, -6.31661390e-10, -6.48481713e-10,
-    -6.64762023e-10, -6.81716905e-10, -6.97726765e-10, -7.14276194e-10,
-    -7.31366523e-10, -7.47917284e-10, -7.64467156e-10, -7.81017029e-10,
-    -7.98106914e-10, -8.14657675e-10, -8.30667091e-10
-    ]
-    ),
-    abs=1e-6
-        # fmt: on
-    )
-
-
-@pytest.fixture
-def mis_rigid():
-    unbalance_magnitudet = np.array([5e-4, 0])
-    unbalance_phaset = np.array([-np.pi / 2, 0])
-
-    misalignment = rotor.run_misalignment(
-        coupling="rigid",
-        dt=0.0001,
-        tI=0,
-        tF=0.005,
-        eCOUP=2e-4,
-        TD=0,
-        TL=0,
-        n1=0,
-        speed=125.66370614359172,
-        unbalance_magnitude=unbalance_magnitudet,
-        unbalance_phase=unbalance_phaset,
-        print_progress=False,
-    )
-
-    return misalignment
-
-
-@pytest.fixture
-def mis_rigid_units():
-    unbalance_magnitudet = Q_(np.array([0.043398083107259365, 0]), "lb*in")
-    unbalance_phaset = Q_(np.array([-90.0, 0.0]), "degrees")
-
-    misalignment = rotor.run_misalignment(
+    return rotor.run_misalignment(
         coupling="rigid",
         dt=0.0001,
         tI=0,
@@ -868,8 +336,6 @@ def mis_rigid_units():
         print_progress=False,
     )
 
-    return misalignment
-
 
 def test_mis_rigid_parameters(mis_rigid):
     assert mis_rigid.dt == 0.0001
@@ -882,192 +348,71 @@ def test_mis_rigid_parameters(mis_rigid):
     assert mis_rigid.speed == 125.66370614359172
 
 
-def test_mis_rigid_parameters_units(mis_rigid_units):
-    assert mis_rigid_units.dt == 0.0001
-    assert mis_rigid_units.tI == 0
-    assert mis_rigid_units.tF == 0.005
-    assert mis_rigid_units.eCOUP == 2e-4
-    assert mis_rigid_units.TD == 0
-    assert mis_rigid_units.TL == 0
-    assert mis_rigid_units.n1 == 0
-    assert mis_rigid_units.speed == 125.66370614359172
-
-
 def test_mis_rigid_forces(mis_rigid):
-    assert mis_rigid.forces[mis_rigid.n1 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [0.00000000e+00, 4.36964689e+00, 1.74771689e+01, 3.93186457e+01,
-       6.98878473e+01, 1.09176370e+02, 1.57173806e+02, 2.13867938e+02,
-       2.79244943e+02, 3.53289607e+02, 4.35985520e+02, 5.27315254e+02,
-       6.27260511e+02, 7.35802235e+02, 8.52920676e+02, 9.78595415e+02,
-       1.11280534e+03, 1.25552860e+03, 1.40674247e+03, 1.56642326e+03,
-       1.73454619e+03, 1.91108521e+03, 2.09601285e+03, 2.28930016e+03,
-       2.49091658e+03, 2.70082986e+03, 2.91900610e+03, 3.14540972e+03,
-       3.38000359e+03, 3.62274912e+03, 3.87360643e+03, 4.13253452e+03,
-       4.39949148e+03, 4.67443471e+03, 4.95732112e+03, 5.24810728e+03,
-       5.54674961e+03, 5.85320444e+03, 6.16742807e+03, 6.48937676e+03,
-       6.81900661e+03, 7.15627347e+03, 7.50113266e+03, 7.85353879e+03,
-       8.21344538e+03, 8.58080462e+03, 8.95556694e+03, 9.33768077e+03,
-       9.72709218e+03, 1.01237446e+04, 1.05275788e+04
-    ]
-    ),
-    abs=1e-2
-        # fmt: on
-    )
-
-    assert mis_rigid.forces[mis_rigid.n1 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [     0.        ,   -695.44989191,  -1390.76098145,  -2085.80340207,
-        -2780.44896672,  -3474.57174737,  -4168.04855579,  -4860.75930785,
-        -5552.58725831,  -6243.41909941,  -6933.14492308,  -7621.65805293,
-        -8308.85475863,  -8994.63387028,  -9678.89631528, -10361.54460307,
-       -11042.48228557, -11721.61342139, -12398.84207127, -13074.07184973,
-       -13747.20555483, -14418.14489306, -15086.79031112, -15753.04094068,
-       -16416.79465562, -17077.94823571, -17736.39762468, -18392.03826576,
-       -19044.76549395, -19694.47496121, -20341.06306924, -20984.42738466,
-       -21624.46701205, -22261.08290347, -22894.17808613, -23523.65779508,
-       -24149.42950266, -24771.40284243, -25389.48943124, -26003.60259834,
-       -26613.65703604, -27219.56839029, -27821.25281318, -28418.62650121,
-       -29011.60524412, -29600.10400851, -30184.03657848, -30763.31527277,
-       -31337.85075352, -31907.55193714, -32472.32601231
-    ]
-    ),
-    abs=1e-2
-        # fmt: on
-    )
-
-    assert mis_rigid.forces[mis_rigid.n2 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [ 0.00000000e+00, -4.36964689e+00, -1.74771689e+01, -3.93186457e+01,
-       -6.98878473e+01, -1.09176370e+02, -1.57173806e+02, -2.13867938e+02,
-       -2.79244943e+02, -3.53289607e+02, -4.35985520e+02, -5.27315254e+02,
-       -6.27260511e+02, -7.35802235e+02, -8.52920676e+02, -9.78595415e+02,
-       -1.11280534e+03, -1.25552860e+03, -1.40674247e+03, -1.56642326e+03,
-       -1.73454619e+03, -1.91108521e+03, -2.09601285e+03, -2.28930016e+03,
-       -2.49091658e+03, -2.70082986e+03, -2.91900610e+03, -3.14540972e+03,
-       -3.38000359e+03, -3.62274912e+03, -3.87360643e+03, -4.13253452e+03,
-       -4.39949148e+03, -4.67443471e+03, -4.95732112e+03, -5.24810728e+03,
-       -5.54674961e+03, -5.85320444e+03, -6.16742807e+03, -6.48937676e+03,
-       -6.81900661e+03, -7.15627347e+03, -7.50113266e+03, -7.85353879e+03,
-       -8.21344538e+03, -8.58080462e+03, -8.95556694e+03, -9.33768077e+03,
-       -9.72709218e+03, -1.01237446e+04, -1.05275788e+04
-    ]
-    ),
-    abs=1e-2
-        # fmt: on
-    )
-
-    assert mis_rigid.forces[mis_rigid.n2 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [    0.        ,   695.44989191,  1390.76098145,  2085.80340207,
-        2780.44896672,  3474.57174737,  4168.04855579,  4860.75930785,
-        5552.58725831,  6243.41909941,  6933.14492308,  7621.65805293,
-        8308.85475863,  8994.63387028,  9678.89631528, 10361.54460307,
-       11042.48228557, 11721.61342139, 12398.84207127, 13074.07184973,
-       13747.20555483, 14418.14489306, 15086.79031112, 15753.04094068,
-       16416.79465562, 17077.94823571, 17736.39762468, 18392.03826576,
-       19044.76549395, 19694.47496121, 20341.06306924, 20984.42738466,
-       21624.46701205, 22261.08290347, 22894.17808613, 23523.65779508,
-       24149.42950266, 24771.40284243, 25389.48943124, 26003.60259834,
-       26613.65703604, 27219.56839029, 27821.25281318, 28418.62650121,
-       29011.60524412, 29600.10400851, 30184.03657848, 30763.31527277,
-       31337.85075352, 31907.55193714, 32472.32601231
-    ]
-    ),
-    abs=1e-2
-        # fmt: on
-    )
-
-
-def test_mis_rigid_forces_units(mis_rigid_units):
-    assert mis_rigid_units.forces[mis_rigid_units.n1 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [0.00000000e+00, 4.36964689e+00, 1.74771689e+01, 3.93186457e+01,
-       6.98878473e+01, 1.09176370e+02, 1.57173806e+02, 2.13867938e+02,
-       2.79244943e+02, 3.53289607e+02, 4.35985520e+02, 5.27315254e+02,
-       6.27260511e+02, 7.35802235e+02, 8.52920676e+02, 9.78595415e+02,
-       1.11280534e+03, 1.25552860e+03, 1.40674247e+03, 1.56642326e+03,
-       1.73454619e+03, 1.91108521e+03, 2.09601285e+03, 2.28930016e+03,
-       2.49091658e+03, 2.70082986e+03, 2.91900610e+03, 3.14540972e+03,
-       3.38000359e+03, 3.62274912e+03, 3.87360643e+03, 4.13253452e+03,
-       4.39949148e+03, 4.67443471e+03, 4.95732112e+03, 5.24810728e+03,
-       5.54674961e+03, 5.85320444e+03, 6.16742807e+03, 6.48937676e+03,
-       6.81900661e+03, 7.15627347e+03, 7.50113266e+03, 7.85353879e+03,
-       8.21344538e+03, 8.58080462e+03, 8.95556694e+03, 9.33768077e+03,
-       9.72709218e+03, 1.01237446e+04, 1.05275788e+04
-    ]
-    ),
-    abs=1e-2
-        # fmt: on
-    )
-
-    assert mis_rigid_units.forces[mis_rigid_units.n1 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [     0.        ,   -695.44989191,  -1390.76098145,  -2085.80340207,
-        -2780.44896672,  -3474.57174737,  -4168.04855579,  -4860.75930785,
-        -5552.58725831,  -6243.41909941,  -6933.14492308,  -7621.65805293,
-        -8308.85475863,  -8994.63387028,  -9678.89631528, -10361.54460307,
-       -11042.48228557, -11721.61342139, -12398.84207127, -13074.07184973,
-       -13747.20555483, -14418.14489306, -15086.79031112, -15753.04094068,
-       -16416.79465562, -17077.94823571, -17736.39762468, -18392.03826576,
-       -19044.76549395, -19694.47496121, -20341.06306924, -20984.42738466,
-       -21624.46701205, -22261.08290347, -22894.17808613, -23523.65779508,
-       -24149.42950266, -24771.40284243, -25389.48943124, -26003.60259834,
-       -26613.65703604, -27219.56839029, -27821.25281318, -28418.62650121,
-       -29011.60524412, -29600.10400851, -30184.03657848, -30763.31527277,
-       -31337.85075352, -31907.55193714, -32472.32601231
-    ]
-    ),
-    abs=1e-2
-        # fmt: on
-    )
-
-    assert mis_rigid_units.forces[mis_rigid_units.n2 * 6, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [ 0.00000000e+00, -4.36964689e+00, -1.74771689e+01, -3.93186457e+01,
-       -6.98878473e+01, -1.09176370e+02, -1.57173806e+02, -2.13867938e+02,
-       -2.79244943e+02, -3.53289607e+02, -4.35985520e+02, -5.27315254e+02,
-       -6.27260511e+02, -7.35802235e+02, -8.52920676e+02, -9.78595415e+02,
-       -1.11280534e+03, -1.25552860e+03, -1.40674247e+03, -1.56642326e+03,
-       -1.73454619e+03, -1.91108521e+03, -2.09601285e+03, -2.28930016e+03,
-       -2.49091658e+03, -2.70082986e+03, -2.91900610e+03, -3.14540972e+03,
-       -3.38000359e+03, -3.62274912e+03, -3.87360643e+03, -4.13253452e+03,
-       -4.39949148e+03, -4.67443471e+03, -4.95732112e+03, -5.24810728e+03,
-       -5.54674961e+03, -5.85320444e+03, -6.16742807e+03, -6.48937676e+03,
-       -6.81900661e+03, -7.15627347e+03, -7.50113266e+03, -7.85353879e+03,
-       -8.21344538e+03, -8.58080462e+03, -8.95556694e+03, -9.33768077e+03,
-       -9.72709218e+03, -1.01237446e+04, -1.05275788e+04
-    ]
-    ),
-    abs=1e-2
-        # fmt: on
-    )
-
-    assert mis_rigid_units.forces[mis_rigid_units.n2 * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array(
-    [    0.        ,   695.44989191,  1390.76098145,  2085.80340207,
-        2780.44896672,  3474.57174737,  4168.04855579,  4860.75930785,
-        5552.58725831,  6243.41909941,  6933.14492308,  7621.65805293,
-        8308.85475863,  8994.63387028,  9678.89631528, 10361.54460307,
-       11042.48228557, 11721.61342139, 12398.84207127, 13074.07184973,
-       13747.20555483, 14418.14489306, 15086.79031112, 15753.04094068,
-       16416.79465562, 17077.94823571, 17736.39762468, 18392.03826576,
-       19044.76549395, 19694.47496121, 20341.06306924, 20984.42738466,
-       21624.46701205, 22261.08290347, 22894.17808613, 23523.65779508,
-       24149.42950266, 24771.40284243, 25389.48943124, 26003.60259834,
-       26613.65703604, 27219.56839029, 27821.25281318, 28418.62650121,
-       29011.60524412, 29600.10400851, 30184.03657848, 30763.31527277,
-       31337.85075352, 31907.55193714, 32472.32601231
-    ]
-    ),
-    abs=1e-2
-        # fmt: on
-    )
+    # fmt: off
+    Fx_n1_rig = np.array([
+            0.     ,     4.31682,    17.26592,    38.84351,    69.04359,
+          107.85803,   155.27676,   211.28794,   275.87817,   349.03266,
+          430.73544,   520.9695 ,   619.71696,   726.95916,   842.67672,
+          966.84955,  1099.45687,  1240.47709,  1389.8878 ,  1547.66558,
+         1713.78589,  1888.22299,  2070.94971,  2261.93744,  2461.15597,
+         2668.57346,  2884.15645,  3107.86985,  3339.67703,  3579.53994,
+         3827.41925,  4083.27451,  4347.06438,  4618.74677,  4898.27908,
+         5185.61834,  5480.72137,  5783.54484,  6094.04533,  6412.17932,
+         6737.90309,  7071.1726 ,  7411.94329,  7760.16982,  8115.80585,
+         8478.80363,  8849.11381,  9226.68502,  9611.46366, 10003.39363,
+        10402.41613
+    ])
+    
+    Fy_n1_rig = np.array([
+             0.     ,   -687.04194,  -1373.94823,  -2060.59157,
+         -2746.84626,  -3432.58871,  -4117.6979 ,  -4802.05574,
+         -5485.54725,  -6168.06069,  -6849.48752,  -7529.72223,
+         -8208.66208,  -8886.20674,  -9562.25785, -10236.71857,
+        -10909.49302, -11580.48583, -12249.60164, -12916.74468,
+        -13581.81844, -14244.72539, -14905.36686, -15563.64295,
+        -16219.4526 , -16872.69375, -17523.26358, -18171.05879,
+        -18815.97606, -19457.91233, -20096.7653 , -20732.43379,
+        -21364.81807, -21993.82016, -22619.34409, -23241.29595,
+        -23859.58398, -24474.11847, -25084.81159, -25691.57716,
+        -26294.33031, -26892.98709, -27487.46407, -28077.67788,
+        -28663.54481, -29244.98043, -29821.89922, -30394.21433,
+        -30961.83738, -31524.67842, -32082.6459
+    ])
+
+    Fx_n2_rig = np.array([
+            0.     ,     -4.31682,    -17.26592,    -38.84351,
+          -69.04359,   -107.85803,   -155.27676,   -211.28794,
+         -275.87817,   -349.03266,   -430.73544,   -520.9695 ,
+         -619.71696,   -726.95916,   -842.67672,   -966.84955,
+        -1099.45687,  -1240.47709,  -1389.8878 ,  -1547.66558,
+        -1713.78589,  -1888.22299,  -2070.94971,  -2261.93744,
+        -2461.15597,  -2668.57346,  -2884.15645,  -3107.86985,
+        -3339.67703,  -3579.53994,  -3827.41925,  -4083.27451,
+        -4347.06438,  -4618.74677,  -4898.27908,  -5185.61834,
+        -5480.72137,  -5783.54484,  -6094.04533,  -6412.17932,
+        -6737.90309,  -7071.1726 ,  -7411.94329,  -7760.16982,
+        -8115.80585,  -8478.80363,  -8849.11381,  -9226.68502,
+        -9611.46366, -10003.39363, -10402.41613
+    ])
+
+    Fy_n2_rig = np.array([
+            0.     ,   687.04194,  1373.94823,  2060.59157,  2746.84626,
+         3432.58871,  4117.6979 ,  4802.05574,  5485.54725,  6168.06069,
+         6849.48752,  7529.72223,  8208.66208,  8886.20674,  9562.25785,
+        10236.71857, 10909.49302, 11580.48583, 12249.60164, 12916.74468,
+        13581.81844, 14244.72539, 14905.36686, 15563.64295, 16219.4526 ,
+        16872.69375, 17523.26358, 18171.05879, 18815.97606, 19457.91233,
+        20096.7653 , 20732.43379, 21364.81807, 21993.82016, 22619.34409,
+        23241.29595, 23859.58398, 24474.11847, 25084.81159, 25691.57716,
+        26294.33031, 26892.98709, 27487.46407, 28077.67788, 28663.54481,
+        29244.98043, 29821.89922, 30394.21433, 30961.83738, 31524.67842,
+        32082.6459
+    ])
+    # fmt: on
+
+    assert_allclose(mis_rigid.forces[mis_rigid.n1 * 6 + 0, :], Fx_n1_rig, rtol=3e-2)
+    assert_allclose(mis_rigid.forces[mis_rigid.n1 * 6 + 1, :], Fy_n1_rig, rtol=3e-2)
+
+    assert_allclose(mis_rigid.forces[mis_rigid.n2 * 6 + 0, :], Fx_n2_rig, rtol=3e-2)
+    assert_allclose(mis_rigid.forces[mis_rigid.n2 * 6 + 1, :], Fy_n2_rig, rtol=3e-2)
diff --git a/ross/tests/test_rotor_assembly.py b/ross/tests/test_rotor_assembly.py
index c9f8380e2..9ef182f40 100644
--- a/ross/tests/test_rotor_assembly.py
+++ b/ross/tests/test_rotor_assembly.py
@@ -1733,16 +1733,12 @@ def rotor_6dof():
 
     shaft_elem = [
         ShaftElement6DoF(
+            l,
+            i_d,
+            o_d,
             material=steel,
-            L=0.25,
-            idl=0,
-            odl=0.05,
-            idr=0,
-            odr=0.05,
-            alpha=0,
-            beta=0,
-            rotary_inertia=False,
-            shear_effects=False,
+            alpha=1,
+            beta=1e-5,
         )
         for l in L
     ]
@@ -1754,44 +1750,16 @@ def rotor_6dof():
         n=4, material=steel, width=0.07, i_d=0.05, o_d=0.28
     )
 
-    kxx = 1e6
-    kyy = 0.8e6
-    kzz = 1e5
-    cxx = 0
-    cyy = 0
-    czz = 0
-    bearing0 = BearingElement6DoF(
-        n=0, kxx=kxx, kyy=kyy, cxx=cxx, cyy=cyy, kzz=kzz, czz=czz
-    )
-    bearing1 = BearingElement6DoF(
-        n=6, kxx=kxx, kyy=kyy, cxx=cxx, cyy=cyy, kzz=kzz, czz=czz
-    )
+    bearing0 = BearingElement6DoF(n=0, kxx=1e6, kyy=8e5, kzz=1e5, cxx=0, cyy=0, czz=0)
+    bearing1 = BearingElement6DoF(n=6, kxx=1e6, kyy=8e5, kzz=1e5, cxx=0, cyy=0, czz=0)
 
     return Rotor(shaft_elem, [disk0, disk1], [bearing0, bearing1])
 
 
 def test_modal_6dof(rotor_6dof):
     modal = rotor_6dof.run_modal(speed=0, sparse=False)
-    wn = np.array(
-        [
-            9.91427121e-05,
-            4.76215566e01,
-            9.17987032e01,
-            9.62914807e01,
-            2.74579882e02,
-            2.96518344e02,
-        ]
-    )
-    wd = np.array(
-        [
-            9.91427121e-05,
-            4.76215566e01,
-            9.17987032e01,
-            9.62914807e01,
-            2.74579882e02,
-            2.96518344e02,
-        ]
-    )
+    wn = np.array([0.0, 47.62138, 91.79647, 96.28891, 274.56591, 296.5005])
+    wd = np.array([0.01079, 47.62156, 91.79656, 96.289, 274.56607, 296.50068])
 
     assert_almost_equal(modal.wn[:6], wn, decimal=2)
     assert_almost_equal(modal.wd[:6], wd, decimal=2)
diff --git a/ross/tests/test_rubbing.py b/ross/tests/test_rubbing.py
index 27e9c5bea..d539a27db 100644
--- a/ross/tests/test_rubbing.py
+++ b/ross/tests/test_rubbing.py
@@ -1,105 +1,75 @@
-import os
-from pathlib import Path
-from tempfile import tempdir
-
 import numpy as np
 import pytest
-from numpy.testing import assert_allclose, assert_almost_equal
+from numpy.testing import assert_allclose
 
 import ross as rs
-from ross.defects.misalignment import MisalignmentFlex
 from ross.units import Q_
 
-steel2 = rs.Material(name="Steel", rho=7850, E=2.17e11, Poisson=0.2992610837438423)
 
-#  Rotor with 6 DoFs, with internal damping, with 10 shaft elements, 2 disks and 2 bearings.
-i_d = 0
-o_d = 0.019
-n = 33
+@pytest.fixture
+def rotor():
+    steel2 = rs.Material(name="Steel", rho=7850, E=2.17e11, Poisson=0.2992610837438423)
+
+    #  Rotor with 6 DoFs, with internal damping, with 10 shaft elements, 2 disks and 2 bearings.
+    i_d = 0
+    o_d = 0.019
 
-# fmt: off
-L = np.array(
-        [0  ,  25,  64, 104, 124, 143, 175, 207, 239, 271,
+    # fmt: off
+    L = np.array([
+        0  ,  25,  64, 104, 124, 143, 175, 207, 239, 271,
         303, 335, 345, 355, 380, 408, 436, 466, 496, 526,
         556, 586, 614, 647, 657, 667, 702, 737, 772, 807,
-        842, 862, 881, 914]
-        )/ 1000
-# fmt: on
-
-L = [L[i] - L[i - 1] for i in range(1, len(L))]
-
-shaft_elem = [
-    rs.ShaftElement6DoF(
-        material=steel2,
-        L=l,
-        idl=i_d,
-        odl=o_d,
-        idr=i_d,
-        odr=o_d,
-        alpha=8.0501,
-        beta=1.0e-5,
-        rotary_inertia=True,
-        shear_effects=True,
-    )
-    for l in L
-]
-
-Id = 0.003844540885417
-Ip = 0.007513248437500
+        842, 862, 881, 914
+    ])/ 1000
+    # fmt: on
 
-disk0 = rs.DiskElement6DoF(n=12, m=2.6375, Id=Id, Ip=Ip)
-disk1 = rs.DiskElement6DoF(n=24, m=2.6375, Id=Id, Ip=Ip)
+    L = [L[i] - L[i - 1] for i in range(1, len(L))]
 
-kxx1 = 4.40e5
-kyy1 = 4.6114e5
-kzz = 0
-cxx1 = 27.4
-cyy1 = 2.505
-czz = 0
-kxx2 = 9.50e5
-kyy2 = 1.09e8
-cxx2 = 50.4
-cyy2 = 100.4553
+    shaft_elem = [
+        rs.ShaftElement6DoF(
+            l,
+            i_d,
+            o_d,
+            material=steel2,
+            alpha=8.0501,
+            beta=1.0e-5,
+        )
+        for l in L
+    ]
 
-bearing0 = rs.BearingElement6DoF(
-    n=4, kxx=kxx1, kyy=kyy1, cxx=cxx1, cyy=cyy1, kzz=kzz, czz=czz
-)
-bearing1 = rs.BearingElement6DoF(
-    n=31, kxx=kxx2, kyy=kyy2, cxx=cxx2, cyy=cyy2, kzz=kzz, czz=czz
-)
+    Id = 0.003844540885417
+    Ip = 0.007513248437500
 
-rotor = rs.Rotor(shaft_elem, [disk0, disk1], [bearing0, bearing1])
+    disk0 = rs.DiskElement6DoF(n=12, m=2.6375, Id=Id, Ip=Ip)
+    disk1 = rs.DiskElement6DoF(n=24, m=2.6375, Id=Id, Ip=Ip)
 
+    kxx1 = 4.40e5
+    kyy1 = 4.6114e5
+    kzz = 0
+    cxx1 = 27.4
+    cyy1 = 2.505
+    czz = 0
+    kxx2 = 9.50e5
+    kyy2 = 1.09e8
+    cxx2 = 50.4
+    cyy2 = 100.4553
 
-@pytest.fixture
-def rub():
-    unbalance_magnitudet = np.array([5e-4, 0])
-    unbalance_phaset = np.array([-np.pi / 2, 0])
-
-    rubbing = rotor.run_rubbing(
-        dt=0.001,
-        tI=0,
-        tF=0.5,
-        deltaRUB=7.95e-5,
-        kRUB=1.1e6,
-        cRUB=40,
-        miRUB=0.3,
-        posRUB=12,
-        speed=125.66370614359172,
-        unbalance_magnitude=unbalance_magnitudet,
-        unbalance_phase=unbalance_phaset,
-        print_progress=True,
+    bearing0 = rs.BearingElement6DoF(
+        n=4, kxx=kxx1, kyy=kyy1, kzz=kzz, cxx=cxx1, cyy=cyy1, czz=czz
+    )
+    bearing1 = rs.BearingElement6DoF(
+        n=31, kxx=kxx2, kyy=kyy2, kzz=kzz, cxx=cxx2, cyy=cyy2, czz=czz
     )
 
-    return rubbing
+    return rs.Rotor(shaft_elem, [disk0, disk1], [bearing0, bearing1])
 
 
 @pytest.fixture
-def rub_units():
-    unbalance_magnitudet = Q_(np.array([0.043398083107259365, 0]), "lb*in")
-    unbalance_phaset = Q_(np.array([-90.0, 0.0]), "degrees")
+def rub(rotor):
+    unbalance_magnitude = np.array([5e-4, 0])
+    unbalance_phase = Q_(np.array([-90.0, 0.0]), "degrees")
 
-    rubbing = rotor.run_rubbing(
+    return rotor.run_rubbing(
         dt=0.001,
         tI=0,
         tF=0.5,
@@ -109,13 +79,11 @@ def rub_units():
         miRUB=0.3,
         posRUB=12,
         speed=Q_(1200, "RPM"),
-        unbalance_magnitude=unbalance_magnitudet,
-        unbalance_phase=unbalance_phaset,
-        print_progress=True,
+        unbalance_magnitude=unbalance_magnitude,
+        unbalance_phase=unbalance_phase,
+        print_progress=False,
     )
 
-    return rubbing
-
 
 def test_rub_parameters(rub):
     assert rub.dt == 0.001
@@ -129,545 +97,182 @@ def test_rub_parameters(rub):
     assert rub.speed == 125.66370614359172
 
 
-def test_rub_parameters_units(rub_units):
-    assert rub_units.dt == 0.001
-    assert rub_units.tI == 0
-    assert rub_units.tF == 0.5
-    assert rub_units.deltaRUB == 7.95e-5
-    assert rub_units.kRUB == 1.1e6
-    assert rub_units.cRUB == 40
-    assert rub_units.miRUB == 0.3
-    assert rub_units.posRUB == 12
-    assert rub_units.speed == 125.66370614359172
-
-
 def test_rub_forces(rub):
-    assert rub.forces_rub[rub.posRUB * 6, :] == pytest.approx(
-        # fmt: off
-    np.array([ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        1.34106430e+00,  2.38612745e+00,  2.49817920e+00,  1.81347173e+00,
-        5.98445473e-01, -1.62547004e+00, -4.23917250e+00, -6.00118806e+00,
-       -6.20373067e+00, -4.73612273e+00, -2.73184699e+00, -2.58971885e+00,
-       -6.09874812e+00, -1.19048343e+01, -1.62930120e+01, -1.64120039e+01,
-       -1.22951018e+01, -6.47561269e+00, -1.92419235e+00,  3.53374623e-02,
-        1.47810258e-01,  6.34772882e-02,  1.83806641e-01,  4.43954908e-01,
-        5.92533172e-01,  5.51738633e-01,  6.18865825e-01,  1.16291085e+00,
-        1.99527847e+00,  2.28657481e+00,  1.21810729e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00, -1.55680614e-01,
-       -1.15087399e+00, -1.31570879e+00, -5.73686616e-01,  1.75075759e-01,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        1.35857558e+00,  4.63547122e+00,  8.26625410e+00,  1.07680988e+01,
-        1.13286366e+01,  9.96961037e+00,  7.44128958e+00,  4.81928761e+00,
-        2.91576005e+00,  1.93968467e+00,  1.61560462e+00,  1.51666598e+00,
-        1.30530274e+00,  8.25200469e-01,  1.14260074e-01, -1.19244337e+00,
-       -2.29090277e+00, -2.84242534e+00, -2.70568998e+00, -1.71797718e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-       -3.62954221e-01, -2.15664186e-02, -5.82573814e-01, -1.27145991e+00,
-       -1.38207772e+00, -8.36333007e-01,  0.00000000e+00, -1.52638197e+00,
-       -3.88459875e+00, -5.84147533e+00, -5.55533169e+00, -2.84156510e+00,
-        0.00000000e+00,  0.00000000e+00, -9.98460245e-01, -5.03924706e+00,
-       -7.75632048e+00, -7.28589841e+00, -4.21050757e+00, -6.78335856e-01,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        4.97669123e-01,  9.43108661e-01,  1.21217301e+00,  1.02059004e+00,
-        5.88846364e-01,  4.25532247e-01,  7.49102363e-01,  1.21248611e+00,
-        1.23719949e+00,  5.26714086e-01,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  8.30579340e-01,  2.26273248e+00,
-        2.80172735e+00,  2.36060364e+00,  1.47121281e+00,  9.05484679e-01,
-        1.10607178e+00,  1.80305085e+00,  2.28352354e+00,  2.03304519e+00,
-        1.10298146e+00,  3.77694650e-02,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00, -3.22380450e+00,
-       -5.68021237e+00, -5.83462167e+00, -3.25195544e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00, -9.90196556e-01,
-       -2.90885128e+00, -2.81169132e+00, -8.81117993e-01,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  2.69024403e-01,  3.11412852e-01,  5.58798236e-01,
-        6.40371175e-01,  3.49164726e-01,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  2.36971006e-01,
-        1.21346945e+00,  2.18249083e+00,  2.84208082e+00,  2.99156571e+00,
-        2.55421246e+00,  1.59111859e+00,  3.14246711e-01,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-       -9.16343507e-02,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00, -1.16296507e+00, -3.23286551e+00,
-       -4.51354090e+00, -3.98257895e+00, -1.81166144e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  2.12219345e-01,  2.21171647e-01,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00, -9.56292801e-01, -1.33271793e+00,
-       -1.03821500e-02,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00]),
-        abs=1e-6
-        # fmt: on
-    )
-
-    assert rub.forces_rub[rub.posRUB * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array([ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-       -2.04165157e+00, -4.37810280e+00, -6.57666995e+00, -8.28694141e+00,
-       -9.37791126e+00, -9.76081035e+00, -9.25783043e+00, -7.68973505e+00,
-       -5.19539491e+00, -2.48084518e+00, -6.09414243e-01, -1.40464287e-01,
-       -2.91868675e-01,  5.94864394e-01,  3.76903307e+00,  7.17867399e+00,
-        8.50615176e+00,  6.92755001e+00,  3.76919148e+00,  1.38886654e+00,
-        1.14793265e+00,  2.47476310e+00,  3.69518317e+00,  3.56151341e+00,
-        2.15845164e+00,  7.22745269e-01,  4.64776971e-01,  1.32361970e+00,
-        2.08542187e+00,  1.74043906e+00,  3.68327123e-01,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  1.40526270e+00,
-        3.91183574e+00,  6.08879987e+00,  6.05701802e+00,  3.35127909e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        1.33684009e+00,  4.74117308e+00,  5.66393572e+00,  3.75144548e+00,
-        2.83459537e-01, -1.61445732e+00, -2.61883874e+00, -2.76818696e+00,
-       -2.47679180e+00, -2.29014373e+00, -2.50794047e+00, -3.08754938e+00,
-       -3.77455164e+00, -4.26560645e+00, -4.33448168e+00, -3.93383052e+00,
-       -3.19736872e+00, -2.30738448e+00, -1.36778475e+00, -4.13019001e-01,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-       -2.30731635e-01, -2.71055349e+00, -3.77777456e+00, -3.19820518e+00,
-       -1.62052825e+00, -1.57786580e-01,  0.00000000e+00, -3.39194075e-01,
-       -1.19486567e+00, -1.40982868e+00, -8.24731765e-01,  1.77245380e-02,
-        0.00000000e+00,  0.00000000e+00,  1.12393272e+00,  3.42883960e+00,
-        6.44392419e+00,  7.74026022e+00,  6.02300041e+00,  2.15021680e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        1.31206244e+00,  2.88224433e+00,  2.70442915e+00,  1.31540273e+00,
-        6.83579157e-02, -1.44657955e-01, -9.74069771e-04,  2.01609711e-01,
-        7.58795510e-02, -1.14712672e-01,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00, -1.16918168e-01, -9.48748823e-02,
-       -3.74701248e-01, -6.68192840e-01, -7.77720853e-01, -8.15718667e-01,
-       -1.19404470e+00, -2.17926682e+00, -3.40215785e+00, -3.99418641e+00,
-       -3.25449483e+00, -1.22032570e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00, -1.72304245e+00,
-       -2.63776333e+00, -2.11110151e+00, -7.22241992e-01,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  1.37988828e-01,
-       -1.36497601e-01,  1.32278945e-01,  5.15571508e-01,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  8.42548433e-01,  2.72091640e+00,  3.24542027e+00,
-        2.03399132e+00,  3.77699186e-02,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00, -1.99799699e-01,
-        5.92260640e-02, -1.02272668e-02, -3.19939917e-01, -7.49813663e-01,
-       -1.06945353e+00, -1.04483902e+00, -5.36963585e-01,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-       -3.23192576e-01,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  1.42596962e-01, -1.71901114e-01,
-       -7.45893810e-02,  3.71940362e-01,  6.24286373e-01,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00, -9.16154067e-02, -8.72890060e-02,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  8.36067934e-01,  1.13854180e+00,
-        4.79721529e-01,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00]),
-        abs=1e-6
-        # fmt: on
-    )
+    # fmt: off
+    Fx_rub = np.array([
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         1.34377,   2.38806,   2.49916,   1.81395,   0.599  ,  -1.62394,
+        -4.23779,  -6.0012 ,  -6.20606,  -4.74044,  -2.73606,  -2.59117,
+        -6.09753, -11.90441, -16.29645, -16.419  , -12.30249,  -6.48021,
+        -1.92542,   0.03546,   0.14775,   0.06315,   0.18367,   0.44434,
+         0.59336,   0.55256,   0.61928,   1.1632 ,   1.99649,   2.28978,
+         1.22352,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,  -0.15535,  -1.14814,  -1.31374,
+        -0.57377,   0.17469,   0.     ,   0.     ,   0.     ,   0.     ,
+         1.36008,   4.63287,   8.26195,  10.76651,  11.33296,   9.97955,
+         7.45322,   4.82873,   2.92041,   1.94037,   1.61479,   1.51651,
+         1.30659,   0.82727,   0.11565,  -1.19318,  -2.29413,  -2.84508,
+        -2.70601,  -1.71827,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,  -0.36415,  -0.02201,  -0.58344,  -1.2737 ,
+        -1.38591,  -0.83953,   0.     ,  -1.51706,  -3.87348,  -5.83944,
+        -5.56822,  -2.86468,   0.     ,   0.     ,  -0.99869,  -5.03332,
+        -7.7501 ,  -7.2838 ,  -4.21365,  -0.68445,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.49957,   0.94559,   1.21312,   1.01922,
+         0.58609,   0.42304,   0.74788,   1.21229,   1.23721,   0.52635,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.8426 ,   2.27205,
+         2.8058 ,   2.35805,   1.46288,   0.89457,   1.09691,   1.79904,
+         2.28564,   2.0398 ,   1.11131,   0.04426,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,  -3.22059,  -5.67507,  -5.82974,
+        -3.24989,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,  -1.01539,
+        -2.92212,  -2.80863,  -0.86681,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.2655 ,   0.31103,   0.5602 ,   0.63978,   0.34391,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.21579,
+         1.19914,   2.17879,   2.84879,   3.0046 ,   2.56763,   1.59991,
+         0.3164 ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+        -0.05146,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,  -1.06699,  -3.18183,  -4.53292,  -4.05864,
+        -1.90024,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.23001,   0.23853,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.     ,   0.     ,  -0.91444,  -1.36307,
+        -0.09154,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ,
+         0.     ,   0.     ,   0.
+    ])
 
+    Fy_rub = np.array([
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+       -2.04429, -4.38016, -6.57714, -8.28573, -9.37587, -9.75933,
+       -9.25813, -7.69222, -5.19938, -2.48476, -0.61176, -0.14117,
+       -0.29231,  0.59298,  3.76722,  7.17877,  8.50843,  6.9303 ,
+        3.77043,  1.38873,  1.14884,  2.47906,  3.70296,  3.57058,
+        2.16593,  0.72693,  0.46617,  1.3243 ,  2.08718,  1.74357,
+        0.372  ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  1.40459,  3.90508,  6.08173,
+        6.05774,  3.36298,  0.     ,  0.     ,  0.     ,  0.     ,
+        1.34137,  4.74269,  5.66655,  3.75628,  0.28809, -1.61396,
+       -2.62114, -2.77201, -2.47978, -2.29098, -2.50753, -3.08839,
+       -3.77908, -4.27479, -4.34706, -3.94654, -3.20655, -2.31137,
+       -1.36806, -0.4128 ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     , -0.22816, -2.7129 , -3.78547, -3.20913,
+       -1.63064, -0.16281,  0.     , -0.33417, -1.19099, -1.40999,
+       -0.82791,  0.01215,  0.     ,  0.     ,  1.12437,  3.42666,
+        6.4417 ,  7.74143,  6.02904,  2.16009,  0.     ,  0.     ,
+        0.     ,  0.     ,  1.33055,  2.8963 ,  2.70812,  1.3118 ,
+        0.0632 , -0.1465 , -0.0018 ,  0.20072,  0.07503, -0.11498,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     , -0.1147 , -0.09399,
+       -0.37461, -0.66764, -0.77524, -0.81072, -1.18813, -2.17631,
+       -3.40636, -4.0075 , -3.27515, -1.24315,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     , -1.72064, -2.63464, -2.10851,
+       -0.72111,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.12732,
+       -0.13905,  0.13201,  0.51607,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.88088,  2.75998,  3.26768,  2.03091,  0.01545,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     , -0.20874,
+        0.05031, -0.0122 , -0.32162, -0.75307, -1.07413, -1.04887,
+       -0.53796,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+       -0.35022,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.18449, -0.16341, -0.07524,  0.3703 ,
+        0.6314 ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     , -0.06506, -0.06606,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.     ,  0.     ,  0.82617,  1.15647,
+        0.53299,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,
+        0.     ,  0.     ,  0.
+    ])
+    # fmt: on
 
-def test_rub_forces_units(rub_units):
-    assert rub_units.forces_rub[rub_units.posRUB * 6, :] == pytest.approx(
-        # fmt: off
-    np.array([ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        1.34106430e+00,  2.38612745e+00,  2.49817920e+00,  1.81347173e+00,
-        5.98445473e-01, -1.62547004e+00, -4.23917250e+00, -6.00118806e+00,
-       -6.20373067e+00, -4.73612273e+00, -2.73184699e+00, -2.58971885e+00,
-       -6.09874812e+00, -1.19048343e+01, -1.62930120e+01, -1.64120039e+01,
-       -1.22951018e+01, -6.47561269e+00, -1.92419235e+00,  3.53374623e-02,
-        1.47810258e-01,  6.34772882e-02,  1.83806641e-01,  4.43954908e-01,
-        5.92533172e-01,  5.51738633e-01,  6.18865825e-01,  1.16291085e+00,
-        1.99527847e+00,  2.28657481e+00,  1.21810729e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00, -1.55680614e-01,
-       -1.15087399e+00, -1.31570879e+00, -5.73686616e-01,  1.75075759e-01,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        1.35857558e+00,  4.63547122e+00,  8.26625410e+00,  1.07680988e+01,
-        1.13286366e+01,  9.96961037e+00,  7.44128958e+00,  4.81928761e+00,
-        2.91576005e+00,  1.93968467e+00,  1.61560462e+00,  1.51666598e+00,
-        1.30530274e+00,  8.25200469e-01,  1.14260074e-01, -1.19244337e+00,
-       -2.29090277e+00, -2.84242534e+00, -2.70568998e+00, -1.71797718e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-       -3.62954221e-01, -2.15664186e-02, -5.82573814e-01, -1.27145991e+00,
-       -1.38207772e+00, -8.36333007e-01,  0.00000000e+00, -1.52638197e+00,
-       -3.88459875e+00, -5.84147533e+00, -5.55533169e+00, -2.84156510e+00,
-        0.00000000e+00,  0.00000000e+00, -9.98460245e-01, -5.03924706e+00,
-       -7.75632048e+00, -7.28589841e+00, -4.21050757e+00, -6.78335856e-01,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        4.97669123e-01,  9.43108661e-01,  1.21217301e+00,  1.02059004e+00,
-        5.88846364e-01,  4.25532247e-01,  7.49102363e-01,  1.21248611e+00,
-        1.23719949e+00,  5.26714086e-01,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  8.30579340e-01,  2.26273248e+00,
-        2.80172735e+00,  2.36060364e+00,  1.47121281e+00,  9.05484679e-01,
-        1.10607178e+00,  1.80305085e+00,  2.28352354e+00,  2.03304519e+00,
-        1.10298146e+00,  3.77694650e-02,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00, -3.22380450e+00,
-       -5.68021237e+00, -5.83462167e+00, -3.25195544e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00, -9.90196556e-01,
-       -2.90885128e+00, -2.81169132e+00, -8.81117993e-01,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  2.69024403e-01,  3.11412852e-01,  5.58798236e-01,
-        6.40371175e-01,  3.49164726e-01,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  2.36971006e-01,
-        1.21346945e+00,  2.18249083e+00,  2.84208082e+00,  2.99156571e+00,
-        2.55421246e+00,  1.59111859e+00,  3.14246711e-01,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-       -9.16343507e-02,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00, -1.16296507e+00, -3.23286551e+00,
-       -4.51354090e+00, -3.98257895e+00, -1.81166144e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  2.12219345e-01,  2.21171647e-01,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00, -9.56292801e-01, -1.33271793e+00,
-       -1.03821500e-02,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00]),
-        abs=1e-6
-        # fmt: on
-    )
-
-    assert rub_units.forces_rub[rub_units.posRUB * 6 + 1, :] == pytest.approx(
-        # fmt: off
-    np.array([ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-       -2.04165157e+00, -4.37810280e+00, -6.57666995e+00, -8.28694141e+00,
-       -9.37791126e+00, -9.76081035e+00, -9.25783043e+00, -7.68973505e+00,
-       -5.19539491e+00, -2.48084518e+00, -6.09414243e-01, -1.40464287e-01,
-       -2.91868675e-01,  5.94864394e-01,  3.76903307e+00,  7.17867399e+00,
-        8.50615176e+00,  6.92755001e+00,  3.76919148e+00,  1.38886654e+00,
-        1.14793265e+00,  2.47476310e+00,  3.69518317e+00,  3.56151341e+00,
-        2.15845164e+00,  7.22745269e-01,  4.64776971e-01,  1.32361970e+00,
-        2.08542187e+00,  1.74043906e+00,  3.68327123e-01,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  1.40526270e+00,
-        3.91183574e+00,  6.08879987e+00,  6.05701802e+00,  3.35127909e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        1.33684009e+00,  4.74117308e+00,  5.66393572e+00,  3.75144548e+00,
-        2.83459537e-01, -1.61445732e+00, -2.61883874e+00, -2.76818696e+00,
-       -2.47679180e+00, -2.29014373e+00, -2.50794047e+00, -3.08754938e+00,
-       -3.77455164e+00, -4.26560645e+00, -4.33448168e+00, -3.93383052e+00,
-       -3.19736872e+00, -2.30738448e+00, -1.36778475e+00, -4.13019001e-01,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-       -2.30731635e-01, -2.71055349e+00, -3.77777456e+00, -3.19820518e+00,
-       -1.62052825e+00, -1.57786580e-01,  0.00000000e+00, -3.39194075e-01,
-       -1.19486567e+00, -1.40982868e+00, -8.24731765e-01,  1.77245380e-02,
-        0.00000000e+00,  0.00000000e+00,  1.12393272e+00,  3.42883960e+00,
-        6.44392419e+00,  7.74026022e+00,  6.02300041e+00,  2.15021680e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        1.31206244e+00,  2.88224433e+00,  2.70442915e+00,  1.31540273e+00,
-        6.83579157e-02, -1.44657955e-01, -9.74069771e-04,  2.01609711e-01,
-        7.58795510e-02, -1.14712672e-01,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00, -1.16918168e-01, -9.48748823e-02,
-       -3.74701248e-01, -6.68192840e-01, -7.77720853e-01, -8.15718667e-01,
-       -1.19404470e+00, -2.17926682e+00, -3.40215785e+00, -3.99418641e+00,
-       -3.25449483e+00, -1.22032570e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00, -1.72304245e+00,
-       -2.63776333e+00, -2.11110151e+00, -7.22241992e-01,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  1.37988828e-01,
-       -1.36497601e-01,  1.32278945e-01,  5.15571508e-01,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  8.42548433e-01,  2.72091640e+00,  3.24542027e+00,
-        2.03399132e+00,  3.77699186e-02,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00, -1.99799699e-01,
-        5.92260640e-02, -1.02272668e-02, -3.19939917e-01, -7.49813663e-01,
-       -1.06945353e+00, -1.04483902e+00, -5.36963585e-01,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-       -3.23192576e-01,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  1.42596962e-01, -1.71901114e-01,
-       -7.45893810e-02,  3.71940362e-01,  6.24286373e-01,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00, -9.16154067e-02, -8.72890060e-02,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  8.36067934e-01,  1.13854180e+00,
-        4.79721529e-01,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-        0.00000000e+00]),
-        abs=1e-6
-        # fmt: on
-    )
+    assert_allclose(rub.forces_rub[rub.posRUB * 6 + 0, :], Fx_rub, rtol=3e-2)
+    assert_allclose(rub.forces_rub[rub.posRUB * 6 + 1, :], Fy_rub, rtol=3e-2)
diff --git a/ross/tests/test_shaft_element.py b/ross/tests/test_shaft_element.py
index 3bed7ea9a..f88b720c6 100644
--- a/ross/tests/test_shaft_element.py
+++ b/ross/tests/test_shaft_element.py
@@ -405,118 +405,124 @@ def test_save_load(tim, tap_tim):
 
 
 @pytest.fixture
-def s6_eb():
-    #  Euler-Bernoulli element
-    """
-    L = 0.1
-    odr = 25.4 * 3 / 4
-    odl = 25.4 * 3 / 4
-    idr = 25.4 * 1 / 2
-    idl = 25.4 * 1 / 2
-    axial_force = 1
-    torque = 1
-    alpha = 2.7
-    beta = 4.8 * 10 ** (-6)
-    """
-
-    return ShaftElement6DoF(
-        L=0.1,
-        idl=25.4 * 1 / 2,
-        odl=25.4 * 3 / 4,
-        idr=25.4 * 1 / 2,
-        odr=25.4 * 3 / 4,
-        material=steel,
-        axial_force=1,
-        torque=1,
-        alpha=2.7,
-        beta=4.8 * 10 ** (-6),
-    )
-
-
-def test_s6_parameters_eb(s6_eb):
-    assert s6_eb.L == 0.1
-    assert s6_eb.idl == 25.4 * 1 / 2
-    assert s6_eb.odl == 25.4 * 3 / 4
-    assert s6_eb.idr == 25.4 * 1 / 2
-    assert s6_eb.odr == 25.4 * 3 / 4
-    assert s6_eb.material.E == 211e9
-    assert s6_eb.material.G_s == 81.2e9
-    assert s6_eb.material.rho == 7810
+def tim_6dof():
+    #  Timoshenko element
+    L = 0.25
+    i_d = 0
+    o_d = 0.05
+    return ShaftElement6DoF(L, i_d, o_d, material=steel, alpha=1, beta=1e-5)
+
+
+def test_parameters_tim_6dof(tim_6dof):
+    assert tim_6dof.L == 0.25
+    assert tim_6dof.idl == 0
+    assert tim_6dof.odl == 0.05
+    assert tim_6dof.idr == 0
+    assert tim_6dof.odr == 0.05
+    assert tim_6dof.material.E == 211e9
+    assert tim_6dof.material.G_s == 81.2e9
+    assert tim_6dof.material.rho == 7810
+    assert_almost_equal(tim_6dof.phi, 0.08795566)
+    assert_almost_equal(tim_6dof.material.Poisson, 0.29926108)
+    assert_almost_equal(tim_6dof.A, 0.00196349)
+    assert_almost_equal(tim_6dof.Ie * 1e7, 3.06796157)
+
+
+def test_mass_matrix_tim_6dof(tim_6dof):
+    # fmt: off
+    M0e_tim = np.array([
+        [ 1.42051,  0.00000,  0.00000,  0.00000,  0.04932,  0.0000,  0.49635,  0.00000,  0.00000,  0.00000, -0.03055,  0.0000],
+        [ 0.00000,  1.42051,  0.00000, -0.04932,  0.00000,  0.0000,  0.00000,  0.49635,  0.00000,  0.03055,  0.00000,  0.0000],
+        [ 0.00000,  0.00000,  1.27791,  0.00000,  0.00000,  0.0000,  0.00000,  0.00000,  0.63895,  0.00000,  0.00000,  0.0000],
+        [ 0.00000, -0.04932,  0.00000,  0.00231,  0.00000,  0.0000,  0.00000, -0.03055,  0.00000, -0.00178,  0.00000,  0.0000],
+        [ 0.04932,  0.00000,  0.00000,  0.00000,  0.00231,  0.0000,  0.03055,  0.00000,  0.00000,  0.00000, -0.00178,  0.0000],
+        [ 0.00000,  0.00000,  0.00000,  0.00000,  0.00000,  0.0004,  0.00000,  0.00000,  0.00000,  0.00000,  0.00000,  0.0002],
+        [ 0.49635,  0.00000,  0.00000,  0.00000,  0.03055,  0.0000,  1.42051,  0.00000,  0.00000,  0.00000, -0.04932,  0.0000],
+        [ 0.00000,  0.49635,  0.00000, -0.03055,  0.00000,  0.0000,  0.00000,  1.42051,  0.00000,  0.04932,  0.00000,  0.0000],
+        [ 0.00000,  0.00000,  0.63895,  0.00000,  0.00000,  0.0000,  0.00000,  0.00000,  1.27791,  0.00000,  0.00000,  0.0000],
+        [ 0.00000,  0.03055,  0.00000, -0.00178,  0.00000,  0.0000,  0.00000,  0.04932,  0.00000,  0.00231,  0.00000,  0.0000],
+        [-0.03055,  0.00000,  0.00000,  0.00000, -0.00178,  0.0000, -0.04932,  0.00000,  0.00000,  0.00000,  0.00231,  0.0000],
+        [ 0.00000,  0.00000,  0.00000,  0.00000,  0.00000,  0.0002,  0.00000,  0.00000,  0.00000,  0.00000,  0.00000,  0.0004],
+])
+    # fmt: on
+    assert_almost_equal(tim_6dof.M(), M0e_tim, decimal=5)
 
 
-def test_s6_M_matrix(s6_eb):
+def test_stiffness_matrix_tim_6dof(tim_6dof):
     # fmt: off
-    M_mat_6DoF = np.array(
-        [[  4.8624e+08,            0,            0,            0,  -4.0523e+06,            0,  -4.8618e+08,            0,            0,            0,  -4.0512e+06,            0],
-         [           0,   4.8624e+08,            0,   4.0523e+06,            0,            0,            0,  -4.8618e+08,            0,   4.0512e+06,            0,            0],
-         [           0,            0,   4.1223e+04,            0,            0,            0,            0,            0,   2.0611e+04,            0,            0,            0],
-         [           0,   4.0523e+06,            0,   5.4023e+05,            0,            0,            0,  -4.0512e+06,            0,  -1.3506e+05,            0,            0],
-         [ -4.0523e+06,            0,            0,            0,   5.4023e+05,            0,   4.0512e+06,            0,            0,            0,  -1.3506e+05,            0],
-         [           0,            0,            0,            0,            0,   1.3505e+06,            0,            0,            0,            0,            0,   6.7527e+05],
-         [ -4.8618e+08,            0,            0,            0,   4.0512e+06,            0,   4.8624e+08,            0,            0,            0,   4.0523e+06,            0],
-         [           0,  -4.8618e+08,            0,  -4.0512e+06,            0,            0,            0,   4.8624e+08,            0,  -4.0523e+06,            0,            0],
-         [           0,            0,   2.0611e+04,            0,            0,            0,            0,            0,   4.1223e+04,            0,            0,            0],
-         [           0,   4.0512e+06,            0,  -1.3506e+05,            0,            0,            0,  -4.0523e+06,            0,   5.4023e+05,            0,            0],
-         [ -4.0512e+06,            0,            0,            0,  -1.3506e+05,            0,   4.0523e+06,            0,            0,            0,   5.4023e+05,            0],
-         [           0,            0,            0,            0,            0,   6.7527e+05,            0,            0,            0,            0,            0,   1.3505e+06]])
+    K0e_tim = np.array([
+        [ 0.45696,  0.00000,   0.0000,  0.00000,  0.05712,  0.00000, -0.45696,  0.00000,   0.0000,  0.00000,  0.05712,  0.00000],
+        [ 0.00000,  0.45696,   0.0000, -0.05712,  0.00000,  0.00000,  0.00000, -0.45696,   0.0000, -0.05712,  0.00000,  0.00000],
+        [ 0.00000,  0.00000,  16.5719,  0.00000,  0.00000,  0.00000,  0.00000,  0.00000, -16.5719,  0.00000,  0.00000,  0.00000],
+        [ 0.00000, -0.05712,   0.0000,  0.00973,  0.00000,  0.00000,  0.00000,  0.05712,   0.0000,  0.00455,  0.00000,  0.00000],
+        [ 0.05712,  0.00000,   0.0000,  0.00000,  0.00973,  0.00000, -0.05712,  0.00000,   0.0000,  0.00000,  0.00455,  0.00000],
+        [ 0.00000,  0.00000,   0.0000,  0.00000,  0.00000,  0.00199,  0.00000,  0.00000,   0.0000,  0.00000,  0.00000, -0.00199],
+        [-0.45696,  0.00000,   0.0000,  0.00000, -0.05712,  0.00000,  0.45696,  0.00000,   0.0000,  0.00000, -0.05712,  0.00000],
+        [ 0.00000, -0.45696,   0.0000,  0.05712,  0.00000,  0.00000,  0.00000,  0.45696,   0.0000,  0.05712,  0.00000,  0.00000],
+        [ 0.00000,  0.00000, -16.5719,  0.00000,  0.00000,  0.00000,  0.00000,  0.00000,  16.5719,  0.00000,  0.00000,  0.00000],
+        [ 0.00000, -0.05712,   0.0000,  0.00455,  0.00000,  0.00000,  0.00000,  0.05712,   0.0000,  0.00973,  0.00000,  0.00000],
+        [ 0.05712,  0.00000,   0.0000,  0.00000,  0.00455,  0.00000, -0.05712,  0.00000,   0.0000,  0.00000,  0.00973,  0.00000],
+        [ 0.00000,  0.00000,   0.0000,  0.00000,  0.00000, -0.00199,  0.00000,  0.00000,   0.0000,  0.00000,  0.00000,  0.00199],
+])
     # fmt: on
-    assert_allclose(s6_eb.M(), M_mat_6DoF, rtol=1e-3)
+    assert_almost_equal(tim_6dof.K() / 1e8, K0e_tim, decimal=5)
 
 
-def test_s6_G_matrix(s6_eb):
+def test_dynamic_stiffness_matrix_tim_6dof(tim_6dof):
     # fmt: off
-    G_mat_6DoF = np.array(
-        [[           0,   3.3495e-02,            0,   2.3779e+02,            0,            0,            0,   -3.3495e-02,           0,   2.3779e+02,            0,            0],
-         [ -3.3495e-02,            0,            0,            0,   2.3779e+02,            0,   3.3495e-02,            0,            0,            0,   2.3779e+02,            0],
-         [           0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
-         [-2.37798e+02,            0,            0,            0,   2.7010e+06,            0,   2.3779e+02,            0,            0,            0,   1.3505e+06,            0],
-         [           0,  -2.3779e+02,            0,  -2.7010e+06,            0,            0,            0,   2.3779e+02,            0,  -1.3505e+06,            0,            0],
-         [           0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
-         [           0,  -3.3495e-02,            0,  -2.3779e+02,            0,            0,            0,   3.3495e-02,            0,  -2.3779e+02,            0,            0],
-         [  3.3495e-02,            0,            0,            0,  -2.3779e+02,            0,  -3.3495e-02,            0,            0,            0,  -2.3779e+02,            0],
-         [           0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
-         [ -2.3779e+02,            0,            0,            0,   1.3505e+06,            0,   2.3779e+02,            0,            0,            0,   2.7010e+06,            0],
-         [           0,  -2.3779e+02,            0,  -1.3505e+06,            0,            0,            0,   2.3779e+02,            0,  -2.7010e+06,            0,            0],
-         [           0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0]])
-
+    Kst0e_tim = np.array([
+        [0.0, -23.00235,  0.0,  0.47922,  0.0,  0.0,  0.0,  23.00235,  0.0,  0.47922,  0.0,  0.0],
+        [0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0,  0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0],
+        [0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0,  0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0],
+        [0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0,  0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0],
+        [0.0 , -0.47922,  0.0,  0.15974,  0.0,  0.0,  0.0,   0.47922,  0.0, -0.03993,  0.0,  0.0],
+        [0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0,  0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0],
+        [0.0,  23.00235,  0.0, -0.47922,  0.0,  0.0,  0.0, -23.00235,  0.0, -0.47922,  0.0,  0.0],
+        [0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0,  0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0],
+        [0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0,  0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0],
+        [0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0,  0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0],
+        [0.0 , -0.47922,  0.0, -0.03993,  0.0,  0.0,  0.0,   0.47922,  0.0,  0.15974,  0.0,  0.0],
+        [0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0,  0.0,   0.00000,  0.0,  0.00000,  0.0,  0.0],
+])
     # fmt: on
-    assert_allclose(s6_eb.G(), G_mat_6DoF, rtol=1e-3)
+    assert_almost_equal(tim_6dof.Kst() * 1e3, Kst0e_tim, decimal=5)
 
 
-def test_s6_K_matrix(s6_eb):
+def test_damping_matrix_tim_6dof(tim_6dof):
     # fmt: off
-    K_mat_6DoF = np.array(
-        [[ 7.7093e+13,            0,            0,  -1.0000e+01,  -3.8547e+12,            0,  -7.7093e+13,            0,            0,   1.0000e+01,  -3.8547e+12,            0],
-         [          0,   7.7093e+13,            0,   3.8547e+12,  -1.0000e+01,            0,            0,  -7.7093e+13,            0,   3.8547e+12,   1.0000e+01,            0],
-         [          0,            0,   3.3411e+14,            0,            0,            0,            0,            0,  -3.3411e+14,            0,            0,            0],
-         [-1.0000e+01,   3.8547e+12,            0,   1.0946e+16,   5.0000e-01,            0,   1.0000e+01,  -3.8547e+12,            0,  -1.0946e+16,   5.0000e-01,            0],
-         [-3.8547e+12,  -1.0000e+01,            0,  -5.0000e-01,   1.0946e+16,            0,   3.8547e+12,   1.0000e+01,            0,  -5.0000e-01,  -1.0946e+16,            0],
-         [          0,            0,            0,            0,            0,   8.4249e+15,            0,            0,            0,            0,            0,  -8.4249e+15],
-         [-7.7093e+13,            0,            0,   1.0000e+01,   3.8547e+12,            0,   7.7093e+13,            0,            0,  -1.0000e+01,   3.8547e+12,            0],
-         [          0,  -7.7093e+13,            0,  -3.8547e+12,   1.0000e+01,            0,            0,   7.7093e+13,            0,  -3.8547e+12,  -1.0000e+01,            0],
-         [          0,            0,  -3.3411e+14,            0,            0,            0,            0,            0,   3.3411e+14,            0,            0,            0],
-         [ 1.0000e+01,   3.8547e+12,            0,  -1.0946e+16,  -5.0000e-01,            0,  -1.0000e+01,  -3.8547e+12,            0,   1.0946e+16,  -5.0000e-01,            0],
-         [-3.8547e+12,   1.0000e+01,            0,   5.0000e-01,  -1.0946e+16,            0,   3.8547e+12,  -1.0000e+01,            0,   5.0000e-01,   1.0946e+16,            0],
-         [          0,            0,            0,            0,            0,  -8.4249e+15,            0,            0,            0,            0,            0,   8.4249e+15]])
-
+    C0e_tim = np.array([
+        [ 0.45838,  0.00000,   0.00000,  0.00000,  0.05717,  0.00000, -0.45647,  0.00000,   0.00000,  0.00000,  0.05709,  0.00000],
+        [ 0.00000,  0.45838,   0.00000, -0.05717,  0.00000,  0.00000,  0.00000, -0.45647,   0.00000, -0.05709,  0.00000,  0.00000],
+        [ 0.00000,  0.00000,  16.57318,  0.00000,  0.00000,  0.00000,  0.00000,  0.00000, -16.57126,  0.00000,  0.00000,  0.00000],
+        [ 0.00000, -0.05717,   0.00000,  0.00973,  0.00000,  0.00000,  0.00000,  0.05709,   0.00000,  0.00455,  0.00000,  0.00000],
+        [ 0.05717,  0.00000,   0.00000,  0.00000,  0.00973,  0.00000, -0.05709,  0.00000,   0.00000,  0.00000,  0.00455,  0.00000],
+        [ 0.00000,  0.00000,   0.00000,  0.00000,  0.00000,  0.00199,  0.00000,  0.00000,   0.00000,  0.00000,  0.00000, -0.00199],
+        [-0.45647,  0.00000,   0.00000,  0.00000, -0.05709,  0.00000,  0.45838,  0.00000,   0.00000,  0.00000, -0.05717,  0.00000],
+        [ 0.00000, -0.45647,   0.00000,  0.05709,  0.00000,  0.00000,  0.00000,  0.45838,   0.00000,  0.05717,  0.00000,  0.00000],
+        [ 0.00000,  0.00000, -16.57126,  0.00000,  0.00000,  0.00000,  0.00000,  0.00000,  16.57318,  0.00000,  0.00000,  0.00000],
+        [ 0.00000, -0.05709,   0.00000,  0.00455,  0.00000,  0.00000,  0.00000,  0.05717,   0.00000,  0.00973,  0.00000,  0.00000],
+        [ 0.05709,  0.00000,   0.00000,  0.00000,  0.00455,  0.00000, -0.05717,  0.00000,   0.00000,  0.00000,  0.00973,  0.00000],
+        [ 0.00000,  0.00000,   0.00000,  0.00000,  0.00000, -0.00199,  0.00000,  0.00000,   0.00000,  0.00000,  0.00000,  0.00199],
+])
     # fmt: on
-    assert_allclose(s6_eb.K(), K_mat_6DoF, rtol=1e-3)
+    assert_almost_equal(tim_6dof.C() / 1e3, C0e_tim, decimal=5)
 
 
-def test_s6_Kst_matrix(s6_eb):
+def test_gyroscopic_matrix_tim_6dof(tim_6dof):
     # fmt: off
-    Kst_mat_6DoF = 0.01 * np.array([
-        [          0,  -9.7239e+10,            0,  -8.1032e+08,            0,            0,            0,   9.7239e+10,            0,  -8.1032e+08,            0,            0],
-        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
-        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
-        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
-        [          0,   8.1032e+08,            0,   1.0804e+08,            0,            0,            0,  -8.1032e+08,            0,  -2.7011e+07,            0,            0],
-        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
-        [          0,   9.7239e+10,            0,   8.1032e+08,            0,            0,            0,  -9.7239e+10,            0,   8.1032e+08,            0,            0],
-        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
-        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
-        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
-        [          0,   8.1032e+08,            0,  -2.7011e+07,            0,            0,            0,  -8.1032e+08,            0,   1.0804e+08,            0,            0],
-        [          0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0]])
+    G0e_tim = np.array([
+        [  0.00000,  19.43344,  0.0, -0.22681,  0.00000,  0.0,   0.00000, -19.43344,  0.0, -0.22681,  0.00000,  0.0],
+        [-19.43344,   0.00000,  0.0,  0.00000, -0.22681,  0.0,  19.43344,   0.00000,  0.0,  0.00000, -0.22681,  0.0],
+        [  0.00000,   0.00000,  0.0,  0.00000,  0.00000,  0.0,   0.00000,   0.00000,  0.0,  0.00000,  0.00000,  0.0],
+        [  0.22681,   0.00000,  0.0,  0.00000,  0.15240,  0.0 , -0.22681,   0.00000,  0.0,  0.00000, -0.04727,  0.0],
+        [  0.00000,   0.22681,  0.0, -0.15240,  0.00000,  0.0,   0.00000 , -0.22681,  0.0,  0.04727,  0.00000,  0.0],
+        [  0.00000,   0.00000,  0.0,  0.00000,  0.00000,  0.0,   0.00000,   0.00000,  0.0,  0.00000,  0.00000,  0.0],
+        [  0.00000, -19.43344,  0.0,  0.22681,  0.00000,  0.0,   0.00000,  19.43344,  0.0,  0.22681,  0.00000,  0.0],
+        [ 19.43344,   0.00000,  0.0,  0.00000,  0.22681,  0.0, -19.43344,   0.00000,  0.0,  0.00000,  0.22681,  0.0],
+        [  0.00000,   0.00000,  0.0,  0.00000,  0.00000,  0.0,   0.00000,   0.00000,  0.0,  0.00000,  0.00000,  0.0],
+        [  0.22681,   0.00000,  0.0,  0.00000, -0.04727,  0.0 , -0.22681,   0.00000,  0.0,  0.00000,  0.15240,  0.0],
+        [  0.00000,   0.22681,  0.0,  0.04727,  0.00000,  0.0,   0.00000 , -0.22681,  0.0, -0.15240,  0.00000,  0.0],
+        [  0.00000,   0.00000,  0.0,  0.00000,  0.00000,  0.0,   0.00000,   0.00000,  0.0,  0.00000,  0.00000,  0.0],
+])
     # fmt: on
-    assert_allclose(s6_eb.Kst(), Kst_mat_6DoF, rtol=1e-3)
+    assert_almost_equal(tim_6dof.G() * 1e3, G0e_tim, decimal=5)