Repository: python-control/python-control
Pull request #955 (2023-12-29 23:56:09)

Fix CI warnings 
This PR makes a bunch of small changes to remove warning messages in the GitHub actions CI runs:
* Access scalar NumPy arrays using `item()` (old usage is newly deprecated in NumPy).
* Change `set_tight_layout(True)` to `fig.set_layout_engine('tight')` (newly deprecated in Matplotlib).
* Update unit tests using deprecated functionality to use current functionality.
* Add warning filters for internal calls to `connect` (which are OK for now).
* Avoid calling tf2zpk with zero in the numerator (SciPy generates warning).
* Fix up some strings that had newly deprecated escape sequences (by using raw strings)
* Turn off warnings for optimal control problems that fail to hit desired precision (but still return valid answers)

Note: it is likely that this PR will conflict with #953.  We should probably merge #953 first and then rebase this PR since the changes here are pretty straightforward.


Patch:
diff --git a/control/flatsys/linflat.py b/control/flatsys/linflat.py
index 9320eec3..e03df514 100644
--- a/control/flatsys/linflat.py
+++ b/control/flatsys/linflat.py
@@ -119,10 +119,10 @@ class LinearFlatSystem(FlatSystem, StateSpace):
         x = np.reshape(x, (-1, 1))
         u = np.reshape(u, (1, -1))
         zflag = [np.zeros(self.nstates + 1)]
-        zflag[0][0] = self.Cf @ x
+        zflag[0][0] = (self.Cf @ x).item()
         H = self.Cf                     # initial state transformation
         for i in range(1, self.nstates + 1):
-            zflag[0][i] = H @ (self.A @ x + self.B @ u)
+            zflag[0][i] = (H @ (self.A @ x + self.B @ u)).item()
             H = H @ self.A       # derivative for next iteration
         return zflag
 
diff --git a/control/freqplot.py b/control/freqplot.py
index 57f24f8d..6c07b178 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 by config.default['freqplot.rcParams'].
+        Default is set up 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
@@ -473,7 +473,7 @@ def bode_plot(
     if ax is None:
         with plt.rc_context(_freqplot_rcParams):
             ax_array = fig.subplots(nrows, ncols, squeeze=False)
-            fig.set_tight_layout(True)
+            fig.set_layout_engine('tight')
             fig.align_labels()
 
         # Set up default sharing of axis limits if not specified
diff --git a/control/grid.py b/control/grid.py
index ef999594..785ec274 100644
--- a/control/grid.py
+++ b/control/grid.py
@@ -1,21 +1,12 @@
-# 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 numpy as np
+from numpy import cos, sin, sqrt, linspace, pi, exp
 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):
@@ -74,15 +65,14 @@ class ModifiedExtremeFinderCycle(angle_helper.ExtremeFinderCycle):
         return lon_min, lon_max, lat_min, lat_max
 
 
-def sgrid(scaling=None):
+def sgrid():
     # 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 degrees
+    # system in degree
     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).
@@ -99,7 +89,6 @@ def sgrid(scaling=None):
         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)
 
@@ -108,20 +97,15 @@ def sgrid(scaling=None):
     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("right")
-
+    axis.set_axis_direction("left")
     ax.axis["wnyneg"] = axis = ax.new_floating_axis(0, 270)
     axis.label.set_visible(False)
     axis.set_axis_direction("left")
@@ -135,41 +119,43 @@ def sgrid(scaling=None):
     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, scaling=scaling)
+    _final_setup(ax)
     return ax, fig
 
 
-# Utility function used by all grid code
-def _final_setup(ax, scaling=None):
+def _final_setup(ax):
     ax.set_xlabel('Real')
     ax.set_ylabel('Imaginary')
-    ax.axhline(y=0, color='black', lw=0.25)
-    ax.axvline(x=0, color='black', lw=0.25)
+    ax.axhline(y=0, color='black', lw=1)
+    ax.axvline(x=0, color='black', lw=1)
+    plt.axis('equal')
 
-    # 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()
 
-# 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)
+    return ax, f
 
-    _final_setup(ax, scaling=scaling)
-    return ax, fig
 
-# 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):
+def zgrid(zetas=None, wns=None, ax=None):
     """Draws discrete damping and frequency grid"""
 
     fig = plt.gcf()
@@ -220,9 +206,5 @@ def zgrid(zetas=None, wns=None, ax=None, scaling=None):
         ax.annotate(r"$\frac{"+num+r"\pi}{T}$", xy=(an_x, an_y),
                     xytext=(an_x, an_y), size=9)
 
-    # 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)
+    _final_setup(ax)
     return ax, fig
diff --git a/control/iosys.py b/control/iosys.py
index fbd5c1db..52262250 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -503,64 +503,42 @@ 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=None, strict=False, dt=None):
+def isdtime(sys, strict=False):
     """
     Check to see if a system is a discrete time system.
 
     Parameters
     ----------
-    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.
+    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
     """
 
-    # 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
+    # Check to see if this is a constant
     if isinstance(sys, (int, float, complex, np.number)):
-        # Constants OK as long as strict checking is off
+        # 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=None, dt=None, strict=False):
+def isctime(sys, strict=False):
     """
     Check to see if a system is a continuous-time system.
 
     Parameters
     ----------
-    sys : I/O system, optional
-        System to be checked.
-    dt : None or number, optional
-        Timebase to be checked.
+    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.
+        If strict is True, make sure that timebase is not None
     """
 
-    # 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
+    # Check to see if this is a constant
     if isinstance(sys, (int, float, complex, np.number)):
-        # Constants OK as long as strict checking is off
+        # 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 0384215a..b63b19c7 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', 'rlocus', 'pzmap', 'dcgain', 'connect']
+__all__ = ['bode', 'nyquist', 'ngrid', '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,106 +195,6 @@ 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()
@@ -354,7 +254,6 @@ 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/optimal.py b/control/optimal.py
index 81372705..ce80eccf 100644
--- a/control/optimal.py
+++ b/control/optimal.py
@@ -956,7 +956,7 @@ def solve_ocp(
         transpose=None, return_states=True, print_summary=True, log=False,
         **kwargs):
 
-    """Compute the solution to an optimal control problem.
+    r"""Compute the solution to an optimal control problem.
 
     The optimal trajectory (states and inputs) is computed so as to
     approximately mimimize a cost function of the following form (for
diff --git a/control/pzmap.py b/control/pzmap.py
index 59f378c3..5ee3d37c 100644
--- a/control/pzmap.py
+++ b/control/pzmap.py
@@ -1,609 +1,133 @@
 # pzmap.py - computations involving poles and zeros
 #
-# Original author: Richard M. Murray
+# Author: Richard M. Murray
 # Date: 7 Sep 2009
 #
 # This file contains functions that compute poles, zeros and related
-# 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.)
+# quantities for a linear system.
 #
-
-import itertools
-import warnings
-from math import pi
-
-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__ = ['pole_zero_map', 'pole_zero_plot', 'pzmap', 'PoleZeroData']
-
-
-# Define default parameter values for this module
-_pzmap_defaults = {
-    '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
-}
-
+# Copyright (c) 2009 by California Institute of Technology
+# All rights reserved.
 #
-# Classes for keeping track of pzmap plots
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
 #
-# The PoleZeroData class keeps track of the information that is on a
-# pole/zero plot.
+# 1. Redistributions of source code must retain the above copyright
+#    notice, this list of conditions and the following disclaimer.
 #
-# 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.
+# 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.
 #
-# 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.
+# 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.
 #
-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.
+from numpy import real, imag, linspace, exp, cos, sin, sqrt
+from math import pi
+from .lti import LTI
+from .iosys import isdtime, isctime
+from .grid import sgrid, zgrid, nogrid
+from . import config
 
-    """
-    # Convert the first argument to a list
-    syslist = sysdata if isinstance(sysdata, (list, tuple)) else [sysdata]
+__all__ = ['pzmap']
 
-    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]
+# Define default parameter values for this module
+_pzmap_defaults = {
+    'pzmap.grid': False,       # Plot omega-damping grid
+    'pzmap.plot': True,        # Generate plot using Matplotlib
+}
 
 
 # TODO: Implement more elegant cross-style axes. See:
-#    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):
+#    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):
     """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
     ----------
-    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,
+    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.
-        If this argument is present, the legacy value of poles and
-        zeros is returned.
+    grid: boolean (default = False)
+        If True plot omega-damping grid.
 
     Returns
     -------
-    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.
+    poles: array
+        The systems poles
+    zeros: array
+        The system's zeros.
 
     Notes
     -----
-    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.
+    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.
 
     """
-    # 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
-
-    # If argument was a singleton, turn it into a tuple
-    if not isinstance(data, (list, tuple)):
-        data = [data]
-
-    # 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
+    # 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')
 
-    # 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)
+    # Make sure there were no extraneous keywords
+    if kwargs:
+        raise TypeError("unrecognized keywords: ", str(kwargs))
 
-        # 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]
+    # Get parameter values
+    plot = config._get_param('pzmap', 'plot', plot, True)
+    grid = config._get_param('pzmap', 'grid', grid, False)
 
-    if plot is False:
-        if len(data) == 1:
-            return poles[0], zeros[0]
-        else:
-            return poles, zeros
+    if not isinstance(sys, LTI):
+        raise TypeError('Argument ``sys``: must be a linear system.')
 
-    # 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]
+    poles = sys.poles()
+    zeros = sys.zeros()
 
-    if len(axs) > 1:
-        # Need to generate a new figure
-        fig, axs = plt.figure(), []
+    if (plot):
+        import matplotlib.pyplot as plt
 
-    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)
+        if grid:
+            if isdtime(sys, strict=True):
+                ax, fig = zgrid()
             else:
-                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)]
+                ax, fig = sgrid()
         else:
-            color = maker_color
+            ax, fig = nogrid()
 
         # Plot the locations of the poles and zeros
         if len(poles) > 0:
-            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)
+            ax.scatter(real(poles), imag(poles), s=50, marker='x',
+                       facecolors='k')
         if len(zeros) > 0:
-            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]
+            ax.scatter(real(zeros), imag(zeros), s=50, marker='o',
+                       facecolors='none', edgecolors='k')
 
+        plt.title(title)
 
-pzmap = pole_zero_plot
+    # Return locations of poles and zeros as a tuple
+    return poles, zeros
diff --git a/control/rlocus.py b/control/rlocus.py
index ea17ae94..41cdec05 100644
--- a/control/rlocus.py
+++ b/control/rlocus.py
@@ -1,6 +1,38 @@
 # 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.
@@ -14,101 +46,48 @@
 # Sawyer B. Fuller (minster@uw.edu) 21 May 2020:
 #   * added compatibility with discrete-time systems.
 #
+# $Id$
 
-import warnings
+# Packages used by this module
 from functools import partial
-
-import matplotlib.pyplot as plt
 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
+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
 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_map', 'root_locus_plot', 'root_locus', 'rlocus']
+__all__ = ['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
 }
 
 
-# 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):
+# 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 plot.
 
-    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.
+    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.
 
     Parameters
     ----------
-    sysdata : PoleZeroMap or LTI object or list
+    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.
@@ -118,83 +97,182 @@ def root_locus_plot(
     ylim : tuple or list, optional
         Set limits of y axis, normally with tuple
         (see :doc:`matplotlib:api/axes_api`).
-    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'].
+    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.
     ax : :class:`matplotlib.axes.Axes`
         Axes on which to create root locus plot
     initial_gain : float, optional
-        Mark the point on the root locus diagram corresponding to the
-        given gain.
+        Used by :func:`sisotool` to indicate initial gain.
 
     Returns
     -------
-    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).
+    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
     -----
-    The root_locus_plot function calls matplotlib.pyplot.axis('equal'), which
+    The root_locus 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.
 
     """
-    from .pzmap import pole_zero_plot
-
-    # Set default parameters
+    # 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)
 
-    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:
-        responses = sysdata
+    # Check for sisotool mode
+    sisotool = kwargs.get('sisotool', False)
 
-    #
-    # 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 siso. sisotool has different requirements
+    if not sys.issiso() and not sisotool:
+        raise ControlMIMONotImplemented(
+            'sys must be single-input single-output (SISO)')
 
-    if plot is False:
-        return responses.loci, responses.gains
+    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
+    else:
+        kvect = np.atleast_1d(kvect)
+        root_array = _RLFindRoots(nump, denp, kvect)
+        root_array = _RLSortRoots(root_array)
+        recompute_on_zoom = False
 
-    # Plot the root loci
-    out = responses.plot(grid=grid, **kwargs)
+    if sisotool:
+        start_roots = _RLFindRoots(nump, denp, initial_gain)
 
-    # Legacy processing: return locations of poles and zeros as a tuple
-    if plot is True:
-        return responses.loci, responses.gains
+    # Make sure there were no extraneous keywords
+    if not sisotool and kwargs:
+        raise TypeError("unrecognized keywords: ", str(kwargs))
 
-    return out
+    # 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))
+
+    return root_array, kvect
 
 
-def _default_gains(num, den, xlim, ylim):
+def _default_gains(num, den, xlim, ylim, zoom_xlim=None, zoom_ylim=None):
     """Unsupervised gains calculation for root locus plot.
 
     References
@@ -203,23 +281,16 @@ def _default_gains(num, den, xlim, ylim):
     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(
@@ -274,7 +345,7 @@ def _default_gains(num, den, xlim, ylim):
         tolerance = x_tolerance
     else:
         tolerance = np.min([x_tolerance, y_tolerance])
-    indexes_too_far = _indexes_filt(root_array, tolerance)
+    indexes_too_far = _indexes_filt(root_array, tolerance, zoom_xlim, zoom_ylim)
 
     # Add more points into the root locus for points that are too far apart
     while len(indexes_too_far) > 0 and kvect.size < 5000:
@@ -286,7 +357,7 @@ def _default_gains(num, den, xlim, ylim):
             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)
+        indexes_too_far = _indexes_filt(root_array, tolerance, zoom_xlim, zoom_ylim)
 
     new_gains = kvect[-1] * np.hstack((np.logspace(0, 3, 4)))
     new_points = _RLFindRoots(num, den, new_gains[1:4])
@@ -296,8 +367,8 @@ def _default_gains(num, den, xlim, ylim):
     return kvect, root_array, xlim, ylim
 
 
-def _indexes_filt(root_array, tolerance):
-    """Calculate the distance between points and return the indices.
+def _indexes_filt(root_array, tolerance, zoom_xlim=None, zoom_ylim=None):
+    """Calculate the distance between points and return the indexes.
 
     Filter the indexes so only the resolution of points within the xlim and
     ylim is improved when zoom is used.
@@ -305,6 +376,48 @@ def _indexes_filt(root_array, tolerance):
     """
     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
 
@@ -445,6 +558,249 @@ def _RLSortRoots(roots):
     return sorted
 
 
-# Alternative ways to call these functions
-root_locus = root_locus_plot
-rlocus = root_locus_plot
+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
diff --git a/control/robust.py b/control/robust.py
index a0e53d19..75930e59 100644
--- a/control/robust.py
+++ b/control/robust.py
@@ -41,6 +41,7 @@
 
 # External packages and modules
 import numpy as np
+import warnings
 from .exception import *
 from .statesp import StateSpace
 from .statefbk import *
@@ -357,7 +358,12 @@ def augw(g, w1=None, w2=None, w3=None):
     # output indices
     oi = np.arange(1, 1 + now1 + now2 + now3 + ny)
 
-    p = connect(sysall, q, ii, oi)
+    # Filter out known warning due to use of connect
+    with warnings.catch_warnings():
+        warnings.filterwarnings(
+            'ignore', message="`connect`", category=DeprecationWarning)
+
+        p = connect(sysall, q, ii, oi)
 
     return p
 
diff --git a/control/sisotool.py b/control/sisotool.py
index 82a989ca..79069162 100644
--- a/control/sisotool.py
+++ b/control/sisotool.py
@@ -1,23 +1,19 @@
 __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 .nlsys import interconnect
-from .statesp import ss, summing_junction
-from .timeresp import step_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
 
 _sisotool_defaults = {
     'sisotool.initial_gain': 1
@@ -90,7 +86,7 @@ def sisotool(sys, initial_gain=None, xlim_rlocus=None, ylim_rlocus=None,
     >>> ct.sisotool(G)                                          # doctest: +SKIP
 
     """
-    from .rlocus import root_locus_map
+    from .rlocus import root_locus
 
     # sys as loop transfer function if SISO
     if not sys.issiso():
@@ -127,51 +123,13 @@ 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)
 
-    # 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()
-
+    # 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)
 
 def _SisotoolUpdate(sys, fig, K, bode_plot_params, tvect=None):
 
@@ -228,7 +186,11 @@ def _SisotoolUpdate(sys, fig, K, bode_plot_params, tvect=None):
         sys_closed = append(sys, -K)
         connects = [[1, 3],
                     [3, 1]]
-        sys_closed = connect(sys_closed, connects, 2, 2)
+        # Filter out known warning due to use of connect
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                'ignore', message="`connect`", category=DeprecationWarning)
+            sys_closed = connect(sys_closed, connects, 2, 2)
     if tvect is None:
         tvect, yout = step_response(sys_closed, T_num=100)
     else:
diff --git a/control/stochsys.py b/control/stochsys.py
index b1bb6ef4..fe11a4fb 100644
--- a/control/stochsys.py
+++ b/control/stochsys.py
@@ -183,14 +183,14 @@ def lqe(*args, **kwargs):
 
 # contributed by Sawyer B. Fuller <minster@uw.edu>
 def dlqe(*args, **kwargs):
-    """dlqe(A, G, C, QN, RN, [, N])
+    r"""dlqe(A, G, C, QN, RN, [, N])
 
     Linear quadratic estimator design (Kalman filter) for discrete-time
     systems. Given the system
 
     .. math::
 
-        x[n+1] &= Ax[n] + Bu[n] + Gw[n] \\\\
+        x[n+1] &= Ax[n] + Bu[n] + Gw[n] \\
         y[n] &= Cx[n] + Du[n] + v[n]
 
     with unbiased process noise w and measurement noise v with covariances
diff --git a/control/timeplot.py b/control/timeplot.py
index 169e7343..58f7d838 100644
--- a/control/timeplot.py
+++ b/control/timeplot.py
@@ -292,7 +292,7 @@ def time_response_plot(
     if ax is None:
         with plt.rc_context(timeplot_rcParams):
             ax_array = fig.subplots(nrows, ncols, sharex=True, squeeze=False)
-            fig.set_tight_layout(True)
+            fig.set_layout_engine('tight')
             fig.align_labels()
 
     else:
diff --git a/control/timeresp.py b/control/timeresp.py
index 3822f13f..f431e6b3 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -678,15 +678,12 @@ class TimeResponseData:
 
         # Create a dict for setting up the data frame
         data = {'time': self.time}
-        if self.ninputs > 0:
-            data.update(
-                {name: self.u[i] for i, name in enumerate(self.input_labels)})
-        if self.noutputs > 0:
-            data.update(
-                {name: self.y[i] for i, name in enumerate(self.output_labels)})
-        if self.nstates > 0:
-            data.update(
-                {name: self.x[i] for i, name in enumerate(self.state_labels)})
+        data.update(
+            {name: self.u[i] for i, name in enumerate(self.input_labels)})
+        data.update(
+            {name: self.y[i] for i, name in enumerate(self.output_labels)})
+        data.update(
+            {name: self.x[i] for i, name in enumerate(self.state_labels)})
 
         return pandas.DataFrame(data)
 
@@ -1087,7 +1084,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
             if U.ndim == 1:
                 U = U.reshape(1, -1)  # pylint: disable=E1103
 
-        # Algorithm: to integrate from time 0 to time dt, with linear
+            # Algorithm: to integrate from time 0 to time dt, with linear
             # interpolation between inputs u(0) = u0 and u(dt) = u1, we solve
             #   xdot = A x + B u,        x(0) = x0
             #   udot = (u1 - u0) / dt,   u(0) = u0.
diff --git a/control/xferfcn.py b/control/xferfcn.py
index 206f0bc9..099f6425 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -468,7 +468,13 @@ class TransferFunction(LTI):
                     numstr = _tf_polynomial_to_string(self.num[no][ni], var=var)
                     denstr = _tf_polynomial_to_string(self.den[no][ni], var=var)
                 elif self.display_format == 'zpk':
-                    z, p, k = tf2zpk(self.num[no][ni], self.den[no][ni])
+                    num = self.num[no][ni]
+                    if num.size == 1 and num.item() == 0:
+                        # Catch a special case that SciPy doesn't handle
+                        z, p, k = tf2zpk([1.], self.den[no][ni])
+                        k = 0
+                    else:
+                        z, p, k = tf2zpk(self.num[no][ni], self.den[no][ni])
                     numstr = _tf_factorized_polynomial_to_string(
                         z, gain=k, var=var)
                     denstr = _tf_factorized_polynomial_to_string(p, var=var)
diff --git a/doc/Makefile b/doc/Makefile
index 71e493f2..a5f7ec5a 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 rlocus-siso_ctime-default.png
+  freqplot-siso_bode-default.png
 classes.pdf: classes.fig
 	fig2dev -Lpdf $< $@
 
@@ -26,9 +26,6 @@ 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 2f6857c3..be7ae7a5 100644
--- a/doc/plotting.rst
+++ b/doc/plotting.rst
@@ -4,15 +4,15 @@
 Plotting data
 *************
 
-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::
+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::
 
   bode_plot(sys)
   nyquist_plot([sys1, sys2])
-  pole_zero_plot(sys)
-  root_locus_plot(sys)
+
+.. root_locus_plot(sys)		# not yet implemented
 
 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()
-  root_locus_map(sys).plot()
+  rootlocus_response(sys).plot()     # implementation pending
 
 The remainder of this chapter provides additional documentation on how
 these response and plotting functions can be customized.
@@ -222,58 +222,6 @@ 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
 ===============================
 
@@ -296,8 +244,6 @@ 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
 
@@ -310,8 +256,6 @@ 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
 
@@ -340,5 +284,4 @@ 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/flatsys_test.py b/control/tests/flatsys_test.py
index 10a5d1bc..a12bf148 100644
--- a/control/tests/flatsys_test.py
+++ b/control/tests/flatsys_test.py
@@ -194,14 +194,17 @@ class TestFlatSys:
         else:
             initial_guess = None
 
-        # Solve the optimal trajectory
-        traj_ocp = fs.solve_flat_ocp(
-            vehicle_flat, timepts, x0, u0,
-            cost=traj_cost, constraints=input_constraints,
-            terminal_cost=terminal_cost, basis=basis,
-            initial_guess=initial_guess,
-            minimize_kwargs={'method': method},
-        )
+        # Solve the optimal trajectory (allow warnings)
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                'ignore', message="unable to solve", category=UserWarning)
+            traj_ocp = fs.solve_flat_ocp(
+                vehicle_flat, timepts, x0, u0,
+                cost=traj_cost, constraints=input_constraints,
+                terminal_cost=terminal_cost, basis=basis,
+                initial_guess=initial_guess,
+                minimize_kwargs={'method': method},
+            )
         xd, ud = traj_ocp.eval(timepts)
 
         if not traj_ocp.success:
diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py
index f452fe7d..dd236769 100644
--- a/control/tests/freqresp_test.py
+++ b/control/tests/freqresp_test.py
@@ -16,7 +16,8 @@ import control as ctrl
 from control.statesp import StateSpace
 from control.xferfcn import TransferFunction
 from control.matlab import ss, tf, bode, rss
-from control.freqplot import bode_plot, nyquist_plot, singular_values_plot
+from control.freqplot import bode_plot, nyquist_plot, nyquist_response, \
+    singular_values_plot, singular_values_response
 from control.tests.conftest import slycotonly
 
 pytestmark = pytest.mark.usefixtures("mplcleanup")
@@ -97,9 +98,8 @@ def test_nyquist_basic(ss_siso):
     tf_siso = tf(ss_siso)
     nyquist_plot(ss_siso)
     nyquist_plot(tf_siso)
-    count, contour = nyquist_plot(
-        tf_siso, plot=False, return_contour=True, omega_num=20)
-    assert len(contour) == 20
+    response = nyquist_response(tf_siso, omega_num=20)
+    assert len(response.contour) == 20
 
     with pytest.warns(UserWarning, match="encirclements was a non-integer"):
         count, contour = nyquist_plot(
@@ -107,9 +107,8 @@ def test_nyquist_basic(ss_siso):
     assert_allclose(contour[0], 1j)
     assert_allclose(contour[-1], 100j)
 
-    count, contour = nyquist_plot(
-        tf_siso, plot=False, omega=np.logspace(-1, 1, 10), return_contour=True)
-    assert len(contour) == 10
+    response = nyquist_response(tf_siso, omega=np.logspace(-1, 1, 10))
+    assert len(response.contour) == 10
 
 
 @pytest.mark.usefixtures("legacy_plot_signature")
@@ -200,7 +199,7 @@ def test_bode_margin(dB, maginfty1, maginfty2, gminv,
     den = [1, 25, 100, 0]
     sys = ctrl.tf(num, den)
     plt.figure()
-    ctrl.bode_plot(sys, margins=True, dB=dB, deg=deg, Hz=Hz)
+    ctrl.bode_plot(sys, display_margins=True, dB=dB, deg=deg, Hz=Hz)
     fig = plt.gcf()
     allaxes = fig.get_axes()
 
@@ -655,7 +654,8 @@ def tsystem(request, ss_mimo_ct, ss_miso_ct, ss_simo_ct, ss_siso_ct, ss_mimo_dt)
 def test_singular_values_plot(tsystem):
     sys = tsystem.sys
     for omega_ref, sigma_ref in zip(tsystem.omegas, tsystem.sigmas):
-        sigma, _ = singular_values_plot(sys, omega_ref, plot=False)
+        response = singular_values_response(sys, omega_ref)
+        sigma = np.real(response.fresp[:, 0, :])
         np.testing.assert_almost_equal(sigma, sigma_ref)
 
 
@@ -663,13 +663,13 @@ def test_singular_values_plot_mpl_base(ss_mimo_ct, ss_mimo_dt):
     sys_ct = ss_mimo_ct.sys
     sys_dt = ss_mimo_dt.sys
     plt.figure()
-    singular_values_plot(sys_ct, plot=True)
+    singular_values_plot(sys_ct)
     fig = plt.gcf()
     allaxes = fig.get_axes()
     assert(len(allaxes) == 1)
     assert(allaxes[0].get_label() == 'control-sigma')
     plt.figure()
-    singular_values_plot([sys_ct, sys_dt], plot=True, Hz=True, dB=True, grid=False)
+    singular_values_plot([sys_ct, sys_dt], Hz=True, dB=True, grid=False)
     fig = plt.gcf()
     allaxes = fig.get_axes()
     assert(len(allaxes) == 1)
@@ -679,10 +679,10 @@ def test_singular_values_plot_mpl_base(ss_mimo_ct, ss_mimo_dt):
 def test_singular_values_plot_mpl_superimpose_nyq(ss_mimo_ct, ss_mimo_dt):
     sys_ct = ss_mimo_ct.sys
     sys_dt = ss_mimo_dt.sys
-    omega_all = np.logspace(-3, 2, 1000)
+    omega_all = np.logspace(-3, int(math.log10(2 * math.pi/sys_dt.dt)), 1000)
     plt.figure()
-    singular_values_plot(sys_ct, omega_all, plot=True)
-    singular_values_plot(sys_dt, omega_all, plot=True)
+    singular_values_plot(sys_ct, omega_all)
+    singular_values_plot(sys_dt, omega_all)
     fig = plt.gcf()
     allaxes = fig.get_axes()
     assert(len(allaxes) == 1)
diff --git a/control/tests/kwargs_test.py b/control/tests/kwargs_test.py
index 53cc0076..0d13e639 100644
--- a/control/tests/kwargs_test.py
+++ b/control/tests/kwargs_test.py
@@ -97,7 +97,6 @@ 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}),
@@ -176,8 +175,6 @@ 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
@@ -196,18 +193,16 @@ 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, TypeError),
-            match="(has no property|unexpected keyword|unrecognized keyword)"):
+    with pytest.raises(AttributeError,
+                       match="(has no property|unexpected 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, TypeError),
-            match="(has no property|unexpected keyword|unrecognized keyword)"):
+    with pytest.raises(AttributeError,
+                       match="(has no property|unexpected keyword)"):
         response.plot(unknown=None)
 
 #
@@ -246,10 +241,8 @@ 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,
@@ -281,7 +274,6 @@ 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 2ba3d5df..e01abcca 100644
--- a/control/tests/matlab_test.py
+++ b/control/tests/matlab_test.py
@@ -424,8 +424,7 @@ class TestMatlab:
     @pytest.mark.parametrize("subsys", ["ss1", "tf1", "tf2"])
     def testRlocus(self, siso, subsys, mplcleanup):
         """Call rlocus()"""
-        rlist, klist = rlocus(getattr(siso, subsys))
-        np.testing.assert_equal(len(rlist), len(klist))
+        rlocus(getattr(siso, subsys))
 
     def testRlocus_list(self, siso, mplcleanup):
         """Test rlocus() with list"""
diff --git a/control/tests/optimal_test.py b/control/tests/optimal_test.py
index d4a2b1d8..f746db7d 100644
--- a/control/tests/optimal_test.py
+++ b/control/tests/optimal_test.py
@@ -745,12 +745,15 @@ def test_optimal_doc(method, npts, initial_guess, fail):
         initial_guess = (state_guess, input_guess)
 
     # Solve the optimal control problem
-    result = opt.solve_ocp(
-        vehicle, timepts, x0, traj_cost, constraints,
-        terminal_cost=term_cost, initial_guess=initial_guess,
-        trajectory_method=method,
-        # minimize_method='COBYLA', # SLSQP',
-    )
+    with warnings.catch_warnings():
+        warnings.filterwarnings(
+            'ignore', message="unable to solve", category=UserWarning)
+        result = opt.solve_ocp(
+            vehicle, timepts, x0, traj_cost, constraints,
+            terminal_cost=term_cost, initial_guess=initial_guess,
+            trajectory_method=method,
+            # minimize_method='COBYLA', # SLSQP',
+        )
 
     if fail == 'xfail':
         assert not result.success
diff --git a/control/tests/pzmap_test.py b/control/tests/pzmap_test.py
index ce8adf6e..8d41807b 100644
--- a/control/tests/pzmap_test.py
+++ b/control/tests/pzmap_test.py
@@ -12,11 +12,9 @@ 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"),
@@ -46,23 +44,20 @@ def test_pzmap(kwargs, setdefaults, dt, editsdefaults, mplcleanup):
 
     pzkwargs = kwargs.copy()
     if setdefaults:
-        for k in ['grid']:
+        for k in ['plot', '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):
-        fig, ax = plt.gcf(), plt.gca()
+        ax = plt.gca()
 
-        assert fig._suptitle.get_text().startswith(
-            kwargs.get('title', 'Pole/zero plot'))
+        assert ax.get_title() == kwargs.get('title', 'Pole Zero Map')
 
         # FIXME: This won't work when zgrid and sgrid are unified
         children = ax.get_children()
@@ -83,43 +78,12 @@ def test_pzmap(kwargs, setdefaults, dt, editsdefaults, mplcleanup):
         assert not plt.get_fignums()
 
 
-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_warns():
+    with pytest.warns(FutureWarning):
+        pzmap(TransferFunction([1], [1, 2]), Plot=True)
 
 
 def test_pzmap_raises():
     with pytest.raises(TypeError):
         # not an LTI system
-        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)
+        pzmap(([1], [1,2]))
diff --git a/control/tests/rlocus_test.py b/control/tests/rlocus_test.py
index 5ee2830c..e61f0c8f 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
+from control.rlocus import root_locus, _RLClickDispatcher
 from control.xferfcn import TransferFunction
 from control.statesp import StateSpace
 from control.bdalg import feedback
@@ -45,7 +45,6 @@ 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]
@@ -56,75 +55,50 @@ class TestRootLocus:
         self.check_cl_poles(sys, roots, klist)
 
         # now check with plotting
-        roots, k_out = root_locus(sys, klist, plot=True)
+        roots, k_out = root_locus(sys, klist)
         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.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
+    @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)
         ax = plt.gca()
-        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
-            assert not isinstance(ax, AA.Axes)
+        n_gridlines = sum([int(line.get_linestyle() in [':', 'dotted',
+                                                        '--', 'dashed'])
+                           for line in ax.lines])
+        if grid is False:
+            assert n_gridlines == 2
         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
+            assert n_gridlines > 2
+        # TODO check validity of grid
 
-        # TODO: check validity of grid
+    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)
 
-    @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]), plot=True)
+            root_locus(TransferFunction([-1, 2], [1, 2]))
 
     # 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, plot=True)
+        root_locus(system)
         fig = plt.gcf()
         ax_rlocus = fig.axes[0]
 
@@ -147,7 +121,6 @@ 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"""
@@ -167,104 +140,4 @@ class TestRootLocus:
         sys = ct.tf(*sp.signal.zpk2tf(
             [-1e-2, 1-1e7j, 1+1e7j], [0, -1e7j, 1e7j], 1))
 
-        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)
+        ct.root_locus(sys)
diff --git a/control/tests/sisotool_test.py b/control/tests/sisotool_test.py
index b5aaeb90..3dd9168a 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.sisotool import _click_dispatcher
+from control.rlocus import _RLClickDispatcher
 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[4].get_data(),
+        assert_array_almost_equal(ax_rlocus.lines[0].get_data(),
                                   initial_point_0, 4)
-        assert_array_almost_equal(ax_rlocus.lines[5].get_data(),
+        assert_array_almost_equal(ax_rlocus.lines[1].get_data(),
                                   initial_point_1, 4)
-        assert_array_almost_equal(ax_rlocus.lines[6].get_data(),
+        assert_array_almost_equal(ax_rlocus.lines[2].get_data(),
                                   initial_point_2, 4)
 
         # Check the step response before moving the point
@@ -79,8 +79,7 @@ class TestSisotool:
             'omega_limits': None,
             'omega_num': None,
             'ax': np.array([[ax_mag], [ax_phase]]),
-            'margins': True,
-            'margin_info': True,
+            'display_margins': 'overlay',
         }
 
         # Check that the xaxes of the bode plot are shared before the rlocus click
@@ -93,8 +92,9 @@ class TestSisotool:
         event = type('test', (object,), {'xdata': 2.31206868287,
                                          'ydata': 15.5983051046,
                                          'inaxes': ax_rlocus.axes})()
-        _click_dispatcher(event=event, sys=tsys, ax=ax_rlocus,
-                          bode_plot_params=bode_plot_params, tvect=None)
+        _RLClickDispatcher(event=event, sys=tsys, fig=fig,
+                           ax_rlocus=ax_rlocus, sisotool=True, plotstr='-',
+                           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.]))
@@ -142,8 +142,9 @@ class TestSisotool:
         event = type('test', (object,), {'xdata': 2.31206868287,
                                          'ydata': 15.5983051046,
                                          'inaxes': ax_rlocus.axes})()
-        _click_dispatcher(event=event, sys=tsys, ax=ax_rlocus,
-                          bode_plot_params=dict(), tvect=tvect)
+        _RLClickDispatcher(event=event, sys=tsys, fig=fig,
+                           ax_rlocus=ax_rlocus, sisotool=True, plotstr='-',
+                           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,
@@ -200,21 +201,3 @@ 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)
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index fb21180b..989c9bd9 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -1230,14 +1230,6 @@ def test_to_pandas():
     np.testing.assert_equal(df['x[0]'], resp.states[0])
     np.testing.assert_equal(df['x[1]'], resp.states[1])
 
-    # System with no states
-    sys = ct.ss([], [], [], 5)
-    resp = ct.input_output_response(sys, timepts, np.sin(timepts), t_eval=T)
-    df = resp.to_pandas()
-    np.testing.assert_equal(df['time'], resp.time)
-    np.testing.assert_equal(df['u[0]'], resp.inputs)
-    np.testing.assert_equal(df['y[0]'], resp.inputs * 5)
-
 
 @pytest.mark.skipif(pandas_check(), reason="pandas installed")
 def test_no_pandas():