Repository: ross-rotordynamics/ross
Pull request #608 (2020-07-01 21:52:21)

Add file to configurate rotordynamic report
- Class "Config" generates an object with several preset parameters to run the rotordynamics analyses. It's a must to check all the options for a correct functioning. The attributes are automatically generated and it's not possible to remove then or add new ones.
- Class "_Dict" is a subclass used to organize nested dictionaries and set each key / value as attribute for the config object.

@raphaeltimbo and @GabrielBachmann I'd like to ask for some ideas for this file. I'm trying to follow a dictionary rule, just like Plotly does, but I couldn't find a pythonic way of how to update the config options. For now the user should change one by one, and it would be exhaustive. 
I'd like to introduce a way for the user input several options at the same time, like we do with Plotly. I've been taking a closer look on Plotly's code, but it seems too complex for what we need to do here.

Here is an example of how we could use `Config`:

```python
>>> options = Config()

# If we want to check the options, we just do:
>>> options
{'rotor_properties': {'rotor_speeds': {'min_speed': 1000.0, 'max_speed': 10000.0, 'oper_speed': 5000.0,
 'trip_speed': 12500.0, 'speed_factor ': 1.25, 'unit': 'rpm'}, 'rotor_id': {'type': 'compressor', 'tag': None}},
 'bearings': {'oper_clearance': None, 'min_clearance': None, 'max_clearance': None}, 'run_campbell': {'speed_range': None,
 'num_modes': 6}, 'plot_ucs': {'stiffness_range': None, 'num': 30, 'num_modes': 16, 'synchronous': False},
 'run_unbalance_response': {'frequency_range': None, 'modes': None, 'cluster_points': False, 'num_modes': 12, 'num_points': 10,
 'rtol': 0.005}, 'stability_level1': {'D': None, 'H': None, 'rated_power': None, 'rho_ratio': None, 'rho_suction': None,
 'rho_discharge': None, 'unit': 'm'}}

# I still need to work on a better repr or str method to display the options for the user

>>> # Displaying the options attributes. The dictionary syntax
>>> options["rotor_properties"]["rotor_id"]
{'type': 'compressor', 'tag': None}

>>> # Displaying the options attributes. The attribute syntax works too
>>> options.rotor_properties.rotor_id
{'type': 'compressor', 'tag': None}

>>> # Changing an option. This is the only syntax acceptable at this moment.
>>> options.rotor_properties.rotor_id.type = "turbine"
>>> options["rotor_properties"]["rotor_id"]
{'type': 'turbine', 'tag': None}

>>> # It's not possible to use the dictionary syntax so far:
>>> options["rotor_properties"]["rotor_id"]["type"] = "turbine" 
Traceback (most recent call last):

  File "<ipython-input-234-1d8cae20ceb8>", line 1, in <module>
    c["rotor_properties"]["rotor_id"]["type"] = "turbine"

TypeError: __setitem__() takes 1 positional argument but 3 were given
```

I wanted to add a `update_config` method to replace the `__setitem__` magic method.


Patch:
diff --git a/requirements.txt b/requirements.txt
index 6b0a1bda..adc57456 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -2,7 +2,7 @@
 numpy
 scipy
 matplotlib
-toml>=0.10.1
+toml
 pandas>=0.23
 plotly>=4.7.0
 bokeh>=1.4.0
diff --git a/ross/__init__.py b/ross/__init__.py
index bed038f3..3093ba63 100644
--- a/ross/__init__.py
+++ b/ross/__init__.py
@@ -1,9 +1,5 @@
 __version__ = "0.3.3"
 
-import plotly.io as _pio
-
-import ross.plotly_theme
-
 from .api_report import *
 from .bearing_seal_element import *
 from .disk_element import *
@@ -12,5 +8,3 @@ from .point_mass import *
 from .rotor_assembly import *
 from .shaft_element import *
 from .utils import visualize_matrix
-
-_pio.templates.default = "ross"
diff --git a/ross/api_report.py b/ross/api_report.py
index dc63f8f1..d82cd33c 100644
--- a/ross/api_report.py
+++ b/ross/api_report.py
@@ -1347,7 +1347,9 @@ class Report:
         )
 
         # Plotting area
-
+        default_layout = dict(
+            width=1200, height=900, plot_bgcolor="white", hoverlabel_align="right"
+        )
         fig1 = go.Figure()
 
         fig1.add_trace(
@@ -1355,6 +1357,7 @@ class Report:
                 x=cross_coupled_Qa,
                 y=log_dec,
                 mode="lines",
+                line=dict(width=3, color=colors1[5]),
                 showlegend=False,
                 hoverinfo="none",
             )
@@ -1364,6 +1367,7 @@ class Report:
                 x=[Qa],
                 y=[log_dec_a],
                 mode="markers",
+                marker=dict(size=10, color=colors1[5]),
                 name="<b>Qa: Anticipated cross-coupling</b>",
                 hoverinfo="none",
             )
@@ -1376,20 +1380,46 @@ class Report:
             xshift=15,
             yshift=15,
             text="<b>Qa</b>",
+            font=dict(size=14),
             showarrow=False,
         )
         fig1.update_xaxes(
             title_text="<b>Applied Cross-Coupled Stiffness, Q (N/m)</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            rangemode="nonnegative",
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig1.update_yaxes(
+            title_text="<b>Log Dec</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
             rangemode="nonnegative",
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
         )
-        fig1.update_yaxes(title_text="<b>Log Dec</b>", rangemode="nonnegative")
         fig1.update_layout(
             title=dict(
                 text=(
                     "<b>Applied Cross-Coupled Stiffness vs. Log Decrement</b><br>"
                     + "<b>(API 684 - SP 6.8.5.10)</b>"
-                )
-            )
+                ),
+                font=dict(size=20),
+            ),
+            legend=dict(
+                font=dict(family="sans-serif", size=14),
+                bgcolor="white",
+                bordercolor="black",
+                borderwidth=2,
+            ),
+            **default_layout,
         )
 
         fig2 = go.Figure()
@@ -1401,6 +1431,7 @@ class Report:
             xshift=40,
             yshift=0,
             text="<b>{}</b>".format(self.tag),
+            font=dict(size=14),
             showarrow=False,
         )
 
@@ -1413,6 +1444,7 @@ class Report:
                 xshift=0,
                 yshift=0,
                 text=f"<b>{text}</b>",
+                font=dict(size=14),
                 opacity=0.4,
                 showarrow=False,
             )
@@ -1422,6 +1454,7 @@ class Report:
                 x=RHO,
                 y=CSR_boundary,
                 mode="lines",
+                line=dict(width=5, color=colors1[5]),
                 showlegend=False,
                 hoverinfo="none",
                 xaxis="x",
@@ -1432,6 +1465,7 @@ class Report:
                 x=0.062428 * RHO,
                 y=CSR_boundary,
                 mode="lines",
+                line=dict(width=5, color=colors1[5]),
                 showlegend=False,
                 hoverinfo="none",
                 xaxis="x2",
@@ -1442,6 +1476,7 @@ class Report:
                 x=[RHO_mean],
                 y=[CSR],
                 mode="markers",
+                marker=dict(size=10, color=colors1[5]),
                 name="<b>CSR: Critical Speed Ratio</b>",
                 hoverinfo="none",
                 xaxis="x",
@@ -1451,19 +1486,38 @@ class Report:
         fig2.update_xaxes(mirror=True)
         fig2.update_yaxes(
             title_text="<b>Maximum Critical Speed Ratio</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
             rangemode="nonnegative",
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
             domain=[0.1, 1],
         )
         fig2.update_layout(
             xaxis=dict(
                 title_text="<b>kg/m³</b>",
+                title_font=dict(family="Arial", size=20),
+                tickfont=dict(size=16),
                 rangemode="nonnegative",
+                gridcolor="lightgray",
+                showline=True,
+                linewidth=2.5,
+                linecolor="black",
                 overlaying="x2",
                 anchor="y",
             ),
             xaxis2=dict(
                 title_text="<b>lb/ft³</b>",
+                title_font=dict(family="Arial", size=20),
+                tickfont=dict(size=16),
                 rangemode="nonnegative",
+                showgrid=False,
+                showline=True,
+                linewidth=2.5,
+                linecolor="black",
                 anchor="free",
                 side="bottom",
                 position=0,
@@ -1472,8 +1526,16 @@ class Report:
                 text=(
                     "<b>CSR vs. Mean Gas Density</b><br>"
                     + "<b>(API 684 - SP 6.8.5.10)</b>"
-                )
+                ),
+                font=dict(size=20),
+            ),
+            legend=dict(
+                font=dict(family="sans-serif", size=14),
+                bgcolor="white",
+                bordercolor="black",
+                borderwidth=2,
             ),
+            **default_layout,
         )
 
         # Level 1 screening criteria - API 684 - SP6.8.5.10
diff --git a/ross/bearing_seal_element.py b/ross/bearing_seal_element.py
index c53715aa..a4321b2e 100644
--- a/ross/bearing_seal_element.py
+++ b/ross/bearing_seal_element.py
@@ -12,12 +12,11 @@ from pathlib import Path
 import numpy as np
 import plotly.graph_objects as go
 import scipy.interpolate as interpolate
-import toml
 
 from ross.element import Element
 from ross.fluid_flow import fluid_flow as flow
 from ross.fluid_flow.fluid_flow_coefficients import (
-    calculate_short_damping_matrix, calculate_short_stiffness_matrix)
+    calculate_damping_matrix, calculate_stiffness_matrix)
 from ross.units import check_units
 from ross.utils import read_table_file
 
@@ -57,8 +56,8 @@ class _Coefficient:
     Examples
     --------
     >>> bearing = bearing_example()
-    >>> bearing.kxx # doctest: +ELLIPSIS
-    [1000000.0...
+    >>> bearing.kxx
+    [1000000.0, 1000000.0...
     """
 
     def __init__(self, coefficient, frequency=None):
@@ -429,10 +428,8 @@ class BearingElement(Element):
         self.n_link = n_link
         self.n_l = n
         self.n_r = n
-        if frequency is not None:
-            self.frequency = np.array(frequency, dtype=np.float64)
-        else:
-            self.frequency = frequency
+
+        self.frequency = np.array(frequency, dtype=np.float64)
         self.tag = tag
         self.color = color
         self.scale_factor = scale_factor
@@ -508,35 +505,114 @@ class BearingElement(Element):
     def __hash__(self):
         return hash(self.tag)
 
-    def save(self, file):
-        try:
-            data = toml.load(file)
-        except FileNotFoundError:
-            data = {}
+    def save(self, file_name=Path(os.getcwd())):
+        """Save 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) / "BearingElement.toml")
 
-        args = {
+        if type(self.frequency) == np.ndarray:
+            try:
+                self.frequency[0]
+                frequency = list(self.frequency)
+            except IndexError:
+                frequency = None
+
+        data["BearingElement"][str(self.n)] = {
             "n": self.n,
-            "kxx": [float(i) for i in self.kxx.coefficient],
-            "cxx": [float(i) for i in self.cxx.coefficient],
-            "kyy": [float(i) for i in self.kyy.coefficient],
-            "kxy": [float(i) for i in self.kxy.coefficient],
-            "kyx": [float(i) for i in self.kyx.coefficient],
-            "cyy": [float(i) for i in self.cyy.coefficient],
-            "cxy": [float(i) for i in self.cxy.coefficient],
-            "cyx": [float(i) for i in self.cyx.coefficient],
+            "kxx": self.kxx.coefficient,
+            "cxx": self.cxx.coefficient,
+            "kyy": self.kyy.coefficient,
+            "kxy": self.kxy.coefficient,
+            "kyx": self.kyx.coefficient,
+            "cyy": self.cyy.coefficient,
+            "cxy": self.cxy.coefficient,
+            "cyx": self.cyx.coefficient,
+            "frequency": frequency,
             "tag": self.tag,
             "n_link": self.n_link,
             "scale_factor": self.scale_factor,
         }
-        if self.frequency is not None:
-            args["frequency"] = [float(i) for i in self.frequency]
-        else:
-            args["frequency"] = self.frequency
+        self.dump_data(data, Path(file_name) / "BearingElement.toml")
+
+    @staticmethod
+    def load(file_name=""):
+        """Load a list of bearing elements saved in a toml format.
+
+        It works as an auxiliary function of the load function in the Rotor class.
+
+        Parameters
+        ----------
+        file_name: string
+            The name of the file of the bearing element to be loaded.
+
+        Returns
+        -------
+        A list of bearing elements.
+
+        Examples
+        --------
+        >>> bearing1 = bearing_example()
+        >>> bearing1.save(os.getcwd())
+        >>> list_of_bearings = BearingElement.load(os.getcwd())
+        >>> bearing1 == list_of_bearings[0]
+        True
+        """
+        bearing_elements = []
+        bearing_elements_dict = BearingElement.get_data(
+            file_name=Path(file_name) / "BearingElement.toml"
+        )
+        for element in bearing_elements_dict["BearingElement"]:
+            bearing = BearingElement(**bearing_elements_dict["BearingElement"][element])
+            bearing.kxx.coefficient = bearing_elements_dict["BearingElement"][element][
+                "kxx"
+            ]
+
+            bearing.kxy.coefficient = bearing_elements_dict["BearingElement"][element][
+                "kxy"
+            ]
+
+            bearing.kyx.coefficient = bearing_elements_dict["BearingElement"][element][
+                "kyx"
+            ]
+
+            bearing.kyy.coefficient = bearing_elements_dict["BearingElement"][element][
+                "kyy"
+            ]
+
+            bearing.cxx.coefficient = bearing_elements_dict["BearingElement"][element][
+                "cxx"
+            ]
 
-        data[f"{self.__class__.__name__}_{self.tag}"] = args
+            bearing.cxy.coefficient = bearing_elements_dict["BearingElement"][element][
+                "cxy"
+            ]
 
-        with open(file, "w") as f:
-            toml.dump(data, f)
+            bearing.cyx.coefficient = bearing_elements_dict["BearingElement"][element][
+                "cyx"
+            ]
+
+            bearing.cyy.coefficient = bearing_elements_dict["BearingElement"][element][
+                "cyy"
+            ]
+
+            bearing_elements.append(bearing)
+        return bearing_elements
 
     def dof_mapping(self):
         """Degrees of freedom mapping.
@@ -1021,8 +1097,8 @@ class BearingElement(Element):
             eccentricity=eccentricity,
             load=load,
         )
-        c = calculate_short_damping_matrix(fluid_flow)
-        k = calculate_short_stiffness_matrix(fluid_flow)
+        c = calculate_damping_matrix(fluid_flow, force_type="short")
+        k = calculate_stiffness_matrix(fluid_flow, force_type="short")
         return cls(
             n,
             kxx=k[0],
@@ -1738,37 +1814,101 @@ class BearingElement6DoF(BearingElement):
             )
         return False
 
-    def save(self, file):
-        try:
-            data = toml.load(file)
-        except FileNotFoundError:
-            data = {}
+    def save(self, file_name=Path(os.getcwd())):
+        """Save 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.
+
+        Examples
+        --------
+        >>> bearing = bearing_6dof_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
 
-        args = {
+        data["BearingElement6DoF"][str(self.n)] = {
             "n": self.n,
-            "kxx": [float(i) for i in self.kxx.coefficient],
-            "cxx": [float(i) for i in self.cxx.coefficient],
-            "kyy": [float(i) for i in self.kyy.coefficient],
-            "kxy": [float(i) for i in self.kxy.coefficient],
-            "kyx": [float(i) for i in self.kyx.coefficient],
-            "kzz": [float(i) for i in self.kzz.coefficient],
-            "cyy": [float(i) for i in self.cyy.coefficient],
-            "cxy": [float(i) for i in self.cxy.coefficient],
-            "cyx": [float(i) for i in self.cyx.coefficient],
-            "czz": [float(i) for i in self.czz.coefficient],
+            "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,
             "n_link": self.n_link,
             "scale_factor": self.scale_factor,
         }
-        if self.frequency is not None:
-            args["frequency"] = [float(i) for i in self.frequency]
-        else:
-            args["frequency"] = self.frequency
+        self.dump_data(data, Path(file_name) / "BearingElement6DoF.toml")
+
+    @staticmethod
+    def load(file_name=""):
+        """Load 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
+        --------
+        >>> bearing1 = bearing_6dof_example()
+        >>> bearing1.save(os.getcwd())
+        >>> list_of_bearings = BearingElement6DoF.load(os.getcwd())
+        >>> bearing1 == list_of_bearings[0]
+        True
+        """
+        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"]
 
-        data[f"{self.__class__.__name__}_{self.tag}"] = args
+            bearing.czz.coefficient = data[element]["czz"]
 
-        with open(file, "w") as f:
-            toml.dump(data, f)
+            bearing_elements.append(bearing)
+        return bearing_elements
 
     def dof_mapping(self):
         """Degrees of freedom mapping.
diff --git a/ross/config.py b/ross/config.py
new file mode 100644
index 00000000..f0024716
--- /dev/null
+++ b/ross/config.py
@@ -0,0 +1,471 @@
+"""Rotordynamic Report Configuration File."""
+import black
+
+__all__ = ["Config"]
+
+
+class _Dict:
+    """Set keys and values as attribute for the Config object.
+
+    Subclass to organize nested dictionaries and set each key / value as attribute for
+    the config object.
+
+    Return
+    ------
+    A dictionary as attribute for the Config() object.
+
+    Examples
+    --------
+    >>> param = _Dict({
+    ...     "stiffness_range": None,
+    ...     "num": 30,
+    ...     "num_modes": 16,
+    ...     "synchronous": False,
+    ... })
+    >>> param.num
+    30
+    """
+
+    def __init__(self, dictionary):
+        for k, v in dictionary.items():
+            if isinstance(v, dict):
+                setattr(self, k, self.__class__(v))
+            else:
+                setattr(self, k, v)
+
+    def __repr__(self):
+        """Return a string representation for the dict attribute.
+
+        This method uses Black code formatting for better visualization.
+
+        Returns
+        -------
+        A string representation for the dictionary options.
+
+        Examples
+        --------
+        >>> param = _Dict({
+        ...     "oper_clearance": None,
+        ...     "min_clearance": None,
+        ...     "max_clearance": None,
+        ... })
+        >>> param # doctest: +ELLIPSIS
+        {"oper_clearance": None, "min_clearance": None, "max_clearance": None}...
+        """
+        return black.format_file_contents(
+            repr(self.__dict__), fast=True, mode=black.FileMode()
+        )
+
+    def __getitem__(self, option):
+        """Return the value for a given option from the dictionary.
+
+        Parameters
+        ----------
+        option : str
+            A dictionary key corresponding to config options as string.
+
+        Raises
+        ------
+        KeyError
+            Raises an error if the parameter doesn't belong to the dictionary.
+
+        Returns
+        -------
+        Return the value for the given key.
+
+        Examples
+        --------
+        >>> param = _Dict({
+        ...     "stiffness_range": None,
+        ...     "num": 30,
+        ...     "num_modes": 16,
+        ...     "synchronous": False,
+        ... })
+        >>> param["num"]
+        30
+        """
+        if option not in self.__dict__.keys():
+            raise KeyError("Option '{}' not found.".format(option))
+
+        return self.__dict__[option]
+
+    def _update(self, **kwargs):
+        """Update the dict values.
+
+        This is an axuliar method for Config.update_config() to set new values for the
+        config dictionary according to kwargs input.
+        The kwargs must respect Config attributes. It's only possible to update
+        existing values from Config dictionary.
+        **See Config attributes reference for infos about the dict options.
+
+        Parameters
+        ----------
+        **kwargs : dict
+            Dictionary with new values for corresponding keys.
+
+        Raises
+        ------
+        KeyError
+            Raises an error if the parameter doesn't belong to the dictionary.
+
+        Examples
+        --------
+        >>> param = _Dict({
+        ...     "stiffness_range": None,
+        ...     "num": 30,
+        ...     "num_modes": 16,
+        ...     "synchronous": False,
+        ... })
+        >>> param._update(num=20, num_modes=10)
+        >>> param # doctest: +ELLIPSIS
+        {"stiffness_range": None, "num": 20, "num_modes": 10, "synchronous": False}...
+        """
+        for k, v in kwargs.items():
+            if k not in self.__dict__.keys():
+                raise KeyError("Option '{}' not found.".format(k))
+
+            if isinstance(v, dict):
+                getattr(self, k)._update(**v)
+            else:
+                self.__dict__[k] = v
+
+
+class Config:
+    """Configurate parameters for rotordynamic report.
+
+    This class generates an object with several preset parameters to run the
+    rotordynamics analyses. It's a must to check all the options for a correct
+    functioning.
+
+    The attributes are automatically generated and it's not possible to remove then or
+    add new ones.
+
+    Attributes
+    ----------
+    rotor_properties : dict
+        Dictionary of rotor properties.
+
+        rotor_speeds : dict
+            Dictionary indicating the operational speeds which will be used in the
+            analyses.
+
+            min_speed : float
+                The machine minimum operational speed.
+            max_speed : float
+                The machine  maximum operational speed.
+            oper_speed : float
+                The machine  nominal operational speed.
+            trip_speed : float
+                The machine overspeed trip.
+            speed_factor : float
+                Multiplicative factor of the speed range - according to API 684.
+                Default is 1.50.
+            unit : str
+                Unit system to speed values. Options: "rad/s", "rpm".
+                Default is "rpm".
+
+        rotor_id : dict
+            Dictionary with rotor identifications.
+
+            type : str
+                Machine type: Options are: "compressor", "turbine", "axial_flow". Each
+                options has it's own considerations (according to API 684). If a different
+                option is input, it will be the software will treat as a "compressor".
+            tag : str
+                Tag for the rotor. If None, it'll copy the tag from rotor object.
+
+    bearings : dict
+        The analyses consider different configurations for the bearings. It should be done
+        for minimum, maximum and the nominal clearance.
+
+        oper_clearance : list
+            List of bearing elements. The coefficients should be calculated for the nominal
+            clearance.
+        min_clearance : list
+            List of bearing elements. The coefficients should be calculated for the minimum
+            clearance.
+        max_clearance : list
+            List of bearing elements. The coefficients should be calculated for the maximum
+            clearance.
+
+    run_campbell : dict
+        Dictionary configurating run_campbell parameters.
+
+        speed_range : list, array
+            Array with the speed range.
+        num_modes : float
+            Number of frequencies that will be calculated.
+            Default is 6.
+
+    plot_ucs : dict
+        Dictionary configurating plot_ucs parameters.
+
+        stiffness_range : tuple, optional
+            Tuple with (start, end) for stiffness range.
+        num : int
+            Number of steps in the range.
+            Default is 30.
+        num_modes : int, optional
+            Number of modes to be calculated.
+            Default is 16.
+        synchronous : bool
+            If True a synchronous analysis is carried out and the frequency of
+            the first forward model will be equal to the speed.
+            Default is False.
+
+    run_unbalance_response : dict
+        Dictionary configurating run_unbalance_response parameters.
+
+        probes : dict
+            Dictionary with the node where the probe is set and its respective
+            orientation angle.
+
+            node : list
+                List with the nodes where probes are located.
+            orientation : list
+                List with the respective orientation angle for the probes.
+                0 or π (rad) corresponds to the X orientation and
+                π / 2 or 3 * π / 2 (rad) corresponds to the Y orientation.
+            unit : str
+                Unit system for the orientation angle. Can be "rad" or "degree".
+                Default is "rad".
+
+        frequency_range : list, array
+            Array with the desired range of frequencies. If None and cluster_points is
+            False, it creates an array from 0 to the max continuos speed times the
+            speed_factor.
+            If None and cluster_points is True, it creates and automatic array based on
+            the number of modes selected.
+            Default is None with cluster_points False.
+        modes : list, optional
+            Modes that will be used to calculate the frequency response
+            (all modes will be used if a list is not given).
+        cluster_points : bool, optional
+            Boolean to activate the automatic frequency spacing method. If True, the
+            method uses _clustering_points() to create an speed_range.
+            Default is False
+        num_points : int, optional
+            The number of points generated per critical speed.
+            The method set the same number of points for slightly less and slightly
+            higher than the natural circular frequency. It means there'll be num_points
+            greater and num_points smaller than a given critical speed.
+            num_points may be between 2 and 12. Anything above this range defaults
+            to 10 and anything below this range defaults to 4.
+            The default is 10.
+        num_modes
+            The number of eigenvalues and eigenvectors to be calculated using ARPACK.
+            It also defines the range for the output array, since the method generates
+            points only for the critical speed calculated by Rotor.run_critical_speed().
+            Default is 12.
+        rtol : float, optional
+            Tolerance (relative) for termination. Applied to scipy.optimize.newton to
+            calculate the approximated critical speeds.
+            Default is 0.005 (0.5%).
+
+    stability_level1 : dict
+        Dictionary configurating stability_level_1 parameters.
+
+        D : list, array
+            Impeller diameter, m (in.) or Blade pitch diameter, m (in.).
+            The disk elements order must be observed to input this list.
+        H : list, array
+            Minimum diffuser width per impeller, m (in.) or Effective blade height, m (in.).
+            The disk elements order must be observed to input this list.
+        rated_power : list
+            Rated power per stage/impeller, W (HP),
+            The disk elements order must be observed to input this list.
+        rho_ratio : list
+            Density ratio between the discharge gas density and the suction
+            gas density per impeller, kg/m3 (lbm/in.3),
+            The disk elements order must be observed to input this list.
+        rho_suction : float
+            Suction gas density in the first stage, kg/m3 (lbm/in.3).
+        rho_discharge : float
+            Discharge gas density in the last stage, kg/m3 (lbm/in.3),
+        unit: str, optional
+            Unit system. Options are "m" (meter) and "in" (inch).
+            Default is "m".
+
+    Returns
+    -------
+    A config object to rotordynamic report.
+
+    Examples
+    --------
+    There are two possible syntax to return the options setup. One is using the
+    object.attribute syntax and the other is the dicionary syntax
+
+    First syntax opion:
+    >>> configs = Config()
+    >>> configs.rotor_properties.rotor_id # doctest: +ELLIPSIS
+    {"type": "compressor", "tag": None}...
+
+    Second syntax opion:
+    >>> configs = Config()
+    >>> configs["rotor_properties"]["rotor_id"] # doctest: +ELLIPSIS
+    {"type": "compressor", "tag": None}...
+    """
+
+    def __init__(self):
+        # fmt: off
+        # Configurating rotor properties
+        self.rotor_properties = _Dict({
+            "rotor_speeds": {
+                "min_speed": None,
+                "max_speed": None,
+                "oper_speed": None,
+                "trip_speed": None,
+                "speed_factor": 1.50,
+                "unit": "rpm",
+            },
+            "rotor_id": {
+                "type": "compressor",
+                "tag": None
+            },
+        })
+
+        # Configurating bearing elements for diferent clearances
+        self.bearings = _Dict({
+            "oper_clearance": None,
+            "min_clearance": None,
+            "max_clearance": None,
+        })
+
+        # Configurating campbell options
+        self.run_campbell = _Dict({
+            "speed_range": None,
+            "num_modes": 6,
+        })
+
+        # Configurating UCS options
+        self.plot_ucs = _Dict({
+            "stiffness_range": None,
+            "num": 30,
+            "num_modes": 16,
+            "synchronous": False,
+        })
+
+        # Configurating unbalance response options
+        self.run_unbalance_response = _Dict({
+            "probes": {
+                "node": None,
+                "orientation": None,
+                "unit": "rad",
+            },
+            "frequency_range": None,
+            "modes": None,
+            "cluster_points": False,
+            "num_modes": 12,
+            "num_points": 10,
+            "rtol": 0.005,
+        })
+
+        # Configurating stability level 1 analysis
+        self.stability_level1 = _Dict({
+            "D": None,
+            "H": None,
+            "rated_power": None,
+            "rho_ratio": None,
+            "rho_suction": None,
+            "rho_discharge": None,
+            "unit": "m",
+        })
+        # fmt: on
+
+    def __repr__(self):
+        """Return a string representation for the config options.
+
+        This method uses Black code formatting for better visualization.
+
+        Returns
+        -------
+        A string representation for the config dictionary.
+
+        Examples
+        --------
+        >>> configs = Config()
+        >>> configs # doctest: +ELLIPSIS
+        {
+            "rotor_properties": {
+                "rotor_speeds": {
+                    "min_speed": None,
+                    "max_speed": None,
+                    "oper_speed": None...
+        """
+        return black.format_file_contents(
+            repr(self.__dict__), fast=True, mode=black.FileMode()
+        )
+
+    def __getitem__(self, option):
+        """Return the value for a given option from config dictionary.
+
+        Parameters
+        ----------
+        option : str
+            A dictionary key corresponding to config options as string.
+
+        Raises
+        ------
+        KeyError
+            Raises an error if the parameter doesn't belong to the dictionary.
+
+        Returns
+        -------
+        Return the value for the given key.
+
+        Examples
+        --------
+        >>> configs = Config()
+        >>> configs["bearings"] # doctest: +ELLIPSIS
+        {"oper_clearance": None, "min_clearance": None, "max_clearance": None}...
+        """
+        if option not in self.__dict__.keys():
+            raise KeyError("Option '{}' not found.".format(option))
+
+        return self.__dict__[option]
+
+    def update_config(self, **kwargs):
+        """Update the config options.
+
+        This method set new values for the config dictionary according to kwargs input.
+        The kwargs must respect Config attributes. It's only possible to update
+        existing values from Config dictionary.
+        **See Config attributes reference for infos about the dict options.
+
+        Parameters
+        ----------
+        **kwargs : dict
+            Dictionary with new values for corresponding keys.
+
+        Raises
+        ------
+        KeyError
+            Raises an error if the parameter doesn't belong to the dictionary.
+
+        Examples
+        --------
+        >>> configs = Config()
+        >>> configs.update_config(
+        ...     rotor_properties=dict(
+        ...         rotor_speeds=dict(min_speed=1000.0, max_speed=10000.0),
+        ...         rotor_id=dict(type="turbine", tag="Model"),
+        ...     )
+        ... )
+        >>> configs.rotor_properties # doctest: +ELLIPSIS
+        {
+            "rotor_speeds": {
+                "min_speed": 1000.0,
+                "max_speed": 10000.0,
+                "oper_speed": None,
+                "trip_speed": None,
+                "speed_factor": 1.5,
+                "unit": "rpm",
+            },
+            "rotor_id": {"type": "turbine", "tag": "Model"},
+        }...
+        """
+        for k, v in kwargs.items():
+            if k not in self.__dict__.keys():
+                raise KeyError("Option '{}' not found.".format(k))
+            else:
+                getattr(self, k)._update(**v)
diff --git a/ross/disk_element.py b/ross/disk_element.py
index 235c56ee..52a9d03e 100644
--- a/ross/disk_element.py
+++ b/ross/disk_element.py
@@ -126,6 +126,66 @@ class DiskElement(Element):
     def __hash__(self):
         return hash(self.tag)
 
+    def save(self, file_name=os.getcwd()):
+        """Save a disk 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 disk element will be saved in.
+
+        Examples
+        --------
+        >>> disk = disk_example()
+        >>> disk.save()
+        """
+        data = self.get_data(Path(file_name) / "DiskElement.toml")
+        data["DiskElement"][str(self.n)] = {
+            "n": self.n,
+            "m": self.m,
+            "Id": self.Id,
+            "Ip": self.Ip,
+            "tag": self.tag,
+            "scale_factor": self.scale_factor,
+            "color": self.color,
+        }
+        self.dump_data(data, Path(file_name) / "DiskElement.toml")
+
+    @staticmethod
+    def load(file_name=os.getcwd()):
+        """Load a list of disk elements saved in a toml format.
+
+        It works as an auxiliary function of the load function in the Rotor class.
+
+        Parameters
+        ----------
+        file_name: str
+            The name of the file of the disk element to be loaded.
+
+        Returns
+        -------
+        disk_elements: list
+            A list of disk elements.
+
+        Examples
+        --------
+        >>> disk1 = disk_example()
+        >>> disk1.save(os.getcwd())
+        >>> list_of_disks = DiskElement.load(os.getcwd())
+        >>> disk1 == list_of_disks[0]
+        True
+        """
+        disk_elements = []
+        with open("DiskElement.toml", "r") as f:
+            disk_elements_dict = toml.load(f)
+            for element in disk_elements_dict["DiskElement"]:
+                disk_elements.append(
+                    DiskElement(**disk_elements_dict["DiskElement"][element])
+                )
+        return disk_elements
+
     def dof_mapping(self):
         """Degrees of freedom mapping.
 
@@ -294,7 +354,12 @@ class DiskElement(Element):
         y_upper.append(None)
         y_pos.extend(y_lower)
 
-        customdata = [self.n, self.Ip, self.Id, self.m]
+        customdata = [
+            self.n,
+            self.Ip,
+            self.Id,
+            self.m,
+        ]
         hovertemplate = (
             f"<b>Disk Node: {customdata[0]}<b><br>"
             + f"<b>Polar Inertia: {customdata[1]:.3e}<b><br>"
@@ -354,7 +419,15 @@ class DiskElement(Element):
     @classmethod
     @check_units
     def from_geometry(
-        cls, n, material, width, i_d, o_d, tag=None, scale_factor=1.0, color="Firebrick"
+        cls,
+        n,
+        material,
+        width,
+        i_d,
+        o_d,
+        tag=None,
+        scale_factor=1.0,
+        color="Firebrick",
     ):
         """Create a disk element from geometry properties.
 
@@ -412,7 +485,9 @@ class DiskElement(Element):
         return cls(n, m, Id, Ip, tag, scale_factor, color)
 
     @classmethod
-    def from_table(cls, file, sheet_name=0, tag=None, scale_factor=None, color=None):
+    def from_table(
+        cls, file, sheet_name=0, tag=None, scale_factor=None, color=None,
+    ):
         """Instantiate one or more disks using inputs from an Excel table.
 
         A header with the names of the columns is required. These names should
diff --git a/ross/element.py b/ross/element.py
index ceeac373..60ab85cb 100644
--- a/ross/element.py
+++ b/ross/element.py
@@ -1,4 +1,3 @@
-import inspect
 from abc import ABC, abstractmethod
 from collections import namedtuple
 from pathlib import Path
@@ -18,16 +17,12 @@ class Element(ABC):
         self.n = n
         self.tag = tag
 
-    def save(self, file):
-        """Save the element in a .toml file.
-
-        This function will save the element to a .toml file.
-        The file will have all the argument's names and values that are needed to
-        reinstantiate the element.
+    def save(self, file_name):
+        """Save the element in a file.
 
         Parameters
         ----------
-        file : str, pathlib.Path
+        file_name: string
             The name of the file the element will be saved in.
 
         Examples
@@ -35,32 +30,18 @@ class Element(ABC):
         >>> # Example using DiskElement
         >>> from ross.disk_element import disk_example
         >>> disk = disk_example()
-        >>> disk.save('/tmp/disk.toml')
+        >>> disk.save()
         """
-        # get __init__ arguments
-        signature = inspect.signature(self.__init__)
-        args_list = list(signature.parameters)
-        args = {arg: getattr(self, arg) for arg in args_list}
-        try:
-            data = toml.load(file)
-        except FileNotFoundError:
-            data = {}
-
-        data[f"{self.__class__.__name__}_{self.tag}"] = args
-        with open(file, "w") as f:
-            toml.dump(data, f)
-
-    @classmethod
-    def read_toml_data(cls, data):
-        """Read and parse data stored in a .toml file.
+        pass
 
-        The data passed to this method needs to be according to the
-        format saved in the .toml file by the .save() method.
+    @staticmethod
+    def load(file_name):
+        """Load elements saved in a file.
 
         Parameters
         ----------
-        data : dict
-            Dictionary obtained from toml.load().
+        file_name: string
+            The name of the file to be loaded.
 
         Returns
         -------
@@ -72,19 +53,12 @@ class Element(ABC):
         >>> from ross.bearing_seal_element import bearing_example
         >>> from ross.bearing_seal_element import BearingElement
         >>> bearing1 = bearing_example()
-        >>> bearing1.save('/tmp/bearing1.toml')
-        >>> bearing1_loaded = BearingElement.load('/tmp/bearing1.toml')
-        >>> bearing1 == bearing1_loaded
+        >>> bearing1.save(Path('.'))
+        >>> list_of_bearings = BearingElement.load()
+        >>> bearing1 == list_of_bearings[0]
         True
         """
-        return cls(**data)
-
-    @classmethod
-    def load(cls, file):
-        data = toml.load(file)
-        # extract single dictionary in the data
-        data = list(data.values())[0]
-        return cls.read_toml_data(data)
+        pass
 
     @abstractmethod
     def M(self):
@@ -204,6 +178,65 @@ class Element(ABC):
         attributes["type"] = self.__class__.__name__
         return pd.Series(attributes)
 
+    @staticmethod
+    def get_data(file_name):
+        """Load elements data saved in a toml format.
+
+        Parameters
+        ----------
+        file_name: string
+            The name of the file to be loaded.
+
+        Returns
+        -------
+        data: dict
+            The element parameters presented as a dictionary.
+
+        Examples
+        --------
+        >>> # Example using BearingElement
+        >>> from ross.bearing_seal_element import bearing_example
+        >>> from ross.bearing_seal_element import BearingElement
+        >>> bearing = bearing_example()
+        >>> bearing.save(Path('.'))
+        >>> BearingElement.get_data('BearingElement.toml') # doctest: +ELLIPSIS
+        {'BearingElement': {'0': {'n': 0, 'kxx': [1000000.0, 1000000.0,...
+        """
+        try:
+            with open(file_name, "r") as f:
+                data = toml.load(f)
+                if data == {"": {}}:
+                    data = {str(file_name.name)[:-5]: {}}
+
+        except FileNotFoundError:
+            data = {str(file_name.name)[:-5]: {}}
+            Element.dump_data(data, file_name)
+        return data
+
+    @staticmethod
+    def dump_data(data, file_name):
+        """Dump element data in a toml file.
+
+        Parameters
+        ----------
+        data: dict
+            The data that should be dumped.
+        file_name: string
+            The name of the file the data will be dumped in.
+
+        Examples
+        --------
+        >>> # Example using BearingElement
+        >>> from ross.bearing_seal_element import bearing_example
+        >>> from ross.bearing_seal_element import BearingElement
+        >>> bearing = bearing_example()
+        >>> bearing.save(Path('.'))
+        >>> data = BearingElement.get_data('BearingElement.toml')
+        >>> BearingElement.dump_data(data, 'BearingElement.toml') # doctest: +ELLIPSIS
+        """
+        with open(file_name, "w") as f:
+            toml.dump(data, f)
+
     @abstractmethod
     def dof_mapping(self):
         """Degrees of freedom mapping.
diff --git a/ross/fluid_flow/fluid_flow.py b/ross/fluid_flow/fluid_flow.py
index e3bb5265..7361fc38 100644
--- a/ross/fluid_flow/fluid_flow.py
+++ b/ross/fluid_flow/fluid_flow.py
@@ -48,8 +48,6 @@ class FluidFlow:
         Output Pressure (Pa).
     load: float
         Load applied to the rotor (N).
-    omegap: float
-        Frequency of the rotor (rad/s).
 
     Geometric data of the problem
     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -124,12 +122,6 @@ class FluidFlow:
         y value of the external radius.
     yri : array of shape (nz, ntheta)
         y value of the internal radius.
-    t : float
-        Time.
-    xp : float
-        Perturbation along x.
-    yp : float
-        Perturbation along y.
     eccentricity : float
         distance between the center of the rotor and the stator.
     radial_clearance: float
@@ -199,7 +191,6 @@ class FluidFlow:
         attitude_angle=None,
         eccentricity=None,
         load=None,
-        omegap=None,
         immediately_calculate_pressure_matrix_numerically=True,
     ):
         if load is None and eccentricity is None:
@@ -246,11 +237,6 @@ class FluidFlow:
             self.eccentricity = (
                 calculate_eccentricity_ratio(modified_s) * self.radial_clearance
             )
-        self.omegap = omegap
-        if self.omegap is None:
-            self.omegap = self.omega
-        else:
-            self.omegap = omegap
         self.eccentricity_ratio = self.eccentricity / self.radial_clearance
         if self.load is None:
             self.load = calculate_rotor_load(
@@ -379,7 +365,7 @@ class FluidFlow:
         self.analytical_pressure_matrix_available = True
         return self.p_mat_analytical
 
-    def calculate_coefficients(self, direction=None):
+    def calculate_coefficients(self):
         """This function calculates the constants that form the Poisson equation
         of the discrete pressure (central differences in the second
         derivatives). It is executed when the class is instantiated.
@@ -430,14 +416,6 @@ class FluidFlow:
                 self.c0w[i][j] = (- w * self.ri[i][j] *
                                   (np.log(self.re[i][j] / self.ri[i][j]) *
                                    (1 + (self.ri[i][j] ** 2) / (self.re[i][j] ** 2 - self.ri[i][j] ** 2)) - 1 / 2))
-                if direction == "x":
-                    a = self.omegap * (self.xp) * np.cos(self.omegap * self.t)
-                    self.c0w[i][j] += self.ri[i][j] * a * np.sin(self.gama[i][j])
-                elif direction == "y":
-                    b = self.omegap * (self.yp) * np.cos(self.omegap * self.t)
-                    self.c0w[i][j] -= self.ri[i][j] * b * np.cos(self.gama[i][j])
-                else:
-                    self.c0w[i][j] += 0
                 # fmt: on
                 if not eccentricity_error:
                     if abs(self.xri[i][j]) > abs(self.xre[i][j]) or abs(
diff --git a/ross/fluid_flow/fluid_flow_coefficients.py b/ross/fluid_flow/fluid_flow_coefficients.py
index e3c20087..e34431e7 100644
--- a/ross/fluid_flow/fluid_flow_coefficients.py
+++ b/ross/fluid_flow/fluid_flow_coefficients.py
@@ -145,121 +145,19 @@ def calculate_oil_film_force(fluid_flow_object, force_type=None):
     return radial_force, tangential_force, force_x, force_y
 
 
-def calculate_stiffness_and_damping_coefficients(fluid_flow_object):
-    """This function calculates the bearing stiffness and damping matrices numerically.
-        Parameters
-        ----------
-        fluid_flow_object: A FluidFlow object.
-        Returns
-        -------
-        Two lists of floats
-            A list of length four including stiffness floats in this order: kxx, kxy, kyx, kyy.
-            And another list of length four including damping floats in this order: cxx, cxy, cyx, cyy.
-            And
-        Examples
-        --------
-        >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
-        >>> my_fluid_flow = fluid_flow_example()
-        >>> calculate_stiffness_and_damping_coefficients(my_fluid_flow)  # doctest: +ELLIPSIS
-        ([428...
-        """
-    N = 6
-    t = np.linspace(0, 2 * np.pi / fluid_flow_object.omegap, N)
-    fluid_flow_object.xp = fluid_flow_object.radial_clearance * 0.0001
-    fluid_flow_object.yp = fluid_flow_object.radial_clearance * 0.0001
-    dx = np.zeros(N)
-    dy = np.zeros(N)
-    xdot = np.zeros(N)
-    ydot = np.zeros(N)
-    radial_force = np.zeros(N)
-    tangential_force = np.zeros(N)
-    force_xx = np.zeros(N)
-    force_yx = np.zeros(N)
-    force_xy = np.zeros(N)
-    force_yy = np.zeros(N)
-    X1 = np.zeros([N, 3])
-    X2 = np.zeros([N, 3])
-    F1 = np.zeros(N)
-    F2 = np.zeros(N)
-    F3 = np.zeros(N)
-    F4 = np.zeros(N)
-
-    for i in range(N):
-        fluid_flow_object.t = t[i]
-
-        delta_x = fluid_flow_object.xp * np.sin(
-            fluid_flow_object.omegap * fluid_flow_object.t
-        )
-        move_rotor_center(fluid_flow_object, delta_x, 0)
-        dx[i] = delta_x
-        xdot[i] = (
-            fluid_flow_object.omegap
-            * fluid_flow_object.xp
-            * np.cos(fluid_flow_object.omegap * fluid_flow_object.t)
-        )
-        fluid_flow_object.calculate_coefficients(direction="x")
-        fluid_flow_object.calculate_pressure_matrix_numerical()
-        [
-            radial_force[i],
-            tangential_force[i],
-            force_xx[i],
-            force_yx[i],
-        ] = calculate_oil_film_force(fluid_flow_object, force_type="numerical")
-
-        delta_y = fluid_flow_object.yp * np.sin(
-            fluid_flow_object.omegap * fluid_flow_object.t
-        )
-        move_rotor_center(fluid_flow_object, -delta_x, 0)
-        move_rotor_center(fluid_flow_object, 0, delta_y)
-        dy[i] = delta_y
-        ydot[i] = (
-            fluid_flow_object.omegap
-            * fluid_flow_object.yp
-            * np.cos(fluid_flow_object.omegap * fluid_flow_object.t)
-        )
-        fluid_flow_object.calculate_coefficients(direction="y")
-        fluid_flow_object.calculate_pressure_matrix_numerical()
-        [
-            radial_force[i],
-            tangential_force[i],
-            force_xy[i],
-            force_yy[i],
-        ] = calculate_oil_film_force(fluid_flow_object, force_type="numerical")
-        move_rotor_center(fluid_flow_object, 0, -delta_y)
-        fluid_flow_object.calculate_coefficients()
-        fluid_flow_object.calculate_pressure_matrix_numerical()
-
-        X1[i] = [1, dx[i], xdot[i]]
-        X2[i] = [1, dy[i], ydot[i]]
-        F1[i] = -force_xx[i]
-        F2[i] = -force_xy[i]
-        F3[i] = -force_yx[i]
-        F4[i] = -force_yy[i]
-
-    P1 = np.dot(
-        np.dot(np.linalg.inv(np.dot(np.transpose(X1), X1)), np.transpose(X1)), F1
-    )
-    P2 = np.dot(
-        np.dot(np.linalg.inv(np.dot(np.transpose(X2), X2)), np.transpose(X2)), F2
-    )
-    P3 = np.dot(
-        np.dot(np.linalg.inv(np.dot(np.transpose(X1), X1)), np.transpose(X1)), F3
-    )
-    P4 = np.dot(
-        np.dot(np.linalg.inv(np.dot(np.transpose(X2), X2)), np.transpose(X2)), F4
-    )
-
-    K = [P1[1], P2[1], P3[1], P4[1]]
-    C = [P1[2], P2[2], P3[2], P4[2]]
-
-    return K, C
-
-
-def calculate_short_stiffness_matrix(fluid_flow_object):
-    """This function calculates the stiffness matrix for the short bearing.
+def calculate_stiffness_matrix(
+    fluid_flow_object, force_type=None, oil_film_force="numerical"
+):
+    """This function calculates the bearing stiffness matrix numerically.
     Parameters
     ----------
     fluid_flow_object: A FluidFlow object.
+    oil_film_force: str
+        If set, calculates the oil film force analytically considering the chosen type: 'short' or 'long'.
+        If set to 'numerical', calculates the oil film force numerically.
+    force_type: str
+        If set, calculates the stiffness matrix analytically considering the chosen type: 'short'.
+        If set to 'numerical', calculates the stiffness matrix numerically.
     Returns
     -------
     list of floats
@@ -268,76 +166,117 @@ def calculate_short_stiffness_matrix(fluid_flow_object):
     --------
     >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
     >>> my_fluid_flow = fluid_flow_example()
-    >>> calculate_short_stiffness_matrix(my_fluid_flow)  # doctest: +ELLIPSIS
+    >>> calculate_stiffness_matrix(my_fluid_flow)  # doctest: +ELLIPSIS
     [417...
     """
-    h0 = 1.0 / (
-        (
-            (np.pi ** 2) * (1 - fluid_flow_object.eccentricity_ratio ** 2)
-            + 16 * fluid_flow_object.eccentricity_ratio ** 2
-        )
-        ** 1.5
-    )
-    a = fluid_flow_object.load / fluid_flow_object.radial_clearance
-    kxx = (
-        a
-        * h0
-        * 4
-        * (
-            (np.pi ** 2) * (2 - fluid_flow_object.eccentricity_ratio ** 2)
-            + 16 * fluid_flow_object.eccentricity_ratio ** 2
-        )
-    )
-    kxy = (
-        a
-        * h0
-        * np.pi
-        * (
-            (np.pi ** 2) * (1 - fluid_flow_object.eccentricity_ratio ** 2) ** 2
-            - 16 * fluid_flow_object.eccentricity_ratio ** 4
-        )
-        / (
-            fluid_flow_object.eccentricity_ratio
-            * np.sqrt(1 - fluid_flow_object.eccentricity_ratio ** 2)
+    if force_type != "numerical" and (
+        force_type == "short" or fluid_flow_object.bearing_type == "short_bearing"
+    ):
+        h0 = 1.0 / (
+            (
+                (np.pi ** 2) * (1 - fluid_flow_object.eccentricity_ratio ** 2)
+                + 16 * fluid_flow_object.eccentricity_ratio ** 2
+            )
+            ** 1.5
         )
-    )
-    kyx = (
-        -a
-        * h0
-        * np.pi
-        * (
-            (np.pi ** 2)
-            * (1 - fluid_flow_object.eccentricity_ratio ** 2)
-            * (1 + 2 * fluid_flow_object.eccentricity_ratio ** 2)
-            + (32 * fluid_flow_object.eccentricity_ratio ** 2)
-            * (1 + fluid_flow_object.eccentricity_ratio ** 2)
+        a = fluid_flow_object.load / fluid_flow_object.radial_clearance
+        kxx = (
+            a
+            * h0
+            * 4
+            * (
+                (np.pi ** 2) * (2 - fluid_flow_object.eccentricity_ratio ** 2)
+                + 16 * fluid_flow_object.eccentricity_ratio ** 2
+            )
         )
-        / (
-            fluid_flow_object.eccentricity_ratio
-            * np.sqrt(1 - fluid_flow_object.eccentricity_ratio ** 2)
+        kxy = (
+            a
+            * h0
+            * np.pi
+            * (
+                (np.pi ** 2) * (1 - fluid_flow_object.eccentricity_ratio ** 2) ** 2
+                - 16 * fluid_flow_object.eccentricity_ratio ** 4
+            )
+            / (
+                fluid_flow_object.eccentricity_ratio
+                * np.sqrt(1 - fluid_flow_object.eccentricity_ratio ** 2)
+            )
         )
-    )
-    kyy = (
-        a
-        * h0
-        * 4
-        * (
-            (np.pi ** 2) * (1 + 2 * fluid_flow_object.eccentricity_ratio ** 2)
-            + (
-                (32 * fluid_flow_object.eccentricity_ratio ** 2)
+        kyx = (
+            -a
+            * h0
+            * np.pi
+            * (
+                (np.pi ** 2)
+                * (1 - fluid_flow_object.eccentricity_ratio ** 2)
+                * (1 + 2 * fluid_flow_object.eccentricity_ratio ** 2)
+                + (32 * fluid_flow_object.eccentricity_ratio ** 2)
                 * (1 + fluid_flow_object.eccentricity_ratio ** 2)
             )
-            / (1 - fluid_flow_object.eccentricity_ratio ** 2)
+            / (
+                fluid_flow_object.eccentricity_ratio
+                * np.sqrt(1 - fluid_flow_object.eccentricity_ratio ** 2)
+            )
         )
-    )
+        kyy = (
+            a
+            * h0
+            * 4
+            * (
+                (np.pi ** 2) * (1 + 2 * fluid_flow_object.eccentricity_ratio ** 2)
+                + (
+                    (32 * fluid_flow_object.eccentricity_ratio ** 2)
+                    * (1 + fluid_flow_object.eccentricity_ratio ** 2)
+                )
+                / (1 - fluid_flow_object.eccentricity_ratio ** 2)
+            )
+        )
+    else:
+
+        [radial_force, tangential_force, force_x, force_y] = calculate_oil_film_force(
+            fluid_flow_object, force_type=oil_film_force
+        )
+        delta = fluid_flow_object.radial_clearance / 100
+
+        move_rotor_center(fluid_flow_object, delta, 0)
+        fluid_flow_object.calculate_coefficients()
+        fluid_flow_object.calculate_pressure_matrix_numerical()
+        [
+            radial_force_x,
+            tangential_force_x,
+            force_x_x,
+            force_y_x,
+        ] = calculate_oil_film_force(fluid_flow_object, force_type=oil_film_force)
+
+        move_rotor_center(fluid_flow_object, -delta, 0)
+        move_rotor_center(fluid_flow_object, 0, delta)
+        fluid_flow_object.calculate_coefficients()
+        fluid_flow_object.calculate_pressure_matrix_numerical()
+        [
+            radial_force_y,
+            tangential_force_y,
+            force_x_y,
+            force_y_y,
+        ] = calculate_oil_film_force(fluid_flow_object, force_type=oil_film_force)
+        move_rotor_center(fluid_flow_object, 0, -delta)
+
+        kxx = (force_x - force_x_x) / delta
+        kyx = (force_y - force_y_x) / delta
+        kxy = (force_x - force_x_y) / delta
+        kyy = (force_y - force_y_y) / delta
+
     return [kxx, kxy, kyx, kyy]
 
 
-def calculate_short_damping_matrix(fluid_flow_object):
-    """This function calculates the damping matrix for the short bearing.
+def calculate_damping_matrix(fluid_flow_object, force_type=None):
+    """Returns the damping matrix calculated analytically.
+    Suitable only for short bearings.
     Parameters
     -------
     fluid_flow_object: A FluidFlow object.
+    force_type: str
+        If set, calculates the stiffness matrix analytically considering the chosen type: 'short'.
+        If set to 'numerical', calculates the stiffness matrix numerically.
     Returns
     -------
     list of floats
@@ -346,23 +285,29 @@ def calculate_short_damping_matrix(fluid_flow_object):
     --------
     >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
     >>> my_fluid_flow = fluid_flow_example()
-    >>> calculate_short_damping_matrix(my_fluid_flow) # doctest: +ELLIPSIS
+    >>> calculate_damping_matrix(my_fluid_flow) # doctest: +ELLIPSIS
     [...
     """
     # fmt: off
-    h0 = 1.0 / (((np.pi ** 2) * (1 - fluid_flow_object.eccentricity_ratio ** 2)
-                 + 16 * fluid_flow_object.eccentricity_ratio ** 2) ** 1.5)
-    a = fluid_flow_object.load / (fluid_flow_object.radial_clearance * fluid_flow_object.omega)
-    cxx = (a * h0 * 2 * np.pi * np.sqrt(1 - fluid_flow_object.eccentricity_ratio ** 2) *
-           ((np.pi ** 2) * (1 + 2 * fluid_flow_object.eccentricity_ratio ** 2)
-            - 16 * fluid_flow_object.eccentricity_ratio ** 2) / fluid_flow_object.eccentricity_ratio)
-    cxy = (-a * h0 * 8 * ((np.pi ** 2) * (1 + 2 * fluid_flow_object.eccentricity_ratio ** 2)
-                          - 16 * fluid_flow_object.eccentricity_ratio ** 2))
-    cyx = cxy
-    cyy = (a * h0 * (2 * np.pi * (
-            (np.pi ** 2) * (1 - fluid_flow_object.eccentricity_ratio ** 2) ** 2
-            + 48 * fluid_flow_object.eccentricity_ratio ** 2)) /
-           (fluid_flow_object.eccentricity_ratio * np.sqrt(1 - fluid_flow_object.eccentricity_ratio ** 2)))
+    if force_type != 'numerical' and (force_type == 'short' or fluid_flow_object.bearing_type == 'short_bearing'):
+        h0 = 1.0 / (((np.pi ** 2) * (1 - fluid_flow_object.eccentricity_ratio ** 2)
+                     + 16 * fluid_flow_object.eccentricity_ratio ** 2) ** 1.5)
+        a = fluid_flow_object.load / (fluid_flow_object.radial_clearance * fluid_flow_object.omega)
+        cxx = (a * h0 * 2 * np.pi * np.sqrt(1 - fluid_flow_object.eccentricity_ratio ** 2) *
+               ((np.pi ** 2) * (1 + 2 * fluid_flow_object.eccentricity_ratio ** 2)
+                - 16 * fluid_flow_object.eccentricity_ratio ** 2) / fluid_flow_object.eccentricity_ratio)
+        cxy = (-a * h0 * 8 * ((np.pi ** 2) * (1 + 2 * fluid_flow_object.eccentricity_ratio ** 2)
+                              - 16 * fluid_flow_object.eccentricity_ratio ** 2))
+        cyx = cxy
+        cyy = (a * h0 * (2 * np.pi * (
+                (np.pi ** 2) * (1 - fluid_flow_object.eccentricity_ratio ** 2) ** 2
+                + 48 * fluid_flow_object.eccentricity_ratio ** 2)) /
+               (fluid_flow_object.eccentricity_ratio * np.sqrt(1 - fluid_flow_object.eccentricity_ratio ** 2)))
+    else:
+        cxx, cxy, cyx, cyy = warnings.warn(
+            "Function calculate_damping_matrix cannot yet be calculated numerically. "
+            "The force_type should be  'short' or 'short_bearing"
+        )
     # fmt: on
     return [cxx, cxy, cyx, cyy]
 
diff --git a/ross/fluid_flow/fluid_flow_graphics.py b/ross/fluid_flow/fluid_flow_graphics.py
index 69d69be7..21969673 100644
--- a/ross/fluid_flow/fluid_flow_graphics.py
+++ b/ross/fluid_flow/fluid_flow_graphics.py
@@ -1,10 +1,11 @@
 import numpy as np
 import plotly.graph_objects as go
+import plotly.io as pio
 
-from ross.plotly_theme import tableau_colors
+pio.renderers.default = "browser"
 
 
-def plot_eccentricity(fluid_flow_object, z=0, fig=None, **kwargs):
+def plot_eccentricity(fluid_flow_object, z=0, **kwargs):
     """Plot the rotor eccentricity.
 
     This function assembles pressure graphic along the z-axis.
@@ -15,8 +16,6 @@ def plot_eccentricity(fluid_flow_object, z=0, fig=None, **kwargs):
     fluid_flow_object: a FluidFlow object
     z: int, optional
         The distance in z where to cut and plot.
-    fig : Plotly graph_objects.Figure()
-        The figure object with the plot.
     kwargs : optional
         Additional key word arguments can be passed to change the plot layout only
         (e.g. width=1000, height=800, ...).
@@ -35,15 +34,29 @@ def plot_eccentricity(fluid_flow_object, z=0, fig=None, **kwargs):
     >>> # to show the plots you can use:
     >>> # show(fig)
     """
-    if fig is None:
-        fig = go.Figure()
+    kwargs_default_values = dict(
+        width=600,
+        height=600,
+        plot_bgcolor="white",
+        hoverlabel_align="right",
+        legend=dict(
+            font=dict(family="sans-serif", size=14),
+            bgcolor="white",
+            bordercolor="black",
+            borderwidth=2,
+        ),
+    )
+    for k, v in kwargs_default_values.items():
+        kwargs.setdefault(k, v)
+
+    fig = go.Figure()
     fig.add_trace(
         go.Scatter(
             x=fluid_flow_object.xre[z],
             y=fluid_flow_object.yre[z],
             mode="markers+lines",
-            marker=dict(color=tableau_colors["red"]),
-            line=dict(color=tableau_colors["red"]),
+            marker=dict(size=10, color="firebrick"),
+            line=dict(width=2.0, color="firebrick"),
             name="Stator",
             legendgroup="Stator",
             hovertemplate=("<b>X: %{x:.3e}</b><br>" + "<b>Y: %{y:.3e}</b>"),
@@ -54,8 +67,8 @@ def plot_eccentricity(fluid_flow_object, z=0, fig=None, **kwargs):
             x=fluid_flow_object.xri[z],
             y=fluid_flow_object.yri[z],
             mode="markers+lines",
-            marker=dict(color=tableau_colors["blue"]),
-            line=dict(color=tableau_colors["blue"]),
+            marker=dict(size=10, color="royalblue"),
+            line=dict(width=2.0, color="royalblue"),
             name="Rotor",
             legendgroup="Rotor",
             hovertemplate=("<b>X: %{x:.3e}</b><br>" + "<b>Y: %{y:.3e}</b>"),
@@ -65,7 +78,7 @@ def plot_eccentricity(fluid_flow_object, z=0, fig=None, **kwargs):
         go.Scatter(
             x=[fluid_flow_object.xi],
             y=[fluid_flow_object.yi],
-            marker=dict(color=tableau_colors["red"]),
+            marker=dict(size=10, color="firebrick"),
             name="Stator",
             legendgroup="Stator",
             showlegend=False,
@@ -76,24 +89,42 @@ def plot_eccentricity(fluid_flow_object, z=0, fig=None, **kwargs):
         go.Scatter(
             x=[0],
             y=[0],
-            marker=dict(color=tableau_colors["blue"]),
+            marker=dict(size=10, color="royalblue"),
             name="Rotor",
             legendgroup="Rotor",
             showlegend=False,
             hovertemplate=("<b>X: %{x:.3e}</b><br>" + "<b>Y: %{y:.3e}</b>"),
         )
     )
-    fig.update_xaxes(title_text="<b>X axis</b>")
-    fig.update_yaxes(title_text="<b>Y axis</b>")
+    fig.update_xaxes(
+        title_text="<b>X axis</b>",
+        title_font=dict(size=16),
+        tickfont=dict(size=14),
+        gridcolor="lightgray",
+        showline=True,
+        linewidth=2.5,
+        linecolor="black",
+        mirror=True,
+    )
+    fig.update_yaxes(
+        title_text="<b>Y axis</b>",
+        title_font=dict(size=16),
+        tickfont=dict(size=14),
+        gridcolor="lightgray",
+        showline=True,
+        linewidth=2.5,
+        linecolor="black",
+        mirror=True,
+    )
     fig.update_layout(
-        title=dict(text="<b>Cut in plane Z={}</b>".format(z)),
-        yaxis=dict(scaleanchor="x", scaleratio=1, **kwargs),
+        title=dict(text="<b>Cut in plane Z={}</b>".format(z), font=dict(size=16)),
+        **kwargs,
     )
 
     return fig
 
 
-def plot_pressure_z(fluid_flow_object, theta=0, fig=None, **kwargs):
+def plot_pressure_z(fluid_flow_object, theta=0, **kwargs):
     """Plot the pressure distribution along the z-axis.
 
     This function assembles pressure graphic along the z-axis for one or both the
@@ -105,8 +136,6 @@ def plot_pressure_z(fluid_flow_object, theta=0, fig=None, **kwargs):
     fluid_flow_object: a FluidFlow object
     theta: int, optional
         The theta to be considered.
-    fig : Plotly graph_objects.Figure()
-        The figure object with the plot.
     kwargs : optional
         Additional key word arguments can be passed to change the plot layout only
         (e.g. width=1000, height=800, ...).
@@ -135,16 +164,29 @@ def plot_pressure_z(fluid_flow_object, theta=0, fig=None, **kwargs):
             "Must calculate the pressure matrix. "
             "Try calling calculate_pressure_matrix_numerical() or calculate_pressure_matrix_analytical() first."
         )
+    kwargs_default_values = dict(
+        width=800,
+        height=600,
+        plot_bgcolor="white",
+        hoverlabel_align="right",
+        legend=dict(
+            font=dict(family="sans-serif", size=14),
+            bgcolor="white",
+            bordercolor="black",
+            borderwidth=2,
+        ),
+    )
+    for k, v in kwargs_default_values.items():
+        kwargs.setdefault(k, v)
 
-    if fig is None:
-        fig = go.Figure()
+    fig = go.Figure()
     if fluid_flow_object.numerical_pressure_matrix_available:
         fig.add_trace(
             go.Scatter(
                 x=fluid_flow_object.z_list,
                 y=fluid_flow_object.p_mat_numerical[:, theta],
                 mode="lines",
-                line=dict(color=tableau_colors["blue"]),
+                line=dict(width=3.0, color="royalblue"),
                 showlegend=True,
                 name="<b>Numerical pressure</b>",
                 hovertemplate=(
@@ -159,7 +201,7 @@ def plot_pressure_z(fluid_flow_object, theta=0, fig=None, **kwargs):
                 x=fluid_flow_object.z_list,
                 y=fluid_flow_object.p_mat_analytical[:, theta],
                 mode="lines",
-                line=dict(color=tableau_colors["red"]),
+                line=dict(width=3.0, color="firebrick"),
                 showlegend=True,
                 name="<b>Analytical pressure</b>",
                 hovertemplate=(
@@ -168,14 +210,33 @@ def plot_pressure_z(fluid_flow_object, theta=0, fig=None, **kwargs):
                 ),
             )
         )
-    fig.update_xaxes(title_text="<b>Axial Length</b>")
-    fig.update_yaxes(title_text="<b>Pressure</b>")
+    fig.update_xaxes(
+        title_text="<b>Axial Length</b>",
+        title_font=dict(size=16),
+        tickfont=dict(size=14),
+        gridcolor="lightgray",
+        showline=True,
+        linewidth=2.5,
+        linecolor="black",
+        mirror=True,
+    )
+    fig.update_yaxes(
+        title_text="<b>Pressure</b>",
+        title_font=dict(size=16),
+        tickfont=dict(size=14),
+        gridcolor="lightgray",
+        showline=True,
+        linewidth=2.5,
+        linecolor="black",
+        mirror=True,
+    )
     fig.update_layout(
         title=dict(
             text=(
                 "<b>Pressure along the flow (axial direction)<b><br>"
                 + "<b>Theta={}</b>".format(theta)
-            )
+            ),
+            font=dict(size=16),
         ),
         **kwargs,
     )
@@ -183,7 +244,7 @@ def plot_pressure_z(fluid_flow_object, theta=0, fig=None, **kwargs):
     return fig
 
 
-def plot_shape(fluid_flow_object, theta=0, fig=None, **kwargs):
+def plot_shape(fluid_flow_object, theta=0, **kwargs):
     """Plot the surface geometry of the rotor.
 
     This function assembles a graphic representing the geometry of the rotor.
@@ -193,8 +254,6 @@ def plot_shape(fluid_flow_object, theta=0, fig=None, **kwargs):
     fluid_flow_object: a FluidFlow object
     theta: int, optional
         The theta to be considered.
-    fig : Plotly graph_objects.Figure()
-        The figure object with the plot.
 
     Returns
     -------
@@ -213,14 +272,28 @@ def plot_shape(fluid_flow_object, theta=0, fig=None, **kwargs):
     >>> # to show the plots you can use:
     >>> # fig.show()
     """
-    if fig is None:
-        fig = go.Figure()
+    kwargs_default_values = dict(
+        width=800,
+        height=600,
+        plot_bgcolor="white",
+        hoverlabel_align="right",
+        legend=dict(
+            font=dict(family="sans-serif", size=14),
+            bgcolor="white",
+            bordercolor="black",
+            borderwidth=2,
+        ),
+    )
+    for k, v in kwargs_default_values.items():
+        kwargs.setdefault(k, v)
+
+    fig = go.Figure()
     fig.add_trace(
         go.Scatter(
             x=fluid_flow_object.z_list,
             y=fluid_flow_object.re[:, theta],
             mode="lines",
-            line=dict(color=tableau_colors["red"]),
+            line=dict(width=3.0, color="firebrick"),
             showlegend=True,
             hoverinfo="none",
             name="<b>Stator</b>",
@@ -231,20 +304,39 @@ def plot_shape(fluid_flow_object, theta=0, fig=None, **kwargs):
             x=fluid_flow_object.z_list,
             y=fluid_flow_object.ri[:, theta],
             mode="lines",
-            line=dict(color=tableau_colors["blue"]),
+            line=dict(width=3.0, color="royalblue"),
             showlegend=True,
             hoverinfo="none",
             name="<b>Rotor</b>",
         )
     )
-    fig.update_xaxes(title_text="<b>Axial Length</b>")
-    fig.update_yaxes(title_text="<b>Radial direction</b>")
+    fig.update_xaxes(
+        title_text="<b>Axial Length</b>",
+        title_font=dict(size=16),
+        tickfont=dict(size=14),
+        gridcolor="lightgray",
+        showline=True,
+        linewidth=2.5,
+        linecolor="black",
+        mirror=True,
+    )
+    fig.update_yaxes(
+        title_text="<b>Radial direction</b>",
+        title_font=dict(size=16),
+        tickfont=dict(size=14),
+        gridcolor="lightgray",
+        showline=True,
+        linewidth=2.5,
+        linecolor="black",
+        mirror=True,
+    )
     fig.update_layout(
         title=dict(
             text=(
                 "<b>Shapes of stator and rotor - Axial direction<b><br>"
                 + "<b>Theta={}</b>".format(theta)
-            )
+            ),
+            font=dict(size=16),
         ),
         **kwargs,
     )
@@ -252,7 +344,7 @@ def plot_shape(fluid_flow_object, theta=0, fig=None, **kwargs):
     return fig
 
 
-def plot_pressure_theta(fluid_flow_object, z=0, fig=None, **kwargs):
+def plot_pressure_theta(fluid_flow_object, z=0, **kwargs):
     """Plot the pressure distribution along theta.
 
     This function assembles pressure graphic in the theta direction for a given z
@@ -264,8 +356,6 @@ def plot_pressure_theta(fluid_flow_object, z=0, fig=None, **kwargs):
     fluid_flow_object: a FluidFlow object
     z: int, optional
         The distance along z-axis to be considered.
-    fig : Plotly graph_objects.Figure()
-        The figure object with the plot.
     kwargs : optional
         Additional key word arguments can be passed to change the plot layout only
         (e.g. width=1000, height=800, ...).
@@ -294,15 +384,29 @@ def plot_pressure_theta(fluid_flow_object, z=0, fig=None, **kwargs):
             "Must calculate the pressure matrix. "
             "Try calling calculate_pressure_matrix_numerical() or calculate_pressure_matrix_analytical() first."
         )
-    if fig is None:
-        fig = go.Figure()
+    kwargs_default_values = dict(
+        width=800,
+        height=600,
+        plot_bgcolor="white",
+        hoverlabel_align="right",
+        legend=dict(
+            font=dict(family="sans-serif", size=14),
+            bgcolor="white",
+            bordercolor="black",
+            borderwidth=2,
+        ),
+    )
+    for k, v in kwargs_default_values.items():
+        kwargs.setdefault(k, v)
+
+    fig = go.Figure()
     if fluid_flow_object.numerical_pressure_matrix_available:
         fig.add_trace(
             go.Scatter(
                 x=fluid_flow_object.gama[z],
                 y=fluid_flow_object.p_mat_numerical[z],
                 mode="lines",
-                line=dict(color=tableau_colors["blue"]),
+                line=dict(width=3.0, color="royalblue"),
                 showlegend=True,
                 name="<b>Numerical pressure</b>",
                 hovertemplate=(
@@ -316,7 +420,7 @@ def plot_pressure_theta(fluid_flow_object, z=0, fig=None, **kwargs):
                 x=fluid_flow_object.gama[z],
                 y=fluid_flow_object.p_mat_analytical[z],
                 mode="lines",
-                line=dict(color=tableau_colors["red"]),
+                line=dict(width=3.0, color="firebrick"),
                 showlegend=True,
                 name="<b>Analytical pressure</b>",
                 hovertemplate=(
@@ -326,17 +430,38 @@ def plot_pressure_theta(fluid_flow_object, z=0, fig=None, **kwargs):
             )
         )
 
-    fig.update_xaxes(title_text="<b>Theta value</b>")
-    fig.update_yaxes(title_text="<b>Pressure</b>")
+    fig.update_xaxes(
+        title_text="<b>Theta value</b>",
+        title_font=dict(size=16),
+        tickfont=dict(size=14),
+        gridcolor="lightgray",
+        showline=True,
+        linewidth=2.5,
+        linecolor="black",
+        mirror=True,
+    )
+    fig.update_yaxes(
+        title_text="<b>Pressure</b>",
+        title_font=dict(size=16),
+        tickfont=dict(size=14),
+        gridcolor="lightgray",
+        showline=True,
+        linewidth=2.5,
+        linecolor="black",
+        mirror=True,
+    )
     fig.update_layout(
-        title=dict(text=("<b>Pressure along Theta | Z={}<b>".format(z))), **kwargs
+        title=dict(
+            text=("<b>Pressure along Theta | Z={}<b>".format(z)), font=dict(size=16)
+        ),
+        **kwargs,
     )
 
     return fig
 
 
 def plot_pressure_theta_cylindrical(
-    fluid_flow_object, z=0, from_numerical=True, fig=None, **kwargs
+    fluid_flow_object, z=0, from_numerical=True, **kwargs
 ):
     """Plot cylindrical pressure graphic in the theta direction.
 
@@ -353,8 +478,6 @@ def plot_pressure_theta_cylindrical(
         If False, takes the analytically calculated one instead.
         If condition cannot be satisfied (matrix not calculated), it will take the one
         that is available and raise a warning.
-    fig : Plotly graph_objects.Figure()
-        The figure object with the plot.
     kwargs : optional
         Additional key word arguments can be passed to change the plot layout only
         (e.g. width=1000, height=800, ...).
@@ -421,9 +544,11 @@ def plot_pressure_theta_cylindrical(
                 continue
             z_matrix[i][j] = pressure_along_theta[i] - min_pressure + 0.01
 
-    if fig is None:
-        fig = go.Figure()
-    fig.add_trace(
+    kwargs_default_values = dict(width=1200, height=900)
+    for k, v in kwargs_default_values.items():
+        kwargs.setdefault(k, v)
+
+    fig = go.Figure(
         go.Barpolar(
             r=r_matrix.ravel(),
             theta=theta_matrix.ravel(),
@@ -433,7 +558,10 @@ def plot_pressure_theta_cylindrical(
                 colorscale="Viridis",
                 cmin=np.amin(z_matrix),
                 cmax=np.amax(z_matrix),
-                colorbar=dict(title=dict(text="<b>Pressure</b>", side="top")),
+                colorbar=dict(
+                    title=dict(text="<b>Pressure</b>", side="top", font=dict(size=16)),
+                    tickfont=dict(size=16),
+                ),
             ),
             thetaunit="degrees",
             name="Pressure",
@@ -448,9 +576,14 @@ def plot_pressure_theta_cylindrical(
     fig.update_layout(
         polar=dict(
             hole=0.5,
+            bgcolor="white",
             bargap=0.0,
+            radialaxis=dict(gridcolor="lightgray", nticks=5),
             angularaxis=dict(
-                rotation=-90 - fluid_flow_object.attitude_angle * 180 / np.pi
+                rotation=-90 - fluid_flow_object.attitude_angle * 180 / np.pi,
+                gridcolor="lightgray",
+                linecolor="black",
+                linewidth=2.5,
             ),
         ),
         **kwargs,
@@ -458,15 +591,12 @@ def plot_pressure_theta_cylindrical(
     return fig
 
 
-def plot_pressure_surface(fluid_flow_object, fig=None, **kwargs):
+def plot_pressure_surface(fluid_flow_object, **kwargs):
     """Assembles pressure surface graphic in the bearing, using Plotly.
 
     Parameters
     ----------
     fluid_flow_object: a FluidFlow object
-    fig : Plotly graph_objects.Figure()
-        The figure object with the plot.
-
     kwargs : optional
         Additional key word arguments can be passed to change the plot layout only
         (e.g. width=1000, height=800, ...).
@@ -495,8 +625,12 @@ def plot_pressure_surface(fluid_flow_object, fig=None, **kwargs):
             "Must calculate the pressure matrix. "
             "Try calling calculate_pressure_matrix_numerical() or calculate_pressure_matrix_analytical() first."
         )
-    if fig is None:
-        fig = go.Figure()
+
+    kwargs_default_values = dict(width=1200, height=900)
+    for k, v in kwargs_default_values.items():
+        kwargs.setdefault(k, v)
+
+    fig = go.Figure()
     if fluid_flow_object.numerical_pressure_matrix_available:
         z, theta = np.meshgrid(fluid_flow_object.z_list, fluid_flow_object.gama[0])
         fig.add_trace(
@@ -507,7 +641,33 @@ def plot_pressure_surface(fluid_flow_object, fig=None, **kwargs):
                 colorscale="Viridis",
                 cmin=np.amin(fluid_flow_object.p_mat_numerical.T),
                 cmax=np.amax(fluid_flow_object.p_mat_numerical.T),
-                colorbar=dict(title=dict(text="<b>Pressure</b>", side="top")),
+                colorbar=dict(
+                    title=dict(text="<b>Pressure</b>", side="top", font=dict(size=16)),
+                    tickfont=dict(size=16),
+                ),
+                name="Pressure",
+                showlegend=False,
+                hovertemplate=(
+                    "<b>Length: %{x:.2e}</b><br>"
+                    + "<b>Angular Position: %{y:.2f}</b><br>"
+                    + "<b>Pressure: %{z:.2f}</b>"
+                ),
+            )
+        )
+    if fluid_flow_object.analytical_pressure_matrix_available:
+        z, theta = np.meshgrid(fluid_flow_object.z_list, fluid_flow_object.gama[0])
+        fig.add_trace(
+            go.Surface(
+                x=z,
+                y=theta,
+                z=fluid_flow_object.p_mat_analytical.T,
+                colorscale="Viridis",
+                cmin=np.amin(fluid_flow_object.p_mat_analytical.T),
+                cmax=np.amax(fluid_flow_object.p_mat_analytical.T),
+                colorbar=dict(
+                    title=dict(text="<b>Pressure</b>", side="top", font=dict(size=16)),
+                    tickfont=dict(size=16),
+                ),
                 name="Pressure",
                 showlegend=False,
                 hovertemplate=(
@@ -521,11 +681,32 @@ def plot_pressure_surface(fluid_flow_object, fig=None, **kwargs):
     fig.update_layout(
         scene=dict(
             bgcolor="white",
-            xaxis=dict(title=dict(text="<b>Rotor Length</b>")),
-            yaxis=dict(title=dict(text="<b>Angular Position</b>")),
-            zaxis=dict(title=dict(text="<b>Pressure</b>")),
+            xaxis=dict(
+                title=dict(text="<b>Rotor Length</b>", font=dict(size=14)),
+                tickfont=dict(size=16),
+                nticks=5,
+                backgroundcolor="lightgray",
+                gridcolor="white",
+                showspikes=False,
+            ),
+            yaxis=dict(
+                title=dict(text="<b>Angular Position</b>", font=dict(size=14),),
+                tickfont=dict(size=16),
+                nticks=5,
+                backgroundcolor="lightgray",
+                gridcolor="white",
+                showspikes=False,
+            ),
+            zaxis=dict(
+                title=dict(text="<b>Pressure</b>", font=dict(size=14),),
+                tickfont=dict(size=16),
+                nticks=5,
+                backgroundcolor="lightgray",
+                gridcolor="white",
+                showspikes=False,
+            ),
         ),
-        title=dict(text="<b>Bearing Pressure Field</b>"),
+        title=dict(text=("<b>Bearing Pressure Field</b>"), font=dict(size=20),),
         **kwargs,
     )
 
diff --git a/ross/materials.py b/ross/materials.py
index 91233fc5..986d057e 100644
--- a/ross/materials.py
+++ b/ross/materials.py
@@ -210,6 +210,7 @@ class Material:
         Examples
         --------
         >>> import ross as rs
+        >>> rs.steel.save_material()
         >>> steel = rs.Material.use_material('Steel')
         """
         data = Material.get_data()
diff --git a/ross/plotly_theme.py b/ross/plotly_theme.py
deleted file mode 100644
index b038bf47..00000000
--- a/ross/plotly_theme.py
+++ /dev/null
@@ -1,377 +0,0 @@
-import plotly.graph_objects as go
-import plotly.io as pio
-
-# tableau colors
-tableau_colors = {
-    "blue": "#1f77b4",
-    "orange": "#ff7f0e",
-    "green": "#2ca02c",
-    "red": "#d62728",
-    "purple": "#9467bd",
-    "brown": "#8c564b",
-    "pink": "#e377c2",
-    "gray": "#7f7f7f",
-    "olive": "#bcbd22",
-    "cyan": "#17becf",
-}
-pio.templates["ross"] = go.layout.Template(
-    layout={
-        "annotationdefaults": {
-            "arrowcolor": "#2a3f5f",
-            "arrowhead": 0,
-            "arrowwidth": 1,
-        },
-        "coloraxis": {"colorbar": {"outlinewidth": 0, "ticks": ""}},
-        "colorscale": {
-            "diverging": [
-                [0, "#8e0152"],
-                [0.1, "#c51b7d"],
-                [0.2, "#de77ae"],
-                [0.3, "#f1b6da"],
-                [0.4, "#fde0ef"],
-                [0.5, "#f7f7f7"],
-                [0.6, "#e6f5d0"],
-                [0.7, "#b8e186"],
-                [0.8, "#7fbc41"],
-                [0.9, "#4d9221"],
-                [1, "#276419"],
-            ],
-            "sequential": [
-                [0.0, "#0d0887"],
-                [0.1111111111111111, "#46039f"],
-                [0.2222222222222222, "#7201a8"],
-                [0.3333333333333333, "#9c179e"],
-                [0.4444444444444444, "#bd3786"],
-                [0.5555555555555556, "#d8576b"],
-                [0.6666666666666666, "#ed7953"],
-                [0.7777777777777778, "#fb9f3a"],
-                [0.8888888888888888, "#fdca26"],
-                [1.0, "#f0f921"],
-            ],
-            "sequentialminus": [
-                [0.0, "#0d0887"],
-                [0.1111111111111111, "#46039f"],
-                [0.2222222222222222, "#7201a8"],
-                [0.3333333333333333, "#9c179e"],
-                [0.4444444444444444, "#bd3786"],
-                [0.5555555555555556, "#d8576b"],
-                [0.6666666666666666, "#ed7953"],
-                [0.7777777777777778, "#fb9f3a"],
-                [0.8888888888888888, "#fdca26"],
-                [1.0, "#f0f921"],
-            ],
-        },
-        "colorway": list(tableau_colors.values()),
-        "font": {"color": "#2a3f5f"},
-        "geo": {
-            "bgcolor": "white",
-            "lakecolor": "white",
-            "landcolor": "white",
-            "showlakes": True,
-            "showland": True,
-            "subunitcolor": "#C8D4E3",
-        },
-        "hoverlabel": {"align": "left"},
-        "hovermode": "closest",
-        "mapbox": {"style": "light"},
-        "paper_bgcolor": "white",
-        "plot_bgcolor": "white",
-        "polar": {
-            "angularaxis": {
-                "gridcolor": "#EBF0F8",
-                "linecolor": "#EBF0F8",
-                "ticks": "",
-            },
-            "bgcolor": "white",
-            "radialaxis": {"gridcolor": "#EBF0F8", "linecolor": "#EBF0F8", "ticks": ""},
-        },
-        "scene": {
-            "xaxis": {
-                "backgroundcolor": "white",
-                "gridcolor": "#DFE8F3",
-                "gridwidth": 2,
-                "linecolor": "#EBF0F8",
-                "showbackground": True,
-                "ticks": "",
-                "zerolinecolor": "#EBF0F8",
-                "showspikes": False,
-            },
-            "yaxis": {
-                "backgroundcolor": "white",
-                "gridcolor": "#DFE8F3",
-                "gridwidth": 2,
-                "linecolor": "#EBF0F8",
-                "showbackground": True,
-                "ticks": "",
-                "zerolinecolor": "#EBF0F8",
-                "showspikes": False,
-            },
-            "zaxis": {
-                "backgroundcolor": "white",
-                "gridcolor": "#DFE8F3",
-                "gridwidth": 2,
-                "linecolor": "#EBF0F8",
-                "showbackground": True,
-                "ticks": "",
-                "zerolinecolor": "#EBF0F8",
-                "showspikes": False,
-            },
-        },
-        "shapedefaults": {"line": {"color": "#2a3f5f"}},
-        "ternary": {
-            "aaxis": {"gridcolor": "#DFE8F3", "linecolor": "#A2B1C6", "ticks": ""},
-            "baxis": {"gridcolor": "#DFE8F3", "linecolor": "#A2B1C6", "ticks": ""},
-            "bgcolor": "white",
-            "caxis": {"gridcolor": "#DFE8F3", "linecolor": "#A2B1C6", "ticks": ""},
-        },
-        "title": {"x": 0.05},
-        "xaxis": {
-            "automargin": True,
-            "gridcolor": "#EBF0F8",
-            "showline": True,
-            "linecolor": "black",
-            "linewidth": 2.0,
-            "mirror": True,
-            "ticks": "",
-            "title": {"standoff": 15},
-            "zeroline": False,
-            "zerolinecolor": "#EBF0F8",
-            "zerolinewidth": 2,
-        },
-        "yaxis": {
-            "automargin": True,
-            "gridcolor": "#EBF0F8",
-            "showline": True,
-            "linecolor": "black",
-            "linewidth": 2.0,
-            "mirror": True,
-            "ticks": "",
-            "title": {"standoff": 15},
-            "zeroline": False,
-            "zerolinecolor": "#EBF0F8",
-            "zerolinewidth": 2,
-        },
-    },
-    data={
-        "bar": [
-            {
-                "error_x": {"color": "#2a3f5f"},
-                "error_y": {"color": "#2a3f5f"},
-                "marker": {"line": {"color": "white", "width": 0.5}},
-                "type": "bar",
-            }
-        ],
-        "barpolar": [
-            {"marker": {"line": {"color": "white", "width": 0.5}}, "type": "barpolar"}
-        ],
-        "carpet": [
-            {
-                "aaxis": {
-                    "endlinecolor": "#2a3f5f",
-                    "gridcolor": "#C8D4E3",
-                    "linecolor": "#C8D4E3",
-                    "minorgridcolor": "#C8D4E3",
-                    "startlinecolor": "#2a3f5f",
-                },
-                "baxis": {
-                    "endlinecolor": "#2a3f5f",
-                    "gridcolor": "#C8D4E3",
-                    "linecolor": "#C8D4E3",
-                    "minorgridcolor": "#C8D4E3",
-                    "startlinecolor": "#2a3f5f",
-                },
-                "type": "carpet",
-            }
-        ],
-        "choropleth": [
-            {"colorbar": {"outlinewidth": 0, "ticks": ""}, "type": "choropleth"}
-        ],
-        "contour": [
-            {
-                "colorbar": {"outlinewidth": 0, "ticks": ""},
-                "colorscale": [
-                    [0.0, "#0d0887"],
-                    [0.1111111111111111, "#46039f"],
-                    [0.2222222222222222, "#7201a8"],
-                    [0.3333333333333333, "#9c179e"],
-                    [0.4444444444444444, "#bd3786"],
-                    [0.5555555555555556, "#d8576b"],
-                    [0.6666666666666666, "#ed7953"],
-                    [0.7777777777777778, "#fb9f3a"],
-                    [0.8888888888888888, "#fdca26"],
-                    [1.0, "#f0f921"],
-                ],
-                "type": "contour",
-            }
-        ],
-        "contourcarpet": [
-            {"colorbar": {"outlinewidth": 0, "ticks": ""}, "type": "contourcarpet"}
-        ],
-        "heatmap": [
-            {
-                "colorbar": {"outlinewidth": 0, "ticks": ""},
-                "colorscale": [
-                    [0.0, "#0d0887"],
-                    [0.1111111111111111, "#46039f"],
-                    [0.2222222222222222, "#7201a8"],
-                    [0.3333333333333333, "#9c179e"],
-                    [0.4444444444444444, "#bd3786"],
-                    [0.5555555555555556, "#d8576b"],
-                    [0.6666666666666666, "#ed7953"],
-                    [0.7777777777777778, "#fb9f3a"],
-                    [0.8888888888888888, "#fdca26"],
-                    [1.0, "#f0f921"],
-                ],
-                "type": "heatmap",
-            }
-        ],
-        "heatmapgl": [
-            {
-                "colorbar": {"outlinewidth": 0, "ticks": ""},
-                "colorscale": [
-                    [0.0, "#0d0887"],
-                    [0.1111111111111111, "#46039f"],
-                    [0.2222222222222222, "#7201a8"],
-                    [0.3333333333333333, "#9c179e"],
-                    [0.4444444444444444, "#bd3786"],
-                    [0.5555555555555556, "#d8576b"],
-                    [0.6666666666666666, "#ed7953"],
-                    [0.7777777777777778, "#fb9f3a"],
-                    [0.8888888888888888, "#fdca26"],
-                    [1.0, "#f0f921"],
-                ],
-                "type": "heatmapgl",
-            }
-        ],
-        "histogram": [
-            {
-                "marker": {"colorbar": {"outlinewidth": 0, "ticks": ""}},
-                "type": "histogram",
-            }
-        ],
-        "histogram2d": [
-            {
-                "colorbar": {"outlinewidth": 0, "ticks": ""},
-                "colorscale": [
-                    [0.0, "#0d0887"],
-                    [0.1111111111111111, "#46039f"],
-                    [0.2222222222222222, "#7201a8"],
-                    [0.3333333333333333, "#9c179e"],
-                    [0.4444444444444444, "#bd3786"],
-                    [0.5555555555555556, "#d8576b"],
-                    [0.6666666666666666, "#ed7953"],
-                    [0.7777777777777778, "#fb9f3a"],
-                    [0.8888888888888888, "#fdca26"],
-                    [1.0, "#f0f921"],
-                ],
-                "type": "histogram2d",
-            }
-        ],
-        "histogram2dcontour": [
-            {
-                "colorbar": {"outlinewidth": 0, "ticks": ""},
-                "colorscale": [
-                    [0.0, "#0d0887"],
-                    [0.1111111111111111, "#46039f"],
-                    [0.2222222222222222, "#7201a8"],
-                    [0.3333333333333333, "#9c179e"],
-                    [0.4444444444444444, "#bd3786"],
-                    [0.5555555555555556, "#d8576b"],
-                    [0.6666666666666666, "#ed7953"],
-                    [0.7777777777777778, "#fb9f3a"],
-                    [0.8888888888888888, "#fdca26"],
-                    [1.0, "#f0f921"],
-                ],
-                "type": "histogram2dcontour",
-            }
-        ],
-        "mesh3d": [{"colorbar": {"outlinewidth": 0, "ticks": ""}, "type": "mesh3d"}],
-        "parcoords": [
-            {
-                "line": {"colorbar": {"outlinewidth": 0, "ticks": ""}},
-                "type": "parcoords",
-            }
-        ],
-        "pie": [{"automargin": True, "type": "pie"}],
-        "scatter": [
-            {
-                "marker": {"colorbar": {"outlinewidth": 0, "ticks": ""}},
-                "type": "scatter",
-            }
-        ],
-        "scatter3d": [
-            {
-                "line": {"colorbar": {"outlinewidth": 0, "ticks": ""}},
-                "marker": {"colorbar": {"outlinewidth": 0, "ticks": ""}},
-                "type": "scatter3d",
-            }
-        ],
-        "scattercarpet": [
-            {
-                "marker": {"colorbar": {"outlinewidth": 0, "ticks": ""}},
-                "type": "scattercarpet",
-            }
-        ],
-        "scattergeo": [
-            {
-                "marker": {"colorbar": {"outlinewidth": 0, "ticks": ""}},
-                "type": "scattergeo",
-            }
-        ],
-        "scattergl": [
-            {
-                "marker": {"colorbar": {"outlinewidth": 0, "ticks": ""}},
-                "type": "scattergl",
-            }
-        ],
-        "scattermapbox": [
-            {
-                "marker": {"colorbar": {"outlinewidth": 0, "ticks": ""}},
-                "type": "scattermapbox",
-            }
-        ],
-        "scatterpolar": [
-            {
-                "marker": {"colorbar": {"outlinewidth": 0, "ticks": ""}},
-                "type": "scatterpolar",
-            }
-        ],
-        "scatterpolargl": [
-            {
-                "marker": {"colorbar": {"outlinewidth": 0, "ticks": ""}},
-                "type": "scatterpolargl",
-            }
-        ],
-        "scatterternary": [
-            {
-                "marker": {"colorbar": {"outlinewidth": 0, "ticks": ""}},
-                "type": "scatterternary",
-            }
-        ],
-        "surface": [
-            {
-                "colorbar": {"outlinewidth": 0, "ticks": ""},
-                "colorscale": [
-                    [0.0, "#0d0887"],
-                    [0.1111111111111111, "#46039f"],
-                    [0.2222222222222222, "#7201a8"],
-                    [0.3333333333333333, "#9c179e"],
-                    [0.4444444444444444, "#bd3786"],
-                    [0.5555555555555556, "#d8576b"],
-                    [0.6666666666666666, "#ed7953"],
-                    [0.7777777777777778, "#fb9f3a"],
-                    [0.8888888888888888, "#fdca26"],
-                    [1.0, "#f0f921"],
-                ],
-                "type": "surface",
-            }
-        ],
-        "table": [
-            {
-                "cells": {"fill": {"color": "#EBF0F8"}, "line": {"color": "white"}},
-                "header": {"fill": {"color": "#C8D4E3"}, "line": {"color": "white"}},
-                "type": "table",
-            }
-        ],
-    },
-)
diff --git a/ross/point_mass.py b/ross/point_mass.py
index e2360c75..43c75112 100644
--- a/ross/point_mass.py
+++ b/ross/point_mass.py
@@ -113,6 +113,64 @@ class PointMass(Element):
             f" my={self.my:{0}.{5}}, tag={self.tag!r})"
         )
 
+    def save(self, file_name=os.getcwd()):
+        """Save a point mass 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 point mass element will be saved in.
+
+        Examples
+        --------
+        >>> point_mass = point_mass_example()
+        >>> point_mass.save()
+        """
+        data = self.get_data(Path(file_name) / "PointMass.toml")
+        data["PointMass"][str(self.n)] = {
+            "n": self.n,
+            "m": self.m,
+            "mx": self.mx,
+            "my": self.my,
+            "tag": self.tag,
+        }
+        self.dump_data(data, Path(file_name) / "PointMass.toml")
+
+    @staticmethod
+    def load(file_name=os.getcwd()):
+        """Load a list of point mass elements saved in a toml format.
+
+        It works as an auxiliary function of the load function in the Rotor class.
+
+        Parameters
+        ----------
+        file_name: str
+            The name of the file of the point mass element to be loaded.
+
+        Returns
+        -------
+        point_mass_elements: list
+            A list of point mass elements.
+
+        Examples
+        --------
+        >>> point_mass1 = point_mass_example()
+        >>> point_mass1.save(os.getcwd())
+        >>> list_of_pointmass = PointMass.load(os.getcwd())
+        >>> point_mass1 == list_of_pointmass[0]
+        True
+        """
+        point_mass_elements = []
+        with open("PointMass.toml", "r") as f:
+            point_mass_dict = toml.load(f)
+            for element in point_mass_dict["PointMass"]:
+                point_mass_elements.append(
+                    PointMass(**point_mass_dict["PointMass"][element])
+                )
+        return point_mass_elements
+
     def M(self):
         """Mass matrix for an instance of a point mass element.
 
@@ -247,7 +305,11 @@ class PointMass(Element):
         zpos, ypos = position
         radius = ypos / 8
 
-        customdata = [self.n, self.mx, self.my]
+        customdata = [
+            self.n,
+            self.mx,
+            self.my,
+        ]
         hovertemplate = (
             f"<b>PointMass Node: {customdata[0]}<b><br>"
             + f"<b>Mass (X): {customdata[1]:.3f}<b><br>"
diff --git a/ross/results.py b/ross/results.py
index 1594dc20..48ffb21d 100644
--- a/ross/results.py
+++ b/ross/results.py
@@ -2,37 +2,15 @@
 
 This module returns graphs for each type of analyses in rotor_assembly.py.
 """
-import copy
-
 import numpy as np
+import plotly.express as px
 import plotly.graph_objects as go
 import scipy.linalg as la
 from plotly.subplots import make_subplots
 from scipy import interpolate
 
-from ross.plotly_theme import tableau_colors
-
-
-class CriticalSpeedResults:
-    """Class used to store results from run_critical_speed() method.
-
-    Parameters
-    ----------
-    wn : array
-        Undamped critical speeds array.
-    wd : array
-        Undamped critical speeds array.
-    log_dec : array
-        Logarithmic decrement for each critical speed.
-    damping_ratio : array
-        Damping ratio for each critical speed.
-    """
-
-    def __init__(self, wn, wd, log_dec, damping_ratio):
-        self.wn = wn
-        self.wd = wd
-        self.log_dec = log_dec
-        self.damping_ratio = damping_ratio
+# set Plotly palette of colors
+colors1 = px.colors.qualitative.Dark24
 
 
 class ModalResults:
@@ -55,7 +33,7 @@ class ModalResults:
     wd : array
         Damped natural frequencies array.
     log_dec : array
-        Logarithmic decrement for each mode.
+        Logarithmic decrement for each .
     damping_ratio : array
         Damping ratio for each mode.
     lti : StateSpaceContinuous
@@ -479,7 +457,7 @@ class ModalResults:
 
         return xn, yn, zn, x_circles, y_circles, z_circles_pos, nn
 
-    def plot_mode_3d(self, mode=None, evec=None, fig=None, **kwargs):
+    def plot_mode3D(self, mode=None, evec=None, **kwargs):
         """Plot (3D view) the mode shapes.
 
         Parameters
@@ -489,8 +467,6 @@ class ModalResults:
             Default is None
         evec : array
             Array containing the system eigenvectors
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
@@ -501,8 +477,7 @@ class ModalResults:
         fig : Plotly graph_objects.Figure()
             The figure object with the plot.
         """
-        if fig is None:
-            fig = go.Figure()
+        fig = go.Figure()
 
         nodes = self.nodes
         kappa_mode = self.kappa_modes[mode]
@@ -514,24 +489,9 @@ class ModalResults:
                     x=zc_pos[10:, node],
                     y=xc[10:, node],
                     z=yc[10:, node],
-                    mode="lines",
-                    line=dict(color=kappa_mode[node]),
-                    name="node {}".format(node),
-                    showlegend=False,
-                    hovertemplate=(
-                        "Nodal Position: %{x:.2f}<br>"
-                        + "X - Relative Displacement: %{y:.2f}<br>"
-                        + "Y - Relative Displacement: %{z:.2f}"
-                    ),
-                )
-            )
-            fig.add_trace(
-                go.Scatter3d(
-                    x=[zc_pos[10, node]],
-                    y=[xc[10, node]],
-                    z=[yc[10, node]],
-                    mode="markers",
-                    marker=dict(color=kappa_mode[node]),
+                    mode="markers+lines",
+                    line=dict(width=4.0, color=kappa_mode[node]),
+                    marker=dict(size=5, color=kappa_mode[node]),
                     name="node {}".format(node),
                     showlegend=False,
                 )
@@ -543,7 +503,7 @@ class ModalResults:
                 y=xn,
                 z=yn,
                 mode="lines",
-                line=dict(color="black", dash="dash"),
+                line=dict(width=3.0, color="black", dash="dash"),
                 name="mode shape",
                 showlegend=False,
             )
@@ -559,43 +519,63 @@ class ModalResults:
                 y=zn_cl * 0,
                 z=zn_cl * 0,
                 mode="lines",
-                line=dict(color="black", dash="dashdot"),
+                line=dict(width=2.0, color="black", dash="dashdot"),
                 hoverinfo="none",
                 showlegend=False,
             )
         )
         fig.update_layout(
+            width=1200,
+            height=900,
             scene=dict(
+                bgcolor="white",
                 xaxis=dict(
-                    title=dict(text="<b>Rotor Length</b>"),
-                    autorange="reversed",
+                    title=dict(text="<b>Rotor Length</b>", font=dict(size=14)),
+                    tickfont=dict(size=16),
+                    range=[zn_cl0 - 0.1, zn_cl1 + 0.1],
                     nticks=5,
+                    backgroundcolor="lightgray",
+                    gridcolor="white",
+                    showspikes=False,
                 ),
                 yaxis=dict(
-                    title=dict(text="<b>Relative Displacement</b>"),
+                    title=dict(
+                        text="<b>Dimensionless deformation</b>", font=dict(size=14)
+                    ),
+                    tickfont=dict(size=16),
                     range=[-2, 2],
                     nticks=5,
+                    backgroundcolor="lightgray",
+                    gridcolor="white",
+                    showspikes=False,
                 ),
                 zaxis=dict(
-                    title=dict(text="<b>Relative Displacement</b>"),
+                    title=dict(
+                        text="<b>Dimensionless deformation</b>", font=dict(size=14)
+                    ),
+                    tickfont=dict(size=16),
                     range=[-2, 2],
                     nticks=5,
+                    backgroundcolor="lightgray",
+                    gridcolor="white",
+                    showspikes=False,
                 ),
             ),
             title=dict(
                 text=(
-                    f"<b>Mode</b> {mode + 1} | "
-                    f"<b>whirl</b>: {self.whirl_direction()[mode]} | "
-                    f"<b>ω<sub>n</sub></b> = {self.wn[mode]:.1f} rad/s | "
+                    f"<b>Mode</b> {mode + 1}<br>"
+                    f"<b>whirl</b>: {self.whirl_direction()[mode]}<br>"
+                    f"<b>ωn</b> = {self.wn[mode]:.1f} rad/s<br>"
                     f"<b>log dec</b> = {self.log_dec[mode]:.1f}"
-                )
+                ),
+                font=dict(size=14),
             ),
             **kwargs,
         )
 
         return fig
 
-    def plot_mode_2d(self, mode=None, evec=None, fig=None, **kwargs):
+    def plot_mode2D(self, mode=None, evec=None, **kwargs):
         """Plot (2D view) the mode shapes.
 
         Parameters
@@ -605,8 +585,6 @@ class ModalResults:
             Default is None
         evec : array
             Array containing the system eigenvectors
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
@@ -633,10 +611,9 @@ class ModalResults:
         zn = np.delete(zn, idx_remove)
         vn = np.delete(vn, idx_remove)
 
-        if fig is None:
-            fig = go.Figure()
+        fig = go.Figure()
 
-        colors = dict(Backward="red", Mixed="black", Forward="blue")
+        colors = dict(Backward="firebrick", Mixed="black", Forward="royalblue")
         whirl_dir = colors[self.whirl_direction()[mode]]
 
         fig.add_trace(
@@ -644,7 +621,7 @@ class ModalResults:
                 x=zn,
                 y=vn,
                 mode="lines",
-                line=dict(color=whirl_dir),
+                line=dict(width=2.5, color=whirl_dir),
                 name="mode shape",
                 showlegend=False,
             )
@@ -655,23 +632,51 @@ class ModalResults:
                 x=nodes_pos,
                 y=np.zeros(len(nodes_pos)),
                 mode="lines",
-                line=dict(color="black", dash="dashdot"),
+                line=dict(width=1.0, color="black", dash="dashdot"),
                 name="centerline",
                 hoverinfo="none",
                 showlegend=False,
             )
         )
 
-        fig.update_xaxes(title_text="<b>Rotor Length</b>")
-        fig.update_yaxes(title_text="<b>Relative Displacement</b>")
+        fig.update_xaxes(
+            title_text="<b>Rotor Length</b>",
+            title_font=dict(size=16),
+            tickfont=dict(size=14),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_yaxes(
+            title_text="<b>Non dimensional deformation</b>",
+            title_font=dict(size=16),
+            tickfont=dict(size=14),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
         fig.update_layout(
+            width=1200,
+            height=900,
+            plot_bgcolor="white",
+            legend=dict(
+                font=dict(family="sans-serif", size=14),
+                bgcolor="white",
+                bordercolor="black",
+                borderwidth=2,
+            ),
             title=dict(
                 text=(
                     f"<b>Mode</b> {mode + 1} | "
                     f"<b>whirl</b>: {self.whirl_direction()[mode]} | "
                     f"<b>ωn</b> = {self.wn[mode]:.1f} rad/s | "
                     f"<b>log dec</b> = {self.log_dec[mode]:.1f}"
-                )
+                ),
+                font=dict(size=16),
             ),
             **kwargs,
         )
@@ -708,7 +713,7 @@ class CampbellResults:
         self.log_dec = log_dec
         self.whirl_values = whirl_values
 
-    def plot(self, harmonics=[1], fig=None, **kwargs):
+    def plot(self, harmonics=[1], **kwargs):
         """Create Campbell Diagram figure using Plotly.
 
         Parameters
@@ -716,8 +721,6 @@ class CampbellResults:
         harmonics: list, optional
             List withe the harmonics to be plotted.
             The default is to plot 1x.
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
@@ -735,13 +738,13 @@ class CampbellResults:
         speed_range = self.speed_range
         log_dec_map = log_dec.flatten()
 
-        if fig is None:
-            fig = go.Figure()
+        fig = go.Figure()
 
         scatter_marker = ["triangle-up", "circle", "triangle-down"]
         for mark, whirl_dir, legend in zip(
             scatter_marker, [0.0, 0.5, 1.0], ["Foward", "Mixed", "Backward"]
         ):
+            num_frequencies = wd.shape[1]
             for i in range(num_frequencies):
                 w_i = wd[:, i]
                 whirl_i = whirl[:, i]
@@ -774,7 +777,7 @@ class CampbellResults:
                                 x=xnew[idx],
                                 y=ynew[idx],
                                 mode="markers",
-                                marker=dict(symbol="x", color="black"),
+                                marker=dict(symbol="x", size=10, color="black"),
                                 name="Crit. Speed",
                                 legendgroup="Crit. Speed",
                                 showlegend=False,
@@ -795,6 +798,7 @@ class CampbellResults:
                             y=w_i[whirl_mask],
                             marker=dict(
                                 symbol=mark,
+                                size=16,
                                 cmax=max(log_dec_map),
                                 cmin=min(log_dec_map),
                                 color=log_dec_i[whirl_mask],
@@ -815,7 +819,7 @@ class CampbellResults:
                     x=speed_range,
                     y=h * speed_range,
                     mode="lines",
-                    line=dict(dash="dashdot"),
+                    line=dict(width=2.5, color=colors1[j], dash="dashdot"),
                     name="{}x speed".format(h),
                     hoverinfo="none",
                 )
@@ -826,36 +830,57 @@ class CampbellResults:
         for mark, legend in zip(scatter_marker, legends):
             fig.add_trace(
                 go.Scatter(
-                    x=[0],
-                    y=[0],
+                    x=[-1000],
+                    y=[-1000],
                     mode="markers",
                     name=legend,
                     legendgroup=legend,
-                    marker=dict(symbol=mark, color="black"),
+                    marker=dict(symbol=mark, size=16, color="black"),
                 )
             )
 
         fig.update_xaxes(
             title_text="<b>Frequency</b>",
-            range=[np.min(speed_range), np.max(speed_range)],
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            range=[0, np.max(speed_range)],
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
         )
         fig.update_yaxes(
-            title_text="<b>Damped Natural Frequencies</b>", range=[0, 1.1 * np.max(wd)]
+            title_text="<b>Damped Natural Frequencies</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            range=[0, 1.1 * np.max(wd)],
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
         )
         fig.update_layout(
+            width=1200,
+            height=900,
+            plot_bgcolor="white",
+            hoverlabel_align="right",
             coloraxis=dict(
                 cmin=min(log_dec_map),
                 cmax=max(log_dec_map),
                 colorscale="rdbu",
-                colorbar=dict(title=dict(text="<b>Log Dec</b>", side="right")),
+                colorbar=dict(
+                    title=dict(text="<b>Log Dec</b>", side="right", font=dict(size=20)),
+                    tickfont=dict(size=16),
+                ),
             ),
             legend=dict(
                 itemsizing="constant",
+                bgcolor="white",
+                borderwidth=2,
+                font=dict(size=14),
                 orientation="h",
-                xanchor="center",
-                x=0.5,
-                yanchor="bottom",
-                y=-0.3,
             ),
             **kwargs,
         )
@@ -891,7 +916,7 @@ class FrequencyResponseResults:
         self.magnitude = magnitude
         self.phase = phase
 
-    def plot_magnitude(self, inp, out, units="mic-pk-pk", fig=None, **mag_kwargs):
+    def plot_magnitude(self, inp, out, units="mic-pk-pk", **mag_kwargs):
         """Plot frequency response (magnitude) using Plotly.
 
         This method plots the frequency response magnitude given an output and
@@ -906,8 +931,6 @@ class FrequencyResponseResults:
         units : str
             Unit system
             Default is "mic-pk-pk"
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         mag_kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
@@ -928,15 +951,20 @@ class FrequencyResponseResults:
         else:
             y_axis_label = "<b>Amplitude (dB)</b>"
 
-        if fig is None:
-            fig = go.Figure()
+        kwargs_default_values = dict(
+            width=1200, height=900, plot_bgcolor="white", hoverlabel_align="right"
+        )
+        for k, v in kwargs_default_values.items():
+            mag_kwargs.setdefault(k, v)
+
+        fig = go.Figure()
 
         fig.add_trace(
             go.Scatter(
                 x=frequency_range,
                 y=mag[inp, out, :],
                 mode="lines",
-                line=dict(color=tableau_colors["blue"]),
+                line=dict(width=3.0, color="royalblue"),
                 name="Amplitude",
                 legendgroup="Amplitude",
                 showlegend=False,
@@ -946,14 +974,38 @@ class FrequencyResponseResults:
 
         fig.update_xaxes(
             title_text="<b>Frequency</b>",
-            range=[np.min(frequency_range), np.max(frequency_range)],
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            range=[0, np.max(frequency_range)],
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_yaxes(
+            title_text=y_axis_label,
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_layout(
+            legend=dict(
+                font=dict(family="sans-serif", size=14),
+                bgcolor="white",
+                bordercolor="black",
+                borderwidth=2,
+            ),
+            **mag_kwargs,
         )
-        fig.update_yaxes(title_text=y_axis_label)
-        fig.update_layout(**mag_kwargs)
 
         return fig
 
-    def plot_phase(self, inp, out, fig=None, **phase_kwargs):
+    def plot_phase(self, inp, out, **phase_kwargs):
         """Plot frequency response (phase) using Plotly.
 
         This method plots the frequency response phase given an output and
@@ -965,8 +1017,6 @@ class FrequencyResponseResults:
             Input.
         out : int
             Output.
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         phase_kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
@@ -980,15 +1030,20 @@ class FrequencyResponseResults:
         frequency_range = self.speed_range
         phase = self.phase
 
-        if fig is None:
-            fig = go.Figure()
+        kwargs_default_values = dict(
+            width=1200, height=900, plot_bgcolor="white", hoverlabel_align="right"
+        )
+        for k, v in kwargs_default_values.items():
+            phase_kwargs.setdefault(k, v)
+
+        fig = go.Figure()
 
         fig.add_trace(
             go.Scatter(
                 x=frequency_range,
                 y=phase[inp, out, :],
                 mode="lines",
-                line=dict(color=tableau_colors["blue"]),
+                line=dict(width=3.0, color="royalblue"),
                 name="Phase",
                 legendgroup="Phase",
                 showlegend=False,
@@ -998,14 +1053,38 @@ class FrequencyResponseResults:
 
         fig.update_xaxes(
             title_text="<b>Frequency</b>",
-            range=[np.min(frequency_range), np.max(frequency_range)],
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            range=[0, np.max(frequency_range)],
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_yaxes(
+            title_text="<b>Phase</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_layout(
+            legend=dict(
+                font=dict(family="sans-serif", size=14),
+                bgcolor="white",
+                bordercolor="black",
+                borderwidth=2,
+            ),
+            **phase_kwargs,
         )
-        fig.update_yaxes(title_text="<b>Phase</b>")
-        fig.update_layout(**phase_kwargs)
 
         return fig
 
-    def plot_polar_bode(self, inp, out, units="mic-pk-pk", fig=None, **polar_kwargs):
+    def plot_polar_bode(self, inp, out, units="mic-pk-pk", **polar_kwargs):
         """Plot frequency response (polar) using Plotly.
 
         This method plots the frequency response (polar graph) given an output and
@@ -1020,8 +1099,6 @@ class FrequencyResponseResults:
         units : str
             Magnitude unit system.
             Default is "mic-pk-pk"
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         polar_kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
@@ -1043,8 +1120,13 @@ class FrequencyResponseResults:
         else:
             r_axis_label = "<b>Amplitude (dB)</b>"
 
-        if fig is None:
-            fig = go.Figure()
+        kwargs_default_values = dict(
+            width=1200, height=900, polar_bgcolor="white", hoverlabel_align="right"
+        )
+        for k, v in kwargs_default_values.items():
+            polar_kwargs.setdefault(k, v)
+
+        fig = go.Figure()
 
         fig.add_trace(
             go.Scatterpolar(
@@ -1053,8 +1135,8 @@ class FrequencyResponseResults:
                 customdata=frequency_range,
                 thetaunit="radians",
                 mode="lines+markers",
-                marker=dict(color=tableau_colors["blue"]),
-                line=dict(color=tableau_colors["blue"]),
+                marker=dict(size=8, color="royalblue"),
+                line=dict(width=3.0, color="royalblue"),
                 name="Polar_plot",
                 legendgroup="Polar",
                 showlegend=False,
@@ -1068,7 +1150,18 @@ class FrequencyResponseResults:
 
         fig.update_layout(
             polar=dict(
-                radialaxis=dict(title=dict(text=r_axis_label), exponentformat="power")
+                radialaxis=dict(
+                    title=dict(text=r_axis_label, font=dict(size=14)),
+                    tickfont=dict(size=14),
+                    gridcolor="lightgray",
+                    exponentformat="power",
+                ),
+                angularaxis=dict(
+                    tickfont=dict(size=14),
+                    gridcolor="lightgray",
+                    linecolor="black",
+                    linewidth=2.5,
+                ),
             ),
             **polar_kwargs,
         )
@@ -1080,10 +1173,10 @@ class FrequencyResponseResults:
         inp,
         out,
         units="mic-pk-pk",
-        mag_kwargs=None,
-        phase_kwargs=None,
-        polar_kwargs=None,
-        subplot_kwargs=None,
+        mag_kwargs={},
+        phase_kwargs={},
+        polar_kwargs={},
+        subplot_kwargs={},
     ):
         """Plot frequency response.
 
@@ -1132,10 +1225,15 @@ class FrequencyResponseResults:
             Plotly figure with Amplitude vs Frequency and Phase vs Frequency and
             polar Amplitude vs Phase plots.
         """
-        mag_kwargs = {} if mag_kwargs is None else copy.copy(mag_kwargs)
-        phase_kwargs = {} if phase_kwargs is None else copy.copy(phase_kwargs)
-        polar_kwargs = {} if polar_kwargs is None else copy.copy(polar_kwargs)
-        subplot_kwargs = {} if subplot_kwargs is None else copy.copy(subplot_kwargs)
+        kwargs_default_values = dict(
+            width=1800,
+            height=900,
+            polar_bgcolor="white",
+            plot_bgcolor="white",
+            hoverlabel_align="right",
+        )
+        for k, v in kwargs_default_values.items():
+            subplot_kwargs.setdefault(k, v)
 
         fig0 = self.plot_magnitude(inp, out, units, **mag_kwargs)
         fig1 = self.plot_phase(inp, out, **phase_kwargs)
@@ -1194,7 +1292,7 @@ class ForcedResponseResults:
     """
 
     def __init__(
-        self, rotor, forced_resp, speed_range, magnitude, phase, unbalance=None
+        self, rotor, forced_resp, speed_range, magnitude, phase, unbalance=None,
     ):
         self.rotor = rotor
         self.forced_resp = forced_resp
@@ -1203,7 +1301,7 @@ class ForcedResponseResults:
         self.phase = phase
         self.unbalance = unbalance
 
-    def plot_magnitude(self, dof, units="m", fig=None, **kwargs):
+    def plot_magnitude(self, dof, units="m", **kwargs):
         """Plot forced response (magnitude) using Plotly.
 
         Parameters
@@ -1213,8 +1311,6 @@ class ForcedResponseResults:
         units : str
             Units to plot the magnitude ('m' or 'mic-pk-pk')
             Default is 'm'
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
@@ -1234,15 +1330,20 @@ class ForcedResponseResults:
             mag = 2 * mag * 1e6
             y_axis_label = "<b>Amplitude (μ pk-pk)</b>"
 
-        if fig is None:
-            fig = go.Figure()
+        kwargs_default_values = dict(
+            width=1200, height=900, plot_bgcolor="white", hoverlabel_align="right"
+        )
+        for k, v in kwargs_default_values.items():
+            kwargs.setdefault(k, v)
+
+        fig = go.Figure()
 
         fig.add_trace(
             go.Scatter(
                 x=frequency_range,
                 y=mag[dof],
                 mode="lines",
-                line=dict(color=tableau_colors["blue"]),
+                line=dict(width=3.0, color="royalblue"),
                 name="Amplitude",
                 legendgroup="Amplitude",
                 showlegend=False,
@@ -1252,22 +1353,45 @@ class ForcedResponseResults:
 
         fig.update_xaxes(
             title_text="<b>Frequency</b>",
-            range=[np.min(frequency_range), np.max(frequency_range)],
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            range=[0, np.max(frequency_range)],
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_yaxes(
+            title_text=y_axis_label,
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            exponentformat="power",
+            mirror=True,
+        )
+        fig.update_layout(
+            legend=dict(
+                font=dict(family="sans-serif", size=14),
+                bgcolor="white",
+                bordercolor="black",
+                borderwidth=2,
+            ),
+            **kwargs,
         )
-        fig.update_yaxes(title_text=y_axis_label, exponentformat="power")
-        fig.update_layout(**kwargs)
 
         return fig
 
-    def plot_phase(self, dof, fig=None, **kwargs):
+    def plot_phase(self, dof, **kwargs):
         """Plot forced response (phase) using Plotly.
 
         Parameters
         ----------
         dof : int
             Degree of freedom.
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
@@ -1281,15 +1405,20 @@ class ForcedResponseResults:
         frequency_range = self.speed_range
         phase = self.phase
 
-        if fig is None:
-            fig = go.Figure()
+        kwargs_default_values = dict(
+            width=1200, height=900, plot_bgcolor="white", hoverlabel_align="right"
+        )
+        for k, v in kwargs_default_values.items():
+            kwargs.setdefault(k, v)
+
+        fig = go.Figure()
 
         fig.add_trace(
             go.Scatter(
                 x=frequency_range,
                 y=phase[dof],
                 mode="lines",
-                line=dict(color=tableau_colors["blue"]),
+                line=dict(width=3.0, color="royalblue"),
                 name="Phase",
                 legendgroup="Phase",
                 showlegend=False,
@@ -1299,14 +1428,38 @@ class ForcedResponseResults:
 
         fig.update_xaxes(
             title_text="<b>Frequency</b>",
-            range=[np.min(frequency_range), np.max(frequency_range)],
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            range=[0, np.max(frequency_range)],
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_yaxes(
+            title_text="<b>Phase Angle</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_layout(
+            legend=dict(
+                font=dict(family="sans-serif", size=14),
+                bgcolor="white",
+                bordercolor="black",
+                borderwidth=2,
+            ),
+            **kwargs,
         )
-        fig.update_yaxes(title_text="<b>Phase Angle</b>")
-        fig.update_layout(**kwargs)
 
         return fig
 
-    def plot_polar_bode(self, dof, units="m", fig=None, **kwargs):
+    def plot_polar_bode(self, dof, units="mic-pk-pk", **kwargs):
         """Plot polar forced response using Plotly.
 
         Parameters
@@ -1316,8 +1469,6 @@ class ForcedResponseResults:
         units : str
             Magnitude unit system.
             Default is "mic-pk-pk"
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
@@ -1339,8 +1490,13 @@ class ForcedResponseResults:
         else:
             r_axis_label = "<b>Amplitude (dB)</b>"
 
-        if fig is None:
-            fig = go.Figure()
+        kwargs_default_values = dict(
+            width=1200, height=900, polar_bgcolor="white", hoverlabel_align="right"
+        )
+        for k, v in kwargs_default_values.items():
+            kwargs.setdefault(k, v)
+
+        fig = go.Figure()
 
         fig.add_trace(
             go.Scatterpolar(
@@ -1349,8 +1505,8 @@ class ForcedResponseResults:
                 customdata=frequency_range,
                 thetaunit="radians",
                 mode="lines+markers",
-                marker=dict(color=tableau_colors["blue"]),
-                line=dict(color=tableau_colors["blue"]),
+                marker=dict(size=8, color="royalblue"),
+                line=dict(width=3.0, color="royalblue"),
                 name="Polar_plot",
                 legendgroup="Polar",
                 showlegend=False,
@@ -1364,7 +1520,18 @@ class ForcedResponseResults:
 
         fig.update_layout(
             polar=dict(
-                radialaxis=dict(title_text=r_axis_label, exponentformat="power")
+                radialaxis=dict(
+                    title_text=r_axis_label,
+                    title_font=dict(family="Arial", size=14),
+                    gridcolor="lightgray",
+                    exponentformat="power",
+                ),
+                angularaxis=dict(
+                    tickfont=dict(size=14),
+                    gridcolor="lightgray",
+                    linecolor="black",
+                    linewidth=2.5,
+                ),
             ),
             **kwargs,
         )
@@ -1374,11 +1541,11 @@ class ForcedResponseResults:
     def plot(
         self,
         dof,
-        units="m",
-        mag_kwargs=None,
-        phase_kwargs=None,
-        polar_kwargs=None,
-        subplot_kwargs=None,
+        units="mic-pk-pk",
+        mag_kwargs={},
+        phase_kwargs={},
+        polar_kwargs={},
+        subplot_kwargs={},
     ):
         """Plot forced response.
 
@@ -1422,14 +1589,19 @@ class ForcedResponseResults:
             Plotly figure with Amplitude vs Frequency and Phase vs Frequency and
             polar Amplitude vs Phase plots.
         """
-        mag_kwargs = {} if mag_kwargs is None else copy.copy(mag_kwargs)
-        phase_kwargs = {} if phase_kwargs is None else copy.copy(phase_kwargs)
-        polar_kwargs = {} if polar_kwargs is None else copy.copy(polar_kwargs)
-        subplot_kwargs = {} if subplot_kwargs is None else copy.copy(subplot_kwargs)
+        kwargs_default_values = dict(
+            width=1800,
+            height=900,
+            polar_bgcolor="white",
+            plot_bgcolor="white",
+            hoverlabel_align="right",
+        )
+        for k, v in kwargs_default_values.items():
+            subplot_kwargs.setdefault(k, v)
 
-        fig0 = self.plot_magnitude(dof, units, **mag_kwargs)
+        fig0 = self.plot_magnitude(dof, **mag_kwargs)
         fig1 = self.plot_phase(dof, **phase_kwargs)
-        fig2 = self.plot_polar_bode(dof, units, **polar_kwargs)
+        fig2 = self.plot_polar_bode(dof, **polar_kwargs)
 
         subplots = make_subplots(
             rows=2, cols=2, specs=[[{}, {"type": "polar", "rowspan": 2}], [{}, None]]
@@ -1564,7 +1736,7 @@ class ForcedResponseResults:
 
         return Mx, My
 
-    def plot_deflected_shape_2d(self, speed, units="mic-pk-pk", fig=None, **kwargs):
+    def plot_deflected_shape_2d(self, speed, units="mic-pk-pk", **kwargs):
         """Plot the 2D deflected shape diagram.
 
         Parameters
@@ -1579,8 +1751,6 @@ class ForcedResponseResults:
                 - "mic-pk-pk" : microns peak to peak
                 - "db" : decibels
             Default is "mic-pk-pk".
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the deflected shape
             plot layout only (e.g. width=1000, height=800, ...).
@@ -1594,16 +1764,22 @@ class ForcedResponseResults:
         if not any(np.isclose(self.speed_range, speed, atol=1e-6)):
             raise ValueError("No data available for this speed value.")
 
+        kwargs_default_values = dict(
+            width=1200, height=900, plot_bgcolor="white", hoverlabel_align="right",
+        )
+        for k, v in kwargs_default_values.items():
+            kwargs.setdefault(k, v)
+
         nodes_pos = self.rotor.nodes_pos
         maj_vect = self._calculate_major_axis(speed=speed)
 
-        if fig is None:
-            fig = go.Figure()
+        fig = go.Figure()
         fig.add_trace(
             go.Scatter(
                 x=nodes_pos,
                 y=maj_vect[4].real,
                 mode="lines",
+                line=dict(width=6.0, color="royalblue"),
                 name="Major Axis",
                 legendgroup="Major_Axis_2d",
                 showlegend=False,
@@ -1618,21 +1794,45 @@ class ForcedResponseResults:
                 x=nodes_pos,
                 y=np.zeros(len(nodes_pos)),
                 mode="lines",
-                line=dict(color="black", dash="dashdot"),
+                line=dict(width=2.0, color="black", dash="dashdot"),
                 showlegend=False,
                 hoverinfo="none",
             )
         )
 
-        fig.update_xaxes(title_text="<b>Rotor Length</b>")
-        fig.update_yaxes(title_text="<b>Major Axis Absolute Amplitude</b>")
-        fig.update_layout(**kwargs)
+        fig.update_xaxes(
+            title_text="<b>Rotor Length</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_yaxes(
+            title_text="<b>Major Axis Absolute Amplitude</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_layout(
+            legend=dict(
+                font=dict(family="sans-serif", size=14),
+                bgcolor="white",
+                bordercolor="black",
+                borderwidth=2,
+            ),
+            **kwargs,
+        )
 
         return fig
 
-    def plot_deflected_shape_3d(
-        self, speed, samples=101, units="mic-pk-pk", fig=None, **kwargs
-    ):
+    def plot_deflected_shape_3d(self, speed, samples=101, units="mic-pk-pk", **kwargs):
         """Plot the 3D deflected shape diagram.
 
         Parameters
@@ -1650,8 +1850,6 @@ class ForcedResponseResults:
                 - "mic-pk-pk" : microns peak to peak
                 - "db" : decibels
             Default is "mic-pk-pk".
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the deflected shape
             plot layout only (e.g. width=1000, height=800, ...).
@@ -1665,6 +1863,10 @@ class ForcedResponseResults:
         if not any(np.isclose(self.speed_range, speed, atol=1e-6)):
             raise ValueError("No data available for this speed value.")
 
+        kwargs_default_values = dict(hoverlabel_align="right")
+        for k, v in kwargs_default_values.items():
+            kwargs.setdefault(k, v)
+
         mag = self.magnitude
         phase = self.phase
         ub = self.unbalance
@@ -1678,8 +1880,7 @@ class ForcedResponseResults:
 
         x_pos = np.repeat(nodes_pos, t.size).reshape(len(nodes_pos), t.size)
 
-        if fig is None:
-            fig = go.Figure()
+        fig = go.Figure()
         for i, n in enumerate(nodes):
             dofx = number_dof * n
             dofy = number_dof * n + 1
@@ -1694,7 +1895,7 @@ class ForcedResponseResults:
                     y=y,
                     z=z,
                     mode="lines",
-                    line=dict(color="royalblue"),
+                    line=dict(width=6.0, color="royalblue"),
                     name="Orbit",
                     legendgroup="Orbit",
                     showlegend=False,
@@ -1715,8 +1916,8 @@ class ForcedResponseResults:
                 y=np.real(maj_vect[3]),
                 z=np.imag(maj_vect[3]),
                 mode="lines+markers",
-                marker=dict(color="black"),
-                line=dict(color="black", dash="dashdot"),
+                marker=dict(size=4.0, color="black"),
+                line=dict(width=6.0, color="black", dash="dashdot"),
                 name="Major Axis",
                 legendgroup="Major_Axis",
                 showlegend=True,
@@ -1736,7 +1937,7 @@ class ForcedResponseResults:
                 y=line,
                 z=line,
                 mode="lines",
-                line=dict(color="black", dash="dashdot"),
+                line=dict(width=2.0, color="black", dash="dashdot"),
                 showlegend=False,
                 hoverinfo="none",
             )
@@ -1751,7 +1952,7 @@ class ForcedResponseResults:
                     y=[0, np.amax(np.real(maj_vect[4])) / 2 * np.cos(p)],
                     z=[0, np.amax(np.real(maj_vect[4])) / 2 * np.sin(p)],
                     mode="lines",
-                    line=dict(color="firebrick"),
+                    line=dict(width=6.0, color="firebrick"),
                     legendgroup="Unbalance",
                     hoverinfo="none",
                     showlegend=False,
@@ -1763,7 +1964,7 @@ class ForcedResponseResults:
                     y=[np.amax(np.real(maj_vect[4])) / 2 * np.cos(p)],
                     z=[np.amax(np.real(maj_vect[4])) / 2 * np.sin(p)],
                     mode="markers",
-                    marker=dict(symbol="diamond", color="firebrick"),
+                    marker=dict(symbol="diamond", size=5.0, color="firebrick"),
                     name="Unbalance",
                     legendgroup="Unbalance",
                     showlegend=True if i == 0 else False,
@@ -1775,19 +1976,51 @@ class ForcedResponseResults:
             i += 1
 
         fig.update_layout(
+            width=1200,
+            height=900,
             scene=dict(
                 bgcolor="white",
-                xaxis=dict(title=dict(text="<b>Rotor Length</b>")),
-                yaxis=dict(title=dict(text="<b>Amplitude - X</b>")),
-                zaxis=dict(title=dict(text="<b>Amplitude - Y</b>")),
+                xaxis=dict(
+                    title=dict(text="<b>Rotor Length</b>", font=dict(size=14)),
+                    tickfont=dict(size=12),
+                    nticks=5,
+                    backgroundcolor="lightgray",
+                    gridcolor="white",
+                    showspikes=False,
+                ),
+                yaxis=dict(
+                    title=dict(text="<b>Amplitude - X</b>", font=dict(size=14)),
+                    tickfont=dict(size=12),
+                    nticks=5,
+                    backgroundcolor="lightgray",
+                    gridcolor="white",
+                    showspikes=False,
+                ),
+                zaxis=dict(
+                    title=dict(text="<b>Amplitude - Y</b>", font=dict(size=14)),
+                    tickfont=dict(size=12),
+                    nticks=5,
+                    backgroundcolor="lightgray",
+                    gridcolor="white",
+                    showspikes=False,
+                ),
+            ),
+            title=dict(
+                text=(f"<b>Deflected Shape</b><br>" f"<b>Speed = {speed}</b>"),
+                font=dict(size=18),
+            ),
+            legend=dict(
+                font=dict(family="sans-serif", size=14),
+                bgcolor="white",
+                bordercolor="black",
+                borderwidth=2,
             ),
-            title=dict(text=(f"<b>Deflected Shape</b><br>" f"<b>Speed = {speed}</b>")),
             **kwargs,
         )
 
         return fig
 
-    def plot_bending_moment(self, speed, units="mic-pk-pk", fig=None, **kwargs):
+    def plot_bending_moment(self, speed, units="mic-pk-pk", **kwargs):
         """Plot the bending moment diagram.
 
         Parameters
@@ -1815,51 +2048,59 @@ class ForcedResponseResults:
         if not any(np.isclose(self.speed_range, speed, atol=1e-6)):
             raise ValueError("No data available for this speed value.")
 
+        kwargs_default_values = dict(
+            width=1200, height=900, plot_bgcolor="white", hoverlabel_align="right",
+        )
+        for k, v in kwargs_default_values.items():
+            kwargs.setdefault(k, v)
+
         Mx, My = self._calculate_bending_moment(speed=speed)
         Mr = np.sqrt(Mx ** 2 + My ** 2)
 
         nodes_pos = self.rotor.nodes_pos
 
-        if fig is None:
-            fig = go.Figure()
+        fig = go.Figure()
         fig.add_trace(
             go.Scatter(
                 x=nodes_pos,
                 y=Mx,
                 mode="lines",
+                line=dict(width=6.0, color=colors1[2]),
                 name="Bending Moment (X dir.)",
                 legendgroup="Mx",
                 showlegend=True,
                 hovertemplate=(
                     "<b>Nodal Position: %{x:.2f}</b><br>" + "<b>Mx: %{y:.2e}</b>"
                 ),
-            )
+            ),
         )
         fig.add_trace(
             go.Scatter(
                 x=nodes_pos,
                 y=My,
                 mode="lines",
+                line=dict(width=6.0, color=colors1[6]),
                 name="Bending Moment (Y dir.)",
                 legendgroup="My",
                 showlegend=True,
                 hovertemplate=(
                     "<b>Nodal Position: %{x:.2f}</b><br>" + "<b>My: %{y:.2e}</b>"
                 ),
-            )
+            ),
         )
         fig.add_trace(
             go.Scatter(
                 x=nodes_pos,
                 y=Mr,
                 mode="lines",
+                line=dict(width=6.0, color=colors1[7]),
                 name="Bending Moment (abs)",
                 legendgroup="Mr",
                 showlegend=True,
                 hovertemplate=(
                     "<b>Nodal Position: %{x:.2f}</b><br>" + "<b>Mr: %{y:.2e}</b>"
                 ),
-            )
+            ),
         )
 
         # plot center line
@@ -1868,15 +2109,41 @@ class ForcedResponseResults:
                 x=nodes_pos,
                 y=np.zeros_like(nodes_pos),
                 mode="lines",
-                line=dict(color="black", dash="dashdot"),
+                line=dict(width=3.0, color="black", dash="dashdot"),
                 showlegend=False,
                 hoverinfo="none",
             )
         )
 
-        fig.update_xaxes(title_text="<b>Rotor Length</b>")
-        fig.update_yaxes(title_text="<b>Bending Moment</b>")
-        fig.update_layout(**kwargs)
+        fig.update_xaxes(
+            title_text="<b>Rotor Length</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_yaxes(
+            title_text="<b>Bending Moment</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_layout(
+            legend=dict(
+                font=dict(family="sans-serif", size=14),
+                bgcolor="white",
+                bordercolor="black",
+                borderwidth=2,
+            ),
+            **kwargs,
+        )
 
         return fig
 
@@ -1885,10 +2152,10 @@ class ForcedResponseResults:
         speed,
         samples=101,
         units="mic-pk-pk",
-        shape2d_kwargs=None,
-        shape3d_kwargs=None,
-        bm_kwargs=None,
-        subplot_kwargs=None,
+        shape2d_kwargs={},
+        shape3d_kwargs={},
+        bm_kwargs={},
+        subplot_kwargs={},
     ):
         """Plot deflected shape diagrams.
 
@@ -1936,10 +2203,11 @@ class ForcedResponseResults:
             Plotly figure with Amplitude vs Frequency and Phase vs Frequency and
             polar Amplitude vs Phase plots.
         """
-        shape2d_kwargs = {} if shape2d_kwargs is None else copy.copy(shape2d_kwargs)
-        shape3d_kwargs = {} if shape3d_kwargs is None else copy.copy(shape3d_kwargs)
-        bm_kwargs = {} if bm_kwargs is None else copy.copy(bm_kwargs)
-        subplot_kwargs = {} if subplot_kwargs is None else copy.copy(subplot_kwargs)
+        kwargs_default_values = dict(
+            width=1800, height=900, plot_bgcolor="white", hoverlabel_align="right",
+        )
+        for k, v in kwargs_default_values.items():
+            subplot_kwargs.setdefault(k, v)
 
         fig0 = self.plot_deflected_shape_2d(speed, units, **shape2d_kwargs)
         fig1 = self.plot_deflected_shape_3d(speed, samples, units, **shape3d_kwargs)
@@ -2029,7 +2297,7 @@ class StaticResults:
         self.nodes_pos = nodes_pos
         self.Vx_axis = Vx_axis
 
-    def plot_deformation(self, fig=None, **kwargs):
+    def plot_deformation(self, **kwargs):
         """Plot the shaft static deformation.
 
         This method plots:
@@ -2037,8 +2305,6 @@ class StaticResults:
 
         Parameters
         ----------
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
@@ -2049,8 +2315,13 @@ class StaticResults:
         fig : Plotly graph_objects.Figure()
             The figure object with the plot.
         """
-        if fig is None:
-            fig = go.Figure()
+        kwargs_default_values = dict(
+            width=1200, height=900, plot_bgcolor="white", hoverlabel_align="right"
+        )
+        for k, v in kwargs_default_values.items():
+            kwargs.setdefault(k, v)
+
+        fig = go.Figure()
 
         shaft_end = max([sublist[-1] for sublist in self.nodes_pos])
 
@@ -2060,7 +2331,7 @@ class StaticResults:
                 x=[-0.01 * shaft_end, 1.01 * shaft_end],
                 y=[0, 0],
                 mode="lines",
-                line=dict(color="black", dash="dashdot"),
+                line=dict(width=3.0, color="black", dash="dashdot"),
                 showlegend=False,
                 hoverinfo="none",
             )
@@ -2082,7 +2353,8 @@ class StaticResults:
                     x=xnew,
                     y=ynew,
                     mode="lines",
-                    name=f"Shaft {count}",
+                    line=dict(width=5.0, color=colors1[count]),
+                    name="Shaft {}".format(count),
                     showlegend=True,
                     hovertemplate=(
                         "Shaft lengh: %{x:.2f}<br>" + "Displacement: %{y:.2e}"
@@ -2091,19 +2363,44 @@ class StaticResults:
             )
             count += 1
 
-        fig.update_xaxes(title_text="<b>Shaft Length</b>")
-        fig.update_yaxes(title_text="<b>Deformation</b>")
-        fig.update_layout(title=dict(text="<b>Static Deformation</b>"), **kwargs)
+        fig.update_xaxes(
+            title_text="<b>Shaft Length</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_yaxes(
+            title_text="<b>Deformation</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_layout(
+            title=dict(text="<b>Static Deformation</b>", font=dict(size=16)),
+            legend=dict(
+                font=dict(family="sans-serif", size=14),
+                bgcolor="white",
+                bordercolor="black",
+                borderwidth=2,
+            ),
+            **kwargs,
+        )
 
         return fig
 
-    def plot_free_body_diagram(self, fig=None, **kwargs):
+    def plot_free_body_diagram(self, **kwargs):
         """Plot the rotor free-body diagram.
 
         Parameters
         ----------
-        subplots : Plotly graph_objects.make_subplots()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. plot_bgcolor="white", ...).
@@ -2114,17 +2411,20 @@ class StaticResults:
         subplots : Plotly graph_objects.make_subplots()
             The figure object with the plot.
         """
+        kwargs_default_values = dict(width=1200, height=900, plot_bgcolor="white")
+        for k, v in kwargs_default_values.items():
+            kwargs.setdefault(k, v)
+
         cols = 1 if len(self.nodes_pos) < 2 else 2
         rows = len(self.nodes_pos) // 2 + len(self.nodes_pos) % 2
-        if fig is None:
-            fig = make_subplots(
-                rows=rows,
-                cols=cols,
-                subplot_titles=[
-                    "<b>Free-Body Diagram - Shaft {}</b>".format(j)
-                    for j in range(len(self.nodes_pos))
-                ],
-            )
+        fig = make_subplots(
+            rows=rows,
+            cols=cols,
+            subplot_titles=[
+                "<b>Free-Body Diagram - Shaft {}</b>".format(j)
+                for j in range(len(self.nodes_pos))
+            ],
+        )
         j = 0
         y_start = 5.0
         for nodes_pos, nodes in zip(self.nodes_pos, self.nodes):
@@ -2136,7 +2436,7 @@ class StaticResults:
                     x=nodes_pos,
                     y=np.zeros(len(nodes_pos)),
                     mode="lines",
-                    line=dict(color="black"),
+                    line=dict(width=10.0, color="black"),
                     hoverinfo="none",
                     showlegend=False,
                 ),
@@ -2148,7 +2448,7 @@ class StaticResults:
                     x=nodes_pos,
                     y=[y_start] * len(nodes_pos),
                     mode="lines",
-                    line=dict(color="black"),
+                    line=dict(width=3.0, color="black"),
                     hoverinfo="none",
                     showlegend=False,
                 ),
@@ -2185,6 +2485,7 @@ class StaticResults:
                 xshift=125,
                 yshift=20,
                 text="<b>Shaft weight = {}N</b>".format(text),
+                font=dict(size=20),
                 align="right",
                 showarrow=False,
             )
@@ -2203,6 +2504,7 @@ class StaticResults:
                         ayref="y{}".format(j + 1),
                         text="<b>Fb = {}N</b>".format(text),
                         textangle=90,
+                        font=dict(size=20),
                         showarrow=True,
                         arrowhead=2,
                         arrowsize=1,
@@ -2227,6 +2529,7 @@ class StaticResults:
                         ayref="y{}".format(j + 1),
                         text="<b>Fd = {}N</b>".format(text),
                         textangle=270,
+                        font=dict(size=20),
                         showarrow=True,
                         arrowhead=2,
                         arrowsize=1,
@@ -2238,7 +2541,17 @@ class StaticResults:
                         col=col,
                     )
 
-            fig.update_xaxes(title_text="<b>Shaft Length</b>", row=row, col=col)
+            fig.update_xaxes(
+                title_text="<b>Shaft Length</b>",
+                title_font=dict(family="Arial", size=20),
+                tickfont=dict(size=16),
+                showgrid=False,
+                showline=True,
+                linewidth=2.5,
+                linecolor="black",
+                row=row,
+                col=col,
+            )
             fig.update_yaxes(
                 visible=False, gridcolor="lightgray", showline=False, row=row, col=col
             )
@@ -2248,7 +2561,7 @@ class StaticResults:
 
         return fig
 
-    def plot_shearing_force(self, fig=None, **kwargs):
+    def plot_shearing_force(self, **kwargs):
         """Plot the rotor shearing force diagram.
 
         This method plots:
@@ -2256,8 +2569,6 @@ class StaticResults:
 
         Parameters
         ----------
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
@@ -2268,8 +2579,13 @@ class StaticResults:
         fig : Plotly graph_objects.Figure()
             The figure object with the plot.
         """
-        if fig is None:
-            fig = go.Figure()
+        kwargs_default_values = dict(
+            width=1200, height=900, plot_bgcolor="white", hoverlabel_align="right"
+        )
+        for k, v in kwargs_default_values.items():
+            kwargs.setdefault(k, v)
+
+        fig = go.Figure()
 
         shaft_end = max([sublist[-1] for sublist in self.nodes_pos])
 
@@ -2279,7 +2595,7 @@ class StaticResults:
                 x=[-0.1 * shaft_end, 1.1 * shaft_end],
                 y=[0, 0],
                 mode="lines",
-                line=dict(color="black", dash="dashdot"),
+                line=dict(width=3.0, color="black", dash="dashdot"),
                 showlegend=False,
                 hoverinfo="none",
             )
@@ -2292,8 +2608,9 @@ class StaticResults:
                     x=Vx_axis,
                     y=Vx,
                     mode="lines",
-                    name=f"Shaft {j}",
-                    legendgroup=f"Shaft {j}",
+                    line=dict(width=5.0, color=colors1[j]),
+                    name="Shaft {}".format(j),
+                    legendgroup="Shaft {}".format(j),
                     showlegend=True,
                     hovertemplate=(
                         "Shaft lengh: %{x:.2f}<br>" + "Shearing Force: %{y:.2f}"
@@ -2303,31 +2620,57 @@ class StaticResults:
             j += 1
 
         fig.update_xaxes(
-            title_text="<b>Shaft Length</b>", range=[-0.1 * shaft_end, 1.1 * shaft_end]
+            title_text="<b>Shaft Length</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            range=[-0.1 * shaft_end, 1.1 * shaft_end],
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_yaxes(
+            title_text="<b>Force</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_layout(
+            title=dict(text="<b>Shearing Force Diagram</b>", font=dict(size=16)),
+            legend=dict(
+                font=dict(family="sans-serif", size=14),
+                bgcolor="white",
+                bordercolor="black",
+                borderwidth=2,
+            ),
+            **kwargs,
         )
-        fig.update_yaxes(title_text="<b>Force</b>")
-        fig.update_layout(title=dict(text="<b>Shearing Force Diagram</b>"), **kwargs)
 
         return fig
 
-    def plot_bending_moment(self, fig=None, **kwargs):
+    def plot_bending_moment(self, **kwargs):
         """Plot the rotor bending moment diagram.
 
         This method plots:
             bending moment diagram.
 
-        Parameters
-        ----------
-        fig : Plotly graph_objects.Figure()
-            Plotly figure with the bending moment diagram plot
-
         Returns
         -------
         fig : Plotly graph_objects.Figure()
             Plotly figure with the bending moment diagram plot
         """
-        if fig is None:
-            fig = go.Figure()
+        kwargs_default_values = dict(
+            width=1200, height=900, plot_bgcolor="white", hoverlabel_align="right",
+        )
+        for k, v in kwargs_default_values.items():
+            kwargs.setdefault(k, v)
+
+        fig = go.Figure()
 
         shaft_end = max([sublist[-1] for sublist in self.nodes_pos])
 
@@ -2337,7 +2680,7 @@ class StaticResults:
                 x=[-0.1 * shaft_end, 1.1 * shaft_end],
                 y=[0, 0],
                 mode="lines",
-                line=dict(color="black", dash="dashdot"),
+                line=dict(width=3.0, color="black", dash="dashdot"),
                 showlegend=False,
                 hoverinfo="none",
             )
@@ -2350,8 +2693,9 @@ class StaticResults:
                     x=nodes_pos,
                     y=Bm,
                     mode="lines",
-                    name=f"Shaft {j}",
-                    legendgroup=f"Shaft {j}",
+                    line=dict(width=5.0, color=colors1[j]),
+                    name="Shaft {}".format(j),
+                    legendgroup="Shaft {}".format(j),
                     showlegend=True,
                     hovertemplate=(
                         "Shaft lengh: %{x:.2f}<br>" + "Bending Moment: %{y:.2f}"
@@ -2360,9 +2704,37 @@ class StaticResults:
             )
             j += 1
 
-        fig.update_xaxes(title_text="<b>Shaft Length</b>")
-        fig.update_yaxes(title_text="<b>Bending Moment</b>")
-        fig.update_layout(title=dict(text="<b>Bending Moment Diagram</b>"), **kwargs)
+        fig.update_xaxes(
+            title_text="<b>Shaft Length</b>",
+            title_font=dict(family="Arial", size=20),
+            range=[-0.1 * shaft_end, 1.1 * shaft_end],
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_yaxes(
+            title_text="<b>Bending Moment</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_layout(
+            title=dict(text="<b>Bending Moment Diagram</b>", font=dict(size=16)),
+            legend=dict(
+                font=dict(family="sans-serif", size=14),
+                bgcolor="white",
+                bordercolor="black",
+                borderwidth=2,
+            ),
+            **kwargs,
+        )
 
         return fig
 
@@ -2603,7 +2975,7 @@ class ConvergenceResults:
         self.eigv_arr = eigv_arr
         self.error_arr = error_arr
 
-    def plot(self, fig=None, **kwargs):
+    def plot(self, **kwargs):
         """Plot convergence results.
 
         This method plots:
@@ -2612,8 +2984,6 @@ class ConvergenceResults:
 
         Parameters
         ----------
-        fig : Plotly graph_objects.make_subplots()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
@@ -2624,15 +2994,19 @@ class ConvergenceResults:
         fig : Plotly graph_objects.make_subplots()
             The figure object with the plot.
         """
-        if fig is None:
-            fig = make_subplots(
-                rows=1,
-                cols=2,
-                subplot_titles=[
-                    "<b>Frequency Evaluation</b>",
-                    "<b>Relative Error Evaluation</b>",
-                ],
-            )
+        kwargs_default_values = dict(plot_bgcolor="white", hoverlabel_align="right")
+
+        for k, v in kwargs_default_values.items():
+            kwargs.setdefault(k, v)
+
+        fig = make_subplots(
+            rows=1,
+            cols=2,
+            subplot_titles=[
+                "<b>Frequency Evaluation</b>",
+                "<b>Relative Error Evaluation</b>",
+            ],
+        )
 
         # plot Frequency vs number of elements
         fig.add_trace(
@@ -2640,6 +3014,8 @@ class ConvergenceResults:
                 x=self.el_num,
                 y=self.eigv_arr,
                 mode="lines+markers",
+                line=dict(width=3.0, color="FireBrick"),
+                marker=dict(size=10, color="FireBrick"),
                 hovertemplate=(
                     "Number of Elements: %{x:.2f}<br>" + "Frequency: %{y:.0f}"
                 ),
@@ -2648,8 +3024,30 @@ class ConvergenceResults:
             row=1,
             col=1,
         )
-        fig.update_xaxes(title_text="<b>Number of Elements</b>", row=1, col=1)
-        fig.update_yaxes(title_text="<b>Frequency</b>", row=1, col=1)
+        fig.update_xaxes(
+            title_text="<b>Number of Elements</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+            row=1,
+            col=1,
+        )
+        fig.update_yaxes(
+            title_text="<b>Frequency</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+            row=1,
+            col=1,
+        )
 
         # plot Error vs number of elements
         fig.add_trace(
@@ -2657,6 +3055,8 @@ class ConvergenceResults:
                 x=self.el_num,
                 y=self.error_arr,
                 mode="lines+markers",
+                line=dict(width=3.0, color="FireBrick"),
+                marker=dict(size=10, color="FireBrick"),
                 hovertemplate=(
                     "Number of Elements: %{x:.2f}<br>" + "Relative Error: %{y:.0f}"
                 ),
@@ -2666,8 +3066,31 @@ class ConvergenceResults:
             col=2,
         )
 
-        fig.update_xaxes(title_text="<b>Number of Elements</b>", row=1, col=2)
-        fig.update_yaxes(title_text="<b>Relative Error (%)</b>", row=1, col=2)
+        fig.update_xaxes(
+            title_text="<b>Number of Elements</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+            row=1,
+            col=2,
+        )
+        fig.update_yaxes(
+            title_text="<b>Relative Error (%)</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            exponentformat="power",
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+            row=1,
+            col=2,
+        )
 
         fig.update_layout(**kwargs)
 
@@ -2719,7 +3142,7 @@ class TimeResponseResults:
         self.nodes_pos = nodes_pos
         self.number_dof = number_dof
 
-    def _plot1d(self, dof, fig=None, **kwargs):
+    def _plot1d(self, dof, **kwargs):
         """Plot time response for a single DoF using Plotly.
 
         This function will take a rotor object and plot its time response using Plotly.
@@ -2728,8 +3151,6 @@ class TimeResponseResults:
         ----------
         dof : int
             Degree of freedom that will be observed.
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
@@ -2749,14 +3170,20 @@ class TimeResponseResults:
         obs_dof = dof % self.number_dof
         obs_dof = dof_dict[str(obs_dof)]
 
-        if fig is None:
-            fig = go.Figure()
+        kwargs_default_values = dict(
+            width=1200, height=900, plot_bgcolor="white", hoverlabel_align="right"
+        )
+        for k, v in kwargs_default_values.items():
+            kwargs.setdefault(k, v)
+
+        fig = go.Figure()
 
         fig.add_trace(
             go.Scatter(
                 x=self.t,
                 y=self.yout[:, dof],
                 mode="lines",
+                line=dict(width=3.0, color="royalblue"),
                 name="Phase",
                 legendgroup="Phase",
                 showlegend=False,
@@ -2764,18 +3191,38 @@ class TimeResponseResults:
             )
         )
 
-        fig.update_xaxes(title_text="<b>Time</b>")
-        fig.update_yaxes(title_text="<b>Amplitude</b>")
+        fig.update_xaxes(
+            title_text="<b>Time</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            range=[0, self.t[-1]],
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_yaxes(
+            title_text="<b>Amplitude</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
         fig.update_layout(
             title=dict(
-                text="<b>Response for node {} - DoF {}</b>".format(dof // 4, obs_dof)
+                text="<b>Response for node {} - DoF {}</b>".format(dof // 4, obs_dof),
+                font=dict(size=20),
             ),
             **kwargs,
         )
 
         return fig
 
-    def _plot2d(self, node, fig=None, **kwargs):
+    def _plot2d(self, node, **kwargs):
         """Plot orbit response (2D).
 
         This function will take a rotor object and plot its orbit response using Plotly.
@@ -2784,8 +3231,6 @@ class TimeResponseResults:
         ----------
         node: int, optional
             Selected node to plot orbit.
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
@@ -2796,14 +3241,20 @@ class TimeResponseResults:
         fig : Plotly graph_objects.Figure()
             The figure object with the plot.
         """
-        if fig is None:
-            fig = go.Figure()
+        kwargs_default_values = dict(
+            width=1200, height=900, plot_bgcolor="white", hoverlabel_align="right"
+        )
+        for k, v in kwargs_default_values.items():
+            kwargs.setdefault(k, v)
+
+        fig = go.Figure()
 
         fig.add_trace(
             go.Scatter(
                 x=self.yout[:, self.number_dof * node],
                 y=self.yout[:, self.number_dof * node + 1],
                 mode="lines",
+                line=dict(width=3.0, color="royalblue"),
                 name="Phase",
                 legendgroup="Phase",
                 showlegend=False,
@@ -2813,23 +3264,42 @@ class TimeResponseResults:
             )
         )
 
-        fig.update_xaxes(title_text="<b>Amplitude - X direction</b>")
-        fig.update_yaxes(title_text="<b>Amplitude - Y direction</b>")
+        fig.update_xaxes(
+            title_text="<b>Amplitude - X direction</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
+        fig.update_yaxes(
+            title_text="<b>Amplitude - Y direction</b>",
+            title_font=dict(family="Arial", size=20),
+            tickfont=dict(size=16),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+        )
         fig.update_layout(
-            title=dict(text="<b>Response for node {}</b>".format(node)), **kwargs
+            title=dict(
+                text="<b>Response for node {}</b>".format(node), font=dict(size=20)
+            ),
+            **kwargs,
         )
 
         return fig
 
-    def _plot3d(self, fig=None, **kwargs):
+    def _plot3d(self, **kwargs):
         """Plot orbit response (3D).
 
         This function will take a rotor object and plot its orbit response using Plotly.
 
         Parameters
         ----------
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. hoverlabel_align="center", ...).
@@ -2840,8 +3310,11 @@ class TimeResponseResults:
         fig : Plotly graph_objects.Figure()
             The figure object with the plot.
         """
-        if fig is None:
-            fig = go.Figure()
+        kwargs_default_values = dict(hoverlabel_align="right")
+        for k, v in kwargs_default_values.items():
+            kwargs.setdefault(k, v)
+
+        fig = go.Figure()
 
         for n in self.nodes_list:
             x_pos = np.ones(self.yout.shape[0]) * self.nodes_pos[n]
@@ -2851,7 +3324,7 @@ class TimeResponseResults:
                     y=self.yout[:, self.number_dof * n],
                     z=self.yout[:, self.number_dof * n + 1],
                     mode="lines",
-                    line=dict(color=tableau_colors["blue"]),
+                    line=dict(width=3.0, color="royalblue"),
                     name="Mean",
                     legendgroup="mean",
                     showlegend=False,
@@ -2873,23 +3346,47 @@ class TimeResponseResults:
                 y=line,
                 z=line,
                 mode="lines",
-                line=dict(color="black", dash="dashdot"),
+                line=dict(width=2.0, color="black", dash="dashdot"),
                 showlegend=False,
             )
         )
 
         fig.update_layout(
+            width=1200,
+            height=900,
             scene=dict(
-                xaxis=dict(title=dict(text="<b>Rotor Length</b>")),
-                yaxis=dict(title=dict(text="<b>Amplitude - X</b>")),
-                zaxis=dict(title=dict(text="<b>Amplitude - Y</b>")),
+                bgcolor="white",
+                xaxis=dict(
+                    title=dict(text="<b>Rotor Length</b>", font=dict(size=14)),
+                    tickfont=dict(size=12),
+                    nticks=5,
+                    backgroundcolor="lightgray",
+                    gridcolor="white",
+                    showspikes=False,
+                ),
+                yaxis=dict(
+                    title=dict(text="<b>Amplitude - X</b>", font=dict(size=14)),
+                    tickfont=dict(size=12),
+                    nticks=5,
+                    backgroundcolor="lightgray",
+                    gridcolor="white",
+                    showspikes=False,
+                ),
+                zaxis=dict(
+                    title=dict(text="<b>Amplitude - Y</b>", font=dict(size=14)),
+                    tickfont=dict(size=12),
+                    nticks=5,
+                    backgroundcolor="lightgray",
+                    gridcolor="white",
+                    showspikes=False,
+                ),
             ),
             **kwargs,
         )
 
         return fig
 
-    def plot(self, plot_type="3d", dof=None, node=None, fig=None, **kwargs):
+    def plot(self, plot_type="3d", dof=None, node=None, **kwargs):
         """Plot time response.
 
         The plot type options are:
@@ -2913,8 +3410,6 @@ class TimeResponseResults:
             Selected node to plot orbit.
             Fill this attribute only when selection plot_type = "2d".
             Default is None
-        fig : Plotly graph_objects.Figure()
-            The figure object with the plot.
         kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py
index 9916e0c2..52af9a2a 100644
--- a/ross/rotor_assembly.py
+++ b/ross/rotor_assembly.py
@@ -18,7 +18,6 @@ import scipy.linalg as la
 import scipy.signal as signal
 import scipy.sparse.linalg as las
 import toml
-from scipy.interpolate import UnivariateSpline
 from scipy.optimize import newton
 
 from ross.bearing_seal_element import (BallBearingElement, BearingElement,
@@ -27,11 +26,10 @@ from ross.bearing_seal_element import (BallBearingElement, BearingElement,
                                        RollerBearingElement, SealElement)
 from ross.disk_element import DiskElement, DiskElement6DoF
 from ross.materials import steel
-from ross.point_mass import PointMass
 from ross.results import (CampbellResults, ConvergenceResults,
-                          CriticalSpeedResults, ForcedResponseResults,
-                          FrequencyResponseResults, ModalResults,
-                          StaticResults, SummaryResults, TimeResponseResults)
+                          ForcedResponseResults, FrequencyResponseResults,
+                          ModalResults, StaticResults, SummaryResults,
+                          TimeResponseResults)
 from ross.shaft_element import ShaftElement, ShaftElement6DoF
 from ross.utils import convert
 
@@ -119,7 +117,11 @@ class Rotor(object):
         tag=None,
     ):
 
-        self.parameters = {"min_w": min_w, "max_w": max_w, "rated_w": rated_w}
+        self.parameters = {
+            "min_w": min_w,
+            "max_w": max_w,
+            "rated_w": rated_w,
+        }
         if tag is None:
             self.tag = "Rotor 0"
 
@@ -191,11 +193,6 @@ class Rotor(object):
             )
         ]
 
-        # check if tags are unique
-        tags_list = [el.tag for el in self.elements]
-        if len(tags_list) != len(set(tags_list)):
-            raise ValueError("Tags should be unique.")
-
         self.number_dof = self._check_number_dof()
 
         ####################################################
@@ -547,28 +544,17 @@ class Rotor(object):
     def run_modal(self, speed, num_modes=12, sparse=True):
         """Run modal analysis.
 
-        Method to calculate eigenvalues and eigvectors for a given rotor system.
-        Tthe natural frequencies and dampings ratios are calculated for a given
-        rotor speed. It means that for each speed input there's a different set of
-        eigenvalues and eigenvectors, hence, different natural frequencies and damping
-        ratios are returned.
+        Method to calculate eigenvalues and eigvectors for a given rotor system
 
         Parameters
         ----------
         speed : float
             Speed at which the eigenvalues and eigenvectors will be calculated.
-        num_modes : int, optional
-            The number of eigenvalues and eigenvectors to be calculated using ARPACK.
-            If sparse=True, it determines the number of eigenvalues and eigenvectors
-            to be calculated. It must be smaller than Rotor.ndof - 1. It is not
-            possible to compute all eigenvectors of a matrix with ARPACK.
-            If sparse=False, num_modes does not have any effect over the method.
+        num_modes : int
+            Number of modes to be calculated. This uses scipy.sparse.eigs method.
             Default is 12.
         sparse : bool, optional
-            If True, ARPACK is used to calculate a desired number (according to
-            num_modes) or eigenvalues and eigenvectors.
-            If False, scipy.linalg.eig() is used to calculate all the eigenvalues and
-            eigenvectors.
+            If sparse, eigenvalues will be calculated with arpack.
             Default is True.
 
         Returns
@@ -587,14 +573,14 @@ class Rotor(object):
         Example
         -------
         >>> rotor = rotor_example()
-        >>> modal = rotor.run_modal(speed=0, sparse=False)
+        >>> modal = rotor.run_modal(speed=0)
         >>> modal.wn[:2]
         array([91.79655318, 96.28899977])
         >>> modal.wd[:2]
         array([91.79655318, 96.28899977])
-        >>> fig = modal.plot_mode_3d(0)
+        >>> fig = modal.plot_mode3D(0)
         """
-        evalues, evectors = self._eigen(speed, num_modes=num_modes, sparse=sparse)
+        evalues, evectors = self._eigen(speed, num_modes=num_modes)
         wn_len = num_modes // 2
         wn = (np.absolute(evalues))[:wn_len]
         wd = (np.imag(evalues))[:wn_len]
@@ -620,81 +606,6 @@ class Rotor(object):
 
         return modal_results
 
-    def run_critical_speed(self, num_modes=12, sparse=True, rtol=0.005):
-        """Calculate the critical speeds and damping ratios for the rotor model.
-
-        This function runs an iterative method over "run_modal()" to minimize
-        (using scipy.optimize.newton) the error between the rotor speed and the rotor
-        critical speeds (rotor speed - critical speed).
-
-        Differently from run_modal(), this function doesn't take a speed input because
-        it iterates over the natural frequencies calculated in the last iteration.
-        The initial value is considered to be the undamped natural frequecies for
-        speed = 0 (no gyroscopic effect).
-
-        Once the error is within an acceptable range defined by "rtol", it returns the
-        approximated critical speed.
-
-        With the critical speeds calculated, the function uses the results to
-        calculate the log dec and damping ratios for each critical speed.
-
-        Parameters
-        ----------
-        num_modes : int, optional
-            The number of eigenvalues and eigenvectors to be calculated using ARPACK.
-            If sparse=True, it determines the number of eigenvalues and eigenvectors
-            to be calculated. It must be smaller than Rotor.ndof - 1. It is not
-            possible to compute all eigenvectors of a matrix with ARPACK.
-            If sparse=False, num_modes does not have any effect over the method.
-            Default is 12.
-        sparse : bool, optional
-            If True, ARPACK is used to calculate a desired number (according to
-            num_modes) or eigenvalues and eigenvectors.
-            If False, scipy.linalg.eig() is used to calculate all the eigenvalues and
-            eigenvectors.
-            Default is True.
-        rtol : float, optional
-            Tolerance (relative) for termination. Applied to scipy.optimize.newton.
-            Default is 0.005 (0.5%).
-
-        Returns
-        -------
-        CriticalSpeedResults : array
-            CriticalSpeedResults.wn : undamped critical speeds.
-            CriticalSpeedResults.wd : damped critical speeds.
-            CriticalSpeedResults.log_dec : log_dec for each critical speed.
-            CriticalSpeedResults.damping_ratio : damping ratio for each critical speed.
-
-        Examples
-        --------
-        >>> rotor = rotor_example()
-        >>> results = rotor.run_critical_speed(num_modes=8)
-        >>> np.round(results.wd)
-        array([ 92.,  96., 271., 300.])
-        >>> np.round(results.wn)
-        array([ 92.,  96., 271., 300.])
-        """
-        _wn = self.run_modal(0, num_modes, sparse).wn
-        wn = np.zeros_like(_wn)
-        wd = np.zeros_like(_wn)
-        log_dec = np.zeros_like(_wn)
-        damping_ratio = np.zeros_like(_wn)
-
-        for i in range(len(wn)):
-            wn_func = lambda s: (s - self.run_modal(s, num_modes, sparse).wn[i])
-            wn[i] = newton(func=wn_func, x0=_wn[i], rtol=rtol)
-
-        for i in range(len(wn)):
-            wd_func = lambda s: (s - self.run_modal(s, num_modes, sparse).wd[i])
-            wd[i] = newton(func=wd_func, x0=wn[i], rtol=rtol)
-
-        for i, s in enumerate(wd):
-            modal = self.run_modal(s, num_modes, sparse)
-            log_dec[i] = modal.log_dec[i]
-            damping_ratio[i] = modal.damping_ratio[i]
-
-        return CriticalSpeedResults(wn, wd, log_dec, damping_ratio)
-
     def convergence(self, n_eigval=0, err_max=1e-02):
         """Run convergence analysis.
 
@@ -1000,76 +911,6 @@ class Rotor(object):
                     break
         # fmt: on
 
-    def _clustering_points(self, num_modes=12, num_points=10, modes=None, rtol=0.005):
-        """Create an array with points clustered close to the natural frequencies.
-
-        This method generates an automatic array to run frequency response analyses.
-        The frequency points are calculated based on the damped natural frequencies and
-        their respective damping ratios. The greater the damping ratio, the more spread
-        the points are. If the damping ratio, for a given critical speed, is smaller
-        than 0.005, it is redefined to be 0.005 (for this method only).
-
-        Parameters
-        ----------
-        num_modes : int, optional
-            The number of eigenvalues and eigenvectors to be calculated using ARPACK.
-            It also defines the range for the output array, since the method generates
-            points only for the critical speed calculated by run_critical_speed().
-            Default is 12.
-        num_points : int, optional
-            The number of points generated for each critical speed.
-            The method set the same number of points for slightly less and slightly
-            higher than the natural circular frequency. It means there'll be num_points
-            greater and num_points smaller than a given critical speed.
-            num_points may be between 2 and 12. Anything above this range defaults
-            to 10 and anything below this range defaults to 4.
-            The default is 10.
-        modes : list, optional
-            Modes that will be used to calculate the frequency response.
-            The possibilities are limited by the num_modes argument.
-            (all modes will be used if a list is not given).
-        rtol : float, optional
-            Tolerance (relative) for termination. Applied to scipy.optimize.newton in
-            run_critical_speed() method.
-            Default is 0.005 (0.5%).
-
-        Returns
-        -------
-        speed_range : array
-            Range of frequencies (or speed).
-
-        Examples
-        --------
-        >>> rotor = rotor_example()
-        >>> speed_range = rotor._clustering_points(num_modes=12, num_points=5)
-        >>> speed_range.shape
-        (60,)
-        """
-        critical_speeds = self.run_critical_speed(num_modes=num_modes, rtol=rtol)
-        omega = critical_speeds.wd
-        damping = critical_speeds.damping_ratio
-        damping = np.array([d if d >= 0.005 else 0.005 for d in damping])
-
-        if num_points > 12:
-            num_points = 10
-        elif num_points < 2:
-            num_points = 4
-
-        if modes is not None:
-            omega = omega[modes]
-            damping = damping[modes]
-
-        a = np.zeros((len(omega), num_points))
-        for i in range(len(omega)):
-            for j in range(num_points):
-                b = 2 * (num_points - j + 1) / (num_points - 1)
-                a[i, j] = 1 + damping[i] ** b
-
-        omega = omega.reshape((len(omega), 1))
-        speed_range = np.sort(np.ravel(np.concatenate((omega / a, omega * a))))
-
-        return speed_range
-
     @staticmethod
     def _index(eigenvalues):
         """Generate indexes to sort eigenvalues and eigenvectors.
@@ -1255,15 +1096,15 @@ class Rotor(object):
         >>> speed = 100.0
         >>> H = rotor.transfer_matrix(speed=speed)
         """
-        lti = self._lti(speed=speed)
-        B = lti.B
-        C = lti.C
-        D = lti.D
+        modal = self.run_modal(speed=speed)
+        B = modal.lti.B
+        C = modal.lti.C
+        D = modal.lti.D
 
         # calculate eigenvalues and eigenvectors using la.eig to get
         # left and right eigenvectors.
 
-        evals, psi = self._eigen(speed=speed, frequency=frequency, sparse=False)
+        evals, psi, = la.eig(self.A(speed, frequency))
 
         psi_inv = la.inv(psi)
 
@@ -1274,6 +1115,7 @@ class Rotor(object):
             idx = np.zeros((2 * m), int)
             idx[0:m] = modes  # modes
             idx[m:] = range(2 * n)[-m:]  # conjugates (see how evalues are ordered)
+
             evals = evals[np.ix_(idx)]
             psi = psi[np.ix_(range(2 * n), idx)]
             psi_inv = psi_inv[np.ix_(idx, range(2 * n))]
@@ -1284,59 +1126,27 @@ class Rotor(object):
 
         return H
 
-    def run_freq_response(
-        self,
-        speed_range=None,
-        modes=None,
-        cluster_points=False,
-        num_modes=12,
-        num_points=10,
-        rtol=0.005,
-    ):
+    def run_freq_response(self, speed_range=None, modes=None):
         """Frequency response for a mdof system.
 
-        This method returns the frequency response for a mdof system given a range of
-        frequencies and the modes that will be used.
+        This method returns the frequency response for a mdof system
+        given a range of frequencies and the modes that will be used.
 
-        General parameters
-        ------------------
+        Parameters
+        ----------
         speed_range : array, optional
-            Array with the desired range of frequencies.
-            Default is 0 to 1.5 x highest damped natural frequency.
+            Array with the desired range of frequencies (the default
+             is 0 to 1.5 x highest damped natural frequency.
         modes : list, optional
             Modes that will be used to calculate the frequency response
             (all modes will be used if a list is not given).
 
-        Frequency spacing parameters
-        ----------------------------
-        cluster_points : bool, optional
-            boolean to activate the automatic frequency spacing method. If True, the
-            method uses _clustering_points() to create an speed_range.
-            Default is False
-        num_points : int, optional
-            The number of points generated per critical speed.
-            The method set the same number of points for slightly less and slightly
-            higher than the natural circular frequency. It means there'll be num_points
-            greater and num_points smaller than a given critical speed.
-            num_points may be between 2 and 12. Anything above this range defaults
-            to 10 and anything below this range defaults to 4.
-            The default is 10.
-        num_modes
-            The number of eigenvalues and eigenvectors to be calculated using ARPACK.
-            It also defines the range for the output array, since the method generates
-            points only for the critical speed calculated by run_critical_speed().
-            Default is 12.
-        rtol : float, optional
-            Tolerance (relative) for termination. Applied to scipy.optimize.newton to
-            calculate the approximated critical speeds.
-            Default is 0.005 (0.5%).
-
         Returns
         -------
         results : array
             Array with the frequencies, magnitude (dB) of the frequency
             response for each pair input/output, and
-            phase of the frequency response for each pair input/output.
+            phase of the frequency response for each pair input/output..
             It will be returned if plot=False.
 
         Examples
@@ -1347,25 +1157,12 @@ class Rotor(object):
         >>> response.magnitude # doctest: +ELLIPSIS
         array([[[1.00000000e-06, 1.00261725e-06, 1.01076952e-06, ...
 
-        Using clustered points option.
-        Set `cluster_points=True` and choose how many modes the method must search and
-        how many points to add just before and after each critical speed.
-
-        >>> response = rotor.run_freq_response(cluster_points=True, num_points=5)
-        >>> response.speed_range.shape
-        (60,)
-
         # plot frequency response function:
         >>> fig = response.plot(inp=13, out=13)
         """
         if speed_range is None:
-            if not cluster_points:
-                modal = self.run_modal(0)
-                speed_range = np.linspace(0, max(modal.evalues.imag) * 1.5, 1000)
-            else:
-                speed_range = self._clustering_points(
-                    num_modes, num_points, modes, rtol
-                )
+            modal = self.run_modal(0)
+            speed_range = np.linspace(0, max(modal.evalues.imag) * 1.5, 1000)
 
         self._check_frequency_array(speed_range)
 
@@ -1384,18 +1181,8 @@ class Rotor(object):
 
         return results
 
-    def forced_response(
-        self,
-        force=None,
-        speed_range=None,
-        modes=None,
-        cluster_points=False,
-        num_modes=12,
-        num_points=10,
-        rtol=0.005,
-        unbalance=None,
-    ):
-        """Forced response for a mdof system.
+    def forced_response(self, force=None, speed_range=None, modes=None, unbalance=None):
+        """Unbalanced response for a mdof system.
 
         This method returns the unbalanced response for a mdof system
         given magnitude and phase of the unbalance, the node where it's
@@ -1403,9 +1190,9 @@ class Rotor(object):
 
         Parameters
         ----------
-        force : list, array
+        force : list
             Unbalance force in each degree of freedom for each value in omega
-        speed_range : list, array
+        speed_range : list, float
             Array with the desired range of frequencies
         modes : list, optional
             Modes that will be used to calculate the frequency response
@@ -1416,30 +1203,6 @@ class Rotor(object):
             response analysis.
             Default is None.
 
-        Frequency spacing parameters
-        ----------------------------
-        cluster_points : bool, optional
-            boolean to activate the automatic frequency spacing method. If True, the
-            method uses _clustering_points() to create an speed_range.
-            Default is False
-        num_points : int, optional
-            The number of points generated per critical speed.
-            The method set the same number of points for slightly less and slightly
-            higher than the natural circular frequency. It means there'll be num_points
-            greater and num_points smaller than a given critical speed.
-            num_points may be between 2 and 12. Anything above this range defaults
-            to 10 and anything below this range defaults to 4.
-            The default is 10.
-        num_modes
-            The number of eigenvalues and eigenvectors to be calculated using ARPACK.
-            It also defines the range for the output array, since the method generates
-            points only for the critical speed calculated by run_critical_speed().
-            Default is 12.
-        rtol : float, optional
-            Tolerance (relative) for termination. Applied to scipy.optimize.newton to
-            calculate the approximated critical speeds.
-            Default is 0.005 (0.5%).
-
         Returns
         -------
         force_resp : array
@@ -1459,26 +1222,8 @@ class Rotor(object):
         >>> resp = rotor.forced_response(force=force, speed_range=speed)
         >>> resp.magnitude # doctest: +ELLIPSIS
         array([[0.00000000e+00, 5.06073311e-04, 2.10044826e-03, ...
-
-        Using clustered points option.
-        Set `cluster_points=True` and choose how many modes the method must search and
-        how many points to add just before and after each critical speed.
-
-        >>> response = rotor.forced_response(
-        ...     force=force, cluster_points=True, num_modes=12, num_points=5
-        ... )
-        >>> response.speed_range.shape
-        (60,)
         """
-        if speed_range is None:
-            if cluster_points:
-                speed_range = self._clustering_points(
-                    num_modes, num_points, modes, rtol
-                )
-
-        freq_resp = self.run_freq_response(
-            speed_range, modes, cluster_points, num_modes, num_points, rtol
-        )
+        freq_resp = self.run_freq_response(speed_range=speed_range, modes=modes)
 
         forced_resp = np.zeros(
             (self.ndof, len(freq_resp.speed_range)), dtype=np.complex
@@ -1546,16 +1291,7 @@ class Rotor(object):
         return F0
 
     def run_unbalance_response(
-        self,
-        node,
-        magnitude,
-        phase,
-        frequency_range=None,
-        modes=None,
-        cluster_points=False,
-        num_modes=12,
-        num_points=10,
-        rtol=0.005,
+        self, node, magnitude, phase, frequency_range=None, modes=None
     ):
         """Unbalanced response for a mdof system.
 
@@ -1577,30 +1313,6 @@ class Rotor(object):
             Modes that will be used to calculate the frequency response
             (all modes will be used if a list is not given).
 
-        Frequency spacing parameters
-        ----------------------------
-        cluster_points : bool, optional
-            boolean to activate the automatic frequency spacing method. If True, the
-            method uses _clustering_points() to create an speed_range.
-            Default is False
-        num_points : int, optional
-            The number of points generated per critical speed.
-            The method set the same number of points for slightly less and slightly
-            higher than the natural circular frequency. It means there'll be num_points
-            greater and num_points smaller than a given critical speed.
-            num_points may be between 2 and 12. Anything above this range defaults
-            to 10 and anything below this range defaults to 4.
-            The default is 10.
-        num_modes
-            The number of eigenvalues and eigenvectors to be calculated using ARPACK.
-            It also defines the range for the output array, since the method generates
-            points only for the critical speed calculated by run_critical_speed().
-            Default is 12.
-        rtol : float, optional
-            Tolerance (relative) for termination. Applied to scipy.optimize.newton to
-            calculate the approximated critical speeds.
-            Default is 0.005 (0.5%).
-
         Returns
         -------
         force_resp : array
@@ -1625,29 +1337,13 @@ class Rotor(object):
         >>> response.magnitude # doctest: +ELLIPSIS
         array([[0.00000000e+00, 5.06073311e-04, 2.10044826e-03, ...
 
-        Using clustered points option.
-        Set `cluster_points=True` and choose how many modes the method must search and
-        how many points to add just before and after each critical speed.
-
-        >>> response2 = rotor.run_unbalance_response(
-        ...     node=3, magnitude=0.01, phase=0.0, cluster_points=True, num_points=5
-        ... )
-        >>> response2.speed_range.shape
-        (60,)
-
-        plot unbalance response:
+        # plot unbalance response:
         >>> fig = response.plot(dof=13)
 
-        plot deflected shape configuration
+        # plot deflected shape configuration
         >>> value = 600
         >>> fig = response.plot_deflected_shape(speed=value)
         """
-        if frequency_range is None:
-            if cluster_points:
-                frequency_range = self._clustering_points(
-                    num_modes, num_points, modes, rtol
-                )
-
         force = np.zeros((self.ndof, len(frequency_range)), dtype=np.complex)
 
         try:
@@ -1656,12 +1352,8 @@ class Rotor(object):
         except TypeError:
             force = self._unbalance_force(node, magnitude, phase, frequency_range)
 
-        # fmt: off
         ub = np.vstack((node, magnitude, phase))
-        forced_response = self.forced_response(
-            force, frequency_range, modes, cluster_points, num_modes, num_points, rtol, ub
-        )
-        # fmt: on
+        forced_response = self.forced_response(force, frequency_range, modes, ub)
 
         return forced_response
 
@@ -1776,10 +1468,11 @@ class Rotor(object):
                 x=x_pos,
                 y=y_pos,
                 text=text,
+                textfont=dict(family="Verdana", size=20, color="black"),
                 mode="markers+text",
                 marker=dict(
                     opacity=0.7,
-                    size=20,
+                    size=30,
                     color="#ffcc99",
                     line=dict(width=1.0, color="black"),
                 ),
@@ -1817,17 +1510,31 @@ class Rotor(object):
 
         fig.update_xaxes(
             title_text="<b>Axial location</b>",
+            title_font=dict(family="Verdana", size=20),
+            tickfont=dict(size=16),
             range=[-0.1 * shaft_end, 1.1 * shaft_end],
             showgrid=False,
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
             mirror=True,
         )
         fig.update_yaxes(
             title_text="<b>Shaft radius</b>",
+            title_font=dict(family="Verdana", size=20),
+            tickfont=dict(size=16),
             range=[-0.3 * shaft_end, 0.3 * shaft_end],
             showgrid=False,
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
             mirror=True,
         )
-        fig.update_layout(title=dict(text="<b>Rotor Model</b>"), **kwargs)
+        fig.update_layout(
+            plot_bgcolor="white",
+            title=dict(text="<b>Rotor Model</b>", font=dict(size=18)),
+            **kwargs,
+        )
 
         return fig
 
@@ -1893,97 +1600,11 @@ class Rotor(object):
 
         return results
 
-    def _calc_ucs(
-        self,
-        stiffness_range=None,
-        num_modes=16,
-        num=20,
-        fig=None,
-        synchronous=False,
-        **kwargs,
-    ):
-
-        if stiffness_range is None:
-            if self.rated_w is not None:
-                bearing = self.bearing_elements[0]
-                k = bearing.kxx.interpolated(self.rated_w)
-                k = int(np.log10(k))
-                stiffness_range = (k - 3, k + 3)
-            else:
-                stiffness_range = (6, 11)
-
-        stiffness_log = np.logspace(*stiffness_range, num=num)
-        rotor_wn = np.zeros((self.number_dof, len(stiffness_log)))
-
-        bearings_elements = []  # exclude the seals
-        for bearing in self.bearing_elements:
-            if not isinstance(bearing, SealElement):
-                bearings_elements.append(bearing)
-
-        for i, k in enumerate(stiffness_log):
-            bearings = [BearingElement(b.n, kxx=k, cxx=0) for b in bearings_elements]
-            rotor = self.__class__(self.shaft_elements, self.disk_elements, bearings)
-            speed = 0
-            if synchronous:
-
-                def wn_diff(x):
-                    """Function to evaluate difference between speed and
-                    natural frequency for the first mode."""
-                    modal = rotor.run_modal(speed=x, num_modes=num_modes)
-                    # get first forward mode
-                    if modal.whirl_direction()[0] == "Forward":
-                        wn0 = modal.wn[0]
-                    else:
-                        wn0 = modal.wn[1]
-
-                    return wn0 - x
-
-                speed = newton(wn_diff, 0)
-            modal = rotor.run_modal(speed=speed, num_modes=num_modes)
-
-            # if sync, select only forward modes
-            if synchronous:
-                rotor_wn[:, i] = modal.wn[modal.whirl_direction() == "Forward"]
-            # if not sync, with speed=0 whirl direction can be confusing, with
-            # two close modes being forward or backward, so we select on mode in
-            # each 2 modes.
-            else:
-                rotor_wn[:, i] = modal.wn[
-                    : int(self.number_dof * 2) : int(self.number_dof / 2)
-                ]
-
-        bearing0 = bearings_elements[0]
-
-        # calculate interception points
-        intersection_points = {"x": [], "y": []}
-
-        # if bearing does not have constant coefficient, check intersection points
-        if not np.isnan(bearing0.frequency).all():
-            for j in range(rotor_wn.shape[0]):
-                for coeff in ["kxx", "kyy"]:
-                    with warnings.catch_warnings():
-                        warnings.simplefilter("ignore")
-                        frequency_curve = UnivariateSpline(rotor_wn[j], stiffness_log)
-                    coeff_curve = getattr(bearing0, coeff).interpolated
-                    # search diff in several points to have a good precision
-                    search_range = np.linspace(
-                        bearing0.frequency[0], bearing0.frequency[-1], 10000
-                    )
-                    intersection_idx = np.argmin(
-                        abs(frequency_curve(search_range) - coeff_curve(search_range))
-                    )
-                    intersection_points["y"].append(search_range[intersection_idx])
-                    intersection_points["x"].append(
-                        frequency_curve(search_range[intersection_idx])
-                    )
-
-        return stiffness_log, rotor_wn, bearing0, intersection_points
-
     def plot_ucs(
         self,
         stiffness_range=None,
         num_modes=16,
-        num=30,
+        num=20,
         fig=None,
         synchronous=False,
         **kwargs,
@@ -2046,49 +1667,66 @@ class Rotor(object):
         >>> rotor = Rotor(shaft_elem, [disk0, disk1], [bearing0, bearing1])
         >>> fig = rotor.plot_ucs()
         """
+        if stiffness_range is None:
+            if self.rated_w is not None:
+                bearing = self.bearing_elements[0]
+                k = bearing.kxx.interpolated(self.rated_w)
+                k = int(np.log10(k))
+                stiffness_range = (k - 3, k + 3)
+            else:
+                stiffness_range = (6, 11)
 
-        stiffness_log, rotor_wn, bearing0, intersection_points = self._calc_ucs(
-            stiffness_range=stiffness_range,
-            num_modes=num_modes,
-            num=num,
-            fig=fig,
-            synchronous=synchronous,
-            **kwargs,
-        )
+        stiffness_log = np.logspace(*stiffness_range, num=num)
+        rotor_wn = np.zeros((self.number_dof, len(stiffness_log)))
 
-        if fig is None:
-            fig = go.Figure()
+        bearings_elements = []  # exclude the seals
+        for bearing in self.bearing_elements:
+            if not isinstance(bearing, SealElement):
+                bearings_elements.append(bearing)
 
-        for j in range(rotor_wn.shape[0]):
-            fig.add_trace(
-                go.Scatter(
-                    x=stiffness_log,
-                    y=rotor_wn[j],
-                    mode="lines",
-                    hoverinfo="none",
-                    showlegend=False,
-                )
-            )
+        for i, k in enumerate(stiffness_log):
+            bearings = [BearingElement(b.n, kxx=k, cxx=0) for b in bearings_elements]
+            rotor = self.__class__(self.shaft_elements, self.disk_elements, bearings)
+            speed = 0
+            if synchronous:
 
-        fig.add_trace(
-            go.Scatter(
-                x=intersection_points["x"],
-                y=intersection_points["y"],
-                mode="markers",
-                marker=dict(symbol="circle-open-dot", color="red", size=8),
-                hovertemplate="Stiffness (N/m): %{x:.2e}<br>Frequency (rad/s): %{y:.2f}",
-                showlegend=False,
-                name="",
-            )
-        )
+                def wn_diff(x):
+                    """Function to evaluate difference between speed and
+                    natural frequency for the first mode."""
+                    modal = rotor.run_modal(speed=x, num_modes=num_modes)
+                    # get first forward mode
+                    if modal.whirl_direction()[0] == "Forward":
+                        wn0 = modal.wn[0]
+                    else:
+                        wn0 = modal.wn[1]
+
+                    return wn0 - x
+
+                speed = newton(wn_diff, 0)
+            modal = rotor.run_modal(speed=speed, num_modes=num_modes)
+
+            # if sync, select only forward modes
+            if synchronous:
+                rotor_wn[:, i] = modal.wn[modal.whirl_direction() == "Forward"]
+            # if not sync, with speed=0 whirl direction can be confusing, with
+            # two close modes being forward or backward, so we select on mode in
+            # each 2 modes.
+            else:
+                rotor_wn[:, i] = modal.wn[
+                    : int(self.number_dof * 2) : int(self.number_dof / 2)
+                ]
+
+        bearing0 = bearings_elements[0]
+
+        if fig is None:
+            fig = go.Figure()
 
         fig.add_trace(
             go.Scatter(
                 x=bearing0.kxx.interpolated(bearing0.frequency),
                 y=bearing0.frequency,
-                mode="lines",
-                line=dict(dash="dashdot"),
-                hoverinfo="none",
+                mode="markers",
+                marker=dict(size=10, symbol="circle", color="#888844"),
                 name="Kxx",
             )
         )
@@ -2096,25 +1734,57 @@ class Rotor(object):
             go.Scatter(
                 x=bearing0.kyy.interpolated(bearing0.frequency),
                 y=bearing0.frequency,
-                mode="lines",
-                line=dict(dash="dashdot"),
-                hoverinfo="none",
+                mode="markers",
+                marker=dict(size=10, symbol="square", color="#888844"),
                 name="Kyy",
             )
         )
 
+        for j in range(rotor_wn.T.shape[1]):
+            fig.add_trace(
+                go.Scatter(
+                    x=stiffness_log,
+                    y=np.transpose(rotor_wn.T)[j],
+                    mode="lines",
+                    line=dict(width=3, color=colors[j]),
+                    hoverinfo="none",
+                    showlegend=False,
+                )
+            )
         fig.update_xaxes(
-            title_text="<b>Bearing Stiffness (N/m)</b>",
+            title_text="<b>Bearing Stiffness</b>",
+            title_font=dict(size=16),
+            tickfont=dict(size=14),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
             type="log",
             exponentformat="power",
         )
         fig.update_yaxes(
-            title_text="<b>Critical Speed (rad/s)</b>",
+            title_text="<b>Critical Speed</b>",
+            title_font=dict(size=16),
+            tickfont=dict(size=14),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
             type="log",
             exponentformat="power",
         )
         fig.update_layout(
-            title=dict(text="<b>Undamped Critical Speed Map</b>"), **kwargs
+            plot_bgcolor="white",
+            legend=dict(
+                font=dict(family="sans-serif", size=14),
+                bgcolor="white",
+                bordercolor="black",
+                borderwidth=2,
+            ),
+            title=dict(text="<b>Undamped Critical Speed Map</b>", font=dict(size=16)),
+            **kwargs,
         )
 
         return fig
@@ -2132,8 +1802,6 @@ class Rotor(object):
             Default is 5.
         stiffness_range : tuple, optional
             Tuple with (start, end) for stiffness range.
-            This will be used to create an evenly numbers spaced evenly on a log scale
-            to create a better visualization (see np.logspace).
         kwargs : optional
             Additional key word arguments can be passed to change the plot layout only
             (e.g. width=1000, height=800, ...).
@@ -2210,11 +1878,39 @@ class Rotor(object):
         )
 
         fig.update_xaxes(
-            title_text="<b>Applied Cross Coupled Stiffness</b>", exponentformat="power"
+            title_text="<b>Applied Cross Coupled Stiffness</b>",
+            title_font=dict(size=16),
+            tickfont=dict(size=14),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+            exponentformat="power",
+        )
+        fig.update_yaxes(
+            title_text="<b>Log Dec</b>",
+            title_font=dict(size=16),
+            tickfont=dict(size=14),
+            gridcolor="lightgray",
+            showline=True,
+            linewidth=2.5,
+            linecolor="black",
+            mirror=True,
+            exponentformat="power",
         )
-        fig.update_yaxes(title_text="<b>Log Dec</b>", exponentformat="power")
         fig.update_layout(
-            title=dict(text="<b>Level 1 stability analysis</b>"), **kwargs
+            width=1200,
+            height=900,
+            plot_bgcolor="white",
+            legend=dict(
+                font=dict(family="sans-serif", size=14),
+                bgcolor="white",
+                bordercolor="black",
+                borderwidth=2,
+            ),
+            title=dict(text="<b>Level 1 stability analysis</b>", font=dict(size=16)),
+            **kwargs,
         )
 
         return fig
@@ -2254,7 +1950,7 @@ class Rotor(object):
         >>> F[:, 4 * node] = 10 * np.cos(2 * t)
         >>> F[:, 4 * node + 1] = 10 * np.sin(2 * t)
         >>> response = rotor.run_time_response(speed, F, t)
-        >>> response.yout[:, dof] # doctest: +ELLIPSIS
+        >>> response.yout[:, dof]
         array([ 0.00000000e+00,  1.86686693e-07,  8.39130663e-07, ...
 
         # plot time response for a single DoF:
@@ -2274,12 +1970,12 @@ class Rotor(object):
 
         return results
 
-    def save_mat(self, file, speed, frequency=None):
+    def save_mat(self, file_path, speed, frequency=None):
         """Save matrices and rotor model to a .mat file.
 
         Parameters
         ----------
-        file : str, pathlib.Path
+        file_path : str
 
         speed: float
             Rotor speed.
@@ -2290,7 +1986,7 @@ class Rotor(object):
         Examples
         --------
         >>> rotor = rotor_example()
-        >>> rotor.save_mat('/tmp/new_matrices', speed=0)
+        >>> rotor.save_mat('new_matrices.mat', speed=0)
         """
         if frequency is None:
             frequency = speed
@@ -2303,33 +1999,53 @@ class Rotor(object):
             "nodes": self.nodes_pos,
         }
 
-        sio.savemat(file, dic)
+        sio.savemat("%s/%s.mat" % (os.getcwd(), file_path), dic)
 
-    def save(self, file):
-        """Save the rotor to a .toml file.
+    def save(self, rotor_name="rotor", file_path=Path(".")):
+        """Save rotor to toml file.
 
         Parameters
         ----------
-        file : str or pathlib.Path
+        file_path : str
 
         Examples
         --------
         >>> rotor = rotor_example()
-        >>> rotor.save('/tmp/rotor.toml')
+        >>> rotor.save('new_rotor')
+        >>> Rotor.remove('new_rotor')
         """
-        with open(file, "w") as f:
+        path_rotor = Path(file_path)
+
+        if os.path.isdir(path_rotor / rotor_name):
+            if int(
+                input(
+                    "There is a rotor with this file_path, do you want to overwrite it? (1 for yes and 0 for no)"
+                )
+            ):
+                shutil.rmtree(path_rotor / rotor_name)
+            else:
+                return "The rotor was not saved."
+
+        os.mkdir(path_rotor / rotor_name)
+        rotor_folder = path_rotor / rotor_name
+        os.mkdir(rotor_folder / "results")
+        os.mkdir(rotor_folder / "elements")
+
+        with open(rotor_folder / "properties.toml", "w") as f:
             toml.dump({"parameters": self.parameters}, f)
-        for el in self.elements:
-            el.save(file)
 
-    @classmethod
-    def load(cls, file):
+        elements_folder = rotor_folder / "elements"
+
+        for element in self.elements:
+            element.save(elements_folder)
+
+    @staticmethod
+    def load(file_path):
         """Load rotor from toml file.
 
         Parameters
         ----------
-        file : str or pathlib.Path
-            String or Path for a .toml file.
+        file_path : str
 
         Returns
         -------
@@ -2338,42 +2054,58 @@ class Rotor(object):
         Example
         -------
         >>> rotor1 = rotor_example()
-        >>> rotor1.save('/tmp/new_rotor1.toml')
-        >>> rotor2 = Rotor.load('/tmp/new_rotor1.toml')
+        >>> rotor1.save(Path('.')/'new_rotor1')
+        >>> rotor2 = Rotor.load(Path('.')/'new_rotor1')
         >>> rotor1 == rotor2
         True
+        >>> Rotor.remove('new_rotor1')
         """
-        data = toml.load(file)
-        parameters = data["parameters"]
+        rotor_path = Path(file_path)
 
-        elements = []
-        for el_name, el_data in data.items():
-            if el_name == "parameters":
-                continue
-            class_name = el_name.split("_")[0]
-            elements.append(globals()[class_name].read_toml_data(el_data))
-
-        shaft_elements = []
-        disk_elements = []
-        bearing_elements = []
-        point_mass_elements = []
-        for el in elements:
-            if isinstance(el, ShaftElement):
-                shaft_elements.append(el)
-            elif isinstance(el, DiskElement):
-                disk_elements.append(el)
-            elif isinstance(el, BearingElement):
-                bearing_elements.append(el)
-            elif isinstance(el, PointMass):
-                point_mass_elements.append(el)
+        if os.path.isdir(rotor_path / "elements"):
+            elements_path = rotor_path / "elements"
+        else:
+            raise FileNotFoundError("Elements folder not found.")
+
+        with open(rotor_path / "properties.toml", "r") as f:
+            parameters = toml.load(f)["parameters"]
+
+        global_elements = {}
+        for el in os.listdir(elements_path):
+            elements = []
+            if ".toml" in el:
+                with open(Path(elements_path) / el, "r") as f:
+                    el_dict = toml.load(f)
+                    element_class = list(el_dict.keys())[0]
+                    for el_number in el_dict[element_class]:
+                        element = (
+                            element_class + f"(**{el_dict[element_class][el_number]})"
+                        )
+                        elements.append(eval(element))
+            global_elements[convert(element_class + "s")] = elements
 
-        return cls(
-            shaft_elements=shaft_elements,
-            disk_elements=disk_elements,
-            bearing_elements=bearing_elements,
-            point_mass_elements=point_mass_elements,
-            **parameters,
-        )
+        return Rotor(**global_elements, **parameters)
+
+    @staticmethod
+    def remove(file_path):
+        """
+        Remove a previously saved rotor in rotors folder.
+
+        Parameters
+        ----------
+        file_path : str
+
+        Example
+        -------
+        >>> rotor = rotor_example()
+        >>> rotor.save('new_rotor2')
+        >>> Rotor.remove('new_rotor2')
+        """
+        try:
+            Rotor.load(file_path)
+            shutil.rmtree(Path(file_path))
+        except:
+            return "This is not a valid rotor."
 
     def run_static(self):
         """Run static analysis.
@@ -2937,7 +2669,11 @@ class CoAxialRotor(Rotor):
         tag=None,
     ):
 
-        self.parameters = {"min_w": min_w, "max_w": max_w, "rated_w": rated_w}
+        self.parameters = {
+            "min_w": min_w,
+            "max_w": max_w,
+            "rated_w": rated_w,
+        }
         if tag is None:
             self.tag = "Rotor 0"
 
diff --git a/ross/shaft_element.py b/ross/shaft_element.py
index 8ba9f7c4..b07119d4 100644
--- a/ross/shaft_element.py
+++ b/ross/shaft_element.py
@@ -3,7 +3,6 @@
 This module defines the ShaftElement classes which will be used to represent the rotor
 shaft. There're 2 options, an element with 8 or 12 degrees of freedom.
 """
-import inspect
 import os
 from pathlib import Path
 
@@ -381,34 +380,70 @@ class ShaftElement(Element):
     def __hash__(self):
         return hash(self.tag)
 
-    def save(self, file):
-        signature = inspect.signature(self.__init__)
-        args_list = list(signature.parameters)
-        args = {arg: getattr(self, arg) for arg in args_list}
-
-        # add material characteristics so that the shaft element can be reconstructed
-        # even if the material is not in the available_materials file.
-        args["material"] = {
-            "name": self.material.name,
-            "rho": self.material.rho,
-            "E": self.material.E,
-            "G_s": self.material.G_s,
-            "color": self.material.color,
+    def save(self, file_name=Path(os.getcwd())):
+        """Save shaft elements to toml file.
+
+        Parameters
+        ----------
+        file_name : str
+
+        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.save()
+        """
+        data = self.get_data(Path(file_name) / "ShaftElement.toml")
+        data["ShaftElement"][str(self.n)] = {
+            "L": self.L,
+            "idl": self.idl,
+            "odl": self.odl,
+            "idr": self.idr,
+            "odr": self.odr,
+            "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,
+            "shear_method_calc": self.shear_method_calc,
         }
+        self.dump_data(data, Path(file_name) / "ShaftElement.toml")
 
-        try:
-            data = toml.load(file)
-        except FileNotFoundError:
-            data = {}
+    @staticmethod
+    def load(file_name="ShaftElement"):
+        """Load a list of shaft elements saved in a toml format.
 
-        data[f"{self.__class__.__name__}_{self.tag}"] = args
-        with open(file, "w") as f:
-            toml.dump(data, f)
+        It works as an auxiliary function of the load function in the Rotor class.
 
-    @classmethod
-    def read_toml_data(cls, data):
-        data["material"] = Material(**data["material"])
-        return cls(**data)
+        Parameters
+        ----------
+        file_name : str, optional
+
+        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.save(os.getcwd())
+        >>> shaft2 = ShaftElement.load(os.getcwd())
+        >>> shaft2 # doctest: +ELLIPSIS
+        [ShaftElement(L=0.25, idl=0.0...
+        """
+        shaft_elements = []
+        with open("ShaftElement.toml", "r") as f:
+            shaft_elements_dict = toml.load(f)
+            for element in shaft_elements_dict["ShaftElement"]:
+                shaft_elements.append(
+                    ShaftElement(**shaft_elements_dict["ShaftElement"][element])
+                )
+        return shaft_elements
 
     @property
     def n(self):
@@ -1453,34 +1488,76 @@ class ShaftElement6DoF(ShaftElement):
             f"\n"
         )
 
-    def save(self, file):
-        signature = inspect.signature(self.__init__)
-        args_list = list(signature.parameters)
-        args = {arg: getattr(self, arg) for arg in args_list}
-
-        # add material characteristics so that the shaft element can be reconstructed
-        # even if the material is not in the available_materials file.
-        args["material"] = {
-            "name": self.material.name,
-            "rho": self.material.rho,
-            "E": self.material.E,
-            "G_s": self.material.G_s,
-            "color": self.material.color,
+    def save(self, file_name=Path(os.getcwd())):
+        """Save shaft elements to toml file.
+
+        Parameters
+        ----------
+        file_name : str
+
+        Examples
+        --------
+        >>> from ross.materials import steel
+        >>> shaft1 = ShaftElement6DoF(
+        ...        L=0.25, idl=0, idr=0, odl=0.05, odr=0.08,
+        ...        material=steel, rotary_inertia=True, shear_effects=True,
+        ... )
+        >>> shaft1.save()
+        """
+        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")
 
-        try:
-            data = toml.load(file)
-        except FileNotFoundError:
-            data = {}
+    @staticmethod
+    def load(file_name="ShaftElement6DoF"):
+        """Load a list of shaft elements saved in a toml format.
 
-        data[f"{self.__class__.__name__}_{self.tag}"] = args
-        with open(file, "w") as f:
-            toml.dump(data, f)
+        It works as an auxiliary function of the load function in the Rotor class.
 
-    @classmethod
-    def read_toml_data(cls, data):
-        data["material"] = Material(**data["material"])
-        return cls(**data)
+        Parameters
+        ----------
+        file_name : str, optional
+
+        Returns
+        -------
+        shaft_elements : list
+            A list of shaft elements
+
+        Examples
+        --------
+        >>> from ross.materials import steel
+        >>> shaft1 = ShaftElement6DoF(
+        ...        L=0.25, idl=0, idr=0, odl=0.05, odr=0.08,
+        ...        material=steel, rotary_inertia=True, shear_effects=True
+        ... )
+        >>> shaft1.save(os.getcwd())
+        >>> shaft2 = ShaftElement6DoF.load(os.getcwd())
+        >>> shaft2 # doctest: +ELLIPSIS
+        [ShaftElement6DoF(L=0.25, idl=0.0...
+        """
+        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
 
     def dof_mapping(self):
         """Degrees of freedom mapping.
diff --git a/ross/stochastic/st_bearing_seal_element.py b/ross/stochastic/st_bearing_seal_element.py
index 88f50286..32db0a05 100644
--- a/ross/stochastic/st_bearing_seal_element.py
+++ b/ross/stochastic/st_bearing_seal_element.py
@@ -8,7 +8,7 @@ import numpy as np
 from ross.bearing_seal_element import BearingElement
 from ross.fluid_flow import fluid_flow as flow
 from ross.fluid_flow.fluid_flow_coefficients import (
-    calculate_short_damping_matrix, calculate_short_stiffness_matrix)
+    calculate_damping_matrix, calculate_stiffness_matrix)
 from ross.stochastic.st_results_elements import plot_histogram
 
 # fmt: on
@@ -512,8 +512,8 @@ class ST_BearingElement:
                 eccentricity=attribute_dict["eccentricity"][i],
                 load=attribute_dict["load"][i],
             )
-            c = calculate_short_damping_matrix(fluid_flow)
-            k = calculate_short_stiffness_matrix(fluid_flow)
+            c = calculate_damping_matrix(fluid_flow, force_type="short")
+            k = calculate_stiffness_matrix(fluid_flow, force_type="short")
             args_dict["kxx"].append(k[0])
             args_dict["kxy"].append(k[1])
             args_dict["kyx"].append(k[2])
diff --git a/ross/utils.py b/ross/utils.py
index bdb986a3..a32c34e3 100644
--- a/ross/utils.py
+++ b/ross/utils.py
@@ -393,13 +393,32 @@ def visualize_matrix(rotor, matrix, frequency=None, **kwargs):
         )
     )
 
+    fig.update_xaxes(
+        tickfont=dict(size=16),
+        showline=True,
+        linewidth=2.5,
+        linecolor="black",
+        mirror=True,
+    )
+    fig.update_yaxes(
+        tickfont=dict(size=16),
+        showline=True,
+        linewidth=2.5,
+        linecolor="black",
+        mirror=True,
+    )
     fig.update_layout(
+        width=1200,
+        height=900,
+        hoverlabel_align="left",
         coloraxis=dict(
             cmin=np.amin(A),
             cmax=np.amax(A),
             colorscale="Rainbow",
             colorbar=dict(
-                title=dict(text="<b>Value</b>", side="top"), exponentformat="power"
+                title=dict(text="<b>Value</b>", side="top", font=dict(size=20)),
+                tickfont=dict(size=16),
+                exponentformat="power",
             ),
         ),
         **kwargs,


Test patch:
diff --git a/ross/tests/test_bearing_seal_element.py b/ross/tests/test_bearing_seal_element.py
index 55ad9fc9..f3e684fe 100644
--- a/ross/tests/test_bearing_seal_element.py
+++ b/ross/tests/test_bearing_seal_element.py
@@ -290,25 +290,19 @@ def test_roller_bearing_element():
     assert_allclose(rollerbearing.G(), G)
 
 
-@pytest.fixture
-def bearing_6dof():
+def test_bearing_6dof():
     bearing_6dof = BearingElement6DoF(
         n=0, kxx=1e6, kyy=0.8e6, kzz=1e5, cxx=2e2, cyy=1.5e2, czz=0.5e2
     )
-
-    return bearing_6dof
-
-
-def test_bearing6(bearing_6dof):
     # fmt: off
     K = np.array(
-        [[1000000., 0., 0.],
-         [0., 800000., 0.],
-         [0., 0., 100000.]])
+        [[1000000.,       0.,       0.],
+         [      0.,  800000.,       0.],
+         [      0.,       0.,  100000.]])
     C = np.array(
-        [[200., 0., 0.],
-         [0., 150., 0.],
-         [0., 0., 50.]])
+        [[200.,   0.,   0.],
+         [  0., 150.,   0.],
+         [  0.,   0.,  50.]])
     M = np.zeros((3, 3))
     G = np.zeros((3, 3))
     # fmt: on
@@ -333,17 +327,3 @@ def test_bearing_6dof_equality():
     assert bearing_6dof_0 == bearing_6dof_1
     assert not bearing_6dof_1 == bearing_6dof_2
     assert not bearing_6dof_0 == bearing_6dof_2
-
-
-def test_save_load(bearing0, bearing_constant, bearing_6dof):
-    bearing0.save("/tmp/bearing0.toml")
-    bearing0_loaded = BearingElement.load("/tmp/bearing0.toml")
-    assert bearing0 == bearing0_loaded
-
-    bearing_constant.save("/tmp/bearing_constant.toml")
-    bearing_constant_loaded = BearingElement.load("/tmp/bearing_constant.toml")
-    assert bearing_constant == bearing_constant_loaded
-
-    bearing_6dof.save("/tmp/bearing_6dof.toml")
-    bearing_6dof_loaded = BearingElement6DoF.load("/tmp/bearing_6dof.toml")
-    assert bearing_6dof == bearing_6dof_loaded
diff --git a/ross/tests/test_disk_element.py b/ross/tests/test_disk_element.py
index 01794e24..f8d0b2fc 100644
--- a/ross/tests/test_disk_element.py
+++ b/ross/tests/test_disk_element.py
@@ -47,12 +47,6 @@ def disk_from_geometry():
     return DiskElement.from_geometry(0, steel, 0.07, 0.05, 0.28)
 
 
-def test_save_load(disk, disk_from_geometry):
-    disk.save("/tmp/disk.toml")
-    disk_loaded = DiskElement.load("/tmp/disk.toml")
-    assert disk == disk_loaded
-
-
 def test_mass_matrix_disk1(disk_from_geometry):
     # fmt: off
     Md1 = np.array([[ 32.58973,   0.     ,   0.     ,   0.     ],
diff --git a/ross/tests/test_fluid_flow.py b/ross/tests/test_fluid_flow.py
index 5ea73416..f2408ad6 100644
--- a/ross/tests/test_fluid_flow.py
+++ b/ross/tests/test_fluid_flow.py
@@ -8,9 +8,8 @@ from numpy.testing import assert_allclose
 
 from ross.fluid_flow import fluid_flow as flow
 from ross.fluid_flow.fluid_flow_coefficients import (
-    calculate_oil_film_force, calculate_short_damping_matrix,
-    calculate_short_stiffness_matrix,
-    calculate_stiffness_and_damping_coefficients, find_equilibrium_position)
+    calculate_damping_matrix, calculate_oil_film_force,
+    calculate_stiffness_matrix, find_equilibrium_position)
 from ross.fluid_flow.fluid_flow_geometry import move_rotor_center
 from ross.fluid_flow.fluid_flow_graphics import (
     plot_eccentricity, plot_pressure_surface, plot_pressure_theta,
@@ -125,7 +124,7 @@ def test_stiffness_matrix():
     Taken from example 5.5.1, page 181 (Dynamics of rotating machine, FRISSWELL)
     """
     bearing = fluid_flow_short_friswell()
-    kxx, kxy, kyx, kyy = calculate_short_stiffness_matrix(bearing)
+    kxx, kxy, kyx, kyy = calculate_stiffness_matrix(bearing, force_type="short")
     assert math.isclose(kxx / 10 ** 6, 12.81, rel_tol=0.01)
     assert math.isclose(kxy / 10 ** 6, 16.39, rel_tol=0.01)
     assert math.isclose(kyx / 10 ** 6, -25.06, rel_tol=0.01)
@@ -141,31 +140,12 @@ def test_stiffness_matrix_numerical(fluid_flow_short_eccentricity):
     """
     bearing = fluid_flow_short_eccentricity
     bearing.calculate_pressure_matrix_numerical()
-    K, C = calculate_stiffness_and_damping_coefficients(bearing)
-    kxx, kxy, kyx, kyy = K[0], K[1], K[2], K[3]
-    k_xx, k_xy, k_yx, k_yy = calculate_short_stiffness_matrix(bearing)
-    assert_allclose(kxx, k_xx, rtol=0.29)
-    assert_allclose(kxy, k_xy, rtol=0.22)
-    assert_allclose(kyx, k_yx, rtol=0.15)
-    assert_allclose(kyy, k_yy, rtol=0.22)
-
-
-def test_damping_matrix_numerical(fluid_flow_short_eccentricity):
-    """
-    This function instantiate a bearing using the fluid flow class and test if it matches the
-    expected results for the damping matrix based on Friswell's book formulas, given the
-    eccentricity ratio.
-    Taken from chapter 5, page 179 (Dynamics of rotating machine, FRISSWELL)
-    """
-    bearing = fluid_flow_short_eccentricity
-    bearing.calculate_pressure_matrix_numerical()
-    K, C = calculate_stiffness_and_damping_coefficients(bearing)
-    cxx, cxy, cyx, cyy = C[0], C[1], C[2], C[3]
-    c_xx, c_xy, c_yx, c_yy = calculate_short_damping_matrix(bearing)
-    assert_allclose(cxx, c_xx, rtol=0.12)
-    assert_allclose(cxy, c_xy, rtol=0.10)
-    assert_allclose(cyx, c_yx, rtol=0.32)
-    assert_allclose(cyy, c_yy, rtol=0.02)
+    kxx, kxy, kyx, kyy = calculate_stiffness_matrix(bearing, force_type="numerical")
+    k_xx, k_xy, k_yx, k_yy = calculate_stiffness_matrix(bearing, force_type="short")
+    assert_allclose(kxx, k_xx, rtol=0.33)
+    assert_allclose(kxy, k_xy, rtol=0.21)
+    assert_allclose(kyx, k_yx, rtol=0.17)
+    assert_allclose(kyy, k_yy, rtol=0.19)
 
 
 def test_damping_matrix():
@@ -175,7 +155,7 @@ def test_damping_matrix():
     Taken from example 5.5.1, page 181 (Dynamics of rotating machine, FRISSWELL)
     """
     bearing = fluid_flow_short_friswell()
-    cxx, cxy, cyx, cyy = calculate_short_damping_matrix(bearing)
+    cxx, cxy, cyx, cyy = calculate_damping_matrix(bearing, force_type="short")
     assert math.isclose(cxx / 10 ** 3, 232.9, rel_tol=0.01)
     assert math.isclose(cxy / 10 ** 3, -81.92, rel_tol=0.01)
     assert math.isclose(cyx / 10 ** 3, -81.92, rel_tol=0.01)
diff --git a/ross/tests/test_point_mass.py b/ross/tests/test_point_mass.py
index 60f3b52f..a5d2dbb0 100644
--- a/ross/tests/test_point_mass.py
+++ b/ross/tests/test_point_mass.py
@@ -27,11 +27,3 @@ def test_local_index():
 
     assert p.dof_local_index().x_0 == 0
     assert p.dof_local_index().y_0 == 1
-
-
-def test_save_load():
-    p = PointMass(n=0, m=10.0, tag="pointmass")
-    p.save("/tmp/point_mass.toml")
-    p_loaded = p.load("/tmp/point_mass.toml")
-
-    assert p == p_loaded
diff --git a/ross/tests/test_rotor_assembly.py b/ross/tests/test_rotor_assembly.py
index c55afa6c..3ba5c66e 100644
--- a/ross/tests/test_rotor_assembly.py
+++ b/ross/tests/test_rotor_assembly.py
@@ -1,4 +1,4 @@
-from pathlib import Path
+import os
 
 import numpy as np
 import pytest
@@ -12,6 +12,8 @@ from ross.rotor_assembly import *
 from ross.rotor_assembly import MAC_modes
 from ross.shaft_element import *
 
+test_dir = os.path.dirname(__file__)
+
 
 @pytest.fixture
 def rotor1():
@@ -1050,67 +1052,6 @@ def test_static_analysis_rotor6(rotor6):
     )
 
 
-def test_run_critical_speed(rotor5, rotor6):
-    results5 = rotor5.run_critical_speed(num_modes=12, rtol=0.005)
-    results6 = rotor6.run_critical_speed(num_modes=12, rtol=0.005)
-
-    wn5 = np.array(
-        [
-            86.10505193,
-            86.60492546,
-            198.93259257,
-            207.97165539,
-            244.95609413,
-            250.53522782,
-        ]
-    )
-    wd5 = np.array(
-        [
-            86.1050519,
-            86.60492544,
-            198.93259256,
-            207.97165539,
-            244.95609413,
-            250.53522782,
-        ]
-    )
-    log_dec5 = np.zeros_like(wd5)
-    damping_ratio5 = np.zeros_like(wd5)
-
-    wd6 = np.array(
-        [
-            61.52110644,
-            63.72862198,
-            117.49491374,
-            118.55829416,
-            233.83724523,
-            236.40346235,
-        ]
-    )
-    wn6 = np.array(
-        [
-            61.52110644,
-            63.72862198,
-            117.49491375,
-            118.55829421,
-            233.83724523,
-            236.40346458,
-        ]
-    )
-    log_dec6 = np.zeros_like(wd6)
-    damping_ratio6 = np.zeros_like(wd6)
-
-    assert_almost_equal(results5.wd, wd5, decimal=4)
-    assert_almost_equal(results5.wn, wn5, decimal=4)
-    assert_almost_equal(results5.log_dec, log_dec5, decimal=4)
-    assert_almost_equal(results5.damping_ratio, damping_ratio5, decimal=4)
-
-    assert_almost_equal(results6.wd, wd6, decimal=4)
-    assert_almost_equal(results6.wn, wn6, decimal=4)
-    assert_almost_equal(results6.log_dec, log_dec6, decimal=4)
-    assert_almost_equal(results6.damping_ratio, damping_ratio6, decimal=4)
-
-
 @pytest.fixture
 def coaxrotor():
     #  Co-axial rotor system with 2 shafts, 4 disks and
@@ -1515,10 +1456,11 @@ def test_H_kappa(rotor7):
 
 
 def test_save_load():
+
     a = rotor_example()
     a.save("teste00000000000000001")
-    b = Rotor.load("teste00000000000000001.rsm")
-    (Path.cwd() / "teste00000000000000001.rsm").unlink()
+    b = Rotor.load("teste00000000000000001")
+    Rotor.remove("teste00000000000000001")
 
     assert a == b
 
@@ -1717,59 +1659,3 @@ def test_modal_6dof(rotor_6dof):
 
     assert_almost_equal(modal.wn[:6], wn, decimal=3)
     assert_almost_equal(modal.wd[:6], wd, decimal=3)
-
-
-@pytest.fixture
-def rotor8():
-    #  Rotor with damping
-    #  Rotor with 6 shaft elements, 2 disks and 2 bearings with frequency dependent coefficients
-    i_d = 0
-    o_d = 0.05
-    n = 6
-    L = [0.25 for _ in range(n)]
-
-    shaft_elem = [
-        ShaftElement(
-            l,
-            i_d,
-            o_d,
-            material=steel,
-            shear_effects=True,
-            rotary_inertia=True,
-            gyroscopic=True,
-        )
-        for l in L
-    ]
-
-    disk0 = DiskElement.from_geometry(2, steel, 0.07, 0.05, 0.28)
-    disk1 = DiskElement.from_geometry(4, steel, 0.07, 0.05, 0.35)
-
-    stfx = [1e7, 1.5e7]
-    stfy = [1e7, 1.5e7]
-    c = [1e3, 1.5e3]
-    frequency = [50, 5000]
-    bearing0 = BearingElement(0, kxx=stfx, kyy=stfy, cxx=c, cyy=c, frequency=frequency)
-    bearing1 = BearingElement(6, kxx=stfx, kyy=stfy, cxx=c, cyy=c, frequency=frequency)
-
-    return Rotor(shaft_elem, [disk0, disk1], [bearing0, bearing1])
-
-
-def test_ucs_calc(rotor8):
-    exp_stiffness_range = np.array([1000000.0, 1832980.710832, 3359818.286284])
-    exp_rotor_wn = np.array([86.658114, 95.660573, 101.868254])
-    exp_intersection_points_x = np.array(
-        [9632232.337582, 9632232.337582, 10333174.417919]
-    )
-    exp_intersection_points_y = np.array([107.920792, 107.920792, 411.881188])
-    stiffness_range, rotor_wn, bearing, intersection_points = rotor8._calc_ucs()
-    assert_allclose(stiffness_range[:3], exp_stiffness_range)
-    assert_allclose(rotor_wn[0, :3], exp_rotor_wn)
-    assert_allclose(intersection_points["x"][:3], exp_intersection_points_x, rtol=1e-3)
-    assert_allclose(intersection_points["y"][:3], exp_intersection_points_y, rtol=1e-3)
-
-
-def test_save_load(rotor8):
-    rotor8.save("/tmp/rotor8.toml")
-    rotor8_loaded = Rotor.load("/tmp/rotor8.toml")
-
-    rotor8 == rotor8_loaded
diff --git a/ross/tests/test_shaft_element.py b/ross/tests/test_shaft_element.py
index d8621908..22392e23 100644
--- a/ross/tests/test_shaft_element.py
+++ b/ross/tests/test_shaft_element.py
@@ -308,15 +308,6 @@ def test_match_gyroscopic_matrix(tap2, tim2):
     assert_almost_equal(G_tap, G_tim, decimal=5)
 
 
-def test_save_load(tim, tap_tim):
-    tim.save("/tmp/tim.toml")
-    tim_loaded = ShaftElement.load("/tmp/tim.toml")
-    assert tim == tim_loaded
-    tap_tim.save("/tmp/tap_tim.toml")
-    tap_tim_loaded = ShaftElement.load("/tmp/tap_tim.toml")
-    assert tap_tim == tap_tim_loaded
-
-
 @pytest.fixture
 def s6_eb():