Repository: ross-rotordynamics/ross
Pull request #270 (2019-09-03 17:50:59)

Add new "Tag" attribute to elements
Disk, bearing and seal elements now get an attribute called "Tag".
This allows the user to name the elements to help refer to the actual equipment (impellers, blades, labyrinth seal...)

If the user do not input a tag, it's named after:
- Disk 1, 2, 3... 
- Bearing 1, 2, 3...
- Seal 1, 2, 3...

As bearings and seals are input in the same list of objects, this situation below may happens:
```

def rotor_example():
    i_d = 0
    o_d = 0.05
    n = 6
    L = [0.25 for _ in range(n)]

    shaft_elem = [
        ShaftElement(
            l, i_d, o_d, steel, shear_effects=True, rotary_inertia=True, gyroscopic=True
        )
        for l in L
    ]

    disk0 = DiskElement.from_geometry(
        n=1, material=steel, width=0.07, i_d=0.05, o_d=0.28
    )
    disk1 = DiskElement.from_geometry(
        n=3, material=steel, width=0.07, i_d=0.05, o_d=0.28, tag="disk_test"
    )
    disk2 = DiskElement.from_geometry(
        n=5, material=steel, width=0.07, i_d=0.05, o_d=0.28, tag="disk_test2"
    )

    stfx = 1e8
    stfy = 1e8
    bearing0 = BearingElement(0, kxx=stfx, kyy=stfy, cxx=1000, tag="brg")
    bearing1 = BearingElement(6, kxx=stfx, kyy=stfy, cxx=1000)
    seal0 = SealElement(2, kxx=1e2, kyy=1e2, cxx=500, cyy=500, tag="seal")
    seal1 = SealElement(4, kxx=1e2, kyy=1e2, cxx=500, cyy=500)

    return Rotor(shaft_elem, [disk0, disk1, disk2], [bearing0, seal0, seal1, bearing1])

>>> rotor = rotor_example()
>>> rotor.df_bearings["tag"]
0    brg
1    seal
2    Seal 2
3    Bearing 3
Name: tag, dtype: object
```
The counting for bearings and seals is not independent. I don't see it as big deal, but I'm open to suggestions.



Patch:
diff --git a/ross/api_report b/ross/api_report
deleted file mode 100644
index 1f3026c5..00000000
--- a/ross/api_report
+++ /dev/null
@@ -1,527 +0,0 @@
-import numpy as np
-from scipy.signal import argrelextrema
-
-from rotor_assembly import Rotor
-import ross as rs
-
-import bokeh.palettes as bp
-from bokeh.plotting import figure
-from bokeh.models import (
-        ColumnDataSource,
-        HoverTool,
-        Span,
-        Label,
-)
-
-# set bokeh palette of colors
-bokeh_colors = bp.RdGy[11]
-
-
-class Report:
-    def __init__(
-            self,
-            rotor,
-            minspeed,
-            maxspeed,
-            speed_units="rpm",
-    ):
-        """Report according to standard analysis.
-
-        - Perform Stability_level1 analysis
-        - Apply Level 1 Screening Criteria
-        - Perform Stability_level2 analysis
-
-        Parameters
-        ----------
-        rotor : object
-            A rotor built from rotor_assembly.
-        maxspeed : float
-            Maximum operation speed.
-        minspeed : float
-            Minimum operation speed.
-        speed_units : str
-            String defining the unit for rotor speed.
-            Default is "rpm".
-            
-        Return
-        ------
-
-        Example
-        -------
-        >>> rotor = rotor_example()
-        >>> report = Report(rotor=rotor,
-                            minspeed=400,
-                            maxspeed=1000) # doctest: +ELLIPSIS
-        [...]
-        """
-        self.rotor = rotor
-        self.maxspeed = maxspeed
-        self.minspeed = minspeed
-        self.speed_units = speed_units
-
-        if speed_units == "rpm":
-            minspeed = minspeed * np.pi / 30
-            maxspeed = maxspeed * np.pi / 30
-
-    @classmethod
-    def from_saved_rotors(cls, path):
-        rotor = rs.Rotor.load(path)
-        return cls(rotor)
-
-    def static_forces(self):
-        """
-        Method to calculate the bearing reaction forces.
-
-        Returns
-        -------
-        Fb : list
-            Bearing reaction forces.
-
-        Example
-        -------
-        >>> rotor = rotor_example()
-        >>> report = Report(rotor=rotor, minspeed=400, maxspeed=1000)
-        >>> report.static_forces() # doctest: +ELLIPSIS
-        [...]
-        """
-        # get reaction forces on bearings
-        Fb = self.rotor.run_static().force_data['Bearings Reaction Forces']
-        Fb = np.array(Fb) / 9.8065
-
-        return Fb
-
-    def unbalance_forces(self, mode):
-        """
-        Method to calculate the unbalance forces.
-
-        Parameters
-        ----------
-        mode : int
-            n'th mode shape.
-
-        Returns
-        -------
-        U : list
-            Unbalancing forces.
-
-        Example
-        -------
-        >>> rotor = rotor_example()
-        >>> report = Report(rotor=rotor, minspeed=400, maxspeed=1000)
-        >>> report.unbalance_forces(mode=0) # doctest: +ELLIPSIS
-        [...]
-        """
-
-        N = 60 * self.maxspeed / (2 * np.pi)
-
-        # get reaction forces on bearings
-        Fb = self.static_forces()
-
-        if mode == 0 or mode == 1:
-            U = [4 * max(6350 * np.sum(Fb) / N, 250e-3 * np.sum(Fb))]
-
-        if mode == 2 or mode == 3:
-            U = 4 * max(6350 * Fb / N, 250e-3 * Fb)
-
-        # get disk masses
-        W3 = 0
-        for disk in self.rotor.disk_elements:
-            if disk.n > self.rotor.df_bearings['n'].iloc[-1]:
-                W3 += disk.m
-                U = [6350 * W3 / N]
-
-        return U
-
-    def unbalance_response(
-            self,
-            clearances,
-            mode,
-    ):
-        """Evaluates the unbalance response for the rotor.
-
-        This analysis takes the critical speeds of interest, calculates the
-        position and weight of the required unbalance and performs the analysis
-        including:
-         - Check if vibration at MCS is below the limit with the applied weight;
-         - Check if the clearances are ok if the vibration deteriorate to the
-         limit level;
-
-        Parameters
-        ----------
-        clearances : dict
-            Dict mapping between node and its clearance in meters.
-            (e.g. clearances = dict(3=0.001, 5=0.002...)
-        mode : int
-            n'th mode shape.
-
-        Returns
-        -------
-        mag_plot : bokeh axes
-            Bokeh axes with unbalance response plot.
-
-        Example
-        -------
-        >>> rotor = rotor_example()
-        >>> report = Report(rotor=rotor)
-        >>> clearances = dict(3=0.001, 5=0.002)
-        >>> report.unbalance_response(clearances=clearances,
-                                      mode=0,
-                                      ) # doctest: +ELLIPSIS
-        [...]
-        """
-        maxspeed = self.maxspeed
-        minspeed = self.minspeed
-        freq_range = np.linspace(0, 1.25 * maxspeed, 201)
-
-        # returns de nodes where forces will be applied
-        node_min, node_max = self.mode_shape(mode)
-
-        nodes = [int(node) for sub_nodes in [node_min, node_max] for node in sub_nodes]
-
-        force = self.unbalance_forces(mode)
-
-        phase = []
-        for node in nodes:
-            phase.append(np.pi)
-
-        response = self.rotor.unbalance_response(nodes, force, phase, freq_range)
-        mag = response.magnitude
-
-        for node in nodes:
-            dof = 4 * node + 1
-            mag_plot = response.plot_magnitude_bokeh(dof)
-
-        magnitude = mag[dof]
-        idx_max = argrelextrema(magnitude, np.greater)[0].tolist()
-        wn = freq_range[idx_max]
-
-        AF_table = []
-        SM_table = []
-        SM_ref_table = []
-
-        for i, peak in enumerate(magnitude[idx_max]):
-            peak_n = 0.707 * peak
-            peak_aux = np.linspace(peak_n, peak_n, len(freq_range))
-
-            idx = np.argwhere(
-                np.diff(np.sign(peak_aux - magnitude))
-            ).flatten()
-            idx = np.sort(np.append(idx, idx_max[i]))
-
-            # if speed range is not long enough to catch the magnitudes
-            try:
-                idx_aux = [
-                    list(idx).index(idx_max[i]) - 1,
-                    list(idx).index(idx_max[i]) + 1,
-                ]
-                idx = idx[idx_aux]
-            except IndexError:
-                idx = [
-                    list(idx).index(idx_max[i]) - 1,
-                    len(freq_range) - 1
-                ]
-
-            # Amplification Factor (AF) - API684 - SP6.8.2.1
-            AF = wn[i] / (
-                freq_range[idx[1]] - freq_range[idx[0]]
-            )
-
-            # Separation Margin (SM) - API684 - SP6.8.2.10
-            if AF > 2.5 and wn[i] < minspeed:
-                SM = min([16, 17 * (1 - 1 / (AF - 1.5))]) / 100
-                SMspeed = wn[i] * (1 + SM)
-                SM_ref = (minspeed - wn[i]) / wn[i]
-                source = ColumnDataSource(
-                    dict(
-                        top=[max(magnitude[idx_max])],
-                        bottom=[0],
-                        left=[wn[i]],
-                        right=[SMspeed],
-                        tag1=[wn[i]],
-                        tag2=[SMspeed],
-                    )
-                )
-
-                mag_plot.quad(
-                    top="top",
-                    bottom="bottom",
-                    left="left",
-                    right="right",
-                    source=source,
-                    line_color=bokeh_colors[8],
-                    line_width=0.8,
-                    fill_alpha=0.2,
-                    fill_color=bokeh_colors[8],
-                    legend="Separation Margin",
-                    name="SM2",
-                )
-                hover = HoverTool(names=["SM2"])
-                hover.tooltips = [
-                    ("Critical Speed :", "@tag1"),
-                    ("Speed at 0.707 x peak amplitude :", "@tag2"),
-                ]
-                mag_plot.add_tools(hover)
-
-            elif AF > 2.5 and wn[i] > maxspeed:
-                SM = min([26, 10 + 17 * (1 - 1 / (AF - 1.5))]) / 100
-                SMspeed = wn[i] * (1 - SM)
-                SM_ref = (wn[i] - maxspeed) / maxspeed
-                source = ColumnDataSource(
-                    dict(
-                        top=[max(magnitude[idx_max])],
-                        bottom=[0],
-                        left=[SMspeed],
-                        right=[wn[i]],
-                        tag1=[wn[i]],
-                        tag2=[SMspeed],
-                    )
-                )
-
-                mag_plot.quad(
-                    top="top",
-                    bottom="bottom",
-                    left="left",
-                    right="right",
-                    source=source,
-                    line_color=bokeh_colors[8],
-                    line_width=0.8,
-                    fill_alpha=0.2,
-                    fill_color=bokeh_colors[8],
-                    legend="Separation Margin",
-                    name="SM2",
-                )
-                hover = HoverTool(names=["SM2"])
-                hover.tooltips = [
-                    ("Critical Speed :", "@tag1"),
-                    ("Speed at 0.707 x peak amplitude :", "@tag2"),
-                ]
-                mag_plot.add_tools(hover)
-
-            else:
-                SM = None
-                SM_ref = None
-                SMspeed = None
-
-        # amplitude limit in micrometers (A1) - API684 - SP6.8.2.11
-        A1 = 25.4 * np.sqrt(12000 / (30 * maxspeed / np.pi))
-
-        # amplitude from mode shape analysis
-        Amax = max(mag[dof])
-
-        # Scale Factor (Scc) - API684 - SP6.8.2.11 / API617 - 4.8.2.11
-        Scc = max(A1 / Amax, 0.5)
-        Scc = min(Scc, 6.0)
-
-        mag_plot.quad(
-            top=max(mag[dof]),
-            bottom=0,
-            left=minspeed,
-            right=maxspeed,
-            line_color="green",
-            line_width=0.8,
-            fill_alpha=0.2,
-            fill_color="green",
-            legend="Operation Speed Range",
-        )
-
-        source = ColumnDataSource(dict(x=freq_range, y=mag[dof]))
-        mag_plot.line(
-            x="x",
-            y="y",
-            source=source,
-            line_color=bokeh_colors[0],
-            line_alpha=1.0,
-            line_width=3,
-        )
-        mag_plot.line(
-            x=[minspeed, maxspeed],
-            y=[A1, A1],
-            line_dash="dotdash",
-            line_width=2.0,
-            line_color=bokeh_colors[1],
-            legend="Av1 - Mechanical test vibration limit",
-        )
-        mag_plot.add_layout(
-            Label(
-                x=(minspeed+maxspeed)/2,
-                y=A1,
-                angle=0,
-                text="Av1",
-                text_font_style="bold",
-                text_font_size="12pt",
-                text_baseline="top",
-                text_align="center",
-                y_offset=20,
-            )
-        )
-        mag_plot.width = 1280
-        mag_plot.height = 720
-        mag_plot.title.text_font_size = "14pt"
-
-        return mag_plot
-
-    def mode_shape(self, mode):
-        """Evaluates the mode shapes for the rotor.
-
-        This analysis presents the vibration mode for each critical speed.
-        The importance is to locate the critical node, where the displacement
-        is the greatest, then apply loads for unbalance response (stability
-        level 1)
-
-        Parameters
-        ----------
-        mode : int
-            the n'th vibration mode
-
-        Returns
-        -------
-        node_min, node_max : list
-            List with nodes where the largest absolute displacements occur
-
-        Example
-        -------
-        >>> rotor = rotor_example()
-        >>> report = Report(rotor=rotor)
-        >>> report.mode_shpae(mode=0) # doctest: +ELLIPSIS
-        [...]
-        """
-        nodes_pos = self.rotor.nodes_pos
-
-        xn, yn, zn, xc, yc, zc_pos, nn = self.rotor.run_mode_shapes().calc_mode_shape(mode=mode)
-
-        # reduce 3D view to 2D view
-        vn = np.zeros(len(zn))
-        for i in range(len(zn)):
-            theta = np.arctan(xn[i] / yn[i])
-            vn[i] = xn[i] * np.sin(theta) + yn[i] * np.cos(theta)
-
-        # remove repetitive values from zn and vn
-        idx_remove = []
-        for i in range(1, len(zn)):
-            if zn[i] == zn[i-1]:
-                idx_remove.append(i)
-        zn = np.delete(zn, idx_remove)
-        vn = np.delete(vn, idx_remove)
-
-        aux_idx_max = argrelextrema(vn, np.greater)[0].tolist()
-        aux_idx_min = argrelextrema(vn, np.less)[0].tolist()
-
-        node_min = []
-        node_max = []
-
-        # verification of rigid modes
-        if len(aux_idx_max) == 0 and len(aux_idx_min) == 0:
-            idx_max = np.argmax(vn)
-            idx_min = np.argmin(vn)
-
-            # corrects the index by the removed points
-            for i in idx_remove:
-                if idx_min > i:
-                    idx_min += 1
-                if idx_max > i:
-                    idx_max += 1
-            node_max = np.round(np.array(idx_max) / nn)
-            node_min = np.round(np.array(idx_min) / nn)
-
-        if len(aux_idx_min) != 0:
-            idx_min = np.where(vn == min(vn[aux_idx_min]))[0].tolist()
-
-            # corrects the index by the removed points
-            for i in idx_remove:
-                if idx_min[0] > i:
-                    idx_min[0] += 1
-            node_min = np.round(np.array(idx_min) / nn)
-
-        if len(aux_idx_max) != 0:
-            idx_max = np.where(vn == max(vn[aux_idx_max]))[0].tolist()
-
-            # corrects the index by the removed points
-            for i in idx_remove:
-                if idx_max[0] > i:
-                    idx_max[0] += 1
-            node_max = np.round(np.array(idx_max) / nn)
-
-        TOOLS = "pan,wheel_zoom,box_zoom,reset,save,box_select"
-        plot = figure(
-            tools=TOOLS,
-            width=1400,
-            height=700,
-            title="Static Analysis",
-            x_axis_label="shaft lenght",
-            y_axis_label="lateral displacement",
-        )
-
-        nodes_pos = np.array(nodes_pos)
-
-        plot.line(zn, vn, line_width=4, line_color="red")
-        plot.line(
-            x=nodes_pos,
-            y=np.zeros(len(nodes_pos)),
-            line_dash="dotdash",
-            line_width=4.0,
-            line_color='black'
-        )
-        plot.circle(
-            x=nodes_pos[self.rotor.df_bearings['n'].tolist()],
-            y=np.zeros(len(self.rotor.df_bearings)),
-            size=12,
-            fill_color='black'
-        )
-        plot.add_layout(
-            Label(
-                x=np.mean(nodes_pos[self.rotor.df_bearings['n'].tolist()]),
-                y=0,
-                angle=0,
-                text="Bearing Span",
-                text_font_style="bold",
-                text_font_size="12pt",
-                text_baseline="top",
-                text_align="center",
-                y_offset=20,
-            )
-        )
-        for node in nodes_pos[self.rotor.df_bearings['n'].tolist()]:
-            plot.add_layout(
-                Span(
-                    location=node,
-                    dimension='height',
-                    line_color='green',
-                    line_dash='dashed',
-                    line_width=3,
-                )
-            )
-
-        return node_min, node_max
-
-    def stability_level_1(self):
-        """Stability analysis level 1.
-
-        This analysis consider a anticipated cross coupling QA based on
-        conditions at the normal operating point and the cross-coupling
-        required to produce a zero log decrement, Q0.
-
-        Components such as seals and impellers are not considered in this
-        analysis.
-
-        Parameters
-        ----------
-        (Check what we need to calculate the applied cross coupling and list
-         them as parameters)
-        """
-        pass
-
-    def stability_level_2(self):
-        """Stability analysis level 2.
-
-        For the level 2 stability analysis additional sources that contribute
-        to the rotor stability shall be considered such as:
-        a)  labyrinth seals;
-        b)  damper seals;
-        c)  impeller/blade flow aerodynamic effects;
-        d)  internal friction.
-
-        Parameters
-        ----------
-        (Check what we need to calculate the applied cross coupling and list
-        them as parameters)
-        """
diff --git a/ross/bearing_seal_element.py b/ross/bearing_seal_element.py
index 6b6bb666..61fdbeec 100644
--- a/ross/bearing_seal_element.py
+++ b/ross/bearing_seal_element.py
@@ -36,20 +36,14 @@ class _Coefficient:
             #  dfitpack.error is not exposed by scipy
             #  so a bare except is used
             except:
-                try:
-                    if len(self.w) in (2, 3):
-                        self.interpolated = interpolate.interp1d(
-                            self.w, self.coefficient, kind=len(self.w) - 1
-                        )
-                except:
-                    raise ValueError(
-                        "Arguments (coefficients and w)" " must have the same dimension"
-                    )
+                raise ValueError(
+                    "Arguments (coefficients and w)" " must have the same dimension"
+                )
         else:
             self.interpolated = lambda x: np.array(self.coefficient[0])
 
     def __eq__(self, other):
-        if np.allclose(self.__dict__["coefficient"], other.__dict__["coefficient"]):
+        if np.allclose(self.__dict__["coefficient"],other.__dict__["coefficient"]):
             return True
         else:
             return False
@@ -125,6 +119,9 @@ class BearingElement(Element):
         (defaults to 0)
     w: array, optional
         Array with the speeds (rad/s).
+     tag : str, optional
+         A tag to name the element
+         Default is None
     Examples
     --------
     >>> # A bearing element located in the first rotor node, with these
@@ -144,7 +141,7 @@ class BearingElement(Element):
     """
 
     def __init__(
-        self, n, kxx, cxx, kyy=None, kxy=0, kyx=0, cyy=None, cxy=0, cyx=0, w=None
+        self, n, kxx, cxx, kyy=None, kxy=0, kyx=0, cyy=None, cxy=0, cyx=0, w=None, tag=None
     ):
 
         args = ["kxx", "kyy", "kxy", "kyx", "cxx", "cyy", "cxy", "cyx"]
@@ -189,6 +186,7 @@ class BearingElement(Element):
         self.n_r = n
 
         self.w = np.array(w, dtype=np.float64)
+        self.tag = tag
         self.color = "#355d7a"
 
     def __repr__(self):
@@ -199,7 +197,7 @@ class BearingElement(Element):
             f" kyx={self.kyx}, kyy={self.kyy},\n"
             f" cxx={self.cxx}, cxy={self.cxy},\n"
             f" cyx={self.cyx}, cyy={self.cyy},\n"
-            f" w={self.w})"
+            f" w={self.w}, tag={self.tag!r})"
         )
 
     def __eq__(self, other):
@@ -243,6 +241,7 @@ class BearingElement(Element):
             "cxy": self.cxy.coefficient,
             "cyx": self.cyx.coefficient,
             "w": w,
+            "tag": self.tag,
         }
         self.dump_data(data, file_name)
 
@@ -338,19 +337,19 @@ class BearingElement(Element):
         Returns
         -------
         """
-        default_values = dict(lw=1.0, alpha=1.0, c="k")
+        default_values = dict(lw=1., alpha=1., c='k')
         for k, v in default_values.items():
             kwargs.setdefault(k, v)
 
         # geometric factors
         zpos, ypos = position
-        coils = 6  # number of points to generate spring
-        step = mean / 5  # spring step
-        a = 2.2  # range(1.8 to 2.8)
-        b = 1.4  # range(1.2 to 1.5)
-        c = mean / 2.9  # defines width
-        d = (coils / a) * step  # modifies damper width
-        n = 5  # number of ground lines
+        coils = 6               # number of points to generate spring
+        step = mean / 5         # spring step
+        a = 2.2                 # range(1.8 to 2.8)
+        b = 1.4                 # range(1.2 to 1.5)
+        c = mean / 2.9          # defines width
+        d = (coils / a)*step    # modifies damper width
+        n = 5                   # number of ground lines
 
         zs0 = zpos - c
         zs1 = zpos + c
@@ -369,7 +368,7 @@ class BearingElement(Element):
             b * ys0 + (coils + 1) * step,
             ys0 + (coils + 1) * step,
             ys0 + (coils + 1) * step,
-            ys0 + (coils + 1) * step,
+            ys0 + (coils + 1) * step
         ]
         yu_top = [-y for y in yl_top]
 
@@ -385,7 +384,7 @@ class BearingElement(Element):
 
         step2 = (zl_g[1] - zl_g[0]) / n
         for i in range(n + 1):
-            zl_g2 = [(zs0 - step) + step2 * (i), (zs0 - step) + step2 * (i + 1)]
+            zl_g2 = [(zs0 - step) + step2*(i), (zs0 - step) + step2*(i+1)]
             yl_g2 = [yl_g[0], 1.1 * yl_g[0]]
             yu_g2 = [-y for y in yl_g2]
             ax.add_line(mlines.Line2D(zl_g2, yl_g2, **kwargs))
@@ -394,7 +393,11 @@ class BearingElement(Element):
         # plot spring
         z_spring = np.array([zs0, zs0, zs0, zs0])
         yl_spring = np.array(
-            [ys0, ys0 + step, ys0 + coils * step, ys0 + (coils + 1) * step]
+            [ys0,
+             ys0 + step,
+             ys0 + coils * step,
+             ys0 + (coils + 1) * step
+             ]
         )
 
         for i in range(coils):
@@ -425,7 +428,7 @@ class BearingElement(Element):
             ys0 + 1.5 * d,
             ys0 + 1.5 * d,
             ys0 + 1.5 * d,
-            ys0 + (coils + 1) * step,
+            ys0 + (coils + 1) * step
         ]
         yu_damper3 = [-y for y in yl_damper3]
 
@@ -454,13 +457,13 @@ class BearingElement(Element):
 
         # geometric factors
         zpos, ypos = position
-        coils = 6  # number of points to generate spring
-        step = mean / 5  # spring step
-        a = 2.2  # range(1.8 to 2.8)
-        b = 1.4  # range(1.2 to 1.5)
-        c = mean / 2.9  # defines width
-        d = (coils / a) * step  # modifies damper width
-        n = 5  # number of ground lines
+        coils = 6               # number of points to generate spring
+        step = mean / 5         # spring step
+        a = 2.2                 # range(1.8 to 2.8)
+        b = 1.4                 # range(1.2 to 1.5)
+        c = mean / 2.9          # defines width
+        d = (coils / a)*step    # modifies damper width
+        n = 5                   # number of ground lines
 
         zs0 = zpos - c
         zs1 = zpos + c
@@ -479,7 +482,7 @@ class BearingElement(Element):
             b * ys0 + (coils + 1) * step,
             ys0 + (coils + 1) * step,
             ys0 + (coils + 1) * step,
-            ys0 + (coils + 1) * step,
+            ys0 + (coils + 1) * step
         ]
         yu_top = [-y for y in yl_top]
 
@@ -494,8 +497,8 @@ class BearingElement(Element):
         bk_ax.line(x=zl_g, y=yu_g, **kwargs)
 
         step2 = (zl_g[1] - zl_g[0]) / n
-        for i in range(n + 1):
-            zl_g2 = [(zs0 - step) + step2 * (i), (zs0 - step) + step2 * (i + 1)]
+        for i in range(n+1):
+            zl_g2 = [(zs0 - step) + step2*(i), (zs0 - step) + step2*(i+1)]
             yl_g2 = [yl_g[0], 1.1 * yl_g[0]]
             yu_g2 = [-y for y in yl_g2]
             bk_ax.line(x=zl_g2, y=yl_g2, **kwargs)
@@ -504,7 +507,11 @@ class BearingElement(Element):
         # plot spring
         z_spring = np.array([zs0, zs0, zs0, zs0])
         yl_spring = np.array(
-            [ys0, ys0 + step, ys0 + coils * step, ys0 + (coils + 1) * step]
+            [ys0,
+             ys0 + step,
+             ys0 + coils * step,
+             ys0 + (coils + 1) * step
+             ]
         )
 
         for i in range(coils):
@@ -535,7 +542,7 @@ class BearingElement(Element):
             ys0 + 1.5 * d,
             ys0 + 1.5 * d,
             ys0 + 1.5 * d,
-            ys0 + (coils + 1) * step,
+            ys0 + (coils + 1) * step
         ]
         yu_damper3 = [-y for y in yl_damper3]
 
@@ -591,18 +598,18 @@ class BearingElement(Element):
         -------
         A bearing object.
         """
-        parameters = read_table_file(file, "bearing", sheet_name, n)
+        parameters = read_table_file(file, 'bearing', sheet_name, n)
         return cls(
-            n=parameters["n"],
-            kxx=parameters["kxx"],
-            cxx=parameters["cxx"],
-            kyy=parameters["kyy"],
-            kxy=parameters["kxy"],
-            kyx=parameters["kyx"],
-            cyy=parameters["cyy"],
-            cxy=parameters["cxy"],
-            cyx=parameters["cyx"],
-            w=parameters["w"],
+            n=parameters['n'],
+            kxx=parameters['kxx'],
+            cxx=parameters['cxx'],
+            kyy=parameters['kyy'],
+            kxy=parameters['kxy'],
+            kyx=parameters['kyx'],
+            cyy=parameters['cyy'],
+            cxy=parameters['cxy'],
+            cyx=parameters['cyx'],
+            w=parameters['w'],
         )
 
     @classmethod
@@ -719,6 +726,7 @@ class SealElement(BearingElement):
         cyx=0,
         w=None,
         seal_leakage=None,
+        tag=None,
     ):
         super().__init__(
             n=n,
@@ -731,6 +739,7 @@ class SealElement(BearingElement):
             cxy=cxy,
             cyx=cyx,
             cyy=cyy,
+            tag=tag,
         )
 
         self.seal_leakage = seal_leakage
@@ -755,22 +764,26 @@ class SealElement(BearingElement):
 
         #  node (x pos), outer diam. (y pos)
         bearing_points_u = [
-            [zpos - step, 1.1 * ypos],  # upper
-            [zpos - step, 1.1 * ypos + 2 * step],
-            [zpos + step, 1.1 * ypos + 2 * step],
-            [zpos + step, 1.1 * ypos],
-            [zpos - step, 1.1 * ypos],
+            [zpos - step, 1.1*ypos],  # upper
+            [zpos - step, 1.1*ypos + 2*step],
+            [zpos + step, 1.1*ypos + 2*step],
+            [zpos + step, 1.1*ypos],
+            [zpos - step, 1.1*ypos],
         ]
         bearing_points_l = [
-            [zpos - step, -1.1 * ypos],  # upper
-            [zpos - step, -1.1 * ypos - 2 * step],
-            [zpos + step, -1.1 * ypos - 2 * step],
-            [zpos + step, -1.1 * ypos],
-            [zpos - step, -1.1 * ypos],
+            [zpos - step, -1.1*ypos],  # upper
+            [zpos - step, -1.1*ypos - 2*step],
+            [zpos + step, -1.1*ypos - 2*step],
+            [zpos + step, -1.1*ypos],
+            [zpos - step, -1.1*ypos],
         ]
 
-        ax.add_patch(mpatches.Polygon(bearing_points_u, color=self.color, picker=True))
-        ax.add_patch(mpatches.Polygon(bearing_points_l, color=self.color, picker=True))
+        ax.add_patch(
+            mpatches.Polygon(bearing_points_u, color=self.color, picker=True)
+        )
+        ax.add_patch(
+            mpatches.Polygon(bearing_points_l, color=self.color, picker=True)
+        )
 
     def bokeh_patch(self, position, length, bk_ax):
         """Seal element patch.
@@ -791,8 +804,8 @@ class SealElement(BearingElement):
 
         # bokeh plot - upper seal visual representation
         bk_ax.quad(
-            top=1.1 * ypos + 2 * step,
-            bottom=1.1 * ypos,
+            top=1.1*ypos + 2*step,
+            bottom=1.1*ypos,
             left=zpos - step,
             right=zpos + step,
             line_color=bokeh_colors[6],
@@ -803,8 +816,8 @@ class SealElement(BearingElement):
         )
         # bokeh plot - lower seal visual representation
         bk_ax.quad(
-            top=-1.1 * ypos,
-            bottom=-1.1 * ypos - 2 * step,
+            top=-1.1*ypos,
+            bottom=-1.1*ypos - 2*step,
             left=zpos - step,
             right=zpos + step,
             line_color=bokeh_colors[6],
diff --git a/ross/disk_element.py b/ross/disk_element.py
index 671797d7..c8f05c6f 100644
--- a/ross/disk_element.py
+++ b/ross/disk_element.py
@@ -24,6 +24,9 @@ class DiskElement(Element):
          Diametral moment of inertia.
      Ip : float
          Polar moment of inertia
+     tag : str, optional
+         A tag to name the element
+         Default is None
      References
      ----------
      .. [1] 'Dynamics of Rotating Machinery' by MI Friswell, JET Penny, SD Garvey
@@ -35,7 +38,7 @@ class DiskElement(Element):
      0.32956362
      """
 
-    def __init__(self, n, m, Id, Ip):
+    def __init__(self, n, m, Id, Ip, tag=None):
         self.n = int(n)
         self.n_l = n
         self.n_r = n
@@ -43,6 +46,7 @@ class DiskElement(Element):
         self.m = m
         self.Id = Id
         self.Ip = Ip
+        self.tag = tag
         self.color = "#bc625b"
 
     def __eq__(self, other):
@@ -70,7 +74,7 @@ class DiskElement(Element):
             f"{self.__class__.__name__}"
             f"(Id={self.Id:{0}.{5}}, Ip={self.Ip:{0}.{5}}, "
             f"m={self.m:{0}.{5}}, color={self.color!r}, "
-            f"n={self.n})"
+            f"n={self.n}, tag={self.tag!r})"
         )
 
     def save(self, file_name):
@@ -80,6 +84,7 @@ class DiskElement(Element):
             "m": self.m,
             "Id": self.Id,
             "Ip": self.Ip,
+            "tag": self.tag,
         }
         self.dump_data(data, file_name)
 
@@ -240,6 +245,7 @@ class DiskElement(Element):
                 IP=[self.Ip],
                 ID=[self.Id],
                 mass=[self.m],
+                tag=[self.tag],
             )
         )
         source_c = ColumnDataSource(
@@ -252,6 +258,7 @@ class DiskElement(Element):
                 IP=[self.Ip],
                 ID=[self.Id],
                 mass=[self.m],
+                tag=[self.tag],
             )
         )
 
@@ -299,13 +306,15 @@ class DiskElement(Element):
             ("Polar Moment of Inertia :", "@IP"),
             ("Diametral Moment of Inertia :", "@ID"),
             ("Disk mass :", "@mass"),
+            ("Tag :", "@tag")
         ]
         hover.mode = "mouse"
 
         return hover
 
+
     @classmethod
-    def from_geometry(cls, n, material, width, i_d, o_d):
+    def from_geometry(cls, n, material, width, i_d, o_d, tag=None):
         """A disk element.
         This class will create a disk element from input data of geometry.
         Parameters
@@ -328,6 +337,9 @@ class DiskElement(Element):
             Diametral moment of inertia.
         Ip : float
             Polar moment of inertia
+        tag : str, optional
+            A tag to name the element
+            Default is None
         References
         ----------
         .. [1] 'Dynamics of Rotating Machinery' by MI Friswell, JET Penny, SD Garvey
@@ -348,7 +360,9 @@ class DiskElement(Element):
         # fmt: on
         Ip = 0.03125 * material.rho * np.pi * width * (o_d ** 4 - i_d ** 4)
 
-        return cls(n, m, Id, Ip)
+        tag = tag
+
+        return cls(n, m, Id, Ip, tag)
 
     @classmethod
     def from_table(cls, file, sheet_name=0):
diff --git a/ross/results.py b/ross/results.py
index 838a2a07..6f8c76c7 100644
--- a/ross/results.py
+++ b/ross/results.py
@@ -1,15 +1,14 @@
 import pickle
 import numpy as np
 from scipy import interpolate
-from scipy.signal import argrelextrema
+
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 import bokeh.palettes as bp
 from mpl_toolkits.mplot3d import Axes3D
-from bokeh.layouts import gridplot, widgetbox
+from bokeh.layouts import gridplot
 from bokeh.plotting import figure, show
 from bokeh.transform import linear_cmap
-from bokeh.models.widgets import DataTable, TableColumn, NumberFormatter
 from bokeh.models import (
         ColumnDataSource,
         ColorBar,
@@ -430,8 +429,8 @@ class FrequencyResponseResults:
         # bokeh plot - create a new plot
         mag_plot = figure(
             tools="pan, box_zoom, wheel_zoom, reset, save",
-            width=1000,
-            height=450,
+            width=900,
+            height=400,
             title="Frequency Response - Magnitude",
             x_axis_label="Frequency",
             y_axis_label=y_axis_label,
@@ -439,143 +438,6 @@ class FrequencyResponseResults:
         mag_plot.xaxis.axis_label_text_font_size = "14pt"
         mag_plot.yaxis.axis_label_text_font_size = "14pt"
 
-        minspeed = 320
-        maxspeed = 450
-
-        idx_max = argrelextrema(mag[inp, out, :], np.greater)
-        wn = frequency_range[idx_max[0]]
-
-        AF_table = []
-        SM_table = []
-        SM_ref_table = []
-
-        for i, peak in enumerate(mag[inp, out, :][idx_max[0]]):
-            peak_n = 0.707 * peak
-            peak_aux = np.linspace(peak_n, peak_n, len(frequency_range))
-
-            idx = np.argwhere(
-                np.diff(np.sign(peak_aux - mag[inp, out, :]))
-            ).flatten()
-            idx = np.sort(np.append(idx, idx_max[0][i]))
-
-            # if speed range is not long enough to catch the magnitudes
-            try:
-                idx_aux = [
-                    list(idx).index(idx_max[0][i]) - 1,
-                    list(idx).index(idx_max[0][i]) + 1,
-                ]
-                idx = idx[idx_aux]
-            except IndexError:
-                idx = [list(idx).index(idx_max[0][i]) - 1,
-                       len(frequency_range) - 1
-                ]
-
-            # Amplification Factor (AF)
-            AF = wn[i] / (
-                frequency_range[idx[1]] - frequency_range[idx[0]]
-            )
-
-            # Separation Margin (SM)
-            if AF > 2.5 and wn[i] < minspeed:
-                SM = min([16, 17 * (1 - 1 / (AF - 1.5))]) / 100
-                SMspeed = wn[i] * (1 + SM)
-                SM_ref = (minspeed - wn[i]) / wn[i]
-                source = ColumnDataSource(
-                    dict(
-                        top=[max(mag[inp, out, :][idx_max[0]])],
-                        bottom=[0],
-                        left=[wn[i]],
-                        right=[SMspeed],
-                        tag1=[wn[i]],
-                        tag2=[SMspeed],
-                    )
-                )
-
-            elif AF > 2.5 and wn[i] > maxspeed:
-                SM = min([26, 10 + 17 * (1 - 1 / (AF - 1.5))]) /100
-                SMspeed = wn[i] * (1 - SM)
-                SM_ref = (wn[i] - maxspeed) / maxspeed
-                source = ColumnDataSource(
-                    dict(
-                        top=[max(mag[inp, out, :][idx_max[0]])],
-                        bottom=[0],
-                        left=[SMspeed],
-                        right=[wn[i]],
-                        tag1=[wn[i]],
-                        tag2=[SMspeed],
-                    )
-                )
-
-            else:
-                SM = None
-                SM_ref = None
-                SMspeed = None
-
-            AF_table.append(AF)
-            SM_table.append(SM)
-            SM_ref_table.append(SM_ref)
-
-            mag_plot.quad(
-                top="top",
-                bottom="bottom",
-                left="left",
-                right="right",
-                source=source,
-                line_color=bokeh_colors[8],
-                line_width=0.8,
-                fill_alpha=0.2,
-                fill_color=bokeh_colors[8],
-                legend="Separation Margin",
-                name="SM2",
-            )
-            hover = HoverTool(names=["SM2"])
-            hover.tooltips = [
-                ("Critical Speed :", "@tag1"),
-                ("Speed at 0.707 x peak amplitude :", "@tag2"),
-            ]
-
-        table_source = ColumnDataSource(
-                dict(
-                    Wd=frequency_range[idx_max[0]],
-                    SM_table=SM_table,
-                    AF_table=AF_table,
-                    SM_ref_table=SM_ref_table,
-                )
-        )
-        form1 = NumberFormatter(format='0.00')
-        form2 = NumberFormatter(format='0.00%')
-        columns = [
-                TableColumn(
-                    field="Wd", title="Critical Speed", formatter=form1
-                ),
-                TableColumn(
-                    field="SM_table", title="Required Separation Margin", formatter=form2
-                ),
-                TableColumn(
-                    field="SM_ref_table", title="Available Separation Margin", formatter=form2
-                ),
-                TableColumn(
-                    field="AF_table", title="Amplification Factor", formatter=form1
-                )
-            ]
-        data_table = DataTable(
-                source=table_source, columns=columns, width=700, height=450
-        )
-        table = widgetbox(data_table)
-
-        mag_plot.quad(
-            top=max(mag[inp, out, :][idx_max[0]]),
-            bottom=0,
-            left=minspeed,
-            right=maxspeed,
-            line_color="green",
-            line_width=0.8,
-            fill_alpha=0.2,
-            fill_color="green",
-            legend="Operation Speed Range",
-        )
-        mag_plot.add_tools(hover)
-
         source = ColumnDataSource(dict(x=frequency_range, y=mag[inp, out, :]))
         mag_plot.line(
             x="x",
@@ -586,7 +448,7 @@ class FrequencyResponseResults:
             line_width=3,
         )
 
-        return mag_plot, table
+        return mag_plot
 
     def plot_phase_matplotlib(self, inp, out, ax=None, **kwargs):
         """Plot frequency response.
@@ -661,8 +523,8 @@ class FrequencyResponseResults:
         # bokeh plot - create a new plot
         phase_plot = figure(
             tools="pan, box_zoom, wheel_zoom, reset, save",
-            width=1000,
-            height=450,
+            width=900,
+            height=400,
             title="Frequency Response - Phase",
             x_axis_label="Frequency",
             y_axis_label="Phase",
@@ -758,11 +620,11 @@ class FrequencyResponseResults:
         --------
         """
         # bokeh plot axes
-        bk_ax0, table = self.plot_magnitude_bokeh(inp, out, ax=ax0)
+        bk_ax0 = self.plot_magnitude_bokeh(inp, out, ax=ax0)
         bk_ax1 = self.plot_phase_bokeh(inp, out, ax=ax1)
 
         # show the bokeh plot results
-        grid_plots = gridplot([[bk_ax0, table], [bk_ax1]])
+        grid_plots = gridplot([[bk_ax0], [bk_ax1]])
         show(grid_plots)
 
         return bk_ax0, bk_ax1
@@ -1131,15 +993,19 @@ class ModeShapeResults:
         self.log_dec = log_dec
         self.kappa_modes = kappa_modes
 
-    def calc_mode_shape(self, mode=None, evec=None):
+    def plot(self, mode=None, evec=None, fig=None, ax=None):
+        if ax is None:
+            fig = plt.figure()
+            ax = fig.gca(projection="3d")
+
         evec0 = self.modes[:, mode]
         nodes = self.nodes
         nodes_pos = self.nodes_pos
+        kappa_modes = self.kappa_modes
         elements_length = self.elements_length
 
         modex = evec0[0::4]
         modey = evec0[1::4]
-
         xmax, ixmax = max(abs(modex)), np.argmax(abs(modex))
         ymax, iymax = max(abs(modey)), np.argmax(abs(modey))
 
@@ -1159,6 +1025,8 @@ class ModeShapeResults:
         y_circles = np.zeros((num_points, len(nodes)))
         z_circles_pos = np.zeros((num_points, len(nodes)))
 
+        kappa_mode = kappa_modes[mode]
+
         for node in nodes:
             x = modex[node] * circle
             x_circles[:, node] = np.real(x)
@@ -1193,35 +1061,23 @@ class ModeShapeResults:
 
             pos0 = nn * n
             pos1 = nn * (n + 1)
-
             xn[pos0:pos1] = Nx @ evec0[xx].real
             yn[pos0:pos1] = Ny @ evec0[yy].real
             zn[pos0:pos1] = (node_pos * onn + Le * zeta).reshape(nn)
 
-        return xn, yn, zn, x_circles, y_circles, z_circles_pos, nn
-
-    def plot(self, mode=None, evec=None, fig=None, ax=None):
-        if ax is None:
-            fig = plt.figure()
-            ax = fig.gca(projection="3d")
-
-        nodes = self.nodes
-        kappa_mode = self.kappa_modes[mode]
-        xn, yn, zn, xc, yc, zc_pos, nn = self.calc_mode_shape(mode=mode, evec=evec)
-
         for node in nodes:
             ax.plot(
-                xc[10:, node],
-                yc[10:, node],
-                zc_pos[10:, node],
+                x_circles[10:, node],
+                y_circles[10:, node],
+                z_circles_pos[10:, node],
                 color=kappa_mode[node],
                 linewidth=0.5,
                 zdir="x",
             )
             ax.scatter(
-                xc[10, node],
-                yc[10, node],
-                zc_pos[10, node],
+                x_circles[10, node],
+                y_circles[10, node],
+                z_circles_pos[10, node],
                 s=5,
                 color=kappa_mode[node],
                 zdir="x",
@@ -1260,7 +1116,6 @@ class StaticResults():
         nodes,
         nodes_pos,
         Vx_axis,
-        force_data,
     ):
 
         self.disp_y = disp_y
@@ -1272,7 +1127,6 @@ class StaticResults():
         self.nodes = nodes
         self.nodes_pos = nodes_pos
         self.Vx_axis = Vx_axis
-        self.force_data = force_data
 
     def plot(self):
         """Plot static analysis graphs.
diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py
index 904f0fb1..d0d801fa 100644
--- a/ross/rotor_assembly.py
+++ b/ross/rotor_assembly.py
@@ -185,6 +185,16 @@ class Rotor(object):
         if bearing_seal_elements is None:
             bearing_seal_elements = []
 
+        for i, disk in enumerate(disk_elements):
+            if disk.tag is None:
+                disk.tag = "Disk " + str(i)
+
+        for i, brg in enumerate(bearing_seal_elements):
+            if brg.__class__.__name__ == "BearingElement" and brg.tag is None:
+                brg.tag = "Bearing " + str(i)
+            if brg.__class__.__name__ == "SealElement" and brg.tag is None:
+                brg.tag = "Seal " + str(i)
+
         self.shaft_elements = shaft_elements
         self.bearing_seal_elements = bearing_seal_elements
         self.disk_elements = disk_elements
@@ -212,6 +222,7 @@ class Rotor(object):
             "rho",
             "volume",
             "m",
+            "tag",
         ]
 
         df_shaft = pd.DataFrame([el.summary() for el in self.shaft_elements])
@@ -1882,7 +1893,8 @@ class Rotor(object):
                     -disp_y[node] * self.df_bearings.loc[i, "kyy"].coefficient[0]
                 )
                 BrgForceToReturn.append(
-                    np.around(-disp_y[node] * self.df_bearings.loc[i, "kyy"].coefficient[0], decimals=1)
+                    "%.1f"
+                    % (-disp_y[node] * self.df_bearings.loc[i, "kyy"].coefficient[0])
                 )
 
             # Disk Forces
@@ -1891,7 +1903,7 @@ class Rotor(object):
                 for i, node in enumerate(self.df_disks["n"]):
                     DskForce[node] = self.df_disks.loc[i, "m"] * -9.8065
                     DskForceToReturn.append(
-                        np.around(self.df_disks.loc[i, "m"] * -9.8065, decimals=1)
+                        "%.1f" % (self.df_disks.loc[i, "m"] * -9.8065)
                     )
 
             # Shaft Weight Forces
@@ -1943,15 +1955,13 @@ class Rotor(object):
                 Bm = np.append(Bm, Bm[i] + Mx[i])
             self.Bm = Bm
 
-            sh_weight = np.around(
-                    sum(self.df_shaft["m"].values) * 9.8065, decimals=1
-            )
+            sh_weight = sum(self.df_shaft["m"].values) * 9.8065
 
             force_data = {
                 "Static displacement vector": disp_y,
                 "Shearing force vector": Vx,
                 "Bending moment vector": Bm,
-                "Shaft Total Weight": sh_weight,
+                "Shaft Total Weight": "%.1f" % sh_weight,
                 "Disks Forces": DskForceToReturn,
                 "Bearings Reaction Forces": BrgForceToReturn,
             }
@@ -1966,7 +1976,6 @@ class Rotor(object):
                 self.nodes,
                 self.nodes_pos,
                 Vx_axis,
-                force_data,
             )
 
         return results
diff --git a/ross/shaft_element.py b/ross/shaft_element.py
index 8cd6519c..6d45226d 100644
--- a/ross/shaft_element.py
+++ b/ross/shaft_element.py
@@ -695,7 +695,6 @@ class ShaftTaperedElement(Element):
     .. math:: [x_0, y_0, \alpha_0, \beta_0, x_1, y_1, \alpha_1, \beta_1]^T
     Where :math:`\alpha_0` and :math:`\alpha_1` are the bending on the yz plane
     and :math:`\beta_0` and :math:`\beta_1` are the bending on the xz plane.
-
     Parameters
     ----------
     L : float
@@ -731,44 +730,6 @@ class ShaftTaperedElement(Element):
     shear_method_calc : string, optional
         Determines which shear calculation method the user will adopt
         Default is 'cowper'
-
-    Returns
-    -------
-    A shaft element
-
-    Attributes
-    ----------
-    Poisson : float
-        Poisson coefficient for the element.
-    A : float
-        Element mid-section area.
-    Ie : float
-        Ie is the second moment of area of the mid cross section about
-        the neutral plane Ie = pi*r**2/4
-    phi : float
-        Constant that is used according to [1] to consider rotary
-        inertia and shear effects. If these are not considered phi=0.
-    References
-    ----------
-    .. [1] 'Dynamics of Rotating Machinery' by MI Friswell, JET Penny, SD Garvey
-       & AW Lees, published by Cambridge University Press, 2010 pp. 158-166.
-    Examples
-    --------
-    >>> from ross.materials import steel
-    >>> Euler_Bernoulli_Element = ShaftTaperedElement(
-    ...                         L=0.5, i_d_l=0.05, i_d_r=0.05, o_d_l=0.1,
-    ...                         o_d_r=0.15, material=steel,
-    ...                         rotary_inertia=False,
-    ...                         shear_effects=False)
-    >>> Euler_Bernoulli_Element.phi
-    0
-    >>> Timoshenko_Element = ShaftTaperedElement(
-    ...                         L=0.5, i_d_l=0.05, i_d_r=0.05, o_d_l=0.1,
-    ...                         o_d_r=0.15, material=steel,
-    ...                         rotary_inertia=True,
-    ...                         shear_effects=True)
-    >>> Timoshenko_Element.phi
-    0.21026260191948362
     """
 
     def __init__(
@@ -809,6 +770,7 @@ class ShaftTaperedElement(Element):
 
         self.L = float(L)
 
+        # diameters
         self.o_d = (float(o_d_l) + float(o_d_r)) / 2
         self.i_d = (float(i_d_l) + float(i_d_r)) / 2
         self.i_d_l = float(i_d_l)
@@ -843,7 +805,7 @@ class ShaftTaperedElement(Element):
         b1 = np.pi * (delta_ro ** 2 - delta_ri ** 2) / A_l
         b2 = 3 * np.pi * (roj ** 2 * delta_ro ** 2 - rij ** 2 * delta_ri ** 2) / (2 * Ie_l)
         gama = np.pi * (roj * delta_ro ** 3 - rij * delta_ri ** 3) / Ie_l
-        delta = np.pi * (delta_ro ** 4 - delta_ri ** 4) / (4 * Ie_l)
+        delta = np.pi * (delta_ro ** 4 - delta_ri ** 4) / 4 * Ie_l
 
         self.a1 = a1
         self.a2 = a2
@@ -873,7 +835,7 @@ class ShaftTaperedElement(Element):
         # hutchinson - kappa as per Hutchinson (2001)
         # cowper - kappa as per Cowper (1996)
         if shear_effects:
-            r = ((i_d_l + i_d_r) / 2) / ((o_d_l + o_d_r) / 2)
+            r = i_d_l / o_d_l
             r2 = r * r
             r12 = (1 + r2) ** 2
             if shear_method_calc == "hutchinson":
@@ -979,11 +941,9 @@ class ShaftTaperedElement(Element):
 
     def M(self):
         r"""Mass matrix for an instance of a shaft tapered element.
-
         Returns
         -------
         Mass matrix for the shaft tapered element.
-
         Examples
         --------
         >>> Timoshenko_Element = ShaftTaperedElement(
@@ -992,10 +952,10 @@ class ShaftTaperedElement(Element):
         ...                         rotary_inertia=True,
         ...                         shear_effects=True)
         >>> Timoshenko_Element.M()[:4, :4]
-        array([[11.36986417,  0.        ,  0.        ,  0.86197637],
-               [ 0.        , 11.36986417, -0.86197637,  0.        ],
-               [ 0.        , -0.86197637,  0.08667495,  0.        ],
-               [ 0.86197637,  0.        ,  0.        ,  0.08667495]])
+        array([[11.35295267,  0.        ,  0.        ,  0.85960486],
+               [ 0.        , 11.35295267, -0.85960486,  0.        ],
+               [ 0.        , -0.85960486,  0.08653475,  0.        ],
+               [ 0.85960486,  0.        ,  0.        ,  0.08653475]])
         """
         phi = self.phi
         L = self.L
@@ -1130,11 +1090,9 @@ class ShaftTaperedElement(Element):
 
     def K(self):
         r"""Stiffness matrix for an instance of a shaft tapered element.
-
         Returns
         -------
         Stiffness matrix for the shaft tapered element.
-
         Examples
         --------
         >>> from ross.materials import steel
@@ -1144,13 +1102,13 @@ class ShaftTaperedElement(Element):
         ...                         rotary_inertia=True,
         ...                         shear_effects=True)
         >>> Timoshenko_Element.K()[:4, :4]/1e6
-        array([[209.25641985,   0.        ,   0.        ,  38.62129051],
-               [  0.        , 209.25641985, -38.62129051,   0.        ],
-               [  0.        , -38.62129051,  11.56619973,   0.        ],
-               [ 38.62129051,   0.        ,   0.        ,  11.56619973]])
+        array([[211.30832227,   0.        ,   0.        ,  54.96237663],
+               [  0.        , 211.30832227, -54.96237663,   0.        ],
+               [  0.        , -54.96237663,  15.71572869,   0.        ],
+               [ 54.96237663,   0.        ,   0.        ,  15.71572869]])
         """
-        L = self.L
         phi = self.phi
+        L = self.L
         a1 = self.a1
         a2 = self.a2
         b1 = self.b1
@@ -1158,9 +1116,10 @@ class ShaftTaperedElement(Element):
         delta = self.delta
         gama = self.gama
         Ie_l = self.Ie_l
+        A_l = self.A_l
 
         # fmt: off
-        k1 = 1260 + 630 * a2 + 504 * b2 + 441 * gama + 396 * delta
+        k1 = 630 * a2 + 504 * b2 + 441 * gama + 396 * delta
         k2 = (
             630
             + 210 * a2
@@ -1202,7 +1161,7 @@ class ShaftTaperedElement(Element):
         k8 = 6 + 3 * a1 + 2 * b1
         k9 = 3 + 1.5 * a1 + b1
 
-        K1 = np.array([
+        K = np.array([
             [  k1,     0,       0,    L*k2,   -k1,    0,       0,    L*k3],
             [   0,    k1,   -L*k2,       0,     0,  -k1,   -L*k3,       0],
             [   0, -L*k2, L**2*k4,       0,     0, L*k2, L**2*k5,       0],
@@ -1213,18 +1172,20 @@ class ShaftTaperedElement(Element):
             [L*k3,     0,       0, L**2*k5, -L*k3,    0,       0, L**2*k6],
         ])
 
+        K = self.material.E * Ie_l * K / (105 * L ** 3 * (1 + phi) ** 2)
+
         K2 = np.array([
             [  k7,     0,       0,    L*k8,   -k7,     0,       0,    L*k8],
             [   0,    k7,   -L*k8,       0,     0,   -k7,   -L*k8,       0],
             [   0, -L*k8, L**2*k9,       0,     0,  L*k8, L**2*k9,       0],
             [L*k8,     0,       0, L**2*k9, -L*k8,     0,       0, L**2*k9],
-            [ -k7,     0,       0,   -L*k8,    k7,     0,       0,   -L*k8],
-            [   0,   -k7,    L*k8,       0,     0,    k7,    L*k8,       0],
-            [   0, -L*k8, L**2*k9,       0,     0,  L*k8, L**2*k9,       0],
-            [L*k8,     0,       0, L**2*k9, -L*k8,     0,       0, L**2*k9],
+            [ -k7,     0,       0,   -L*k8,    k7,     0,       0,    L*k8],
+            [   0,   -k7,    L*k8,       0,     0,    k7,   -L*k8,       0],
+            [   0, -L*k8, L**2*k9,       0,     0, -L*k8, L**2*k9,       0],
+            [L*k8,     0,       0, L**2*k9,  L*k8,     0,       0, L**2*k9],
         ])
-
-        K = self.material.E * Ie_l / (105 * L ** 3 * (1 + phi) ** 2) * (K1 + 105 * phi * K2)
+        K2 = self.material.G_s * A_l * phi ** 2 * K2 / (12 * self.kappa * L * (1 + phi) ** 2)
+        K = K + K2
         # fmt: on
 
         return K
@@ -1239,17 +1200,6 @@ class ShaftTaperedElement(Element):
 
         Examples
         --------
-        >>> from ross.materials import steel
-        >>> Timoshenko_Element = ShaftTaperedElement(
-        ...                         0.5, 0.05, 0.05, 0.1,
-        ...                         0.15, steel,
-        ...                         rotary_inertia=True,
-        ...                         shear_effects=True)
-        >>> Timoshenko_Element.C()[:4, :4]
-        array([[0., 0., 0., 0.],
-               [0., 0., 0., 0.],
-               [0., 0., 0., 0.],
-               [0., 0., 0., 0.]])
         """
         C = np.zeros((8, 8))
 
@@ -1257,11 +1207,9 @@ class ShaftTaperedElement(Element):
 
     def G(self):
         """Gyroscopic matrix for an instance of a shaft tapred element.
-
         Returns
         -------
         Gyroscopic matrix for the shaft tapered element.
- 
         Examples
         --------
         >>> from ross.materials import steel
@@ -1272,10 +1220,10 @@ class ShaftTaperedElement(Element):
         ...                         rotary_inertia=True,
         ...                         shear_effects=True)
         >>> Timoshenko_Element.G()[:4, :4]
-        array([[ 0.        ,  0.30940809, -0.01085902,  0.        ],
-               [-0.30940809,  0.        ,  0.        , -0.01085902],
-               [ 0.01085902,  0.        ,  0.        ,  0.0067206 ],
-               [ 0.        ,  0.01085902, -0.0067206 ,  0.        ]])
+        array([[ 0.        ,  0.2998902 , -0.00955058,  0.        ],
+               [-0.2998902 ,  0.        ,  0.        , -0.00955058],
+               [ 0.00955058,  0.        ,  0.        ,  0.00663842],
+               [ 0.        ,  0.00955058, -0.00663842,  0.        ]])
         """
 
         G = np.zeros((8, 8))


Test patch:
diff --git a/ross/tests/test_bearing_seal_element.py b/ross/tests/test_bearing_seal_element.py
index 2390fbb2..a698b62d 100644
--- a/ross/tests/test_bearing_seal_element.py
+++ b/ross/tests/test_bearing_seal_element.py
@@ -119,7 +119,7 @@ def test_bearing1_matrices(bearing1):
     assert_allclose(bearing1.C(314.2), C)
 
 
-def test_bearing_error_speed_not_given():
+def test_bearing_error1():
     speed = np.linspace(0, 10000, 5)
     kx = 1e8 * speed
     cx = 1e8 * speed
@@ -157,38 +157,6 @@ def test_bearing_constant(bearing_constant):
     assert_allclose(bearing_constant.cxx.interpolated(300.9), 0, rtol=1e5)
 
 
-def test_bearing_len_2():
-    bearing = BearingElement(
-        n=0,
-        kxx=[481, 4810],
-        cxx=[3.13, 10.81],
-        kyy=[481, 4810],
-        kxy=[194, 2078],
-        kyx=[-194, -2078],
-        cyy=[3.13, 10.81],
-        cxy=[0.276, 0.69],
-        cyx=[-0.276, -0.69],
-        w=[115.19, 345.575],
-    )
-    assert_allclose(bearing.kxx.interpolated(115.19), 481, rtol=1e5)
-
-
-def test_bearing_len_3():
-    bearing = BearingElement(
-        n=0,
-        kxx=[481, 4810, 18810],
-        cxx=[3.13, 10.81, 22.99],
-        kyy=[481, 4810, 18810],
-        kxy=[194, 2078, 8776],
-        kyx=[-194, -2078, -8776],
-        cyy=[3.13, 10.81, 22.99],
-        cxy=[0.276, 0.69, 1.19],
-        cyx=[-0.276, -0.69, -1.19],
-        w=[115.19, 345.575, 691.15],
-    )
-    assert_allclose(bearing.kxx.interpolated(115.19), 481, rtol=1e5)
-
-
 def test_equality(bearing0, bearing1, bearing_constant):
     assert bearing0 == bearing0
     assert bearing0 == bearing1
@@ -197,9 +165,7 @@ def test_equality(bearing0, bearing1, bearing_constant):
 
 
 def test_from_table():
-    bearing_file = (
-        os.path.dirname(os.path.realpath(__file__)) + "/data/bearing_seal.xls"
-    )
+    bearing_file = os.path.dirname(os.path.realpath(__file__)) + "/data/bearing_seal.xls"
     bearing = BearingElement.from_table(0, bearing_file)
     assert bearing.n == 0
     assert math.isclose(bearing.w[2], 523.6, abs_tol=0.1)
diff --git a/ross/tests/test_shaft_element.py b/ross/tests/test_shaft_element.py
index 76183a51..ee1731d3 100644
--- a/ross/tests/test_shaft_element.py
+++ b/ross/tests/test_shaft_element.py
@@ -153,8 +153,6 @@ def test_from_table():
         assert shaft[9].n == 8
         assert_allclose(shaft[8].material.E, 206842718795.05, rtol=3e-06)
         assert_allclose(shaft[9].material.E, 6894.75, atol=0.008)
-
-
 # Shaft Tapered Element tests
 @pytest.fixture
 def tap_tim():
@@ -194,51 +192,49 @@ def test_parameters_tap_tim(tap_tim):
     assert tap_tim.material.rho == 7810
     assert_almost_equal(tap_tim.material.Poisson, 0.29926108)
     assert_almost_equal(tap_tim.A, 0.024052818754)
-    assert_almost_equal(tap_tim.Ie * 1e5,  4.60385984)
-    assert_almost_equal(tap_tim.phi, 0.4208816002)
+    assert_almost_equal(tap_tim.Ie * 1e5, 4.44854429)
+    assert_almost_equal(tap_tim.phi, 0.40668276)
 
 
 def test_mass_matrix_tap_tim(tap_tim):
     # fmt: off
-    M0e_tim = np.array([
-            [41.84146162,  0.        ,  0.        ,  1.85735573, 10.27465893,  0.        ,  0.        , -1.16147981],
-            [ 0.        , 41.84146162, -1.85735573,  0.        ,  0.        , 10.27465893,  1.16147981,  0.        ],
-            [ 0.        , -1.85735573,  0.18686024,  0.        ,  0.        , -1.04154454, -0.10311193,  0.        ],
-            [ 1.85735573,  0.        ,  0.        ,  0.18686024,  1.04154454,  0.        ,  0.        , -0.10311193],
-            [10.27465893,  0.        ,  0.        ,  1.04154454, 17.35069604,  0.        ,  0.        , -1.1330391 ],
-            [ 0.        , 10.27465893, -1.04154454,  0.        ,  0.        , 17.35069604,  1.1330391 ,  0.        ],
-            [ 0.        ,  1.16147981, -0.10311193,  0.        ,  0.        ,  1.1330391 ,  0.09851808,  0.        ],
-            [-1.16147981,  0.        ,  0.        , -0.10311193, -1.1330391 ,  0.        ,  0.        ,  0.09851808]])
+    M0e_tim = np.array([[41.85725794,  0.        ,  0.        ,  1.8603991 , 10.27411395,  0.        ,  0.        , -1.14850813],
+                        [ 0.        , 41.85725794, -1.8603991 ,  0.        ,  0.        , 10.27411395,  1.14850813,  0.        ],
+                        [ 0.        , -1.8603991 ,  0.18670638,  0.        ,  0.        , -1.04155143, -0.10238917,  0.        ],
+                        [ 1.8603991 ,  0.        ,  0.        ,  0.18670638,  1.04155143,  0.        ,  0.        , -0.10238917],
+                        [10.27411395,  0.        ,  0.        ,  1.04155143, 17.33598969,  0.        ,  0.        , -1.14296052],
+                        [ 0.        , 10.27411395, -1.04155143,  0.        ,  0.        , 17.33598969,  1.14296052,  0.        ],
+                        [ 0.        ,  1.14850813, -0.10238917,  0.        ,  0.        ,  1.14296052,  0.09198448,  0.        ],
+                        [-1.14850813,  0.        ,  0.        , -0.10238917, -1.14296052,  0.        ,  0.        ,  0.09198448]])
     # fmt: on
     assert_almost_equal(tap_tim.M(), M0e_tim, decimal=5)
 
 
+@pytest.mark.skip(reason="inconsistencies in stiffness matrix")
 def test_stiffness_matrix_tap_tim(tap_tim):
     # fmt: off
-    K0e_tim = np.array([
-            [ 23.3411494 ,   0.        ,   0.        ,   6.2946873 , -23.3411494 ,   0.        ,   0.        ,   3.04177246],
-            [  0.        ,  23.3411494 ,  -6.2946873 ,   0.        ,   0.        , -23.3411494 ,  -3.04177246,   0.        ],
-            [  0.        ,  -6.2946873 ,   1.91793266,   0.        ,   0.        ,   6.2946873 ,   0.59994226,   0.        ],
-            [  6.2946873 ,   0.        ,   0.        ,   1.91793266,  -6.2946873 ,   0.        ,   0.        ,   0.59994226],
-            [-23.3411494 ,   0.        ,   0.        ,  -6.2946873 ,  23.3411494 ,   0.        ,   0.        ,  -3.04177246],
-            [  0.        , -23.3411494 ,   6.2946873 ,   0.        ,   0.        ,  23.3411494 ,   3.04177246,   0.        ],
-            [  0.        ,  -3.04177246,   0.59994226,   0.        ,   0.        ,   3.04177246,   0.61676672,   0.        ],
-            [  3.04177246,   0.        ,   0.        ,   0.59994226,  -3.04177246,   0.        ,   0.        ,   0.61676672]])
+    K0e_tim = np.array([[-19.5879458 ,   0.        ,   0.        ,   5.57911856,  19.5879458 ,   0.        ,   0.        ,   1.92061626],
+                        [  0.        , -19.5879458 ,  -5.57911856,   0.        ,   0.        ,  19.5879458 ,  -1.92061626,   0.        ],
+                        [  0.        ,  -5.57911856,   1.78916039,   0.        ,   0.        ,   5.57911856,   0.44248703,   0.        ],
+                        [  5.57911856,   0.        ,   0.        ,   1.78916039,  -5.57911856,   0.        ,   0.        ,   0.44248703],
+                        [ 19.5879458 ,   0.        ,   0.        ,  -5.57911856, -19.5879458 ,   0.        ,   0.        ,   0.03401395],
+                        [  0.        ,  19.5879458 ,   5.57911856,   0.        ,   0.        , -19.5879458 ,  -0.03401395,   0.        ],
+                        [  0.        ,  -1.92061626,   0.44248703,   0.        ,   0.        ,  -0.03401395,   0.32575947,   0.        ],
+                        [  1.92061626,   0.        ,   0.        ,   0.44248703,   0.03401395,   0.        ,   0.        ,   0.32575947]])
     # fmt: on
     assert_almost_equal(tap_tim.K() / 1e8, K0e_tim, decimal=5)
 
 
 def test_gyroscopic_matrix_tap_tim(tap_tim):
     # fmt: off
-    G0e_tim = np.array([
-            [ 0.        ,  1.23853234,  0.15544731,  0.        ,  0.        , -1.23853234, -0.03173553,  0.        ],
-            [-1.23853234,  0.        ,  0.        ,  0.15544731,  1.23853234,  0.        ,  0.        , -0.03173553],
-            [-0.15544731,  0.        ,  0.        ,  0.11850408,  0.15544731,  0.        ,  0.        , -0.01563679],
-            [ 0.        , -0.15544731, -0.11850408,  0.        ,  0.        ,  0.15544731,  0.01563679,  0.        ],
-            [ 0.        , -1.23853234, -0.15544731,  0.        ,  0.        ,  1.23853234,  0.03173553,  0.        ],
-            [ 1.23853234,  0.        ,  0.        , -0.15544731, -1.23853234,  0.        ,  0.        ,  0.03173553],
-            [ 0.03173553,  0.        ,  0.        , -0.01563679, -0.03173553,  0.        ,  0.        ,  0.01089192],
-            [ 0.        ,  0.03173553,  0.01563679,  0.        ,  0.        , -0.03173553, -0.01089192,  0.        ]])
+    G0e_tim = np.array([[ 0.        ,  1.19360236,  0.15864756,  0.        ,  0.        , -1.19360236, -0.05013491,  0.        ],
+                        [-1.19360236,  0.        ,  0.        ,  0.15864756,  1.19360236,  0.        ,  0.        , -0.05013491],
+                        [-0.15864756,  0.        ,  0.        ,  0.11754246,  0.15864756,  0.        ,  0.        , -0.01449657],
+                        [ 0.        , -0.15864756, -0.11754246,  0.        ,  0.        ,  0.15864756,  0.01449657,  0.        ],
+                        [ 0.        , -1.19360236, -0.15864756,  0.        ,  0.        ,  1.19360236,  0.05013491,  0.        ],
+                        [ 1.19360236,  0.        ,  0.        , -0.15864756, -1.19360236,  0.        ,  0.        ,  0.05013491],
+                        [ 0.05013491,  0.        ,  0.        , -0.01449657, -0.05013491,  0.        ,  0.        , -0.00213197],
+                        [ 0.        ,  0.05013491,  0.01449657,  0.        ,  0.        , -0.05013491,  0.00213197,  0.        ]])
     # fmt: on
     assert_almost_equal(tap_tim.G(), G0e_tim, decimal=5)
 
@@ -274,6 +270,7 @@ def test_match_mass_matrix(tap2, tim2):
     assert_almost_equal(M_tap, M_tim, decimal=5)
 
 
+@pytest.mark.skip(reason="inconsistencies in stiffness matrix")
 def test_match_stiffness_matrix(tap2, tim2):
     K_tap = tap2.K()
     K_tim = tim2.K()