diff --git a/CHANGELOG.txt b/CHANGELOG.txt
index c9cc2509a84..ed65d0c1014 100644
--- a/CHANGELOG.txt
+++ b/CHANGELOG.txt
@@ -8,6 +8,9 @@ Changes:
  - obspy.clients.seishub:
    * submodule removed completely, since it is outdated and not even test
      servers have been running for years (see #2994)
+ - obspy.taup:
+   * add option "indicate_wave_type" to distinguish S waves in ray paths
+     plot by using wiggly lines for shear waves (see #3047)
 
 maintenance_1.3.x
 =================
diff --git a/misc/docs/source/tutorial/code_snippets/travel_time_body_waves.py b/misc/docs/source/tutorial/code_snippets/travel_time_body_waves.py
index 2fbdbc98241..e21ad498a7e 100644
--- a/misc/docs/source/tutorial/code_snippets/travel_time_body_waves.py
+++ b/misc/docs/source/tutorial/code_snippets/travel_time_body_waves.py
@@ -37,7 +37,7 @@
     ax = arrivals.plot_rays(plot_type='spherical',
                             legend=False, label_arrivals=True,
                             plot_all=True,
-                            show=False, ax=ax)
+                            show=False, ax=ax, indicate_wave_type=True)
 
 # Annotate regions
 ax.text(0, 0, 'Solid\ninner\ncore',
diff --git a/obspy/CONTRIBUTORS.txt b/obspy/CONTRIBUTORS.txt
index 620a6a06a55..fca4ce73efe 100644
--- a/obspy/CONTRIBUTORS.txt
+++ b/obspy/CONTRIBUTORS.txt
@@ -89,6 +89,8 @@ Reyes, Celso
 Ringler, Adam
 Rothenhäusler, Nicolas
 Russo, Emiliano
+Rudge, John F.
+Russell, Stuart
 Sales de Andrade, Elliott
 Satriano, Claudio
 Saul, Joachim
diff --git a/obspy/taup/tau.py b/obspy/taup/tau.py
index aa21298ab85..c5804687b1c 100644
--- a/obspy/taup/tau.py
+++ b/obspy/taup/tau.py
@@ -17,7 +17,7 @@
 from .taup_pierce import TauPPierce
 from .taup_time import TauPTime
 from .taup_geo import calc_dist, add_geo_to_arrivals
-from .utils import parse_phase_list
+from .utils import parse_phase_list, split_ray_path
 import obspy.geodetics.base as geodetics
 
 # Pretty paired colors. Reorder to have saturated colors first and remove
@@ -246,7 +246,8 @@ def plot_times(self, phase_list=None, plot_all=True, legend=False,
 
     def plot_rays(self, phase_list=None, plot_type="spherical",
                   plot_all=True, legend=False, label_arrivals=False,
-                  show=True, fig=None, ax=None):
+                  show=True, fig=None, ax=None,
+                  indicate_wave_type=False):
         """
         Plot ray paths if any have been calculated.
 
@@ -284,6 +285,9 @@ def plot_rays(self, phase_list=None, plot_type="spherical",
             will be created. Must be a polar axes for the spherical plot and
             a regular one for the Cartesian plot.
         :type ax: :class:`matplotlib.axes.Axes`
+        :param indicate_wave_type: Distinguish between p and s waves when
+            plotting ray path. s waves indicated by wiggly lines.
+        :type indicate_wave_type: bool
         :returns: Matplotlib axes with the plot
         :rtype: :class:`matplotlib.axes.Axes`
         """
@@ -312,7 +316,7 @@ def plot_rays(self, phase_list=None, plot_type="spherical",
             distance = arrival.distance
             if distance < 0:
                 distance = (distance % 360)
-            if abs(dist - distance) / dist > 1E-5:
+            if abs(dist - distance) > 1E-5 * dist:
                 if plot_all is False:
                     continue
                 # Mirror on axis.
@@ -362,10 +366,28 @@ def plot_rays(self, phase_list=None, plot_type="spherical",
             radius = self.model.radius_of_planet
             for ray in arrivals:
                 color = colors.get(ray.name, 'k')
-                # Requires interpolation,or diffracted phases look funny.
-                ax.plot(intp(ray.path["dist"], 100),
-                        radius - intp(ray.path["depth"], 100),
-                        color=color, label=ray.name, lw=2.0)
+
+                # Requires interpolation, or diffracted phases look funny.
+                if indicate_wave_type:
+                    # Plot p, s, and diff phases separately
+                    paths, waves = split_ray_path(ray.path, self.model)
+                    for path, wave in zip(paths, waves):
+                        if wave == "s":
+                            # Make s waves wiggly
+                            with mpl.rc_context({'path.sketch': (2, 10, 1)}):
+                                ax.plot(intp(path["dist"], 100),
+                                        radius - intp(path["depth"], 100),
+                                        color=color, label=ray.name, lw=1.5)
+                        else:
+                            # p and diff waves
+                            ax.plot(intp(path["dist"], 100),
+                                    radius - intp(path["depth"], 100),
+                                    color=color, label=ray.name, lw=2.0)
+                else:
+                    # Plot each ray as a single path
+                    ax.plot(intp(ray.path["dist"], 100),
+                            radius - intp(ray.path["depth"], 100),
+                            color=color, label=ray.name, lw=2.0)
                 ax.set_yticks(radius - discons)
                 ax.xaxis.set_major_formatter(plt.NullFormatter())
                 ax.yaxis.set_major_formatter(plt.NullFormatter())
@@ -442,8 +464,25 @@ def plot_rays(self, phase_list=None, plot_type="spherical",
             # Plot the ray paths:
             for ray in arrivals:
                 color = colors.get(ray.name, 'k')
-                ax.plot(np.rad2deg(ray.path["dist"]), ray.path["depth"],
-                        color=color, label=ray.name, lw=2.0)
+                if indicate_wave_type:
+                    # Plot p, s, and diff phases separately
+                    paths, waves = split_ray_path(ray.path, self.model)
+                    for path, wave in zip(paths, waves):
+                        if wave == "s":
+                            # Make s waves wiggly
+                            with mpl.rc_context({'path.sketch': (2, 10, 1)}):
+                                ax.plot(np.rad2deg(path["dist"]),
+                                        path["depth"],
+                                        color=color, label=ray.name, lw=1.5)
+                        else:
+                            # p and diff waves
+                            ax.plot(np.rad2deg(path["dist"]), path["depth"],
+                                    color=color, label=ray.name, lw=2.0)
+                else:
+                    # Plot each ray as a single path
+                    ax.plot(np.rad2deg(ray.path["dist"]), ray.path["depth"],
+                            color=color,
+                            label=ray.name, lw=2.0)
 
             # Pretty station marker:
             ms = 14
@@ -975,7 +1014,7 @@ def plot_ray_paths(source_depth, min_degrees=0, max_degrees=360, npoints=10,
                    plot_type='spherical', phase_list=['P', 'S', 'PP'],
                    model='iasp91', plot_all=True, legend=False,
                    label_arrivals=False, verbose=False, fig=None, show=True,
-                   ax=None):
+                   ax=None, indicate_wave_type=False):
     """
     Plot ray paths for seismic phases.
 
@@ -1018,6 +1057,9 @@ def plot_ray_paths(source_depth, min_degrees=0, max_degrees=360, npoints=10,
     :param ax: Axes to plot in. If not given, a new figure with an axes
         will be created.
     :type ax: :class:`matplotlib.axes.Axes`
+    :param indicate_wave_type: Distinguish between p and s waves when
+        plotting ray path. s waves indicated by wiggly lines.
+    :type indicate_wave_type: bool
     :returns: Matplotlib axes with the plot
     :rtype: :class:`matplotlib.axes.Axes`
 
@@ -1056,7 +1098,8 @@ def plot_ray_paths(source_depth, min_degrees=0, max_degrees=360, npoints=10,
                                            phase_list=phase_list)
             ax = arrivals.plot_rays(phase_list=phase_list, show=False,
                                     ax=ax, plot_type=plot_type,
-                                    plot_all=plot_all, legend=False)
+                                    plot_all=plot_all, legend=False,
+                                    indicate_wave_type=indicate_wave_type)
             plotted = True
         except ValueError:
             norays.append(degree)
diff --git a/obspy/taup/utils.py b/obspy/taup/utils.py
index 2a936646038..b98d731e4e4 100644
--- a/obspy/taup/utils.py
+++ b/obspy/taup/utils.py
@@ -5,6 +5,8 @@
 import inspect
 from pathlib import Path
 
+import numpy as np
+
 ROOT = Path(Path(inspect.getfile(inspect.currentframe())).resolve()).parent
 
 
@@ -54,3 +56,126 @@ def get_phase_names(phase_name):
         names.append(phase_name)
 
     return names
+
+
+def split_ray_path(path, model):
+    """
+    Split and label ray path according to type of wave.
+
+    :param path: Path taken by ray.
+    :type path: :class:`~numpy.ndarray`
+        (dtype = :class:`~obspy.taup.helper_classes.TimeDist`)
+    :param model: The model used to calculate the ray path.
+    :type model: :class:`obspy.taup.tau_model.TauModel`
+    :returns: A list of paths and a list of wave types.
+        Wave types are either "p", "s", or "diff".
+    :rtype: list[:class:`~numpy.ndarray`]
+        and list[str]
+
+    The ray path is split at all discontinuities in the model
+    and labelled according to wave type.
+    """
+
+    # Get discontinuity depths in the model in km
+    discs = model.s_mod.v_mod.get_discontinuity_depths()[:-1]
+
+    # Split path at discontinuity depths
+    depths = path["depth"]
+    is_disc = np.isin(depths, discs)
+    is_disc[0] = False  # Don't split on first point
+    idx = np.where(is_disc)[0]
+
+    # Split ray path, including start and end points in each segment
+    splitted = np.split(path, idx)
+    paths = [
+        np.append(s, splitted[i + 1][0]) for i, s in enumerate(splitted[:-1])
+    ]
+
+    # Classify the waves as p, s, or diff
+    wave_types = [_classify_path(p, model) for p in paths]
+
+    return paths, wave_types
+
+
+def _expected_delay_time(ray_param, depth0, depth1, wave, model):
+    """
+    Expected delay time between two depths for a given wave type (p or s).
+    """
+
+    # Convert depths to radii
+    radius0 = model.radius_of_planet - depth0
+    radius1 = model.radius_of_planet - depth1
+
+    # Velocity model from TauModel
+    v_mod = model.s_mod.v_mod
+
+    # Get velocities
+    if depth1 >= depth0:
+        v0 = v_mod.evaluate_below(depth0, wave)[0]
+        v1 = v_mod.evaluate_above(depth1, wave)[0]
+    else:
+        v0 = v_mod.evaluate_above(depth0, wave)[0]
+        v1 = v_mod.evaluate_below(depth1, wave)[0]
+
+    # Calculate time for segment if velocity non-zero
+    # - if velocity zero then return zero time
+    if v0 > 0.0:
+
+        eta0 = radius0 / v0
+        eta1 = radius1 / v1
+
+        def vertical_slowness(eta, p):
+            y = eta**2 - p**2
+            return np.sqrt(y * (y > 0))  # in s
+
+        n0 = vertical_slowness(eta0, ray_param)
+        n1 = vertical_slowness(eta1, ray_param)
+
+        if ray_param == 0.0:
+            return 0.5 * ((1.0 / v0) + (1.0 / v1)) * abs(radius1 - radius0)
+        return 0.5 * (n0 + n1) * abs(np.log(radius1 / radius0))
+
+    return 0.0
+
+
+def _classify_path(path, model):
+    """
+    Determine whether we have a p or s-wave path by comparing delay times.
+    """
+
+    # Examine just the first two points near the shallowest part of the path
+    if path[0]["depth"] < path[-1]["depth"]:
+        point0 = path[0]
+        point1 = path[1]
+    else:
+        point0 = path[-2]
+        point1 = path[-1]
+
+    # Ray parameter
+    ray_param = point0["p"]
+
+    # Depths
+    depth0 = point0["depth"]
+    depth1 = point1["depth"]
+
+    # If no change in depth then this is a diffracted/head wave segment
+    if depth0 == depth1:
+        return "diff"
+
+    # Delay time for this segment from ray path
+    travel_time = point1["time"] - point0["time"]
+    distance = abs(point0["dist"] - point1["dist"])
+    delay_time = travel_time - ray_param * distance
+
+    # Get the expected delay time for each wave type
+    delay_p = _expected_delay_time(ray_param, depth0, depth1, "p", model)
+    delay_s = _expected_delay_time(ray_param, depth0, depth1, "s", model)
+
+    # Difference between predictions and given delay times
+    error_p = (delay_p / delay_time) - 1.0
+    error_s = (delay_s / delay_time) - 1.0
+
+    # Check which wave type matches the given delay time the best
+    if abs(error_p) < abs(error_s):
+        return "p"
+    return "s"