From 1a6d80668b3ea81329f0e474400be1a84205fd86 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Fri, 14 Aug 2020 13:45:03 -0700
Subject: [PATCH 01/32] eliminate _evalfr in lti system classes, replaced with
 __call__, which evaluates at many frequencies at once

---
 control/bench/time_freqresp.py      |   4 +-
 control/frdata.py                   | 168 +++++++++++---------------
 control/freqplot.py                 |  14 +--
 control/lti.py                      |  74 ++++++++++--
 control/margins.py                  |  16 +--
 control/nichols.py                  |   2 +-
 control/statesp.py                  | 147 ++++++++---------------
 control/tests/frd_test.py           | 179 ++++++++++++----------------
 control/tests/freqresp_test.py      |  12 +-
 control/tests/statesp_array_test.py |  23 +---
 control/tests/statesp_test.py       |  39 ++----
 control/tests/xferfcn_test.py       |  38 ++----
 control/xferfcn.py                  | 110 +++--------------
 doc/control.rst                     |   1 +
 examples/robust_mimo.py             |   2 +-
 examples/robust_siso.py             |   8 +-
 16 files changed, 333 insertions(+), 504 deletions(-)

diff --git a/control/bench/time_freqresp.py b/control/bench/time_freqresp.py
index 1945cbc24..3ae837082 100644
--- a/control/bench/time_freqresp.py
+++ b/control/bench/time_freqresp.py
@@ -8,7 +8,7 @@
 sys_tf = tf(sys)
 w = logspace(-1,1,50)
 ntimes = 1000
-time_ss = timeit("sys.freqresp(w)", setup="from __main__ import sys, w", number=ntimes)
-time_tf = timeit("sys_tf.freqresp(w)", setup="from __main__ import sys_tf, w", number=ntimes)
+time_ss = timeit("sys.freqquency_response(w)", setup="from __main__ import sys, w", number=ntimes)
+time_tf = timeit("sys_tf.frequency_response(w)", setup="from __main__ import sys_tf, w", number=ntimes)
 print("State-space model on %d states: %f" % (nstates, time_ss))
 print("Transfer-function model on %d states: %f" % (nstates, time_tf))
diff --git a/control/frdata.py b/control/frdata.py
index 8ca9dfd9d..084f4aba3 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -47,8 +47,8 @@
 # External function declarations
 from warnings import warn
 import numpy as np
-from numpy import angle, array, empty, ones, \
-    real, imag, absolute, eye, linalg, where, dot
+from numpy import angle, array, empty, ones, isin, \
+    real, imag, absolute, eye, linalg, where, dot, sort
 from scipy.interpolate import splprep, splev
 from .lti import LTI
 
@@ -107,17 +107,10 @@ def __init__(self, *args, **kwargs):
                 # not an FRD, but still a system, second argument should be
                 # the frequency range
                 otherlti = args[0]
-                self.omega = array(args[1], dtype=float)
-                self.omega.sort()
+                self.omega = sort(np.asarray(args[1], dtype=float))
                 numfreq = len(self.omega)
-
                 # calculate frequency response at my points
-                self.fresp = empty(
-                    (otherlti.outputs, otherlti.inputs, numfreq),
-                    dtype=complex)
-                for k, w in enumerate(self.omega):
-                    self.fresp[:, :, k] = otherlti._evalfr(w)
-
+                self.fresp = otherlti(1j * self.omega, squeeze=False)
             else:
                 # The user provided a response and a freq vector
                 self.fresp = array(args[0], dtype=complex)
@@ -141,7 +134,7 @@ def __init__(self, *args, **kwargs):
             self.omega = args[0].omega
             self.fresp = args[0].fresp
         else:
-            raise ValueError("Needs 1 or 2 arguments; receivd %i." % len(args))
+            raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
 
         # create interpolation functions
         if smooth:
@@ -163,7 +156,7 @@ def __str__(self):
         mimo = self.inputs > 1 or self.outputs > 1
         outstr = ['Frequency response data']
 
-        mt, pt, wt = self.freqresp(self.omega)
+        #mt, pt, wt = self.frequency_response(self.omega)
         for i in range(self.inputs):
             for j in range(self.outputs):
                 if mimo:
@@ -173,7 +166,7 @@ def __str__(self):
                 outstr.extend(
                     ['%12.3f  %10.4g%+10.4gj' % (w, m, p)
                      for m, p, w in zip(real(self.fresp[j, i, :]),
-                                        imag(self.fresp[j, i, :]), wt)])
+                                        imag(self.fresp[j, i, :]), self.omega)])
 
         return '\n'.join(outstr)
 
@@ -342,28 +335,14 @@ def __pow__(self, other):
             return (FRD(ones(self.fresp.shape), self.omega) / self) * \
                 (self**(other+1))
 
-    def evalfr(self, omega):
-        """Evaluate a transfer function at a single angular frequency.
-
-        self._evalfr(omega) returns the value of the frequency response
-        at frequency omega.
-
-        Note that a "normal" FRD only returns values for which there is an
-        entry in the omega vector. An interpolating FRD can return
-        intermediate values.
-
-        """
-        warn("FRD.evalfr(omega) will be deprecated in a future release "
-             "of python-control; use sys.eval(omega) instead",
-             PendingDeprecationWarning)         # pragma: no coverage
-        return self._evalfr(omega)
-
     # Define the `eval` function to evaluate an FRD at a given (real)
     # frequency.  Note that we choose to use `eval` instead of `evalfr` to
     # avoid confusion with :func:`evalfr`, which takes a complex number as its
     # argument.  Similarly, we don't use `__call__` to avoid confusion between
     # G(s) for a transfer function and G(omega) for an FRD object.
-    def eval(self, omega):
+    # update Sawyer B. Fuller 2020.08.14: __call__ added to provide a uniform 
+    # interface to systems in general and the lti.frequency_response method
+    def eval(self, omega, squeeze=True):
         """Evaluate a transfer function at a single angular frequency.
 
         self.evalfr(omega) returns the value of the frequency response
@@ -371,81 +350,69 @@ def eval(self, omega):
 
         Note that a "normal" FRD only returns values for which there is an
         entry in the omega vector. An interpolating FRD can return
-        intermediate values.
+        intermediate values.     
+        
+        squeeze: bool, optional (default=True)
+            If True and sys is single input, single output (SISO), return a 
+            1D array or scalar depending on omega's length.
 
         """
-        return self._evalfr(omega)
-
-    # Internal function to evaluate the frequency responses
-    def _evalfr(self, omega):
-        """Evaluate a transfer function at a single angular frequency."""
-        # Preallocate the output.
-        if getattr(omega, '__iter__', False):
-            out = empty((self.outputs, self.inputs, len(omega)), dtype=complex)
-        else:
-            out = empty((self.outputs, self.inputs), dtype=complex)
+        omega_array = np.array(omega, ndmin=1) # array-like version of omega 
+        if any(omega_array.imag > 0): 
+            raise ValueError("FRD.eval can only accept real-valued omega")
 
         if self.ifunc is None:
-            try:
-                out = self.fresp[:, :, where(self.omega == omega)[0][0]]
-            except Exception:
+            elements = isin(self.omega, omega)
+            if sum(elements) < len(omega_array):
                 raise ValueError(
-                    "Frequency %f not in frequency list, try an interpolating"
-                    " FRD if you want additional points" % omega)
-        else:
-            if getattr(omega, '__iter__', False):
-                for i in range(self.outputs):
-                    for j in range(self.inputs):
-                        for k, w in enumerate(omega):
-                            frraw = splev(w, self.ifunc[i, j], der=0)
-                            out[i, j, k] = frraw[0] + 1.0j * frraw[1]
+                    "not all frequencies omega are in frequency list of FRD "
+                    "system. Try an interpolating FRD for additional points.")
             else:
-                for i in range(self.outputs):
-                    for j in range(self.inputs):
-                        frraw = splev(omega, self.ifunc[i, j], der=0)
-                        out[i, j] = frraw[0] + 1.0j * frraw[1]
-
+                out = self.fresp[:, :, elements]
+        else:
+            out = empty((self.outputs, self.inputs, len(omega_array)))
+            for i in range(self.outputs):
+                for j in range(self.inputs):
+                    for k, w in enumerate(omega_array):
+                        frraw = splev(w, self.ifunc[i, j], der=0)
+                        out[i, j, k] = frraw[0] + 1.0j * frraw[1]
+        if not hasattr(omega, '__len__'):
+            # omega is a scalar, squeeze down array along last dim
+            out = np.squeeze(out, axis=2)
+        if squeeze and self.issiso():
+            out = out[0][0]
         return out
 
-    # Method for generating the frequency response of the system
-    def freqresp(self, omega):
-        """Evaluate the frequency response at a list of angular frequencies.
-
-        Reports the value of the magnitude, phase, and angular frequency of
-        the requency response evaluated at omega, where omega is a list of
-        angular frequencies, and is a sorted version of the input omega.
-
-        Parameters
-        ----------
-        omega : array_like
-            A list of frequencies in radians/sec at which the system should be
-            evaluated. The list can be either a python list or a numpy array
-            and will be sorted before evaluation.
-
-        Returns
-        -------
-        mag : (self.outputs, self.inputs, len(omega)) ndarray
-            The magnitude (absolute value, not dB or log10) of the system
-            frequency response.
-        phase : (self.outputs, self.inputs, len(omega)) ndarray
-            The wrapped phase in radians of the system frequency response.
-        omega : ndarray or list or tuple
-            The list of sorted frequencies at which the response was
-            evaluated.
-        """
-        # Preallocate outputs.
-        numfreq = len(omega)
-        mag = empty((self.outputs, self.inputs, numfreq))
-        phase = empty((self.outputs, self.inputs, numfreq))
-
-        omega.sort()
+    def __call__(self, s, squeeze=True):
+        """Evaluate the system's transfer function at complex frequencies s.
 
-        for k, w in enumerate(omega):
-            fresp = self._evalfr(w)
-            mag[:, :, k] = abs(fresp)
-            phase[:, :, k] = angle(fresp)
+        For a SISO system, returns the complex value of the
+        transfer function.  For a MIMO transfer fuction, returns a
+        matrix of values.
+        
+        Raises an error if s is not purely imaginary. 
+        
+        squeeze: bool, optional (default=True)
+            If True and sys is single input, single output (SISO), return a 
+            1D array or scalar depending on omega's length.
 
-        return mag, phase, omega
+        """
+        if any(abs(np.array(s, ndmin=1).real) > 0):
+            raise ValueError("__call__: FRD systems can only accept"
+                            "purely imaginary frequencies")
+        # need to preserve array or scalar status
+        if hasattr(s, '__len__'):
+            return self.eval(np.asarray(s).imag, squeeze=squeeze)
+        else:
+            return self.eval(complex(s).imag, squeeze=squeeze)
+        
+    def freqresp(self, omega):
+        warn("FrequencyResponseData.freqresp(omega) will be deprecated in a "
+             "future release of python-control; use "
+             "FrequencyResponseData.frequency_response(sys, omega), or "
+             "freqresp(sys, omega) in the MATLAB compatibility module "
+             "instead", PendingDeprecationWarning)        
+        return self.frequency_response(omega)
 
     def feedback(self, other=1, sign=-1):
         """Feedback interconnection between two FRD objects."""
@@ -515,11 +482,10 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1):
             "Frequency ranges of FRD do not match, conversion not implemented")
 
     elif isinstance(sys, LTI):
-        omega.sort()
-        fresp = empty((sys.outputs, sys.inputs, len(omega)), dtype=complex)
-        for k, w in enumerate(omega):
-            fresp[:, :, k] = sys._evalfr(w)
-
+        omega = np.sort(omega)
+        fresp = sys(1j * omega)
+        if len(fresp.shape) == 1: 
+            fresp = fresp[np.newaxis, np.newaxis, :]
         return FRD(fresp, omega, smooth=True)
 
     elif isinstance(sys, (int, float, complex, np.number)):
diff --git a/control/freqplot.py b/control/freqplot.py
index 7b296c111..ef63dd2fb 100644
--- a/control/freqplot.py
+++ b/control/freqplot.py
@@ -136,7 +136,7 @@ def bode_plot(syslist, omega=None,
     Notes
     -----
     1. Alternatively, you may use the lower-level method (mag, phase, freq)
-    = sys.freqresp(freq) to generate the frequency response for a system,
+    = sys.frequency_response(freq) to generate the frequency response for a system,
     but it returns a MIMO response.
 
     2. If a discrete time model is given, the frequency response is plotted
@@ -207,7 +207,7 @@ def bode_plot(syslist, omega=None,
             else:
                 nyquistfrq = None
             # Get the magnitude and phase of the system
-            mag_tmp, phase_tmp, omega_sys = sys.freqresp(omega_sys)
+            mag_tmp, phase_tmp, omega_sys = sys.frequency_response(omega_sys)
             mag = np.atleast_1d(np.squeeze(mag_tmp))
             phase = np.atleast_1d(np.squeeze(phase_tmp))
             phase = unwrap(phase)
@@ -524,7 +524,7 @@ def nyquist_plot(syslist, omega=None, plot=True, label_freq=0,
                 "Nyquist is currently only implemented for SISO systems.")
         else:
             # Get the magnitude and phase of the system
-            mag_tmp, phase_tmp, omega = sys.freqresp(omega)
+            mag_tmp, phase_tmp, omega = sys.frequency_response(omega)
             mag = np.squeeze(mag_tmp)
             phase = np.squeeze(phase_tmp)
 
@@ -654,7 +654,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     omega_plot = omega / (2. * math.pi) if Hz else omega
 
     # TODO: Need to add in the mag = 1 lines
-    mag_tmp, phase_tmp, omega = S.freqresp(omega)
+    mag_tmp, phase_tmp, omega = S.frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['s'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -664,7 +664,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     plot_axes['s'].tick_params(labelbottom=False)
     plot_axes['s'].grid(grid, which='both')
 
-    mag_tmp, phase_tmp, omega = (P * S).freqresp(omega)
+    mag_tmp, phase_tmp, omega = (P * S).frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['ps'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -674,7 +674,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     plot_axes['ps'].set_ylabel("$|PS|$" + " (dB)" if dB else "")
     plot_axes['ps'].grid(grid, which='both')
 
-    mag_tmp, phase_tmp, omega = (C * S).freqresp(omega)
+    mag_tmp, phase_tmp, omega = (C * S).frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['cs'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -685,7 +685,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     plot_axes['cs'].set_ylabel("$|CS|$" + " (dB)" if dB else "")
     plot_axes['cs'].grid(grid, which='both')
 
-    mag_tmp, phase_tmp, omega = T.freqresp(omega)
+    mag_tmp, phase_tmp, omega = T.frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['t'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
diff --git a/control/lti.py b/control/lti.py
index 8db14794b..6d39955bb 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -13,7 +13,8 @@
 """
 
 import numpy as np
-from numpy import absolute, real
+from numpy import absolute, real, angle, abs
+from warnings import warn
 
 __all__ = ['issiso', 'timebase', 'timebaseEqual', 'isdtime', 'isctime',
            'pole', 'zero', 'damp', 'evalfr', 'freqresp', 'dcgain']
@@ -109,10 +110,55 @@ def damp(self):
         Z = -real(splane_poles)/wn
         return wn, Z, poles
 
+    def frequency_response(self, omega, squeeze=True):
+        """Evaluate the linear time-invariant system at an array of angular 
+        frequencies.
+
+        Reports the frequency response of the system,
+
+             G(j*omega) = mag*exp(j*phase)
+
+        for continuous time systems. For discrete time systems, the response is
+        evaluated around the unit circle such that
+
+             G(exp(j*omega*dt)) = mag*exp(j*phase).
+
+        Parameters
+        ----------
+        omega : array_like or float
+            A list, tuple, array, or scalar value of frequencies in 
+            radians/sec at which the system will be evaluated. 
+        squeeze: bool, optional (default=True)
+            If True and sys is single input, single output (SISO), return a 
+            1D array or scalar depending on omega's length.
+
+        Returns
+        -------
+        mag : (self.outputs, self.inputs, len(omega)) ndarray
+            The magnitude (absolute value, not dB or log10) of the system
+            frequency response.
+        phase : (self.outputs, self.inputs, len(omega)) ndarray
+            The wrapped phase in radians of the system frequency response.
+        omega : ndarray
+            The (sorted) frequencies at which the response was evaluated.
+        
+        """        
+        omega = np.sort(np.array(omega, ndmin=1))
+        if isdtime(self, strict=True):
+            # Convert the frequency to discrete time
+            if np.any(omega * self.dt > np.pi):
+                warn("__call__: evaluation above Nyquist frequency")
+            s = np.exp(1j * omega * self.dt)
+        else:
+            s = 1j * omega
+        response = self.__call__(s, squeeze=squeeze)
+        return abs(response), angle(response), omega
+
     def dcgain(self):
         """Return the zero-frequency gain"""
         raise NotImplementedError("dcgain not implemented for %s objects" %
                                   str(self.__class__))
+    
 
 # Test to see if a system is SISO
 def issiso(sys, strict=False):
@@ -379,20 +425,22 @@ def damp(sys, doprint=True):
                       (p.real, p.imag, d, w))
     return wn, damping, poles
 
-def evalfr(sys, x):
+def evalfr(sys, x, squeeze=True):
     """
-    Evaluate the transfer function of an LTI system for a single complex
-    number x.
+    Evaluate the transfer function of an LTI system for complex frequency x.
 
     To evaluate at a frequency, enter x = omega*j, where omega is the
-    frequency in radians
+    frequency in radians per second
 
     Parameters
     ----------
     sys: StateSpace or TransferFunction
         Linear system
-    x: scalar
+    x: scalar or array-like
         Complex number
+    squeeze: bool, optional (default=True)
+        If True and sys is single input, single output (SISO), return a 
+        1D array or scalar depending on omega's length.
 
     Returns
     -------
@@ -417,12 +465,12 @@ def evalfr(sys, x):
 
     .. todo:: Add example with MIMO system
     """
-    if issiso(sys):
+    if squeeze and issiso(sys):
         return sys.horner(x)[0][0]
     return sys.horner(x)
 
 
-def freqresp(sys, omega):
+def freqresp(sys, omega, squeeze=True):
     """
     Frequency response of an LTI system at multiple angular frequencies.
 
@@ -434,6 +482,9 @@ def freqresp(sys, omega):
         A list of frequencies in radians/sec at which the system should be
         evaluated. The list can be either a python list or a numpy array
         and will be sorted before evaluation.
+    squeeze: bool, optional (default=True)
+        If True and sys is single input, single output (SISO), return a 
+        1D array or scalar depending on omega's length.
 
     Returns
     -------
@@ -453,9 +504,8 @@ def freqresp(sys, omega):
 
     Notes
     -----
-    This function is a wrapper for StateSpace.freqresp and
-    TransferFunction.freqresp.  The output omega is a sorted version of the
-    input omega.
+    This function is a wrapper for StateSpace.frequency_response and
+    TransferFunction.frequency_response.  
 
     Examples
     --------
@@ -480,7 +530,7 @@ def freqresp(sys, omega):
         #>>> # frequency response from the 1st input to the 2nd output, for
         #>>> # s = 0.1i, i, 10i.
     """
-    return sys.freqresp(omega)
+    return sys.frequency_response(omega, squeeze=squeeze)
 
 
 def dcgain(sys):
diff --git a/control/margins.py b/control/margins.py
index 7bdcf6caa..4094b43bf 100644
--- a/control/margins.py
+++ b/control/margins.py
@@ -213,15 +213,15 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
         # a bit coarse, have the interpolated frd evaluated again
         def _mod(w):
             """Calculate |G(jw)| - 1"""
-            return np.abs(sys._evalfr(w)[0][0]) - 1
+            return np.abs(sys(1j * w)[0][0]) - 1
 
         def _arg(w):
             """Calculate the phase angle at -180 deg"""
-            return np.angle(-sys._evalfr(w)[0][0])
+            return np.angle(-sys(1j * w)[0][0])
 
         def _dstab(w):
             """Calculate the distance from -1 point"""
-            return np.abs(sys._evalfr(w)[0][0] + 1.)
+            return np.abs(sys(1j * w)[0][0] + 1.)
 
         # Find all crossings, note that this depends on omega having
         # a correct range
@@ -232,7 +232,7 @@ def _dstab(w):
 
         # find the phase crossings ang(H(jw) == -180
         widx = np.where(np.diff(np.sign(_arg(sys.omega))))[0]
-        widx = widx[np.real(sys._evalfr(sys.omega[widx])[0][0]) <= 0]
+        widx = widx[np.real(sys(1j * sys.omega[widx])[0][0]) <= 0]
         w_180 = np.array(
             [sp.optimize.brentq(_arg, sys.omega[i], sys.omega[i+1])
              for i in widx])
@@ -249,10 +249,10 @@ def _dstab(w):
     # margins, as iterables, converted frdata and xferfcn calculations to
     # vector for this
     with np.errstate(all='ignore'):
-        gain_w_180 = np.abs(sys._evalfr(w_180)[0][0])
+        gain_w_180 = np.abs(sys(1j * w_180))
         GM = 1.0/gain_w_180
-    SM = np.abs(sys._evalfr(wstab)[0][0]+1)
-    PM = np.remainder(np.angle(sys._evalfr(wc)[0][0], deg=True), 360.0) - 180.0
+    SM = np.abs(sys(1j * wstab)+1)
+    PM = np.remainder(np.angle(sys(1j * wc), deg=True), 360.0) - 180.0
 
     if returnall:
         return GM, PM, SM, w_180, wc, wstab
@@ -313,7 +313,7 @@ def phase_crossover_frequencies(sys):
 
     # using real() to avoid rounding errors and results like 1+0j
     # it would be nice to have a vectorized version of self.evalfr here
-    gain = np.real(np.asarray([tf._evalfr(f)[0][0] for f in realposfreq]))
+    gain = np.real(np.asarray([tf(1j * f)[0][0] for f in realposfreq]))
 
     return realposfreq, gain
 
diff --git a/control/nichols.py b/control/nichols.py
index c8a98ed5e..c4ac4de0c 100644
--- a/control/nichols.py
+++ b/control/nichols.py
@@ -96,7 +96,7 @@ def nichols_plot(sys_list, omega=None, grid=None):
 
     for sys in sys_list:
         # Get the magnitude and phase of the system
-        mag_tmp, phase_tmp, omega = sys.freqresp(omega)
+        mag_tmp, phase_tmp, omega = sys.frequency_response(omega)
         mag = np.squeeze(mag_tmp)
         phase = np.squeeze(phase_tmp)
 
diff --git a/control/statesp.py b/control/statesp.py
index 522d187a9..a4f3aec71 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -54,7 +54,7 @@
 import math
 import numpy as np
 from numpy import any, array, asarray, concatenate, cos, delete, \
-    dot, empty, exp, eye, isinf, ones, pad, sin, zeros, squeeze
+    dot, empty, exp, eye, isinf, ones, pad, sin, zeros, squeeze, pi
 from numpy.random import rand, randn
 from numpy.linalg import solve, eigvals, matrix_rank
 from numpy.linalg.linalg import LinAlgError
@@ -437,99 +437,27 @@ def __rdiv__(self, other):
 
         raise NotImplementedError("StateSpace.__rdiv__ is not implemented yet.")
 
-    def evalfr(self, omega):
-        """Evaluate a SS system's transfer function at a single frequency.
+    def __call__(self, s, squeeze=True):
+        """Evaluate system's transfer function at complex frequencies s/z.
+        
+        If squeeze is True (default) and sys is single input, single output 
+        (SISO), return a 1D array or scalar depending on s.
 
-        self._evalfr(omega) returns the value of the transfer function matrix
-        with input value s = i * omega.
-
-        """
-        warn("StateSpace.evalfr(omega) will be deprecated in a future "
-             "release of python-control; use evalfr(sys, omega*1j) instead",
-             PendingDeprecationWarning)
-        return self._evalfr(omega)
-
-    def _evalfr(self, omega):
-        """Evaluate a SS system's transfer function at a single frequency"""
-        # Figure out the point to evaluate the transfer function
-        if isdtime(self, strict=True):
-            dt = timebase(self)
-            s = exp(1.j * omega * dt)
-            if omega * dt > math.pi:
-                warn("_evalfr: frequency evaluation above Nyquist frequency")
-        else:
-            s = omega * 1.j
-
-        return self.horner(s)
-
-    def horner(self, s):
-        """Evaluate the systems's transfer function for a complex variable
-
-        Returns a matrix of values evaluated at complex variable s.
-        """
-        resp = np.dot(self.C, solve(s * eye(self.states) - self.A,
-                                    self.B)) + self.D
-        return array(resp)
-
-    def freqresp(self, omega):
-        """Evaluate the system's transfer function at a list of frequencies
-
-        Reports the frequency response of the system,
-
-             G(j*omega) = mag*exp(j*phase)
-
-        for continuous time. For discrete time systems, the response is
-        evaluated around the unit circle such that
-
-             G(exp(j*omega*dt)) = mag*exp(j*phase).
-
-        Parameters
-        ----------
-        omega : array_like
-            A list of frequencies in radians/sec at which the system should be
-            evaluated. The list can be either a python list or a numpy array
-            and will be sorted before evaluation.
-
-        Returns
-        -------
-        mag : (self.outputs, self.inputs, len(omega)) ndarray
-            The magnitude (absolute value, not dB or log10) of the system
-            frequency response.
-        phase : (self.outputs, self.inputs, len(omega)) ndarray
-            The wrapped phase in radians of the system frequency response.
-        omega : ndarray
-            The list of sorted frequencies at which the response was
-            evaluated.
-        """
-        # In case omega is passed in as a list, rather than a proper array.
-        omega = np.asarray(omega)
-
-        numFreqs = len(omega)
-        Gfrf = np.empty((self.outputs, self.inputs, numFreqs),
-                        dtype=np.complex128)
-
-        # Sort frequency and calculate complex frequencies on either imaginary
-        # axis (continuous time) or unit circle (discrete time).
-        omega.sort()
-        if isdtime(self, strict=True):
-            dt = timebase(self)
-            cmplx_freqs = exp(1.j * omega * dt)
-            if max(np.abs(omega)) * dt > math.pi:
-                warn("freqresp: frequency evaluation above Nyquist frequency")
-        else:
-            cmplx_freqs = omega * 1.j
-
-        # Do the frequency response evaluation. Use TB05AD from Slycot
-        # if it's available, otherwise use the built-in horners function.
+        """        
+        # Use TB05AD from Slycot if available
         try:
             from slycot import tb05ad
 
-            n = np.shape(self.A)[0]
+            # preallocate
+            s_arr = np.array(s, ndmin=1) # array-like version of s
+            out = np.empty((self.outputs, self.inputs, len(s_arr)), 
+                           dtype=complex)
+            n = self.states
             m = self.inputs
             p = self.outputs
             # The first call both evaluates C(sI-A)^-1 B and also returns
             # Hessenberg transformed matrices at, bt, ct.
-            result = tb05ad(n, m, p, cmplx_freqs[0], self.A,
+            result = tb05ad(n, m, p, s_arr[0], self.A,
                             self.B, self.C, job='NG')
             # When job='NG', result = (at, bt, ct, g_i, hinvb, info)
             at = result[0]
@@ -537,27 +465,54 @@ def freqresp(self, omega):
             ct = result[2]
 
             # TB05AD frequency evaluation does not include direct feedthrough.
-            Gfrf[:, :, 0] = result[3] + self.D
+            out[:, :, 0] = result[3] + self.D
 
             # Now, iterate through the remaining frequencies using the
             # transformed state matrices, at, bt, ct.
 
             # Start at the second frequency, already have the first.
-            for kk, cmplx_freqs_kk in enumerate(cmplx_freqs[1:numFreqs]):
-                result = tb05ad(n, m, p, cmplx_freqs_kk, at,
-                                bt, ct, job='NH')
+            for kk, s_kk in enumerate(s_arr[1:len(s_arr)]):
+                result = tb05ad(n, m, p, s_kk, at, bt, ct, job='NH')
                 # When job='NH', result = (g_i, hinvb, info)
 
                 # kk+1 because enumerate starts at kk = 0.
                 # but zero-th spot is already filled.
-                Gfrf[:, :, kk+1] = result[0] + self.D
-
+                out[:, :, kk+1] = result[0] + self.D
+            
+            if not hasattr(s, '__len__'): 
+                # received a scalar s, squeeze down the array
+                out = np.squeeze(out, axis=2)
         except ImportError:  # Slycot unavailable. Fall back to horner.
-            for kk, cmplx_freqs_kk in enumerate(cmplx_freqs):
-                Gfrf[:, :, kk] = self.horner(cmplx_freqs_kk)
+            out = self.horner(s)
+        if squeeze and self.issiso():
+            return out[0][0]
+        else:
+            return out
 
-        #      mag           phase           omega
-        return np.abs(Gfrf), np.angle(Gfrf), omega
+    def horner(self, s):
+        """Evaluate systems's transfer function at complex frequencies s.
+        """
+        s_arr = np.array(s, ndmin=1) # force to be an array
+        # Preallocate 
+        out = empty((self.outputs, self.inputs, len(s_arr)), dtype=complex)
+        
+        for idx, s_idx in enumerate(s_arr):
+            out[:,:,idx] = \
+                np.dot(self.C, 
+                       solve(s_idx * eye(self.states) - self.A, self.B)) \
+                + self.D
+        if not hasattr(s, '__len__'): 
+            # received a scalar s, squeeze down the array along last dim
+            out = np.squeeze(out, axis=2)
+        return out
+
+    def freqresp(self, omega):
+        warn("StateSpace.freqresp(omega) will be deprecated in a "
+             "future release of python-control; use "
+             "StateSpace.frequency_response(sys, omega), or "
+             "freqresp(sys, omega) in the MATLAB compatibility module "
+             "instead", PendingDeprecationWarning)        
+        return self.frequency_response(omega)
 
     # Compute poles and zeros
     def pole(self):
diff --git a/control/tests/frd_test.py b/control/tests/frd_test.py
index fcbc10263..3d57753bc 100644
--- a/control/tests/frd_test.py
+++ b/control/tests/frd_test.py
@@ -37,7 +37,7 @@ def testSISOtf(self):
         assert isinstance(frd, FRD)
 
         np.testing.assert_array_almost_equal(
-            frd.freqresp([1.0]), h.freqresp([1.0]))
+            frd.frequency_response([1.0]), h.frequency_response([1.0]))
 
     def testOperators(self):
         # get two SISO transfer functions
@@ -48,39 +48,39 @@ def testOperators(self):
         f2 = FRD(h2, omega)
 
         np.testing.assert_array_almost_equal(
-            (f1 + f2).freqresp([0.1, 1.0, 10])[0],
-            (h1 + h2).freqresp([0.1, 1.0, 10])[0])
+            (f1 + f2).frequency_response([0.1, 1.0, 10])[0],
+            (h1 + h2).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1 + f2).freqresp([0.1, 1.0, 10])[1],
-            (h1 + h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 + f2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 + h2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1 - f2).freqresp([0.1, 1.0, 10])[0],
-            (h1 - h2).freqresp([0.1, 1.0, 10])[0])
+            (f1 - f2).frequency_response([0.1, 1.0, 10])[0],
+            (h1 - h2).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1 - f2).freqresp([0.1, 1.0, 10])[1],
-            (h1 - h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 - f2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 - h2).frequency_response([0.1, 1.0, 10])[1])
 
         # multiplication and division
         np.testing.assert_array_almost_equal(
-            (f1 * f2).freqresp([0.1, 1.0, 10])[1],
-            (h1 * h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 * f2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 * h2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1 / f2).freqresp([0.1, 1.0, 10])[1],
-            (h1 / h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 / f2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 / h2).frequency_response([0.1, 1.0, 10])[1])
 
         # with default conversion from scalar
         np.testing.assert_array_almost_equal(
-            (f1 * 1.5).freqresp([0.1, 1.0, 10])[1],
-            (h1 * 1.5).freqresp([0.1, 1.0, 10])[1])
+            (f1 * 1.5).frequency_response([0.1, 1.0, 10])[1],
+            (h1 * 1.5).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1 / 1.7).freqresp([0.1, 1.0, 10])[1],
-            (h1 / 1.7).freqresp([0.1, 1.0, 10])[1])
+            (f1 / 1.7).frequency_response([0.1, 1.0, 10])[1],
+            (h1 / 1.7).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (2.2 * f2).freqresp([0.1, 1.0, 10])[1],
-            (2.2 * h2).freqresp([0.1, 1.0, 10])[1])
+            (2.2 * f2).frequency_response([0.1, 1.0, 10])[1],
+            (2.2 * h2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (1.3 / f2).freqresp([0.1, 1.0, 10])[1],
-            (1.3 / h2).freqresp([0.1, 1.0, 10])[1])
+            (1.3 / f2).frequency_response([0.1, 1.0, 10])[1],
+            (1.3 / h2).frequency_response([0.1, 1.0, 10])[1])
 
     def testOperatorsTf(self):
         # get two SISO transfer functions
@@ -92,24 +92,24 @@ def testOperatorsTf(self):
         f2  # reference to avoid pyflakes error
 
         np.testing.assert_array_almost_equal(
-            (f1 + h2).freqresp([0.1, 1.0, 10])[0],
-            (h1 + h2).freqresp([0.1, 1.0, 10])[0])
+            (f1 + h2).frequency_response([0.1, 1.0, 10])[0],
+            (h1 + h2).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1 + h2).freqresp([0.1, 1.0, 10])[1],
-            (h1 + h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 + h2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 + h2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1 - h2).freqresp([0.1, 1.0, 10])[0],
-            (h1 - h2).freqresp([0.1, 1.0, 10])[0])
+            (f1 - h2).frequency_response([0.1, 1.0, 10])[0],
+            (h1 - h2).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1 - h2).freqresp([0.1, 1.0, 10])[1],
-            (h1 - h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 - h2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 - h2).frequency_response([0.1, 1.0, 10])[1])
         # multiplication and division
         np.testing.assert_array_almost_equal(
-            (f1 * h2).freqresp([0.1, 1.0, 10])[1],
-            (h1 * h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 * h2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 * h2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1 / h2).freqresp([0.1, 1.0, 10])[1],
-            (h1 / h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 / h2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 / h2).frequency_response([0.1, 1.0, 10])[1])
         # the reverse does not work
 
     def testbdalg(self):
@@ -121,45 +121,45 @@ def testbdalg(self):
         f2 = FRD(h2, omega)
 
         np.testing.assert_array_almost_equal(
-            (bdalg.series(f1, f2)).freqresp([0.1, 1.0, 10])[0],
-            (bdalg.series(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+            (bdalg.series(f1, f2)).frequency_response([0.1, 1.0, 10])[0],
+            (bdalg.series(h1, h2)).frequency_response([0.1, 1.0, 10])[0])
 
         np.testing.assert_array_almost_equal(
-            (bdalg.parallel(f1, f2)).freqresp([0.1, 1.0, 10])[0],
-            (bdalg.parallel(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+            (bdalg.parallel(f1, f2)).frequency_response([0.1, 1.0, 10])[0],
+            (bdalg.parallel(h1, h2)).frequency_response([0.1, 1.0, 10])[0])
 
         np.testing.assert_array_almost_equal(
-            (bdalg.feedback(f1, f2)).freqresp([0.1, 1.0, 10])[0],
-            (bdalg.feedback(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+            (bdalg.feedback(f1, f2)).frequency_response([0.1, 1.0, 10])[0],
+            (bdalg.feedback(h1, h2)).frequency_response([0.1, 1.0, 10])[0])
 
         np.testing.assert_array_almost_equal(
-            (bdalg.negate(f1)).freqresp([0.1, 1.0, 10])[0],
-            (bdalg.negate(h1)).freqresp([0.1, 1.0, 10])[0])
+            (bdalg.negate(f1)).frequency_response([0.1, 1.0, 10])[0],
+            (bdalg.negate(h1)).frequency_response([0.1, 1.0, 10])[0])
 
 #       append() and connect() not implemented for FRD objects
 #        np.testing.assert_array_almost_equal(
-#            (bdalg.append(f1, f2)).freqresp([0.1, 1.0, 10])[0],
-#            (bdalg.append(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+#            (bdalg.append(f1, f2)).frequency_response([0.1, 1.0, 10])[0],
+#            (bdalg.append(h1, h2)).frequency_response([0.1, 1.0, 10])[0])
 #
 #        f3 = bdalg.append(f1, f2, f2)
 #        h3 = bdalg.append(h1, h2, h2)
 #        Q = np.mat([ [1, 2], [2, -1] ])
 #        np.testing.assert_array_almost_equal(
-#           (bdalg.connect(f3, Q, [2], [1])).freqresp([0.1, 1.0, 10])[0],
-#            (bdalg.connect(h3, Q, [2], [1])).freqresp([0.1, 1.0, 10])[0])
+#           (bdalg.connect(f3, Q, [2], [1])).frequency_response([0.1, 1.0, 10])[0],
+#            (bdalg.connect(h3, Q, [2], [1])).frequency_response([0.1, 1.0, 10])[0])
 
     def testFeedback(self):
         h1 = TransferFunction([1], [1, 2, 2])
         omega = np.logspace(-1, 2, 10)
         f1 = FRD(h1, omega)
         np.testing.assert_array_almost_equal(
-            f1.feedback(1).freqresp([0.1, 1.0, 10])[0],
-            h1.feedback(1).freqresp([0.1, 1.0, 10])[0])
+            f1.feedback(1).frequency_response([0.1, 1.0, 10])[0],
+            h1.feedback(1).frequency_response([0.1, 1.0, 10])[0])
 
         # Make sure default argument also works
         np.testing.assert_array_almost_equal(
-            f1.feedback().freqresp([0.1, 1.0, 10])[0],
-            h1.feedback().freqresp([0.1, 1.0, 10])[0])
+            f1.feedback().frequency_response([0.1, 1.0, 10])[0],
+            h1.feedback().frequency_response([0.1, 1.0, 10])[0])
 
     def testFeedback2(self):
         h2 = StateSpace([[-1.0, 0], [0, -2.0]], [[0.4], [0.1]],
@@ -192,11 +192,11 @@ def testMIMO(self):
         omega = np.logspace(-1, 2, 10)
         f1 = FRD(sys, omega)
         np.testing.assert_array_almost_equal(
-            sys.freqresp([0.1, 1.0, 10])[0],
-            f1.freqresp([0.1, 1.0, 10])[0])
+            sys.frequency_response([0.1, 1.0, 10])[0],
+            f1.frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            sys.freqresp([0.1, 1.0, 10])[1],
-            f1.freqresp([0.1, 1.0, 10])[1])
+            sys.frequency_response([0.1, 1.0, 10])[1],
+            f1.frequency_response([0.1, 1.0, 10])[1])
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
     def testMIMOfb(self):
@@ -208,11 +208,11 @@ def testMIMOfb(self):
         f1 = FRD(sys, omega).feedback([[0.1, 0.3], [0.0, 1.0]])
         f2 = FRD(sys.feedback([[0.1, 0.3], [0.0, 1.0]]), omega)
         np.testing.assert_array_almost_equal(
-            f1.freqresp([0.1, 1.0, 10])[0],
-            f2.freqresp([0.1, 1.0, 10])[0])
+            f1.frequency_response([0.1, 1.0, 10])[0],
+            f2.frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            f1.freqresp([0.1, 1.0, 10])[1],
-            f2.freqresp([0.1, 1.0, 10])[1])
+            f1.frequency_response([0.1, 1.0, 10])[1],
+            f2.frequency_response([0.1, 1.0, 10])[1])
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
     def testMIMOfb2(self):
@@ -224,11 +224,11 @@ def testMIMOfb2(self):
         f1 = FRD(sys, omega).feedback(K)
         f2 = FRD(sys.feedback(K), omega)
         np.testing.assert_array_almost_equal(
-            f1.freqresp([0.1, 1.0, 10])[0],
-            f2.freqresp([0.1, 1.0, 10])[0])
+            f1.frequency_response([0.1, 1.0, 10])[0],
+            f2.frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            f1.freqresp([0.1, 1.0, 10])[1],
-            f2.freqresp([0.1, 1.0, 10])[1])
+            f1.frequency_response([0.1, 1.0, 10])[1],
+            f2.frequency_response([0.1, 1.0, 10])[1])
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
     def testMIMOMult(self):
@@ -240,11 +240,11 @@ def testMIMOMult(self):
         f1 = FRD(sys, omega)
         f2 = FRD(sys, omega)
         np.testing.assert_array_almost_equal(
-            (f1*f2).freqresp([0.1, 1.0, 10])[0],
-            (sys*sys).freqresp([0.1, 1.0, 10])[0])
+            (f1*f2).frequency_response([0.1, 1.0, 10])[0],
+            (sys*sys).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1*f2).freqresp([0.1, 1.0, 10])[1],
-            (sys*sys).freqresp([0.1, 1.0, 10])[1])
+            (f1*f2).frequency_response([0.1, 1.0, 10])[1],
+            (sys*sys).frequency_response([0.1, 1.0, 10])[1])
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
     def testMIMOSmooth(self):
@@ -257,14 +257,14 @@ def testMIMOSmooth(self):
         f1 = FRD(sys, omega, smooth=True)
         f2 = FRD(sys2, omega, smooth=True)
         np.testing.assert_array_almost_equal(
-            (f1*f2).freqresp([0.1, 1.0, 10])[0],
-            (sys*sys2).freqresp([0.1, 1.0, 10])[0])
+            (f1*f2).frequency_response([0.1, 1.0, 10])[0],
+            (sys*sys2).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1*f2).freqresp([0.1, 1.0, 10])[1],
-            (sys*sys2).freqresp([0.1, 1.0, 10])[1])
+            (f1*f2).frequency_response([0.1, 1.0, 10])[1],
+            (sys*sys2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1*f2).freqresp([0.1, 1.0, 10])[2],
-            (sys*sys2).freqresp([0.1, 1.0, 10])[2])
+            (f1*f2).frequency_response([0.1, 1.0, 10])[2],
+            (sys*sys2).frequency_response([0.1, 1.0, 10])[2])
 
     def testAgainstOctave(self):
         # with data from octave:
@@ -277,8 +277,8 @@ def testAgainstOctave(self):
         omega = np.logspace(-1, 2, 10)
         f1 = FRD(sys, omega)
         np.testing.assert_array_almost_equal(
-            (f1.freqresp([1.0])[0] *
-             np.exp(1j*f1.freqresp([1.0])[1])).reshape(3, 2),
+            (f1.frequency_response([1.0])[0] *
+             np.exp(1j*f1.frequency_response([1.0])[1])).reshape(3, 2),
             np.matrix('0.4-0.2j 0; 0 0.1-0.2j; 0 0.3-0.1j'))
 
     def test_string_representation(self):
@@ -383,36 +383,13 @@ def test_operator_conversion(self):
     def test_eval(self):
         sys_tf = ct.tf([1], [1, 2, 1])
         frd_tf = FRD(sys_tf, np.logspace(-1, 1, 3))
-        np.testing.assert_almost_equal(sys_tf.evalfr(1), frd_tf.eval(1))
+        np.testing.assert_almost_equal(sys_tf(1j), frd_tf.eval(1))
+        np.testing.assert_almost_equal(sys_tf(1j), frd_tf(1j))
 
         # Should get an error if we evaluate at an unknown frequency
-        self.assertRaises(ValueError, frd_tf.eval, 2)
-
-    # This test only works in Python 3 due to a conflict with the same
-    # warning type in other test modules (frd_test.py).  See
-    # https://bugs.python.org/issue4180 for more details
-    @unittest.skipIf(pysys.version_info < (3, 0), "test requires Python 3+")
-    def test_evalfr_deprecated(self):
-        sys_tf = ct.tf([1], [1, 2, 1])
-        frd_tf = FRD(sys_tf, np.logspace(-1, 1, 3))
-
-        # Deprecated version of the call (should generate warning)
-        import warnings
-        with warnings.catch_warnings():
-            # Make warnings generate an exception
-            warnings.simplefilter('error')
-
-            # Make sure that we get a pending deprecation warning
-            self.assertRaises(PendingDeprecationWarning, frd_tf.evalfr, 1.)
-
-        # FRD.evalfr() is being deprecated
-        import warnings
-        with warnings.catch_warnings():
-            # Make warnings generate an exception
-            warnings.simplefilter('error')
-
-            # Make sure that we get a pending deprecation warning
-            self.assertRaises(PendingDeprecationWarning, frd_tf.evalfr, 1.)
+        self.assertRaises(Exception, frd_tf.eval(2))
+        # Should get an error if we use __call__ at real-valued frequency
+        self.assertRaises(ValueError, frd_tf(2))
 
     def test_repr_str(self):
         # repr printing
diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py
index 9d59a1972..991a4cacd 100644
--- a/control/tests/freqresp_test.py
+++ b/control/tests/freqresp_test.py
@@ -18,7 +18,7 @@
 from control.exception import slycot_check
 
 
-class TestFreqresp(unittest.TestCase):
+class TestFrequencyResponse(unittest.TestCase):
    def setUp(self):
       self.A = np.matrix('1,1;0,1')
       self.C = np.matrix('1,0')
@@ -30,7 +30,7 @@ def test_siso(self):
       sys = StateSpace(self.A,B,self.C,D)
 
       # test frequency response
-      frq=sys.freqresp(self.omega)
+      frq=sys.frequency_response(self.omega)
 
       # test bode plot
       bode(sys)
@@ -84,7 +84,7 @@ def test_superimpose(self):
       self.assertEqual(len(ax.get_lines()), 2)
 
    def test_doubleint(self):
-      # 30 May 2016, RMM: added to replicate typecast bug in freqresp.py
+      # 30 May 2016, RMM: added to replicate typecast bug in frequency_response.py
       A = np.matrix('0, 1; 0, 0');
       B = np.matrix('0; 1');
       C = np.matrix('1, 0');
@@ -99,7 +99,7 @@ def test_mimo(self):
       D = np.matrix('0,0')
       sysMIMO = ss(self.A,B,self.C,D)
 
-      frqMIMO = sysMIMO.freqresp(self.omega)
+      frqMIMO = sysMIMO.frequency_response(self.omega)
       tfMIMO = tf(sysMIMO)
 
       #bode(sysMIMO) # - should throw not implemented exception
@@ -167,7 +167,7 @@ def test_discrete(self):
          omega_ok = np.linspace(10e-4,0.99,100) * np.pi/sys.dt
 
          # Test frequency response
-         ret = sys.freqresp(omega_ok)
+         ret = sys.frequency_response(omega_ok)
 
          # Check for warning if frequency is out of range
          import warnings
@@ -179,7 +179,7 @@ def test_discrete(self):
 
             # Look for a warning about sampling above Nyquist frequency
             omega_bad = np.linspace(10e-4,1.1,10) * np.pi/sys.dt
-            ret = sys.freqresp(omega_bad)
+            ret = sys.frequency_response(omega_bad)
             print("len(w) =", len(w))
             self.assertEqual(len(w), 1)
             self.assertIn("above", str(w[-1].message))
diff --git a/control/tests/statesp_array_test.py b/control/tests/statesp_array_test.py
index a2d034075..02e5d69d6 100644
--- a/control/tests/statesp_array_test.py
+++ b/control/tests/statesp_array_test.py
@@ -180,7 +180,7 @@ def test_multiply_ss(self):
         np.testing.assert_array_almost_equal(sys.C, C)
         np.testing.assert_array_almost_equal(sys.D, D)
 
-    def test_evalfr(self):
+    def test_call(self):
         """Evaluate the frequency response at one frequency."""
 
         A = [[-2, 0.5], [0.5, -0.3]]
@@ -196,22 +196,7 @@ def test_evalfr(self):
 
         # Correct versions of the call
         np.testing.assert_almost_equal(evalfr(sys, 1j), resp)
-        np.testing.assert_almost_equal(sys._evalfr(1.), resp)
-
-        # Deprecated version of the call (should generate warning)
-        import warnings
-        with warnings.catch_warnings(record=True) as w:
-            # Set up warnings filter to only show warnings in control module
-            warnings.filterwarnings("ignore")
-            warnings.filterwarnings("always", module="control")
-
-            # Make sure that we get a pending deprecation warning
-            sys.evalfr(1.)
-            assert len(w) == 1
-            assert issubclass(w[-1].category, PendingDeprecationWarning)
-
-            # Leave the warnings filter like we found it
-            warnings.resetwarnings()
+        np.testing.assert_almost_equal(sys(1.j), resp)
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
     def test_freq_resp(self):
@@ -233,7 +218,7 @@ def test_freq_resp(self):
                       [-0.438157380501337, -1.40720969147217]]]
         true_omega = [0.1, 10.]
 
-        mag, phase, omega = sys.freqresp(true_omega)
+        mag, phase, omega = sys.frequency_response(true_omega)
 
         np.testing.assert_almost_equal(mag, true_mag)
         np.testing.assert_almost_equal(phase, true_phase)
@@ -518,7 +503,7 @@ def test_horner(self):
         self.sys322.horner(1.+1.j)
 
         # Make sure result agrees with frequency response
-        mag, phase, omega = self.sys322.freqresp([1])
+        mag, phase, omega = self.sys322.frequency_response([1])
         np.testing.assert_array_almost_equal(
             self.sys322.horner(1.j),
             mag[:,:,0] * np.exp(1.j * phase[:,:,0]))
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 96404d79f..77a6d3ec2 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -189,7 +189,7 @@ def test_multiply_ss(self):
         np.testing.assert_array_almost_equal(sys.C, C)
         np.testing.assert_array_almost_equal(sys.D, D)
 
-    def test_evalfr(self):
+    def test_call(self):
         """Evaluate the frequency response at one frequency."""
 
         A = [[-2, 0.5], [0.5, -0.3]]
@@ -205,7 +205,14 @@ def test_evalfr(self):
 
         # Correct versions of the call
         np.testing.assert_almost_equal(evalfr(sys, 1j), resp)
-        np.testing.assert_almost_equal(sys._evalfr(1.), resp)
+        np.testing.assert_almost_equal(sys(1.j), resp)
+
+    def test_freqresp_deprecated(self):
+        A = [[-2, 0.5], [0.5, -0.3]]
+        B = [[0.3, -1.3], [0.1, 0.]]
+        C = [[0., 0.1], [-0.3, -0.2]]
+        D = [[0., -0.8], [-0.3, 0.]]
+        sys = StateSpace(A, B, C, D)
 
         # Deprecated version of the call (should generate warning)
         import warnings
@@ -215,8 +222,7 @@ def test_evalfr(self):
             warnings.filterwarnings("always", module="control")
 
             # Make sure that we get a pending deprecation warning
-            sys.evalfr(1.)
-            assert len(w) == 1
+            sys.freqresp(1.)
             assert issubclass(w[-1].category, PendingDeprecationWarning)
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
@@ -239,12 +245,12 @@ def test_freq_resp(self):
                       [-0.438157380501337, -1.40720969147217]]]
         true_omega = [0.1, 10.]
 
-        mag, phase, omega = sys.freqresp(true_omega)
+        mag, phase, omega = sys.frequency_response(true_omega)
 
         np.testing.assert_almost_equal(mag, true_mag)
         np.testing.assert_almost_equal(phase, true_phase)
         np.testing.assert_equal(omega, true_omega)
-
+        
     @unittest.skipIf(not slycot_check(), "slycot not installed")
     def test_minreal(self):
         """Test a minreal model reduction."""
@@ -652,26 +658,5 @@ def test_copy_constructor(self):
         # Change the A matrix for the original system
         linsys.A[0, 0] = -3
         np.testing.assert_array_equal(cpysys.A, [[-1]]) # original value
-
-    def test_sample_system_prewarping(self): 
-        """test that prewarping works when converting from cont to discrete time system"""
-        A = np.array([
-            [ 0.00000000e+00,  1.00000000e+00,  0.00000000e+00, 0.00000000e+00],
-            [-3.81097561e+01, -1.12500000e+00,  0.00000000e+00, 0.00000000e+00],
-            [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00, 1.00000000e+00],
-            [ 0.00000000e+00,  0.00000000e+00, -1.66356135e+04, -1.34748470e+01]])
-        B = np.array([
-            [    0.        ], [   38.1097561 ],[    0.     ],[16635.61352143]])
-        C = np.array([[0.90909091, 0.        , 0.09090909, 0.       ],])
-        wwarp = 50
-        Ts = 0.025
-        plant = StateSpace(A,B,C,0)
-        plant_d_warped = plant.sample(Ts, 'bilinear', prewarp_frequency=wwarp)
-        np.testing.assert_array_almost_equal(
-            evalfr(plant, wwarp*1j), 
-            evalfr(plant_d_warped, np.exp(wwarp*1j*Ts)), 
-            decimal=4)
-
-
 if __name__ == "__main__":
     unittest.main()
diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py
index 02e6c2b37..4ee314df1 100644
--- a/control/tests/xferfcn_test.py
+++ b/control/tests/xferfcn_test.py
@@ -386,7 +386,7 @@ def test_slice(self):
         self.assertEqual((sys1.inputs, sys1.outputs), (2, 1))
         self.assertEqual(sys1.dt, 0.5)
 
-    def test_evalfr_siso(self):
+    def test_call_siso(self):
         """Evaluate the frequency response of a SISO system at one frequency."""
 
         sys = TransferFunction([1., 3., 5], [1., 6., 2., -1])
@@ -400,38 +400,22 @@ def test_evalfr_siso(self):
         # Test call version as well
         np.testing.assert_almost_equal(sys(1.j), -0.5 - 0.5j)
         np.testing.assert_almost_equal(
-            sys(32.j), 0.00281959302585077 - 0.030628473607392j)
+            sys(32j), 0.00281959302585077 - 0.030628473607392j)
 
         # Test internal version (with real argument)
         np.testing.assert_array_almost_equal(
-            sys._evalfr(1.), np.array([[-0.5 - 0.5j]]))
+            sys(1j), np.array([[-0.5 - 0.5j]]))
         np.testing.assert_array_almost_equal(
-            sys._evalfr(32.),
+            sys(32j),
             np.array([[0.00281959302585077 - 0.030628473607392j]]))
 
-    # This test only works in Python 3 due to a conflict with the same
-    # warning type in other test modules (frd_test.py).  See
-    # https://bugs.python.org/issue4180 for more details
     @unittest.skipIf(pysys.version_info < (3, 0), "test requires Python 3+")
-    def test_evalfr_deprecated(self):
-        sys = TransferFunction([1., 3., 5], [1., 6., 2., -1])
-
-        # Deprecated version of the call (should generate warning)
-        import warnings
-        with warnings.catch_warnings():
-            # Make warnings generate an exception
-            warnings.simplefilter('error')
-
-            # Make sure that we get a pending deprecation warning
-            self.assertRaises(PendingDeprecationWarning, sys.evalfr, 1.)
-
-    @unittest.skipIf(pysys.version_info < (3, 0), "test requires Python 3+")
-    def test_evalfr_dtime(self):
+    def test_call_dtime(self):
         sys = TransferFunction([1., 3., 5], [1., 6., 2., -1], 0.1)
         np.testing.assert_array_almost_equal(sys(1j), -0.5 - 0.5j)
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
-    def test_evalfr_mimo(self):
+    def test_call_mimo(self):
         """Evaluate the frequency response of a MIMO system at one frequency."""
 
         num = [[[1., 2.], [0., 3.], [2., -1.]],
@@ -443,12 +427,12 @@ def test_evalfr_mimo(self):
                 [-0.083333333333333, -0.188235294117647 - 0.847058823529412j,
                  -1. - 8.j]]
 
-        np.testing.assert_array_almost_equal(sys._evalfr(2.), resp)
+        np.testing.assert_array_almost_equal(evalfr(sys, 2j), resp)
 
         # Test call version as well
         np.testing.assert_array_almost_equal(sys(2.j), resp)
 
-    def test_freqresp_siso(self):
+    def test_frequency_response_siso(self):
         """Evaluate the magnitude and phase of a SISO system at
         multiple frequencies."""
 
@@ -459,14 +443,14 @@ def test_freqresp_siso(self):
                        -1.32655885133871]]]
         trueomega = [0.1, 1., 10.]
 
-        mag, phase, omega = sys.freqresp(trueomega)
+        mag, phase, omega = sys.frequency_response(trueomega, squeeze=False)
 
         np.testing.assert_array_almost_equal(mag, truemag)
         np.testing.assert_array_almost_equal(phase, truephase)
         np.testing.assert_array_almost_equal(omega, trueomega)
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
-    def test_freqresp_mimo(self):
+    def test_frequency_response_mimo(self):
         """Evaluate the magnitude and phase of a MIMO system at
         multiple frequencies."""
 
@@ -489,7 +473,7 @@ def test_freqresp_mimo(self):
                        [-1.66852323, -1.89254688, -1.62050658],
                        [-0.13298964, -1.10714871, -2.75046720]]]
 
-        mag, phase, omega = sys.freqresp(true_omega)
+        mag, phase, omega = sys.frequency_response(true_omega)
 
         np.testing.assert_array_almost_equal(mag, true_mag)
         np.testing.assert_array_almost_equal(phase, true_phase)
diff --git a/control/xferfcn.py b/control/xferfcn.py
index f50d5141d..bb342c732 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -204,18 +204,21 @@ def __init__(self, *args):
 
         self._truncatecoeff()
 
-    def __call__(self, s):
-        """Evaluate the system's transfer function for a complex variable
+    def __call__(self, s, squeeze=True):
+        """Evaluate the system's transfer function at the complex frequency s/z.
 
-        For a SISO transfer function, returns the value of the
+        For a SISO transfer function, returns the complex value of the
         transfer function.  For a MIMO transfer fuction, returns a
-        matrix of values evaluated at complex variable s."""
+        matrix of values.
+        
+        If squeeze is True (default) and sys is single input, single output 
+        (SISO), return a 1D array or scalar depending on s.
 
-        if self.issiso():
-            # return a scalar
+        """    
+        if squeeze and self.issiso():
+            # return a scalar/1d array of outputs
             return self.horner(s)[0][0]
         else:
-            # return a matrix
             return self.horner(s)
 
     def _truncatecoeff(self):
@@ -608,40 +611,10 @@ def __getitem__(self, key):
         else:
             return TransferFunction(num, den, self.dt)
 
-    def evalfr(self, omega):
-        """Evaluate a transfer function at a single angular frequency.
-
-        self._evalfr(omega) returns the value of the transfer function
-        matrix with input value s = i * omega.
-
-        """
-        warn("TransferFunction.evalfr(omega) will be deprecated in a "
-             "future release of python-control; use evalfr(sys, omega*1j) "
-             "instead", PendingDeprecationWarning)
-        return self._evalfr(omega)
-
-    def _evalfr(self, omega):
-        """Evaluate a transfer function at a single angular frequency."""
-        # TODO: implement for discrete time systems
-        if isdtime(self, strict=True):
-            # Convert the frequency to discrete time
-            dt = timebase(self)
-            s = exp(1.j * omega * dt)
-            if np.any(omega * dt > pi):
-                warn("_evalfr: frequency evaluation above Nyquist frequency")
-        else:
-            s = 1.j * omega
-
-        return self.horner(s)
-
     def horner(self, s):
-        """Evaluate the systems's transfer function for a complex variable
-
-        Returns a matrix of values evaluated at complex variable s.
-        """
-
+        "Evaluate system's transfer function at complex frequencies s."
         # Preallocate the output.
-        if getattr(s, '__iter__', False):
+        if hasattr(s, '__len__'):
             out = empty((self.outputs, self.inputs, len(s)), dtype=complex)
         else:
             out = empty((self.outputs, self.inputs), dtype=complex)
@@ -654,59 +627,12 @@ def horner(self, s):
         return out
 
     def freqresp(self, omega):
-        """Evaluate the transfer function at a list of angular frequencies.
-
-        Reports the frequency response of the system,
-
-             G(j*omega) = mag*exp(j*phase)
-
-        for continuous time. For discrete time systems, the response is
-        evaluated around the unit circle such that
-
-             G(exp(j*omega*dt)) = mag*exp(j*phase).
-
-        Parameters
-        ----------
-        omega : array_like
-            A list of frequencies in radians/sec at which the system should be
-            evaluated. The list can be either a python list or a numpy array
-            and will be sorted before evaluation.
-
-        Returns
-        -------
-        mag : (self.outputs, self.inputs, len(omega)) ndarray
-            The magnitude (absolute value, not dB or log10) of the system
-            frequency response.
-        phase : (self.outputs, self.inputs, len(omega)) ndarray
-            The wrapped phase in radians of the system frequency response.
-        omega : ndarray or list or tuple
-            The list of sorted frequencies at which the response was
-            evaluated.
-        """
-        # Preallocate outputs.
-        numfreq = len(omega)
-        mag = empty((self.outputs, self.inputs, numfreq))
-        phase = empty((self.outputs, self.inputs, numfreq))
-
-        # Figure out the frequencies
-        omega.sort()
-        if isdtime(self, strict=True):
-            dt = timebase(self)
-            slist = np.array([exp(1.j * w * dt) for w in omega])
-            if max(omega) * dt > pi:
-                warn("freqresp: frequency evaluation above Nyquist frequency")
-        else:
-            slist = np.array([1j * w for w in omega])
-
-        # Compute frequency response for each input/output pair
-        for i in range(self.outputs):
-            for j in range(self.inputs):
-                fresp = (polyval(self.num[i][j], slist) /
-                         polyval(self.den[i][j], slist))
-                mag[i, j, :] = abs(fresp)
-                phase[i, j, :] = angle(fresp)
-
-        return mag, phase, omega
+        warn("TransferFunction.freqresp(omega) will be deprecated in a "
+             "future release of python-control; use "
+             "TransferFunction.frequency_response(sys, omega), or "
+             "freqresp(sys, omega) in the MATLAB compatibility module "
+             "instead", PendingDeprecationWarning)        
+        return self.frequency_response(omega)
 
     def pole(self):
         """Compute the poles of a transfer function."""
diff --git a/doc/control.rst b/doc/control.rst
index 57d64b1eb..0705e1b9f 100644
--- a/doc/control.rst
+++ b/doc/control.rst
@@ -172,6 +172,7 @@ Utility functions and conversions
     tfdata
     timebase
     timebaseEqual
+    common_timebase
     unwrap
     use_fbs_defaults
     use_matlab_defaults
diff --git a/examples/robust_mimo.py b/examples/robust_mimo.py
index d4e1335e6..9cc5a1c3b 100644
--- a/examples/robust_mimo.py
+++ b/examples/robust_mimo.py
@@ -43,7 +43,7 @@ def triv_sigma(g, w):
     g - LTI object, order n
     w - frequencies, length m
     s - (m,n) array of singular values of g(1j*w)"""
-    m, p, _ = g.freqresp(w)
+    m, p, _ = g.frequency_response(w)
     sjw = (m*np.exp(1j*p)).transpose(2, 0, 1)
     sv = np.linalg.svd(sjw, compute_uv=False)
     return sv
diff --git a/examples/robust_siso.py b/examples/robust_siso.py
index 87fcdb707..17ce10927 100644
--- a/examples/robust_siso.py
+++ b/examples/robust_siso.py
@@ -50,10 +50,10 @@
 
 # frequency response
 omega = np.logspace(-2, 2, 101)
-ws1mag, _, _ = ws1.freqresp(omega)
-s1mag, _, _ = s1.freqresp(omega)
-ws2mag, _, _ = ws2.freqresp(omega)
-s2mag, _, _ = s2.freqresp(omega)
+ws1mag, _, _ = ws1.frequency_response(omega)
+s1mag, _, _ = s1.frequency_response(omega)
+ws2mag, _, _ = ws2.frequency_response(omega)
+s2mag, _, _ = s2.frequency_response(omega)
 
 plt.figure(1)
 # text uses log-scaled absolute, but dB are probably more familiar to most control engineers

From b05492c08e73ebb29739acdf6924df6a4fe73764 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Sun, 16 Aug 2020 00:51:38 -0700
Subject: [PATCH 02/32] fixed remaining failing unit tests and cleanup of
 docfiles

---
 control/frdata.py                   |  13 +--
 control/lti.py                      |   9 +-
 control/margins.py                  |  17 ++--
 control/rlocus.py                   |  10 +--
 control/statesp.py                  | 122 +++++++++++++++++-----------
 control/tests/frd_test.py           |   4 +-
 control/tests/margin_test.py        |   8 +-
 control/tests/statesp_array_test.py |   2 +-
 control/xferfcn.py                  |  63 ++++++++------
 9 files changed, 147 insertions(+), 101 deletions(-)

diff --git a/control/frdata.py b/control/frdata.py
index 084f4aba3..6016e1460 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -47,7 +47,7 @@
 # External function declarations
 from warnings import warn
 import numpy as np
-from numpy import angle, array, empty, ones, isin, \
+from numpy import angle, array, empty, ones, \
     real, imag, absolute, eye, linalg, where, dot, sort
 from scipy.interpolate import splprep, splev
 from .lti import LTI
@@ -354,7 +354,7 @@ def eval(self, omega, squeeze=True):
         
         squeeze: bool, optional (default=True)
             If True and sys is single input, single output (SISO), return a 
-            1D array or scalar depending on omega's length.
+            1D array or scalar the same length as omega.
 
         """
         omega_array = np.array(omega, ndmin=1) # array-like version of omega 
@@ -362,15 +362,16 @@ def eval(self, omega, squeeze=True):
             raise ValueError("FRD.eval can only accept real-valued omega")
 
         if self.ifunc is None:
-            elements = isin(self.omega, omega)
+            elements = np.isin(self.omega, omega) # binary array
             if sum(elements) < len(omega_array):
                 raise ValueError(
-                    "not all frequencies omega are in frequency list of FRD "
-                    "system. Try an interpolating FRD for additional points.")
+                    '''not all frequencies omega are in frequency list of FRD 
+                    system. Try an interpolating FRD for additional points.''')
             else:
                 out = self.fresp[:, :, elements]
         else:
-            out = empty((self.outputs, self.inputs, len(omega_array)))
+            out = empty((self.outputs, self.inputs, len(omega_array)), 
+                         dtype=complex)
             for i in range(self.outputs):
                 for j in range(self.inputs):
                     for k, w in enumerate(omega_array):
diff --git a/control/lti.py b/control/lti.py
index 6d39955bb..c6bf318a9 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -465,9 +465,14 @@ def evalfr(sys, x, squeeze=True):
 
     .. todo:: Add example with MIMO system
     """
+    out = sys.horner(x)
+    if not hasattr(x, '__len__'): 
+            # received a scalar x, squeeze down the array along last dim
+            out = np.squeeze(out, axis=2)
     if squeeze and issiso(sys):
-        return sys.horner(x)[0][0]
-    return sys.horner(x)
+        return out[0][0]
+    else: 
+        return out
 
 
 def freqresp(sys, omega, squeeze=True):
diff --git a/control/margins.py b/control/margins.py
index 4094b43bf..cda5095aa 100644
--- a/control/margins.py
+++ b/control/margins.py
@@ -213,15 +213,15 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
         # a bit coarse, have the interpolated frd evaluated again
         def _mod(w):
             """Calculate |G(jw)| - 1"""
-            return np.abs(sys(1j * w)[0][0]) - 1
+            return np.abs(sys(1j * w)) - 1
 
         def _arg(w):
             """Calculate the phase angle at -180 deg"""
-            return np.angle(-sys(1j * w)[0][0])
+            return np.angle(-sys(1j * w))
 
         def _dstab(w):
             """Calculate the distance from -1 point"""
-            return np.abs(sys(1j * w)[0][0] + 1.)
+            return np.abs(sys(1j * w) + 1.)
 
         # Find all crossings, note that this depends on omega having
         # a correct range
@@ -232,7 +232,7 @@ def _dstab(w):
 
         # find the phase crossings ang(H(jw) == -180
         widx = np.where(np.diff(np.sign(_arg(sys.omega))))[0]
-        widx = widx[np.real(sys(1j * sys.omega[widx])[0][0]) <= 0]
+        widx = widx[np.real(sys(1j * sys.omega[widx])) <= 0]
         w_180 = np.array(
             [sp.optimize.brentq(_arg, sys.omega[i], sys.omega[i+1])
              for i in widx])
@@ -296,11 +296,10 @@ def phase_crossover_frequencies(sys):
     (array([ 1.73205081,  0.        ]), array([-0.5 ,  0.25]))
     """
 
+    if not sys.issiso(): 
+        raise ValueError("MIMO systems not yet implemented.")
     # Convert to a transfer function
     tf = xferfcn._convert_to_transfer_function(sys)
-
-    # if not siso, fall back to (0,0) element
-    #! TODO: should add a check and warning here
     num = tf.num[0][0]
     den = tf.den[0][0]
 
@@ -313,7 +312,9 @@ def phase_crossover_frequencies(sys):
 
     # using real() to avoid rounding errors and results like 1+0j
     # it would be nice to have a vectorized version of self.evalfr here
-    gain = np.real(np.asarray([tf(1j * f)[0][0] for f in realposfreq]))
+    # update Sawyer B. Fuller 2020.08.15: your wish is my command.
+    #gain = np.real(np.asarray([tf(1j * f) for f in realposfreq]))
+    gain = np.real(tf(1j * realposfreq))
 
     return realposfreq, gain
 
diff --git a/control/rlocus.py b/control/rlocus.py
index 9f7ff4568..39c2fcfb6 100644
--- a/control/rlocus.py
+++ b/control/rlocus.py
@@ -493,7 +493,7 @@ def _RLFindRoots(nump, denp, kvect):
     """Find the roots for the root locus."""
     # Convert numerator and denominator to polynomials if they aren't
     roots = []
-    for k in kvect:
+    for k in np.array(kvect, ndmin=1):
         curpoly = denp + k * nump
         curroots = curpoly.r
         if len(curroots) < denp.order:
@@ -577,10 +577,10 @@ def _RLFeedbackClicksPoint(event, sys, fig, ax_rlocus, sisotool=False):
     # 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.horner(s)
-        K_xlim = -1. / sys.horner(
+        K = -1. / sys(s)
+        K_xlim = -1. / sys(
             complex(event.xdata + 0.05 * abs(xlim[1] - xlim[0]), event.ydata))
-        K_ylim = -1. / sys.horner(
+        K_ylim = -1. / sys(
             complex(event.xdata, event.ydata + 0.05 * abs(ylim[1] - ylim[0])))
 
     except TypeError:
@@ -621,7 +621,7 @@ def _RLFeedbackClicksPoint(event, sys, fig, ax_rlocus, sisotool=False):
             ax_rlocus.plot(s.real, s.imag, 'k.', marker='s', markersize=8,
                            zorder=20, label='gain_point')
 
-        return K.real[0][0]
+        return K.real
 
 
 def _removeLine(label, ax):
diff --git a/control/statesp.py b/control/statesp.py
index a4f3aec71..40748c9fa 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -437,60 +437,91 @@ def __rdiv__(self, other):
 
         raise NotImplementedError("StateSpace.__rdiv__ is not implemented yet.")
 
-    def __call__(self, s, squeeze=True):
-        """Evaluate system's transfer function at complex frequencies s/z.
+    def __call__(self, x, squeeze=True):
+        """Evaluate system's transfer function at complex frequencies.
+        
+        Evaluates at complex frequency x, where x is s or z dependong on 
+        whether the system is continuous or discrete-time.
         
         If squeeze is True (default) and sys is single input, single output 
-        (SISO), return a 1D array or scalar depending on s.
+        (SISO), return a 1D array or scalar depending on the size of x. 
+        For a MIMO system, returns an (n_outputs, n_inputs, n_x) array.
 
         """        
-        # Use TB05AD from Slycot if available
+        # Use Slycot if available
         try:
-            from slycot import tb05ad
-
-            # preallocate
-            s_arr = np.array(s, ndmin=1) # array-like version of s
-            out = np.empty((self.outputs, self.inputs, len(s_arr)), 
-                           dtype=complex)
-            n = self.states
-            m = self.inputs
-            p = self.outputs
-            # The first call both evaluates C(sI-A)^-1 B and also returns
-            # Hessenberg transformed matrices at, bt, ct.
-            result = tb05ad(n, m, p, s_arr[0], self.A,
-                            self.B, self.C, job='NG')
-            # When job='NG', result = (at, bt, ct, g_i, hinvb, info)
-            at = result[0]
-            bt = result[1]
-            ct = result[2]
-
-            # TB05AD frequency evaluation does not include direct feedthrough.
-            out[:, :, 0] = result[3] + self.D
-
-            # Now, iterate through the remaining frequencies using the
-            # transformed state matrices, at, bt, ct.
-
-            # Start at the second frequency, already have the first.
-            for kk, s_kk in enumerate(s_arr[1:len(s_arr)]):
-                result = tb05ad(n, m, p, s_kk, at, bt, ct, job='NH')
-                # When job='NH', result = (g_i, hinvb, info)
-
-                # kk+1 because enumerate starts at kk = 0.
-                # but zero-th spot is already filled.
-                out[:, :, kk+1] = result[0] + self.D
-            
-            if not hasattr(s, '__len__'): 
-                # received a scalar s, squeeze down the array
-                out = np.squeeze(out, axis=2)
-        except ImportError:  # Slycot unavailable. Fall back to horner.
-            out = self.horner(s)
+            out = self.slycot_horner(x)
+        except ImportError:  # Slycot unavailable. use built-in horner.
+            out = self.horner(x)
+        if not hasattr(x, '__len__'): 
+            # received a scalar x, squeeze down the array along last dim
+            out = np.squeeze(out, axis=2)
         if squeeze and self.issiso():
             return out[0][0]
         else:
             return out
 
+    def slycot_horner(self, s):
+        """Evaluate system's transfer function at complex frequencies s 
+        using Horner's method from Slycot.
+
+        Expects inputs and outputs to be formatted correctly. Use __call__
+        for a more user-friendly interface. 
+
+        Parameters
+            s : array-like
+
+        Returns
+            output : array of size (outputs, inputs, len(s))
+        
+        """
+        from slycot import tb05ad
+
+        # preallocate
+        s_arr = np.array(s, ndmin=1) # array-like version of s
+        out = np.empty((self.outputs, self.inputs, len(s_arr)), 
+                        dtype=complex)
+        n = self.states
+        m = self.inputs
+        p = self.outputs
+        # The first call both evaluates C(sI-A)^-1 B and also returns
+        # Hessenberg transformed matrices at, bt, ct.
+        result = tb05ad(n, m, p, s_arr[0], self.A,
+                        self.B, self.C, job='NG')
+        # When job='NG', result = (at, bt, ct, g_i, hinvb, info)
+        at = result[0]
+        bt = result[1]
+        ct = result[2]
+
+        # TB05AD frequency evaluation does not include direct feedthrough.
+        out[:, :, 0] = result[3] + self.D
+
+        # Now, iterate through the remaining frequencies using the
+        # transformed state matrices, at, bt, ct.
+
+        # Start at the second frequency, already have the first.
+        for kk, s_kk in enumerate(s_arr[1:len(s_arr)]):
+            result = tb05ad(n, m, p, s_kk, at, bt, ct, job='NH')
+            # When job='NH', result = (g_i, hinvb, info)
+
+            # kk+1 because enumerate starts at kk = 0.
+            # but zero-th spot is already filled.
+            out[:, :, kk+1] = result[0] + self.D
+        return out
+
     def horner(self, s):
-        """Evaluate systems's transfer function at complex frequencies s.
+        """Evaluate system's transfer function at complex frequencies s 
+        using Horner's method.
+
+        Expects inputs and outputs to be formatted correctly. Use __call__
+        for a more user-friendly interface. 
+
+        Parameters
+            s : array-like
+
+        Returns
+            output : array of size (outputs, inputs, len(s))
+        
         """
         s_arr = np.array(s, ndmin=1) # force to be an array
         # Preallocate 
@@ -501,9 +532,6 @@ def horner(self, s):
                 np.dot(self.C, 
                        solve(s_idx * eye(self.states) - self.A, self.B)) \
                 + self.D
-        if not hasattr(s, '__len__'): 
-            # received a scalar s, squeeze down the array along last dim
-            out = np.squeeze(out, axis=2)
         return out
 
     def freqresp(self, omega):
@@ -879,7 +907,7 @@ def dcgain(self):
             if self.isctime():
                 gain = np.asarray(self.D-self.C.dot(np.linalg.solve(self.A, self.B)))
             else:
-                gain = self.horner(1)
+                gain = np.squeeze(self.horner(1))
         except LinAlgError:
             # eigenvalue at DC
             gain = np.tile(np.nan, (self.outputs, self.inputs))
diff --git a/control/tests/frd_test.py b/control/tests/frd_test.py
index 3d57753bc..5335a473e 100644
--- a/control/tests/frd_test.py
+++ b/control/tests/frd_test.py
@@ -387,9 +387,9 @@ def test_eval(self):
         np.testing.assert_almost_equal(sys_tf(1j), frd_tf(1j))
 
         # Should get an error if we evaluate at an unknown frequency
-        self.assertRaises(Exception, frd_tf.eval(2))
+        self.assertRaises(ValueError, frd_tf.eval, 2)
         # Should get an error if we use __call__ at real-valued frequency
-        self.assertRaises(ValueError, frd_tf(2))
+        self.assertRaises(ValueError, frd_tf, 2)
 
     def test_repr_str(self):
         # repr printing
diff --git a/control/tests/margin_test.py b/control/tests/margin_test.py
index 2f60c7bc6..70d6b9b5b 100755
--- a/control/tests/margin_test.py
+++ b/control/tests/margin_test.py
@@ -76,7 +76,7 @@ def test_stability_margins_omega(sys, refout, refoutall):
 def test_stability_margins_3input(sys, refout, refoutall):
     """Test stability_margins() function with mag, phase, omega input"""
     omega = np.logspace(-2, 2, 2000)
-    mag, phase, omega_ = sys.freqresp(omega)
+    mag, phase, omega_ = sys.frequency_response(omega)
     out = stability_margins((mag, phase*180/np.pi, omega_))
     assert_allclose(out, refout, atol=1.5e-3)
 
@@ -92,7 +92,7 @@ def test_margin_sys(sys, refout, refoutall):
 def test_margin_3input(sys, refout, refoutall):
     """Test margin() function with mag, phase, omega input"""
     omega = np.logspace(-2, 2, 2000)
-    mag, phase, omega_ = sys.freqresp(omega)
+    mag, phase, omega_ = sys.frequency_response(omega)
     out = margin((mag, phase*180/np.pi, omega_))
     assert_allclose(out, np.array(refout)[[0, 1, 3, 4]], atol=1.5e-3)
 
@@ -113,7 +113,7 @@ def test_phase_crossover_frequencies():
                            [[3], [4]]],
                           [[[1, 2, 3, 4], [1, 1]],
                            [[1, 1], [1, 1]]])
-    omega, gain = phase_crossover_frequencies(tf)
+    omega, gain = phase_crossover_frequencies(tf[0,0])
     assert_allclose(omega, [1.73205,  0.], atol=1.5e-3)
     assert_allclose(gain, [-0.5,  0.25], atol=1.5e-3)
 
@@ -123,7 +123,7 @@ def test_mag_phase_omega():
     sys = TransferFunction(15, [1, 6, 11, 6])
     out = stability_margins(sys)
     omega = np.logspace(-2, 2, 1000)
-    mag, phase, omega = sys.freqresp(omega)
+    mag, phase, omega = sys.frequency_response(omega)
     out2 = stability_margins((mag, phase*180/np.pi, omega))
     ind = [0, 1, 3, 4]   # indices of gm, pm, wg, wp -- ignore sm
     marg1 = np.array(out)[ind]
diff --git a/control/tests/statesp_array_test.py b/control/tests/statesp_array_test.py
index 02e5d69d6..ffc89ab97 100644
--- a/control/tests/statesp_array_test.py
+++ b/control/tests/statesp_array_test.py
@@ -505,7 +505,7 @@ def test_horner(self):
         # Make sure result agrees with frequency response
         mag, phase, omega = self.sys322.frequency_response([1])
         np.testing.assert_array_almost_equal(
-            self.sys322.horner(1.j),
+            np.squeeze(self.sys322.horner(1.j)),
             mag[:,:,0] * np.exp(1.j * phase[:,:,0]))
 
     def tearDown(self):
diff --git a/control/xferfcn.py b/control/xferfcn.py
index bb342c732..f43f282b4 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -204,22 +204,48 @@ def __init__(self, *args):
 
         self._truncatecoeff()
 
-    def __call__(self, s, squeeze=True):
-        """Evaluate the system's transfer function at the complex frequency s/z.
-
-        For a SISO transfer function, returns the complex value of the
-        transfer function.  For a MIMO transfer fuction, returns a
-        matrix of values.
+    def __call__(self, x, squeeze=True):
+        """Evaluate system's transfer function at complex frequencies.
+        
+        Evaluates at complex frequencies x, where x is s or z dependong on 
+        whether the system is continuous or discrete-time.
         
         If squeeze is True (default) and sys is single input, single output 
-        (SISO), return a 1D array or scalar depending on s.
-
-        """    
+        (SISO), return a 1D array or scalar depending on the shape of x. 
+        For a MIMO system, returns an (n_outputs, n_inputs, n_x) array.
+
+        """        
+        out = self.horner(x)
+        if not hasattr(x, '__len__'): 
+            # received a scalar x, squeeze down the array along last dim
+            out = np.squeeze(out, axis=2)
         if squeeze and self.issiso():
             # return a scalar/1d array of outputs
-            return self.horner(s)[0][0]
+            return out[0][0]
         else:
-            return self.horner(s)
+            return out
+
+    def horner(self, s):
+        """Evaluate system's transfer function at complex frequencies s 
+        using Horner's method.
+
+        Expects inputs and outputs to be formatted correctly. Use __call__
+        for a more user-friendly interface. 
+
+        Parameters
+            s : array-like
+
+        Returns
+            output : array of size (outputs, inputs, len(s))
+        
+        """
+        s_arr = np.array(s, ndmin=1) # force to be an array
+        out = empty((self.outputs, self.inputs, len(s_arr)), dtype=complex)
+        for i in range(self.outputs):
+            for j in range(self.inputs):
+                out[i][j] = (polyval(self.num[i][j], s) /
+                             polyval(self.den[i][j], s))
+        return out
 
     def _truncatecoeff(self):
         """Remove extraneous zero coefficients from num and den.
@@ -611,21 +637,6 @@ def __getitem__(self, key):
         else:
             return TransferFunction(num, den, self.dt)
 
-    def horner(self, s):
-        "Evaluate system's transfer function at complex frequencies s."
-        # Preallocate the output.
-        if hasattr(s, '__len__'):
-            out = empty((self.outputs, self.inputs, len(s)), dtype=complex)
-        else:
-            out = empty((self.outputs, self.inputs), dtype=complex)
-
-        for i in range(self.outputs):
-            for j in range(self.inputs):
-                out[i][j] = (polyval(self.num[i][j], s) /
-                             polyval(self.den[i][j], s))
-
-        return out
-
     def freqresp(self, omega):
         warn("TransferFunction.freqresp(omega) will be deprecated in a "
              "future release of python-control; use "

From 2c4ac62e5b69d6a8511dc2173c32338d98b43e01 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Sun, 16 Aug 2020 00:59:08 -0700
Subject: [PATCH 03/32] minor fix following suggested change in frd str method

---
 control/frdata.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/control/frdata.py b/control/frdata.py
index 6016e1460..cd0680f8d 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -156,7 +156,6 @@ def __str__(self):
         mimo = self.inputs > 1 or self.outputs > 1
         outstr = ['Frequency response data']
 
-        #mt, pt, wt = self.frequency_response(self.omega)
         for i in range(self.inputs):
             for j in range(self.outputs):
                 if mimo:
@@ -164,9 +163,10 @@ def __str__(self):
                 outstr.append('Freq [rad/s]  Response')
                 outstr.append('------------  ---------------------')
                 outstr.extend(
-                    ['%12.3f  %10.4g%+10.4gj' % (w, m, p)
-                     for m, p, w in zip(real(self.fresp[j, i, :]),
-                                        imag(self.fresp[j, i, :]), self.omega)])
+                    ['%12.3f  %10.4g%+10.4gj' % (w, re, im)
+                     for w, re, im in zip(self.omega,
+                                          real(self.fresp[j, i, :]),
+                                          imag(self.fresp[j, i, :]))])
 
         return '\n'.join(outstr)
 

From 6213a543debba42c270d4cf7cad88776ade554b9 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Sun, 16 Aug 2020 01:11:03 -0700
Subject: [PATCH 04/32] fixed a few bugs that slipped through

---
 control/matlab/__init__.py  | 4 ++--
 control/statesp.py          | 1 +
 control/tests/iosys_test.py | 4 ++--
 3 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/control/matlab/__init__.py b/control/matlab/__init__.py
index 413dc6d86..a2ec5efc1 100644
--- a/control/matlab/__init__.py
+++ b/control/matlab/__init__.py
@@ -224,8 +224,8 @@
 \*  :func:`~control.nichols`    Nichols plot
 \*  :func:`margin`              gain and phase margins
 \   lti/allmargin               all crossover frequencies and margins
-\*  :func:`freqresp`            frequency response over a frequency grid
-\*  :func:`evalfr`              frequency response at single frequency
+\*  :func:`freqresp`            frequency response 
+\*  :func:`evalfr`              frequency response at complex frequency s
 ==  ==========================  ============================================
 
 
diff --git a/control/statesp.py b/control/statesp.py
index 40748c9fa..1e502323e 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -527,6 +527,7 @@ def horner(self, s):
         # Preallocate 
         out = empty((self.outputs, self.inputs, len(s_arr)), dtype=complex)
         
+        #TODO: can this be vectorized? 
         for idx, s_idx in enumerate(s_arr):
             out[:,:,idx] = \
                 np.dot(self.C, 
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index 0738e8b18..fdbcc6088 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -860,8 +860,8 @@ def test_lineariosys_statespace(self):
         np.testing.assert_array_equal(
             iosys_siso.pole(), self.siso_linsys.pole())
         omega = np.logspace(.1, 10, 100)
-        mag_io, phase_io, omega_io = iosys_siso.freqresp(omega)
-        mag_ss, phase_ss, omega_ss = self.siso_linsys.freqresp(omega)
+        mag_io, phase_io, omega_io = iosys_siso.frequency_response(omega)
+        mag_ss, phase_ss, omega_ss = self.siso_linsys.frequency_response(omega)
         np.testing.assert_array_equal(mag_io, mag_ss)
         np.testing.assert_array_equal(phase_io, phase_ss)
         np.testing.assert_array_equal(omega_io, omega_ss)

From 8486b7a46cb4d518614d0845b3c6f76f8bb20887 Mon Sep 17 00:00:00 2001
From: sawyerbfuller <58706249+sawyerbfuller@users.noreply.github.com>
Date: Mon, 17 Aug 2020 10:49:14 -0700
Subject: [PATCH 05/32] Update control/frdata.py

Co-authored-by: Ben Greiner <code@bnavigator.de>
---
 control/frdata.py | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/control/frdata.py b/control/frdata.py
index cd0680f8d..686ae8823 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -391,12 +391,24 @@ def __call__(self, s, squeeze=True):
         transfer function.  For a MIMO transfer fuction, returns a
         matrix of values.
         
-        Raises an error if s is not purely imaginary. 
-        
+        Parameters
+        ----------
+        s : scalar or array_like
+            Complex frequencies
         squeeze: bool, optional (default=True)
             If True and sys is single input, single output (SISO), return a 
             1D array or scalar depending on omega's length.
 
+        Returns
+        -------
+        ndarray or scalar
+            Frequency response
+
+        Raises
+        ------
+        ValueError
+            If `s` is not purely imaginary. 
+
         """
         if any(abs(np.array(s, ndmin=1).real) > 0):
             raise ValueError("__call__: FRD systems can only accept"

From c8bf92f15729a973d011566b08b8159ad95f3841 Mon Sep 17 00:00:00 2001
From: sawyerbfuller <58706249+sawyerbfuller@users.noreply.github.com>
Date: Mon, 17 Aug 2020 10:49:51 -0700
Subject: [PATCH 06/32] Update control/frdata.py

Co-authored-by: Ben Greiner <code@bnavigator.de>
---
 control/frdata.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/control/frdata.py b/control/frdata.py
index 686ae8823..1fd359d67 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -352,10 +352,11 @@ def eval(self, omega, squeeze=True):
         entry in the omega vector. An interpolating FRD can return
         intermediate values.     
         
+        Parameters
+        ----------
         squeeze: bool, optional (default=True)
             If True and sys is single input, single output (SISO), return a 
             1D array or scalar the same length as omega.
-
         """
         omega_array = np.array(omega, ndmin=1) # array-like version of omega 
         if any(omega_array.imag > 0): 

From 371986c50d1319af713095d19a509faeb2115a5b Mon Sep 17 00:00:00 2001
From: sawyerbfuller <58706249+sawyerbfuller@users.noreply.github.com>
Date: Mon, 17 Aug 2020 10:50:43 -0700
Subject: [PATCH 07/32] Update control/lti.py

doc fix

Co-authored-by: Ben Greiner <code@bnavigator.de>
---
 control/lti.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/control/lti.py b/control/lti.py
index c6bf318a9..e89e81825 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -436,7 +436,7 @@ def evalfr(sys, x, squeeze=True):
     ----------
     sys: StateSpace or TransferFunction
         Linear system
-    x: scalar or array-like
+    x: scalar or array_like
         Complex number
     squeeze: bool, optional (default=True)
         If True and sys is single input, single output (SISO), return a 

From d8fbaa0c42b26aba42376060926768b4d0b9c892 Mon Sep 17 00:00:00 2001
From: sawyerbfuller <58706249+sawyerbfuller@users.noreply.github.com>
Date: Mon, 17 Aug 2020 11:04:50 -0700
Subject: [PATCH 08/32] Update control/statesp.py

suggested doctoring fix for statesspace.slycot_horner

Co-authored-by: Ben Greiner <code@bnavigator.de>
---
 control/statesp.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/control/statesp.py b/control/statesp.py
index 1e502323e..debe58d70 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -469,7 +469,9 @@ def slycot_horner(self, s):
         for a more user-friendly interface. 
 
         Parameters
-            s : array-like
+        ----------
+        s : array_like
+            Complex frequency
 
         Returns
             output : array of size (outputs, inputs, len(s))

From 17bc2a818da112881e49e9daa48a8b16c08ebace Mon Sep 17 00:00:00 2001
From: sawyerbfuller <58706249+sawyerbfuller@users.noreply.github.com>
Date: Mon, 17 Aug 2020 11:08:13 -0700
Subject: [PATCH 09/32] Apply suggestions from code review

doctoring fixes to adhere to numpydoc convention

Co-authored-by: Ben Greiner <code@bnavigator.de>
---
 control/statesp.py | 14 +++++++++-----
 control/xferfcn.py |  9 ++++++---
 2 files changed, 15 insertions(+), 8 deletions(-)

diff --git a/control/statesp.py b/control/statesp.py
index debe58d70..08eab1beb 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -474,8 +474,9 @@ def slycot_horner(self, s):
             Complex frequency
 
         Returns
-            output : array of size (outputs, inputs, len(s))
-        
+        -------
+        output : (outputs, inputs, len(s)) complex ndarray
+            Frequency response
         """
         from slycot import tb05ad
 
@@ -519,11 +520,14 @@ def horner(self, s):
         for a more user-friendly interface. 
 
         Parameters
-            s : array-like
+        ----------
+        s : array_like
+            Complex frequencies
 
         Returns
-            output : array of size (outputs, inputs, len(s))
-        
+        -------
+        output : (outputs, inputs, len(s)) complex ndarray
+            Frequency response
         """
         s_arr = np.array(s, ndmin=1) # force to be an array
         # Preallocate 
diff --git a/control/xferfcn.py b/control/xferfcn.py
index f43f282b4..c8de0b265 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -233,11 +233,14 @@ def horner(self, s):
         for a more user-friendly interface. 
 
         Parameters
-            s : array-like
+        ----------
+        s : array_like
+            Complex frequencies
 
         Returns
-            output : array of size (outputs, inputs, len(s))
-        
+        -------
+        output : (outputs, inputs, len(s)) complex ndarray
+            Frequency response
         """
         s_arr = np.array(s, ndmin=1) # force to be an array
         out = empty((self.outputs, self.inputs, len(s_arr)), dtype=complex)

From 4ad33db513b13f916608442712d449001918d11d Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Mon, 17 Aug 2020 11:13:06 -0700
Subject: [PATCH 10/32] suggested changes in docstrings

---
 control/frdata.py  | 4 ++--
 control/margins.py | 6 ++++--
 control/rlocus.py  | 2 +-
 3 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/control/frdata.py b/control/frdata.py
index cd0680f8d..9ece4d068 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -365,8 +365,8 @@ def eval(self, omega, squeeze=True):
             elements = np.isin(self.omega, omega) # binary array
             if sum(elements) < len(omega_array):
                 raise ValueError(
-                    '''not all frequencies omega are in frequency list of FRD 
-                    system. Try an interpolating FRD for additional points.''')
+                    "not all frequencies omega are in frequency list of FRD "
+                    "system. Try an interpolating FRD for additional points.")
             else:
                 out = self.fresp[:, :, elements]
         else:
diff --git a/control/margins.py b/control/margins.py
index cda5095aa..35178831d 100644
--- a/control/margins.py
+++ b/control/margins.py
@@ -56,6 +56,8 @@
 from . import xferfcn
 from .lti import issiso
 from . import frdata
+from .exception import ControlMIMONotImplemented
+
 
 __all__ = ['stability_margins', 'phase_crossover_frequencies', 'margin']
 
@@ -155,7 +157,7 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
 
         # check for siso
         if not issiso(sys):
-            raise ValueError("Can only do margins for SISO system")
+            raise ControlMIMONotImplemented()
 
         # real and imaginary part polynomials in omega:
         rnum, inum = _polyimsplit(sys.num[0][0])
@@ -297,7 +299,7 @@ def phase_crossover_frequencies(sys):
     """
 
     if not sys.issiso(): 
-        raise ValueError("MIMO systems not yet implemented.")
+        raise ControlMIMONotImplemented()
     # Convert to a transfer function
     tf = xferfcn._convert_to_transfer_function(sys)
     num = tf.num[0][0]
diff --git a/control/rlocus.py b/control/rlocus.py
index 39c2fcfb6..40eb0b8c8 100644
--- a/control/rlocus.py
+++ b/control/rlocus.py
@@ -472,7 +472,7 @@ def _systopoly1d(sys):
         sys = _convert_to_transfer_function(sys)
 
         # Make sure we have a SISO system
-        if (sys.inputs > 1 or sys.outputs > 1):
+        if not sys.issiso():
             raise ControlMIMONotImplemented()
 
         # Start by extracting the numerator and denominator from system object

From 60433e028d1d7feff6ed97b456b8d93c229a9d7e Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Mon, 17 Aug 2020 12:54:04 -0700
Subject: [PATCH 11/32] better docstrings conforming to numpydoc

---
 control/frdata.py  | 29 +++++++++++++++++------------
 control/lti.py     | 36 +++++++++++++++++-------------------
 control/statesp.py | 26 ++++++++++++++++++++------
 control/xferfcn.py | 32 +++++++++++++++++++++++---------
 4 files changed, 77 insertions(+), 46 deletions(-)

diff --git a/control/frdata.py b/control/frdata.py
index 0a332d44e..9f36ab8f8 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -100,6 +100,7 @@ def __init__(self, *args, **kwargs):
         object, other than an FRD, call FRD(sys, omega)
 
         """
+        # TODO: discrete-time FRD systems? 
         smooth = kwargs.get('smooth', False)
 
         if len(args) == 2:
@@ -386,29 +387,33 @@ def eval(self, omega, squeeze=True):
         return out
 
     def __call__(self, s, squeeze=True):
-        """Evaluate the system's transfer function at complex frequencies s.
+        """Evaluate system's transfer function at complex frequencies.
+        
+        Returns the complex frequency response `sys(x)` where `x` is `s` for 
+        continuous-time systems and `x` is `z` for discrete-time systems. 
 
-        For a SISO system, returns the complex value of the
-        transfer function.  For a MIMO transfer fuction, returns a
-        matrix of values.
+        To evaluate at a frequency omega in radians per second, enter 
+        x = omega*j.
         
         Parameters
         ----------
-        s : scalar or array_like
-            Complex frequencies
+        x: complex scalar or array_like 
+            Complex frequency(s) 
         squeeze: bool, optional (default=True)
-            If True and sys is single input, single output (SISO), return a 
-            1D array or scalar depending on omega's length.
-
+            If True and sys is single input single output (SISO), returns a 
+            1D array or scalar depending on the length of x. 
+            
         Returns
         -------
-        ndarray or scalar
-            Frequency response
+        fresp : (num_outputs, num_inputs, len(x)) or len(x) complex ndarray
+            The frequency response of the system. Array is len(x) if and only if
+            system is SISO and squeeze=True.
 
         Raises
         ------
         ValueError
-            If `s` is not purely imaginary. 
+            If `s` is not purely imaginary, because FrequencyDomainData systems
+            are only defined at imaginary frequency values. 
 
         """
         if any(abs(np.array(s, ndmin=1).real) > 0):
diff --git a/control/lti.py b/control/lti.py
index e89e81825..8820e2a36 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -428,23 +428,29 @@ def damp(sys, doprint=True):
 def evalfr(sys, x, squeeze=True):
     """
     Evaluate the transfer function of an LTI system for complex frequency x.
+    
+    Returns the complex frequency response `sys(x)` where `x` is `s` for 
+    continuous-time systems and `x` is `z` for discrete-time systems. 
 
-    To evaluate at a frequency, enter x = omega*j, where omega is the
-    frequency in radians per second
+    To evaluate at a frequency omega in radians per second, enter x = omega*j, 
+    for continuous-time systems, or x = exp(j*omega*dt) for discrete-time 
+    systems. 
 
     Parameters
     ----------
     sys: StateSpace or TransferFunction
         Linear system
-    x: scalar or array_like
-        Complex number
+    x: complex scalar or array_like 
+        Complex frequency(s) 
     squeeze: bool, optional (default=True)
-        If True and sys is single input, single output (SISO), return a 
-        1D array or scalar depending on omega's length.
-
+        If True and sys is single input single output (SISO), returns a 
+        1D array or scalar depending on the length of x. 
+        
     Returns
     -------
-    fresp: ndarray
+    fresp : (num_outputs, num_inputs, len(x)) or len(x) complex ndarray 
+        The frequency response of the system. If system is SISO and squeezeThe size of the array depends on 
+        whether system is SISO and squeeze keyword.
 
     See Also
     --------
@@ -453,8 +459,8 @@ def evalfr(sys, x, squeeze=True):
 
     Notes
     -----
-    This function is a wrapper for StateSpace.evalfr and
-    TransferFunction.evalfr.
+    This function is a wrapper for StateSpace.__call__ and
+    TransferFunction.__call__.
 
     Examples
     --------
@@ -465,15 +471,7 @@ def evalfr(sys, x, squeeze=True):
 
     .. todo:: Add example with MIMO system
     """
-    out = sys.horner(x)
-    if not hasattr(x, '__len__'): 
-            # received a scalar x, squeeze down the array along last dim
-            out = np.squeeze(out, axis=2)
-    if squeeze and issiso(sys):
-        return out[0][0]
-    else: 
-        return out
-
+    return sys.__call__(x, squeeze=squeeze)
 
 def freqresp(sys, omega, squeeze=True):
     """
diff --git a/control/statesp.py b/control/statesp.py
index 08eab1beb..c5685e713 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -440,12 +440,26 @@ def __rdiv__(self, other):
     def __call__(self, x, squeeze=True):
         """Evaluate system's transfer function at complex frequencies.
         
-        Evaluates at complex frequency x, where x is s or z dependong on 
-        whether the system is continuous or discrete-time.
-        
-        If squeeze is True (default) and sys is single input, single output 
-        (SISO), return a 1D array or scalar depending on the size of x. 
-        For a MIMO system, returns an (n_outputs, n_inputs, n_x) array.
+        Returns the complex frequency response `sys(x)` where `x` is `s` for 
+        continuous-time systems and `x` is `z` for discrete-time systems. 
+
+        To evaluate at a frequency omega in radians per second, enter 
+        x = omega*j, for continuous-time systems, or x = exp(j*omega*dt) for 
+        discrete-time systems. 
+
+        Parameters
+        ----------
+        x: complex scalar or array_like 
+            Complex frequency(s) 
+        squeeze: bool, optional (default=True)
+            If True and sys is single input single output (SISO), returns a 
+            1D array or scalar depending on the length of x. 
+            
+        Returns
+        -------
+        fresp : (num_outputs, num_inputs, len(x)) or len(x) complex ndarray
+            The frequency response of the system. Array is len(x) if and only if
+            system is SISO and squeeze=True.
 
         """        
         # Use Slycot if available
diff --git a/control/xferfcn.py b/control/xferfcn.py
index c8de0b265..ad527325b 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -207,14 +207,28 @@ def __init__(self, *args):
     def __call__(self, x, squeeze=True):
         """Evaluate system's transfer function at complex frequencies.
         
-        Evaluates at complex frequencies x, where x is s or z dependong on 
-        whether the system is continuous or discrete-time.
-        
-        If squeeze is True (default) and sys is single input, single output 
-        (SISO), return a 1D array or scalar depending on the shape of x. 
-        For a MIMO system, returns an (n_outputs, n_inputs, n_x) array.
+        Returns the complex frequency response `sys(x)` where `x` is `s` for 
+        continuous-time systems and `x` is `z` for discrete-time systems. 
+
+        To evaluate at a frequency omega in radians per second, enter 
+        x = omega*j, for continuous-time systems, or x = exp(j*omega*dt) for 
+        discrete-time systems. 
+
+        Parameters
+        ----------
+        x: complex scalar or array_like 
+            Complex frequency(s) 
+        squeeze: bool, optional (default=True)
+            If True and sys is single input single output (SISO), returns a 
+            1D array or scalar depending on the length of x. 
+            
+        Returns
+        -------
+        fresp : (num_outputs, num_inputs, len(x)) or len(x) complex ndarray
+            The frequency response of the system. Array is len(x) if and only if
+            system is SISO and squeeze=True.
 
-        """        
+        """     
         out = self.horner(x)
         if not hasattr(x, '__len__'): 
             # received a scalar x, squeeze down the array along last dim
@@ -226,7 +240,7 @@ def __call__(self, x, squeeze=True):
             return out
 
     def horner(self, s):
-        """Evaluate system's transfer function at complex frequencies s 
+        """Evaluates system's transfer function at complex frequencies s 
         using Horner's method.
 
         Expects inputs and outputs to be formatted correctly. Use __call__
@@ -234,7 +248,7 @@ def horner(self, s):
 
         Parameters
         ----------
-        s : array_like
+        s : array_like or scalar
             Complex frequencies
 
         Returns

From 262fde685bfd5a5ba07578142d46ac42cfcd6e3d Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Mon, 17 Aug 2020 13:32:00 -0700
Subject: [PATCH 12/32] converted statesp.horner to take care of trying to use
 slycot_horner and doing fallback

---
 control/statesp.py | 45 +++++++++++++++++++++++++++------------------
 1 file changed, 27 insertions(+), 18 deletions(-)

diff --git a/control/statesp.py b/control/statesp.py
index c5685e713..3ea597424 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -463,10 +463,7 @@ def __call__(self, x, squeeze=True):
 
         """        
         # Use Slycot if available
-        try:
-            out = self.slycot_horner(x)
-        except ImportError:  # Slycot unavailable. use built-in horner.
-            out = self.horner(x)
+        out = self.horner(x)
         if not hasattr(x, '__len__'): 
             # received a scalar x, squeeze down the array along last dim
             out = np.squeeze(out, axis=2)
@@ -526,33 +523,45 @@ def slycot_horner(self, s):
             out[:, :, kk+1] = result[0] + self.D
         return out
 
-    def horner(self, s):
-        """Evaluate system's transfer function at complex frequencies s 
-        using Horner's method.
+    def horner(self, x):
+        """Evaluate system's transfer function at complex frequencies x
+        using Horner's method. 
 
+        Evaluates sys(x) where `x` is `s` for continuous-time systems and `z` 
+        for discrete-time systems. 
+        
         Expects inputs and outputs to be formatted correctly. Use __call__
         for a more user-friendly interface. 
 
         Parameters
         ----------
-        s : array_like
+        x : array_like
             Complex frequencies
 
         Returns
         -------
         output : (outputs, inputs, len(s)) complex ndarray
             Frequency response
+
+        Notes
+        -----
+            Attempts to use Slycot library, with a fall-back to python code.
         """
-        s_arr = np.array(s, ndmin=1) # force to be an array
-        # Preallocate 
-        out = empty((self.outputs, self.inputs, len(s_arr)), dtype=complex)
-        
-        #TODO: can this be vectorized? 
-        for idx, s_idx in enumerate(s_arr):
-            out[:,:,idx] = \
-                np.dot(self.C, 
-                       solve(s_idx * eye(self.states) - self.A, self.B)) \
-                + self.D
+        try:
+            out = self.slycot_horner(x)
+        except (ImportError, Exception):  
+            # Fall back because either Slycot unavailable or cannot handle 
+            # certain cases.
+            x_arr = np.array(x, ndmin=1) # force to be an array    
+            # Preallocate 
+            out = empty((self.outputs, self.inputs, len(x_arr)), dtype=complex)
+            
+            #TODO: can this be vectorized? 
+            for idx, x_idx in enumerate(x_arr):
+                out[:,:,idx] = \
+                    np.dot(self.C, 
+                        solve(x_idx * eye(self.states) - self.A, self.B)) \
+                    + self.D
         return out
 
     def freqresp(self, omega):

From 4c3ed0b7ad8c6b1a2f41378c8800263b32eb7483 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Mon, 17 Aug 2020 14:40:30 -0700
Subject: [PATCH 13/32] renamed slycot_horner to slycot_laub

---
 control/statesp.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/control/statesp.py b/control/statesp.py
index 3ea597424..0671ed1dc 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -472,9 +472,9 @@ def __call__(self, x, squeeze=True):
         else:
             return out
 
-    def slycot_horner(self, s):
+    def slycot_laub(self, s):
         """Evaluate system's transfer function at complex frequencies s 
-        using Horner's method from Slycot.
+        using Laub's method from Slycot.
 
         Expects inputs and outputs to be formatted correctly. Use __call__
         for a more user-friendly interface. 
@@ -525,7 +525,7 @@ def slycot_horner(self, s):
 
     def horner(self, x):
         """Evaluate system's transfer function at complex frequencies x
-        using Horner's method. 
+        using Laub's or Horner's method. 
 
         Evaluates sys(x) where `x` is `s` for continuous-time systems and `z` 
         for discrete-time systems. 
@@ -545,10 +545,11 @@ def horner(self, x):
 
         Notes
         -----
-            Attempts to use Slycot library, with a fall-back to python code.
+            Attempts to use Laub's method from Slycot library, with a 
+            fall-back to python code.
         """
         try:
-            out = self.slycot_horner(x)
+            out = self.slycot_laub(x)
         except (ImportError, Exception):  
             # Fall back because either Slycot unavailable or cannot handle 
             # certain cases.

From b37c16b6d53e9ec0f7d63b2645cb3ef0a75b8ae7 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Mon, 17 Aug 2020 15:11:14 -0700
Subject: [PATCH 14/32] docs more consistent

---
 control/frdata.py  | 31 ++++++++++++++++++-------------
 control/lti.py     | 26 +++++++++++++-------------
 control/statesp.py | 40 +++++++++++++++++++---------------------
 control/xferfcn.py | 24 ++++++++++++++----------
 4 files changed, 64 insertions(+), 57 deletions(-)

diff --git a/control/frdata.py b/control/frdata.py
index 9f36ab8f8..61b70a8b6 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -344,10 +344,7 @@ def __pow__(self, other):
     # update Sawyer B. Fuller 2020.08.14: __call__ added to provide a uniform 
     # interface to systems in general and the lti.frequency_response method
     def eval(self, omega, squeeze=True):
-        """Evaluate a transfer function at a single angular frequency.
-
-        self.evalfr(omega) returns the value of the frequency response
-        at frequency omega.
+        """Evaluate a transfer function at angular frequency omega.
 
         Note that a "normal" FRD only returns values for which there is an
         entry in the omega vector. An interpolating FRD can return
@@ -355,10 +352,19 @@ def eval(self, omega, squeeze=True):
         
         Parameters
         ----------
+        omega: float or array_like
+            Frequencies in radians per second
         squeeze: bool, optional (default=True)
-            If True and sys is single input, single output (SISO), return a 
-            1D array or scalar the same length as omega.
-        """
+            If True and sys is single input single output (SISO), returns a 
+            1D array or scalar depending on the length of omega 
+            
+        Returns
+        -------
+        fresp : (num_outputs, num_inputs, len(x)) or len(x) complex ndarray
+            The frequency response of the system. Array is len(x) if and only if
+            system is SISO and squeeze=True.
+
+        """     
         omega_array = np.array(omega, ndmin=1) # array-like version of omega 
         if any(omega_array.imag > 0): 
             raise ValueError("FRD.eval can only accept real-valued omega")
@@ -389,16 +395,15 @@ def eval(self, omega, squeeze=True):
     def __call__(self, s, squeeze=True):
         """Evaluate system's transfer function at complex frequencies.
         
-        Returns the complex frequency response `sys(x)` where `x` is `s` for 
-        continuous-time systems and `x` is `z` for discrete-time systems. 
+        Returns the complex frequency response `sys(s)`.
 
         To evaluate at a frequency omega in radians per second, enter 
-        x = omega*j.
-        
+        x = omega*1j or use FRD.eval(omega)
+
         Parameters
         ----------
-        x: complex scalar or array_like 
-            Complex frequency(s) 
+        s: complex scalar or array_like 
+            Complex frequencies 
         squeeze: bool, optional (default=True)
             If True and sys is single input single output (SISO), returns a 
             1D array or scalar depending on the length of x. 
diff --git a/control/lti.py b/control/lti.py
index 8820e2a36..b0b9d7060 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -134,10 +134,10 @@ def frequency_response(self, omega, squeeze=True):
 
         Returns
         -------
-        mag : (self.outputs, self.inputs, len(omega)) ndarray
+        mag : (self.outputs, self.inputs, len(omega)) or len(omega) ndarray
             The magnitude (absolute value, not dB or log10) of the system
             frequency response.
-        phase : (self.outputs, self.inputs, len(omega)) ndarray
+        phase : (self.outputs, self.inputs, len(omega)) or len(omega) ndarray
             The wrapped phase in radians of the system frequency response.
         omega : ndarray
             The (sorted) frequencies at which the response was evaluated.
@@ -430,10 +430,10 @@ def evalfr(sys, x, squeeze=True):
     Evaluate the transfer function of an LTI system for complex frequency x.
     
     Returns the complex frequency response `sys(x)` where `x` is `s` for 
-    continuous-time systems and `x` is `z` for discrete-time systems. 
+    continuous-time systems and `z` for discrete-time systems. 
 
-    To evaluate at a frequency omega in radians per second, enter x = omega*j, 
-    for continuous-time systems, or x = exp(j*omega*dt) for discrete-time 
+    To evaluate at a frequency omega in radians per second, enter x = omega*1j, 
+    for continuous-time systems, or x = exp(1j*omega*dt) for discrete-time 
     systems. 
 
     Parameters
@@ -448,9 +448,9 @@ def evalfr(sys, x, squeeze=True):
         
     Returns
     -------
-    fresp : (num_outputs, num_inputs, len(x)) or len(x) complex ndarray 
-        The frequency response of the system. If system is SISO and squeezeThe size of the array depends on 
-        whether system is SISO and squeeze keyword.
+    fresp : (sys.outputs, sys.inputs, len(x)) or len(x) complex ndarray 
+        The frequency response of the system. Array is len(x) if and only if
+        system is SISO and squeeze=True.
 
     See Also
     --------
@@ -481,22 +481,22 @@ def freqresp(sys, omega, squeeze=True):
     ----------
     sys: StateSpace or TransferFunction
         Linear system
-    omega: array_like
+    omega: float or array_like
         A list of frequencies in radians/sec at which the system should be
         evaluated. The list can be either a python list or a numpy array
         and will be sorted before evaluation.
     squeeze: bool, optional (default=True)
-        If True and sys is single input, single output (SISO), return a 
+        If True and sys is single input, single output (SISO), returns  
         1D array or scalar depending on omega's length.
 
     Returns
     -------
-    mag : (self.outputs, self.inputs, len(omega)) ndarray
+    mag : (sys.outputs, sys.inputs, len(omega)) or len(omega) ndarray 
         The magnitude (absolute value, not dB or log10) of the system
         frequency response.
-    phase : (self.outputs, self.inputs, len(omega)) ndarray
+    phase : (sys.outputs, sys.inputs, len(omega)) or len(omega) ndarray
         The wrapped phase in radians of the system frequency response.
-    omega : ndarray or list or tuple
+    omega : ndarray
         The list of sorted frequencies at which the response was
         evaluated.
 
diff --git a/control/statesp.py b/control/statesp.py
index 0671ed1dc..2b9ef016a 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -438,26 +438,26 @@ def __rdiv__(self, other):
         raise NotImplementedError("StateSpace.__rdiv__ is not implemented yet.")
 
     def __call__(self, x, squeeze=True):
-        """Evaluate system's transfer function at complex frequencies.
+        """Evaluate system's transfer function at complex frequency.
         
         Returns the complex frequency response `sys(x)` where `x` is `s` for 
-        continuous-time systems and `x` is `z` for discrete-time systems. 
+        continuous-time systems and `z` for discrete-time systems. 
 
         To evaluate at a frequency omega in radians per second, enter 
-        x = omega*j, for continuous-time systems, or x = exp(j*omega*dt) for 
+        x = omega*1j, for continuous-time systems, or x = exp(1j*omega*dt) for 
         discrete-time systems. 
 
         Parameters
         ----------
-        x: complex scalar or array_like 
-            Complex frequency(s) 
+        x: complex or complex array_like 
+            Complex frequencies
         squeeze: bool, optional (default=True)
             If True and sys is single input single output (SISO), returns a 
             1D array or scalar depending on the length of x. 
             
         Returns
         -------
-        fresp : (num_outputs, num_inputs, len(x)) or len(x) complex ndarray
+        fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
             The frequency response of the system. Array is len(x) if and only if
             system is SISO and squeeze=True.
 
@@ -472,8 +472,8 @@ def __call__(self, x, squeeze=True):
         else:
             return out
 
-    def slycot_laub(self, s):
-        """Evaluate system's transfer function at complex frequencies s 
+    def slycot_laub(self, x):
+        """Evaluate system's transfer function at complex frequency
         using Laub's method from Slycot.
 
         Expects inputs and outputs to be formatted correctly. Use __call__
@@ -481,27 +481,25 @@ def slycot_laub(self, s):
 
         Parameters
         ----------
-        s : array_like
+        x : complex array_like or complex
             Complex frequency
 
         Returns
         -------
-        output : (outputs, inputs, len(s)) complex ndarray
+        output : (number_outputs, number_inputs, len(x)) complex ndarray
             Frequency response
         """
         from slycot import tb05ad
 
         # preallocate
-        s_arr = np.array(s, ndmin=1) # array-like version of s
-        out = np.empty((self.outputs, self.inputs, len(s_arr)), 
-                        dtype=complex)
+        x_arr = np.array(x, ndmin=1) # array-like version of s
         n = self.states
         m = self.inputs
         p = self.outputs
+        out = np.empty((p, m, len(x_arr)), dtype=complex)
         # The first call both evaluates C(sI-A)^-1 B and also returns
         # Hessenberg transformed matrices at, bt, ct.
-        result = tb05ad(n, m, p, s_arr[0], self.A,
-                        self.B, self.C, job='NG')
+        result = tb05ad(n, m, p, x_arr[0], self.A, self.B, self.C, job='NG')
         # When job='NG', result = (at, bt, ct, g_i, hinvb, info)
         at = result[0]
         bt = result[1]
@@ -514,8 +512,8 @@ def slycot_laub(self, s):
         # transformed state matrices, at, bt, ct.
 
         # Start at the second frequency, already have the first.
-        for kk, s_kk in enumerate(s_arr[1:len(s_arr)]):
-            result = tb05ad(n, m, p, s_kk, at, bt, ct, job='NH')
+        for kk, x_kk in enumerate(x_arr[1:len(x_arr)]):
+            result = tb05ad(n, m, p, x_kk, at, bt, ct, job='NH')
             # When job='NH', result = (g_i, hinvb, info)
 
             # kk+1 because enumerate starts at kk = 0.
@@ -524,10 +522,10 @@ def slycot_laub(self, s):
         return out
 
     def horner(self, x):
-        """Evaluate system's transfer function at complex frequencies x
+        """Evaluate system's transfer function at complex frequency
         using Laub's or Horner's method. 
 
-        Evaluates sys(x) where `x` is `s` for continuous-time systems and `z` 
+        Evaluates `sys(x)` where `x` is `s` for continuous-time systems and `z` 
         for discrete-time systems. 
         
         Expects inputs and outputs to be formatted correctly. Use __call__
@@ -535,12 +533,12 @@ def horner(self, x):
 
         Parameters
         ----------
-        x : array_like
+        x : complex array_like or complex
             Complex frequencies
 
         Returns
         -------
-        output : (outputs, inputs, len(s)) complex ndarray
+        output : (number_outputs, number_inputs, len(x)) complex ndarray
             Frequency response
 
         Notes
diff --git a/control/xferfcn.py b/control/xferfcn.py
index ad527325b..3c898a320 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -208,23 +208,23 @@ def __call__(self, x, squeeze=True):
         """Evaluate system's transfer function at complex frequencies.
         
         Returns the complex frequency response `sys(x)` where `x` is `s` for 
-        continuous-time systems and `x` is `z` for discrete-time systems. 
+        continuous-time systems and `z` for discrete-time systems. 
 
         To evaluate at a frequency omega in radians per second, enter 
-        x = omega*j, for continuous-time systems, or x = exp(j*omega*dt) for 
+        x = omega*1j, for continuous-time systems, or x = exp(1j*omega*dt) for 
         discrete-time systems. 
 
         Parameters
         ----------
-        x: complex scalar or array_like 
-            Complex frequency(s) 
+        x: complex array_like or complex
+            Complex frequencies
         squeeze: bool, optional (default=True)
             If True and sys is single input single output (SISO), returns a 
             1D array or scalar depending on the length of x. 
             
         Returns
         -------
-        fresp : (num_outputs, num_inputs, len(x)) or len(x) complex ndarray
+        fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
             The frequency response of the system. Array is len(x) if and only if
             system is SISO and squeeze=True.
 
@@ -239,22 +239,26 @@ def __call__(self, x, squeeze=True):
         else:
             return out
 
-    def horner(self, s):
-        """Evaluates system's transfer function at complex frequencies s 
-        using Horner's method.
+    def horner(self, x):
+        """Evaluate system's transfer function at complex frequency
+        using Horner's method. 
 
+        Evaluates `sys(x)` where `x` is `s` for continuous-time systems and `z`
+        for discrete-time systems. 
+        
         Expects inputs and outputs to be formatted correctly. Use __call__
         for a more user-friendly interface. 
 
         Parameters
         ----------
-        s : array_like or scalar
+        x : complex array_like or complex
             Complex frequencies
 
         Returns
         -------
-        output : (outputs, inputs, len(s)) complex ndarray
+        output : (self.outputs, self.inputs, len(x)) complex ndarray
             Frequency response
+
         """
         s_arr = np.array(s, ndmin=1) # force to be an array
         out = empty((self.outputs, self.inputs, len(s_arr)), dtype=complex)

From 4e4b97f05d0bb17b7682f3c238deada9037cb3a1 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Mon, 17 Aug 2020 15:51:06 -0700
Subject: [PATCH 15/32] forgot to change a variable name everywhere in a
 function

---
 control/xferfcn.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/control/xferfcn.py b/control/xferfcn.py
index 3c898a320..abd6338dd 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -260,12 +260,12 @@ def horner(self, x):
             Frequency response
 
         """
-        s_arr = np.array(s, ndmin=1) # force to be an array
+        s_arr = np.array(x, ndmin=1) # force to be an array
         out = empty((self.outputs, self.inputs, len(s_arr)), dtype=complex)
         for i in range(self.outputs):
             for j in range(self.inputs):
-                out[i][j] = (polyval(self.num[i][j], s) /
-                             polyval(self.den[i][j], s))
+                out[i][j] = (polyval(self.num[i][j], x) /
+                             polyval(self.den[i][j], x))
         return out
 
     def _truncatecoeff(self):

From 5626a3a2744744c2342346a73991b1dab933ea5b Mon Sep 17 00:00:00 2001
From: Ben Greiner <code@bnavigator.de>
Date: Tue, 18 Aug 2020 12:47:06 +0200
Subject: [PATCH 16/32] revert doc/control.rst

common_timebase is no longer existent
---
 doc/control.rst | 1 -
 1 file changed, 1 deletion(-)

diff --git a/doc/control.rst b/doc/control.rst
index 0705e1b9f..57d64b1eb 100644
--- a/doc/control.rst
+++ b/doc/control.rst
@@ -172,7 +172,6 @@ Utility functions and conversions
     tfdata
     timebase
     timebaseEqual
-    common_timebase
     unwrap
     use_fbs_defaults
     use_matlab_defaults

From 5caffa414afbab1dafee8537bf54a26787bb5eaf Mon Sep 17 00:00:00 2001
From: sawyerbfuller <58706249+sawyerbfuller@users.noreply.github.com>
Date: Wed, 19 Aug 2020 09:38:49 -0700
Subject: [PATCH 17/32] Update control/xferfcn.py

numpydoc convention docstrings

Co-authored-by: Ben Greiner <code@bnavigator.de>
---
 control/xferfcn.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/control/xferfcn.py b/control/xferfcn.py
index abd6338dd..1f7fd1243 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -211,8 +211,8 @@ def __call__(self, x, squeeze=True):
         continuous-time systems and `z` for discrete-time systems. 
 
         To evaluate at a frequency omega in radians per second, enter 
-        x = omega*1j, for continuous-time systems, or x = exp(1j*omega*dt) for 
-        discrete-time systems. 
+        ``x = omega*1j``, for continuous-time systems, or ``x = exp(1j*omega*dt)``
+        for discrete-time systems.
 
         Parameters
         ----------

From c888b6ee86554e94a8393bb1afcf676a24dc91a6 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Wed, 19 Aug 2020 10:27:57 -0700
Subject: [PATCH 18/32] numpydoc convention updates

---
 control/frdata.py  | 40 ++++++++++++++++++++++++----------------
 control/lti.py     |  7 ++++---
 control/statesp.py | 41 ++++++++++++++++++++++++-----------------
 control/xferfcn.py | 33 ++++++++++++++++++++-------------
 4 files changed, 72 insertions(+), 49 deletions(-)

diff --git a/control/frdata.py b/control/frdata.py
index 61b70a8b6..ffb83c73e 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -355,14 +355,14 @@ def eval(self, omega, squeeze=True):
         omega: float or array_like
             Frequencies in radians per second
         squeeze: bool, optional (default=True)
-            If True and sys is single input single output (SISO), returns a 
-            1D array or scalar depending on the length of omega 
+            If True and `sys` is single input single output (SISO), returns a 
+            1D array or scalar depending on the length of `omega` 
             
         Returns
         -------
-        fresp : (num_outputs, num_inputs, len(x)) or len(x) complex ndarray
-            The frequency response of the system. Array is len(x) if and only if
-            system is SISO and squeeze=True.
+        fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
+            The frequency response of the system. Array is ``len(x)`` if and only
+            if system is SISO and ``squeeze=True``.
 
         """     
         omega_array = np.array(omega, ndmin=1) # array-like version of omega 
@@ -398,27 +398,28 @@ def __call__(self, s, squeeze=True):
         Returns the complex frequency response `sys(s)`.
 
         To evaluate at a frequency omega in radians per second, enter 
-        x = omega*1j or use FRD.eval(omega)
+        ``x = omega * 1j`` or use ``sys.eval(omega)``
 
         Parameters
         ----------
         s: complex scalar or array_like 
             Complex frequencies 
         squeeze: bool, optional (default=True)
-            If True and sys is single input single output (SISO), returns a 
-            1D array or scalar depending on the length of x. 
+            If True and `sys` is single input single output (SISO), returns a 
+            1D array or scalar depending on the length of x`. 
             
         Returns
         -------
-        fresp : (num_outputs, num_inputs, len(x)) or len(x) complex ndarray
-            The frequency response of the system. Array is len(x) if and only if
-            system is SISO and squeeze=True.
+        fresp : (self.outputs, self.inputs, len(s)) or len(s) complex ndarray
+            The frequency response of the system. Array is ``len(s)`` if and 
+            only if system is SISO and ``squeeze=True``.
 
         Raises
         ------
         ValueError
-            If `s` is not purely imaginary, because FrequencyDomainData systems
-            are only defined at imaginary frequency values. 
+            If `s` is not purely imaginary, because 
+            :class:`FrequencyDomainData` systems are only defined at imaginary 
+            frequency values. 
 
         """
         if any(abs(np.array(s, ndmin=1).real) > 0):
@@ -431,11 +432,18 @@ def __call__(self, s, squeeze=True):
             return self.eval(complex(s).imag, squeeze=squeeze)
         
     def freqresp(self, omega):
-        warn("FrequencyResponseData.freqresp(omega) will be deprecated in a "
+        """(deprecated) Evaluate transfer function at complex frequencies. 
+        
+        .. deprecated::0.9.0 
+            Method has been given the more pythonic name 
+            :meth:`FrequencyResponseData.frequency_response`. Or use 
+            :func:`freqresp` in the MATLAB compatibility module.
+        """
+        warn("FrequencyResponseData.freqresp(omega) will be removed in a "
              "future release of python-control; use "
-             "FrequencyResponseData.frequency_response(sys, omega), or "
+             "FrequencyResponseData.frequency_response(omega), or "
              "freqresp(sys, omega) in the MATLAB compatibility module "
-             "instead", PendingDeprecationWarning)        
+             "instead", DeprecationWarning)        
         return self.frequency_response(omega)
 
     def feedback(self, other=1, sign=-1):
diff --git a/control/lti.py b/control/lti.py
index b0b9d7060..4af9dfc86 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -432,9 +432,10 @@ def evalfr(sys, x, squeeze=True):
     Returns the complex frequency response `sys(x)` where `x` is `s` for 
     continuous-time systems and `z` for discrete-time systems. 
 
-    To evaluate at a frequency omega in radians per second, enter x = omega*1j, 
-    for continuous-time systems, or x = exp(1j*omega*dt) for discrete-time 
-    systems. 
+    To evaluate at a frequency omega in radians per second, enter 
+    ``x = omega * 1j`` for continuous-time systems, or 
+    ``x = exp(1j * omega * dt)`` for discrete-time systems, or use 
+    ``freqresp(sys, omega)``. 
 
     Parameters
     ----------
diff --git a/control/statesp.py b/control/statesp.py
index 2b9ef016a..d131e1aca 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -444,22 +444,23 @@ def __call__(self, x, squeeze=True):
         continuous-time systems and `z` for discrete-time systems. 
 
         To evaluate at a frequency omega in radians per second, enter 
-        x = omega*1j, for continuous-time systems, or x = exp(1j*omega*dt) for 
-        discrete-time systems. 
+        ``x = omega * 1j``, for continuous-time systems, or 
+        ``x = exp(1j * omega * dt)`` for discrete-time systems. Or use 
+        :meth:`StateSpace.frequency_response`. 
 
         Parameters
         ----------
         x: complex or complex array_like 
             Complex frequencies
         squeeze: bool, optional (default=True)
-            If True and sys is single input single output (SISO), returns a 
-            1D array or scalar depending on the length of x. 
+            If True and `sys` is single input single output (SISO), returns a 
+            1D array or scalar depending on the length of `x`. 
             
         Returns
         -------
         fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
-            The frequency response of the system. Array is len(x) if and only if
-            system is SISO and squeeze=True.
+            The frequency response of the system. Array is ``len(x)`` if and 
+            only if system is SISO and ``squeeze=True``.
 
         """        
         # Use Slycot if available
@@ -476,7 +477,7 @@ def slycot_laub(self, x):
         """Evaluate system's transfer function at complex frequency
         using Laub's method from Slycot.
 
-        Expects inputs and outputs to be formatted correctly. Use __call__
+        Expects inputs and outputs to be formatted correctly. Use ``sys(x)``
         for a more user-friendly interface. 
 
         Parameters
@@ -492,7 +493,7 @@ def slycot_laub(self, x):
         from slycot import tb05ad
 
         # preallocate
-        x_arr = np.array(x, ndmin=1) # array-like version of s
+        x_arr = np.atleast_1d(x) # array-like version of x
         n = self.states
         m = self.inputs
         p = self.outputs
@@ -528,7 +529,7 @@ def horner(self, x):
         Evaluates `sys(x)` where `x` is `s` for continuous-time systems and `z` 
         for discrete-time systems. 
         
-        Expects inputs and outputs to be formatted correctly. Use __call__
+        Expects inputs and outputs to be formatted correctly. Use ``sys(x)``
         for a more user-friendly interface. 
 
         Parameters
@@ -538,20 +539,20 @@ def horner(self, x):
 
         Returns
         -------
-        output : (number_outputs, number_inputs, len(x)) complex ndarray
+        output : (self.outputs, self.inputs, len(x)) complex ndarray
             Frequency response
 
         Notes
         -----
-            Attempts to use Laub's method from Slycot library, with a 
-            fall-back to python code.
+        Attempts to use Laub's method from Slycot library, with a 
+        fall-back to python code.
         """
         try:
             out = self.slycot_laub(x)
         except (ImportError, Exception):  
             # Fall back because either Slycot unavailable or cannot handle 
             # certain cases.
-            x_arr = np.array(x, ndmin=1) # force to be an array    
+            x_arr = np.atleast_1d(x) # force to be an array    
             # Preallocate 
             out = empty((self.outputs, self.inputs, len(x_arr)), dtype=complex)
             
@@ -564,11 +565,17 @@ def horner(self, x):
         return out
 
     def freqresp(self, omega):
-        warn("StateSpace.freqresp(omega) will be deprecated in a "
+        """(deprecated) Evaluate transfer function at complex frequencies. 
+        
+        .. deprecated::0.9.0 
+            Method has been given the more pythonic name 
+            :meth:`StateSpace.frequency_response`. Or use 
+            :func:`freqresp` in the MATLAB compatibility module.
+        """
+        warn("StateSpace.freqresp(omega) will be removed in a "
              "future release of python-control; use "
-             "StateSpace.frequency_response(sys, omega), or "
-             "freqresp(sys, omega) in the MATLAB compatibility module "
-             "instead", PendingDeprecationWarning)        
+             "sys.frequency_response(omega), or freqresp(sys, omega) in the "
+             "MATLAB compatibility module instead", DeprecationWarning)        
         return self.frequency_response(omega)
 
     # Compute poles and zeros
diff --git a/control/xferfcn.py b/control/xferfcn.py
index abd6338dd..9613fa18c 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -211,22 +211,23 @@ def __call__(self, x, squeeze=True):
         continuous-time systems and `z` for discrete-time systems. 
 
         To evaluate at a frequency omega in radians per second, enter 
-        x = omega*1j, for continuous-time systems, or x = exp(1j*omega*dt) for 
-        discrete-time systems. 
+        ``x = omega * 1j``, for continuous-time systems, or 
+        ``x = exp(1j * omega * dt)`` for discrete-time systems. Or use 
+        :meth:`TransferFunction.frequency_response`. 
 
         Parameters
         ----------
         x: complex array_like or complex
             Complex frequencies
         squeeze: bool, optional (default=True)
-            If True and sys is single input single output (SISO), returns a 
-            1D array or scalar depending on the length of x. 
+            If True and `sys` is single input single output (SISO), returns a 
+            1D array or scalar depending on the length of `x`. 
             
         Returns
         -------
         fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
-            The frequency response of the system. Array is len(x) if and only if
-            system is SISO and squeeze=True.
+            The frequency response of the system. Array is `len(x)` if and 
+            only if system is SISO and ``squeeze=True``.
 
         """     
         out = self.horner(x)
@@ -246,7 +247,7 @@ def horner(self, x):
         Evaluates `sys(x)` where `x` is `s` for continuous-time systems and `z`
         for discrete-time systems. 
         
-        Expects inputs and outputs to be formatted correctly. Use __call__
+        Expects inputs and outputs to be formatted correctly. Use ``sys(x)``
         for a more user-friendly interface. 
 
         Parameters
@@ -260,8 +261,8 @@ def horner(self, x):
             Frequency response
 
         """
-        s_arr = np.array(x, ndmin=1) # force to be an array
-        out = empty((self.outputs, self.inputs, len(s_arr)), dtype=complex)
+        x_arr = np.atleast_1d(x) # force to be an array
+        out = empty((self.outputs, self.inputs, len(x_arr)), dtype=complex)
         for i in range(self.outputs):
             for j in range(self.inputs):
                 out[i][j] = (polyval(self.num[i][j], x) /
@@ -659,11 +660,17 @@ def __getitem__(self, key):
             return TransferFunction(num, den, self.dt)
 
     def freqresp(self, omega):
-        warn("TransferFunction.freqresp(omega) will be deprecated in a "
+        """(deprecated) Evaluate transfer function at complex frequencies. 
+        
+        .. deprecated::0.9.0 
+            Method has been given the more pythonic name 
+            :meth:`TransferFunction.frequency_response`. Or use 
+            :func:`freqresp` in the MATLAB compatibility module.
+        """
+        warn("TransferFunction.freqresp(omega) will be removed in a "
              "future release of python-control; use "
-             "TransferFunction.frequency_response(sys, omega), or "
-             "freqresp(sys, omega) in the MATLAB compatibility module "
-             "instead", PendingDeprecationWarning)        
+             "sys.frequency_response(omega), or freqresp(sys, omega) in the "
+             "MATLAB compatibility module instead", DeprecationWarning)        
         return self.frequency_response(omega)
 
     def pole(self):

From 3cf80afbbb3bd9ccf4ce7a761709f7c2535756a0 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Wed, 19 Aug 2020 10:35:17 -0700
Subject: [PATCH 19/32] small numpydoc change in trasnferfunction to resolve
 merge

---
 control/xferfcn.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/control/xferfcn.py b/control/xferfcn.py
index cfe963e11..9613fa18c 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -211,8 +211,9 @@ def __call__(self, x, squeeze=True):
         continuous-time systems and `z` for discrete-time systems. 
 
         To evaluate at a frequency omega in radians per second, enter 
-        ``x = omega*1j``, for continuous-time systems, or ``x = exp(1j*omega*dt)``
-        for discrete-time systems.
+        ``x = omega * 1j``, for continuous-time systems, or 
+        ``x = exp(1j * omega * dt)`` for discrete-time systems. Or use 
+        :meth:`TransferFunction.frequency_response`. 
 
         Parameters
         ----------

From 330e517002114e00cd07fc00f978a1be3c6b2fc4 Mon Sep 17 00:00:00 2001
From: Sawyer Fuller <minster@uw.edu>
Date: Wed, 19 Aug 2020 10:44:32 -0700
Subject: [PATCH 20/32] unit tests now test for deprecation of sys.freqresp

---
 control/tests/statesp_test.py |  2 +-
 control/tests/xferfcn_test.py | 13 +++++++++++++
 2 files changed, 14 insertions(+), 1 deletion(-)

diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 77a6d3ec2..e96a91014 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -223,7 +223,7 @@ def test_freqresp_deprecated(self):
 
             # Make sure that we get a pending deprecation warning
             sys.freqresp(1.)
-            assert issubclass(w[-1].category, PendingDeprecationWarning)
+            assert issubclass(w[-1].category, DeprecationWarning)
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
     def test_freq_resp(self):
diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py
index 4ee314df1..140979a79 100644
--- a/control/tests/xferfcn_test.py
+++ b/control/tests/xferfcn_test.py
@@ -432,6 +432,19 @@ def test_call_mimo(self):
         # Test call version as well
         np.testing.assert_array_almost_equal(sys(2.j), resp)
 
+    def test_freqresp_deprecated(self):
+        sys = TransferFunction([1., 3., 5], [1., 6., 2., -1.])
+        # Deprecated version of the call (should generate warning)
+        import warnings
+        with warnings.catch_warnings(record=True) as w:
+            # Set up warnings filter to only show warnings in control module
+            warnings.filterwarnings("ignore")
+            warnings.filterwarnings("always", module="control")
+
+            # Make sure that we get a pending deprecation warning
+            sys.freqresp(1.)
+            assert issubclass(w[-1].category, DeprecationWarning)
+
     def test_frequency_response_siso(self):
         """Evaluate the magnitude and phase of a SISO system at
         multiple frequencies."""

From 195bec3effa50052518ebf4cd042c846eecdd964 Mon Sep 17 00:00:00 2001
From: Richard Murray <murray@cds.caltech.edu>
Date: Tue, 5 Jan 2021 22:20:28 -0800
Subject: [PATCH 21/32] rebase sawyerbfuller-call-method against master

---
 control/bench/time_freqresp.py |   4 +-
 control/frdata.py              | 192 ++++++++++++++--------------
 control/freqplot.py            |  24 ++--
 control/lti.py                 | 141 ++++++++++++++-------
 control/margins.py             |  14 +-
 control/matlab/__init__.py     |   4 +-
 control/nichols.py             |   2 +-
 control/rlocus.py              |  12 +-
 control/statesp.py             | 225 ++++++++++++++++++---------------
 control/tests/frd_test.py      | 158 +++++++++++------------
 control/tests/freqresp_test.py |  12 +-
 control/tests/iosys_test.py    |   4 +-
 control/tests/margin_test.py   |  12 +-
 control/tests/statesp_test.py  |  45 +++++--
 control/tests/xferfcn_test.py  |  84 +++++++-----
 control/xferfcn.py             | 175 +++++++++++--------------
 examples/robust_mimo.py        |   2 +-
 examples/robust_siso.py        |   8 +-
 18 files changed, 594 insertions(+), 524 deletions(-)

diff --git a/control/bench/time_freqresp.py b/control/bench/time_freqresp.py
index 1945cbc24..3ae837082 100644
--- a/control/bench/time_freqresp.py
+++ b/control/bench/time_freqresp.py
@@ -8,7 +8,7 @@
 sys_tf = tf(sys)
 w = logspace(-1,1,50)
 ntimes = 1000
-time_ss = timeit("sys.freqresp(w)", setup="from __main__ import sys, w", number=ntimes)
-time_tf = timeit("sys_tf.freqresp(w)", setup="from __main__ import sys_tf, w", number=ntimes)
+time_ss = timeit("sys.freqquency_response(w)", setup="from __main__ import sys, w", number=ntimes)
+time_tf = timeit("sys_tf.frequency_response(w)", setup="from __main__ import sys_tf, w", number=ntimes)
 print("State-space model on %d states: %f" % (nstates, time_ss))
 print("Transfer-function model on %d states: %f" % (nstates, time_tf))
diff --git a/control/frdata.py b/control/frdata.py
index 8ca9dfd9d..95c84e0ba 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -48,7 +48,7 @@
 from warnings import warn
 import numpy as np
 from numpy import angle, array, empty, ones, \
-    real, imag, absolute, eye, linalg, where, dot
+    real, imag, absolute, eye, linalg, where, dot, sort
 from scipy.interpolate import splprep, splev
 from .lti import LTI
 
@@ -100,6 +100,7 @@ def __init__(self, *args, **kwargs):
         object, other than an FRD, call FRD(sys, omega)
 
         """
+        # TODO: discrete-time FRD systems?
         smooth = kwargs.get('smooth', False)
 
         if len(args) == 2:
@@ -107,17 +108,10 @@ def __init__(self, *args, **kwargs):
                 # not an FRD, but still a system, second argument should be
                 # the frequency range
                 otherlti = args[0]
-                self.omega = array(args[1], dtype=float)
-                self.omega.sort()
+                self.omega = sort(np.asarray(args[1], dtype=float))
                 numfreq = len(self.omega)
-
                 # calculate frequency response at my points
-                self.fresp = empty(
-                    (otherlti.outputs, otherlti.inputs, numfreq),
-                    dtype=complex)
-                for k, w in enumerate(self.omega):
-                    self.fresp[:, :, k] = otherlti._evalfr(w)
-
+                self.fresp = otherlti(1j * self.omega, squeeze=False)
             else:
                 # The user provided a response and a freq vector
                 self.fresp = array(args[0], dtype=complex)
@@ -141,7 +135,7 @@ def __init__(self, *args, **kwargs):
             self.omega = args[0].omega
             self.fresp = args[0].fresp
         else:
-            raise ValueError("Needs 1 or 2 arguments; receivd %i." % len(args))
+            raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
 
         # create interpolation functions
         if smooth:
@@ -163,7 +157,6 @@ def __str__(self):
         mimo = self.inputs > 1 or self.outputs > 1
         outstr = ['Frequency response data']
 
-        mt, pt, wt = self.freqresp(self.omega)
         for i in range(self.inputs):
             for j in range(self.outputs):
                 if mimo:
@@ -171,9 +164,10 @@ def __str__(self):
                 outstr.append('Freq [rad/s]  Response')
                 outstr.append('------------  ---------------------')
                 outstr.extend(
-                    ['%12.3f  %10.4g%+10.4gj' % (w, m, p)
-                     for m, p, w in zip(real(self.fresp[j, i, :]),
-                                        imag(self.fresp[j, i, :]), wt)])
+                    ['%12.3f  %10.4g%+10.4gj' % (w, re, im)
+                     for w, re, im in zip(self.omega,
+                                          real(self.fresp[j, i, :]),
+                                          imag(self.fresp[j, i, :]))])
 
         return '\n'.join(outstr)
 
@@ -342,110 +336,115 @@ def __pow__(self, other):
             return (FRD(ones(self.fresp.shape), self.omega) / self) * \
                 (self**(other+1))
 
-    def evalfr(self, omega):
-        """Evaluate a transfer function at a single angular frequency.
-
-        self._evalfr(omega) returns the value of the frequency response
-        at frequency omega.
-
-        Note that a "normal" FRD only returns values for which there is an
-        entry in the omega vector. An interpolating FRD can return
-        intermediate values.
-
-        """
-        warn("FRD.evalfr(omega) will be deprecated in a future release "
-             "of python-control; use sys.eval(omega) instead",
-             PendingDeprecationWarning)         # pragma: no coverage
-        return self._evalfr(omega)
-
     # Define the `eval` function to evaluate an FRD at a given (real)
     # frequency.  Note that we choose to use `eval` instead of `evalfr` to
     # avoid confusion with :func:`evalfr`, which takes a complex number as its
     # argument.  Similarly, we don't use `__call__` to avoid confusion between
     # G(s) for a transfer function and G(omega) for an FRD object.
-    def eval(self, omega):
-        """Evaluate a transfer function at a single angular frequency.
-
-        self.evalfr(omega) returns the value of the frequency response
-        at frequency omega.
+    # update Sawyer B. Fuller 2020.08.14: __call__ added to provide a uniform
+    # interface to systems in general and the lti.frequency_response method
+    def eval(self, omega, squeeze=True):
+        """Evaluate a transfer function at angular frequency omega.
 
         Note that a "normal" FRD only returns values for which there is an
         entry in the omega vector. An interpolating FRD can return
         intermediate values.
 
+        Parameters
+        ----------
+        omega: float or array_like
+            Frequencies in radians per second
+        squeeze: bool, optional (default=True)
+            If True and `sys` is single input single output (SISO), returns a
+            1D array or scalar depending on the length of `omega`
+
+        Returns
+        -------
+        fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
+            The frequency response of the system. Array is ``len(x)`` if and
+            only if system is SISO and ``squeeze=True``.
+
         """
-        return self._evalfr(omega)
-
-    # Internal function to evaluate the frequency responses
-    def _evalfr(self, omega):
-        """Evaluate a transfer function at a single angular frequency."""
-        # Preallocate the output.
-        if getattr(omega, '__iter__', False):
-            out = empty((self.outputs, self.inputs, len(omega)), dtype=complex)
-        else:
-            out = empty((self.outputs, self.inputs), dtype=complex)
+        omega_array = np.array(omega, ndmin=1) # array-like version of omega
+        if any(omega_array.imag > 0):
+            raise ValueError("FRD.eval can only accept real-valued omega")
 
         if self.ifunc is None:
-            try:
-                out = self.fresp[:, :, where(self.omega == omega)[0][0]]
-            except Exception:
+            elements = np.isin(self.omega, omega) # binary array
+            if sum(elements) < len(omega_array):
                 raise ValueError(
-                    "Frequency %f not in frequency list, try an interpolating"
-                    " FRD if you want additional points" % omega)
-        else:
-            if getattr(omega, '__iter__', False):
-                for i in range(self.outputs):
-                    for j in range(self.inputs):
-                        for k, w in enumerate(omega):
-                            frraw = splev(w, self.ifunc[i, j], der=0)
-                            out[i, j, k] = frraw[0] + 1.0j * frraw[1]
+                    "not all frequencies omega are in frequency list of FRD "
+                    "system. Try an interpolating FRD for additional points.")
             else:
-                for i in range(self.outputs):
-                    for j in range(self.inputs):
-                        frraw = splev(omega, self.ifunc[i, j], der=0)
-                        out[i, j] = frraw[0] + 1.0j * frraw[1]
-
+                out = self.fresp[:, :, elements]
+        else:
+            out = empty((self.outputs, self.inputs, len(omega_array)),
+                         dtype=complex)
+            for i in range(self.outputs):
+                for j in range(self.inputs):
+                    for k, w in enumerate(omega_array):
+                        frraw = splev(w, self.ifunc[i, j], der=0)
+                        out[i, j, k] = frraw[0] + 1.0j * frraw[1]
+        if not hasattr(omega, '__len__'):
+            # omega is a scalar, squeeze down array along last dim
+            out = np.squeeze(out, axis=2)
+        if squeeze and self.issiso():
+            out = out[0][0]
         return out
 
-    # Method for generating the frequency response of the system
-    def freqresp(self, omega):
-        """Evaluate the frequency response at a list of angular frequencies.
+    def __call__(self, s, squeeze=True):
+        """Evaluate system's transfer function at complex frequencies.
 
-        Reports the value of the magnitude, phase, and angular frequency of
-        the requency response evaluated at omega, where omega is a list of
-        angular frequencies, and is a sorted version of the input omega.
+        Returns the complex frequency response `sys(s)`.
+
+        To evaluate at a frequency omega in radians per second, enter
+        ``x = omega * 1j`` or use ``sys.eval(omega)``
 
         Parameters
         ----------
-        omega : array_like
-            A list of frequencies in radians/sec at which the system should be
-            evaluated. The list can be either a python list or a numpy array
-            and will be sorted before evaluation.
+        s: complex scalar or array_like
+            Complex frequencies
+        squeeze: bool, optional (default=True)
+            If True and `sys` is single input single output (SISO), returns a
+            1D array or scalar depending on the length of x`.
 
         Returns
         -------
-        mag : (self.outputs, self.inputs, len(omega)) ndarray
-            The magnitude (absolute value, not dB or log10) of the system
-            frequency response.
-        phase : (self.outputs, self.inputs, len(omega)) ndarray
-            The wrapped phase in radians of the system frequency response.
-        omega : ndarray or list or tuple
-            The list of sorted frequencies at which the response was
-            evaluated.
-        """
-        # Preallocate outputs.
-        numfreq = len(omega)
-        mag = empty((self.outputs, self.inputs, numfreq))
-        phase = empty((self.outputs, self.inputs, numfreq))
+        fresp : (self.outputs, self.inputs, len(s)) or len(s) complex ndarray
+            The frequency response of the system. Array is ``len(s)`` if and
+            only if system is SISO and ``squeeze=True``.
 
-        omega.sort()
+        Raises
+        ------
+        ValueError
+            If `s` is not purely imaginary, because
+            :class:`FrequencyDomainData` systems are only defined at imaginary
+            frequency values.
+
+        """
+        if any(abs(np.array(s, ndmin=1).real) > 0):
+            raise ValueError("__call__: FRD systems can only accept"
+                            "purely imaginary frequencies")
+        # need to preserve array or scalar status
+        if hasattr(s, '__len__'):
+            return self.eval(np.asarray(s).imag, squeeze=squeeze)
+        else:
+            return self.eval(complex(s).imag, squeeze=squeeze)
 
-        for k, w in enumerate(omega):
-            fresp = self._evalfr(w)
-            mag[:, :, k] = abs(fresp)
-            phase[:, :, k] = angle(fresp)
+    def freqresp(self, omega):
+        """(deprecated) Evaluate transfer function at complex frequencies.
 
-        return mag, phase, omega
+        .. deprecated::0.9.0
+            Method has been given the more pythonic name
+            :meth:`FrequencyResponseData.frequency_response`. Or use
+            :func:`freqresp` in the MATLAB compatibility module.
+        """
+        warn("FrequencyResponseData.freqresp(omega) will be removed in a "
+             "future release of python-control; use "
+             "FrequencyResponseData.frequency_response(omega), or "
+             "freqresp(sys, omega) in the MATLAB compatibility module "
+             "instead", DeprecationWarning)
+        return self.frequency_response(omega)
 
     def feedback(self, other=1, sign=-1):
         """Feedback interconnection between two FRD objects."""
@@ -515,11 +514,10 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1):
             "Frequency ranges of FRD do not match, conversion not implemented")
 
     elif isinstance(sys, LTI):
-        omega.sort()
-        fresp = empty((sys.outputs, sys.inputs, len(omega)), dtype=complex)
-        for k, w in enumerate(omega):
-            fresp[:, :, k] = sys._evalfr(w)
-
+        omega = np.sort(omega)
+        fresp = sys(1j * omega)
+        if len(fresp.shape) == 1:
+            fresp = fresp[np.newaxis, np.newaxis, :]
         return FRD(fresp, omega, smooth=True)
 
     elif isinstance(sys, (int, float, complex, np.number)):
diff --git a/control/freqplot.py b/control/freqplot.py
index 5a03694db..38b525aec 100644
--- a/control/freqplot.py
+++ b/control/freqplot.py
@@ -151,14 +151,14 @@ def bode_plot(syslist, omega=None,
 
     Notes
     -----
-    1. Alternatively, you may use the lower-level method
-       ``(mag, phase, freq) = sys.freqresp(freq)`` to generate the frequency
-       response for a system, but it returns a MIMO response.
+    1. Alternatively, you may use the lower-level methods
+       :meth:`LTI.frequency_response` or ``sys(s)`` or ``sys(z)`` or to
+       generate the frequency response for a single system.
 
     2. If a discrete time model is given, the frequency response is plotted
-       along the upper branch of the unit circle, using the mapping z = exp(j
-       \\omega dt) where omega ranges from 0 to pi/dt and dt is the discrete
-       timebase.  If not timebase is specified (dt = True), dt is set to 1.
+       along the upper branch of the unit circle, using the mapping ``z = exp(1j
+       * omega * dt)`` where `omega` ranges from 0 to `pi/dt` and `dt` is the discrete
+       timebase.  If timebase not specified (``dt=True``), `dt` is set to 1.
 
     Examples
     --------
@@ -228,7 +228,7 @@ def bode_plot(syslist, omega=None,
                 nyquistfrq = None
 
             # Get the magnitude and phase of the system
-            mag_tmp, phase_tmp, omega_sys = sys.freqresp(omega_sys)
+            mag_tmp, phase_tmp, omega_sys = sys.frequency_response(omega_sys)
             mag = np.atleast_1d(np.squeeze(mag_tmp))
             phase = np.atleast_1d(np.squeeze(phase_tmp))
 
@@ -588,7 +588,7 @@ def nyquist_plot(syslist, omega=None, plot=True, label_freq=0,
                 "Nyquist is currently only implemented for SISO systems.")
         else:
             # Get the magnitude and phase of the system
-            mag_tmp, phase_tmp, omega = sys.freqresp(omega)
+            mag_tmp, phase_tmp, omega = sys.frequency_response(omega)
             mag = np.squeeze(mag_tmp)
             phase = np.squeeze(phase_tmp)
 
@@ -718,7 +718,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     omega_plot = omega / (2. * math.pi) if Hz else omega
 
     # TODO: Need to add in the mag = 1 lines
-    mag_tmp, phase_tmp, omega = S.freqresp(omega)
+    mag_tmp, phase_tmp, omega = S.frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['s'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -728,7 +728,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     plot_axes['s'].tick_params(labelbottom=False)
     plot_axes['s'].grid(grid, which='both')
 
-    mag_tmp, phase_tmp, omega = (P * S).freqresp(omega)
+    mag_tmp, phase_tmp, omega = (P * S).frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['ps'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -738,7 +738,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     plot_axes['ps'].set_ylabel("$|PS|$" + " (dB)" if dB else "")
     plot_axes['ps'].grid(grid, which='both')
 
-    mag_tmp, phase_tmp, omega = (C * S).freqresp(omega)
+    mag_tmp, phase_tmp, omega = (C * S).frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['cs'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -749,7 +749,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     plot_axes['cs'].set_ylabel("$|CS|$" + " (dB)" if dB else "")
     plot_axes['cs'].grid(grid, which='both')
 
-    mag_tmp, phase_tmp, omega = T.freqresp(omega)
+    mag_tmp, phase_tmp, omega = T.frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['t'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
diff --git a/control/lti.py b/control/lti.py
index e41fe416b..eff14be2a 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -13,11 +13,11 @@
 """
 
 import numpy as np
-from numpy import absolute, real
+from numpy import absolute, real, angle, abs
 from warnings import warn
 
-__all__ = ['issiso', 'timebase', 'common_timebase', 'timebaseEqual', 
-           'isdtime', 'isctime', 'pole', 'zero', 'damp', 'evalfr', 
+__all__ = ['issiso', 'timebase', 'common_timebase', 'timebaseEqual',
+           'isdtime', 'isctime', 'pole', 'zero', 'damp', 'evalfr',
            'freqresp', 'dcgain']
 
 class LTI:
@@ -111,11 +111,56 @@ def damp(self):
         Z = -real(splane_poles)/wn
         return wn, Z, poles
 
+    def frequency_response(self, omega, squeeze=True):
+        """Evaluate the linear time-invariant system at an array of angular
+        frequencies.
+
+        Reports the frequency response of the system,
+
+             G(j*omega) = mag*exp(j*phase)
+
+        for continuous time systems. For discrete time systems, the response is
+        evaluated around the unit circle such that
+
+             G(exp(j*omega*dt)) = mag*exp(j*phase).
+
+        Parameters
+        ----------
+        omega : array_like or float
+            A list, tuple, array, or scalar value of frequencies in
+            radians/sec at which the system will be evaluated.
+        squeeze: bool, optional (default=True)
+            If True and sys is single input, single output (SISO), return a
+            1D array or scalar depending on omega's length.
+
+        Returns
+        -------
+        mag : (self.outputs, self.inputs, len(omega)) or len(omega) ndarray
+            The magnitude (absolute value, not dB or log10) of the system
+            frequency response.
+        phase : (self.outputs, self.inputs, len(omega)) or len(omega) ndarray
+            The wrapped phase in radians of the system frequency response.
+        omega : ndarray
+            The (sorted) frequencies at which the response was evaluated.
+
+        """
+        omega = np.sort(np.array(omega, ndmin=1))
+        if isdtime(self, strict=True):
+            # Convert the frequency to discrete time
+            if np.any(omega * self.dt > np.pi):
+                warn("__call__: evaluation above Nyquist frequency")
+            s = np.exp(1j * omega * self.dt)
+        else:
+            s = 1j * omega
+        response = self.__call__(s, squeeze=squeeze)
+        return abs(response), angle(response), omega
+
     def dcgain(self):
         """Return the zero-frequency gain"""
         raise NotImplementedError("dcgain not implemented for %s objects" %
                                   str(self.__class__))
 
+
 # Test to see if a system is SISO
 def issiso(sys, strict=False):
     """
@@ -162,50 +207,50 @@ def timebase(sys, strict=True):
 def common_timebase(dt1, dt2):
     """
     Find the common timebase when interconnecting systems
-    
+
     Parameters
     ----------
-    dt1, dt2: number or system with a 'dt' attribute (e.g. TransferFunction 
+    dt1, dt2: number or system with a 'dt' attribute (e.g. TransferFunction
         or StateSpace system)
 
     Returns
     -------
     dt: number
-        The common timebase of dt1 and dt2, as specified in 
-        :ref:`conventions-ref`. 
-        
+        The common timebase of dt1 and dt2, as specified in
+        :ref:`conventions-ref`.
+
     Raises
     ------
     ValueError
         when no compatible time base can be found
     """
-    # explanation: 
+    # explanation:
     # if either dt is None, they are compatible with anything
-    # if either dt is True (discrete with unspecified time base), 
+    # if either dt is True (discrete with unspecified time base),
     #   use the timebase of the other, if it is also discrete
-    # otherwise both dts must be equal 
+    # otherwise both dts must be equal
     if hasattr(dt1, 'dt'):
         dt1 = dt1.dt
     if hasattr(dt2, 'dt'):
         dt2 = dt2.dt
 
-    if dt1 is None: 
+    if dt1 is None:
         return dt2
-    elif dt2 is None: 
+    elif dt2 is None:
         return dt1
-    elif dt1 is True: 
+    elif dt1 is True:
         if dt2 > 0:
             return dt2
-        else: 
+        else:
             raise ValueError("Systems have incompatible timebases")
-    elif dt2 is True: 
-        if dt1 > 0: 
+    elif dt2 is True:
+        if dt1 > 0:
             return dt1
-        else: 
+        else:
             raise ValueError("Systems have incompatible timebases")
     elif np.isclose(dt1, dt2):
         return dt1
-    else: 
+    else:
         raise ValueError("Systems have incompatible timebases")
 
 # Check to see if two timebases are equal
@@ -221,9 +266,9 @@ def timebaseEqual(sys1, sys2):
     timebase (dt > 0) then their timebases must be equal.
     """
     warn("timebaseEqual will be deprecated in a future release of "
-         "python-control; use :func:`common_timebase` instead", 
+         "python-control; use :func:`common_timebase` instead",
          PendingDeprecationWarning)
-        
+
     if (type(sys1.dt) == bool or type(sys2.dt) == bool):
         # Make sure both are unspecified discrete timebases
         return type(sys1.dt) == type(sys2.dt) and sys1.dt == sys2.dt
@@ -413,24 +458,33 @@ def damp(sys, doprint=True):
                       (p.real, p.imag, d, w))
     return wn, damping, poles
 
-def evalfr(sys, x):
+def evalfr(sys, x, squeeze=True):
     """
-    Evaluate the transfer function of an LTI system for a single complex
-    number x.
+    Evaluate the transfer function of an LTI system for complex frequency x.
+
+    Returns the complex frequency response `sys(x)` where `x` is `s` for
+    continuous-time systems and `z` for discrete-time systems.
 
-    To evaluate at a frequency, enter x = omega*j, where omega is the
-    frequency in radians
+    To evaluate at a frequency omega in radians per second, enter
+    ``x = omega * 1j`` for continuous-time systems, or
+    ``x = exp(1j * omega * dt)`` for discrete-time systems, or use
+    ``freqresp(sys, omega)``.
 
     Parameters
     ----------
     sys: StateSpace or TransferFunction
         Linear system
-    x: scalar
-        Complex number
+    x: complex scalar or array_like
+        Complex frequency(s)
+    squeeze: bool, optional (default=True)
+        If True and sys is single input single output (SISO), returns a
+        1D array or scalar depending on the length of x.
 
     Returns
     -------
-    fresp: ndarray
+    fresp : (sys.outputs, sys.inputs, len(x)) or len(x) complex ndarray
+        The frequency response of the system. Array is len(x) if and only if
+        system is SISO and squeeze=True.
 
     See Also
     --------
@@ -439,8 +493,8 @@ def evalfr(sys, x):
 
     Notes
     -----
-    This function is a wrapper for StateSpace.evalfr and
-    TransferFunction.evalfr.
+    This function is a wrapper for StateSpace.__call__ and
+    TransferFunction.__call__.
 
     Examples
     --------
@@ -451,12 +505,9 @@ def evalfr(sys, x):
 
     .. todo:: Add example with MIMO system
     """
-    if issiso(sys):
-        return sys.horner(x)[0][0]
-    return sys.horner(x)
-
+    return sys.__call__(x, squeeze=squeeze)
 
-def freqresp(sys, omega):
+def freqresp(sys, omega, squeeze=True):
     """
     Frequency response of an LTI system at multiple angular frequencies.
 
@@ -464,19 +515,22 @@ def freqresp(sys, omega):
     ----------
     sys: StateSpace or TransferFunction
         Linear system
-    omega: array_like
+    omega: float or array_like
         A list of frequencies in radians/sec at which the system should be
         evaluated. The list can be either a python list or a numpy array
         and will be sorted before evaluation.
+    squeeze: bool, optional (default=True)
+        If True and sys is single input, single output (SISO), returns
+        1D array or scalar depending on omega's length.
 
     Returns
     -------
-    mag : (self.outputs, self.inputs, len(omega)) ndarray
+    mag : (sys.outputs, sys.inputs, len(omega)) or len(omega) ndarray
         The magnitude (absolute value, not dB or log10) of the system
         frequency response.
-    phase : (self.outputs, self.inputs, len(omega)) ndarray
+    phase : (sys.outputs, sys.inputs, len(omega)) or len(omega) ndarray
         The wrapped phase in radians of the system frequency response.
-    omega : ndarray or list or tuple
+    omega : ndarray
         The list of sorted frequencies at which the response was
         evaluated.
 
@@ -487,9 +541,8 @@ def freqresp(sys, omega):
 
     Notes
     -----
-    This function is a wrapper for StateSpace.freqresp and
-    TransferFunction.freqresp.  The output omega is a sorted version of the
-    input omega.
+    This function is a wrapper for StateSpace.frequency_response and
+    TransferFunction.frequency_response.
 
     Examples
     --------
@@ -514,7 +567,7 @@ def freqresp(sys, omega):
         #>>> # frequency response from the 1st input to the 2nd output, for
         #>>> # s = 0.1i, i, 10i.
     """
-    return sys.freqresp(omega)
+    return sys.frequency_response(omega, squeeze=squeeze)
 
 
 def dcgain(sys):
diff --git a/control/margins.py b/control/margins.py
index 03e78352f..20da2a879 100644
--- a/control/margins.py
+++ b/control/margins.py
@@ -339,24 +339,24 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
         # a bit coarse, have the interpolated frd evaluated again
         def _mod(w):
             """Calculate |G(jw)| - 1"""
-            return np.abs(sys._evalfr(w)[0][0]) - 1
+            return np.abs(sys(1j * w)) - 1
 
         def _arg(w):
             """Calculate the phase angle at -180 deg"""
-            return np.angle(-sys._evalfr(w)[0][0])
+            return np.angle(-sys(1j * w))
 
         def _dstab(w):
             """Calculate the distance from -1 point"""
-            return np.abs(sys._evalfr(w)[0][0] + 1.)
+            return np.abs(sys(1j * w) + 1.)
 
         # find the phase crossings ang(H(jw) == -180
         widx = np.where(np.diff(np.sign(_arg(sys.omega))))[0]
-        widx = widx[np.real(sys._evalfr(sys.omega[widx])[0][0]) <= 0]
+        widx = widx[np.real(sys(1j * sys.omega[widx])) <= 0]
         w_180 = np.array(
             [sp.optimize.brentq(_arg, sys.omega[i], sys.omega[i+1])
              for i in widx])
         # TODO: replace by evalfr(sys, 1J*w) or sys(1J*w), (needs gh-449)
-        w180_resp = sys._evalfr(w_180)[0][0]
+        w180_resp = sys(1j * w_180)
 
         # Find all crossings, note that this depends on omega having
         # a correct range
@@ -364,7 +364,7 @@ def _dstab(w):
         wc = np.array(
             [sp.optimize.brentq(_mod, sys.omega[i], sys.omega[i+1])
              for i in widx])
-        wc_resp = sys._evalfr(wc)[0][0]
+        wc_resp = sys(1j * wc)
 
         # find all stab margins?
         widx, = np.where(np.diff(np.sign(np.diff(_dstab(sys.omega)))) > 0)
@@ -374,7 +374,7 @@ def _dstab(w):
                                          ).x
              for i in widx])
         wstab = wstab[(wstab >= sys.omega[0]) * (wstab <= sys.omega[-1])]
-        ws_resp = sys._evalfr(wstab)[0][0]
+        ws_resp = sys(1j * wstab)
 
     with np.errstate(all='ignore'):  # |G|=0 is okay and yields inf
         GM = 1. / np.abs(w180_resp)
diff --git a/control/matlab/__init__.py b/control/matlab/__init__.py
index 413dc6d86..6b88214c6 100644
--- a/control/matlab/__init__.py
+++ b/control/matlab/__init__.py
@@ -224,8 +224,8 @@
 \*  :func:`~control.nichols`    Nichols plot
 \*  :func:`margin`              gain and phase margins
 \   lti/allmargin               all crossover frequencies and margins
-\*  :func:`freqresp`            frequency response over a frequency grid
-\*  :func:`evalfr`              frequency response at single frequency
+\*  :func:`freqresp`            frequency response
+\*  :func:`evalfr`              frequency response at complex frequency s
 ==  ==========================  ============================================
 
 
diff --git a/control/nichols.py b/control/nichols.py
index ca0505957..c1d8ff9b6 100644
--- a/control/nichols.py
+++ b/control/nichols.py
@@ -95,7 +95,7 @@ def nichols_plot(sys_list, omega=None, grid=None):
 
     for sys in sys_list:
         # Get the magnitude and phase of the system
-        mag_tmp, phase_tmp, omega = sys.freqresp(omega)
+        mag_tmp, phase_tmp, omega = sys.frequency_response(omega)
         mag = np.squeeze(mag_tmp)
         phase = np.squeeze(phase_tmp)
 
diff --git a/control/rlocus.py b/control/rlocus.py
index 479a833ab..3a3c1c2b6 100644
--- a/control/rlocus.py
+++ b/control/rlocus.py
@@ -476,7 +476,7 @@ def _systopoly1d(sys):
         sys = _convert_to_transfer_function(sys)
 
         # Make sure we have a SISO system
-        if (sys.inputs > 1 or sys.outputs > 1):
+        if not sys.issiso():
             raise ControlMIMONotImplemented()
 
         # Start by extracting the numerator and denominator from system object
@@ -497,7 +497,7 @@ def _RLFindRoots(nump, denp, kvect):
     """Find the roots for the root locus."""
     # Convert numerator and denominator to polynomials if they aren't
     roots = []
-    for k in kvect:
+    for k in np.array(kvect, ndmin=1):
         curpoly = denp + k * nump
         curroots = curpoly.r
         if len(curroots) < denp.order:
@@ -581,10 +581,10 @@ def _RLFeedbackClicksPoint(event, sys, fig, ax_rlocus, sisotool=False):
     # 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.horner(s)
-        K_xlim = -1. / sys.horner(
+        K = -1. / sys(s)
+        K_xlim = -1. / sys(
             complex(event.xdata + 0.05 * abs(xlim[1] - xlim[0]), event.ydata))
-        K_ylim = -1. / sys.horner(
+        K_ylim = -1. / sys(
             complex(event.xdata, event.ydata + 0.05 * abs(ylim[1] - ylim[0])))
 
     except TypeError:
@@ -625,7 +625,7 @@ def _RLFeedbackClicksPoint(event, sys, fig, ax_rlocus, sisotool=False):
             ax_rlocus.plot(s.real, s.imag, 'k.', marker='s', markersize=8,
                            zorder=20, label='gain_point')
 
-        return K.real[0][0]
+        return K.real
 
 
 def _removeLine(label, ax):
diff --git a/control/statesp.py b/control/statesp.py
index b41f18c7d..8bf4c0d99 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -54,7 +54,7 @@
 import math
 import numpy as np
 from numpy import any, array, asarray, concatenate, cos, delete, \
-    dot, empty, exp, eye, isinf, ones, pad, sin, zeros, squeeze
+    dot, empty, exp, eye, isinf, ones, pad, sin, zeros, squeeze, pi
 from numpy.random import rand, randn
 from numpy.linalg import solve, eigvals, matrix_rank
 from numpy.linalg.linalg import LinAlgError
@@ -640,125 +640,146 @@ def __rdiv__(self, other):
         raise NotImplementedError(
             "StateSpace.__rdiv__ is not implemented yet.")
 
-    def evalfr(self, omega):
-        """Evaluate a SS system's transfer function at a single frequency.
+    def __call__(self, x, squeeze=True):
+        """Evaluate system's transfer function at complex frequency.
 
-        self._evalfr(omega) returns the value of the transfer function matrix
-        with input value s = i * omega.
+        Returns the complex frequency response `sys(x)` where `x` is `s` for
+        continuous-time systems and `z` for discrete-time systems.
 
-        """
-        warn("StateSpace.evalfr(omega) will be deprecated in a future "
-             "release of python-control; use evalfr(sys, omega*1j) instead",
-             PendingDeprecationWarning)
-        return self._evalfr(omega)
-
-    def _evalfr(self, omega):
-        """Evaluate a SS system's transfer function at a single frequency"""
-        # Figure out the point to evaluate the transfer function
-        if isdtime(self, strict=True):
-            s = exp(1.j * omega * self.dt)
-            if omega * self.dt > math.pi:
-                warn("_evalfr: frequency evaluation above Nyquist frequency")
-        else:
-            s = omega * 1.j
+        To evaluate at a frequency omega in radians per second, enter
+        ``x = omega * 1j``, for continuous-time systems, or
+        ``x = exp(1j * omega * dt)`` for discrete-time systems. Or use
+        :meth:`StateSpace.frequency_response`.
 
-        return self.horner(s)
+        Parameters
+        ----------
+        x: complex or complex array_like
+            Complex frequencies
+        squeeze: bool, optional (default=True)
+            If True and `sys` is single input single output (SISO), returns a
+            1D array or scalar depending on the length of `x`.
 
-    def horner(self, s):
-        """Evaluate the systems's transfer function for a complex variable
+        Returns
+        -------
+        fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
+            The frequency response of the system. Array is ``len(x)`` if and
+            only if system is SISO and ``squeeze=True``.
 
-        Returns a matrix of values evaluated at complex variable s.
         """
-        resp = np.dot(self.C, solve(s * eye(self.states) - self.A,
-                                    self.B)) + self.D
-        return array(resp)
+        # Use Slycot if available
+        out = self.horner(x)
+        if not hasattr(x, '__len__'):
+            # received a scalar x, squeeze down the array along last dim
+            out = np.squeeze(out, axis=2)
+        if squeeze and self.issiso():
+            return out[0][0]
+        else:
+            return out
 
-    def freqresp(self, omega):
-        """Evaluate the system's transfer function at a list of frequencies
+    def slycot_laub(self, x):
+        """Evaluate system's transfer function at complex frequency
+        using Laub's method from Slycot.
+
+        Expects inputs and outputs to be formatted correctly. Use ``sys(x)``
+        for a more user-friendly interface.
+
+        Parameters
+        ----------
+        x : complex array_like or complex
+            Complex frequency
 
-        Reports the frequency response of the system,
+        Returns
+        -------
+        output : (number_outputs, number_inputs, len(x)) complex ndarray
+            Frequency response
+        """
+        from slycot import tb05ad
+
+        # preallocate
+        x_arr = np.atleast_1d(x) # array-like version of x
+        n = self.states
+        m = self.inputs
+        p = self.outputs
+        out = np.empty((p, m, len(x_arr)), dtype=complex)
+        # The first call both evaluates C(sI-A)^-1 B and also returns
+        # Hessenberg transformed matrices at, bt, ct.
+        result = tb05ad(n, m, p, x_arr[0], self.A, self.B, self.C, job='NG')
+        # When job='NG', result = (at, bt, ct, g_i, hinvb, info)
+        at = result[0]
+        bt = result[1]
+        ct = result[2]
+
+        # TB05AD frequency evaluation does not include direct feedthrough.
+        out[:, :, 0] = result[3] + self.D
+
+        # Now, iterate through the remaining frequencies using the
+        # transformed state matrices, at, bt, ct.
+
+        # Start at the second frequency, already have the first.
+        for kk, x_kk in enumerate(x_arr[1:len(x_arr)]):
+            result = tb05ad(n, m, p, x_kk, at, bt, ct, job='NH')
+            # When job='NH', result = (g_i, hinvb, info)
+
+            # kk+1 because enumerate starts at kk = 0.
+            # but zero-th spot is already filled.
+            out[:, :, kk+1] = result[0] + self.D
+        return out
 
-             G(j*omega) = mag*exp(j*phase)
+    def horner(self, x):
+        """Evaluate system's transfer function at complex frequency
+        using Laub's or Horner's method.
 
-        for continuous time. For discrete time systems, the response is
-        evaluated around the unit circle such that
+        Evaluates `sys(x)` where `x` is `s` for continuous-time systems and `z`
+        for discrete-time systems.
 
-             G(exp(j*omega*dt)) = mag*exp(j*phase).
+        Expects inputs and outputs to be formatted correctly. Use ``sys(x)``
+        for a more user-friendly interface.
 
         Parameters
         ----------
-        omega : array_like
-            A list of frequencies in radians/sec at which the system should be
-            evaluated. The list can be either a python list or a numpy array
-            and will be sorted before evaluation.
+        x : complex array_like or complex
+            Complex frequencies
 
         Returns
         -------
-        mag : (self.outputs, self.inputs, len(omega)) ndarray
-            The magnitude (absolute value, not dB or log10) of the system
-            frequency response.
-        phase : (self.outputs, self.inputs, len(omega)) ndarray
-            The wrapped phase in radians of the system frequency response.
-        omega : ndarray
-            The list of sorted frequencies at which the response was
-            evaluated.
-        """
-        # In case omega is passed in as a list, rather than a proper array.
-        omega = np.asarray(omega)
-
-        numFreqs = len(omega)
-        Gfrf = np.empty((self.outputs, self.inputs, numFreqs),
-                        dtype=np.complex128)
-
-        # Sort frequency and calculate complex frequencies on either imaginary
-        # axis (continuous time) or unit circle (discrete time).
-        omega.sort()
-        if isdtime(self, strict=True):
-            cmplx_freqs = exp(1.j * omega * self.dt)
-            if max(np.abs(omega)) * self.dt > math.pi:
-                warn("freqresp: frequency evaluation above Nyquist frequency")
-        else:
-            cmplx_freqs = omega * 1.j
+        output : (self.outputs, self.inputs, len(x)) complex ndarray
+            Frequency response
 
-        # Do the frequency response evaluation. Use TB05AD from Slycot
-        # if it's available, otherwise use the built-in horners function.
+        Notes
+        -----
+        Attempts to use Laub's method from Slycot library, with a
+        fall-back to python code.
+        """
         try:
-            from slycot import tb05ad
-
-            n = np.shape(self.A)[0]
-            m = self.inputs
-            p = self.outputs
-            # The first call both evaluates C(sI-A)^-1 B and also returns
-            # Hessenberg transformed matrices at, bt, ct.
-            result = tb05ad(n, m, p, cmplx_freqs[0], self.A,
-                            self.B, self.C, job='NG')
-            # When job='NG', result = (at, bt, ct, g_i, hinvb, info)
-            at = result[0]
-            bt = result[1]
-            ct = result[2]
-
-            # TB05AD frequency evaluation does not include direct feedthrough.
-            Gfrf[:, :, 0] = result[3] + self.D
-
-            # Now, iterate through the remaining frequencies using the
-            # transformed state matrices, at, bt, ct.
-
-            # Start at the second frequency, already have the first.
-            for kk, cmplx_freqs_kk in enumerate(cmplx_freqs[1:numFreqs]):
-                result = tb05ad(n, m, p, cmplx_freqs_kk, at,
-                                bt, ct, job='NH')
-                # When job='NH', result = (g_i, hinvb, info)
-
-                # kk+1 because enumerate starts at kk = 0.
-                # but zero-th spot is already filled.
-                Gfrf[:, :, kk+1] = result[0] + self.D
-
-        except ImportError:  # Slycot unavailable. Fall back to horner.
-            for kk, cmplx_freqs_kk in enumerate(cmplx_freqs):
-                Gfrf[:, :, kk] = self.horner(cmplx_freqs_kk)
-
-        #      mag           phase           omega
-        return np.abs(Gfrf), np.angle(Gfrf), omega
+            out = self.slycot_laub(x)
+        except (ImportError, Exception):
+            # Fall back because either Slycot unavailable or cannot handle
+            # certain cases.
+            x_arr = np.atleast_1d(x) # force to be an array
+            # Preallocate
+            out = empty((self.outputs, self.inputs, len(x_arr)), dtype=complex)
+
+            #TODO: can this be vectorized?
+            for idx, x_idx in enumerate(x_arr):
+                out[:,:,idx] = \
+                    np.dot(self.C,
+                        solve(x_idx * eye(self.states) - self.A, self.B)) \
+                    + self.D
+        return out
+
+    def freqresp(self, omega):
+        """(deprecated) Evaluate transfer function at complex frequencies.
+
+        .. deprecated::0.9.0
+            Method has been given the more pythonic name
+            :meth:`StateSpace.frequency_response`. Or use
+            :func:`freqresp` in the MATLAB compatibility module.
+        """
+        warn("StateSpace.freqresp(omega) will be removed in a "
+             "future release of python-control; use "
+             "sys.frequency_response(omega), or freqresp(sys, omega) in the "
+             "MATLAB compatibility module instead", DeprecationWarning)
+        return self.frequency_response(omega)
 
     # Compute poles and zeros
     def pole(self):
@@ -1148,7 +1169,7 @@ def dcgain(self):
                 gain = np.asarray(self.D -
                                   self.C.dot(np.linalg.solve(self.A, self.B)))
             else:
-                gain = self.horner(1)
+                gain = np.squeeze(self.horner(1))
         except LinAlgError:
             # eigenvalue at DC
             gain = np.tile(np.nan, (self.outputs, self.inputs))
diff --git a/control/tests/frd_test.py b/control/tests/frd_test.py
index 18f2f17b1..5343112fe 100644
--- a/control/tests/frd_test.py
+++ b/control/tests/frd_test.py
@@ -39,8 +39,8 @@ def testSISOtf(self):
         frd = FRD(h, omega)
         assert isinstance(frd, FRD)
 
-        mag1, phase1, omega1 = frd.freqresp([1.0])
-        mag2, phase2, omega2 = h.freqresp([1.0])
+        mag1, phase1, omega1 = frd.frequency_response([1.0])
+        mag2, phase2, omega2 = h.frequency_response([1.0])
         np.testing.assert_array_almost_equal(mag1, mag2)
         np.testing.assert_array_almost_equal(phase1, phase2)
         np.testing.assert_array_almost_equal(omega1, omega2)
@@ -54,39 +54,39 @@ def testOperators(self):
         f2 = FRD(h2, omega)
 
         np.testing.assert_array_almost_equal(
-            (f1 + f2).freqresp([0.1, 1.0, 10])[0],
-            (h1 + h2).freqresp([0.1, 1.0, 10])[0])
+            (f1 + f2).frequency_response([0.1, 1.0, 10])[0],
+            (h1 + h2).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1 + f2).freqresp([0.1, 1.0, 10])[1],
-            (h1 + h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 + f2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 + h2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1 - f2).freqresp([0.1, 1.0, 10])[0],
-            (h1 - h2).freqresp([0.1, 1.0, 10])[0])
+            (f1 - f2).frequency_response([0.1, 1.0, 10])[0],
+            (h1 - h2).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1 - f2).freqresp([0.1, 1.0, 10])[1],
-            (h1 - h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 - f2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 - h2).frequency_response([0.1, 1.0, 10])[1])
 
         # multiplication and division
         np.testing.assert_array_almost_equal(
-            (f1 * f2).freqresp([0.1, 1.0, 10])[1],
-            (h1 * h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 * f2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 * h2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1 / f2).freqresp([0.1, 1.0, 10])[1],
-            (h1 / h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 / f2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 / h2).frequency_response([0.1, 1.0, 10])[1])
 
         # with default conversion from scalar
         np.testing.assert_array_almost_equal(
-            (f1 * 1.5).freqresp([0.1, 1.0, 10])[1],
-            (h1 * 1.5).freqresp([0.1, 1.0, 10])[1])
+            (f1 * 1.5).frequency_response([0.1, 1.0, 10])[1],
+            (h1 * 1.5).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1 / 1.7).freqresp([0.1, 1.0, 10])[1],
-            (h1 / 1.7).freqresp([0.1, 1.0, 10])[1])
+            (f1 / 1.7).frequency_response([0.1, 1.0, 10])[1],
+            (h1 / 1.7).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (2.2 * f2).freqresp([0.1, 1.0, 10])[1],
-            (2.2 * h2).freqresp([0.1, 1.0, 10])[1])
+            (2.2 * f2).frequency_response([0.1, 1.0, 10])[1],
+            (2.2 * h2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (1.3 / f2).freqresp([0.1, 1.0, 10])[1],
-            (1.3 / h2).freqresp([0.1, 1.0, 10])[1])
+            (1.3 / f2).frequency_response([0.1, 1.0, 10])[1],
+            (1.3 / h2).frequency_response([0.1, 1.0, 10])[1])
 
     def testOperatorsTf(self):
         # get two SISO transfer functions
@@ -98,24 +98,24 @@ def testOperatorsTf(self):
         f2  # reference to avoid pyflakes error
 
         np.testing.assert_array_almost_equal(
-            (f1 + h2).freqresp([0.1, 1.0, 10])[0],
-            (h1 + h2).freqresp([0.1, 1.0, 10])[0])
+            (f1 + h2).frequency_response([0.1, 1.0, 10])[0],
+            (h1 + h2).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1 + h2).freqresp([0.1, 1.0, 10])[1],
-            (h1 + h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 + h2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 + h2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1 - h2).freqresp([0.1, 1.0, 10])[0],
-            (h1 - h2).freqresp([0.1, 1.0, 10])[0])
+            (f1 - h2).frequency_response([0.1, 1.0, 10])[0],
+            (h1 - h2).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1 - h2).freqresp([0.1, 1.0, 10])[1],
-            (h1 - h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 - h2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 - h2).frequency_response([0.1, 1.0, 10])[1])
         # multiplication and division
         np.testing.assert_array_almost_equal(
-            (f1 * h2).freqresp([0.1, 1.0, 10])[1],
-            (h1 * h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 * h2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 * h2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1 / h2).freqresp([0.1, 1.0, 10])[1],
-            (h1 / h2).freqresp([0.1, 1.0, 10])[1])
+            (f1 / h2).frequency_response([0.1, 1.0, 10])[1],
+            (h1 / h2).frequency_response([0.1, 1.0, 10])[1])
         # the reverse does not work
 
     def testbdalg(self):
@@ -127,45 +127,45 @@ def testbdalg(self):
         f2 = FRD(h2, omega)
 
         np.testing.assert_array_almost_equal(
-            (bdalg.series(f1, f2)).freqresp([0.1, 1.0, 10])[0],
-            (bdalg.series(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+            (bdalg.series(f1, f2)).frequency_response([0.1, 1.0, 10])[0],
+            (bdalg.series(h1, h2)).frequency_response([0.1, 1.0, 10])[0])
 
         np.testing.assert_array_almost_equal(
-            (bdalg.parallel(f1, f2)).freqresp([0.1, 1.0, 10])[0],
-            (bdalg.parallel(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+            (bdalg.parallel(f1, f2)).frequency_response([0.1, 1.0, 10])[0],
+            (bdalg.parallel(h1, h2)).frequency_response([0.1, 1.0, 10])[0])
 
         np.testing.assert_array_almost_equal(
-            (bdalg.feedback(f1, f2)).freqresp([0.1, 1.0, 10])[0],
-            (bdalg.feedback(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+            (bdalg.feedback(f1, f2)).frequency_response([0.1, 1.0, 10])[0],
+            (bdalg.feedback(h1, h2)).frequency_response([0.1, 1.0, 10])[0])
 
         np.testing.assert_array_almost_equal(
-            (bdalg.negate(f1)).freqresp([0.1, 1.0, 10])[0],
-            (bdalg.negate(h1)).freqresp([0.1, 1.0, 10])[0])
+            (bdalg.negate(f1)).frequency_response([0.1, 1.0, 10])[0],
+            (bdalg.negate(h1)).frequency_response([0.1, 1.0, 10])[0])
 
 #       append() and connect() not implemented for FRD objects
 #        np.testing.assert_array_almost_equal(
-#            (bdalg.append(f1, f2)).freqresp([0.1, 1.0, 10])[0],
-#            (bdalg.append(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+#            (bdalg.append(f1, f2)).frequency_response([0.1, 1.0, 10])[0],
+#            (bdalg.append(h1, h2)).frequency_response([0.1, 1.0, 10])[0])
 #
 #        f3 = bdalg.append(f1, f2, f2)
 #        h3 = bdalg.append(h1, h2, h2)
 #        Q = np.mat([ [1, 2], [2, -1] ])
 #        np.testing.assert_array_almost_equal(
-#           (bdalg.connect(f3, Q, [2], [1])).freqresp([0.1, 1.0, 10])[0],
-#            (bdalg.connect(h3, Q, [2], [1])).freqresp([0.1, 1.0, 10])[0])
+#           (bdalg.connect(f3, Q, [2], [1])).frequency_response([0.1, 1.0, 10])[0],
+#            (bdalg.connect(h3, Q, [2], [1])).frequency_response([0.1, 1.0, 10])[0])
 
     def testFeedback(self):
         h1 = TransferFunction([1], [1, 2, 2])
         omega = np.logspace(-1, 2, 10)
         f1 = FRD(h1, omega)
         np.testing.assert_array_almost_equal(
-            f1.feedback(1).freqresp([0.1, 1.0, 10])[0],
-            h1.feedback(1).freqresp([0.1, 1.0, 10])[0])
+            f1.feedback(1).frequency_response([0.1, 1.0, 10])[0],
+            h1.feedback(1).frequency_response([0.1, 1.0, 10])[0])
 
         # Make sure default argument also works
         np.testing.assert_array_almost_equal(
-            f1.feedback().freqresp([0.1, 1.0, 10])[0],
-            h1.feedback().freqresp([0.1, 1.0, 10])[0])
+            f1.feedback().frequency_response([0.1, 1.0, 10])[0],
+            h1.feedback().frequency_response([0.1, 1.0, 10])[0])
 
     def testFeedback2(self):
         h2 = StateSpace([[-1.0, 0], [0, -2.0]], [[0.4], [0.1]],
@@ -198,11 +198,11 @@ def testMIMO(self):
         omega = np.logspace(-1, 2, 10)
         f1 = FRD(sys, omega)
         np.testing.assert_array_almost_equal(
-            sys.freqresp([0.1, 1.0, 10])[0],
-            f1.freqresp([0.1, 1.0, 10])[0])
+            sys.frequency_response([0.1, 1.0, 10])[0],
+            f1.frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            sys.freqresp([0.1, 1.0, 10])[1],
-            f1.freqresp([0.1, 1.0, 10])[1])
+            sys.frequency_response([0.1, 1.0, 10])[1],
+            f1.frequency_response([0.1, 1.0, 10])[1])
 
     @slycotonly
     def testMIMOfb(self):
@@ -214,11 +214,11 @@ def testMIMOfb(self):
         f1 = FRD(sys, omega).feedback([[0.1, 0.3], [0.0, 1.0]])
         f2 = FRD(sys.feedback([[0.1, 0.3], [0.0, 1.0]]), omega)
         np.testing.assert_array_almost_equal(
-            f1.freqresp([0.1, 1.0, 10])[0],
-            f2.freqresp([0.1, 1.0, 10])[0])
+            f1.frequency_response([0.1, 1.0, 10])[0],
+            f2.frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            f1.freqresp([0.1, 1.0, 10])[1],
-            f2.freqresp([0.1, 1.0, 10])[1])
+            f1.frequency_response([0.1, 1.0, 10])[1],
+            f2.frequency_response([0.1, 1.0, 10])[1])
 
     @slycotonly
     def testMIMOfb2(self):
@@ -232,11 +232,11 @@ def testMIMOfb2(self):
         f1 = FRD(sys, omega).feedback(K)
         f2 = FRD(sys.feedback(K), omega)
         np.testing.assert_array_almost_equal(
-            f1.freqresp([0.1, 1.0, 10])[0],
-            f2.freqresp([0.1, 1.0, 10])[0])
+            f1.frequency_response([0.1, 1.0, 10])[0],
+            f2.frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            f1.freqresp([0.1, 1.0, 10])[1],
-            f2.freqresp([0.1, 1.0, 10])[1])
+            f1.frequency_response([0.1, 1.0, 10])[1],
+            f2.frequency_response([0.1, 1.0, 10])[1])
 
     @slycotonly
     def testMIMOMult(self):
@@ -248,11 +248,11 @@ def testMIMOMult(self):
         f1 = FRD(sys, omega)
         f2 = FRD(sys, omega)
         np.testing.assert_array_almost_equal(
-            (f1*f2).freqresp([0.1, 1.0, 10])[0],
-            (sys*sys).freqresp([0.1, 1.0, 10])[0])
+            (f1*f2).frequency_response([0.1, 1.0, 10])[0],
+            (sys*sys).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1*f2).freqresp([0.1, 1.0, 10])[1],
-            (sys*sys).freqresp([0.1, 1.0, 10])[1])
+            (f1*f2).frequency_response([0.1, 1.0, 10])[1],
+            (sys*sys).frequency_response([0.1, 1.0, 10])[1])
 
     @slycotonly
     def testMIMOSmooth(self):
@@ -265,14 +265,14 @@ def testMIMOSmooth(self):
         f1 = FRD(sys, omega, smooth=True)
         f2 = FRD(sys2, omega, smooth=True)
         np.testing.assert_array_almost_equal(
-            (f1*f2).freqresp([0.1, 1.0, 10])[0],
-            (sys*sys2).freqresp([0.1, 1.0, 10])[0])
+            (f1*f2).frequency_response([0.1, 1.0, 10])[0],
+            (sys*sys2).frequency_response([0.1, 1.0, 10])[0])
         np.testing.assert_array_almost_equal(
-            (f1*f2).freqresp([0.1, 1.0, 10])[1],
-            (sys*sys2).freqresp([0.1, 1.0, 10])[1])
+            (f1*f2).frequency_response([0.1, 1.0, 10])[1],
+            (sys*sys2).frequency_response([0.1, 1.0, 10])[1])
         np.testing.assert_array_almost_equal(
-            (f1*f2).freqresp([0.1, 1.0, 10])[2],
-            (sys*sys2).freqresp([0.1, 1.0, 10])[2])
+            (f1*f2).frequency_response([0.1, 1.0, 10])[2],
+            (sys*sys2).frequency_response([0.1, 1.0, 10])[2])
 
     def testAgainstOctave(self):
         # with data from octave:
@@ -285,8 +285,8 @@ def testAgainstOctave(self):
         omega = np.logspace(-1, 2, 10)
         f1 = FRD(sys, omega)
         np.testing.assert_array_almost_equal(
-            (f1.freqresp([1.0])[0] *
-             np.exp(1j * f1.freqresp([1.0])[1])).reshape(3, 2),
+            (f1.frequency_response([1.0])[0] *
+             np.exp(1j * f1.frequency_response([1.0])[1])).reshape(3, 2),
             np.array([[0.4 - 0.2j, 0], [0, 0.1 - 0.2j], [0, 0.3 - 0.1j]]))
 
     def test_string_representation(self, capsys):
@@ -407,13 +407,9 @@ def test_eval(self):
         with pytest.raises(ValueError):
             frd_tf.eval(2)
 
-
-    def test_evalfr_deprecated(self):
-        sys_tf = ct.tf([1], [1, 2, 1])
-        frd_tf = FRD(sys_tf, np.logspace(-1, 1, 3))
-
-        with pytest.deprecated_call():
-            frd_tf.evalfr(1.)
+        # Should get an error if we use __call__ at real-valued frequency
+        with pytest.raises(ValueError):
+            frd_tf(2)
 
     def test_repr_str(self):
         # repr printing
diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py
index 111d4296c..e6a6934e8 100644
--- a/control/tests/freqresp_test.py
+++ b/control/tests/freqresp_test.py
@@ -18,7 +18,6 @@
 from control.matlab import ss, tf, bode, rss
 from control.tests.conftest import slycotonly
 
-
 pytestmark = pytest.mark.usefixtures("mplcleanup")
 
 
@@ -32,7 +31,7 @@ def test_siso():
     omega = np.linspace(10e-2, 10e2, 1000)
 
     # test frequency response
-    sys.freqresp(omega)
+    sys.frequency_response(omega)
 
     # test bode plot
     bode(sys)
@@ -97,7 +96,7 @@ def test_superimpose():
 def test_doubleint():
     """Test typcast bug with double int
 
-    30 May 2016, RMM: added to replicate typecast bug in freqresp.py
+    30 May 2016, RMM: added to replicate typecast bug in frequency_response.py
     """
     A = np.array([[0, 1], [0, 0]])
     B = np.array([[0], [1]])
@@ -117,7 +116,7 @@ def test_mimo():
     omega = np.linspace(10e-2, 10e2, 1000)
     sysMIMO = ss(A, B, C, D)
 
-    sysMIMO.freqresp(omega)
+    sysMIMO.frequency_response(omega)
     tf(sysMIMO)
 
 
@@ -215,13 +214,13 @@ def test_discrete(dsystem_type):
     omega_ok = np.linspace(10e-4, 0.99, 100) * np.pi / dsys.dt
 
     # Test frequency response
-    dsys.freqresp(omega_ok)
+    dsys.frequency_response(omega_ok)
 
     # Check for warning if frequency is out of range
     with pytest.warns(UserWarning, match="above.*Nyquist"):
         # Look for a warning about sampling above Nyquist frequency
         omega_bad = np.linspace(10e-4, 1.1, 10) * np.pi / dsys.dt
-        dsys.freqresp(omega_bad)
+        dsys.frequency_response(omega_bad)
 
     # Test bode plots (currently only implemented for SISO)
     if (dsys.inputs == 1 and dsys.outputs == 1):
@@ -332,6 +331,7 @@ def test_initial_phase(TF, initial_phase, default_phase, expected_phase):
      pytest.param(ctrl.tf([1], [1, 0, 0, 0, 0, 0]),
                   -270, -3*math.pi/2, math.pi/2,    id="order5, -270"),
     ])
+
 def test_phase_wrap(TF, wrap_phase, min_phase, max_phase):
     mag, phase, omega = ctrl.bode(TF, wrap_phase=wrap_phase)
     assert(min(phase) >= min_phase)
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index 4a8e09930..0dcbf3325 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -1132,8 +1132,8 @@ def test_lineariosys_statespace(self, tsys):
         np.testing.assert_array_equal(
             iosys_siso.pole(), tsys.siso_linsys.pole())
         omega = np.logspace(.1, 10, 100)
-        mag_io, phase_io, omega_io = iosys_siso.freqresp(omega)
-        mag_ss, phase_ss, omega_ss = tsys.siso_linsys.freqresp(omega)
+        mag_io, phase_io, omega_io = iosys_siso.frequency_response(omega)
+        mag_ss, phase_ss, omega_ss = tsys.siso_linsys.frequency_response(omega)
         np.testing.assert_array_equal(mag_io, mag_ss)
         np.testing.assert_array_equal(phase_io, phase_ss)
         np.testing.assert_array_equal(omega_io, omega_ss)
diff --git a/control/tests/margin_test.py b/control/tests/margin_test.py
index 80916da1b..81a4b68b4 100644
--- a/control/tests/margin_test.py
+++ b/control/tests/margin_test.py
@@ -93,7 +93,7 @@ def test_stability_margins_3input(tsys):
     sys, refout, refoutall = tsys
     """Test stability_margins() function with mag, phase, omega input"""
     omega = np.logspace(-2, 2, 2000)
-    mag, phase, omega_ = sys.freqresp(omega)
+    mag, phase, omega_ = sys.frequency_response(omega)
     out = stability_margins((mag, phase*180/np.pi, omega_))
     assert_allclose(out, refout, atol=1.5e-3)
 
@@ -109,7 +109,7 @@ def test_margin_3input(tsys):
     sys, refout, refoutall = tsys
     """Test margin() function with mag, phase, omega input"""
     omega = np.logspace(-2, 2, 2000)
-    mag, phase, omega_ = sys.freqresp(omega)
+    mag, phase, omega_ = sys.frequency_response(omega)
     out = margin((mag, phase*180/np.pi, omega_))
     assert_allclose(out, np.array(refout)[[0, 1, 3, 4]], atol=1.5e-3)
 
@@ -136,8 +136,10 @@ def test_phase_crossover_frequencies_mimo():
                            [[3], [4]]],
                           [[[1, 2, 3, 4], [1, 1]],
                            [[1, 1], [1, 1]]])
-    with pytest.raises(ControlMIMONotImplemented):
-        omega, gain = phase_crossover_frequencies(tf)
+
+    omega, gain = phase_crossover_frequencies(tf[0,0])
+    assert_allclose(omega, [1.73205,  0.], atol=1.5e-3)
+    assert_allclose(gain, [-0.5,  0.25], atol=1.5e-3)
 
 
 def test_mag_phase_omega():
@@ -145,7 +147,7 @@ def test_mag_phase_omega():
     sys = TransferFunction(15, [1, 6, 11, 6])
     out = stability_margins(sys)
     omega = np.logspace(-2, 2, 1000)
-    mag, phase, omega = sys.freqresp(omega)
+    mag, phase, omega = sys.frequency_response(omega)
     out2 = stability_margins((mag, phase*180/np.pi, omega))
     ind = [0, 1, 3, 4]   # indices of gm, pm, wg, wp -- ignore sm
     marg1 = np.array(out)[ind]
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 02fac1776..7bf90f4c6 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -341,18 +341,14 @@ def test_evalfr(self, dt, omega, resp):
         else:
             s = omega * 1j
 
-        # Correct version of the call
+        # Correct versions of the call
         np.testing.assert_allclose(evalfr(sys, s), resp, atol=1e-3)
-        # Deprecated version of the call (should generate warning)
-        with pytest.deprecated_call():
+        np.testing.assert_allclose(sys(s), resp, atol=1e-3)
+
+        # Deprecated version of the call (should generate error)
+        with pytest.raises(AttributeError):
             np.testing.assert_allclose(sys.evalfr(omega), resp, atol=1e-3)
 
-        # call above nyquist frequency
-        if dt:
-            with pytest.warns(UserWarning):
-                np.testing.assert_allclose(sys._evalfr(omega + 2 * np.pi / dt),
-                                           resp,
-                                           atol=1e-3)
 
     @slycotonly
     def test_freq_resp(self):
@@ -374,7 +370,7 @@ def test_freq_resp(self):
                       [-0.438157380501337, -1.40720969147217]]]
         true_omega = [0.1, 10.]
 
-        mag, phase, omega = sys.freqresp(true_omega)
+        mag, phase, omega = sys.frequency_response(true_omega)
 
         np.testing.assert_almost_equal(mag, true_mag)
         np.testing.assert_almost_equal(phase, true_phase)
@@ -739,9 +735,9 @@ def test_horner(self, sys322):
         sys322.horner(1. + 1.j)
 
         # Make sure result agrees with frequency response
-        mag, phase, omega = sys322.freqresp([1])
+        mag, phase, omega = sys322.frequency_response([1])
         np.testing.assert_array_almost_equal(
-            sys322.horner(1.j),
+            np.squeeze(sys322.horner(1.j)),
             mag[:, :, 0] * np.exp(1.j * phase[:, :, 0]))
 
 class TestRss:
@@ -902,7 +898,6 @@ def test_statespace_defaults(self, matarrayout):
 refkey_n = {None: 'p3', '.3g': 'p3', '.5g': 'p5'}
 refkey_r = {None: 'p', 'partitioned': 'p', 'separate': 's'}
 
-
 @pytest.mark.parametrize(" gmats,  ref",
                          [(LTX_G1, LTX_G1_REF),
                           (LTX_G2, LTX_G2_REF)])
@@ -928,3 +923,27 @@ def test_latex_repr(gmats, ref, repr_type, num_format, editsdefaults):
     refkey = "{}_{}".format(refkey_n[num_format], refkey_r[repr_type])
     assert g._repr_latex_() == ref[refkey]
 
+    def test_pole_static(self):
+        """Regression: pole() of static gain is empty array."""
+        np.testing.assert_array_equal(np.array([]),
+                                      StateSpace([], [], [], [[1]]).pole())
+
+    def test_copy_constructor(self):
+        # Create a set of matrices for a simple linear system
+        A = np.array([[-1]])
+        B = np.array([[1]])
+        C = np.array([[1]])
+        D = np.array([[0]])
+
+        # Create the first linear system and a copy
+        linsys = StateSpace(A, B, C, D)
+        cpysys = StateSpace(linsys)
+
+        # Change the original A matrix
+        A[0, 0] = -2
+        np.testing.assert_array_equal(linsys.A, [[-1]]) # original value
+        np.testing.assert_array_equal(cpysys.A, [[-1]]) # original value
+
+        # Change the A matrix for the original system
+        linsys.A[0, 0] = -3
+        np.testing.assert_array_equal(cpysys.A, [[-1]]) # original value
diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py
index c0b5e227f..1ec596a6e 100644
--- a/control/tests/xferfcn_test.py
+++ b/control/tests/xferfcn_test.py
@@ -404,34 +404,6 @@ def test_slice(self):
         assert (sys1.inputs, sys1.outputs) == (2, 1)
         assert sys1.dt == 0.5
 
-    @pytest.mark.parametrize("omega, resp",
-                             [(1, np.array([[-0.5 - 0.5j]])),
-                              (32, np.array([[0.002819593 - 0.03062847j]]))])
-    @pytest.mark.parametrize("dt", [None, 0, 1e-3])
-    def test_evalfr_siso(self, dt, omega, resp):
-        """Evaluate the frequency response at single frequencies"""
-        sys = TransferFunction([1., 3., 5], [1., 6., 2., -1])
-
-        if dt:
-            sys = sample_system(sys, dt)
-            s = np.exp(omega * 1j * dt)
-        else:
-            s = omega * 1j
-
-        # Correct versions of the call
-        np.testing.assert_allclose(evalfr(sys, s), resp, atol=1e-3)
-        np.testing.assert_allclose(sys(s), resp, atol=1e-3)
-        # Deprecated version of the call (should generate warning)
-        with pytest.deprecated_call():
-            np.testing.assert_allclose(sys.evalfr(omega), resp, atol=1e-3)
-
-        # call above nyquist frequency
-        if dt:
-            with pytest.warns(UserWarning):
-                np.testing.assert_allclose(sys._evalfr(omega + 2 * np.pi / dt),
-                                           resp,
-                                           atol=1e-3)
-
     def test_is_static_gain(self):
         numstatic = 1.1
         denstatic = 1.2
@@ -454,10 +426,37 @@ def test_is_static_gain(self):
         assert not TransferFunction(numdynamicmimo,
                                     denstaticmimo).is_static_gain()
 
+    @pytest.mark.parametrize("omega, resp",
+                             [(1, np.array([[-0.5 - 0.5j]])),
+                              (32, np.array([[0.002819593 - 0.03062847j]]))])
+    @pytest.mark.parametrize("dt", [None, 0, 1e-3])
+    def test_call_siso(self, dt, omega, resp):
+        """Evaluate the frequency response of a SISO system at one frequency."""
+        sys = TransferFunction([1., 3., 5], [1., 6., 2., -1])
+
+        if dt:
+            sys = sample_system(sys, dt)
+            s = np.exp(omega * 1j * dt)
+        else:
+            s = omega * 1j
+
+        # Correct versions of the call
+        np.testing.assert_allclose(evalfr(sys, s), resp, atol=1e-3)
+        np.testing.assert_allclose(sys(s), resp, atol=1e-3)
+        # Deprecated version of the call (should generate exception)
+        with pytest.raises(AttributeError):
+            np.testing.assert_allclose(sys.evalfr(omega), resp, atol=1e-3)
+
+
+    @nopython2
+    def test_call_dtime(self):
+        sys = TransferFunction([1., 3., 5], [1., 6., 2., -1], 0.1)
+        np.testing.assert_array_almost_equal(sys(1j), -0.5 - 0.5j)
 
     @slycotonly
-    def test_evalfr_mimo(self):
-        """Evaluate the frequency response of a MIMO system at a freq"""
+    def test_call_mimo(self):
+        """Evaluate the frequency response of a MIMO system at one frequency."""
+
         num = [[[1., 2.], [0., 3.], [2., -1.]],
                [[1.], [4., 0.], [1., -4., 3.]]]
         den = [[[-3., 2., 4.], [1., 0., 0.], [2., -1.]],
@@ -467,13 +466,28 @@ def test_evalfr_mimo(self):
                 [-0.083333333333333, -0.188235294117647 - 0.847058823529412j,
                  -1. - 8.j]]
 
-        np.testing.assert_array_almost_equal(sys._evalfr(2.), resp)
+        np.testing.assert_array_almost_equal(evalfr(sys, 2j), resp)
 
         # Test call version as well
         np.testing.assert_array_almost_equal(sys(2.j), resp)
 
-    def test_freqresp_siso(self):
-        """Evaluate the SISO magnitude and phase at multiple frequencies"""
+    def test_freqresp_deprecated(self):
+        sys = TransferFunction([1., 3., 5], [1., 6., 2., -1.])
+        # Deprecated version of the call (should generate warning)
+        import warnings
+        with warnings.catch_warnings(record=True) as w:
+            # Set up warnings filter to only show warnings in control module
+            warnings.filterwarnings("ignore")
+            warnings.filterwarnings("always", module="control")
+
+            # Make sure that we get a pending deprecation warning
+            sys.freqresp(1.)
+            assert issubclass(w[-1].category, DeprecationWarning)
+
+    def test_frequency_response_siso(self):
+        """Evaluate the magnitude and phase of a SISO system at
+        multiple frequencies."""
+
         sys = TransferFunction([1., 3., 5], [1., 6., 2., -1])
 
         truemag = [[[4.63507337473906, 0.707106781186548, 0.0866592803995351]]]
@@ -481,7 +495,7 @@ def test_freqresp_siso(self):
                        -1.32655885133871]]]
         trueomega = [0.1, 1., 10.]
 
-        mag, phase, omega = sys.freqresp(trueomega)
+        mag, phase, omega = sys.frequency_response(trueomega, squeeze=False)
 
         np.testing.assert_array_almost_equal(mag, truemag)
         np.testing.assert_array_almost_equal(phase, truephase)
@@ -509,7 +523,7 @@ def test_freqresp_mimo(self):
                        [-1.66852323, -1.89254688, -1.62050658],
                        [-0.13298964, -1.10714871, -2.75046720]]]
 
-        mag, phase, omega = sys.freqresp(true_omega)
+        mag, phase, omega = sys.frequency_response(true_omega)
 
         np.testing.assert_array_almost_equal(mag, true_mag)
         np.testing.assert_array_almost_equal(phase, true_phase)
diff --git a/control/xferfcn.py b/control/xferfcn.py
index 9a70e36b6..a0941b74a 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -231,19 +231,70 @@ def __init__(self, *args, **kwargs):
                     dt = config.defaults['control.default_dt']
         self.dt = dt
 
-    def __call__(self, s):
-        """Evaluate the system's transfer function for a complex variable
+    def __call__(self, x, squeeze=True):
+        """Evaluate system's transfer function at complex frequencies.
 
-        For a SISO transfer function, returns the value of the
-        transfer function.  For a MIMO transfer fuction, returns a
-        matrix of values evaluated at complex variable s."""
+        Returns the complex frequency response `sys(x)` where `x` is `s` for
+        continuous-time systems and `z` for discrete-time systems.
 
-        if self.issiso():
-            # return a scalar
-            return self.horner(s)[0][0]
+        To evaluate at a frequency omega in radians per second, enter
+        ``x = omega * 1j``, for continuous-time systems, or
+        ``x = exp(1j * omega * dt)`` for discrete-time systems. Or use
+        :meth:`TransferFunction.frequency_response`.
+
+        Parameters
+        ----------
+        x: complex array_like or complex
+            Complex frequencies
+        squeeze: bool, optional (default=True)
+            If True and `sys` is single input single output (SISO), returns a
+            1D array or scalar depending on the length of `x`.
+
+        Returns
+        -------
+        fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
+            The frequency response of the system. Array is `len(x)` if and
+            only if system is SISO and ``squeeze=True``.
+
+        """
+        out = self.horner(x)
+        if not hasattr(x, '__len__'):
+            # received a scalar x, squeeze down the array along last dim
+            out = np.squeeze(out, axis=2)
+        if squeeze and self.issiso():
+            # return a scalar/1d array of outputs
+            return out[0][0]
         else:
-            # return a matrix
-            return self.horner(s)
+            return out
+
+    def horner(self, x):
+        """Evaluate system's transfer function at complex frequency
+        using Horner's method.
+
+        Evaluates `sys(x)` where `x` is `s` for continuous-time systems and `z`
+        for discrete-time systems.
+
+        Expects inputs and outputs to be formatted correctly. Use ``sys(x)``
+        for a more user-friendly interface.
+
+        Parameters
+        ----------
+        x : complex array_like or complex
+            Complex frequencies
+
+        Returns
+        -------
+        output : (self.outputs, self.inputs, len(x)) complex ndarray
+            Frequency response
+
+        """
+        x_arr = np.atleast_1d(x) # force to be an array
+        out = empty((self.outputs, self.inputs, len(x_arr)), dtype=complex)
+        for i in range(self.outputs):
+            for j in range(self.inputs):
+                out[i][j] = (polyval(self.num[i][j], x) /
+                             polyval(self.den[i][j], x))
+        return out
 
     def _truncatecoeff(self):
         """Remove extraneous zero coefficients from num and den.
@@ -607,103 +658,19 @@ def __getitem__(self, key):
         else:
             return TransferFunction(num, den, self.dt)
 
-    def evalfr(self, omega):
-        """Evaluate a transfer function at a single angular frequency.
-
-        self._evalfr(omega) returns the value of the transfer function
-        matrix with input value s = i * omega.
-
-        """
-        warn("TransferFunction.evalfr(omega) will be deprecated in a "
-             "future release of python-control; use evalfr(sys, omega*1j) "
-             "instead", PendingDeprecationWarning)
-        return self._evalfr(omega)
-
-    def _evalfr(self, omega):
-        """Evaluate a transfer function at a single angular frequency."""
-        # TODO: implement for discrete time systems
-        if isdtime(self, strict=True):
-            # Convert the frequency to discrete time
-            s = exp(1.j * omega * self.dt)
-            if np.any(omega * self.dt > pi):
-                warn("_evalfr: frequency evaluation above Nyquist frequency")
-        else:
-            s = 1.j * omega
-
-        return self.horner(s)
-
-    def horner(self, s):
-        """Evaluate the systems's transfer function for a complex variable
-
-        Returns a matrix of values evaluated at complex variable s.
-        """
-
-        # Preallocate the output.
-        if getattr(s, '__iter__', False):
-            out = empty((self.outputs, self.inputs, len(s)), dtype=complex)
-        else:
-            out = empty((self.outputs, self.inputs), dtype=complex)
-
-        for i in range(self.outputs):
-            for j in range(self.inputs):
-                out[i][j] = (polyval(self.num[i][j], s) /
-                             polyval(self.den[i][j], s))
-
-        return out
-
     def freqresp(self, omega):
-        """Evaluate the transfer function at a list of angular frequencies.
-
-        Reports the frequency response of the system,
+        """(deprecated) Evaluate transfer function at complex frequencies.
 
-             G(j*omega) = mag*exp(j*phase)
-
-        for continuous time. For discrete time systems, the response is
-        evaluated around the unit circle such that
-
-             G(exp(j*omega*dt)) = mag*exp(j*phase).
-
-        Parameters
-        ----------
-        omega : array_like
-            A list of frequencies in radians/sec at which the system should be
-            evaluated. The list can be either a python list or a numpy array
-            and will be sorted before evaluation.
-
-        Returns
-        -------
-        mag : (self.outputs, self.inputs, len(omega)) ndarray
-            The magnitude (absolute value, not dB or log10) of the system
-            frequency response.
-        phase : (self.outputs, self.inputs, len(omega)) ndarray
-            The wrapped phase in radians of the system frequency response.
-        omega : ndarray or list or tuple
-            The list of sorted frequencies at which the response was
-            evaluated.
+        .. deprecated::0.9.0
+            Method has been given the more pythonic name
+            :meth:`TransferFunction.frequency_response`. Or use
+            :func:`freqresp` in the MATLAB compatibility module.
         """
-        # Preallocate outputs.
-        numfreq = len(omega)
-        mag = empty((self.outputs, self.inputs, numfreq))
-        phase = empty((self.outputs, self.inputs, numfreq))
-
-        # Figure out the frequencies
-        omega.sort()
-        if isdtime(self, strict=True):
-            slist = np.array([exp(1.j * w * self.dt) for w in omega])
-            if max(omega) * self.dt > pi:
-                warn("freqresp: frequency evaluation above Nyquist frequency")
-        else:
-            slist = np.array([1j * w for w in omega])
-
-        # Compute frequency response for each input/output pair
-        for i in range(self.outputs):
-            for j in range(self.inputs):
-                fresp = (polyval(self.num[i][j], slist) /
-                         polyval(self.den[i][j], slist))
-                mag[i, j, :] = abs(fresp)
-                phase[i, j, :] = angle(fresp)
-
-        return mag, phase, omega
+        warn("TransferFunction.freqresp(omega) will be removed in a "
+             "future release of python-control; use "
+             "sys.frequency_response(omega), or freqresp(sys, omega) in the "
+             "MATLAB compatibility module instead", DeprecationWarning)
+        return self.frequency_response(omega)
 
     def pole(self):
         """Compute the poles of a transfer function."""
diff --git a/examples/robust_mimo.py b/examples/robust_mimo.py
index d4e1335e6..9cc5a1c3b 100644
--- a/examples/robust_mimo.py
+++ b/examples/robust_mimo.py
@@ -43,7 +43,7 @@ def triv_sigma(g, w):
     g - LTI object, order n
     w - frequencies, length m
     s - (m,n) array of singular values of g(1j*w)"""
-    m, p, _ = g.freqresp(w)
+    m, p, _ = g.frequency_response(w)
     sjw = (m*np.exp(1j*p)).transpose(2, 0, 1)
     sv = np.linalg.svd(sjw, compute_uv=False)
     return sv
diff --git a/examples/robust_siso.py b/examples/robust_siso.py
index 87fcdb707..17ce10927 100644
--- a/examples/robust_siso.py
+++ b/examples/robust_siso.py
@@ -50,10 +50,10 @@
 
 # frequency response
 omega = np.logspace(-2, 2, 101)
-ws1mag, _, _ = ws1.freqresp(omega)
-s1mag, _, _ = s1.freqresp(omega)
-ws2mag, _, _ = ws2.freqresp(omega)
-s2mag, _, _ = s2.freqresp(omega)
+ws1mag, _, _ = ws1.frequency_response(omega)
+s1mag, _, _ = s1.frequency_response(omega)
+ws2mag, _, _ = ws2.frequency_response(omega)
+s2mag, _, _ = s2.frequency_response(omega)
 
 plt.figure(1)
 # text uses log-scaled absolute, but dB are probably more familiar to most control engineers

From bddc7926f752a9e64b0e0beb939fde515abdf847 Mon Sep 17 00:00:00 2001
From: bnavigator <code@bnavigator.de>
Date: Wed, 6 Jan 2021 23:16:38 +0100
Subject: [PATCH 22/32] fix merge strategy 'errors'

---
 control/frdata.py | 26 +++++---------------------
 control/lti.py    |  1 -
 2 files changed, 5 insertions(+), 22 deletions(-)

diff --git a/control/frdata.py b/control/frdata.py
index 168163d28..32d3d2c00 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -348,36 +348,21 @@ def eval(self, omega, squeeze=True):
 
         Note that a "normal" FRD only returns values for which there is an
         entry in the omega vector. An interpolating FRD can return
-        intermediate values.     
-        
-        Parameters
-        ----------
-        omega: float or array_like
-            Frequencies in radians per second
-        squeeze: bool, optional (default=True)
-            If True and `sys` is single input single output (SISO), returns a 
-            1D array or scalar depending on the length of `omega` 
-            
-        Returns
-        -------
-        fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
-            The frequency response of the system. Array is ``len(x)`` if and only
-            if system is SISO and ``squeeze=True``.
+        intermediate values.
 
         Parameters
         ----------
-        omega: float or array_like
+        omega : float or array_like
             Frequencies in radians per second
-        squeeze: bool, optional (default=True)
+        squeeze : bool, optional (default=True)
             If True and `sys` is single input single output (SISO), returns a
             1D array or scalar depending on the length of `omega`
 
         Returns
         -------
         fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
-            The frequency response of the system. Array is ``len(x)`` if and
-            only if system is SISO and ``squeeze=True``.
-
+            The frequency response of the system. Array is ``len(x)`` if and only
+            if system is SISO and ``squeeze=True``.
         """
         omega_array = np.array(omega, ndmin=1) # array-like version of omega
         if any(omega_array.imag > 0):
@@ -434,7 +419,6 @@ def __call__(self, s, squeeze=True):
             If `s` is not purely imaginary, because
             :class:`FrequencyDomainData` systems are only defined at imaginary
             frequency values.
-
         """
         if any(abs(np.array(s, ndmin=1).real) > 0):
             raise ValueError("__call__: FRD systems can only accept"
diff --git a/control/lti.py b/control/lti.py
index 70fbcdc7f..eff14be2a 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -159,7 +159,6 @@ def dcgain(self):
         """Return the zero-frequency gain"""
         raise NotImplementedError("dcgain not implemented for %s objects" %
                                   str(self.__class__))
-    
 
 
 # Test to see if a system is SISO

From e9bff6a4071012f208f1577336450790978e584c Mon Sep 17 00:00:00 2001
From: bnavigator <code@bnavigator.de>
Date: Wed, 6 Jan 2021 23:29:02 +0100
Subject: [PATCH 23/32] restore test_phase_crossover_frequencies_mimo

---
 control/tests/margin_test.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/control/tests/margin_test.py b/control/tests/margin_test.py
index 81a4b68b4..4965bbe89 100644
--- a/control/tests/margin_test.py
+++ b/control/tests/margin_test.py
@@ -136,10 +136,8 @@ def test_phase_crossover_frequencies_mimo():
                            [[3], [4]]],
                           [[[1, 2, 3, 4], [1, 1]],
                            [[1, 1], [1, 1]]])
-
-    omega, gain = phase_crossover_frequencies(tf[0,0])
-    assert_allclose(omega, [1.73205,  0.], atol=1.5e-3)
-    assert_allclose(gain, [-0.5,  0.25], atol=1.5e-3)
+    with pytest.raises(ControlMIMONotImplemented):
+        omega, gain = phase_crossover_frequencies(tf)
 
 
 def test_mag_phase_omega():

From 73fc6a6056aba8051959df9705cd3f7f88844e13 Mon Sep 17 00:00:00 2001
From: bnavigator <code@bnavigator.de>
Date: Wed, 6 Jan 2021 23:43:04 +0100
Subject: [PATCH 24/32] reinstate second check in frd_test test_eval()

---
 control/tests/frd_test.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/control/tests/frd_test.py b/control/tests/frd_test.py
index 5343112fe..581ee9b25 100644
--- a/control/tests/frd_test.py
+++ b/control/tests/frd_test.py
@@ -401,7 +401,8 @@ def test_operator_conversion(self):
     def test_eval(self):
         sys_tf = ct.tf([1], [1, 2, 1])
         frd_tf = FRD(sys_tf, np.logspace(-1, 1, 3))
-        np.testing.assert_almost_equal(evalfr(sys_tf, 1J), frd_tf.eval(1))
+        np.testing.assert_almost_equal(sys_tf(1j), frd_tf.eval(1))
+        np.testing.assert_almost_equal(sys_tf(1j), frd_tf(1j))
 
         # Should get an error if we evaluate at an unknown frequency
         with pytest.raises(ValueError):

From 7acfc779f10806610db0680c3e305ab9e26beb08 Mon Sep 17 00:00:00 2001
From: bnavigator <code@bnavigator.de>
Date: Thu, 7 Jan 2021 00:01:35 +0100
Subject: [PATCH 25/32] statesp_test: rename from test_evalfr to test_call

---
 control/tests/statesp_test.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 7bf90f4c6..041792377 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -327,7 +327,7 @@ def test_multiply_ss(self, sys222, sys322):
                                          [-3.00137118e-01+3.42881660e-03j,
                                           6.32015038e-04-1.21462255e-02j]]))])
     @pytest.mark.parametrize("dt", [None, 0, 1e-3])
-    def test_evalfr(self, dt, omega, resp):
+    def test_call(self, dt, omega, resp):
         """Evaluate the frequency response at single frequencies"""
         A = [[-2, 0.5], [0.5, -0.3]]
         B = [[0.3, -1.3], [0.1, 0.]]
@@ -345,9 +345,9 @@ def test_evalfr(self, dt, omega, resp):
         np.testing.assert_allclose(evalfr(sys, s), resp, atol=1e-3)
         np.testing.assert_allclose(sys(s), resp, atol=1e-3)
 
-        # Deprecated version of the call (should generate error)
+        # Deprecated name of the call (should generate error)
         with pytest.raises(AttributeError):
-            np.testing.assert_allclose(sys.evalfr(omega), resp, atol=1e-3)
+            sys.evalfr(omega)
 
 
     @slycotonly

From a26520d6b1d124b5f13c6083082b3fde5971099c Mon Sep 17 00:00:00 2001
From: bnavigator <code@bnavigator.de>
Date: Thu, 7 Jan 2021 00:02:08 +0100
Subject: [PATCH 26/32] statesp_test: remove duplicate test_pole_static and
 test_copy_constructor

---
 control/tests/statesp_test.py | 25 -------------------------
 1 file changed, 25 deletions(-)

diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 041792377..27c37e247 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -922,28 +922,3 @@ def test_latex_repr(gmats, ref, repr_type, num_format, editsdefaults):
     g = StateSpace(*gmats)
     refkey = "{}_{}".format(refkey_n[num_format], refkey_r[repr_type])
     assert g._repr_latex_() == ref[refkey]
-
-    def test_pole_static(self):
-        """Regression: pole() of static gain is empty array."""
-        np.testing.assert_array_equal(np.array([]),
-                                      StateSpace([], [], [], [[1]]).pole())
-
-    def test_copy_constructor(self):
-        # Create a set of matrices for a simple linear system
-        A = np.array([[-1]])
-        B = np.array([[1]])
-        C = np.array([[1]])
-        D = np.array([[0]])
-
-        # Create the first linear system and a copy
-        linsys = StateSpace(A, B, C, D)
-        cpysys = StateSpace(linsys)
-
-        # Change the original A matrix
-        A[0, 0] = -2
-        np.testing.assert_array_equal(linsys.A, [[-1]]) # original value
-        np.testing.assert_array_equal(cpysys.A, [[-1]]) # original value
-
-        # Change the A matrix for the original system
-        linsys.A[0, 0] = -3
-        np.testing.assert_array_equal(cpysys.A, [[-1]]) # original value

From dae8d4292bab4a9b12a6433bfd5ea4f48eec4d18 Mon Sep 17 00:00:00 2001
From: sawyerbfuller <58706249+sawyerbfuller@users.noreply.github.com>
Date: Wed, 6 Jan 2021 15:17:26 -0800
Subject: [PATCH 27/32] Update control/xferfcn.py

Co-authored-by: Ben Greiner <code@bnavigator.de>
---
 control/xferfcn.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/control/xferfcn.py b/control/xferfcn.py
index a0941b74a..4c0d26a07 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -279,7 +279,7 @@ def horner(self, x):
 
         Parameters
         ----------
-        x : complex array_like or complex
+        x : complex array_like or complex scalar
             Complex frequencies
 
         Returns

From 3b10af001bf07644062c0c2da0067bc501cd3eff Mon Sep 17 00:00:00 2001
From: bnavigator <code@bnavigator.de>
Date: Thu, 7 Jan 2021 00:27:32 +0100
Subject: [PATCH 28/32] check for warning the pytest way

---
 control/tests/xferfcn_test.py | 9 +--------
 1 file changed, 1 insertion(+), 8 deletions(-)

diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py
index 1ec596a6e..eb8755f82 100644
--- a/control/tests/xferfcn_test.py
+++ b/control/tests/xferfcn_test.py
@@ -474,15 +474,8 @@ def test_call_mimo(self):
     def test_freqresp_deprecated(self):
         sys = TransferFunction([1., 3., 5], [1., 6., 2., -1.])
         # Deprecated version of the call (should generate warning)
-        import warnings
-        with warnings.catch_warnings(record=True) as w:
-            # Set up warnings filter to only show warnings in control module
-            warnings.filterwarnings("ignore")
-            warnings.filterwarnings("always", module="control")
-
-            # Make sure that we get a pending deprecation warning
+        with pytest.warns(DeprecationWarning):
             sys.freqresp(1.)
-            assert issubclass(w[-1].category, DeprecationWarning)
 
     def test_frequency_response_siso(self):
         """Evaluate the magnitude and phase of a SISO system at

From 1bf445bbc885c3ccf1fb7c9ead3538262afc983a Mon Sep 17 00:00:00 2001
From: Richard Murray <murray@cds.caltech.edu>
Date: Wed, 6 Jan 2021 21:42:12 -0800
Subject: [PATCH 29/32] add unit tests for freqresp/evalfr warnings/errors

---
 control/frdata.py             | 2 +-
 control/tests/frd_test.py     | 8 ++++++--
 control/tests/statesp_test.py | 6 ++++++
 3 files changed, 13 insertions(+), 3 deletions(-)

diff --git a/control/frdata.py b/control/frdata.py
index 32d3d2c00..837bc1dac 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -421,7 +421,7 @@ def __call__(self, s, squeeze=True):
             frequency values.
         """
         if any(abs(np.array(s, ndmin=1).real) > 0):
-            raise ValueError("__call__: FRD systems can only accept"
+            raise ValueError("__call__: FRD systems can only accept "
                             "purely imaginary frequencies")
         # need to preserve array or scalar status
         if hasattr(s, '__len__'):
diff --git a/control/tests/frd_test.py b/control/tests/frd_test.py
index 581ee9b25..7418dddda 100644
--- a/control/tests/frd_test.py
+++ b/control/tests/frd_test.py
@@ -405,11 +405,15 @@ def test_eval(self):
         np.testing.assert_almost_equal(sys_tf(1j), frd_tf(1j))
 
         # Should get an error if we evaluate at an unknown frequency
-        with pytest.raises(ValueError):
+        with pytest.raises(ValueError, match="not .* in frequency list"):
             frd_tf.eval(2)
 
+        # Should get an error if we evaluate at an complex number
+        with pytest.raises(ValueError, match="can only accept real-valued"):
+            frd_tf.eval(2 + 1j)
+
         # Should get an error if we use __call__ at real-valued frequency
-        with pytest.raises(ValueError):
+        with pytest.raises(ValueError, match="only accept purely imaginary"):
             frd_tf(2)
 
     def test_repr_str(self):
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 27c37e247..25cfe6f47 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -376,6 +376,12 @@ def test_freq_resp(self):
         np.testing.assert_almost_equal(phase, true_phase)
         np.testing.assert_equal(omega, true_omega)
 
+        # Deprecated version of the call (should return warning)
+        with pytest.warns(DeprecationWarning, match="will be removed"):
+            from control import freqresp
+            mag, phase, omega = sys.freqresp(true_omega)
+            np.testing.assert_almost_equal(mag, true_mag)
+
     def test_is_static_gain(self):
         A0 = np.zeros((2,2))
         A1 = A0.copy()

From 3dbaacd7d9d03da4558ad1a8b88dcdcbe47e9261 Mon Sep 17 00:00:00 2001
From: bnavigator <code@bnavigator.de>
Date: Thu, 7 Jan 2021 13:03:57 +0100
Subject: [PATCH 30/32] test frd freqresp deprecation

---
 control/tests/frd_test.py      |  6 ++++
 control/tests/freqresp_test.py | 54 +++++++++++++++++++---------------
 control/tests/statesp_test.py  |  1 -
 3 files changed, 37 insertions(+), 24 deletions(-)

diff --git a/control/tests/frd_test.py b/control/tests/frd_test.py
index 7418dddda..f8ee3eb20 100644
--- a/control/tests/frd_test.py
+++ b/control/tests/frd_test.py
@@ -416,6 +416,12 @@ def test_eval(self):
         with pytest.raises(ValueError, match="only accept purely imaginary"):
             frd_tf(2)
 
+    def test_freqresp_deprecated(self):
+        sys_tf = ct.tf([1], [1, 2, 1])
+        frd_tf = FRD(sys_tf, np.logspace(-1, 1, 3))
+        with pytest.warns(DeprecationWarning):
+            frd_tf.freqresp(1.)
+
     def test_repr_str(self):
         # repr printing
         array = np.array
diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py
index e6a6934e8..752dacb79 100644
--- a/control/tests/freqresp_test.py
+++ b/control/tests/freqresp_test.py
@@ -21,24 +21,46 @@
 pytestmark = pytest.mark.usefixtures("mplcleanup")
 
 
-def test_siso():
-    """Test SISO frequency response"""
+@pytest.fixture
+def ss_siso():
     A = np.array([[1, 1], [0, 1]])
     B = np.array([[0], [1]])
     C = np.array([[1, 0]])
     D = 0
-    sys = StateSpace(A, B, C, D)
+    return StateSpace(A, B, C, D)
+
+
+@pytest.fixture
+def ss_mimo():
+    A = np.array([[1, 1], [0, 1]])
+    B = np.array([[1, 0], [0, 1]])
+    C = np.array([[1, 0]])
+    D = np.array([[0, 0]])
+    return StateSpace(A, B, C, D)
+
+def test_freqresp_siso(ss_siso):
+    """Test SISO frequency response"""
     omega = np.linspace(10e-2, 10e2, 1000)
 
     # test frequency response
-    sys.frequency_response(omega)
+    ctrl.freqresp(ss_siso, omega)
 
-    # test bode plot
-    bode(sys)
 
-    # Convert to transfer function and test bode
-    systf = tf(sys)
-    bode(systf)
+@slycotonly
+def test_freqresp_mimo(ss_mimo):
+    """Test MIMO frequency response calls"""
+    omega = np.linspace(10e-2, 10e2, 1000)
+    ctrl.freqresp(ss_mimo, omega)
+    tf_mimo = tf(ss_mimo)
+    ctrl.freqresp(tf_mimo, omega)
+
+
+def test_bode_basic(ss_siso):
+    """Test bode plot call (Very basic)"""
+    # TODO: proper test
+    tf_siso = tf(ss_siso)
+    bode(ss_siso)
+    bode(tf_siso)
 
 
 @pytest.mark.filterwarnings("ignore:.*non-positive left xlim:UserWarning")
@@ -106,20 +128,6 @@ def test_doubleint():
     bode(sys)
 
 
-@slycotonly
-def test_mimo():
-    """Test MIMO frequency response calls"""
-    A = np.array([[1, 1], [0, 1]])
-    B = np.array([[1, 0], [0, 1]])
-    C = np.array([[1, 0]])
-    D = np.array([[0, 0]])
-    omega = np.linspace(10e-2, 10e2, 1000)
-    sysMIMO = ss(A, B, C, D)
-
-    sysMIMO.frequency_response(omega)
-    tf(sysMIMO)
-
-
 @pytest.mark.parametrize(
     "Hz, Wcp, Wcg",
     [pytest.param(False, 6.0782869, 10., id="omega"),
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 25cfe6f47..48f27a3b5 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -378,7 +378,6 @@ def test_freq_resp(self):
 
         # Deprecated version of the call (should return warning)
         with pytest.warns(DeprecationWarning, match="will be removed"):
-            from control import freqresp
             mag, phase, omega = sys.freqresp(true_omega)
             np.testing.assert_almost_equal(mag, true_mag)
 

From 178bfa0e6cc6bd24b6f2f6aff906f940ea3b9183 Mon Sep 17 00:00:00 2001
From: sawyerbfuller <58706249+sawyerbfuller@users.noreply.github.com>
Date: Thu, 7 Jan 2021 11:40:01 -0800
Subject: [PATCH 31/32] doc updates to specify frequency response outputs more
 precisely

Co-authored-by: Ben Greiner <code@bnavigator.de>
---
 control/frdata.py  | 27 ++++++++++---------
 control/lti.py     | 67 +++++++++++++++++++++++++++-------------------
 control/statesp.py | 14 ++++++----
 3 files changed, 64 insertions(+), 44 deletions(-)

diff --git a/control/frdata.py b/control/frdata.py
index 837bc1dac..22dbb298f 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -356,11 +356,11 @@ def eval(self, omega, squeeze=True):
             Frequencies in radians per second
         squeeze : bool, optional (default=True)
             If True and `sys` is single input single output (SISO), returns a
-            1D array or scalar depending on the length of `omega`
+            1D array rather than a 3D array.
 
         Returns
         -------
-        fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
+        fresp : (self.outputs, self.inputs, len(x)) or (len(x), ) complex ndarray
             The frequency response of the system. Array is ``len(x)`` if and only
             if system is SISO and ``squeeze=True``.
         """
@@ -393,25 +393,28 @@ def eval(self, omega, squeeze=True):
 
     def __call__(self, s, squeeze=True):
         """Evaluate system's transfer function at complex frequencies.
-
-        Returns the complex frequency response `sys(s)`.
+        
+        Returns the complex frequency response `sys(s)` of system `sys` with
+        `m = sys.inputs` number of inputs and `p = sys.outputs` number of
+        outputs.
 
         To evaluate at a frequency omega in radians per second, enter
-        ``x = omega * 1j`` or use ``sys.eval(omega)``
+        ``s = omega * 1j`` or use ``sys.eval(omega)``
 
         Parameters
         ----------
-        s: complex scalar or array_like
+        s : complex scalar or array_like
             Complex frequencies
-        squeeze: bool, optional (default=True)
-            If True and `sys` is single input single output (SISO), returns a
-            1D array or scalar depending on the length of x`.
+        squeeze : bool, optional (default=True)
+            If True and `sys` is single input single output (SISO), i.e. `m=1`,
+            `p=1`, return a 1D array rather than a 3D array.
 
         Returns
         -------
-        fresp : (self.outputs, self.inputs, len(s)) or len(s) complex ndarray
-            The frequency response of the system. Array is ``len(s)`` if and
-            only if system is SISO and ``squeeze=True``.
+        fresp : (p, m, len(s)) complex ndarray or (len(s),) complex ndarray 
+            The frequency response of the system. Array is ``(len(s), )`` if
+            and only if system is SISO and ``squeeze=True``.
+
 
         Raises
         ------
diff --git a/control/lti.py b/control/lti.py
index eff14be2a..152c5c73b 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -124,25 +124,30 @@ def frequency_response(self, omega, squeeze=True):
 
              G(exp(j*omega*dt)) = mag*exp(j*phase).
 
+        In general the system may be multiple input, multiple output (MIMO), where
+        `m = self.inputs` number of inputs and `p = self.outputs` number of
+        outputs.
+
         Parameters
         ----------
-        omega : array_like or float
+        omega : float or array_like
             A list, tuple, array, or scalar value of frequencies in
             radians/sec at which the system will be evaluated.
-        squeeze: bool, optional (default=True)
-            If True and sys is single input, single output (SISO), return a
-            1D array or scalar depending on omega's length.
+        squeeze : bool, optional (default=True)
+            If True and the system is single input single output (SISO), i.e. `m=1`,
+            `p=1`, return a 1D array rather than a 3D array.
 
         Returns
         -------
-        mag : (self.outputs, self.inputs, len(omega)) or len(omega) ndarray
+        mag : (p, m, len(omega)) ndarray or (len(omega),) ndarray 
             The magnitude (absolute value, not dB or log10) of the system
-            frequency response.
-        phase : (self.outputs, self.inputs, len(omega)) or len(omega) ndarray
+            frequency response. Array is ``(len(omega), )`` if
+            and only if system is SISO and ``squeeze=True``.
+        phase : (p, m, len(omega)) ndarray or (len(omega),) ndarray
             The wrapped phase in radians of the system frequency response.
         omega : ndarray
             The (sorted) frequencies at which the response was evaluated.
-
+            
         """
         omega = np.sort(np.array(omega, ndmin=1))
         if isdtime(self, strict=True):
@@ -463,7 +468,9 @@ def evalfr(sys, x, squeeze=True):
     Evaluate the transfer function of an LTI system for complex frequency x.
 
     Returns the complex frequency response `sys(x)` where `x` is `s` for
-    continuous-time systems and `z` for discrete-time systems.
+    continuous-time systems and `z` for discrete-time systems, with
+    `m = sys.inputs` number of inputs and `p = sys.outputs` number of
+    outputs.
 
     To evaluate at a frequency omega in radians per second, enter
     ``x = omega * 1j`` for continuous-time systems, or
@@ -474,17 +481,17 @@ def evalfr(sys, x, squeeze=True):
     ----------
     sys: StateSpace or TransferFunction
         Linear system
-    x: complex scalar or array_like
+    x : complex scalar or array_like
         Complex frequency(s)
-    squeeze: bool, optional (default=True)
-        If True and sys is single input single output (SISO), returns a
-        1D array or scalar depending on the length of x.
+    squeeze : bool, optional (default=True)
+        If True and `sys` is single input single output (SISO), i.e. `m=1`,
+        `p=1`, return a 1D array rather than a 3D array.
 
     Returns
     -------
-    fresp : (sys.outputs, sys.inputs, len(x)) or len(x) complex ndarray
-        The frequency response of the system. Array is len(x) if and only if
-        system is SISO and squeeze=True.
+    fresp : (p, m, len(x)) complex ndarray or (len(x),) complex ndarray 
+        The frequency response of the system. Array is ``(len(x), )`` if
+        and only if system is SISO and ``squeeze=True``.
 
     See Also
     --------
@@ -493,8 +500,8 @@ def evalfr(sys, x, squeeze=True):
 
     Notes
     -----
-    This function is a wrapper for StateSpace.__call__ and
-    TransferFunction.__call__.
+    This function is a wrapper for :meth:`StateSpace.__call__` and
+    :meth:`TransferFunction.__call__`.
 
     Examples
     --------
@@ -510,25 +517,31 @@ def evalfr(sys, x, squeeze=True):
 def freqresp(sys, omega, squeeze=True):
     """
     Frequency response of an LTI system at multiple angular frequencies.
+    
+    In general the system may be multiple input, multiple output (MIMO), where
+    `m = sys.inputs` number of inputs and `p = sys.outputs` number of
+    outputs.
 
     Parameters
     ----------
     sys: StateSpace or TransferFunction
         Linear system
-    omega: float or array_like
+    omega : float or array_like
         A list of frequencies in radians/sec at which the system should be
         evaluated. The list can be either a python list or a numpy array
         and will be sorted before evaluation.
-    squeeze: bool, optional (default=True)
-        If True and sys is single input, single output (SISO), returns
-        1D array or scalar depending on omega's length.
+    squeeze : bool, optional (default=True)
+        If True and `sys` is single input, single output (SISO), returns
+        1D array rather than a 3D array.
 
     Returns
     -------
-    mag : (sys.outputs, sys.inputs, len(omega)) or len(omega) ndarray
+    mag : (p, m, len(omega)) ndarray or (len(omega),) ndarray 
         The magnitude (absolute value, not dB or log10) of the system
-        frequency response.
-    phase : (sys.outputs, sys.inputs, len(omega)) or len(omega) ndarray
+        frequency response. Array is ``(len(omega), )`` if and only if system 
+        is SISO and ``squeeze=True``.
+
+    phase : (p, m, len(omega)) ndarray or (len(omega),) ndarray 
         The wrapped phase in radians of the system frequency response.
     omega : ndarray
         The list of sorted frequencies at which the response was
@@ -541,8 +554,8 @@ def freqresp(sys, omega, squeeze=True):
 
     Notes
     -----
-    This function is a wrapper for StateSpace.frequency_response and
-    TransferFunction.frequency_response.
+    This function is a wrapper for :meth:`StateSpace.frequency_response` and
+    :meth:`TransferFunction.frequency_response`.
 
     Examples
     --------
diff --git a/control/statesp.py b/control/statesp.py
index 8bf4c0d99..ffd229108 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -645,6 +645,10 @@ def __call__(self, x, squeeze=True):
 
         Returns the complex frequency response `sys(x)` where `x` is `s` for
         continuous-time systems and `z` for discrete-time systems.
+        
+        In general the system may be multiple input, multiple output (MIMO), where
+        `m = self.inputs` number of inputs and `p = self.outputs` number of
+        outputs.
 
         To evaluate at a frequency omega in radians per second, enter
         ``x = omega * 1j``, for continuous-time systems, or
@@ -653,15 +657,15 @@ def __call__(self, x, squeeze=True):
 
         Parameters
         ----------
-        x: complex or complex array_like
+        x : complex or complex array_like
             Complex frequencies
-        squeeze: bool, optional (default=True)
-            If True and `sys` is single input single output (SISO), returns a
-            1D array or scalar depending on the length of `x`.
+        squeeze : bool, optional (default=True)
+            If True and `self` is single input single output (SISO), returns a
+            1D array rather than a 3D array.
 
         Returns
         -------
-        fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
+        fresp : (p, m, len(x)) complex ndarray or (len(x),) complex ndarray
             The frequency response of the system. Array is ``len(x)`` if and
             only if system is SISO and ``squeeze=True``.
 

From a6310988374348128f038420d480d593f7995765 Mon Sep 17 00:00:00 2001
From: sawyerbfuller <58706249+sawyerbfuller@users.noreply.github.com>
Date: Thu, 7 Jan 2021 11:57:47 -0800
Subject: [PATCH 32/32] missed a couple of doc updates in xferfcn

Co-authored-by: Ben Greiner <code@bnavigator.de>
---
 control/xferfcn.py | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/control/xferfcn.py b/control/xferfcn.py
index 4c0d26a07..fba674efe 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -237,6 +237,10 @@ def __call__(self, x, squeeze=True):
         Returns the complex frequency response `sys(x)` where `x` is `s` for
         continuous-time systems and `z` for discrete-time systems.
 
+        In general the system may be multiple input, multiple output (MIMO), where
+        `m = self.inputs` number of inputs and `p = self.outputs` number of
+        outputs.
+        
         To evaluate at a frequency omega in radians per second, enter
         ``x = omega * 1j``, for continuous-time systems, or
         ``x = exp(1j * omega * dt)`` for discrete-time systems. Or use
@@ -244,15 +248,15 @@ def __call__(self, x, squeeze=True):
 
         Parameters
         ----------
-        x: complex array_like or complex
+        x : complex array_like or complex
             Complex frequencies
-        squeeze: bool, optional (default=True)
+        squeeze : bool, optional (default=True)
             If True and `sys` is single input single output (SISO), returns a
-            1D array or scalar depending on the length of `x`.
+            1D array rather than a 3D array.
 
         Returns
         -------
-        fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
+        fresp : (p, m, len(x)) complex ndarray or or (len(x), ) complex ndarray
             The frequency response of the system. Array is `len(x)` if and
             only if system is SISO and ``squeeze=True``.