Repository: ross-rotordynamics/ross
Pull request #535 (2020-04-24 20:30:11)

Plotting stochastic rotor orbit
This PR adds methods to plot stochastic orbits using `.run_time_response()` method from `ST_Rotor` class.

EDIT:
Just to clarify how it works, The user runs a stochastic time response analysis and it's given 3 options for plotting the results:

1. plotting 1D: amplitude vs time for a given degree of freedom.
2. plotting 2D: orbit of a node from the rotor.
3. plotting 3D: orbits of all nodes from the rotor. 

Example:
```
        >>> import ross.stochastic as srs
        >>> rotors = srs.st_rotor_example()

        # Running Time Response and saving the results

        >>> size = 1000
        >>> ndof = rotors.ndof
        >>> node = 3 # node where the force is applied

        >>> dof = 9
        >>> speed = 250.0
        >>> t = np.linspace(0, 10, size)
        >>> F = np.zeros((size, ndof))
        >>> F[:, 4 * node] = 10 * np.cos(2 * t)
        >>> F[:, 4 * node + 1] = 10 * np.sin(2 * t)
        >>> results = rotors.run_time_response(speed, F, t)

        # Plotting Time vs Amplitude (the user adds a dof)
        >>> results.plot(plot_type="1d", dof=dof, conf_interval=[90])

        # Plotting orbit of a rotor node (the user adds a node)
        >>> results.plot(plot_type="2d", node=node, conf_interval=[90])

        # Plotting orbit of all nodes
        >>> results.plot(plot_type="3d", conf_interval=[90])
        <matplotlib...
```


Patch:
diff --git a/.gitignore b/.gitignore
index 4e8193c7..81034dd3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -111,23 +111,8 @@ Thumbs.db
 notebooks/
 
 # docs #
-########
 /docs/examples/*_files
 /docs/generated/
 /paper/whedon/
 /paper/paper.pdf
 /paper/Makefile
-
-### VirtualEnv ###
-##################
-.Python
-pyvenv.cfg
-.env
-.venv
-env/
-venv/
-ENV/
-env.bak/
-venv.bak/
-pip-selfcheck.json
-
diff --git a/docs/download.rst b/docs/download.rst
index 0e46bf69..a31b1dd5 100644
--- a/docs/download.rst
+++ b/docs/download.rst
@@ -12,4 +12,3 @@ Here you can download the tutorial and examples notebooks.
 | :download:`Example 6 </examples/example_6.ipynb>`
 | :download:`Example 7 </examples/example_7.ipynb>`
 | :download:`Example 8 </examples/example_8.ipynb>`
-| :download:`Stochastic ROSS - Tutorial </examples/example_stochastic_ross.ipynb>`
diff --git a/docs/examples.rst b/docs/examples.rst
index 642ab9c2..717b8e93 100644
--- a/docs/examples.rst
+++ b/docs/examples.rst
@@ -11,5 +11,4 @@ Examples
    examples/example_5
    examples/example_6
    examples/example_7
-   examples/example_8
-   examples/example_stochastic_ross
\ No newline at end of file
+   examples/example_8
\ No newline at end of file
diff --git a/docs/examples/example_stochastic_ross.ipynb b/docs/examples/tutorial_stochastic_ross.ipynb
similarity index 100%
rename from docs/examples/example_stochastic_ross.ipynb
rename to docs/examples/tutorial_stochastic_ross.ipynb
diff --git a/docs/refs.bib b/docs/refs.bib
index 1290c2b4..6d2cd068 100644
--- a/docs/refs.bib
+++ b/docs/refs.bib
@@ -27,13 +27,6 @@
   publisher={Cambridge university press}
 }
 
-@book{lalanne2001,
-  title={Rotordynamics Prediction in Engineering},
-  author={M Lalanne and G Ferraris},
-  year={1998},
-  publisher={John Wiley and Sons
-}
-
 @article{zang2005review,
   title={A review of robust optimal design and its application in dynamics},
   author={Zang, C and Friswell, MI and Mottershead, JE},
@@ -634,10 +627,3 @@ year={2011}
 	year={2012}
 }
 
-@article{Hutchingson2001,
-	title={Shear Coefficients for Timoshenko Beam Theory},
-	author={Hutchingson, J. R.},
-	journal={J. Appl. Mech.},
-	volume={68},
-	year={2001}
-}
\ No newline at end of file
diff --git a/requirements.txt b/requirements.txt
index 95fdafc1..f1c7b782 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,3 @@
-
 numpy
 scipy
 matplotlib
@@ -6,4 +5,4 @@ toml
 pandas>=0.23
 bokeh
 xlrd
-pint
+pint
\ No newline at end of file
diff --git a/ross/__init__.py b/ross/__init__.py
index 3093ba63..97d27422 100644
--- a/ross/__init__.py
+++ b/ross/__init__.py
@@ -1,4 +1,4 @@
-__version__ = "0.3.3"
+__version__ = "0.3.2"
 
 from .api_report import *
 from .bearing_seal_element import *
diff --git a/ross/bearing_seal_element.py b/ross/bearing_seal_element.py
index 33524cac..25981fb4 100644
--- a/ross/bearing_seal_element.py
+++ b/ross/bearing_seal_element.py
@@ -21,7 +21,6 @@ __all__ = [
     "SealElement",
     "BallBearingElement",
     "RollerBearingElement",
-    "BearingElement6DoF",
     "MagneticBearingElement",
 ]
 # fmt: on
@@ -419,7 +418,7 @@ class BearingElement(Element):
         return dict(x_0=0, y_0=1)
 
     def dof_global_index(self):
-        """Get the global index for an element specific degree of freedom."""
+        """Get the global index for a element specific degree of freedom."""
         global_index = super().dof_global_index()
 
         if self.n_link is not None:
@@ -1454,357 +1453,6 @@ class MagneticBearingElement(BearingElement):
         )
 
 
-class BearingElement6DoF(BearingElement):
-    """A generalistic 6 DoF bearing element. This class will create a bearing
-    element based on the user supplied stiffness and damping coefficients. These
-    are determined alternatively, via purposefully built codes.
-
-    Parameters
-    ----------
-    kxx: float, optional
-        Direct stiffness in the x direction.
-        Default is None
-    kxy: float, optional
-        Cross stiffness between xy directions.
-        Default is None
-    kyx: float, optional
-        Cross stiffness between yx directions.
-        Default is None
-    kyy: float, optional
-        Direct stiffness in the y direction.
-        Default is None
-    kzz: float, optional
-        Direct stiffness in the z direction.
-        Default is None
-    cxx: float, optional
-        Direct damping in the x direction.
-        Default is None
-    cxy: float, optional
-        Cross damping between xy directions.
-        Default is None
-    cyx: float, optional
-        Cross damping between yx directions.
-        Default is None
-    cyy: float, optional
-        Direct damping in the y direction.
-        Default is None
-    czz: float, optional
-        Direct damping in the z direction.
-        Default is None
-    tag : str, optional
-        A tag to name the element
-        Default is None
-    Examples
-    --------
-    >>> n = 0
-    >>> kxx = 1.0e7
-    >>> kyy = 1.5e7
-    >>> kzz = 5.0e5
-    >>> bearing = BearingElement6DoF(n=n, kxx=kxx, kyy=kyy, kzz=kzz,
-    ...                              cxx=0, cyy=0)
-    >>> bearing.K(0)
-    array([[10000000.,        0.,        0.],
-           [       0., 15000000.,        0.],
-           [       0.,        0.,   500000.]])
-    """
-
-    def __init__(
-        self,
-        n,
-        kxx,
-        kyy,
-        cxx,
-        cyy,
-        kxy=0.0,
-        kyx=0.0,
-        kzz=0.0,
-        cxy=0.0,
-        cyx=0.0,
-        czz=0.0,
-        tag=None,
-        scale_factor=1,
-        color="#355d7a",
-    ):
-        super().__init__(
-            n=n,
-            kxx=kxx,
-            cxx=cxx,
-            kyy=kyy,
-            kxy=kxy,
-            kyx=kyx,
-            cyy=cyy,
-            cxy=cxy,
-            cyx=cyx,
-            frequency=None,
-            tag=tag,
-            n_link=None,
-            scale_factor=scale_factor,
-            color=color,
-        )
-
-        new_args = ["kzz", "czz"]
-
-        args_dict = locals()
-        coefficients = {}
-
-        if kzz is None:
-            args_dict["kzz"] = (
-                kxx * 0.6
-            )  # NSK manufacturer sugestion for deep groove ball bearings
-        if czz is None:
-            args_dict["czz"] = cxx
-
-        for arg in new_args:
-            if arg[0] == "k":
-                coefficients[arg] = _Stiffness_Coefficient(
-                    coefficient=args_dict[arg], w=None
-                )
-            else:
-                coefficients[arg] = _Damping_Coefficient(args_dict[arg], None)
-
-        coefficients_len = [len(v.coefficient) for v in coefficients.values()]
-
-        for c in coefficients_len:
-            if c != 1:
-                raise ValueError(
-                    "Arguments (coefficients and frequency)"
-                    " must have the same dimension"
-                )
-
-        for k, v in coefficients.items():
-            setattr(self, k, v)
-
-    def __repr__(self):
-        """This function returns a string representation of a bearing element.
-        Parameters
-        ----------
-
-        Returns
-        -------
-        A string representation of a bearing object.
-        Examples
-        --------
-        >>> bearing = bearing_example()
-        >>> bearing # doctest: +ELLIPSIS
-        BearingElement(n=0, n_link=None,
-         kxx=[...
-        """
-        return (
-            f"{self.__class__.__name__}"
-            f"(n={self.n}, n_link={self.n_link},\n"
-            f" kxx={self.kxx}, kxy={self.kxy},\n"
-            f" kyx={self.kyx}, kyy={self.kyy},\n"
-            f" kzz={self.kzz}, cxx={self.cxx},\n"
-            f" cxy={self.cxy}, cyx={self.cyx},\n"
-            f" cyy={self.cyy}, czz={self.czz},\n"
-            f" frequency={self.frequency}, tag={self.tag!r})"
-        )
-
-    def __eq__(self, other):
-        """This function allows bearing elements to be compared.
-        Parameters
-        ----------
-        other: object
-            The second object to be compared with.
-
-        Returns
-        -------
-        bool
-            True if the comparison is true; False otherwise.
-        Examples
-        --------
-        >>> bearing1 = bearing_example()
-        >>> bearing2 = bearing_example()
-        >>> bearing1 == bearing2
-        True
-        """
-        compared_attributes = [
-            "kxx",
-            "kyy",
-            "kxy",
-            "kyx",
-            "cxx",
-            "cyy",
-            "cxy",
-            "cyx",
-            "kzz",
-            "czz",
-            "frequency",
-            "n",
-        ]
-        if isinstance(other, self.__class__):
-            return all(
-                (
-                    np.array(getattr(self, attr)).all()
-                    == np.array(getattr(other, attr)).all()
-                    for attr in compared_attributes
-                )
-            )
-        return False
-
-    def save(self, file_name=Path(os.getcwd())):
-        """Saves a bearing element in a toml format. It works as an auxiliary
-        function of the save function in the Rotor class.
-
-        Parameters
-        ----------
-        file_name: string
-            The name of the file the bearing element will be saved in.
-
-        Returns
-        -------
-        None
-
-        Examples
-        --------
-        >>> bearing = bearing_example()
-        >>> bearing.save(Path(os.getcwd()))
-        """
-        data = self.get_data(Path(file_name) / "BearingElement6DoF.toml")
-
-        if type(self.frequency) == np.ndarray:
-            try:
-                self.frequency[0]
-                frequency = list(self.frequency)
-            except IndexError:
-                frequency = None
-
-        data["BearingElement6DoF"][str(self.n)] = {
-            "n": self.n,
-            "kxx": self.kxx.coefficient,
-            "cxx": self.cxx.coefficient,
-            "kyy": self.kyy.coefficient,
-            "kxy": self.kxy.coefficient,
-            "kyx": self.kyx.coefficient,
-            "kzz": self.kzz.coefficient,
-            "cyy": self.cyy.coefficient,
-            "cxy": self.cxy.coefficient,
-            "cyx": self.cyx.coefficient,
-            "czz": self.czz.coefficient,
-            "frequency": frequency,
-            "tag": self.tag,
-        }
-        self.dump_data(data, Path(file_name) / "BearingElement6DoF.toml")
-
-    @staticmethod
-    def load(file_name=""):
-        """Loads a list of bearing elements saved in a toml format.
-
-        Parameters
-        ----------
-        file_name: string
-            The name of the file of the bearing element to be loaded.
-
-        Returns
-        -------
-        A list of bearing elements.
-
-        Examples
-        --------
-
-        """
-        bearing_elements = []
-        bearing_elements_dict = BearingElement.get_data(
-            file_name=Path(file_name) / "BearingElement6DoF.toml"
-        )
-        for element in bearing_elements_dict["BearingElement6DoF"]:
-            bearing = BearingElement6DoF(
-                **bearing_elements_dict["BearingElement6DoF"][element]
-            )
-            data = bearing_elements_dict["BearingElement6DoF"]
-            bearing.kxx.coefficient = data[element]["kxx"]
-
-            bearing.kxy.coefficient = data[element]["kxy"]
-
-            bearing.kyx.coefficient = data[element]["kyx"]
-
-            bearing.kyy.coefficient = data[element]["kyy"]
-
-            bearing.kzz.coefficient = data[element]["kzz"]
-
-            bearing.cxx.coefficient = data[element]["cxx"]
-
-            bearing.cxy.coefficient = data[element]["cxy"]
-
-            bearing.cyx.coefficient = data[element]["cyx"]
-
-            bearing.cyy.coefficient = data[element]["cyy"]
-
-            bearing.czz.coefficient = data[element]["czz"]
-
-            bearing_elements.append(bearing)
-        return bearing_elements
-
-    def dof_mapping(self):
-        return dict(x_0=0, y_0=1, z_0=2)
-
-    def dof_global_index(self):
-        """Get the global index for an element specific degree of freedom for the 6DoF element."""
-        global_index = super().dof_global_index()
-
-        if self.n_link is not None:
-            global_index = global_index._asdict()
-            global_index[f"x_{self.n_link}"] = 6 * self.n_link
-            global_index[f"y_{self.n_link}"] = 6 * self.n_link + 1
-            global_index[f"z_{self.n_link}"] = 6 * self.n_link + 2
-            dof_tuple = namedtuple("GlobalIndex", global_index)
-            global_index = dof_tuple(**global_index)
-
-        return global_index
-
-    def K(self, frequency):
-        """Returns the stiffness matrix for a given excitation frequency.
-
-        Returns
-        -------
-        K: np.ndarray
-            A 3x3 matrix of floats containing the kxx, kxy, kyx, kyy and kzz values.
-
-        Examples
-        --------
-        >>> bearing = bearing_6dof_example()
-        >>> bearing.K(0)
-        array([[1000000.,       0.,       0.],
-               [      0.,  800000.,       0.],
-               [      0.,       0.,  100000.]])
-        """
-        kxx = self.kxx.interpolated(frequency)
-        kyy = self.kyy.interpolated(frequency)
-        kxy = self.kxy.interpolated(frequency)
-        kyx = self.kyx.interpolated(frequency)
-        kzz = self.kzz.interpolated(frequency)
-
-        K = np.array([[kxx, kxy, 0], [kyx, kyy, 0], [0, 0, kzz]])
-
-        return K
-
-    def C(self, frequency):
-        """Returns the damping matrix for a given excitation frequency.
-
-        Returns
-        -------
-        C: np.ndarray
-            A 3x3 matrix of floats containing the cxx, cxy, cyx, cyy, and czz values.
-
-        Examples
-        --------
-        >>> bearing = bearing_6dof_example()
-        >>> bearing.C(0)
-        array([[200.,   0.,   0.],
-               [  0., 150.,   0.],
-               [  0.,   0.,  50.]])
-        """
-        cxx = self.cxx.interpolated(frequency)
-        cyy = self.cyy.interpolated(frequency)
-        cxy = self.cxy.interpolated(frequency)
-        cyx = self.cyx.interpolated(frequency)
-        czz = self.czz.interpolated(frequency)
-
-        C = np.array([[cxx, cxy, 0], [cyx, cyy, 0], [0, 0, czz]])
-
-        return C
-
-
 def bearing_example():
     """This function returns an instance of a simple bearing.
     The purpose is to make available a simple model
@@ -1857,33 +1505,3 @@ def seal_example():
     w = np.linspace(0, 200, 11)
     seal = SealElement(n=0, kxx=kxx, kyy=kyy, cxx=cxx, cyy=cyy, frequency=w)
     return seal
-
-
-def bearing_6dof_example():
-    """This function returns an instance of a simple bearing.
-    The purpose is to make available a simple model
-    so that doctest can be written using it.
-
-    Parameters
-    ----------
-
-    Returns
-    -------
-    An instance of a bearing object.
-
-    Examples
-    --------
-    >>> bearing = bearing_example()
-    >>> bearing.frequency[0]
-    0.0
-    """
-    kxx = 1e6
-    kyy = 0.8e6
-    kzz = 1e5
-    cxx = 2e2
-    cyy = 1.5e2
-    czz = 0.5e2
-    bearing = BearingElement6DoF(
-        n=0, kxx=kxx, kyy=kyy, cxx=cxx, cyy=cyy, kzz=kzz, czz=czz
-    )
-    return bearing
diff --git a/ross/disk_element.py b/ross/disk_element.py
index 0363e8df..eb869324 100644
--- a/ross/disk_element.py
+++ b/ross/disk_element.py
@@ -10,7 +10,7 @@ from bokeh.models import ColumnDataSource, HoverTool
 from ross.element import Element
 from ross.utils import read_table_file
 
-__all__ = ["DiskElement", "DiskElement6DoF"]
+__all__ = ["DiskElement"]
 bokeh_colors = bp.RdGy[11]
 
 
@@ -557,227 +557,3 @@ def disk_example():
     """
     disk = DiskElement(0, 32.589_727_65, 0.178_089_28, 0.329_563_62)
     return disk
-
-
-class DiskElement6DoF(DiskElement):
-    """A disk element for 6 DoFs.
-
-     This class will create a disk element with 6 DoF from input data of inertia and mass.
-
-     Parameters
-     ----------
-     n: int
-         Node in which the disk will be inserted.
-     m : float
-         Mass of the disk element.
-     Id : float
-         Diametral moment of inertia.
-     Ip : float
-         Polar moment of inertia
-     tag : str, optional
-         A tag to name the element
-         Default is None
-
-     Examples
-     --------
-     >>> disk = DiskElement6DoF(n=0, m=32, Id=0.2, Ip=0.3)
-     >>> disk.Ip
-     0.3
-     """
-
-    def dof_mapping(self):
-        """6DoFs 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,u_0 - horizontal translation, 
-        y_0,v_0 - vertical translation, 
-        z_0,w_0 - axial translation, 
-        theta_0 - rotation around horizontal, 
-        psi_0   - rotation around vertical, 
-        phi_0   - torsion around axial,
-        """
-        return dict(u_0=0, v_0=1, w_0=2, theta_0=3, psi_0=4, phi_0=5,)
-
-    def M(self):
-        """6DoFs mass matrix.
-
-        This method will return the mass matrix for an instance of a disk
-        element with 6DoFs.
-
-        Returns
-        -------
-        Mass matrix for the 6DoFs disk element.
-
-        Examples
-        --------
-        >>> disk = DiskElement6DoF(0, 32.58972765, 0.17808928, 0.32956362)
-        >>> disk.M().round(2)
-        array([[32.59,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
-               [ 0.  , 32.59,  0.  ,  0.  ,  0.  ,  0.  ],
-               [ 0.  ,  0.  , 32.59,  0.  ,  0.  ,  0.  ],
-               [ 0.  ,  0.  ,  0.  ,  0.18,  0.  ,  0.  ],
-               [ 0.  ,  0.  ,  0.  ,  0.  ,  0.18,  0.  ],
-               [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.33]])
-        """
-        m = self.m
-        Id = self.Id
-        Ip = self.Ip
-        # fmt: off
-        M = np.array([
-            [  m,  0,  0,  0,  0,  0],
-            [  0,  m,  0,  0,  0,  0],
-            [  0,  0,  m,  0,  0,  0],
-            [  0,  0,  0, Id,  0,  0],
-            [  0,  0,  0,  0, Id,  0],
-            [  0,  0,  0,  0,  0, Ip]])
-        # fmt: on
-        return M
-
-    def K(self):
-        """6DoFs stiffness matrix.
-
-        This method will return the stiffness matrix for an instance of a disk
-        element with 6DoFs.
-
-        Returns
-        -------
-        K: np.ndarray
-            A matrix of floats containing the values of the stiffness matrix.
-
-        Examples
-        --------
-        >>> disk = disk_example_6dof()
-        >>> disk.K().round(2)
-        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.33, 0.  , 0.  ]])
-        """
-        Ip = self.Ip
-        # fmt: off
-        K = 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, Ip, 0, 0]])
-        # fmt: on
-        return K
-
-    def C(self):
-        """6DoFs damping matrix.
-
-        Returns
-        -------
-        C: np.ndarray
-            A matrix of floats containing the values of the damping matrix.
-
-        Examples
-        --------
-        >>> disk = disk_example_6dof()
-        >>> disk.C()
-        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.]])
-        """
-        # fmt: off
-        C = np.zeros((6, 6))
-        # fmt: on
-        return C
-
-    def G(self):
-        """6DoFs gyroscopic matrix.
-
-        This method will return the gyroscopic matrix for an instance of a disk
-        element.
-
-        Returns
-        -------
-        G: np.ndarray
-            Gyroscopic matrix for the disk element.
-
-        Examples
-        --------
-        >>> disk = DiskElement6DoF(0, 32.58972765, 0.17808928, 0.32956362)
-        >>> disk.G().round(2)
-        array([[ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
-               [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
-               [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
-               [ 0.  ,  0.  ,  0.  ,  0.  ,  0.33,  0.  ],
-               [ 0.  ,  0.  ,  0.  , -0.33,  0.  ,  0.  ],
-               [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ]])
-        """
-        Ip = self.Ip
-        # fmt: off
-        G = np.array([
-            [ 0, 0, 0,   0,  0, 0],
-            [ 0, 0, 0,   0,  0, 0],
-            [ 0, 0, 0,   0,  0, 0],
-            [ 0, 0, 0,   0, Ip, 0],
-            [ 0, 0, 0, -Ip,  0, 0],
-            [ 0, 0, 0,   0,  0, 0]])
-        # fmt: on
-        return G
-
-
-def disk_example():
-    """This function returns an instance of a simple disk.
-    The purpose is to make available a simple model
-    so that doctest can be written using it.
-
-    Parameters
-    ----------
-
-    Returns
-    -------
-    An instance of a disk object.
-
-    Examples
-    --------
-    >>> disk = disk_example()
-    >>> disk.Ip
-    0.32956362
-    """
-    disk = DiskElement(0, 32.589_727_65, 0.178_089_28, 0.329_563_62)
-    return disk
-
-
-def disk_example_6dof():
-    """This function returns an instance of a simple disk.
-    The purpose is to make available a simple 6DoFs model
-    so that doctest can be written using it.
-
-    Parameters
-    ----------
-
-    Returns
-    -------
-    An instance of a disk object.
-
-    Examples
-    --------
-    >>> disk = disk_example()
-    >>> disk.Ip
-    0.32956362
-    """
-    disk = DiskElement6DoF(0, 32.589_727_65, 0.178_089_28, 0.329_563_62)
-    return disk
diff --git a/ross/shaft_element.py b/ross/shaft_element.py
index 5eec886b..d305f3ad 100644
--- a/ross/shaft_element.py
+++ b/ross/shaft_element.py
@@ -12,7 +12,7 @@ from ross.materials import Material, steel
 from ross.units import check_units
 from ross.utils import read_table_file
 
-__all__ = ["ShaftElement", "ShaftElement6DoF"]
+__all__ = ["ShaftElement"]
 bokeh_colors = bp.RdGy[11]
 
 
@@ -25,7 +25,7 @@ class ShaftElement(Element):
 
     .. math::
 
-        [x_0, y_0, \alpha_0, \beta_0, x_1, y_1, \alpha_1, \beta_1]**T
+        [x_0, y_0, \alpha_0, \beta_0, x_1, y_1, \alpha_1, \beta_1]^T
     Where :math:`\alpha_0` and :math:`\alpha_1` are the bending on the yz plane
     and :math:`\beta_0` and :math:`\beta_1` are the bending on the xz plane.
 
@@ -246,7 +246,7 @@ class ShaftElement(Element):
         self.beam_cg = L * c1 / (4 * c2)
         self.axial_cg_pos = None
 
-        # Slenderness ratio of beam elements (G*A*L**2) / (E*I)
+        # 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
 
@@ -1191,953 +1191,3 @@ class ShaftElement(Element):
         ]
 
         return elements
-
-
-class ShaftElement6DoF(ShaftElement):
-
-    r"""A 6 Degrees of Freedom shaft element.
-    =========================================
-    This class will create a shaft element that takes into
-    account shear stress, rotary inertia and gyroscopic effects.
-    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.
-
-    Parameters
-    ----------
-    L : float, pint.Quantity
-        Element length.
-    idl : float, pint.Quantity
-        Inner diameter of the element at the left node.
-    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.
-    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 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;
-        Default is zero.
-    torque : float, optional
-        Torque moment;
-        Default is zero.
-    shear_effects : bool, optional
-        Determine if shear effects are taken into account;
-        Default is True.
-    rotary_inertia : bool, optional
-        Determine if rotary_inertia effects are taken into account;
-        Default is True.
-    gyroscopic : bool, optional
-        Determine if gyroscopic effects are taken into account;
-        Default is True.
-    tag : str, optional
-        Element tag;
-        Default is None.
-
-    Returns
-    -------
-    A 6 degrees of freedom shaft element, with available gyroscopic, shear and rotary inertia effects.
-
-    Attributes
-    ----------
-    Poisson : float
-        Poisson coefficient for the element.
-    kappa : float
-        Shear coefficient for the element, determined from :cite:`Hutchingson2001` formulation.
-
-    References
-    ----------
-    .. bibliography:: ../../../docs/refs.bib
-
-    Examples
-    --------
-    >>> from ross.materials import steel
-    >>> eixo1 = 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)
-    >>> Euler_Bernoulli_Element = ShaftElement6DoF(
-    ...                         material=steel, L=0.5, idl=0.05, odl=0.1,
-    ...                         idr=0.05, odr=0.15,
-    ...                         alpha=0.01, beta=100,
-    ...                         rotary_inertia=False,
-    ...                         shear_effects=False)
-    >>> Euler_Bernoulli_Element.kappa
-    0.7099387976608923
-    """
-
-    @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:
-            # os.chdir(Path(os.path.dirname(ross.__file__)))
-            self.material = Material.use_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
-
-        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
-
-    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 save(self, file_name=Path(os.getcwd())):
-        """Save shaft elements to toml file.
-
-        Parameters
-        ----------
-        file_name : str
-
-        Examples
-        --------
-
-        """
-
-        data = self.get_data(Path(file_name) / "ShaftElement6DoF.toml")
-        data["ShaftElement6DoF"][str(self.n)] = {
-            "L": self.L,
-            "idl": self.idl,
-            "odl": self.odl,
-            "idr": self.idr,
-            "odr": self.odr,
-            "alpha": self.alpha,
-            "beta": self.beta,
-            "material": self.material.name,
-            "n": self.n,
-            "axial_force": self.axial_force,
-            "torque": self.torque,
-            "shear_effects": self.shear_effects,
-            "rotary_inertia": self.rotary_inertia,
-            "gyroscopic": self.gyroscopic,
-        }
-        self.dump_data(data, Path(file_name) / "ShaftElement6DoF.toml")
-
-    @staticmethod
-    def load(file_name="ShaftElement6DoF"):
-        """Load previously saved shaft elements from toml file.
-
-        Parameters
-        ----------
-        file_name : str, optional
-
-        Examples
-        --------
-
-        """
-        shaft_elements = []
-        with open("ShaftElement6DoF.toml", "r") as f:
-            shaft_elements_dict = toml.load(f)
-            for element in shaft_elements_dict["ShaftElement6DoF"]:
-                shaft_elements.append(
-                    ShaftElement6DoF(**shaft_elements_dict["ShaftElement6DoF"][element])
-                )
-        return shaft_elements
-
-    @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
-
-    def dof_mapping(self):
-        """
-        Method to map the element's 12 degrees of freedom
-
-        Parameters
-        ----------
-        Returns
-        -------
-        The numbering of the degrees of freedom for each element node. 
-        
-        Being the following their ordering for an element:
-
-        x_0,u_0 - horizontal translation, 
-        y_0,v_0 - vertical translation, 
-        z_0,w_0 - axial translation, 
-        theta_0 - rotation around horizontal, 
-        psi_0   - rotation around vertical, 
-        phi_0   - torsion around axial
-        x_1,u_1 - horizontal translation, 
-        y_1,v_1 - vertical translation, 
-        z_1,w_1 - axial translation, 
-        theta_1 - rotation around horizontal, 
-        psi_1   - rotation around vertical, 
-        phi_1   - torsion around axial
-        """
-        return dict(
-            u_0=0,
-            v_0=1,
-            w_0=2,
-            theta_0=3,
-            psi_0=4,
-            phi_0=5,
-            u_1=6,
-            v_1=7,
-            w_1=8,
-            theta_1=9,
-            psi_1=10,
-            phi_1=11,
-        )
-
-    def __repr__(self):
-        """This function returns a string representation of a 6 DoF shaft element.
-
-        Parameters
-        ----------
-        Returns
-        -------
-        A string representation of a 6 DoF shaft object.
-
-        Examples
-        --------
-        >>> from ross.materials import steel
-        >>> shaft1 = ShaftElement(
-        ...        L=0.25, idl=0, idr=0, odl=0.05, odr=0.08,
-        ...        material=steel, rotary_inertia=True, shear_effects=True
-        ... )
-        >>> shaft1 # doctest: +ELLIPSIS
-        ShaftElement(L=0.25, idl=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})"
-        )
-
-    def __str__(self):
-        """
-        Method to convert object into string
-        
-        Parameters
-        ----------
-        Returns
-        -------
-        The object's parameters translated to strings
-        """
-        return (
-            f"\nElem. N:    {self.n}"
-            f"\nLenght:     {self.L:{10}.{5}}"
-            f"\nLeft Int. Diam.: {self.idl:{10}.{5}}"
-            f"\nLeft Out. Diam.: {self.odl:{10}.{5}}"
-            f"\nRight Int. Diam.: {self.idr:{10}.{5}}"
-            f"\nRight Out. Diam.: {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}"
-            f"\n"
-        )
-
-    def M(self):
-        r"""Mass matrix for an instance of a 6 DoF shaft element.
-
-        Returns
-        -------
-        M: np.ndarray
-            Mass matrix for the 6 DoF shaft element.
-
-        Examples
-        --------
-        >>> Timoshenko_Element = ShaftElement6DoF(0.25, 0, 0.05, material=steel)
-        >>> Timoshenko_Element.M().shape
-        (12, 12)
-        """
-
-        # temporary material and geometrical constants
-        L = self.L
-        tempG = self.material.E / (2 * (1 + self.material.Poisson))
-        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
-            )
-        )
-
-        # element level matrix declaration
-
-        aux1 = self.material.rho * tempS * L / 420
-        aux2 = self.material.rho * tempJ * L / 6
-        # fmt: off
-
-        # 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],
-        ])
-
-        # 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],
-        ])
-
-        # Axial terms inertia matrix
-        Ma = self.material.rho * tempS * L / 6 * 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],
-            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-            [ 0, 0, 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, 0, 0, 0],
-            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-            [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-        ])
-
-        # Torsional terms inertias matrix
-        Mr = self.material.rho * tempI * L / 6 * 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, 2, 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, 2],
-        ])
-
-        # fmt: on
-
-        M = M + Ms + Ma + Mr
-
-        return M
-
-    def K(self):
-        r"""Stiffness matrix for an instance of a 6 DoF shaft element.
-
-        Returns
-        -------
-        K: np.ndarray
-            Omega independent stiffness matrix for the 6 DoF shaft element.
-
-        Examples
-        --------
-        >>> from ross.materials import steel
-        >>> Timoshenko_Element = ShaftElement6DoF(0.25, 0, 0.05, material=steel)
-        >>> 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 taperad 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
-            )
-        )
-
-        # temporary variables
-        A = (
-            12
-            * self.material.E
-            * tempI
-            / (self.material.G_s * self.kappa * tempS * L ** 2)
-        )
-
-        # auxiliary variables
-        a1 = self.material.E * tempI / ((1 + A) * L ** 3)
-        a2 = self.material.G_s * tempJ / L
-        a3 = self.material.E * tempS / L
-
-        # fmt: off
-        # pure stiffness matrix [Kc], added to the axial loads stiffness matrix [Ka],
-        # torsional stiffnesses matrix [Kr] and Tinoshenko 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],
-        ])
-
-        # 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],
-        ])
-
-        # 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],
-        ])
-
-        # fmt: on
-        # Dynamic stiffness matrix is added independently in "def Kst"
-        # Kst = self.material.rho*tempI/(15*L)*np.array(12,12)
-
-        K = Kc_plus + Kf + Kt
-
-        return K
-
-    def Kst(self):
-        r"""Dynamic stiffness matrix for an instance of a 6 DoF shaft element.
-
-        Returns
-        -------
-        Kst: np.ndarray
-            Dynamic stiffness matrix for the 6 DoF shaft element. This is
-            directly dependent on the rotation speed Omega. It needs to be
-            multiplied by the adequate Omega value when used in time depen-
-            dent analyses. The matrix multiplier term is: 
-            
-            [(Iz*Omega*rho)/(15*L)] * [Kst]
-
-            and here the Omega value has been suppressed and must be added
-            in the adequate analyses.
-
-        Examples
-        --------
-        >>> from ross.materials import steel
-        >>> Timoshenko_Element = ShaftElement6DoF(0.25, 0, 0.05, material=steel)
-        >>> 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
-            )
-        )
-
-        # 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],
-        ])
-        # fmt: on
-
-        return Kst
-
-    def C(self):
-        r"""Proportional damping matrix for an instance of a 6 DoF shaft element.
-
-        Returns
-        -------
-        C: np.ndarray
-            Proportional damping matrix for the 6 DoF shaft element.
-
-        Examples
-        --------
-        >>> from ross.materials import steel
-        >>> Timoshenko_Element = ShaftElement6DoF(L=0.25, idl=0, odl=0.05, material=steel)
-        >>> Timoshenko_Element.C().shape
-        (12, 12)
-        """
-
-        # proportional damping matrix
-        C = self.alpha * self.M() + self.beta * self.K()
-
-        return C
-
-    def G(self):
-        r"""Gyroscopic matrix for an instance of a 6 DoFs shaft element.
-
-        Returns
-        -------
-        G: np.ndarray
-            Gyroscopic matrix for the 6 DoF shaft element. Similar to the Kst
-            stiffness matrix, this Gyro matrix is also multiplied by the value
-            of the rotating speed Omega. It is omitted from this and must be 
-            added in the respective analyses.
-
-        Examples
-        --------
-        >>> from ross.materials import steel
-        >>> Timoshenko_Element = ShaftElement6DoF(0.25, 0, 0.05, material=steel)
-        >>> Timoshenko_Element.G().shape
-        (12, 12)
-        """
-
-        if self.gyroscopic:
-            # temporary material and geometrical constants, determined as mean values
-            # from the left and right radii of the tapered 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
-                )
-            )
-
-            # fmt: off
-            # Gyroscopic effect matrix
-            G = (self.material.rho * tempI / (15 * L)) * np.array([
-                [   0, -36, 0,  -3*L,      0, 0,    0,   36, 0,  -3*L,      0, 0],
-                [  36,   0, 0,     0,   -3*L, 0,  -36,    0, 0,     0,   -3*L, 0],
-                [   0,   0, 0,     0,      0, 0,    0,    0, 0,     0,      0, 0],
-                [ 3*L,   0, 0,     0,-4*L**2, 0, -3*L,    0, 0,     0,   L**2, 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],
-                [ -36,   0, 0,     0,    3*L, 0,   36,    0, 0,     0,    3*L, 0],
-                [   0,   0, 0,     0,      0, 0,    0,    0, 0,     0,      0, 0],
-                [ 3*L,   0, 0,     0,   L**2, 0, -3*L,    0, 0,     0,-4*L**2, 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],
-            ])
-            # fmt: on
-
-        return G
-
-    def bokeh_patch(self, position, check_sld, bk_ax):
-        """Shaft element patch.
-        Patch that will be used to draw the shaft element.
-        Parameters
-        ----------
-        position : float
-            Position in which the patch will be drawn.
-        check_sld : bool
-            If True, HoverTool displays only the slenderness ratio and color
-            the elements in yellow if slenderness ratio < 1.6
-        bk_ax : bokeh plotting axes, optional
-            Axes in which the plot will be drawn.
-        Returns
-        -------
-        hover : Bokeh HoverTool
-            Bokeh HoverTool axes
-        """
-        if check_sld is True and self.slenderness_ratio < 1.6:
-            bk_color = "yellow"
-            legend = "Shaft - Slenderness Ratio < 1.6"
-        else:
-            bk_color = self.material.color
-            legend = "Shaft"
-
-        # bokeh plot - plot the shaft
-        z_upper = [position, position, position + self.L, position + self.L]
-        y_upper = [self.idl / 2, self.odl / 2, self.odr / 2, self.idr / 2]
-        z_lower = [position, position, position + self.L, position + self.L]
-        y_lower = [-self.idl / 2, -self.odl / 2, -self.odr / 2, -self.idr / 2]
-
-        source = ColumnDataSource(
-            dict(
-                z_l=[z_lower],
-                y_l=[y_lower],
-                z_u=[z_upper],
-                y_u=[y_upper],
-                sld=[self.slenderness_ratio],
-                elnum=[self.n],
-                out_d_l=[self.odl],
-                out_d_r=[self.odr],
-                in_d_l=[self.idl],
-                in_d_r=[self.idr],
-                alpha_factor=[self.alpha],
-                beta_factor=[self.beta],
-                length=[self.L],
-                mat=[self.material.name],
-            )
-        )
-
-        bk_ax.patches(
-            xs="z_u",
-            ys="y_u",
-            source=source,
-            line_color=bokeh_colors[0],
-            line_width=1,
-            fill_alpha=0.5,
-            fill_color=bk_color,
-            legend_label=legend,
-            name="u_shaft",
-        )
-        bk_ax.patches(
-            xs="z_l",
-            ys="y_l",
-            source=source,
-            line_color=bokeh_colors[0],
-            line_width=1,
-            fill_alpha=0.5,
-            fill_color=bk_color,
-            legend_label=legend,
-            name="l_shaft",
-        )
-
-        hover = HoverTool(names=["l_shaft", "u_shaft"])
-        if check_sld:
-            hover.tooltips = [
-                ("Element Number :", "@elnum"),
-                ("Slenderness Ratio :", "@sld"),
-            ]
-        else:
-            hover.tooltips = [
-                ("Element Number :", "@elnum"),
-                ("Left Outer Diameter :", "@out_d_l"),
-                ("Left Inner Diameter :", "@in_d_l"),
-                ("Right Outer Diameter :", "@out_d_r"),
-                ("Right Inner Diameter :", "@in_d_r"),
-                ("Alpha Damp. Factor :", "@alpha_factor"),
-                ("Beta Damp. Factor :", "@beta_factor"),
-                ("Element Length :", "@length"),
-                ("Material :", "@mat"),
-            ]
-        hover.mode = "mouse"
-
-        return hover
-
-    @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(
-                le,
-                (s_idr - s_idl) * i * le / L + s_idl,
-                (s_odr - s_odl) * i * le / L + s_odl,
-                (s_idr - s_idl) * (i + 1) * le / L + s_idl,
-                (s_odr - s_odl) * (i + 1) * le / L + s_odl,
-                alpha,
-                beta,
-                material,
-                n,
-                gyroscopic,
-            )
-            for i in range(ne)
-        ]
-
-        return elements
diff --git a/ross/stochastic/st_results.py b/ross/stochastic/st_results.py
index 03b7eb5d..addf19e8 100644
--- a/ross/stochastic/st_results.py
+++ b/ross/stochastic/st_results.py
@@ -3,9 +3,11 @@
 This module returns graphs for each type of analyses in st_rotor_assembly.py.
 """
 import bokeh.palettes as bp
+import matplotlib.pyplot as plt
 import numpy as np
 from bokeh.layouts import gridplot
 from bokeh.plotting import figure
+from matplotlib import cm
 
 # set bokeh palette of colors
 colors1 = bp.Category10[10]
@@ -73,8 +75,8 @@ class ST_CampbellResults:
             kwargs.setdefault(k, v)
 
         fig = figure(
-            height=900,
-            width=900,
+            width=640,
+            height=480,
             tools="pan, box_zoom, wheel_zoom, reset, save",
             title="Campbell Diagram",
             y_axis_label="Damped Natural Frequencies",
@@ -172,8 +174,8 @@ class ST_CampbellResults:
             kwargs.setdefault(k, v)
 
         fig = figure(
-            height=900,
-            width=900,
+            width=640,
+            height=480,
             tools="pan, box_zoom, wheel_zoom, reset, save",
             title="Campbell Diagram",
             y_axis_label="Log Dec",
@@ -356,8 +358,8 @@ class ST_FrequencyResponseResults:
                     x=self.speed_range,
                     y=mag_percentile,
                     line_color=colors1[i],
-                    line_alpha=1.0,
-                    line_width=3,
+                    line_alpha=0.6,
+                    line_width=2.5,
                     muted_color=colors1[i],
                     muted_alpha=0.1,
                     legend_label="percentile: {}%".format(p),
@@ -371,8 +373,8 @@ class ST_FrequencyResponseResults:
                     x=self.speed_range,
                     y=mag_conf1,
                     line_color=colors1[i],
-                    line_alpha=1.0,
-                    line_width=3,
+                    line_alpha=0.6,
+                    line_width=2.5,
                     muted_color=colors1[i],
                     muted_alpha=0.1,
                     legend_label="confidence interval: {}%".format(p),
@@ -381,8 +383,8 @@ class ST_FrequencyResponseResults:
                     x=self.speed_range,
                     y=mag_conf2,
                     line_color=colors1[i],
-                    line_alpha=1.0,
-                    line_width=3,
+                    line_alpha=0.6,
+                    line_width=2.5,
                     muted_color=colors1[i],
                     muted_alpha=0.1,
                     legend_label="confidence interval: {}%".format(p),
@@ -449,8 +451,8 @@ class ST_FrequencyResponseResults:
                     x=self.speed_range,
                     y=phs_percentile,
                     line_color=colors1[i],
-                    line_alpha=1.0,
-                    line_width=3,
+                    line_alpha=0.6,
+                    line_width=2.5,
                     muted_color=colors1[i],
                     muted_alpha=0.1,
                     legend_label="percentile: {}%".format(p),
@@ -463,8 +465,8 @@ class ST_FrequencyResponseResults:
                 x=self.speed_range,
                 y=phs_conf1,
                 line_color=colors1[i],
-                line_alpha=1.0,
-                line_width=3,
+                line_alpha=0.6,
+                line_width=2.5,
                 muted_color=colors1[i],
                 muted_alpha=0.1,
                 legend_label="confidence interval: {}%".format(p),
@@ -473,8 +475,8 @@ class ST_FrequencyResponseResults:
                 x=self.speed_range,
                 y=phs_conf2,
                 line_color=colors1[i],
-                line_alpha=1.0,
-                line_width=3,
+                line_alpha=0.6,
+                line_width=2.5,
                 legend_label="confidence interval: {}%".format(p),
             )
 
@@ -533,7 +535,7 @@ class ST_FrequencyResponseResults:
 
 
 class ST_TimeResponseResults:
-    """Store stochastic results and provide plots for Time Response.
+    """Store stochastic results and provide plots for Time Response and Orbit Response.
 
     Parameters
     ----------
@@ -543,28 +545,37 @@ class ST_TimeResponseResults:
         System response.
     xout : array
         Time evolution of the state vector.
-    dof : int
-        Degree of freedom that will be observ
+    nodes_list: array
+        list with nodes from a rotor model.
+    nodes_pos: array
+        Rotor nodes axial positions.
 
     Returns
     -------
     fig : bokeh figure
         Returns the bokeh axes object with the plot
+    ax : matplotlib.axes
+        Matplotlib axes with orbit response plot.
     """
 
-    def __init__(self, time_range, yout, xout, dof):
+    def __init__(self, time_range, yout, xout, nodes_list=[], nodes_pos=[]):
         self.time_range = time_range
         self.yout = yout
         self.xout = xout
-        self.dof = dof
+        self.nodes_list = nodes_list
+        self.nodes_pos = nodes_pos
 
-    def plot(self, percentile=[], conf_interval=[], *args, **kwargs):
+    def _plot_time_response(
+        self, dof, percentile=[], conf_interval=[], *args, **kwargs
+    ):
         """Plot time response.
 
         This method plots the time response given.
 
         Parameters
         ----------
+        dof : int
+            Degree of freedom that will be observed.
         percentile : list, optional
             Sequence of percentiles to compute, which must be
             between 0 and 100 inclusive.
@@ -582,21 +593,21 @@ class ST_TimeResponseResults:
         grid_plots : bokeh figure
             Bokeh figure with time response plot.
         """
-        if self.dof % 4 == 0:
+        if dof % 4 == 0:
             obs_dof = "x"
-        elif self.dof % 4 == 1:
+        elif dof % 4 == 1:
             obs_dof = "y"
-        elif self.dof % 4 == 2:
+        elif dof % 4 == 2:
             obs_dof = "α"
-        elif self.dof % 4 == 3:
+        elif dof % 4 == 3:
             obs_dof = "β"
 
         fig = figure(
             tools="pan, box_zoom, wheel_zoom, reset, save",
-            width=900,
-            height=900,
+            width=640,
+            height=480,
             title="Response for node {} and degree of freedom {}".format(
-                self.dof // 4, obs_dof
+                dof // 4, obs_dof
             ),
             x_axis_label="Time (s)",
             y_axis_label="Amplitude",
@@ -606,29 +617,28 @@ class ST_TimeResponseResults:
         fig.axis.major_label_text_font_size = "14pt"
         fig.title.text_font_size = "14pt"
 
-        if len(percentile):
-            for i, p in enumerate(percentile):
-                phs_percentile = np.percentile(self.yout[..., self.dof], p, axis=0)
-                fig.line(
-                    x=self.time_range,
-                    y=phs_percentile,
-                    line_color=colors1[i],
-                    line_alpha=1.0,
-                    line_width=3,
-                    muted_color=colors1[i],
-                    muted_alpha=0.1,
-                    legend_label="percentile: {}%".format(p),
-                )
+        for i, p in enumerate(percentile):
+            tr_percentile = np.percentile(self.yout[..., dof], p, axis=0)
+            fig.line(
+                x=self.time_range,
+                y=tr_percentile,
+                line_color=colors1[i],
+                line_alpha=0.6,
+                line_width=2.5,
+                muted_color=colors1[i],
+                muted_alpha=0.1,
+                legend_label="percentile: {}%".format(p),
+            )
 
         for i, p in enumerate(conf_interval):
-            conf1 = np.percentile(self.yout[..., self.dof], 50 + p / 2, axis=0)
-            conf2 = np.percentile(self.yout[..., self.dof], 50 - p / 2, axis=0)
+            conf1 = np.percentile(self.yout[..., dof], 50 + p / 2, axis=0)
+            conf2 = np.percentile(self.yout[..., dof], 50 - p / 2, axis=0)
             fig.line(
                 x=self.time_range,
                 y=conf1,
                 line_color=colors1[i],
-                line_alpha=1.0,
-                line_width=3,
+                line_alpha=0.6,
+                line_width=2.5,
                 muted_color=colors1[i],
                 muted_alpha=0.1,
                 legend_label="confidence interval: {}%".format(p),
@@ -637,12 +647,12 @@ class ST_TimeResponseResults:
                 x=self.time_range,
                 y=conf2,
                 line_color=colors1[i],
-                line_alpha=1.0,
-                line_width=3,
+                line_alpha=0.6,
+                line_width=2.5,
                 legend_label="confidence interval: {}%".format(p),
             )
 
-        t_mean = np.mean(self.yout[..., self.dof], axis=0)
+        t_mean = np.mean(self.yout[..., dof], axis=0)
 
         fig.line(
             x=self.time_range,
@@ -660,6 +670,250 @@ class ST_TimeResponseResults:
 
         return fig
 
+    def _plot_orbit_2d(self, node, percentile=[], conf_interval=[], *args, **kwargs):
+        """Plot orbit response (2D).
+
+        This function plots orbits for a given node on the rotor system in a 2D view.
+
+        Parameters
+        ----------
+        node : int
+            Select the node to display the respective orbit response.
+        percentile : list, optional
+            Sequence of percentiles to compute, which must be
+            between 0 and 100 inclusive.
+        conf_interval : list, optional
+            Sequence of confidence intervals to compute, which must be
+            between 0 and 100 inclusive.
+        args : optional
+            Additional plot axes
+        kwargs : optional
+            Additional key word arguments can be passed to change
+            the plot (e.g. color="blue").
+
+        Returns
+        -------
+        ax : matplotlib.axes
+            Matplotlib axes with orbit response plot.
+        """
+        fig = figure(
+            tools="pan, box_zoom, wheel_zoom, reset, save",
+            width=640,
+            height=480,
+            title="Rotor Orbit: node {}".format(node),
+            x_axis_label="Amplitude",
+            y_axis_label="Amplitude",
+        )
+        fig.xaxis.axis_label_text_font_size = "14pt"
+        fig.yaxis.axis_label_text_font_size = "14pt"
+        fig.axis.major_label_text_font_size = "14pt"
+        fig.title.text_font_size = "14pt"
+
+        for i, p in enumerate(percentile):
+            fig.line(
+                np.percentile(self.yout[..., 4 * node], 50 + p / 2, axis=0),
+                np.percentile(self.yout[..., 4 * node + 1], 50 + p / 2, axis=0),
+                line_color=colors2[i],
+                line_alpha=0.6,
+                line_width=2.5,
+                muted_color=colors2[i],
+                muted_alpha=0.1,
+                legend_label="percentile: {}%".format(p),
+            )
+
+        for i, p in enumerate(conf_interval):
+            fig.line(
+                np.percentile(self.yout[..., 4 * node], 50 + p / 2, axis=0),
+                np.percentile(self.yout[..., 4 * node + 1], 50 + p / 2, axis=0),
+                line_color=colors1[i],
+                line_alpha=0.6,
+                line_width=2.5,
+                muted_color=colors1[i],
+                muted_alpha=0.1,
+                legend_label="confidence interval: {}%".format(p),
+            )
+            fig.line(
+                np.percentile(self.yout[..., 4 * node], 50 - p / 2, axis=0),
+                np.percentile(self.yout[..., 4 * node + 1], 50 - p / 2, axis=0),
+                line_color=colors1[i],
+                line_alpha=0.6,
+                line_width=2.5,
+                muted_color=colors1[i],
+                muted_alpha=0.1,
+                legend_label="confidence interval: {}%".format(p),
+            )
+
+        fig.line(
+            np.mean(self.yout[..., 4 * node], axis=0),
+            np.mean(self.yout[..., 4 * node + 1], axis=0),
+            line_color="black",
+            line_alpha=1.0,
+            line_width=3,
+            muted_color="black",
+            muted_alpha=0.1,
+            legend_label="Mean",
+        )
+        fig.legend.background_fill_alpha = 0.1
+        fig.legend.click_policy = "mute"
+        fig.legend.label_text_font_size = "10pt"
+
+        return fig
+
+    def _plot_orbit_3d(self, percentile=[], conf_interval=[], *args, **kwargs):
+        """Plot orbit response (3D).
+
+        This function plots orbits for each node on the rotor system in a 3D view.
+
+        Parameters
+        ----------
+        percentile : list, optional
+            Sequence of percentiles to compute, which must be
+            between 0 and 100 inclusive.
+        conf_interval : list, optional
+            Sequence of confidence intervals to compute, which must be
+            between 0 and 100 inclusive.
+        args : optional
+            Additional plot axes
+        kwargs : optional
+            Additional key word arguments can be passed to change
+            the plot (e.g. color="blue").
+
+        Returns
+        -------
+        ax : matplotlib.axes
+            Matplotlib axes with orbit response plot.
+        """
+        from mpl_toolkits.mplot3d import Axes3D
+
+        fig = plt.figure()
+        ax = fig.gca(projection="3d")
+
+        # plot center line
+        line = np.zeros(len(self.nodes_pos))
+        ax.plot(line, line, self.nodes_pos, "k-.", linewidth=1.5, zdir="x")
+
+        for i, p in enumerate(percentile):
+            for n in self.nodes_list:
+                ax.plot(
+                    np.percentile(self.yout[..., 4 * n], p, axis=0),
+                    np.percentile(self.yout[..., 4 * n + 1], p, axis=0),
+                    self.nodes_pos[n],
+                    color=cm.get_cmap("tab20c")(i),
+                    zdir="x",
+                    label="percentile: {}%".format(p),
+                )
+
+        for i, p in enumerate(conf_interval):
+            for n in self.nodes_list:
+                ax.plot(
+                    np.percentile(self.yout[..., 4 * n], 50 + p / 2, axis=0),
+                    np.percentile(self.yout[..., 4 * n + 1], 50 + p / 2, axis=0),
+                    self.nodes_pos[n],
+                    color=cm.get_cmap("tab10")(i),
+                    zdir="x",
+                    label="confidence interval: {}%".format(p),
+                )
+                ax.plot(
+                    np.percentile(self.yout[..., 4 * n], 50 - p / 2, axis=0),
+                    np.percentile(self.yout[..., 4 * n + 1], 50 - p / 2, axis=0),
+                    self.nodes_pos[n],
+                    color=cm.get_cmap("tab10")(i),
+                    zdir="x",
+                    label="confidence interval: {}%".format(p),
+                )
+
+        for n in self.nodes_list:
+            ax.plot(
+                np.mean(self.yout[..., 4 * n], axis=0),
+                np.mean(self.yout[..., 4 * n + 1], axis=0),
+                self.nodes_pos[n],
+                color="k",
+                zdir="x",
+                label="Mean".format(p),
+            )
+
+        ax.set_xlabel("Rotor length (m)", labelpad=16, fontsize=14)
+        ax.set_ylabel("Amplitude - X direction (m)", labelpad=16, fontsize=14)
+        ax.set_zlabel("Amplitude - Y direction (m)", labelpad=16, fontsize=14)
+        ax.set_title("Rotor Orbits", fontsize=16)
+        ax.tick_params(axis="both", which="major", labelsize=14)
+        ax.tick_params(axis="both", which="minor", labelsize=14)
+
+        return ax
+
+    def plot(
+        self,
+        plot_type="1d",
+        node=None,
+        dof=None,
+        percentile=[],
+        conf_interval=[],
+        *args,
+        **kwargs
+    ):
+        """Plot stochastic time or orbit response.
+
+        This function plots calls the auxiliary methods to plot the stochastic orbit
+        response. The plot type options are:
+            - 1d: plot time response for a given degree of freedom of a rotor system.
+            - 2d: plot orbit of a selected node of a rotor system.
+            - 3d: plot orbits for each node on the rotor system in a 3D view.
+
+        If plot_type = "1d": input a dof.
+        If plot_type = "2d": input a node.
+        if plot_type = "3d": no need to input a dof or node.
+
+        Parameters
+        ----------
+        plot_type : str, optional
+            Defines with type of plot to display.
+            Options are: "2d" or "3d"
+            Default is "3d".
+        node : int, optional
+            Select the node to display the respective orbit response.
+            Default is None.
+        dof : int
+            Degree of freedom that will be observed.
+        percentile : list, optional
+            Sequence of percentiles to compute, which must be
+            between 0 and 100 inclusive.
+        conf_interval : list, optional
+            Sequence of confidence intervals to compute, which must be
+            between 0 and 100 inclusive.
+        args : optional
+            Additional plot axes
+        kwargs : optional
+            Additional key word arguments can be passed to change
+            the plot (e.g. color="blue").
+
+        Raise
+        -----
+        ValueError
+            Error raised if no node is specified when plot_type = '2d'.
+        ValueError
+            Error raised if an odd string is specified to plot_type
+
+        Returns
+        -------
+        ax : matplotlib.axes
+            Matplotlib axes with orbit response plot.
+        """
+        if plot_type == "1d":
+            return self._plot_time_response(
+                dof, percentile, conf_interval, *args, **kwargs
+            )
+        elif plot_type == "2d":
+            if node not in self.nodes_list:
+                raise ValueError("Please insert a valid node.")
+            else:
+                return self._plot_orbit_2d(
+                    node, percentile, conf_interval, *args, **kwargs
+                )
+        elif plot_type == "3d":
+            return self._plot_orbit_3d(percentile, conf_interval, *args, **kwargs)
+        else:
+            raise ValueError("Plot type not supported. Choose between '2d' or '3d'.")
+
 
 class ST_ForcedResponseResults:
     """Store stochastic results and provide plots for Forced Response.
@@ -743,8 +997,8 @@ class ST_ForcedResponseResults:
                     x=self.frequency_range,
                     y=mag_percentile,
                     line_color=colors1[i],
-                    line_alpha=1.0,
-                    line_width=3,
+                    line_alpha=0.6,
+                    line_width=2.5,
                     muted_color=colors1[i],
                     muted_alpha=0.1,
                     legend_label="percentile: {}%".format(p),
@@ -758,8 +1012,8 @@ class ST_ForcedResponseResults:
                     x=self.frequency_range,
                     y=mag_conf1,
                     line_color=colors1[i],
-                    line_alpha=1.0,
-                    line_width=3,
+                    line_alpha=0.6,
+                    line_width=2.5,
                     muted_color=colors1[i],
                     muted_alpha=0.1,
                     legend_label="confidence interval: {}%".format(p),
@@ -768,8 +1022,8 @@ class ST_ForcedResponseResults:
                     x=self.frequency_range,
                     y=mag_conf2,
                     line_color=colors1[i],
-                    line_alpha=1.0,
-                    line_width=3,
+                    line_alpha=0.6,
+                    line_width=2.5,
                     muted_color=colors1[i],
                     muted_alpha=0.1,
                     legend_label="confidence interval: {}%".format(p),
@@ -838,8 +1092,8 @@ class ST_ForcedResponseResults:
                     x=self.frequency_range,
                     y=phs_percentile,
                     line_color=colors1[i],
-                    line_alpha=1.0,
-                    line_width=3,
+                    line_alpha=0.6,
+                    line_width=2.5,
                     muted_color=colors1[i],
                     muted_alpha=0.1,
                     legend_label="percentile: {}%".format(p),
@@ -852,8 +1106,8 @@ class ST_ForcedResponseResults:
                 x=self.frequency_range,
                 y=phs_conf1,
                 line_color=colors1[i],
-                line_alpha=1.0,
-                line_width=3,
+                line_alpha=0.6,
+                line_width=2.5,
                 muted_color=colors1[i],
                 muted_alpha=0.1,
                 legend_label="confidence interval: {}%".format(p),
@@ -862,8 +1116,10 @@ class ST_ForcedResponseResults:
                 x=self.frequency_range,
                 y=phs_conf2,
                 line_color=colors1[i],
-                line_alpha=1.0,
-                line_width=3,
+                line_alpha=0.6,
+                line_width=2.5,
+                muted_color=colors1[i],
+                muted_alpha=0.1,
                 legend_label="confidence interval: {}%".format(p),
             )
 
diff --git a/ross/stochastic/st_rotor_assembly.py b/ross/stochastic/st_rotor_assembly.py
index 3668c6f1..06b9e6df 100644
--- a/ross/stochastic/st_rotor_assembly.py
+++ b/ross/stochastic/st_rotor_assembly.py
@@ -201,13 +201,11 @@ class ST_Rotor(object):
         self.is_random = is_random
         self.attribute_dict = attribute_dict
 
-        # iter_break allows the iteration to stop at once to collect ndof
-        self.iter_break = True
-        ndof = list(iter(self))[0].ndof
-
         # common parameters
         self.RV_size = RV_size
-        self.ndof = ndof
+
+        # collect a series of attributes from a rotor instance
+        self._get_rotor_args()
 
     def __iter__(self):
         """Return an iterator for the container.
@@ -285,8 +283,37 @@ class ST_Rotor(object):
             raise KeyError("Object does not have parameter: {}.".format(key))
         self.attribute_dict[key] = value
 
+    def _get_rotor_args(self):
+        """Get relevant attributes from a rotor system.
+
+        This auxiliary funtion get some relevant attributes from a rotor system, such as
+        the nodes numbers, nodes positions and number of degrees of freedom, and add it
+        to the stochastic rotor as attribute. If an attribute is somehow afected by a
+        random variable, the function returns its mean.
+        """
+        self.iter_break = True
+        aux_rotor = list(iter(self))[0]
+
+        self.ndof = aux_rotor.ndof
+        self.nodes = aux_rotor.nodes
+
+        if "shaft_elements" in self.is_random:
+            if any(
+                "L" in sh.is_random
+                for sh in self.attribute_dict["shaft_elements"]
+                if isinstance(sh, ST_ShaftElement)
+            ):
+                nodes_pos_matrix = np.zeros(len(self.nodes), self.RV_size)
+                for i, rotor in enumerate(iter(self)):
+                    nodes_pos_matrix[:, i] = rotor.nodes_pos
+                self.nodes_pos = np.mean(nodes_pos_matrix, axis=1)
+            else:
+                self.nodes_pos = aux_rotor.nodes_pos
+        else:
+            self.nodes_pos = aux_rotor.nodes_pos
+
     @staticmethod
-    def get_args(idx, *args):
+    def _get_args(idx, *args):
         """Build new list of arguments from a random list of arguments.
 
         This funtion takes a list with random values or with lists of random
@@ -311,7 +338,7 @@ class ST_Rotor(object):
         >>> rotors = srs.st_rotor_example()
         >>> old_list = [1, 2, [3, 4], 5]
         >>> index = [0, 1]
-        >>> new_list = [rotors.get_args(idx, old_list) for idx in index]
+        >>> new_list = [rotors._get_args(idx, old_list) for idx in index]
         >>> new_list
         [[1, 2, 3, 5], [1, 2, 4, 5]]
         """
@@ -364,7 +391,7 @@ class ST_Rotor(object):
             arg = []
             for key, value in args_dict.items():
                 if key in is_random and key in self.is_random:
-                    arg.append(self.get_args(i, value))
+                    arg.append(self._get_args(i, value))
                 elif key in is_random and key not in self.is_random:
                     arg.append(value[i])
                 else:
@@ -514,15 +541,16 @@ class ST_Rotor(object):
             magnitude[:, i] = results.magnitude[inp, out, :]
             phase[:, i] = results.phase[inp, out, :]
 
-        results = ST_FrequencyResponseResults(speed_range, magnitude, phase,)
+        results = ST_FrequencyResponseResults(speed_range, magnitude, phase)
 
         return results
 
-    def run_time_response(self, speed, force, time_range, dof, ic=None):
+    def run_time_response(self, speed, force, time_range, ic=None):
         """Stochastic time response for multiples rotor systems.
 
         This function will take a rotor object and plot its time response
-        given a force and a time.
+        given a force and a time. This method displays the amplitude vs time or the
+        rotor orbits.
         The force and ic parameters can be passed as random variables.
 
         Parameters
@@ -537,8 +565,6 @@ class ST_Rotor(object):
             ST_Rotor.RV_size
         time_range : 1-dimensional array
             Time array.
-        dof : int
-            Degree of freedom that will be observed.
         ic : 1-dimensional array, 2-dimensional array, optional
             The initial conditions on the state vector (zero by default).
             Inputing a 2-dimensional array, the method considers the
@@ -568,12 +594,16 @@ class ST_Rotor(object):
         >>> F = np.zeros((size, ndof))
         >>> F[:, 4 * node] = 10 * np.cos(2 * t)
         >>> F[:, 4 * node + 1] = 10 * np.sin(2 * t)
-        >>> results = rotors.run_time_response(speed, F, t, dof)
+        >>> results = rotors.run_time_response(speed, F, t)
 
-        # Plotting Time Response with bokeh
+        # Plotting Time Response 1D, 2D and 3D
 
-        >>> results.plot(conf_interval=[90]) # doctest: +ELLIPSIS
+        >>> results.plot(plot_type="1d", dof=dof, conf_interval=[90]) # doctest: +ELLIPSIS
+        Figure...
+        >>> results.plot(plot_type="2d", node=node, conf_interval=[90]) # doctest: +ELLIPSIS
         Figure...
+        >>> results.plot(plot_type="3d", conf_interval=[90]) # doctest: +ELLIPSIS
+        <matplotlib...
         """
         t_size = len(time_range)
         RV_size = self.RV_size
@@ -600,7 +630,9 @@ class ST_Rotor(object):
                 yout[i] = y
                 i += 1
 
-        results = ST_TimeResponseResults(time_range, yout, xout, dof)
+        results = ST_TimeResponseResults(
+            time_range, yout, xout, self.nodes, self.nodes_pos
+        )
 
         return results
 


Test patch:
diff --git a/ross/tests/test_bearing_seal_element.py b/ross/tests/test_bearing_seal_element.py
index f3e684fe..ccc6e6e4 100644
--- a/ross/tests/test_bearing_seal_element.py
+++ b/ross/tests/test_bearing_seal_element.py
@@ -6,7 +6,6 @@ import pytest
 from numpy.testing import assert_allclose
 
 from ross.bearing_seal_element import (BallBearingElement, BearingElement,
-                                       BearingElement6DoF,
                                        RollerBearingElement)
 
 # fmt: on
@@ -288,42 +287,3 @@ def test_roller_bearing_element():
     assert_allclose(rollerbearing.K(0), K)
     assert_allclose(rollerbearing.C(0), C)
     assert_allclose(rollerbearing.G(), G)
-
-
-def test_bearing_6dof():
-    bearing_6dof = BearingElement6DoF(
-        n=0, kxx=1e6, kyy=0.8e6, kzz=1e5, cxx=2e2, cyy=1.5e2, czz=0.5e2
-    )
-    # fmt: off
-    K = np.array(
-        [[1000000.,       0.,       0.],
-         [      0.,  800000.,       0.],
-         [      0.,       0.,  100000.]])
-    C = np.array(
-        [[200.,   0.,   0.],
-         [  0., 150.,   0.],
-         [  0.,   0.,  50.]])
-    M = np.zeros((3, 3))
-    G = np.zeros((3, 3))
-    # fmt: on
-
-    assert_allclose(bearing_6dof.K(0), K, rtol=1e-3)
-    assert_allclose(bearing_6dof.C(0), C, rtol=1e-3)
-    assert_allclose(bearing_6dof.M(), M, rtol=1e-3)
-    assert_allclose(bearing_6dof.G(), G, rtol=1e-3)
-
-
-def test_bearing_6dof_equality():
-    bearing_6dof_0 = BearingElement6DoF(
-        n=0, kxx=1e6, kyy=0.8e6, kzz=1e5, cxx=2e2, cyy=1.5e2, czz=0.5e2
-    )
-    bearing_6dof_1 = BearingElement6DoF(
-        n=0, kxx=1e6, kyy=0.8e6, kzz=1e5, cxx=2e2, cyy=1.5e2, czz=0.5e2
-    )
-    bearing_6dof_2 = BearingElement6DoF(
-        n=0, kxx=2e6, kyy=0.8e6, kzz=1e5, cxx=2e2, cyy=1.5e2, czz=0.5e2
-    )
-
-    assert bearing_6dof_0 == bearing_6dof_1
-    assert not bearing_6dof_1 == bearing_6dof_2
-    assert not bearing_6dof_0 == bearing_6dof_2
diff --git a/ross/tests/test_disk_element.py b/ross/tests/test_disk_element.py
index c46a57aa..c1dcc4a1 100644
--- a/ross/tests/test_disk_element.py
+++ b/ross/tests/test_disk_element.py
@@ -4,7 +4,7 @@ import numpy as np
 import pytest
 from numpy.testing import assert_allclose, assert_almost_equal
 
-from ross.disk_element import DiskElement, DiskElement6DoF
+from ross.disk_element import DiskElement
 from ross.materials import steel
 
 
@@ -75,85 +75,3 @@ def test_from_table():
         assert_allclose(disks[1].m, 6.90999178227835)
         assert_allclose(disks[1].Ip, 0.0469996988106328, atol=1.6e-07)
         assert_allclose(disks[1].Id, 0.0249998397928898, atol=1.6e-07)
-
-
-@pytest.fixture
-def disk_6dof():
-    return DiskElement6DoF(0, 32.58973, 0.17809, 0.32956)
-
-
-def test_index_6dof(disk_6dof):
-    assert disk_6dof.dof_local_index().u_0 == 0
-    assert disk_6dof.dof_local_index().v_0 == 1
-    assert disk_6dof.dof_local_index().w_0 == 2
-    assert disk_6dof.dof_local_index().theta_0 == 3
-    assert disk_6dof.dof_local_index().psi_0 == 4
-    assert disk_6dof.dof_local_index().phi_0 == 5
-
-
-def test_mass_matrix_disk_6dof(disk_6dof):
-    # fmt: off
-    Md1 = np.array([[32.58973, 0.,       0.,       0.,      0.,      0.   ],
-                    [0.,       32.58973, 0.,       0.,      0.,      0.   ],
-                    [0.,       0.,       32.58973, 0.,      0.,      0.   ],
-                    [0.,       0.,       0.,       0.17809, 0.,      0.   ],
-                    [0.,       0.,       0.,       0.,      0.17809, 0.   ],
-                    [0.,       0.,       0.,       0,       0.,      0.32956]])
-    # fmt: on
-    assert_almost_equal(disk_6dof.M(), Md1, decimal=5)
-
-
-def test_gyroscopic_matrix_disk_6dof(disk_6dof):
-    # fmt: off
-    Gd1 = 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.32956, 0],
-            [ 0, 0, 0, -0.32956,   0,       0],
-            [ 0, 0, 0,  0,         0,       0]])
-    # fmt: on
-
-    assert_almost_equal(disk_6dof.G(), Gd1, decimal=5)
-
-
-@pytest.fixture
-def disk_from_geometry_6dof():
-    return DiskElement6DoF.from_geometry(0, steel, 0.07, 0.05, 0.28)
-
-
-def test_mass_matrix_disk1_6dof(disk_from_geometry_6dof):
-    # fmt: off
-    Md1 = np.array([[32.58973, 0.,       0.,       0.,      0.,      0.   ],
-                    [0.,       32.58973, 0.,       0.,      0.,      0.   ],
-                    [0.,       0.,       32.58973, 0.,      0.,      0.   ],
-                    [0.,       0.,       0.,       0.17809, 0.,      0.   ],
-                    [0.,       0.,       0.,       0.,      0.17809, 0.   ],
-                    [0.,       0.,       0.,       0,       0.,      0.32956]])
-    # fmt: on
-
-    M_class = disk_from_geometry_6dof.M()
-    assert_almost_equal(M_class, Md1, decimal=5)
-
-
-def test_gyroscopic_matrix_disk1_6dof(disk_from_geometry_6dof):
-    # fmt: off
-    Gd1 = 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.32956, 0],
-            [ 0, 0, 0, -0.32956,   0,       0],
-            [ 0, 0, 0,  0,         0,       0]])
-    # fmt: on
-    assert_almost_equal(disk_from_geometry_6dof.G(), Gd1, decimal=5)
-
-
-def test_from_table_6dof():
-    for file_name in ["/data/shaft_us.xls", "/data/shaft_si.xls"]:
-        file_name = os.path.dirname(os.path.realpath(__file__)) + file_name
-        disks = DiskElement6DoF.from_table(file_name, sheet_name="More")
-
-        assert_allclose(disks[1].m, 6.90999178227835)
-        assert_allclose(disks[1].Ip, 0.0469996988106328, atol=1.6e-07)
-        assert_allclose(disks[1].Id, 0.0249998397928898, atol=1.6e-07)
diff --git a/ross/tests/test_shaft_element.py b/ross/tests/test_shaft_element.py
index 22392e23..75c611a0 100644
--- a/ross/tests/test_shaft_element.py
+++ b/ross/tests/test_shaft_element.py
@@ -5,7 +5,7 @@ import pytest
 from numpy.testing import assert_allclose, assert_almost_equal
 
 from ross.materials import steel
-from ross.shaft_element import ShaftElement, ShaftElement6DoF
+from ross.shaft_element import ShaftElement
 
 test_dir = os.path.dirname(__file__)
 
@@ -306,121 +306,3 @@ def test_match_gyroscopic_matrix(tap2, tim2):
     G_tap = tap2.G()
     G_tim = tim2.G()
     assert_almost_equal(G_tap, G_tim, decimal=5)
-
-
-@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 test_s6_M_matrix(s6_eb):
-    # 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]])
-    # fmt: on
-    assert_allclose(s6_eb.M(), M_mat_6DoF, rtol=1e-3)
-
-
-def test_s6_G_matrix(s6_eb):
-    # fmt: off
-    G_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],
-         [  9.7239e+10,            0,            0,            0,  -8.1032e+08,            0,  -9.7239e+10,            0,            0,            0,  -8.1032e+08,            0],
-         [           0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
-         [  8.1032e+08,            0,            0,            0,  -1.0804e+08,            0,  -8.1032e+08,            0,            0,            0,   2.7011e+07,            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],
-         [ -9.7239e+10,            0,            0,            0,   8.1032e+08,            0,   9.7239e+10,            0,            0,            0,   8.1032e+08,            0],
-         [           0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0],
-         [  8.1032e+08,            0,            0,            0,   2.7011e+07,            0,  -8.1032e+08,            0,            0,            0,  -1.0804e+08,            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]])
-    # fmt: on
-    assert_allclose(s6_eb.G(), G_mat_6DoF, rtol=1e-3)
-
-
-def test_s6_K_matrix(s6_eb):
-    # 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]])
-
-    # fmt: on
-    assert_allclose(s6_eb.K(), K_mat_6DoF, rtol=1e-3)
-
-
-def test_s6_Kst_matrix(s6_eb):
-    # 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]])
-    # fmt: on
-    assert_allclose(s6_eb.Kst(), Kst_mat_6DoF, rtol=1e-3)
diff --git a/ross/tests/test_stochastic_rotor_assembly.py b/ross/tests/test_stochastic_rotor_assembly.py
index 8687df5f..ef8e5ab3 100644
--- a/ross/tests/test_stochastic_rotor_assembly.py
+++ b/ross/tests/test_stochastic_rotor_assembly.py
@@ -420,7 +420,6 @@ def test_time_response(rotor1):
     ndof = rotor1.ndof
     RV_size = rotor1.RV_size
     node = 3
-    dof = 9
     speed = 250.0
     t = np.linspace(0, 10, size)
     F = np.zeros((RV_size, size, ndof))
@@ -429,7 +428,7 @@ def test_time_response(rotor1):
     for i, f in enumerate(force):
         F[i, :, 4 * node] = f * np.cos(2 * t)
         F[i, :, 4 * node + 1] = f * np.sin(2 * t)
-    results = rotor1.run_time_response(speed, F, t, dof)
+    results = rotor1.run_time_response(speed, F, t)
 
     # fmt: off
     yout = np.array([