Repository: python-control/python-control
Pull request #953 (2023-12-29 18:34:35)

Update pole/zero and root locus plots to use _map/_plot pattern
This PR updates pole/zero and root locus plots to use the _map/_plot calling pattern described in #645. The following patterns will work:
```
rldata = ct.root_locus_map(syslist[, klist])
lines = rldata.plot(plot_options)

lines = ct.root_locus_plot(response[, klist, plot_options])

ct.root_locus(sys)
ct.rlocus(sys)

ct.pole_zero_map(sys).plot()
ct.pole_zero_plot(sys)
ct.pzmap(sys)
```

Everything is mostly backwards compatible except that the outputs from a _plot function are now an array of Line2D objects instead of the roots and gains.  You can still get the original data using the _map function and there is some legacy processing if you used the `plot` keyword to try let some code work without changes.

Summary of changes:
* All pole/zero and root locus plots now have a `_map` function and a `_plot` function.  You can access the latter via the `.plot()` method on the response.
* The `_plot` functions accept either the output of the `_map` function or a list of systems (in which case the `_map` function is called internally).  This allows the common pattern of `pole_zero_plot(sys)`, `root_locus_plot(sys)` to work as expected.
* The interactive mode for root locus plots is on by default (it used to be there only within `sisotool`).  Clicking on a portion of the root locus diagram will generate markers at the locations on the loci corresponding to that gain and add a message above the plot giving the frequency and damping ratio for the point that was clicked.
* The `print_gain` keyword is replaced with `interactive`.
* The (deprecated) `Plot` keyword is now `plot` (and generates a warning, since it triggers the legacy return values).
* Recalculation of the root locus plot when zooming no longer works (you can still zoom in and out, you just don't get a recalculated curve).
* The short versions `pzmap` and `rlocus` are still there.
* The `control.matlab` version of `rlocus` returns `roots, gains` (compatible with MATLAB)
* Added unit tests plus user documentation.
* Some other small fixes, code streamlining, etc along the way.

**This PR is going to break existing code.** It would be great if a few people could try this out so that we can make sure we are OK with the changes here.

**Examples** (from the user documentation):

Pole/zero maps and root locus diagrams provide insights into system response based on the locations of system poles and zeros in the complex plane. The `pole_zero_map()` function returns the poles and zeros and can be used to generate a pole/zero plot:
```
sys = ct.tf([1, 2], [1, 2, 3], name='SISO transfer function')
response = ct.pole_zero_map(sys)
ct.pole_zero_plot(response)
```
![pzmap-siso_ctime-default](https://github.com/python-control/python-control/assets/293362/147e2d14-ae3c-4cd7-a00d-2929a7d39d60)

A root locus plot shows the location of the closed loop poles of a system as a function of the loop gain:
```
ct.root_locus_map(sys).plot()
```

![rlocus-siso_ctime-default](https://github.com/python-control/python-control/assets/293362/1464cd83-68e5-4703-a4d5-4de589e02fce)

The grid in the left hand plane shows lines of constant damping ratio as well as arcs corresponding to the frequency of the complex pole. The grid can be turned off using the `grid` keyword. Setting grid to `Fals`e will turn off the grid but show the real and imaginary axis. To completely remove all lines except the root loci, use `grid=’empty’`.

On systems that support interactive plots, clicking on a location on the root locus diagram will mark the pole locations on all branches of the diagram and display the gain and damping ratio below the plot title:

![rlocus-siso_ctime-clicked](https://github.com/python-control/python-control/assets/293362/4e7cc157-58b8-4f76-aca4-28e892eacc52)

Root locus diagrams are also supported for discrete time systems, in which case the grid is show inside the unit circle:
```
sysd = sys.sample(0.1)
ct.root_locus_plot(sysd)
```

![rlocus-siso_dtime-default](https://github.com/python-control/python-control/assets/293362/8bd17e0e-ff3d-4e7f-b0ea-0f073761dbbe)

Lists of systems can also be given, in which case the root locus diagram for each system is plotted in different colors:
```
sys1 = ct.tf([1], [1, 2, 1], name='sys1')
sys2 = ct.tf([1, 0.2], [1, 1, 3, 1, 1], name='sys2')
ct.root_locus_plot([sys1, sys2], grid=False)
```

![rlocus-siso_multiple-nogrid](https://github.com/python-control/python-control/assets/293362/623ca08c-2665-4657-8085-8ee9d645cbc7)



Patch:
diff --git a/control/freqplot.py b/control/freqplot.py
index 53351541..57f24f8d 100644
--- a/control/freqplot.py
+++ b/control/freqplot.py
@@ -162,7 +162,7 @@ def bode_plot(
         values with no plot.
     rcParams : dict
         Override the default parameters used for generating plots.
-        Default is set up config.default['freqplot.rcParams'].
+        Default is set by config.default['freqplot.rcParams'].
     wrap_phase : bool or float
         If wrap_phase is `False` (default), then the phase will be unwrapped
         so that it is continuously increasing or decreasing.  If wrap_phase is
diff --git a/control/grid.py b/control/grid.py
index 785ec274..ef999594 100644
--- a/control/grid.py
+++ b/control/grid.py
@@ -1,12 +1,21 @@
-import numpy as np
-from numpy import cos, sin, sqrt, linspace, pi, exp
+# grid.py - code to add gridlines to root locus and pole-zero diagrams
+#
+# This code generates grids for pole-zero diagrams (including root locus
+# diagrams).  Rather than just draw a grid in place, it uses the AxisArtist
+# package to generate a custom grid that will scale with the figure.
+#
+
 import matplotlib.pyplot as plt
-from mpl_toolkits.axisartist import SubplotHost
-from mpl_toolkits.axisartist.grid_helper_curvelinear \
-    import GridHelperCurveLinear
 import mpl_toolkits.axisartist.angle_helper as angle_helper
+import numpy as np
 from matplotlib.projections import PolarAxes
 from matplotlib.transforms import Affine2D
+from mpl_toolkits.axisartist import SubplotHost
+from mpl_toolkits.axisartist.grid_helper_curvelinear import \
+    GridHelperCurveLinear
+from numpy import cos, exp, linspace, pi, sin, sqrt
+
+from .iosys import isdtime
 
 
 class FormatterDMS(object):
@@ -65,14 +74,15 @@ class ModifiedExtremeFinderCycle(angle_helper.ExtremeFinderCycle):
         return lon_min, lon_max, lat_min, lat_max
 
 
-def sgrid():
+def sgrid(scaling=None):
     # From matplotlib demos:
     # https://matplotlib.org/gallery/axisartist/demo_curvelinear_grid.html
     # https://matplotlib.org/gallery/axisartist/demo_floating_axis.html
 
     # PolarAxes.PolarTransform takes radian. However, we want our coordinate
-    # system in degree
+    # system in degrees
     tr = Affine2D().scale(np.pi/180., 1.) + PolarAxes.PolarTransform()
+
     # polar projection, which involves cycle, and also has limits in
     # its coordinates, needs a special method to find the extremes
     # (min, max of the coordinate within the view).
@@ -89,6 +99,7 @@ def sgrid():
         tr, extreme_finder=extreme_finder, grid_locator1=grid_locator1,
         tick_formatter1=tick_formatter1)
 
+    # Set up an axes with a specialized grid helper
     fig = plt.gcf()
     ax = SubplotHost(fig, 1, 1, 1, grid_helper=grid_helper)
 
@@ -97,15 +108,20 @@ def sgrid():
     ax.axis[:].major_ticklabels.set_visible(visible)
     ax.axis[:].major_ticks.set_visible(False)
     ax.axis[:].invert_ticklabel_direction()
+    ax.axis[:].major_ticklabels.set_color('gray')
 
+    # Set up internal tickmarks and labels along the real/imag axes
     ax.axis["wnxneg"] = axis = ax.new_floating_axis(0, 180)
     axis.set_ticklabel_direction("-")
     axis.label.set_visible(False)
+
     ax.axis["wnxpos"] = axis = ax.new_floating_axis(0, 0)
     axis.label.set_visible(False)
+
     ax.axis["wnypos"] = axis = ax.new_floating_axis(0, 90)
     axis.label.set_visible(False)
-    axis.set_axis_direction("left")
+    axis.set_axis_direction("right")
+
     ax.axis["wnyneg"] = axis = ax.new_floating_axis(0, 270)
     axis.label.set_visible(False)
     axis.set_axis_direction("left")
@@ -119,43 +135,41 @@ def sgrid():
     ax.axis["bottom"].get_helper().nth_coord_ticks = 0
 
     fig.add_subplot(ax)
-
-    # RECTANGULAR X Y AXES WITH SCALE
-    # par2 = ax.twiny()
-    # par2.axis["top"].toggle(all=False)
-    # par2.axis["right"].toggle(all=False)
-    # new_fixed_axis = par2.get_grid_helper().new_fixed_axis
-    # par2.axis["left"] = new_fixed_axis(loc="left",
-    #                                   axes=par2,
-    #                                   offset=(0, 0))
-    # par2.axis["bottom"] = new_fixed_axis(loc="bottom",
-    #                                     axes=par2,
-    #                                     offset=(0, 0))
-    # FINISH RECTANGULAR
-
     ax.grid(True, zorder=0, linestyle='dotted')
 
-    _final_setup(ax)
+    _final_setup(ax, scaling=scaling)
     return ax, fig
 
 
-def _final_setup(ax):
+# Utility function used by all grid code
+def _final_setup(ax, scaling=None):
     ax.set_xlabel('Real')
     ax.set_ylabel('Imaginary')
-    ax.axhline(y=0, color='black', lw=1)
-    ax.axvline(x=0, color='black', lw=1)
-    plt.axis('equal')
+    ax.axhline(y=0, color='black', lw=0.25)
+    ax.axvline(x=0, color='black', lw=0.25)
 
+    # Set up the scaling for the axes
+    scaling = 'equal' if scaling is None else scaling
+    plt.axis(scaling)
 
-def nogrid():
-    f = plt.gcf()
-    ax = plt.axes()
 
-    _final_setup(ax)
-    return ax, f
+# If not grid is given, at least separate stable/unstable regions
+def nogrid(dt=None, ax=None, scaling=None):
+    fig = plt.gcf()
+    if ax is None:
+        ax = fig.gca()
+
+    # Draw the unit circle for discrete time systems
+    if isdtime(dt=dt, strict=True):
+        s = np.linspace(0, 2*pi, 100)
+        ax.plot(np.cos(s), np.sin(s), 'k--', lw=0.5, dashes=(5, 5))
 
+    _final_setup(ax, scaling=scaling)
+    return ax, fig
 
-def zgrid(zetas=None, wns=None, ax=None):
+# Grid for discrete time system (drawn, not rendered by AxisArtist)
+# TODO (at some point): think about using customized grid generator?
+def zgrid(zetas=None, wns=None, ax=None, scaling=None):
     """Draws discrete damping and frequency grid"""
 
     fig = plt.gcf()
@@ -206,5 +220,9 @@ def zgrid(zetas=None, wns=None, ax=None):
         ax.annotate(r"$\frac{"+num+r"\pi}{T}$", xy=(an_x, an_y),
                     xytext=(an_x, an_y), size=9)
 
-    _final_setup(ax)
+    # Set default axes to allow some room around the unit circle
+    ax.set_xlim([-1.1, 1.1])
+    ax.set_ylim([-1.1, 1.1])
+
+    _final_setup(ax, scaling=scaling)
     return ax, fig
diff --git a/control/iosys.py b/control/iosys.py
index 52262250..fbd5c1db 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -503,42 +503,64 @@ def common_timebase(dt1, dt2):
         raise ValueError("Systems have incompatible timebases")
 
 # Check to see if a system is a discrete time system
-def isdtime(sys, strict=False):
+def isdtime(sys=None, strict=False, dt=None):
     """
     Check to see if a system is a discrete time system.
 
     Parameters
     ----------
-    sys : I/O or LTI system
-        System to be checked
-    strict: bool (default = False)
-        If strict is True, make sure that timebase is not None
+    sys : I/O system, optional
+        System to be checked.
+    dt : None or number, optional
+        Timebase to be checked.
+    strict: bool, default=False
+        If strict is True, make sure that timebase is not None.
     """
 
-    # Check to see if this is a constant
+    # See if we were passed a timebase instead of a system
+    if sys is None:
+        if dt is None:
+            return True if not strict else False
+        else:
+            return dt > 0
+    elif dt is not None:
+        raise TypeError("passing both system and timebase not allowed")
+
+    # Check timebase of the system
     if isinstance(sys, (int, float, complex, np.number)):
-        # OK as long as strict checking is off
+        # Constants OK as long as strict checking is off
         return True if not strict else False
     else:
         return sys.isdtime(strict)
 
 
 # Check to see if a system is a continuous time system
-def isctime(sys, strict=False):
+def isctime(sys=None, dt=None, strict=False):
     """
     Check to see if a system is a continuous-time system.
 
     Parameters
     ----------
-    sys : I/O or LTI system
-        System to be checked
+    sys : I/O system, optional
+        System to be checked.
+    dt : None or number, optional
+        Timebase to be checked.
     strict: bool (default = False)
-        If strict is True, make sure that timebase is not None
+        If strict is True, make sure that timebase is not None.
     """
 
-    # Check to see if this is a constant
+    # See if we were passed a timebase instead of a system
+    if sys is None:
+        if dt is None:
+            return True if not strict else False
+        else:
+            return dt == 0
+    elif dt is not None:
+        raise TypeError("passing both system and timebase not allowed")
+
+    # Check timebase of the system
     if isinstance(sys, (int, float, complex, np.number)):
-        # OK as long as strict checking is off
+        # Constants OK as long as strict checking is off
         return True if not strict else False
     else:
         return sys.isctime(strict)
diff --git a/control/matlab/wrappers.py b/control/matlab/wrappers.py
index b63b19c7..0384215a 100644
--- a/control/matlab/wrappers.py
+++ b/control/matlab/wrappers.py
@@ -12,14 +12,14 @@ from ..xferfcn import tf
 from ..lti import LTI
 from ..exception import ControlArgument
 
-__all__ = ['bode', 'nyquist', 'ngrid', 'dcgain', 'connect']
+__all__ = ['bode', 'nyquist', 'ngrid', 'rlocus', 'pzmap', 'dcgain', 'connect']
 
 def bode(*args, **kwargs):
     """bode(syslist[, omega, dB, Hz, deg, ...])
 
     Bode plot of the frequency response.
 
-    Plots a bode gain and phase diagram
+    Plots a bode gain and phase diagram.
 
     Parameters
     ----------
@@ -195,6 +195,106 @@ def _parse_freqplot_args(*args):
     return syslist, omega, plotstyle, other
 
 
+# TODO: rewrite to call root_locus_map, without using legacy plot keyword
+def rlocus(*args, **kwargs):
+    """rlocus(sys[, klist, xlim, ylim, ...])
+
+    Root locus diagram.
+
+    Calculate the root locus by finding the roots of 1 + k * G(s) where G
+    is a linear system with transfer function num(s)/den(s) and each k is
+    an element of kvect.
+
+    Parameters
+    ----------
+    sys : LTI object
+        Linear input/output systems (SISO only, for now).
+    kvect : array_like, optional
+        Gains to use in computing plot of closed-loop poles.
+    xlim : tuple or list, optional
+        Set limits of x axis, normally with tuple
+        (see :doc:`matplotlib:api/axes_api`).
+    ylim : tuple or list, optional
+        Set limits of y axis, normally with tuple
+        (see :doc:`matplotlib:api/axes_api`).
+
+    Returns
+    -------
+    roots : ndarray
+        Closed-loop root locations, arranged in which each row corresponds
+        to a gain in gains.
+    gains : ndarray
+        Gains used.  Same as kvect keyword argument if provided.
+
+    Notes
+    -----
+    This function is a wrapper for :func:`~control.root_locus_plot`,
+    with legacy return arguments.
+
+    """
+    from ..rlocus import root_locus_plot
+
+    # Use the plot keyword to get legacy behavior
+    kwargs = dict(kwargs)       # make a copy since we modify this
+    if 'plot' not in kwargs:
+        kwargs['plot'] = True
+
+    # Turn off deprecation warning
+    with warnings.catch_warnings():
+        warnings.filterwarnings(
+            'ignore', message='.* return values of .* is deprecated',
+            category=DeprecationWarning)
+        retval = root_locus_plot(*args, **kwargs)
+
+    return retval
+
+
+# TODO: rewrite to call pole_zero_map, without using legacy plot keyword
+def pzmap(*args, **kwargs):
+    """pzmap(sys[, grid, plot])
+
+    Plot a pole/zero map for a linear system.
+
+    Parameters
+    ----------
+    sys: LTI (StateSpace or TransferFunction)
+        Linear system for which poles and zeros are computed.
+    plot: bool, optional
+        If ``True`` a graph is generated with Matplotlib,
+        otherwise the poles and zeros are only computed and returned.
+    grid: boolean (default = False)
+        If True plot omega-damping grid.
+
+    Returns
+    -------
+    poles: array
+        The system's poles.
+    zeros: array
+        The system's zeros.
+
+    Notes
+    -----
+    This function is a wrapper for :func:`~control.pole_zero_plot`,
+    with legacy return arguments.
+
+    """
+    from ..pzmap import pole_zero_plot
+
+    # Use the plot keyword to get legacy behavior
+    kwargs = dict(kwargs)       # make a copy since we modify this
+    if 'plot' not in kwargs:
+        kwargs['plot'] = True
+
+    # Turn off deprecation warning
+    with warnings.catch_warnings():
+        warnings.filterwarnings(
+            'ignore', message='.* return values of .* is deprecated',
+            category=DeprecationWarning)
+        retval = pole_zero_plot(*args, **kwargs)
+
+    return retval
+
+
 from ..nichols import nichols_grid
 def ngrid():
     return nichols_grid()
@@ -254,6 +354,7 @@ def dcgain(*args):
 
 from ..bdalg import connect as ct_connect
 def connect(*args):
+
     """Index-based interconnection of an LTI system.
 
     The system `sys` is a system typically constructed with `append`, with
diff --git a/control/pzmap.py b/control/pzmap.py
index 5ee3d37c..59f378c3 100644
--- a/control/pzmap.py
+++ b/control/pzmap.py
@@ -1,133 +1,609 @@
 # pzmap.py - computations involving poles and zeros
 #
-# Author: Richard M. Murray
+# Original author: Richard M. Murray
 # Date: 7 Sep 2009
 #
 # This file contains functions that compute poles, zeros and related
-# quantities for a linear system.
-#
-# Copyright (c) 2009 by California Institute of Technology
-# All rights reserved.
-#
-# Redistribution and use in source and binary forms, with or without
-# modification, are permitted provided that the following conditions
-# are met:
-#
-# 1. Redistributions of source code must retain the above copyright
-#    notice, this list of conditions and the following disclaimer.
-#
-# 2. Redistributions in binary form must reproduce the above copyright
-#    notice, this list of conditions and the following disclaimer in the
-#    documentation and/or other materials provided with the distribution.
-#
-# 3. Neither the name of the California Institute of Technology nor
-#    the names of its contributors may be used to endorse or promote
-#    products derived from this software without specific prior
-#    written permission.
-#
-# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
-# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
-# FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CALTECH
-# OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
-# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
-# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
-# USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
-# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
-# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
-# OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
-# SUCH DAMAGE.
+# quantities for a linear system, as well as the main functions for
+# storing and plotting pole/zero and root locus diagrams.  (The actual
+# computation of root locus diagrams is in rlocus.py.)
 #
 
-from numpy import real, imag, linspace, exp, cos, sin, sqrt
+import itertools
+import warnings
 from math import pi
-from .lti import LTI
-from .iosys import isdtime, isctime
-from .grid import sgrid, zgrid, nogrid
+
+import matplotlib.pyplot as plt
+import numpy as np
+from numpy import cos, exp, imag, linspace, real, sin, sqrt
+
 from . import config
+from .freqplot import _freqplot_defaults, _get_line_labels
+from .grid import nogrid, sgrid, zgrid
+from .iosys import isctime, isdtime
+from .lti import LTI
+from .statesp import StateSpace
+from .xferfcn import TransferFunction
 
-__all__ = ['pzmap']
+__all__ = ['pole_zero_map', 'pole_zero_plot', 'pzmap', 'PoleZeroData']
 
 
 # Define default parameter values for this module
 _pzmap_defaults = {
-    'pzmap.grid': False,       # Plot omega-damping grid
-    'pzmap.plot': True,        # Generate plot using Matplotlib
+    'pzmap.grid': None,                 # Plot omega-damping grid
+    'pzmap.marker_size': 6,             # Size of the markers
+    'pzmap.marker_width': 1.5,          # Width of the markers
+    'pzmap.expansion_factor': 1.8,      # Amount to scale plots beyond features
+    'pzmap.buffer_factor': 1.05,        # Buffer to leave around plot peaks
 }
 
+#
+# Classes for keeping track of pzmap plots
+#
+# The PoleZeroData class keeps track of the information that is on a
+# pole/zero plot.
+#
+# In addition to the locations of poles and zeros, you can also save a set
+# of gains and loci for use in generating a root locus plot.  The gain
+# variable is a 1D array consisting of a list of increasing gains.  The
+# loci variable is a 2D array indexed by [gain_idx, root_idx] that can be
+# plotted using the `pole_zero_plot` function.
+#
+# The PoleZeroList class is used to return a list of pole/zero plots.  It
+# is a lightweight wrapper on the built-in list class that includes a
+# `plot` method, allowing plotting a set of root locus diagrams.
+#
+class PoleZeroData:
+    """Pole/zero data object.
+
+    This class is used as the return type for computing pole/zero responses
+    and root locus diagrams.  It contains information on the location of
+    system poles and zeros, as well as the gains and loci for root locus
+    diagrams.
+
+    Attributes
+    ----------
+    poles : ndarray
+        1D array of system poles.
+    zeros : ndarray
+        1D array of system zeros.
+    gains : ndarray, optional
+        1D array of gains for root locus plots.
+    loci : ndarray, optiona
+        2D array of poles, with each row corresponding to a gain.
+    sysname : str, optional
+        System name.
+    sys : StateSpace or TransferFunction
+        System corresponding to the data.
+
+    """
+    def __init__(
+            self, poles, zeros, gains=None, loci=None, dt=None, sysname=None,
+            sys=None):
+        """Create a pole/zero map object.
+
+        Parameters
+        ----------
+        poles : ndarray
+            1D array of system poles.
+        zeros : ndarray
+            1D array of system zeros.
+        gains : ndarray, optional
+            1D array of gains for root locus plots.
+        loci : ndarray, optiona
+            2D array of poles, with each row corresponding to a gain.
+        sysname : str, optional
+            System name.
+        sys : StateSpace or TransferFunction
+            System corresponding to the data.
+
+        """
+        self.poles = poles
+        self.zeros = zeros
+        self.gains = gains
+        self.loci = loci
+        self.dt = dt
+        self.sysname = sysname
+        self.sys = sys
+
+    # Implement functions to allow legacy assignment to tuple
+    def __iter__(self):
+        return iter((self.poles, self.zeros))
+
+    def plot(self, *args, **kwargs):
+        """Plot the pole/zero data.
+
+        See :func:`~control.pole_zero_plot` for description of arguments
+        and keywords.
+
+        """
+        # If this is a root locus plot, use rlocus defaults for grid
+        if self.loci is not None:
+            from .rlocus import _rlocus_defaults
+            kwargs = kwargs.copy()
+            kwargs['grid'] = config._get_param(
+                'rlocus', 'grid', kwargs.get('grid', None), _rlocus_defaults)
+
+        return pole_zero_plot(self, *args, **kwargs)
+
+
+class PoleZeroList(list):
+    """List of PoleZeroData objects."""
+    def plot(self, *args, **kwargs):
+        """Plot pole/zero data.
+
+        See :func:`~control.pole_zero_plot` for description of arguments
+        and keywords.
+
+        """
+        return pole_zero_plot(self, *args, **kwargs)
+
+
+# Pole/zero map
+def pole_zero_map(sysdata):
+    """Compute the pole/zero map for an LTI system.
+
+    Parameters
+    ----------
+    sys : LTI system (StateSpace or TransferFunction)
+        Linear system for which poles and zeros are computed.
+
+    Returns
+    -------
+    pzmap_data : PoleZeroMap
+        Pole/zero map containing the poles and zeros of the system.  Use
+        `pzmap_data.plot()` or `pole_zero_plot(pzmap_data)` to plot the
+        pole/zero map.
+
+    """
+    # Convert the first argument to a list
+    syslist = sysdata if isinstance(sysdata, (list, tuple)) else [sysdata]
+
+    responses = []
+    for idx, sys in enumerate(syslist):
+        responses.append(
+            PoleZeroData(
+                sys.poles(), sys.zeros(), dt=sys.dt, sysname=sys.name))
+
+    if isinstance(sysdata, (list, tuple)):
+        return PoleZeroList(responses)
+    else:
+        return responses[0]
+
 
 # TODO: Implement more elegant cross-style axes. See:
-#    http://matplotlib.sourceforge.net/examples/axes_grid/demo_axisline_style.html
-#    http://matplotlib.sourceforge.net/examples/axes_grid/demo_curvelinear_grid.html
-def pzmap(sys, plot=None, grid=None, title='Pole Zero Map', **kwargs):
+#    https://matplotlib.org/2.0.2/examples/axes_grid/demo_axisline_style.html
+#    https://matplotlib.org/2.0.2/examples/axes_grid/demo_curvelinear_grid.html
+def pole_zero_plot(
+        data, plot=None, grid=None, title=None, marker_color=None,
+        marker_size=None, marker_width=None, legend_loc='upper right',
+        xlim=None, ylim=None, interactive=None, ax=None, scaling=None,
+        initial_gain=None, **kwargs):
     """Plot a pole/zero map for a linear system.
 
+    If the system data include root loci, a root locus diagram for the
+    system is plotted.  When the root locus for a single system is plotted,
+    clicking on a location on the root locus will mark the gain on all
+    branches of the diagram and show the system gain and damping for the
+    given pole in the axes title.  Set to False to turn off this behavior.
+
     Parameters
     ----------
-    sys: LTI (StateSpace or TransferFunction)
-        Linear system for which poles and zeros are computed.
-    plot: bool, optional
-        If ``True`` a graph is generated with Matplotlib,
+    sysdata : List of PoleZeroData objects or LTI systems
+        List of pole/zero response data objects generated by pzmap_response()
+        or rootlocus_response() that are to be plotted.  If a list of systems
+        is given, the poles and zeros of those systems will be plotted.
+    grid : bool or str, optional
+        If `True` plot omega-damping grid, if `False` show imaginary axis
+        for continuous time systems, unit circle for discrete time systems.
+        If `empty`, do not draw any additonal lines.  Default value is set
+        by config.default['pzmap.grid'] or config.default['rlocus.grid'].
+    plot : bool, optional
+        (legacy) If ``True`` a graph is generated with Matplotlib,
         otherwise the poles and zeros are only computed and returned.
-    grid: boolean (default = False)
-        If True plot omega-damping grid.
+        If this argument is present, the legacy value of poles and
+        zeros is returned.
 
     Returns
     -------
-    poles: array
-        The systems poles
-    zeros: array
-        The system's zeros.
+    lines : array of list of Line2D
+        Array of Line2D objects for each set of markers in the plot. The
+        shape of the array is given by (nsys, 2) where nsys is the number
+        of systems or responses passed to the function.  The second index
+        specifies the pzmap object type:
+
+        * lines[idx, 0]: poles
+        * lines[idx, 1]: zeros
+
+    poles, zeros: list of arrays
+        (legacy) If the `plot` keyword is given, the system poles and zeros
+        are returned.
+
+    Other Parameters
+    ----------------
+    scaling : str or list, optional
+        Set the type of axis scaling.  Can be 'equal' (default), 'auto', or
+        a list of the form [xmin, xmax, ymin, ymax].
+    title : str, optional
+        Set the title of the plot.  Defaults plot type and system name(s).
+    marker_color : str, optional
+        Set the color of the markers used for poles and zeros.
+    marker_color : int, optional
+        Set the size of the markers used for poles and zeros.
+    marker_width : int, optional
+        Set the line width of the markers used for poles and zeros.
+    legend_loc : str, optional
+        For plots with multiple lines, a legend will be included in the
+        given location.  Default is 'center right'.  Use False to supress.
+    xlim : list, optional
+        Set the limits for the x axis.
+    ylim : list, optional
+        Set the limits for the y axis.
+    interactive : bool, optional
+        Turn off interactive mode for root locus plots.
+    initial_gain : float, optional
+        If given, the specified system gain will be marked on the plot.
 
     Notes
     -----
-    The pzmap function calls matplotlib.pyplot.axis('equal'), which means
-    that trying to reset the axis limits may not behave as expected.  To
-    change the axis limits, use matplotlib.pyplot.gca().axis('auto') and
-    then set the axis limits to the desired values.
+    By default, the pzmap function calls matplotlib.pyplot.axis('equal'),
+    which means that trying to reset the axis limits may not behave as
+    expected.  To change the axis limits, use the `scaling` keyword of use
+    matplotlib.pyplot.gca().axis('auto') and then set the axis limits to
+    the desired values.
 
     """
-    # Check to see if legacy 'Plot' keyword was used
-    if 'Plot' in kwargs:
-        import warnings
-        warnings.warn("'Plot' keyword is deprecated in pzmap; use 'plot'",
-                      FutureWarning)
-        plot = kwargs.pop('Plot')
+    # Get parameter values
+    grid = config._get_param('pzmap', 'grid', grid, _pzmap_defaults)
+    marker_size = config._get_param('pzmap', 'marker_size', marker_size, 6)
+    marker_width = config._get_param('pzmap', 'marker_width', marker_width, 1.5)
+    xlim_user, ylim_user = xlim, ylim
+    freqplot_rcParams = config._get_param(
+        'freqplot', 'rcParams', kwargs, _freqplot_defaults,
+        pop=True, last=True)
+    user_ax = ax
 
-    # Make sure there were no extraneous keywords
-    if kwargs:
-        raise TypeError("unrecognized keywords: ", str(kwargs))
+    # If argument was a singleton, turn it into a tuple
+    if not isinstance(data, (list, tuple)):
+        data = [data]
 
-    # Get parameter values
-    plot = config._get_param('pzmap', 'plot', plot, True)
-    grid = config._get_param('pzmap', 'grid', grid, False)
+    # If we are passed a list of systems, compute response first
+    if all([isinstance(
+            sys, (StateSpace, TransferFunction)) for sys in data]):
+        # Get the response, popping off keywords used there
+        pzmap_responses = pole_zero_map(data)
+    elif all([isinstance(d, PoleZeroData) for d in data]):
+        pzmap_responses = data
+    else:
+        raise TypeError("unknown system data type")
+
+    # Decide whether we are plotting any root loci
+    rlocus_plot = any([resp.loci is not None for resp in pzmap_responses])
+
+    # Turn on interactive mode by default, if allowed
+    if interactive is None and rlocus_plot and len(pzmap_responses) == 1 \
+       and pzmap_responses[0].sys is not None:
+        interactive = True
+
+    # Legacy return value processing
+    if plot is not None:
+        warnings.warn(
+            "`pole_zero_plot` return values of poles, zeros is deprecated; "
+            "use pole_zero_map()", DeprecationWarning)
 
-    if not isinstance(sys, LTI):
-        raise TypeError('Argument ``sys``: must be a linear system.')
+        # Extract out the values that we will eventually return
+        poles = [response.poles for response in pzmap_responses]
+        zeros = [response.zeros for response in pzmap_responses]
 
-    poles = sys.poles()
-    zeros = sys.zeros()
+    if plot is False:
+        if len(data) == 1:
+            return poles[0], zeros[0]
+        else:
+            return poles, zeros
+
+    # Initialize the figure
+    # TODO: turn into standard utility function (from plotutil.py?)
+    if user_ax is None:
+        fig = plt.gcf()
+        axs = fig.get_axes()
+    else:
+        fig = ax.figure
+        axs = [ax]
 
-    if (plot):
-        import matplotlib.pyplot as plt
+    if len(axs) > 1:
+        # Need to generate a new figure
+        fig, axs = plt.figure(), []
 
-        if grid:
-            if isdtime(sys, strict=True):
-                ax, fig = zgrid()
+    with plt.rc_context(freqplot_rcParams):
+        if grid and grid != 'empty':
+            plt.clf()
+            if all([isctime(dt=response.dt) for response in data]):
+                ax, fig = sgrid(scaling=scaling)
+            elif all([isdtime(dt=response.dt) for response in data]):
+                ax, fig = zgrid(scaling=scaling)
             else:
-                ax, fig = sgrid()
+                raise ValueError(
+                    "incompatible time bases; don't know how to grid")
+            # Store the limits for later use
+            xlim, ylim = ax.get_xlim(), ax.get_ylim()
+        elif len(axs) == 0:
+            if grid == 'empty':
+                # Leave off grid entirely
+                ax = plt.axes()
+                xlim = ylim = [0, 0]    # use data to set limits
+            else:
+                # draw stability boundary; use first response timebase
+                ax, fig = nogrid(data[0].dt, scaling=scaling)
+                xlim, ylim = ax.get_xlim(), ax.get_ylim()
+        else:
+            # Use the existing axes and any grid that is there
+            ax = axs[0]
+
+            # Store the limits for later use
+            xlim, ylim = ax.get_xlim(), ax.get_ylim()
+
+            # Issue a warning if the user tried to set the grid type
+            if grid:
+                warnings.warn("axis already exists; grid keyword ignored")
+
+    # Handle color cycle manually as all root locus segments
+    # of the same system are expected to be of the same color
+    color_cycle = plt.rcParams['axes.prop_cycle'].by_key()['color']
+    color_offset = 0
+    if len(ax.lines) > 0:
+        last_color = ax.lines[-1].get_color()
+        if last_color in color_cycle:
+            color_offset = color_cycle.index(last_color) + 1
+
+    # Create a list of lines for the output
+    out = np.empty(
+        (len(pzmap_responses), 3 if rlocus_plot else 2), dtype=object)
+    for i, j in itertools.product(range(out.shape[0]), range(out.shape[1])):
+        out[i, j] = []          # unique list in each element
+
+    # Plot the responses (and keep track of axes limits)
+    for idx, response in enumerate(pzmap_responses):
+        poles = response.poles
+        zeros = response.zeros
+
+        # Get the color to use for this system
+        if marker_color is None:
+            color = color_cycle[(color_offset + idx) % len(color_cycle)]
         else:
-            ax, fig = nogrid()
+            color = maker_color
 
         # Plot the locations of the poles and zeros
         if len(poles) > 0:
-            ax.scatter(real(poles), imag(poles), s=50, marker='x',
-                       facecolors='k')
+            label = response.sysname if response.loci is None else None
+            out[idx, 0] = ax.plot(
+                real(poles), imag(poles), marker='x', linestyle='',
+                markeredgecolor=color, markerfacecolor=color,
+                markersize=marker_size, markeredgewidth=marker_width,
+                label=label)
         if len(zeros) > 0:
-            ax.scatter(real(zeros), imag(zeros), s=50, marker='o',
-                       facecolors='none', edgecolors='k')
+            out[idx, 1] = ax.plot(
+                real(zeros), imag(zeros), marker='o', linestyle='',
+                markeredgecolor=color, markerfacecolor='none',
+                markersize=marker_size, markeredgewidth=marker_width)
+
+        # Plot the loci, if present
+        if response.loci is not None:
+            for locus in response.loci.transpose():
+                out[idx, 2] += ax.plot(
+                    real(locus), imag(locus), color=color,
+                    label=response.sysname)
+
+            # Compute the axis limits to use based on the response
+            resp_xlim, resp_ylim = _compute_root_locus_limits(response)
+
+            # Keep track of the current limits
+            xlim = [min(xlim[0], resp_xlim[0]), max(xlim[1], resp_xlim[1])]
+            ylim = [min(ylim[0], resp_ylim[0]), max(ylim[1], resp_ylim[1])]
+
+            # Plot the initial gain, if given
+            if initial_gain is not None:
+                _mark_root_locus_gain(ax, response.sys, initial_gain)
+
+            # TODO: add arrows to root loci (reuse Nyquist arrow code?)
+
+    # Set the axis limits to something reasonable
+    if rlocus_plot:
+        # Set up the limits for the plot using information from loci
+        ax.set_xlim(xlim if xlim_user is None else xlim_user)
+        ax.set_ylim(ylim if ylim_user is None else ylim_user)
+    else:
+        # No root loci => only set axis limits if users specified them
+        if xlim_user is not None:
+            ax.set_xlim(xlim_user)
+        if ylim_user is not None:
+            ax.set_ylim(ylim_user)
+
+    # List of systems that are included in this plot
+    lines, labels = _get_line_labels(ax)
+
+    # Add legend if there is more than one system plotted
+    if len(labels) > 1 and legend_loc is not False:
+        if response.loci is None:
+            # Use "x o" for the system label, via matplotlib tuple handler
+            from matplotlib.legend_handler import HandlerTuple
+            from matplotlib.lines import Line2D
+
+            line_tuples = []
+            for pole_line in lines:
+                zero_line = Line2D(
+                    [0], [0], marker='o', linestyle='',
+                    markeredgecolor=pole_line.get_markerfacecolor(),
+                    markerfacecolor='none', markersize=marker_size,
+                    markeredgewidth=marker_width)
+            handle = (pole_line, zero_line)
+            line_tuples.append(handle)
+
+            with plt.rc_context(freqplot_rcParams):
+                ax.legend(
+                    line_tuples, labels, loc=legend_loc,
+                    handler_map={tuple: HandlerTuple(ndivide=None)})
+        else:
+            # Regular legend, with lines
+            with plt.rc_context(freqplot_rcParams):
+                ax.legend(lines, labels, loc=legend_loc)
+
+    # Add the title
+    if title is None:
+        title = "Pole/zero plot for " + ", ".join(labels)
+    if user_ax is None:
+        with plt.rc_context(freqplot_rcParams):
+            fig.suptitle(title)
+
+    # Add dispather to handle choosing a point on the diagram
+    if interactive:
+        if len(pzmap_responses) > 1:
+            raise NotImplementedError(
+                "interactive mode only allowed for single system")
+        elif pzmap_responses[0].sys == None:
+            raise SystemError("missing system information")
+        else:
+            sys = pzmap_responses[0].sys
+
+        # Define function to handle mouse clicks
+        def _click_dispatcher(event):
+            # Find the gain corresponding to the clicked point
+            K, s = _find_root_locus_gain(event, sys, ax)
+
+            if K is not None:
+                # Mark the gain on the root locus diagram
+                _mark_root_locus_gain(ax, sys, K)
+
+                # Display the parameters in the axes title
+                with plt.rc_context(freqplot_rcParams):
+                    ax.set_title(_create_root_locus_label(sys, K, s))
+
+            ax.figure.canvas.draw()
+
+        fig.canvas.mpl_connect('button_release_event', _click_dispatcher)
+
+    # Legacy processing: return locations of poles and zeros as a tuple
+    if plot is True:
+        if len(data) == 1:
+            return poles, zeros
+        else:
+            TypeError("system lists not supported with legacy return values")
+
+    return out
+
+
+# Utility function to find gain corresponding to a click event
+def _find_root_locus_gain(event, sys, ax):
+    # Get the current axis limits to set various thresholds
+    xlim, ylim = ax.get_xlim(), ax.get_ylim()
+
+    # Catch type error when event click is in the figure but not in an axis
+    try:
+        s = complex(event.xdata, event.ydata)
+        K = -1. / sys(s)
+        K_xlim = -1. / sys(
+            complex(event.xdata + 0.05 * abs(xlim[1] - xlim[0]), event.ydata))
+        K_ylim = -1. / sys(
+            complex(event.xdata, event.ydata + 0.05 * abs(ylim[1] - ylim[0])))
+
+    except TypeError:
+        K = float('inf')
+        K_xlim = float('inf')
+        K_ylim = float('inf')
+
+    #
+    # Compute tolerances for deciding if we clicked on the root locus
+    #
+    # This is a bit of black magic that sets some limits for how close we
+    # need to be to the root locus in order to consider it a click on the
+    # actual curve.  Otherwise, we will just ignore the click.
+
+    x_tolerance = 0.1 * abs((xlim[1] - xlim[0]))
+    y_tolerance = 0.1 * abs((ylim[1] - ylim[0]))
+    gain_tolerance = np.mean([x_tolerance, y_tolerance]) * 0.1 + \
+        0.1 * max([abs(K_ylim.imag/K_ylim.real), abs(K_xlim.imag/K_xlim.real)])
+
+    # Decide whether to pay attention to this event
+    if abs(K.real) > 1e-8 and abs(K.imag / K.real) < gain_tolerance and \
+       event.inaxes == ax.axes and K.real > 0.:
+        return K.real, s
+
+    else:
+        return None, s
+
+
+# Mark points corresponding to a given gain on root locus plot
+def _mark_root_locus_gain(ax, sys, K):
+    from .rlocus import _RLFindRoots, _systopoly1d
+
+    # Remove any previous gain points
+    for line in reversed(ax.lines):
+        if line.get_label() == '_gain_point':
+            line.remove()
+            del line
+
+    # Visualise clicked point, displaying all roots
+    # TODO: allow marker parameters to be set
+    nump, denp = _systopoly1d(sys)
+    root_array = _RLFindRoots(nump, denp, K.real)
+    ax.plot(
+        [root.real for root in root_array], [root.imag for root in root_array],
+        marker='s', markersize=6, zorder=20, label='_gain_point', color='k')
+
+
+# Return a string identifying a clicked point
+def _create_root_locus_label(sys, K, s):
+    # Figure out the damping ratio
+    if isdtime(sys, strict=True):
+        zeta = -np.cos(np.angle(np.log(s)))
+    else:
+        zeta = -1 * s.real / abs(s)
+
+    return "Clicked at: %.4g%+.4gj   gain = %.4g  damping = %.4g" % \
+        (s.real, s.imag, K.real, zeta)
+
+
+# Utility function to compute limits for root loci
+def _compute_root_locus_limits(response):
+    loci = response.loci
+
+    # Start with information about zeros, if present
+    if response.sys is not None and response.sys.zeros().size > 0:
+        xlim = [
+            min(0, np.min(response.sys.zeros().real)),
+            max(0, np.max(response.sys.zeros().real))
+        ]
+        ylim = max(0, np.max(response.sys.zeros().imag))
+    else:
+        xlim, ylim = [0, 0], 0
+
+    # Go through each locus and look for features
+    rho = config._get_param('pzmap', 'buffer_factor')
+    for locus in loci.transpose():
+        # Include all starting points
+        xlim = [min(xlim[0], locus[0].real), max(xlim[1], locus[0].real)]
+        ylim = max(ylim, locus[0].imag)
+
+        # Find the local maxima of root locus curve
+        xpeaks = np.where(
+            np.diff(np.abs(locus.real)) < 0, locus.real[0:-1], 0)
+        xlim = [
+            min(xlim[0], np.min(xpeaks) * rho),
+            max(xlim[1], np.max(xpeaks) * rho)
+        ]
+
+        ypeaks = np.where(
+            np.diff(np.abs(locus.imag)) < 0, locus.imag[0:-1], 0)
+        ylim = max(ylim, np.max(ypeaks) * rho)
+
+    if isctime(dt=response.dt):
+        # Adjust the limits to include some space around features
+        # TODO: use _k_max and project out to max k for all value?
+        rho = config._get_param('pzmap', 'expansion_factor')
+        xlim[0] = rho * xlim[0] if xlim[0] < 0 else 0
+        xlim[1] = rho * xlim[1] if xlim[1] > 0 else 0
+        ylim = rho * ylim if ylim > 0 else np.max(np.abs(xlim))
+
+    return xlim, [-ylim, ylim]
 
-        plt.title(title)
 
-    # Return locations of poles and zeros as a tuple
-    return poles, zeros
+pzmap = pole_zero_plot
diff --git a/control/rlocus.py b/control/rlocus.py
index 41cdec05..ea17ae94 100644
--- a/control/rlocus.py
+++ b/control/rlocus.py
@@ -1,38 +1,6 @@
 # rlocus.py - code for computing a root locus plot
 # Code contributed by Ryan Krauss, 2010
 #
-# Copyright (c) 2010 by Ryan Krauss
-# All rights reserved.
-#
-# Redistribution and use in source and binary forms, with or without
-# modification, are permitted provided that the following conditions
-# are met:
-#
-# 1. Redistributions of source code must retain the above copyright
-#    notice, this list of conditions and the following disclaimer.
-#
-# 2. Redistributions in binary form must reproduce the above copyright
-#    notice, this list of conditions and the following disclaimer in the
-#    documentation and/or other materials provided with the distribution.
-#
-# 3. Neither the name of the California Institute of Technology nor
-#    the names of its contributors may be used to endorse or promote
-#    products derived from this software without specific prior
-#    written permission.
-#
-# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
-# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
-# FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CALTECH
-# OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
-# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
-# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
-# USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
-# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
-# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
-# OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
-# SUCH DAMAGE.
-#
 # RMM, 17 June 2010: modified to be a standalone piece of code
 #   * Added BSD copyright info to file (per Ryan)
 #   * Added code to convert (num, den) to poly1d's if they aren't already.
@@ -46,48 +14,101 @@
 # Sawyer B. Fuller (minster@uw.edu) 21 May 2020:
 #   * added compatibility with discrete-time systems.
 #
-# $Id$
 
-# Packages used by this module
+import warnings
 from functools import partial
-import numpy as np
-import matplotlib as mpl
+
 import matplotlib.pyplot as plt
-from numpy import array, poly1d, row_stack, zeros_like, real, imag
-import scipy.signal             # signal processing toolbox
+import numpy as np
+import scipy.signal  # signal processing toolbox
+from numpy import array, imag, poly1d, real, row_stack, zeros_like
+
+from . import config
+from .exception import ControlMIMONotImplemented
 from .iosys import isdtime
+from .lti import LTI
 from .xferfcn import _convert_to_transfer_function
-from .exception import ControlMIMONotImplemented
-from .sisotool import _SisotoolUpdate
-from .grid import sgrid, zgrid
-from . import config
-import warnings
 
-__all__ = ['root_locus', 'rlocus']
+__all__ = ['root_locus_map', 'root_locus_plot', 'root_locus', 'rlocus']
 
 # Default values for module parameters
 _rlocus_defaults = {
     'rlocus.grid': True,
-    'rlocus.plotstr': 'b' if int(mpl.__version__[0]) == 1 else 'C0',
-    'rlocus.print_gain': True,
-    'rlocus.plot': True
 }
 
 
-# Main function: compute a root locus diagram
-def root_locus(sys, kvect=None, xlim=None, ylim=None,
-               plotstr=None, plot=True, print_gain=None, grid=None, ax=None,
-               initial_gain=None, **kwargs):
+# Root locus map
+def root_locus_map(sysdata, gains=None):
+    """Compute the root locus map for an LTI system.
+
+    Calculate the root locus by finding the roots of 1 + k * G(s) where G
+    is a linear system with transfer function num(s)/den(s) and each k is
+    an element of kvect.
+
+    Parameters
+    ----------
+    sys : LTI system or list of LTI systems
+        Linear input/output systems (SISO only, for now).
+    kvect : array_like, optional
+        Gains to use in computing plot of closed-loop poles.
+
+    Returns
+    -------
+    rldata : PoleZeroData or list of PoleZeroData
+        Root locus data object(s) corresponding to the .  The loci of
+        the root locus diagram are available in the array
+        `rldata.loci`, indexed by the gain index and the locus index,
+        and the gains are in the array `rldata.gains`.
+
+    Notes
+    -----
+    For backward compatibility, the `rldata` return object can be
+    assigned to the tuple `roots, gains`.
+
+    """
+    from .pzmap import PoleZeroData, PoleZeroList
+
+    # Convert the first argument to a list
+    syslist = sysdata if isinstance(sysdata, (list, tuple)) else [sysdata]
+
+    responses = []
+    for idx, sys in enumerate(syslist):
+        if not sys.issiso():
+            raise ControlMIMONotImplemented(
+                "sys must be single-input single-output (SISO)")
+
+        # Convert numerator and denominator to polynomials if they aren't
+        nump, denp = _systopoly1d(sys[0, 0])
+
+        if gains is None:
+            kvect, root_array, _, _ = _default_gains(nump, denp, None, None)
+        else:
+            kvect = np.atleast_1d(gains)
+            root_array = _RLFindRoots(nump, denp, kvect)
+            root_array = _RLSortRoots(root_array)
+
+        responses.append(PoleZeroData(
+            sys.poles(), sys.zeros(), kvect, root_array,
+            dt=sys.dt, sysname=sys.name, sys=sys))
+
+    if isinstance(sysdata, (list, tuple)):
+        return PoleZeroList(responses)
+    else:
+        return responses[0]
+
+
+def root_locus_plot(
+        sysdata, kvect=None, grid=None, plot=None, **kwargs):
 
     """Root locus plot.
 
-    Calculate the root locus by finding the roots of 1+k*TF(s)
-    where TF is self.num(s)/self.den(s) and each k is an element
-    of kvect.
+    Calculate the root locus by finding the roots of 1 + k * G(s) where G
+    is a linear system with transfer function num(s)/den(s) and each k is
+    an element of kvect.
 
     Parameters
     ----------
-    sys : LTI object
+    sysdata : PoleZeroMap or LTI object or list
         Linear input/output systems (SISO only, for now).
     kvect : array_like, optional
         Gains to use in computing plot of closed-loop poles.
@@ -97,182 +118,83 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None,
     ylim : tuple or list, optional
         Set limits of y axis, normally with tuple
         (see :doc:`matplotlib:api/axes_api`).
-    plotstr : :func:`matplotlib.pyplot.plot` format string, optional
-        plotting style specification
-    plot : boolean, optional
-        If True (default), plot root locus diagram.
-    print_gain : bool
-        If True (default), report mouse clicks when close to the root locus
-        branches, calculate gain, damping and print.
-    grid : bool
-        If True plot omega-damping grid.  Default is False.
+    plot : bool, optional
+        (legacy) If given, `root_locus_plot` returns the legacy return values
+        of roots and gains.  If False, just return the values with no plot.
+    grid : bool or str, optional
+        If `True` plot omega-damping grid, if `False` show imaginary axis
+        for continuous time systems, unit circle for discrete time systems.
+        If `empty`, do not draw any additonal lines.  Default value is set
+        by config.default['rlocus.grid'].
     ax : :class:`matplotlib.axes.Axes`
         Axes on which to create root locus plot
     initial_gain : float, optional
-        Used by :func:`sisotool` to indicate initial gain.
+        Mark the point on the root locus diagram corresponding to the
+        given gain.
 
     Returns
     -------
-    roots : ndarray
-        Closed-loop root locations, arranged in which each row corresponds
-        to a gain in gains
-    gains : ndarray
-        Gains used.  Same as kvect keyword argument if provided.
+    lines : array of list of Line2D
+        Array of Line2D objects for each set of markers in the plot. The
+        shape of the array is given by (nsys, 3) where nsys is the number
+        of systems or responses passed to the function.  The second index
+        specifies the object type:
+
+        * lines[idx, 0]: poles
+        * lines[idx, 1]: zeros
+        * lines[idx, 2]: loci
+
+    roots, gains : ndarray
+        (legacy) If the `plot` keyword is given, returns the
+        closed-loop root locations, arranged such that each row
+        corresponds to a gain in gains, and the array of gains (ame as
+        kvect keyword argument if provided).
 
     Notes
     -----
-    The root_locus function calls matplotlib.pyplot.axis('equal'), which
+    The root_locus_plot function calls matplotlib.pyplot.axis('equal'), which
     means that trying to reset the axis limits may not behave as expected.
     To change the axis limits, use matplotlib.pyplot.gca().axis('auto') and
     then set the axis limits to the desired values.
 
     """
-    # Check to see if legacy 'Plot' keyword was used
-    if 'Plot' in kwargs:
-        warnings.warn("'Plot' keyword is deprecated in root_locus; "
-                      "use 'plot'", FutureWarning)
-        # Map 'Plot' keyword to 'plot' keyword
-        plot = kwargs.pop('Plot')
-
-    # Check to see if legacy 'PrintGain' keyword was used
-    if 'PrintGain' in kwargs:
-        warnings.warn("'PrintGain' keyword is deprecated in root_locus; "
-                      "use 'print_gain'", FutureWarning)
-        # Map 'PrintGain' keyword to 'print_gain' keyword
-        print_gain = kwargs.pop('PrintGain')
-
-    # Get parameter values
-    plotstr = config._get_param('rlocus', 'plotstr', plotstr, _rlocus_defaults)
-    grid = config._get_param('rlocus', 'grid', grid, _rlocus_defaults)
-    print_gain = config._get_param(
-        'rlocus', 'print_gain', print_gain, _rlocus_defaults)
+    from .pzmap import pole_zero_plot
 
-    # Check for sisotool mode
-    sisotool = kwargs.get('sisotool', False)
-
-    # make sure siso. sisotool has different requirements
-    if not sys.issiso() and not sisotool:
-        raise ControlMIMONotImplemented(
-            'sys must be single-input single-output (SISO)')
+    # Set default parameters
+    grid = config._get_param('rlocus', 'grid', grid, _rlocus_defaults)
 
-    sys_loop = sys[0,0]
-    # Convert numerator and denominator to polynomials if they aren't
-    (nump, denp) = _systopoly1d(sys_loop)
-
-    # if discrete-time system and if xlim and ylim are not given,
-    #  that we a view of the unit circle
-    if xlim is None and isdtime(sys, strict=True):
-        xlim = (-1.2, 1.2)
-    if ylim is None and isdtime(sys, strict=True):
-        xlim = (-1.3, 1.3)
-
-    if kvect is None:
-        kvect, root_array, xlim, ylim = _default_gains(nump, denp, xlim, ylim)
-        recompute_on_zoom = True
+    if isinstance(sysdata, list) and all(
+            [isinstance(sys, LTI) for sys in sysdata]) or \
+            isinstance(sysdata, LTI):
+        responses = root_locus_map(sysdata, gains=kvect)
     else:
-        kvect = np.atleast_1d(kvect)
-        root_array = _RLFindRoots(nump, denp, kvect)
-        root_array = _RLSortRoots(root_array)
-        recompute_on_zoom = False
+        responses = sysdata
 
-    if sisotool:
-        start_roots = _RLFindRoots(nump, denp, initial_gain)
+    #
+    # Process `plot` keyword
+    #
+    # See bode_plot for a description of how this keyword is handled to
+    # support legacy implementatoins of root_locus.
+    #
+    if plot is not None:
+        warnings.warn(
+            "`root_locus` return values of roots, gains is deprecated; "
+            "use root_locus_map()", DeprecationWarning)
 
-    # Make sure there were no extraneous keywords
-    if not sisotool and kwargs:
-        raise TypeError("unrecognized keywords: ", str(kwargs))
+    if plot is False:
+        return responses.loci, responses.gains
 
-    # Create the Plot
-    if plot:
-        if sisotool:
-            fig = kwargs['fig']
-            ax = fig.axes[1]
-        else:
-            if ax is None:
-                ax = plt.gca()
-            fig = ax.figure
-            ax.set_title('Root Locus')
-
-        if print_gain and not sisotool:
-            fig.canvas.mpl_connect(
-                'button_release_event',
-                partial(_RLClickDispatcher, sys=sys, fig=fig,
-                        ax_rlocus=fig.axes[0], plotstr=plotstr))
-        elif sisotool:
-            fig.axes[1].plot(
-                [root.real for root in start_roots],
-                [root.imag for root in start_roots],
-                marker='s', markersize=6, zorder=20, color='k', label='gain_point')
-            s = start_roots[0][0]
-            if isdtime(sys, strict=True):
-                zeta = -np.cos(np.angle(np.log(s)))
-            else:
-                zeta = -1 * s.real / abs(s)
-            fig.suptitle(
-                "Clicked at: %10.4g%+10.4gj  gain: %10.4g  damp: %10.4g" %
-                (s.real, s.imag, initial_gain, zeta),
-                fontsize=12 if int(mpl.__version__[0]) == 1 else 10)
-            fig.canvas.mpl_connect(
-                'button_release_event',
-                partial(_RLClickDispatcher, sys=sys, fig=fig,
-                        ax_rlocus=fig.axes[1], plotstr=plotstr,
-                        sisotool=sisotool,
-                        bode_plot_params=kwargs['bode_plot_params'],
-                        tvect=kwargs['tvect']))
-
-
-        if recompute_on_zoom:
-            # update gains and roots when xlim/ylim change. Only then are
-            # data on available. I.e., cannot combine with _RLClickDispatcher
-            dpfun = partial(
-                _RLZoomDispatcher, sys=sys, ax_rlocus=ax, plotstr=plotstr)
-            # TODO: the next too lines seem to take a long time to execute
-            # TODO: is there a way to speed them up?  (RMM, 6 Jun 2019)
-            ax.callbacks.connect('xlim_changed', dpfun)
-            ax.callbacks.connect('ylim_changed', dpfun)
-
-        # plot open loop poles
-        poles = array(denp.r)
-        ax.plot(real(poles), imag(poles), 'x')
-
-        # plot open loop zeros
-        zeros = array(nump.r)
-        if zeros.size > 0:
-            ax.plot(real(zeros), imag(zeros), 'o')
-
-        # Now plot the loci
-        for index, col in enumerate(root_array.T):
-            ax.plot(real(col), imag(col), plotstr, label='rootlocus')
-
-        # Set up plot axes and labels
-        ax.set_xlabel('Real')
-        ax.set_ylabel('Imaginary')
-
-        # Set up the limits for the plot
-        # Note: need to do this before computing grid lines
-        if xlim:
-            ax.set_xlim(xlim)
-        if ylim:
-            ax.set_ylim(ylim)
-
-        # Draw the grid
-        if grid:
-            if isdtime(sys, strict=True):
-                zgrid(ax=ax)
-            else:
-                _sgrid_func(ax)
-        else:
-            ax.axhline(0., linestyle=':', color='k', linewidth=.75, zorder=-20)
-            ax.axvline(0., linestyle=':', color='k', linewidth=.75, zorder=-20)
-            if isdtime(sys, strict=True):
-                ax.add_patch(plt.Circle(
-                    (0, 0), radius=1.0, linestyle=':', edgecolor='k',
-                    linewidth=0.75, fill=False, zorder=-20))
+    # Plot the root loci
+    out = responses.plot(grid=grid, **kwargs)
+
+    # Legacy processing: return locations of poles and zeros as a tuple
+    if plot is True:
+        return responses.loci, responses.gains
 
-    return root_array, kvect
+    return out
 
 
-def _default_gains(num, den, xlim, ylim, zoom_xlim=None, zoom_ylim=None):
+def _default_gains(num, den, xlim, ylim):
     """Unsupervised gains calculation for root locus plot.
 
     References
@@ -281,16 +203,23 @@ def _default_gains(num, den, xlim, ylim, zoom_xlim=None, zoom_ylim=None):
     Saddle River, NJ : New Delhi: Prentice Hall..
 
     """
+    # Compute the break points on the real axis for the root locus plot
     k_break, real_break = _break_points(num, den)
+
+    # Decide on the maximum gain to use and create the gain vector
     kmax = _k_max(num, den, real_break, k_break)
     kvect = np.hstack((np.linspace(0, kmax, 50), np.real(k_break)))
     kvect.sort()
 
+    # Find the roots for all of the gains and sort them
     root_array = _RLFindRoots(num, den, kvect)
     root_array = _RLSortRoots(root_array)
+
+    # Keep track of the open loop poles and zeros
     open_loop_poles = den.roots
     open_loop_zeros = num.roots
 
+    # ???
     if open_loop_zeros.size != 0 and \
        open_loop_zeros.size < open_loop_poles.size:
         open_loop_zeros_xl = np.append(
@@ -345,7 +274,7 @@ def _default_gains(num, den, xlim, ylim, zoom_xlim=None, zoom_ylim=None):
         tolerance = x_tolerance
     else:
         tolerance = np.min([x_tolerance, y_tolerance])
-    indexes_too_far = _indexes_filt(root_array, tolerance, zoom_xlim, zoom_ylim)
+    indexes_too_far = _indexes_filt(root_array, tolerance)
 
     # Add more points into the root locus for points that are too far apart
     while len(indexes_too_far) > 0 and kvect.size < 5000:
@@ -357,7 +286,7 @@ def _default_gains(num, den, xlim, ylim, zoom_xlim=None, zoom_ylim=None):
             root_array = np.insert(root_array, index + 1, new_points, axis=0)
 
         root_array = _RLSortRoots(root_array)
-        indexes_too_far = _indexes_filt(root_array, tolerance, zoom_xlim, zoom_ylim)
+        indexes_too_far = _indexes_filt(root_array, tolerance)
 
     new_gains = kvect[-1] * np.hstack((np.logspace(0, 3, 4)))
     new_points = _RLFindRoots(num, den, new_gains[1:4])
@@ -367,8 +296,8 @@ def _default_gains(num, den, xlim, ylim, zoom_xlim=None, zoom_ylim=None):
     return kvect, root_array, xlim, ylim
 
 
-def _indexes_filt(root_array, tolerance, zoom_xlim=None, zoom_ylim=None):
-    """Calculate the distance between points and return the indexes.
+def _indexes_filt(root_array, tolerance):
+    """Calculate the distance between points and return the indices.
 
     Filter the indexes so only the resolution of points within the xlim and
     ylim is improved when zoom is used.
@@ -376,48 +305,6 @@ def _indexes_filt(root_array, tolerance, zoom_xlim=None, zoom_ylim=None):
     """
     distance_points = np.abs(np.diff(root_array, axis=0))
     indexes_too_far = list(np.unique(np.where(distance_points > tolerance)[0]))
-
-    if zoom_xlim is not None and zoom_ylim is not None:
-        x_tolerance_zoom = 0.05 * (zoom_xlim[1] - zoom_xlim[0])
-        y_tolerance_zoom = 0.05 * (zoom_ylim[1] - zoom_ylim[0])
-        tolerance_zoom = np.min([x_tolerance_zoom, y_tolerance_zoom])
-        indexes_too_far_zoom = list(
-            np.unique(np.where(distance_points > tolerance_zoom)[0]))
-        indexes_too_far_filtered = []
-
-        for index in indexes_too_far_zoom:
-            for point in root_array[index]:
-                if (zoom_xlim[0] <= point.real <= zoom_xlim[1]) and \
-                   (zoom_ylim[0] <= point.imag <= zoom_ylim[1]):
-                    indexes_too_far_filtered.append(index)
-                    break
-
-        # Check if zoom box is not overshot & insert points where neccessary
-        if len(indexes_too_far_filtered) == 0 and len(root_array) < 500:
-            limits = [zoom_xlim[0], zoom_xlim[1], zoom_ylim[0], zoom_ylim[1]]
-            for index, limit in enumerate(limits):
-                if index <= 1:
-                    asign = np.sign(real(root_array)-limit)
-                else:
-                    asign = np.sign(imag(root_array) - limit)
-                signchange = ((np.roll(asign, 1, axis=0)
-                               - asign) != 0).astype(int)
-                signchange[0] = np.zeros((len(root_array[0])))
-                if len(np.where(signchange == 1)[0]) > 0:
-                    indexes_too_far_filtered.append(
-                        np.where(signchange == 1)[0][0]-1)
-
-        if len(indexes_too_far_filtered) > 0:
-            if indexes_too_far_filtered[0] != 0:
-                indexes_too_far_filtered.insert(
-                    0, indexes_too_far_filtered[0]-1)
-            if not indexes_too_far_filtered[-1] + 1 >= len(root_array) - 2:
-                indexes_too_far_filtered.append(
-                    indexes_too_far_filtered[-1] + 1)
-
-        indexes_too_far.extend(indexes_too_far_filtered)
-
-    indexes_too_far = list(np.unique(indexes_too_far))
     indexes_too_far.sort()
     return indexes_too_far
 
@@ -558,249 +445,6 @@ def _RLSortRoots(roots):
     return sorted
 
 
-def _RLZoomDispatcher(event, sys, ax_rlocus, plotstr):
-    """Rootlocus plot zoom dispatcher"""
-    sys_loop = sys[0,0]
-    nump, denp = _systopoly1d(sys_loop)
-    xlim, ylim = ax_rlocus.get_xlim(), ax_rlocus.get_ylim()
-
-    kvect, root_array, xlim, ylim = _default_gains(
-        nump, denp, xlim=None, ylim=None, zoom_xlim=xlim, zoom_ylim=ylim)
-    _removeLine('rootlocus', ax_rlocus)
-
-    for i, col in enumerate(root_array.T):
-        ax_rlocus.plot(real(col), imag(col), plotstr, label='rootlocus',
-                       scalex=False, scaley=False)
-
-
-def _RLClickDispatcher(event, sys, fig, ax_rlocus, plotstr, sisotool=False,
-                       bode_plot_params=None, tvect=None):
-    """Rootlocus plot click dispatcher"""
-
-    # Zoom is handled by specialized callback above, only do gain plot
-    if event.inaxes == ax_rlocus.axes and \
-       plt.get_current_fig_manager().toolbar.mode not in \
-       {'zoom rect', 'pan/zoom'}:
-
-        # if a point is clicked on the rootlocus plot visually emphasize it
-        K = _RLFeedbackClicksPoint(event, sys, fig, ax_rlocus, sisotool)
-        if sisotool and K is not None:
-            _SisotoolUpdate(sys, fig, K, bode_plot_params, tvect)
-
-    # Update the canvas
-    fig.canvas.draw()
-
-
-def _RLFeedbackClicksPoint(event, sys, fig, ax_rlocus, sisotool=False):
-    """Display root-locus gain feedback point for clicks on root-locus plot"""
-    sys_loop = sys[0,0]
-    (nump, denp) = _systopoly1d(sys_loop)
-
-    xlim = ax_rlocus.get_xlim()
-    ylim = ax_rlocus.get_ylim()
-    x_tolerance = 0.1 * abs((xlim[1] - xlim[0]))
-    y_tolerance = 0.1 * abs((ylim[1] - ylim[0]))
-    gain_tolerance = np.mean([x_tolerance, y_tolerance])*0.1
-
-    # Catch type error when event click is in the figure but not in an axis
-    try:
-        s = complex(event.xdata, event.ydata)
-        K = -1. / sys_loop(s)
-        K_xlim = -1. / sys_loop(
-            complex(event.xdata + 0.05 * abs(xlim[1] - xlim[0]), event.ydata))
-        K_ylim = -1. / sys_loop(
-            complex(event.xdata, event.ydata + 0.05 * abs(ylim[1] - ylim[0])))
-
-    except TypeError:
-        K = float('inf')
-        K_xlim = float('inf')
-        K_ylim = float('inf')
-
-    gain_tolerance += 0.1 * max([abs(K_ylim.imag/K_ylim.real),
-                                 abs(K_xlim.imag/K_xlim.real)])
-
-    if abs(K.real) > 1e-8 and abs(K.imag / K.real) < gain_tolerance and \
-       event.inaxes == ax_rlocus.axes and K.real > 0.:
-
-        if isdtime(sys, strict=True):
-            zeta = -np.cos(np.angle(np.log(s)))
-        else:
-            zeta = -1 * s.real / abs(s)
-
-        # Display the parameters in the output window and figure
-        print("Clicked at %10.4g%+10.4gj gain %10.4g damp %10.4g" %
-              (s.real, s.imag, K.real, zeta))
-        fig.suptitle(
-            "Clicked at: %10.4g%+10.4gj  gain: %10.4g  damp: %10.4g" %
-            (s.real, s.imag, K.real, zeta),
-            fontsize=12 if int(mpl.__version__[0]) == 1 else 10)
-
-        # Remove the previous line
-        _removeLine(label='gain_point', ax=ax_rlocus)
-
-        # Visualise clicked point, display all roots for sisotool mode
-        if sisotool:
-            root_array = _RLFindRoots(nump, denp, K.real)
-            ax_rlocus.plot(
-                [root.real for root in root_array],
-                [root.imag for root in root_array],
-                marker='s', markersize=6, zorder=20, label='gain_point', color='k')
-        else:
-            ax_rlocus.plot(s.real, s.imag, 'k.', marker='s', markersize=8,
-                           zorder=20, label='gain_point')
-
-        return K.real
-
-
-def _removeLine(label, ax):
-    """Remove a line from the ax when a label is specified"""
-    for line in reversed(ax.lines):
-        if line.get_label() == label:
-            line.remove()
-            del line
-
-
-def _sgrid_func(ax, zeta=None, wn=None):
-    # Get locator function for x-axis, y-axis tick marks
-    xlocator = ax.get_xaxis().get_major_locator()
-    ylocator = ax.get_yaxis().get_major_locator()
-
-    # Decide on the location for the labels (?)
-    ylim = ax.get_ylim()
-    ytext_pos_lim = ylim[1] - (ylim[1] - ylim[0]) * 0.03
-    xlim = ax.get_xlim()
-    xtext_pos_lim = xlim[0] + (xlim[1] - xlim[0]) * 0.0
-
-    # Create a list of damping ratios, if needed
-    if zeta is None:
-        zeta = _default_zetas(xlim, ylim)
-
-    # Figure out the angles for the different damping ratios
-    angles = []
-    for z in zeta:
-        if (z >= 1e-4) and (z <= 1):
-            angles.append(np.pi/2 + np.arcsin(z))
-        else:
-            zeta.remove(z)
-    y_over_x = np.tan(angles)
-
-    # zeta-constant lines
-    for index, yp in enumerate(y_over_x):
-        ax.plot([0, xlocator()[0]], [0, yp * xlocator()[0]], color='gray',
-                linestyle='dashed', linewidth=0.5)
-        ax.plot([0, xlocator()[0]], [0, -yp * xlocator()[0]], color='gray',
-                linestyle='dashed', linewidth=0.5)
-        an = "%.2f" % zeta[index]
-        if yp < 0:
-            xtext_pos = 1/yp * ylim[1]
-            ytext_pos = yp * xtext_pos_lim
-            if np.abs(xtext_pos) > np.abs(xtext_pos_lim):
-                xtext_pos = xtext_pos_lim
-            else:
-                ytext_pos = ytext_pos_lim
-            ax.annotate(an, textcoords='data', xy=[xtext_pos, ytext_pos],
-                        fontsize=8)
-    ax.plot([0, 0], [ylim[0], ylim[1]],
-            color='gray', linestyle='dashed', linewidth=0.5)
-
-    # omega-constant lines
-    angles = np.linspace(-90, 90, 20) * np.pi/180
-    if wn is None:
-        wn = _default_wn(xlocator(), ylocator())
-
-    for om in wn:
-        if om < 0:
-            # Generate the lines for natural frequency curves
-            yp = np.sin(angles) * np.abs(om)
-            xp = -np.cos(angles) * np.abs(om)
-
-            # Plot the natural frequency contours
-            ax.plot(xp, yp, color='gray', linestyle='dashed', linewidth=0.5)
-
-            # Annotate the natural frequencies by listing on x-axis
-            # Note: need to filter values for proper plotting in Jupyter
-            if (om > xlim[0]):
-                an = "%.2f" % -om
-                ax.annotate(an, textcoords='data', xy=[om, 0], fontsize=8)
-
-
-def _default_zetas(xlim, ylim):
-    """Return default list of damping coefficients
-
-    This function computes a list of damping coefficients based on the limits
-    of the graph.  A set of 4 damping coefficients are computed for the x-axis
-    and a set of three damping coefficients are computed for the y-axis
-    (corresponding to the normal 4:3 plot aspect ratio in `matplotlib`?).
-
-    Parameters
-    ----------
-    xlim : array_like
-        List of x-axis limits [min, max]
-    ylim : array_like
-        List of y-axis limits [min, max]
-
-    Returns
-    -------
-    zeta : list
-        List of default damping coefficients for the plot
-
-    """
-    # Damping coefficient lines that intersect the x-axis
-    sep1 = -xlim[0] / 4
-    ang1 = [np.arctan((sep1*i)/ylim[1]) for i in np.arange(1, 4, 1)]
-
-    # Damping coefficient lines that intersection the y-axis
-    sep2 = ylim[1] / 3
-    ang2 = [np.arctan(-xlim[0]/(ylim[1]-sep2*i)) for i in np.arange(1, 3, 1)]
-
-    # Put the lines together and add one at -pi/2 (negative real axis)
-    angles = np.concatenate((ang1, ang2))
-    angles = np.insert(angles, len(angles), np.pi/2)
-
-    # Return the damping coefficients corresponding to these angles
-    zeta = np.sin(angles)
-    return zeta.tolist()
-
-
-def _default_wn(xloc, yloc, max_lines=7):
-    """Return default wn for root locus plot
-
-    This function computes a list of natural frequencies based on the grid
-    parameters of the graph.
-
-    Parameters
-    ----------
-    xloc : array_like
-        List of x-axis tick values
-    ylim : array_like
-        List of y-axis limits [min, max]
-    max_lines : int, optional
-        Maximum number of frequencies to generate (default = 7)
-
-    Returns
-    -------
-    wn : list
-        List of default natural frequencies for the plot
-
-    """
-    sep = xloc[1]-xloc[0]       # separation between x-ticks
-
-    # Decide whether to use the x or y axis for determining wn
-    if yloc[-1] / sep > max_lines*10:
-        # y-axis scale >> x-axis scale
-        wn = yloc                   # one frequency per y-axis tick mark
-    else:
-        wn = xloc                   # one frequency per x-axis tick mark
-
-        # Insert additional frequencies to span the y-axis
-        while np.abs(wn[0]) < yloc[-1]:
-            wn = np.insert(wn, 0, wn[0]-sep)
-
-    # If there are too many values, cut them in half
-    while len(wn) > max_lines:
-        wn = wn[0:-1:2]
-
-    return wn
-
-
-rlocus = root_locus
+# Alternative ways to call these functions
+root_locus = root_locus_plot
+rlocus = root_locus_plot
diff --git a/control/sisotool.py b/control/sisotool.py
index 9af5268b..82a989ca 100644
--- a/control/sisotool.py
+++ b/control/sisotool.py
@@ -1,19 +1,23 @@
 __all__ = ['sisotool', 'rootlocus_pid_designer']
 
+import warnings
+from functools import partial
+
+import matplotlib.pyplot as plt
+import numpy as np
+
 from control.exception import ControlMIMONotImplemented
+from control.statesp import _convert_to_statespace
+
+from . import config
+from .bdalg import append, connect
 from .freqplot import bode_plot
-from .timeresp import step_response
 from .iosys import common_timebase, isctime, isdtime
 from .lti import frequency_response
-from .xferfcn import tf
-from .statesp import ss, summing_junction
-from .bdalg import append, connect
 from .nlsys import interconnect
-from control.statesp import _convert_to_statespace
-from . import config
-import numpy as np
-import matplotlib.pyplot as plt
-import warnings
+from .statesp import ss, summing_junction
+from .timeresp import step_response
+from .xferfcn import tf
 
 _sisotool_defaults = {
     'sisotool.initial_gain': 1
@@ -86,7 +90,7 @@ def sisotool(sys, initial_gain=None, xlim_rlocus=None, ylim_rlocus=None,
     >>> ct.sisotool(G)                                          # doctest: +SKIP
 
     """
-    from .rlocus import root_locus
+    from .rlocus import root_locus_map
 
     # sys as loop transfer function if SISO
     if not sys.issiso():
@@ -123,13 +127,51 @@ def sisotool(sys, initial_gain=None, xlim_rlocus=None, ylim_rlocus=None,
     initial_gain = config._get_param('sisotool', 'initial_gain',
             initial_gain, _sisotool_defaults)
 
-    # First time call to setup the bode and step response plots
+    # First time call to setup the Bode and step response plots
     _SisotoolUpdate(sys, fig, initial_gain, bode_plot_params)
 
-    # Setup the root-locus plot window
-    root_locus(sys, initial_gain=initial_gain, xlim=xlim_rlocus,
-        ylim=ylim_rlocus, plotstr=plotstr_rlocus, grid=rlocus_grid,
-        fig=fig, bode_plot_params=bode_plot_params, tvect=tvect, sisotool=True)
+    # root_locus(
+    #     sys[0, 0], initial_gain=initial_gain, xlim=xlim_rlocus,
+    #     ylim=ylim_rlocus, plotstr=plotstr_rlocus, grid=rlocus_grid,
+    #     ax=fig.axes[1])
+    ax_rlocus = fig.axes[1]
+    root_locus_map(sys[0, 0]).plot(
+        xlim=xlim_rlocus, ylim=ylim_rlocus, grid=rlocus_grid,
+        initial_gain=initial_gain, ax=ax_rlocus)
+    if rlocus_grid is False:
+        # Need to generate grid manually, since root_locus_plot() won't
+        from .grid import nogrid
+        nogrid(sys.dt, ax=ax_rlocus)
+
+    # Reset the button release callback so that we can update all plots
+    fig.canvas.mpl_connect(
+        'button_release_event', partial(
+            _click_dispatcher, sys=sys, ax=fig.axes[1],
+            bode_plot_params=bode_plot_params, tvect=tvect))
+
+
+def _click_dispatcher(event, sys, ax, bode_plot_params, tvect):
+    # Zoom handled by specialized callback in rlocus, only handle gain plot
+    if event.inaxes == ax.axes and \
+       plt.get_current_fig_manager().toolbar.mode not in \
+       {'zoom rect', 'pan/zoom'}:
+        fig = ax.figure
+
+        # if a point is clicked on the rootlocus plot visually emphasize it
+        # K = _RLFeedbackClicksPoint(
+        #     event, sys, fig, ax_rlocus, show_clicked=True)
+        from .pzmap import _create_root_locus_label, _find_root_locus_gain, \
+            _mark_root_locus_gain
+
+        K, s = _find_root_locus_gain(event, sys, ax)
+        if K is not None:
+            _mark_root_locus_gain(ax, sys, K)
+            fig.suptitle(_create_root_locus_label(sys, K, s), fontsize=10)
+            _SisotoolUpdate(sys, fig, K, bode_plot_params, tvect)
+
+        # Update the canvas
+        fig.canvas.draw()
+
 
 def _SisotoolUpdate(sys, fig, K, bode_plot_params, tvect=None):
 
diff --git a/doc/Makefile b/doc/Makefile
index a5f7ec5a..71e493f2 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -16,7 +16,7 @@ help:
 
 # Rules to create figures
 FIGS = classes.pdf timeplot-mimo_step-default.png \
-  freqplot-siso_bode-default.png
+  freqplot-siso_bode-default.png rlocus-siso_ctime-default.png
 classes.pdf: classes.fig
 	fig2dev -Lpdf $< $@
 
@@ -26,6 +26,9 @@ timeplot-mimo_step-default.png: ../control/tests/timeplot_test.py
 freqplot-siso_bode-default.png: ../control/tests/freqplot_test.py
 	PYTHONPATH=.. python $<
 
+rlocus-siso_ctime-default.png: ../control/tests/rlocus_test.py
+	PYTHONPATH=.. python $<
+
 # Catch-all target: route all unknown targets to Sphinx using the new
 # "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
 html pdf clean doctest: Makefile $(FIGS)
diff --git a/doc/plotting.rst b/doc/plotting.rst
index be7ae7a5..2f6857c3 100644
--- a/doc/plotting.rst
+++ b/doc/plotting.rst
@@ -4,15 +4,15 @@
 Plotting data
 *************
 
-The Python Control Toolbox contains a number of functions for plotting
-input/output responses in the time and frequency domain, root locus
-diagrams, and other standard charts used in control system analysis, for
-example::
+The Python Control Systems Toolbox contains a number of functions for
+plotting input/output responses in the time and frequency domain, root
+locus diagrams, and other standard charts used in control system analysis,
+for example::
 
   bode_plot(sys)
   nyquist_plot([sys1, sys2])
-
-.. root_locus_plot(sys)		# not yet implemented
+  pole_zero_plot(sys)
+  root_locus_plot(sys)
 
 While plotting functions can be called directly, the standard pattern used
 in the toolbox is to provide a function that performs the basic computation
@@ -35,7 +35,7 @@ analysis object, allowing the following type of calls::
   step_response(sys).plot()
   frequency_response(sys).plot()
   nyquist_response(sys).plot()
-  rootlocus_response(sys).plot()     # implementation pending
+  root_locus_map(sys).plot()
 
 The remainder of this chapter provides additional documentation on how
 these response and plotting functions can be customized.
@@ -222,6 +222,58 @@ sensitivity functions for a feedback control system in standard form::
 .. image:: freqplot-gangof4.png
 
 
+Pole/zero data
+==============
+
+Pole/zero maps and root locus diagrams provide insights into system
+response based on the locations of system poles and zeros in the complex
+plane.  The :func:`~control.pole_zero_map` function returns the poles and
+zeros and can be used to generate a pole/zero plot::
+
+  sys = ct.tf([1, 2], [1, 2, 3], name='SISO transfer function')
+  response = ct.pole_zero_map(sys)
+  ct.pole_zero_plot(response)
+
+.. image:: pzmap-siso_ctime-default.png
+
+A root locus plot shows the location of the closed loop poles of a system
+as a function of the loop gain::
+
+  ct.root_locus_map(sys).plot()
+
+.. image:: rlocus-siso_ctime-default.png
+
+The grid in the left hand plane shows lines of constant damping ratio as
+well as arcs corresponding to the frequency of the complex pole.  The grid
+can be turned off using the `grid` keyword.  Setting `grid` to `False` will
+turn off the grid but show the real and imaginary axis.  To completely
+remove all lines except the root loci, use `grid='empty'`.
+
+On systems that support interactive plots, clicking on a location on the
+root locus diagram will mark the pole locations on all branches of the
+diagram and display the gain and damping ratio for the clicked point below
+the plot title:
+
+.. image:: rlocus-siso_ctime-clicked.png
+
+Root locus diagrams are also supported for discrete time systems, in which
+case the grid is show inside the unit circle::
+
+  sysd = sys.sample(0.1)
+  ct.root_locus_plot(sysd)
+
+.. image:: rlocus-siso_dtime-default.png
+
+Lists of systems can also be given, in which case the root locus diagram
+for each system is plotted in different colors::
+
+  sys1 = ct.tf([1], [1, 2, 1], name='sys1')
+  sys2 = ct.tf([1, 0.2], [1, 1, 3, 1, 1], name='sys2')
+  ct.root_locus_plot([sys1, sys2], grid=False)
+
+.. image:: rlocus-siso_multiple-nogrid.png
+
+
 Response and plotting functions
 ===============================
 
@@ -244,6 +296,8 @@ number of encirclements for a Nyquist plot) as well as plotting (via the
    ~control.initial_response
    ~control.input_output_response
    ~control.nyquist_response
+   ~control.pole_zero_map
+   ~control.root_locus_map
    ~control.singular_values_response
    ~control.step_response
 
@@ -256,6 +310,8 @@ Plotting functions
    ~control.bode_plot
    ~control.describing_function_plot
    ~control.nichols_plot
+   ~control.pole_zero_plot
+   ~control.root_locus_plot
    ~control.singular_values_plot
    ~control.time_response_plot
 
@@ -284,4 +340,5 @@ The following classes are used in generating response data.
    ~control.DescribingFunctionResponse
    ~control.FrequencyResponseData
    ~control.NyquistResponseData
+   ~control.PoleZeroData
    ~control.TimeResponseData


Test patch:
diff --git a/control/tests/kwargs_test.py b/control/tests/kwargs_test.py
index 0d13e639..53cc0076 100644
--- a/control/tests/kwargs_test.py
+++ b/control/tests/kwargs_test.py
@@ -97,6 +97,7 @@ def test_kwarg_search(module, prefix):
      (control.pzmap, 1, 0, (), {}),
      (control.rlocus, 0, 1, (), {}),
      (control.root_locus, 0, 1, (), {}),
+     (control.root_locus_plot, 0, 1, (), {}),
      (control.rss, 0, 0, (2, 1, 1), {}),
      (control.set_defaults, 0, 0, ('control',), {'default_dt': True}),
      (control.ss, 0, 0, (0, 0, 0, 0), {'dt': 1}),
@@ -175,6 +176,8 @@ def test_matplotlib_kwargs(function, nsysargs, moreargs, kwargs, mplcleanup):
         (control.frequency_response, control.bode, True),
         (control.frequency_response, control.bode_plot, True),
         (control.nyquist_response, control.nyquist_plot, False),
+        (control.pole_zero_map, control.pole_zero_plot, False),
+        (control.root_locus_map, control.root_locus_plot, False),
     ])
 def test_response_plot_kwargs(data_fcn, plot_fcn, mimo):
     # Create a system for testing
@@ -193,16 +196,18 @@ def test_response_plot_kwargs(data_fcn, plot_fcn, mimo):
     plot_fcn(response)
 
     # Now add an unrecognized keyword and make sure there is an error
-    with pytest.raises(AttributeError,
-                       match="(has no property|unexpected keyword)"):
+    with pytest.raises(
+            (AttributeError, TypeError),
+            match="(has no property|unexpected keyword|unrecognized keyword)"):
         plot_fcn(response, unknown=None)
 
     # Call the plotting function via the response and make sure it works
     response.plot()
 
     # Now add an unrecognized keyword and make sure there is an error
-    with pytest.raises(AttributeError,
-                       match="(has no property|unexpected keyword)"):
+    with pytest.raises(
+            (AttributeError, TypeError),
+            match="(has no property|unexpected keyword|unrecognized keyword)"):
         response.plot(unknown=None)
 
 #
@@ -241,8 +246,10 @@ kwarg_unittest = {
     'nyquist_response': test_response_plot_kwargs,
     'nyquist_plot': test_matplotlib_kwargs,
     'pzmap': test_unrecognized_kwargs,
+    'pole_zero_plot': test_unrecognized_kwargs,
     'rlocus': test_unrecognized_kwargs,
     'root_locus': test_unrecognized_kwargs,
+    'root_locus_plot': test_unrecognized_kwargs,
     'rss': test_unrecognized_kwargs,
     'set_defaults': test_unrecognized_kwargs,
     'singular_values_plot': test_matplotlib_kwargs,
@@ -274,6 +281,7 @@ kwarg_unittest = {
     'flatsys.LinearFlatSystem.__init__': test_unrecognized_kwargs,
     'NonlinearIOSystem.linearize': test_unrecognized_kwargs,
     'NyquistResponseData.plot': test_response_plot_kwargs,
+    'PoleZeroData.plot': test_response_plot_kwargs,
     'InterconnectedSystem.__init__':
         interconnect_test.test_interconnect_exceptions,
     'StateSpace.__init__':
diff --git a/control/tests/matlab_test.py b/control/tests/matlab_test.py
index e01abcca..2ba3d5df 100644
--- a/control/tests/matlab_test.py
+++ b/control/tests/matlab_test.py
@@ -424,7 +424,8 @@ class TestMatlab:
     @pytest.mark.parametrize("subsys", ["ss1", "tf1", "tf2"])
     def testRlocus(self, siso, subsys, mplcleanup):
         """Call rlocus()"""
-        rlocus(getattr(siso, subsys))
+        rlist, klist = rlocus(getattr(siso, subsys))
+        np.testing.assert_equal(len(rlist), len(klist))
 
     def testRlocus_list(self, siso, mplcleanup):
         """Test rlocus() with list"""
diff --git a/control/tests/pzmap_test.py b/control/tests/pzmap_test.py
index 8d41807b..ce8adf6e 100644
--- a/control/tests/pzmap_test.py
+++ b/control/tests/pzmap_test.py
@@ -12,9 +12,11 @@ import pytest
 from matplotlib import pyplot as plt
 from mpl_toolkits.axisartist import Axes as mpltAxes
 
+import control as ct
 from control import TransferFunction, config, pzmap
 
 
+@pytest.mark.filterwarnings("ignore:.*return values.*:DeprecationWarning")
 @pytest.mark.parametrize("kwargs",
                          [pytest.param(dict(), id="default"),
                           pytest.param(dict(plot=False), id="plot=False"),
@@ -44,20 +46,23 @@ def test_pzmap(kwargs, setdefaults, dt, editsdefaults, mplcleanup):
 
     pzkwargs = kwargs.copy()
     if setdefaults:
-        for k in ['plot', 'grid']:
+        for k in ['grid']:
             if k in pzkwargs:
                 v = pzkwargs.pop(k)
                 config.set_defaults('pzmap', **{k: v})
 
+    if kwargs.get('plot', None) is None:
+        pzkwargs['plot'] = True         # use to get legacy return values
     P, Z = pzmap(T, **pzkwargs)
 
     np.testing.assert_allclose(P, Pref, rtol=1e-3)
     np.testing.assert_allclose(Z, Zref, rtol=1e-3)
 
     if kwargs.get('plot', True):
-        ax = plt.gca()
+        fig, ax = plt.gcf(), plt.gca()
 
-        assert ax.get_title() == kwargs.get('title', 'Pole Zero Map')
+        assert fig._suptitle.get_text().startswith(
+            kwargs.get('title', 'Pole/zero plot'))
 
         # FIXME: This won't work when zgrid and sgrid are unified
         children = ax.get_children()
@@ -78,12 +83,43 @@ def test_pzmap(kwargs, setdefaults, dt, editsdefaults, mplcleanup):
         assert not plt.get_fignums()
 
 
-def test_pzmap_warns():
-    with pytest.warns(FutureWarning):
-        pzmap(TransferFunction([1], [1, 2]), Plot=True)
+def test_polezerodata():
+    sys = ct.rss(4, 1, 1)
+    pzdata = ct.pole_zero_map(sys)
+    np.testing.assert_equal(pzdata.poles, sys.poles())
+    np.testing.assert_equal(pzdata.zeros, sys.zeros())
+
+    # Extract data from PoleZeroData
+    poles, zeros = pzdata
+    np.testing.assert_equal(poles, sys.poles())
+    np.testing.assert_equal(zeros, sys.zeros())
+
+    # Legacy return format
+    for plot in [True, False]:
+        with pytest.warns(DeprecationWarning, match=".* values .* deprecated"):
+            poles, zeros = ct.pole_zero_plot(pzdata, plot=False)
+        np.testing.assert_equal(poles, sys.poles())
+        np.testing.assert_equal(zeros, sys.zeros())
 
 
 def test_pzmap_raises():
     with pytest.raises(TypeError):
         # not an LTI system
-        pzmap(([1], [1,2]))
+        pzmap(([1], [1, 2]))
+
+    sys1 = ct.rss(2, 1, 1)
+    sys2 = sys1.sample(0.1)
+    with pytest.raises(ValueError, match="incompatible time bases"):
+        pzdata = ct.pole_zero_plot([sys1, sys2], grid=True)
+
+    with pytest.warns(UserWarning, match="axis already exists"):
+        fig, ax = plt.figure(), plt.axes()
+        ct.pole_zero_plot(sys1, ax=ax, grid='empty')
+
+
+def test_pzmap_limits():
+    sys = ct.tf([1, 2], [1, 2, 3])
+    out = ct.pole_zero_plot(sys, xlim=[-1, 1], ylim=[-1, 1])
+    ax = ct.get_plot_axes(out)[0, 0]
+    assert ax.get_xlim() == (-1, 1)
+    assert ax.get_ylim() == (-1, 1)
diff --git a/control/tests/rlocus_test.py b/control/tests/rlocus_test.py
index e61f0c8f..5ee2830c 100644
--- a/control/tests/rlocus_test.py
+++ b/control/tests/rlocus_test.py
@@ -9,7 +9,7 @@ from numpy.testing import assert_array_almost_equal
 import pytest
 
 import control as ct
-from control.rlocus import root_locus, _RLClickDispatcher
+from control.rlocus import root_locus
 from control.xferfcn import TransferFunction
 from control.statesp import StateSpace
 from control.bdalg import feedback
@@ -45,6 +45,7 @@ class TestRootLocus:
             poles = np.sort(poles)
             np.testing.assert_array_almost_equal(poles, poles_expected)
 
+    @pytest.mark.filterwarnings("ignore:.*return values.*:DeprecationWarning")
     def testRootLocus(self, sys):
         """Basic root locus (no plot)"""
         klist = [-1, 0, 1]
@@ -55,50 +56,75 @@ class TestRootLocus:
         self.check_cl_poles(sys, roots, klist)
 
         # now check with plotting
-        roots, k_out = root_locus(sys, klist)
+        roots, k_out = root_locus(sys, klist, plot=True)
         np.testing.assert_equal(len(roots), len(klist))
         np.testing.assert_allclose(klist, k_out)
         self.check_cl_poles(sys, roots, klist)
 
+    @pytest.mark.filterwarnings("ignore:.*return values.*:DeprecationWarning")
     def test_without_gains(self, sys):
         roots, kvect = root_locus(sys, plot=False)
         self.check_cl_poles(sys, roots, kvect)
 
-    @pytest.mark.slow
-    @pytest.mark.parametrize('grid', [None, True, False])
-    def test_root_locus_plot_grid(self, sys, grid):
-        rlist, klist = root_locus(sys, grid=grid)
+    @pytest.mark.parametrize("grid", [None, True, False, 'empty'])
+    @pytest.mark.parametrize("method", ['plot', 'map', 'response', 'pzmap'])
+    def test_root_locus_plot_grid(self, sys, grid, method):
+        import mpl_toolkits.axisartist as AA
+
+        # Generate the root locus plot
+        plt.clf()
+        if method == 'plot':
+            ct.root_locus_plot(sys, grid=grid)
+        elif method == 'map':
+            ct.root_locus_map(sys).plot(grid=grid)
+        elif method == 'response':
+            response = ct.root_locus_map(sys)
+            ct.root_locus_plot(response, grid=grid)
+        elif method == 'pzmap':
+            response = ct.root_locus_map(sys)
+            ct.pole_zero_plot(response, grid=grid)
+
+        # Count the number of dotted/dashed lines in the plot
         ax = plt.gca()
-        n_gridlines = sum([int(line.get_linestyle() in [':', 'dotted',
-                                                        '--', 'dashed'])
-                           for line in ax.lines])
-        if grid is False:
-            assert n_gridlines == 2
-        else:
+        n_gridlines = sum([int(
+            line.get_linestyle() in [':', 'dotted', '--', 'dashed'] or
+            line.get_linewidth() < 1
+        ) for line in ax.lines])
+
+        # Make sure they line up with what we expect
+        if grid == 'empty':
+            assert n_gridlines == 0
+            assert not isinstance(ax, AA.Axes)
+        elif grid is False or method == 'pzmap' and grid is None:
+            assert n_gridlines == 2 if sys.isctime() else 3
+            assert not isinstance(ax, AA.Axes)
+        elif sys.isdtime(strict=True):
             assert n_gridlines > 2
-        # TODO check validity of grid
+            assert not isinstance(ax, AA.Axes)
+        else:
+            # Continuous time, with grid => check that AxisArtist was used
+            assert isinstance(ax, AA.Axes)
+            for spine in ['wnxneg', 'wnxpos', 'wnyneg', 'wnypos']:
+                assert spine in ax.axis
 
-    def test_root_locus_warnings(self):
-        sys = TransferFunction([1000], [1, 25, 100, 0])
-        with pytest.warns(FutureWarning, match="Plot.*deprecated"):
-            rlist, klist = root_locus(sys, Plot=True)
-        with pytest.warns(FutureWarning, match="PrintGain.*deprecated"):
-            rlist, klist = root_locus(sys, PrintGain=True)
+        # TODO: check validity of grid
 
+    @pytest.mark.filterwarnings("ignore:.*return values.*:DeprecationWarning")
     def test_root_locus_neg_false_gain_nonproper(self):
         """ Non proper TranferFunction with negative gain: Not implemented"""
         with pytest.raises(ValueError, match="with equal order"):
-            root_locus(TransferFunction([-1, 2], [1, 2]))
+            root_locus(TransferFunction([-1, 2], [1, 2]), plot=True)
 
     # TODO: cover and validate negative false_gain branch in _default_gains()
 
+    @pytest.mark.skip("Zooming functionality no longer implemented")
     @pytest.mark.skipif(plt.get_current_fig_manager().toolbar is None,
                         reason="Requires the zoom toolbar")
     def test_root_locus_zoom(self):
         """Check the zooming functionality of the Root locus plot"""
         system = TransferFunction([1000], [1, 25, 100, 0])
         plt.figure()
-        root_locus(system)
+        root_locus(system, plot=True)
         fig = plt.gcf()
         ax_rlocus = fig.axes[0]
 
@@ -121,6 +147,7 @@ class TestRootLocus:
         assert_array_almost_equal(zoom_x, zoom_x_valid)
         assert_array_almost_equal(zoom_y, zoom_y_valid)
 
+    @pytest.mark.filterwarnings("ignore:.*return values.*:DeprecationWarning")
     @pytest.mark.timeout(2)
     def test_rlocus_default_wn(self):
         """Check that default wn calculation works properly"""
@@ -140,4 +167,104 @@ class TestRootLocus:
         sys = ct.tf(*sp.signal.zpk2tf(
             [-1e-2, 1-1e7j, 1+1e7j], [0, -1e7j, 1e7j], 1))
 
-        ct.root_locus(sys)
+        ct.root_locus(sys, plot=True)
+
+
+@pytest.mark.parametrize(
+    "sys, grid, xlim, ylim, interactive", [
+        (ct.tf([1], [1, 2, 1]), None, None, None, False),
+    ])
+def test_root_locus_plots(sys, grid, xlim, ylim, interactive):
+    ct.root_locus_map(sys).plot(
+        grid=grid, xlim=xlim, ylim=ylim, interactive=interactive)
+    # TODO: add tests to make sure everything "looks" OK
+
+
+# Generate plots used in documentation
+def test_root_locus_documentation(savefigs=False):
+    plt.figure()
+    sys = ct.tf([1, 2], [1, 2, 3], name='SISO transfer function')
+    response = ct.pole_zero_map(sys)
+    ct.pole_zero_plot(response)
+    plt.savefig('pzmap-siso_ctime-default.png')
+
+    plt.figure()
+    ct.root_locus_map(sys).plot()
+    plt.savefig('rlocus-siso_ctime-default.png')
+
+    # TODO: generate event in order to generate real title
+    plt.figure()
+    out = ct.root_locus_map(sys).plot(initial_gain=3.506)
+    ax = ct.get_plot_axes(out)[0, 0]
+    freqplot_rcParams = ct.config._get_param('freqplot', 'rcParams')
+    with plt.rc_context(freqplot_rcParams):
+        ax.set_title(
+            "Clicked at: -2.729+1.511j  gain = 3.506  damping = 0.8748")
+    plt.savefig('rlocus-siso_ctime-clicked.png')
+
+    plt.figure()
+    sysd = sys.sample(0.1)
+    ct.root_locus_plot(sysd)
+    plt.savefig('rlocus-siso_dtime-default.png')
+
+    plt.figure()
+    sys1 = ct.tf([1, 2], [1, 2, 3], name='sys1')
+    sys2 = ct.tf([1, 0.2], [1, 1, 3, 1, 1], name='sys2')
+    ct.root_locus_plot([sys1, sys2], grid=False)
+    plt.savefig('rlocus-siso_multiple-nogrid.png')
+
+
+if __name__ == "__main__":
+    #
+    # Interactive mode: generate plots for manual viewing
+    #
+    # Running this script in python (or better ipython) will show a
+    # collection of figures that should all look OK on the screeen.
+    #
+
+    # In interactive mode, turn on ipython interactive graphics
+    plt.ion()
+
+    # Start by clearing existing figures
+    plt.close('all')
+
+    # Define systems to be tested
+    sys_secord = ct.tf([1], [1, 1, 1], name="2P")
+    sys_seczero = ct.tf([1, 0, -1], [1, 1, 1], name="2P, 2Z")
+    sys_fbs_a = ct.tf([1, 1], [1, 0, 0], name="FBS 12_19a")
+    sys_fbs_b = ct.tf(
+        ct.tf([1, 1], [1, 2, 0]) * ct.tf([1], [1, 2 ,4]), name="FBS 12_19b")
+    sys_fbs_c = ct.tf([1, 1], [1, 0, 1, 0], name="FBS 12_19c")
+    sys_fbs_d = ct.tf([1, 2, 2], [1, 0, 1, 0], name="FBS 12_19d")
+    sys_poles = sys_fbs_d.poles()
+    sys_zeros = sys_fbs_d.zeros()
+    sys_discrete = ct.zpk(
+        sys_zeros / 3, sys_poles / 3, 1, dt=True, name="discrete")
+
+    # Run through a large number of test cases
+    test_cases = [
+        # sys          grid   xlim      ylim      inter
+        (sys_secord,   None,  None,     None,     None),
+        (sys_seczero,  None,  None,     None,     None),
+        (sys_fbs_a,    None,  None,     None,     None),
+        (sys_fbs_b,    None,  None,     None,     False),
+        (sys_fbs_c,    None,  None,     None,     None),
+        (sys_fbs_c,    None,  None,     [-2, 2],  None),
+        (sys_fbs_c,    True,  [-3, 3],  None,     None),
+        (sys_fbs_d,    None,  None,     None,     None),
+        (ct.zpk(sys_zeros * 10, sys_poles * 10, 1, name="12_19d * 10"),
+                       None,  None,     None,     None),
+        (ct.zpk(sys_zeros / 10, sys_poles / 10, 1, name="12_19d / 10"),
+                       True,  None,     None,     None),
+        (sys_discrete, None,  None,     None,     None),
+        (sys_discrete, True,  None,     None,     None),
+        (sys_fbs_d,    True,  None,     None,     True),
+    ]
+
+    for sys, grid, xlim, ylim, interactive in test_cases:
+        plt.figure()
+        test_root_locus_plots(
+            sys, grid=grid, xlim=xlim, ylim=ylim, interactive=interactive)
+
+    # Run tests that generate plots for the documentation
+    test_root_locus_documentation(savefigs=True)
diff --git a/control/tests/sisotool_test.py b/control/tests/sisotool_test.py
index 5a86c73d..b5aaeb90 100644
--- a/control/tests/sisotool_test.py
+++ b/control/tests/sisotool_test.py
@@ -7,7 +7,7 @@ from numpy.testing import assert_array_almost_equal
 import pytest
 
 from control.sisotool import sisotool, rootlocus_pid_designer
-from control.rlocus import _RLClickDispatcher
+from control.sisotool import _click_dispatcher
 from control.xferfcn import TransferFunction
 from control.statesp import StateSpace
 from control import c2d
@@ -57,11 +57,11 @@ class TestSisotool:
         initial_point_0 = (np.array([-22.53155977]), np.array([0.]))
         initial_point_1 = (np.array([-1.23422011]), np.array([-6.54667031]))
         initial_point_2 = (np.array([-1.23422011]), np.array([6.54667031]))
-        assert_array_almost_equal(ax_rlocus.lines[0].get_data(),
+        assert_array_almost_equal(ax_rlocus.lines[4].get_data(),
                                   initial_point_0, 4)
-        assert_array_almost_equal(ax_rlocus.lines[1].get_data(),
+        assert_array_almost_equal(ax_rlocus.lines[5].get_data(),
                                   initial_point_1, 4)
-        assert_array_almost_equal(ax_rlocus.lines[2].get_data(),
+        assert_array_almost_equal(ax_rlocus.lines[6].get_data(),
                                   initial_point_2, 4)
 
         # Check the step response before moving the point
@@ -93,9 +93,8 @@ class TestSisotool:
         event = type('test', (object,), {'xdata': 2.31206868287,
                                          'ydata': 15.5983051046,
                                          'inaxes': ax_rlocus.axes})()
-        _RLClickDispatcher(event=event, sys=tsys, fig=fig,
-                           ax_rlocus=ax_rlocus, sisotool=True, plotstr='-',
-                           bode_plot_params=bode_plot_params, tvect=None)
+        _click_dispatcher(event=event, sys=tsys, ax=ax_rlocus,
+                          bode_plot_params=bode_plot_params, tvect=None)
 
         # Check the moved root locus plot points
         moved_point_0 = (np.array([-29.91742755]), np.array([0.]))
@@ -143,9 +142,8 @@ class TestSisotool:
         event = type('test', (object,), {'xdata': 2.31206868287,
                                          'ydata': 15.5983051046,
                                          'inaxes': ax_rlocus.axes})()
-        _RLClickDispatcher(event=event, sys=tsys, fig=fig,
-                           ax_rlocus=ax_rlocus, sisotool=True, plotstr='-',
-                           bode_plot_params=dict(), tvect=tvect)
+        _click_dispatcher(event=event, sys=tsys, ax=ax_rlocus,
+                          bode_plot_params=dict(), tvect=tvect)
         assert_array_almost_equal(tvect, ax_step.lines[0].get_data()[0])
 
     @pytest.mark.skipif(plt.get_current_fig_manager().toolbar is None,
@@ -202,3 +200,21 @@ class TestPidDesigner:
     def test_pid_designer_2(self, plant, kwargs):
         rootlocus_pid_designer(plant, **kwargs)
 
+
+if __name__ == "__main__":
+    #
+    # Interactive mode: generate plots for manual viewing
+    #
+    # Running this script in python (or better ipython) will show a
+    # collection of figures that should all look OK on the screeen.
+    #
+    import control as ct
+
+    # In interactive mode, turn on ipython interactive graphics
+    plt.ion()
+
+    # Start by clearing existing figures
+    plt.close('all')
+
+    tsys = ct.tf([1000], [1, 25, 100, 0])
+    ct.sisotool(tsys)