Repository: ross-rotordynamics/ross
Pull request #582 (2020-06-11 21:30:26)

Add plotly theme
This PR is related to #580.
We want to create a plotly theme to avoid having to set every style parameter for each plot, saving us lots of lines of code, decreasing burden in maintaining the code and keeping the style consistent across different plots.


Patch:
diff --git a/ross/__init__.py b/ross/__init__.py
index 3093ba63..bed038f3 100644
--- a/ross/__init__.py
+++ b/ross/__init__.py
@@ -1,5 +1,9 @@
 __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 *
@@ -8,3 +12,5 @@ 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/fluid_flow/fluid_flow_graphics.py b/ross/fluid_flow/fluid_flow_graphics.py
index 21969673..69d69be7 100644
--- a/ross/fluid_flow/fluid_flow_graphics.py
+++ b/ross/fluid_flow/fluid_flow_graphics.py
@@ -1,11 +1,10 @@
 import numpy as np
 import plotly.graph_objects as go
-import plotly.io as pio
 
-pio.renderers.default = "browser"
+from ross.plotly_theme import tableau_colors
 
 
-def plot_eccentricity(fluid_flow_object, z=0, **kwargs):
+def plot_eccentricity(fluid_flow_object, z=0, fig=None, **kwargs):
     """Plot the rotor eccentricity.
 
     This function assembles pressure graphic along the z-axis.
@@ -16,6 +15,8 @@ def plot_eccentricity(fluid_flow_object, z=0, **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, ...).
@@ -34,29 +35,15 @@ def plot_eccentricity(fluid_flow_object, z=0, **kwargs):
     >>> # to show the plots you can use:
     >>> # show(fig)
     """
-    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()
+    if fig is None:
+        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(size=10, color="firebrick"),
-            line=dict(width=2.0, color="firebrick"),
+            marker=dict(color=tableau_colors["red"]),
+            line=dict(color=tableau_colors["red"]),
             name="Stator",
             legendgroup="Stator",
             hovertemplate=("<b>X: %{x:.3e}</b><br>" + "<b>Y: %{y:.3e}</b>"),
@@ -67,8 +54,8 @@ def plot_eccentricity(fluid_flow_object, z=0, **kwargs):
             x=fluid_flow_object.xri[z],
             y=fluid_flow_object.yri[z],
             mode="markers+lines",
-            marker=dict(size=10, color="royalblue"),
-            line=dict(width=2.0, color="royalblue"),
+            marker=dict(color=tableau_colors["blue"]),
+            line=dict(color=tableau_colors["blue"]),
             name="Rotor",
             legendgroup="Rotor",
             hovertemplate=("<b>X: %{x:.3e}</b><br>" + "<b>Y: %{y:.3e}</b>"),
@@ -78,7 +65,7 @@ def plot_eccentricity(fluid_flow_object, z=0, **kwargs):
         go.Scatter(
             x=[fluid_flow_object.xi],
             y=[fluid_flow_object.yi],
-            marker=dict(size=10, color="firebrick"),
+            marker=dict(color=tableau_colors["red"]),
             name="Stator",
             legendgroup="Stator",
             showlegend=False,
@@ -89,42 +76,24 @@ def plot_eccentricity(fluid_flow_object, z=0, **kwargs):
         go.Scatter(
             x=[0],
             y=[0],
-            marker=dict(size=10, color="royalblue"),
+            marker=dict(color=tableau_colors["blue"]),
             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>",
-        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_xaxes(title_text="<b>X axis</b>")
+    fig.update_yaxes(title_text="<b>Y axis</b>")
     fig.update_layout(
-        title=dict(text="<b>Cut in plane Z={}</b>".format(z), font=dict(size=16)),
-        **kwargs,
+        title=dict(text="<b>Cut in plane Z={}</b>".format(z)),
+        yaxis=dict(scaleanchor="x", scaleratio=1, **kwargs),
     )
 
     return fig
 
 
-def plot_pressure_z(fluid_flow_object, theta=0, **kwargs):
+def plot_pressure_z(fluid_flow_object, theta=0, fig=None, **kwargs):
     """Plot the pressure distribution along the z-axis.
 
     This function assembles pressure graphic along the z-axis for one or both the
@@ -136,6 +105,8 @@ def plot_pressure_z(fluid_flow_object, theta=0, **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, ...).
@@ -164,29 +135,16 @@ def plot_pressure_z(fluid_flow_object, theta=0, **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)
 
-    fig = go.Figure()
+    if fig is None:
+        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(width=3.0, color="royalblue"),
+                line=dict(color=tableau_colors["blue"]),
                 showlegend=True,
                 name="<b>Numerical pressure</b>",
                 hovertemplate=(
@@ -201,7 +159,7 @@ def plot_pressure_z(fluid_flow_object, theta=0, **kwargs):
                 x=fluid_flow_object.z_list,
                 y=fluid_flow_object.p_mat_analytical[:, theta],
                 mode="lines",
-                line=dict(width=3.0, color="firebrick"),
+                line=dict(color=tableau_colors["red"]),
                 showlegend=True,
                 name="<b>Analytical pressure</b>",
                 hovertemplate=(
@@ -210,33 +168,14 @@ def plot_pressure_z(fluid_flow_object, theta=0, **kwargs):
                 ),
             )
         )
-    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_xaxes(title_text="<b>Axial Length</b>")
+    fig.update_yaxes(title_text="<b>Pressure</b>")
     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,
     )
@@ -244,7 +183,7 @@ def plot_pressure_z(fluid_flow_object, theta=0, **kwargs):
     return fig
 
 
-def plot_shape(fluid_flow_object, theta=0, **kwargs):
+def plot_shape(fluid_flow_object, theta=0, fig=None, **kwargs):
     """Plot the surface geometry of the rotor.
 
     This function assembles a graphic representing the geometry of the rotor.
@@ -254,6 +193,8 @@ def plot_shape(fluid_flow_object, theta=0, **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
     -------
@@ -272,28 +213,14 @@ def plot_shape(fluid_flow_object, theta=0, **kwargs):
     >>> # to show the plots you can use:
     >>> # fig.show()
     """
-    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 fig is None:
+        fig = go.Figure()
     fig.add_trace(
         go.Scatter(
             x=fluid_flow_object.z_list,
             y=fluid_flow_object.re[:, theta],
             mode="lines",
-            line=dict(width=3.0, color="firebrick"),
+            line=dict(color=tableau_colors["red"]),
             showlegend=True,
             hoverinfo="none",
             name="<b>Stator</b>",
@@ -304,39 +231,20 @@ def plot_shape(fluid_flow_object, theta=0, **kwargs):
             x=fluid_flow_object.z_list,
             y=fluid_flow_object.ri[:, theta],
             mode="lines",
-            line=dict(width=3.0, color="royalblue"),
+            line=dict(color=tableau_colors["blue"]),
             showlegend=True,
             hoverinfo="none",
             name="<b>Rotor</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_xaxes(title_text="<b>Axial Length</b>")
+    fig.update_yaxes(title_text="<b>Radial direction</b>")
     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,
     )
@@ -344,7 +252,7 @@ def plot_shape(fluid_flow_object, theta=0, **kwargs):
     return fig
 
 
-def plot_pressure_theta(fluid_flow_object, z=0, **kwargs):
+def plot_pressure_theta(fluid_flow_object, z=0, fig=None, **kwargs):
     """Plot the pressure distribution along theta.
 
     This function assembles pressure graphic in the theta direction for a given z
@@ -356,6 +264,8 @@ def plot_pressure_theta(fluid_flow_object, z=0, **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, ...).
@@ -384,29 +294,15 @@ def plot_pressure_theta(fluid_flow_object, z=0, **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)
-
-    fig = go.Figure()
+    if fig is None:
+        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(width=3.0, color="royalblue"),
+                line=dict(color=tableau_colors["blue"]),
                 showlegend=True,
                 name="<b>Numerical pressure</b>",
                 hovertemplate=(
@@ -420,7 +316,7 @@ def plot_pressure_theta(fluid_flow_object, z=0, **kwargs):
                 x=fluid_flow_object.gama[z],
                 y=fluid_flow_object.p_mat_analytical[z],
                 mode="lines",
-                line=dict(width=3.0, color="firebrick"),
+                line=dict(color=tableau_colors["red"]),
                 showlegend=True,
                 name="<b>Analytical pressure</b>",
                 hovertemplate=(
@@ -430,38 +326,17 @@ def plot_pressure_theta(fluid_flow_object, z=0, **kwargs):
             )
         )
 
-    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_xaxes(title_text="<b>Theta value</b>")
+    fig.update_yaxes(title_text="<b>Pressure</b>")
     fig.update_layout(
-        title=dict(
-            text=("<b>Pressure along Theta | Z={}<b>".format(z)), font=dict(size=16)
-        ),
-        **kwargs,
+        title=dict(text=("<b>Pressure along Theta | Z={}<b>".format(z))), **kwargs
     )
 
     return fig
 
 
 def plot_pressure_theta_cylindrical(
-    fluid_flow_object, z=0, from_numerical=True, **kwargs
+    fluid_flow_object, z=0, from_numerical=True, fig=None, **kwargs
 ):
     """Plot cylindrical pressure graphic in the theta direction.
 
@@ -478,6 +353,8 @@ 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, ...).
@@ -544,11 +421,9 @@ def plot_pressure_theta_cylindrical(
                 continue
             z_matrix[i][j] = pressure_along_theta[i] - min_pressure + 0.01
 
-    kwargs_default_values = dict(width=1200, height=900)
-    for k, v in kwargs_default_values.items():
-        kwargs.setdefault(k, v)
-
-    fig = go.Figure(
+    if fig is None:
+        fig = go.Figure()
+    fig.add_trace(
         go.Barpolar(
             r=r_matrix.ravel(),
             theta=theta_matrix.ravel(),
@@ -558,10 +433,7 @@ 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", font=dict(size=16)),
-                    tickfont=dict(size=16),
-                ),
+                colorbar=dict(title=dict(text="<b>Pressure</b>", side="top")),
             ),
             thetaunit="degrees",
             name="Pressure",
@@ -576,14 +448,9 @@ 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,
-                gridcolor="lightgray",
-                linecolor="black",
-                linewidth=2.5,
+                rotation=-90 - fluid_flow_object.attitude_angle * 180 / np.pi
             ),
         ),
         **kwargs,
@@ -591,12 +458,15 @@ def plot_pressure_theta_cylindrical(
     return fig
 
 
-def plot_pressure_surface(fluid_flow_object, **kwargs):
+def plot_pressure_surface(fluid_flow_object, fig=None, **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, ...).
@@ -625,12 +495,8 @@ def plot_pressure_surface(fluid_flow_object, **kwargs):
             "Must calculate the pressure matrix. "
             "Try calling calculate_pressure_matrix_numerical() or calculate_pressure_matrix_analytical() first."
         )
-
-    kwargs_default_values = dict(width=1200, height=900)
-    for k, v in kwargs_default_values.items():
-        kwargs.setdefault(k, v)
-
-    fig = go.Figure()
+    if fig is None:
+        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(
@@ -641,33 +507,7 @@ def plot_pressure_surface(fluid_flow_object, **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", 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),
-                ),
+                colorbar=dict(title=dict(text="<b>Pressure</b>", side="top")),
                 name="Pressure",
                 showlegend=False,
                 hovertemplate=(
@@ -681,32 +521,11 @@ def plot_pressure_surface(fluid_flow_object, **kwargs):
     fig.update_layout(
         scene=dict(
             bgcolor="white",
-            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,
-            ),
+            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>")),
         ),
-        title=dict(text=("<b>Bearing Pressure Field</b>"), font=dict(size=20),),
+        title=dict(text="<b>Bearing Pressure Field</b>"),
         **kwargs,
     )
 
diff --git a/ross/plotly_theme.py b/ross/plotly_theme.py
new file mode 100644
index 00000000..b038bf47
--- /dev/null
+++ b/ross/plotly_theme.py
@@ -0,0 +1,377 @@
+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/results.py b/ross/results.py
index 4228aa58..1594dc20 100644
--- a/ross/results.py
+++ b/ross/results.py
@@ -2,15 +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
 
-# set Plotly palette of colors
-colors1 = px.colors.qualitative.Dark24
+from ross.plotly_theme import tableau_colors
 
 
 class CriticalSpeedResults:
@@ -479,7 +479,7 @@ class ModalResults:
 
         return xn, yn, zn, x_circles, y_circles, z_circles_pos, nn
 
-    def plot_mode3D(self, mode=None, evec=None, **kwargs):
+    def plot_mode_3d(self, mode=None, evec=None, fig=None, **kwargs):
         """Plot (3D view) the mode shapes.
 
         Parameters
@@ -489,6 +489,8 @@ 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, ...).
@@ -499,7 +501,8 @@ class ModalResults:
         fig : Plotly graph_objects.Figure()
             The figure object with the plot.
         """
-        fig = go.Figure()
+        if fig is None:
+            fig = go.Figure()
 
         nodes = self.nodes
         kappa_mode = self.kappa_modes[mode]
@@ -511,9 +514,24 @@ class ModalResults:
                     x=zc_pos[10:, node],
                     y=xc[10:, node],
                     z=yc[10:, node],
-                    mode="markers+lines",
-                    line=dict(width=4.0, color=kappa_mode[node]),
-                    marker=dict(size=5, color=kappa_mode[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]),
                     name="node {}".format(node),
                     showlegend=False,
                 )
@@ -525,7 +543,7 @@ class ModalResults:
                 y=xn,
                 z=yn,
                 mode="lines",
-                line=dict(width=3.0, color="black", dash="dash"),
+                line=dict(color="black", dash="dash"),
                 name="mode shape",
                 showlegend=False,
             )
@@ -541,63 +559,43 @@ class ModalResults:
                 y=zn_cl * 0,
                 z=zn_cl * 0,
                 mode="lines",
-                line=dict(width=2.0, color="black", dash="dashdot"),
+                line=dict(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>", font=dict(size=14)),
-                    tickfont=dict(size=16),
-                    range=[zn_cl0 - 0.1, zn_cl1 + 0.1],
+                    title=dict(text="<b>Rotor Length</b>"),
+                    autorange="reversed",
                     nticks=5,
-                    backgroundcolor="lightgray",
-                    gridcolor="white",
-                    showspikes=False,
                 ),
                 yaxis=dict(
-                    title=dict(
-                        text="<b>Dimensionless deformation</b>", font=dict(size=14)
-                    ),
-                    tickfont=dict(size=16),
+                    title=dict(text="<b>Relative Displacement</b>"),
                     range=[-2, 2],
                     nticks=5,
-                    backgroundcolor="lightgray",
-                    gridcolor="white",
-                    showspikes=False,
                 ),
                 zaxis=dict(
-                    title=dict(
-                        text="<b>Dimensionless deformation</b>", font=dict(size=14)
-                    ),
-                    tickfont=dict(size=16),
+                    title=dict(text="<b>Relative Displacement</b>"),
                     range=[-2, 2],
                     nticks=5,
-                    backgroundcolor="lightgray",
-                    gridcolor="white",
-                    showspikes=False,
                 ),
             ),
             title=dict(
                 text=(
-                    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>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>log dec</b> = {self.log_dec[mode]:.1f}"
-                ),
-                font=dict(size=14),
+                )
             ),
             **kwargs,
         )
 
         return fig
 
-    def plot_mode2D(self, mode=None, evec=None, **kwargs):
+    def plot_mode_2d(self, mode=None, evec=None, fig=None, **kwargs):
         """Plot (2D view) the mode shapes.
 
         Parameters
@@ -607,6 +605,8 @@ 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,9 +633,10 @@ class ModalResults:
         zn = np.delete(zn, idx_remove)
         vn = np.delete(vn, idx_remove)
 
-        fig = go.Figure()
+        if fig is None:
+            fig = go.Figure()
 
-        colors = dict(Backward="firebrick", Mixed="black", Forward="royalblue")
+        colors = dict(Backward="red", Mixed="black", Forward="blue")
         whirl_dir = colors[self.whirl_direction()[mode]]
 
         fig.add_trace(
@@ -643,7 +644,7 @@ class ModalResults:
                 x=zn,
                 y=vn,
                 mode="lines",
-                line=dict(width=2.5, color=whirl_dir),
+                line=dict(color=whirl_dir),
                 name="mode shape",
                 showlegend=False,
             )
@@ -654,51 +655,23 @@ class ModalResults:
                 x=nodes_pos,
                 y=np.zeros(len(nodes_pos)),
                 mode="lines",
-                line=dict(width=1.0, color="black", dash="dashdot"),
+                line=dict(color="black", dash="dashdot"),
                 name="centerline",
                 hoverinfo="none",
                 showlegend=False,
             )
         )
 
-        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_xaxes(title_text="<b>Rotor Length</b>")
+        fig.update_yaxes(title_text="<b>Relative Displacement</b>")
         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,
         )
@@ -735,7 +708,7 @@ class CampbellResults:
         self.log_dec = log_dec
         self.whirl_values = whirl_values
 
-    def plot(self, harmonics=[1], **kwargs):
+    def plot(self, harmonics=[1], fig=None, **kwargs):
         """Create Campbell Diagram figure using Plotly.
 
         Parameters
@@ -743,6 +716,8 @@ 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, ...).
@@ -760,13 +735,13 @@ class CampbellResults:
         speed_range = self.speed_range
         log_dec_map = log_dec.flatten()
 
-        fig = go.Figure()
+        if fig is None:
+            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]
@@ -799,7 +774,7 @@ class CampbellResults:
                                 x=xnew[idx],
                                 y=ynew[idx],
                                 mode="markers",
-                                marker=dict(symbol="x", size=10, color="black"),
+                                marker=dict(symbol="x", color="black"),
                                 name="Crit. Speed",
                                 legendgroup="Crit. Speed",
                                 showlegend=False,
@@ -820,7 +795,6 @@ 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],
@@ -841,7 +815,7 @@ class CampbellResults:
                     x=speed_range,
                     y=h * speed_range,
                     mode="lines",
-                    line=dict(width=2.5, color=colors1[j], dash="dashdot"),
+                    line=dict(dash="dashdot"),
                     name="{}x speed".format(h),
                     hoverinfo="none",
                 )
@@ -852,55 +826,36 @@ class CampbellResults:
         for mark, legend in zip(scatter_marker, legends):
             fig.add_trace(
                 go.Scatter(
-                    x=[-1000],
-                    y=[-1000],
+                    x=[0],
+                    y=[0],
                     mode="markers",
                     name=legend,
                     legendgroup=legend,
-                    marker=dict(symbol=mark, size=16, color="black"),
+                    marker=dict(symbol=mark, color="black"),
                 )
             )
 
         fig.update_xaxes(
             title_text="<b>Frequency</b>",
-            title_font=dict(family="Arial", size=20),
-            tickfont=dict(size=16),
             range=[np.min(speed_range), 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>",
-            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,
+            title_text="<b>Damped Natural Frequencies</b>", range=[0, 1.1 * np.max(wd)]
         )
         fig.update_layout(
-            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", font=dict(size=20)),
-                    tickfont=dict(size=16),
-                ),
+                colorbar=dict(title=dict(text="<b>Log Dec</b>", side="right")),
             ),
             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,
         )
@@ -936,7 +891,7 @@ class FrequencyResponseResults:
         self.magnitude = magnitude
         self.phase = phase
 
-    def plot_magnitude(self, inp, out, units="mic-pk-pk", **mag_kwargs):
+    def plot_magnitude(self, inp, out, units="mic-pk-pk", fig=None, **mag_kwargs):
         """Plot frequency response (magnitude) using Plotly.
 
         This method plots the frequency response magnitude given an output and
@@ -951,6 +906,8 @@ 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, ...).
@@ -971,20 +928,15 @@ class FrequencyResponseResults:
         else:
             y_axis_label = "<b>Amplitude (dB)</b>"
 
-        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()
+        if fig is None:
+            fig = go.Figure()
 
         fig.add_trace(
             go.Scatter(
                 x=frequency_range,
                 y=mag[inp, out, :],
                 mode="lines",
-                line=dict(width=3.0, color="royalblue"),
+                line=dict(color=tableau_colors["blue"]),
                 name="Amplitude",
                 legendgroup="Amplitude",
                 showlegend=False,
@@ -994,38 +946,14 @@ class FrequencyResponseResults:
 
         fig.update_xaxes(
             title_text="<b>Frequency</b>",
-            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,
+            range=[np.min(frequency_range), np.max(frequency_range)],
         )
+        fig.update_yaxes(title_text=y_axis_label)
+        fig.update_layout(**mag_kwargs)
 
         return fig
 
-    def plot_phase(self, inp, out, **phase_kwargs):
+    def plot_phase(self, inp, out, fig=None, **phase_kwargs):
         """Plot frequency response (phase) using Plotly.
 
         This method plots the frequency response phase given an output and
@@ -1037,6 +965,8 @@ 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, ...).
@@ -1050,20 +980,15 @@ class FrequencyResponseResults:
         frequency_range = self.speed_range
         phase = self.phase
 
-        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()
+        if fig is None:
+            fig = go.Figure()
 
         fig.add_trace(
             go.Scatter(
                 x=frequency_range,
                 y=phase[inp, out, :],
                 mode="lines",
-                line=dict(width=3.0, color="royalblue"),
+                line=dict(color=tableau_colors["blue"]),
                 name="Phase",
                 legendgroup="Phase",
                 showlegend=False,
@@ -1073,38 +998,14 @@ class FrequencyResponseResults:
 
         fig.update_xaxes(
             title_text="<b>Frequency</b>",
-            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,
+            range=[np.min(frequency_range), np.max(frequency_range)],
         )
+        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", **polar_kwargs):
+    def plot_polar_bode(self, inp, out, units="mic-pk-pk", fig=None, **polar_kwargs):
         """Plot frequency response (polar) using Plotly.
 
         This method plots the frequency response (polar graph) given an output and
@@ -1119,6 +1020,8 @@ 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, ...).
@@ -1140,13 +1043,8 @@ class FrequencyResponseResults:
         else:
             r_axis_label = "<b>Amplitude (dB)</b>"
 
-        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()
+        if fig is None:
+            fig = go.Figure()
 
         fig.add_trace(
             go.Scatterpolar(
@@ -1155,8 +1053,8 @@ class FrequencyResponseResults:
                 customdata=frequency_range,
                 thetaunit="radians",
                 mode="lines+markers",
-                marker=dict(size=8, color="royalblue"),
-                line=dict(width=3.0, color="royalblue"),
+                marker=dict(color=tableau_colors["blue"]),
+                line=dict(color=tableau_colors["blue"]),
                 name="Polar_plot",
                 legendgroup="Polar",
                 showlegend=False,
@@ -1170,18 +1068,7 @@ class FrequencyResponseResults:
 
         fig.update_layout(
             polar=dict(
-                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,
-                ),
+                radialaxis=dict(title=dict(text=r_axis_label), exponentformat="power")
             ),
             **polar_kwargs,
         )
@@ -1193,10 +1080,10 @@ class FrequencyResponseResults:
         inp,
         out,
         units="mic-pk-pk",
-        mag_kwargs={},
-        phase_kwargs={},
-        polar_kwargs={},
-        subplot_kwargs={},
+        mag_kwargs=None,
+        phase_kwargs=None,
+        polar_kwargs=None,
+        subplot_kwargs=None,
     ):
         """Plot frequency response.
 
@@ -1245,15 +1132,10 @@ class FrequencyResponseResults:
             Plotly figure with Amplitude vs Frequency and Phase vs Frequency and
             polar Amplitude vs Phase plots.
         """
-        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)
+        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)
 
         fig0 = self.plot_magnitude(inp, out, units, **mag_kwargs)
         fig1 = self.plot_phase(inp, out, **phase_kwargs)
@@ -1312,7 +1194,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
@@ -1321,7 +1203,7 @@ class ForcedResponseResults:
         self.phase = phase
         self.unbalance = unbalance
 
-    def plot_magnitude(self, dof, units="m", **kwargs):
+    def plot_magnitude(self, dof, units="m", fig=None, **kwargs):
         """Plot forced response (magnitude) using Plotly.
 
         Parameters
@@ -1331,6 +1213,8 @@ 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, ...).
@@ -1350,20 +1234,15 @@ class ForcedResponseResults:
             mag = 2 * mag * 1e6
             y_axis_label = "<b>Amplitude (μ pk-pk)</b>"
 
-        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()
+        if fig is None:
+            fig = go.Figure()
 
         fig.add_trace(
             go.Scatter(
                 x=frequency_range,
                 y=mag[dof],
                 mode="lines",
-                line=dict(width=3.0, color="royalblue"),
+                line=dict(color=tableau_colors["blue"]),
                 name="Amplitude",
                 legendgroup="Amplitude",
                 showlegend=False,
@@ -1373,45 +1252,22 @@ class ForcedResponseResults:
 
         fig.update_xaxes(
             title_text="<b>Frequency</b>",
-            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,
+            range=[np.min(frequency_range), np.max(frequency_range)],
         )
+        fig.update_yaxes(title_text=y_axis_label, exponentformat="power")
+        fig.update_layout(**kwargs)
 
         return fig
 
-    def plot_phase(self, dof, **kwargs):
+    def plot_phase(self, dof, fig=None, **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, ...).
@@ -1425,20 +1281,15 @@ class ForcedResponseResults:
         frequency_range = self.speed_range
         phase = self.phase
 
-        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()
+        if fig is None:
+            fig = go.Figure()
 
         fig.add_trace(
             go.Scatter(
                 x=frequency_range,
                 y=phase[dof],
                 mode="lines",
-                line=dict(width=3.0, color="royalblue"),
+                line=dict(color=tableau_colors["blue"]),
                 name="Phase",
                 legendgroup="Phase",
                 showlegend=False,
@@ -1448,38 +1299,14 @@ class ForcedResponseResults:
 
         fig.update_xaxes(
             title_text="<b>Frequency</b>",
-            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,
+            range=[np.min(frequency_range), np.max(frequency_range)],
         )
+        fig.update_yaxes(title_text="<b>Phase Angle</b>")
+        fig.update_layout(**kwargs)
 
         return fig
 
-    def plot_polar_bode(self, dof, units="m", **kwargs):
+    def plot_polar_bode(self, dof, units="m", fig=None, **kwargs):
         """Plot polar forced response using Plotly.
 
         Parameters
@@ -1489,6 +1316,8 @@ 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, ...).
@@ -1510,13 +1339,8 @@ class ForcedResponseResults:
         else:
             r_axis_label = "<b>Amplitude (dB)</b>"
 
-        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()
+        if fig is None:
+            fig = go.Figure()
 
         fig.add_trace(
             go.Scatterpolar(
@@ -1525,8 +1349,8 @@ class ForcedResponseResults:
                 customdata=frequency_range,
                 thetaunit="radians",
                 mode="lines+markers",
-                marker=dict(size=8, color="royalblue"),
-                line=dict(width=3.0, color="royalblue"),
+                marker=dict(color=tableau_colors["blue"]),
+                line=dict(color=tableau_colors["blue"]),
                 name="Polar_plot",
                 legendgroup="Polar",
                 showlegend=False,
@@ -1540,18 +1364,7 @@ class ForcedResponseResults:
 
         fig.update_layout(
             polar=dict(
-                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,
-                ),
+                radialaxis=dict(title_text=r_axis_label, exponentformat="power")
             ),
             **kwargs,
         )
@@ -1562,10 +1375,10 @@ class ForcedResponseResults:
         self,
         dof,
         units="m",
-        mag_kwargs={},
-        phase_kwargs={},
-        polar_kwargs={},
-        subplot_kwargs={},
+        mag_kwargs=None,
+        phase_kwargs=None,
+        polar_kwargs=None,
+        subplot_kwargs=None,
     ):
         """Plot forced response.
 
@@ -1609,15 +1422,10 @@ class ForcedResponseResults:
             Plotly figure with Amplitude vs Frequency and Phase vs Frequency and
             polar Amplitude vs Phase plots.
         """
-        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)
+        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)
 
         fig0 = self.plot_magnitude(dof, units, **mag_kwargs)
         fig1 = self.plot_phase(dof, **phase_kwargs)
@@ -1756,7 +1564,7 @@ class ForcedResponseResults:
 
         return Mx, My
 
-    def plot_deflected_shape_2d(self, speed, units="mic-pk-pk", **kwargs):
+    def plot_deflected_shape_2d(self, speed, units="mic-pk-pk", fig=None, **kwargs):
         """Plot the 2D deflected shape diagram.
 
         Parameters
@@ -1771,6 +1579,8 @@ 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, ...).
@@ -1784,22 +1594,16 @@ 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)
 
-        fig = go.Figure()
+        if fig is None:
+            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,
@@ -1814,45 +1618,21 @@ class ForcedResponseResults:
                 x=nodes_pos,
                 y=np.zeros(len(nodes_pos)),
                 mode="lines",
-                line=dict(width=2.0, color="black", dash="dashdot"),
+                line=dict(color="black", dash="dashdot"),
                 showlegend=False,
                 hoverinfo="none",
             )
         )
 
-        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,
-        )
+        fig.update_xaxes(title_text="<b>Rotor Length</b>")
+        fig.update_yaxes(title_text="<b>Major Axis Absolute Amplitude</b>")
+        fig.update_layout(**kwargs)
 
         return fig
 
-    def plot_deflected_shape_3d(self, speed, samples=101, units="mic-pk-pk", **kwargs):
+    def plot_deflected_shape_3d(
+        self, speed, samples=101, units="mic-pk-pk", fig=None, **kwargs
+    ):
         """Plot the 3D deflected shape diagram.
 
         Parameters
@@ -1870,6 +1650,8 @@ 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, ...).
@@ -1883,10 +1665,6 @@ 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
@@ -1900,7 +1678,8 @@ class ForcedResponseResults:
 
         x_pos = np.repeat(nodes_pos, t.size).reshape(len(nodes_pos), t.size)
 
-        fig = go.Figure()
+        if fig is None:
+            fig = go.Figure()
         for i, n in enumerate(nodes):
             dofx = number_dof * n
             dofy = number_dof * n + 1
@@ -1915,7 +1694,7 @@ class ForcedResponseResults:
                     y=y,
                     z=z,
                     mode="lines",
-                    line=dict(width=6.0, color="royalblue"),
+                    line=dict(color="royalblue"),
                     name="Orbit",
                     legendgroup="Orbit",
                     showlegend=False,
@@ -1936,8 +1715,8 @@ class ForcedResponseResults:
                 y=np.real(maj_vect[3]),
                 z=np.imag(maj_vect[3]),
                 mode="lines+markers",
-                marker=dict(size=4.0, color="black"),
-                line=dict(width=6.0, color="black", dash="dashdot"),
+                marker=dict(color="black"),
+                line=dict(color="black", dash="dashdot"),
                 name="Major Axis",
                 legendgroup="Major_Axis",
                 showlegend=True,
@@ -1957,7 +1736,7 @@ class ForcedResponseResults:
                 y=line,
                 z=line,
                 mode="lines",
-                line=dict(width=2.0, color="black", dash="dashdot"),
+                line=dict(color="black", dash="dashdot"),
                 showlegend=False,
                 hoverinfo="none",
             )
@@ -1972,7 +1751,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(width=6.0, color="firebrick"),
+                    line=dict(color="firebrick"),
                     legendgroup="Unbalance",
                     hoverinfo="none",
                     showlegend=False,
@@ -1984,7 +1763,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", size=5.0, color="firebrick"),
+                    marker=dict(symbol="diamond", color="firebrick"),
                     name="Unbalance",
                     legendgroup="Unbalance",
                     showlegend=True if i == 0 else False,
@@ -1996,51 +1775,19 @@ class ForcedResponseResults:
             i += 1
 
         fig.update_layout(
-            width=1200,
-            height=900,
             scene=dict(
                 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,
-                ),
-            ),
-            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,
+                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>")),
             ),
+            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", **kwargs):
+    def plot_bending_moment(self, speed, units="mic-pk-pk", fig=None, **kwargs):
         """Plot the bending moment diagram.
 
         Parameters
@@ -2068,59 +1815,51 @@ 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
 
-        fig = go.Figure()
+        if fig is None:
+            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
@@ -2129,41 +1868,15 @@ class ForcedResponseResults:
                 x=nodes_pos,
                 y=np.zeros_like(nodes_pos),
                 mode="lines",
-                line=dict(width=3.0, color="black", dash="dashdot"),
+                line=dict(color="black", dash="dashdot"),
                 showlegend=False,
                 hoverinfo="none",
             )
         )
 
-        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,
-        )
+        fig.update_xaxes(title_text="<b>Rotor Length</b>")
+        fig.update_yaxes(title_text="<b>Bending Moment</b>")
+        fig.update_layout(**kwargs)
 
         return fig
 
@@ -2172,10 +1885,10 @@ class ForcedResponseResults:
         speed,
         samples=101,
         units="mic-pk-pk",
-        shape2d_kwargs={},
-        shape3d_kwargs={},
-        bm_kwargs={},
-        subplot_kwargs={},
+        shape2d_kwargs=None,
+        shape3d_kwargs=None,
+        bm_kwargs=None,
+        subplot_kwargs=None,
     ):
         """Plot deflected shape diagrams.
 
@@ -2223,11 +1936,10 @@ class ForcedResponseResults:
             Plotly figure with Amplitude vs Frequency and Phase vs Frequency and
             polar Amplitude vs Phase plots.
         """
-        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)
+        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)
 
         fig0 = self.plot_deflected_shape_2d(speed, units, **shape2d_kwargs)
         fig1 = self.plot_deflected_shape_3d(speed, samples, units, **shape3d_kwargs)
@@ -2317,7 +2029,7 @@ class StaticResults:
         self.nodes_pos = nodes_pos
         self.Vx_axis = Vx_axis
 
-    def plot_deformation(self, **kwargs):
+    def plot_deformation(self, fig=None, **kwargs):
         """Plot the shaft static deformation.
 
         This method plots:
@@ -2325,6 +2037,8 @@ 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, ...).
@@ -2335,13 +2049,8 @@ class StaticResults:
         fig : Plotly graph_objects.Figure()
             The figure object with the plot.
         """
-        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()
+        if fig is None:
+            fig = go.Figure()
 
         shaft_end = max([sublist[-1] for sublist in self.nodes_pos])
 
@@ -2351,7 +2060,7 @@ class StaticResults:
                 x=[-0.01 * shaft_end, 1.01 * shaft_end],
                 y=[0, 0],
                 mode="lines",
-                line=dict(width=3.0, color="black", dash="dashdot"),
+                line=dict(color="black", dash="dashdot"),
                 showlegend=False,
                 hoverinfo="none",
             )
@@ -2373,8 +2082,7 @@ class StaticResults:
                     x=xnew,
                     y=ynew,
                     mode="lines",
-                    line=dict(width=5.0, color=colors1[count]),
-                    name="Shaft {}".format(count),
+                    name=f"Shaft {count}",
                     showlegend=True,
                     hovertemplate=(
                         "Shaft lengh: %{x:.2f}<br>" + "Displacement: %{y:.2e}"
@@ -2383,44 +2091,19 @@ class StaticResults:
             )
             count += 1
 
-        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,
-        )
+        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)
 
         return fig
 
-    def plot_free_body_diagram(self, **kwargs):
+    def plot_free_body_diagram(self, fig=None, **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", ...).
@@ -2431,20 +2114,17 @@ 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
-        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))
-            ],
-        )
+        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))
+                ],
+            )
         j = 0
         y_start = 5.0
         for nodes_pos, nodes in zip(self.nodes_pos, self.nodes):
@@ -2456,7 +2136,7 @@ class StaticResults:
                     x=nodes_pos,
                     y=np.zeros(len(nodes_pos)),
                     mode="lines",
-                    line=dict(width=10.0, color="black"),
+                    line=dict(color="black"),
                     hoverinfo="none",
                     showlegend=False,
                 ),
@@ -2468,7 +2148,7 @@ class StaticResults:
                     x=nodes_pos,
                     y=[y_start] * len(nodes_pos),
                     mode="lines",
-                    line=dict(width=3.0, color="black"),
+                    line=dict(color="black"),
                     hoverinfo="none",
                     showlegend=False,
                 ),
@@ -2505,7 +2185,6 @@ class StaticResults:
                 xshift=125,
                 yshift=20,
                 text="<b>Shaft weight = {}N</b>".format(text),
-                font=dict(size=20),
                 align="right",
                 showarrow=False,
             )
@@ -2524,7 +2203,6 @@ 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,
@@ -2549,7 +2227,6 @@ 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,
@@ -2561,17 +2238,7 @@ class StaticResults:
                         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_xaxes(title_text="<b>Shaft Length</b>", row=row, col=col)
             fig.update_yaxes(
                 visible=False, gridcolor="lightgray", showline=False, row=row, col=col
             )
@@ -2581,7 +2248,7 @@ class StaticResults:
 
         return fig
 
-    def plot_shearing_force(self, **kwargs):
+    def plot_shearing_force(self, fig=None, **kwargs):
         """Plot the rotor shearing force diagram.
 
         This method plots:
@@ -2589,6 +2256,8 @@ 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, ...).
@@ -2599,13 +2268,8 @@ class StaticResults:
         fig : Plotly graph_objects.Figure()
             The figure object with the plot.
         """
-        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()
+        if fig is None:
+            fig = go.Figure()
 
         shaft_end = max([sublist[-1] for sublist in self.nodes_pos])
 
@@ -2615,7 +2279,7 @@ class StaticResults:
                 x=[-0.1 * shaft_end, 1.1 * shaft_end],
                 y=[0, 0],
                 mode="lines",
-                line=dict(width=3.0, color="black", dash="dashdot"),
+                line=dict(color="black", dash="dashdot"),
                 showlegend=False,
                 hoverinfo="none",
             )
@@ -2628,9 +2292,8 @@ class StaticResults:
                     x=Vx_axis,
                     y=Vx,
                     mode="lines",
-                    line=dict(width=5.0, color=colors1[j]),
-                    name="Shaft {}".format(j),
-                    legendgroup="Shaft {}".format(j),
+                    name=f"Shaft {j}",
+                    legendgroup=f"Shaft {j}",
                     showlegend=True,
                     hovertemplate=(
                         "Shaft lengh: %{x:.2f}<br>" + "Shearing Force: %{y:.2f}"
@@ -2640,57 +2303,31 @@ class StaticResults:
             j += 1
 
         fig.update_xaxes(
-            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,
+            title_text="<b>Shaft Length</b>", range=[-0.1 * shaft_end, 1.1 * shaft_end]
         )
+        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, **kwargs):
+    def plot_bending_moment(self, fig=None, **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
         """
-        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()
+        if fig is None:
+            fig = go.Figure()
 
         shaft_end = max([sublist[-1] for sublist in self.nodes_pos])
 
@@ -2700,7 +2337,7 @@ class StaticResults:
                 x=[-0.1 * shaft_end, 1.1 * shaft_end],
                 y=[0, 0],
                 mode="lines",
-                line=dict(width=3.0, color="black", dash="dashdot"),
+                line=dict(color="black", dash="dashdot"),
                 showlegend=False,
                 hoverinfo="none",
             )
@@ -2713,9 +2350,8 @@ class StaticResults:
                     x=nodes_pos,
                     y=Bm,
                     mode="lines",
-                    line=dict(width=5.0, color=colors1[j]),
-                    name="Shaft {}".format(j),
-                    legendgroup="Shaft {}".format(j),
+                    name=f"Shaft {j}",
+                    legendgroup=f"Shaft {j}",
                     showlegend=True,
                     hovertemplate=(
                         "Shaft lengh: %{x:.2f}<br>" + "Bending Moment: %{y:.2f}"
@@ -2724,37 +2360,9 @@ class StaticResults:
             )
             j += 1
 
-        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,
-        )
+        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)
 
         return fig
 
@@ -2995,7 +2603,7 @@ class ConvergenceResults:
         self.eigv_arr = eigv_arr
         self.error_arr = error_arr
 
-    def plot(self, **kwargs):
+    def plot(self, fig=None, **kwargs):
         """Plot convergence results.
 
         This method plots:
@@ -3004,6 +2612,8 @@ 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, ...).
@@ -3014,19 +2624,15 @@ class ConvergenceResults:
         fig : Plotly graph_objects.make_subplots()
             The figure object with the plot.
         """
-        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>",
-            ],
-        )
+        if fig is None:
+            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(
@@ -3034,8 +2640,6 @@ 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}"
                 ),
@@ -3044,30 +2648,8 @@ class ConvergenceResults:
             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,
-        )
+        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)
 
         # plot Error vs number of elements
         fig.add_trace(
@@ -3075,8 +2657,6 @@ 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}"
                 ),
@@ -3086,31 +2666,8 @@ class ConvergenceResults:
             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_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_layout(**kwargs)
 
@@ -3162,7 +2719,7 @@ class TimeResponseResults:
         self.nodes_pos = nodes_pos
         self.number_dof = number_dof
 
-    def _plot1d(self, dof, **kwargs):
+    def _plot1d(self, dof, fig=None, **kwargs):
         """Plot time response for a single DoF using Plotly.
 
         This function will take a rotor object and plot its time response using Plotly.
@@ -3171,6 +2728,8 @@ 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, ...).
@@ -3190,20 +2749,14 @@ class TimeResponseResults:
         obs_dof = dof % self.number_dof
         obs_dof = dof_dict[str(obs_dof)]
 
-        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()
+        if fig is None:
+            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,
@@ -3211,38 +2764,18 @@ class TimeResponseResults:
             )
         )
 
-        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_xaxes(title_text="<b>Time</b>")
+        fig.update_yaxes(title_text="<b>Amplitude</b>")
         fig.update_layout(
             title=dict(
-                text="<b>Response for node {} - DoF {}</b>".format(dof // 4, obs_dof),
-                font=dict(size=20),
+                text="<b>Response for node {} - DoF {}</b>".format(dof // 4, obs_dof)
             ),
             **kwargs,
         )
 
         return fig
 
-    def _plot2d(self, node, **kwargs):
+    def _plot2d(self, node, fig=None, **kwargs):
         """Plot orbit response (2D).
 
         This function will take a rotor object and plot its orbit response using Plotly.
@@ -3251,6 +2784,8 @@ 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, ...).
@@ -3261,20 +2796,14 @@ class TimeResponseResults:
         fig : Plotly graph_objects.Figure()
             The figure object with the plot.
         """
-        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()
+        if fig is None:
+            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,
@@ -3284,42 +2813,23 @@ class TimeResponseResults:
             )
         )
 
-        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_xaxes(title_text="<b>Amplitude - X direction</b>")
+        fig.update_yaxes(title_text="<b>Amplitude - Y direction</b>")
         fig.update_layout(
-            title=dict(
-                text="<b>Response for node {}</b>".format(node), font=dict(size=20)
-            ),
-            **kwargs,
+            title=dict(text="<b>Response for node {}</b>".format(node)), **kwargs
         )
 
         return fig
 
-    def _plot3d(self, **kwargs):
+    def _plot3d(self, fig=None, **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", ...).
@@ -3330,11 +2840,8 @@ class TimeResponseResults:
         fig : Plotly graph_objects.Figure()
             The figure object with the plot.
         """
-        kwargs_default_values = dict(hoverlabel_align="right")
-        for k, v in kwargs_default_values.items():
-            kwargs.setdefault(k, v)
-
-        fig = go.Figure()
+        if fig is None:
+            fig = go.Figure()
 
         for n in self.nodes_list:
             x_pos = np.ones(self.yout.shape[0]) * self.nodes_pos[n]
@@ -3344,7 +2851,7 @@ class TimeResponseResults:
                     y=self.yout[:, self.number_dof * n],
                     z=self.yout[:, self.number_dof * n + 1],
                     mode="lines",
-                    line=dict(width=3.0, color="royalblue"),
+                    line=dict(color=tableau_colors["blue"]),
                     name="Mean",
                     legendgroup="mean",
                     showlegend=False,
@@ -3366,47 +2873,23 @@ class TimeResponseResults:
                 y=line,
                 z=line,
                 mode="lines",
-                line=dict(width=2.0, color="black", dash="dashdot"),
+                line=dict(color="black", dash="dashdot"),
                 showlegend=False,
             )
         )
 
         fig.update_layout(
-            width=1200,
-            height=900,
             scene=dict(
-                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,
-                ),
+                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>")),
             ),
             **kwargs,
         )
 
         return fig
 
-    def plot(self, plot_type="3d", dof=None, node=None, **kwargs):
+    def plot(self, plot_type="3d", dof=None, node=None, fig=None, **kwargs):
         """Plot time response.
 
         The plot type options are:
@@ -3430,6 +2913,8 @@ 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 d66f59b7..0b8b2ad9 100644
--- a/ross/rotor_assembly.py
+++ b/ross/rotor_assembly.py
@@ -18,6 +18,7 @@ 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,
@@ -117,11 +118,7 @@ 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"
 
@@ -589,7 +586,7 @@ class Rotor(object):
         array([91.79655318, 96.28899977])
         >>> modal.wd[:2]
         array([91.79655318, 96.28899977])
-        >>> fig = modal.plot_mode3D(0)
+        >>> fig = modal.plot_mode_3d(0)
         """
         evalues, evectors = self._eigen(speed, num_modes=num_modes, sparse=sparse)
         wn_len = num_modes // 2
@@ -1773,11 +1770,10 @@ 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=30,
+                    size=20,
                     color="#ffcc99",
                     line=dict(width=1.0, color="black"),
                 ),
@@ -1815,31 +1811,17 @@ 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(
-            plot_bgcolor="white",
-            title=dict(text="<b>Rotor Model</b>", font=dict(size=18)),
-            **kwargs,
-        )
+        fig.update_layout(title=dict(text="<b>Rotor Model</b>"), **kwargs)
 
         return fig
 
@@ -1905,7 +1887,7 @@ class Rotor(object):
 
         return results
 
-    def plot_ucs(
+    def _calc_ucs(
         self,
         stiffness_range=None,
         num_modes=16,
@@ -1914,6 +1896,92 @@ class Rotor(object):
         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,
+        fig=None,
+        synchronous=False,
+        **kwargs,
+    ):
         """Plot undamped critical speed map.
 
         This method will plot the undamped critical speed map for a given range
@@ -1972,66 +2040,49 @@ 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 = 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]
+        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,
+        )
 
         if fig is None:
             fig = go.Figure()
 
+        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,
+                )
+            )
+
+        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="",
+            )
+        )
+
         fig.add_trace(
             go.Scatter(
                 x=bearing0.kxx.interpolated(bearing0.frequency),
                 y=bearing0.frequency,
-                mode="markers",
-                marker=dict(size=10, symbol="circle", color="#888844"),
+                mode="lines",
+                line=dict(dash="dashdot"),
+                hoverinfo="none",
                 name="Kxx",
             )
         )
@@ -2039,57 +2090,25 @@ class Rotor(object):
             go.Scatter(
                 x=bearing0.kyy.interpolated(bearing0.frequency),
                 y=bearing0.frequency,
-                mode="markers",
-                marker=dict(size=10, symbol="square", color="#888844"),
+                mode="lines",
+                line=dict(dash="dashdot"),
+                hoverinfo="none",
                 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</b>",
-            title_font=dict(size=16),
-            tickfont=dict(size=14),
-            gridcolor="lightgray",
-            showline=True,
-            linewidth=2.5,
-            linecolor="black",
-            mirror=True,
+            title_text="<b>Bearing Stiffness (N/m)</b>",
             type="log",
             exponentformat="power",
         )
         fig.update_yaxes(
-            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,
+            title_text="<b>Critical Speed (rad/s)</b>",
             type="log",
             exponentformat="power",
         )
         fig.update_layout(
-            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,
+            title=dict(text="<b>Undamped Critical Speed Map</b>"), **kwargs
         )
 
         return fig
@@ -2107,6 +2126,8 @@ 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, ...).
@@ -2974,11 +2995,7 @@ 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"
 


Test patch:
diff --git a/ross/tests/test_rotor_assembly.py b/ross/tests/test_rotor_assembly.py
index af851263..9d332b68 100644
--- a/ross/tests/test_rotor_assembly.py
+++ b/ross/tests/test_rotor_assembly.py
@@ -1720,3 +1720,52 @@ 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)