diff --git a/README.rst b/README.rst
index 48b398bc8..7a30f2cb8 100644
--- a/README.rst
+++ b/README.rst
@@ -24,7 +24,7 @@ Features
 Links
 =====
 
-- Project home page: http://python-control.sourceforge.net
+- Project home page: http://python-control.org
 - Source code repository: https://github.com/python-control/python-control
 - Documentation: http://python-control.readthedocs.org/
 - Issue tracker: https://github.com/python-control/python-control/issues
@@ -46,7 +46,7 @@ https://github.com/python-control/Slycot
 Installation
 ============
 
-The package may be installed using pip or distutils.
+The package may be installed using pip, conda, or distutils.
 
 Pip
 ---
@@ -56,6 +56,14 @@ To install using pip::
   pip install slycot   # optional
   pip install control
 
+conda-forge
+-----------
+
+Binaries are available from conda-forge for selected platforms (Linux and
+MacOS).  Install using
+
+  conda install -c conda-forge control
+
 Distutils
 ---------
 
diff --git a/control/canonical.py b/control/canonical.py
index 4acd78ff8..c0244d75f 100644
--- a/control/canonical.py
+++ b/control/canonical.py
@@ -6,8 +6,8 @@
 from .statesp import StateSpace
 from .statefbk import ctrb, obsv
 
-from numpy import zeros, shape, poly
-from numpy.linalg import solve, matrix_rank
+from numpy import zeros, shape, poly, iscomplex, hstack
+from numpy.linalg import solve, matrix_rank, eig
 
 __all__ = ['canonical_form', 'reachable_form', 'observable_form']
 
@@ -22,7 +22,7 @@ def canonical_form(xsys, form='reachable'):
         Canonical form for transformation.  Chosen from:
           * 'reachable' - reachable canonical form
           * 'observable' - observable canonical form
-          * 'modal' - modal canonical form [not implemented]
+          * 'modal' - modal canonical form
 
     Returns
     -------
@@ -37,6 +37,8 @@ def canonical_form(xsys, form='reachable'):
         return reachable_form(xsys)
     elif form == 'observable':
         return observable_form(xsys)
+    elif form == 'modal':
+        return modal_form(xsys)
     else:
         raise ControlNotImplemented(
             "Canonical form '%s' not yet implemented" % form)
@@ -142,3 +144,69 @@ def observable_form(xsys):
     zsys.B = Tzx * xsys.B
 
     return zsys, Tzx
+
+def modal_form(xsys):
+    """Convert a system into modal canonical form
+
+    Parameters
+    ----------
+    xsys : StateSpace object
+        System to be transformed, with state `x`
+
+    Returns
+    -------
+    zsys : StateSpace object
+        System in modal canonical form, with state `z`
+    T : matrix
+        Coordinate transformation: z = T * x
+    """
+    # Check to make sure we have a SISO system
+    if not issiso(xsys):
+        raise ControlNotImplemented(
+            "Canonical forms for MIMO systems not yet supported")
+
+    # Create a new system, starting with a copy of the old one
+    zsys = StateSpace(xsys)
+
+    # Calculate eigenvalues and matrix of eigenvectors Tzx,
+    eigval, eigvec = eig(xsys.A)
+
+    # Eigenvalues and according eigenvectors are not sorted,
+    # thus modal transformation is ambiguous
+    # Sorting eigenvalues and respective vectors by largest to smallest eigenvalue
+    idx = eigval.argsort()[::-1]
+    eigval = eigval[idx]
+    eigvec = eigvec[:,idx]
+
+    # If all eigenvalues are real, the matrix of eigenvectors is Tzx directly
+    if not iscomplex(eigval).any():
+        Tzx = eigvec
+    else:
+        # A is an arbitrary semisimple matrix
+
+        # Keep track of complex conjugates (need only one)
+        lst_conjugates = []
+        Tzx = None
+        for val, vec in zip(eigval, eigvec.T):
+            if iscomplex(val):
+                if val not in lst_conjugates:
+                    lst_conjugates.append(val.conjugate())
+                    if Tzx is not None:
+                        Tzx = hstack((Tzx, hstack((vec.real.T, vec.imag.T))))
+                    else:
+                        Tzx = hstack((vec.real.T, vec.imag.T))
+                else:
+                    # if conjugate has already been seen, skip this eigenvalue
+                    lst_conjugates.remove(val)
+            else:
+                if Tzx is not None:
+                    Tzx = hstack((Tzx, vec.real.T))
+                else:
+                    Tzx = vec.real.T
+
+    # Generate the system matrices for the desired canonical form
+    zsys.A = solve(Tzx, xsys.A).dot(Tzx)
+    zsys.B = solve(Tzx, xsys.B)
+    zsys.C = xsys.C.dot(Tzx)
+
+    return zsys, Tzx
diff --git a/control/config.py b/control/config.py
index 98198aae7..10ab8a1ed 100644
--- a/control/config.py
+++ b/control/config.py
@@ -18,6 +18,8 @@
 def use_matlab_defaults():
     """
     Use MATLAB compatible configuration settings
+
+    The following conventions are used:
         * Bode plots plot gain in dB, phase in degrees, frequency in Hertz
     """
     # Bode plot defaults
@@ -28,7 +30,9 @@ def use_matlab_defaults():
 # Set defaults to match FBS (Astrom and Murray)
 def use_fbs_defaults():
     """
-    Use `Astrom and Murray <http://fbsbook.org>`_ compatible settings
+    Use `Feedback Systems <http://fbsbook.org>`_ (FBS) compatible settings
+
+    The following conventions are used:
         * Bode plots plot gain in powers of ten, phase in degrees,
           frequency in Hertz
     """
diff --git a/control/frdata.py b/control/frdata.py
index 8c4e1f2f6..d34200455 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -452,8 +452,8 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1):
     scalar, then the number of inputs and outputs can be specified
     manually, as in:
 
-    >>> sys = _convertToFRD(3.) # Assumes inputs = outputs = 1
-    >>> sys = _convertToFRD(1., inputs=3, outputs=2)
+    >>> frd = _convertToFRD(3., omega) # Assumes inputs = outputs = 1
+    >>> frd = _convertToFRD(1., omegs, inputs=3, outputs=2)
 
     In the latter example, sys's matrix transfer function is [[1., 1., 1.]
                                                               [1., 1., 1.]].
diff --git a/control/freqplot.py b/control/freqplot.py
index 33cf0ad61..f5761e282 100644
--- a/control/freqplot.py
+++ b/control/freqplot.py
@@ -70,8 +70,8 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
     ----------
     syslist : linsys
         List of linear input/output systems (single system is OK)
-    omega : freq_range
-        Range of frequencies in rad/sec
+    omega : list
+        List of frequencies in rad/sec to be used for frequency response
     dB : boolean
         If True, plot result in dB
     Hz : boolean
@@ -106,7 +106,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
     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
-    time base.  If not timebase is specified (dt = True), dt is set to 1.
+    timebase.  If not timebase is specified (dt = True), dt is set to 1.
 
     Examples
     --------
diff --git a/control/lti.py b/control/lti.py
index 0bd05e902..5950d9d58 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -82,9 +82,21 @@ def isctime(self, strict=False):
         return self.dt == 0
 
     def issiso(self):
+        '''Check to see if a system is single input, single output'''
         return self.inputs == 1 and self.outputs == 1
 
     def damp(self):
+        '''Natural frequency, damping ratio of system poles
+
+        Returns
+        -------
+        wn : array
+            Natural frequencies for each system pole
+        zeta : array
+            Damping ratio for each system pole
+        poles : array
+            Array of system poles
+        '''
         poles = self.pole()
 
         if isdtime(self, strict=True):
@@ -102,6 +114,16 @@ def dcgain(self):
 
 # Test to see if a system is SISO
 def issiso(sys, strict=False):
+    """
+    Check to see if a system is single input, single output
+
+    Parameters
+    ----------
+    sys : LTI system
+        System to be checked
+    strict: bool (default = False)
+        If strict is True, do not treat scalars as SISO
+    """
     if isinstance(sys, (int, float, complex, np.number)) and not strict:
         return True
     elif not isinstance(sys, LTI):
diff --git a/control/statefbk.py b/control/statefbk.py
index a2d75c30d..0fb377a47 100644
--- a/control/statefbk.py
+++ b/control/statefbk.py
@@ -113,11 +113,11 @@ def place(A, B, p):
     return K
 
 
-def place_varga(A, B, p):
+def place_varga(A, B, p, dtime=False, alpha=None):
     """Place closed loop eigenvalues
-    K = place_varga(A, B, p)
+    K = place_varga(A, B, p, dtime=False, alpha=None)
 
-    Parameters
+    Required Parameters
     ----------
     A : 2-d array
         Dynamics matrix
@@ -125,6 +125,20 @@ def place_varga(A, B, p):
         Input matrix
     p : 1-d list
         Desired eigenvalue locations
+
+    Optional Parameters
+    ---------------
+    dtime: False for continuous time pole placement or True for discrete time.
+            The default is dtime=False.
+    alpha: double scalar
+           If DICO='C', then place_varga will leave the eigenvalues with real
+           real part less than alpha untouched.
+           If DICO='D', the place_varga will leave eigenvalues with modulus
+           less than alpha untouched.
+
+           By default (alpha=None), place_varga computes alpha such that all
+           poles will be placed.
+
     Returns
     -------
     K : 2-d array
@@ -146,7 +160,7 @@ def place_varga(A, B, p):
     --------
     >>> A = [[-1, -1], [0, 1]]
     >>> B = [[0], [1]]
-    >>> K = place(A, B, [-2, -5])
+    >>> K = place_varga(A, B, [-2, -5])
 
     See Also:
     --------
@@ -160,24 +174,46 @@ def place_varga(A, B, p):
         raise ControlSlycot("can't find slycot module 'sb01bd'")
 
     # Convert the system inputs to NumPy arrays
-    A_mat = np.array(A);
-    B_mat = np.array(B);
+    A_mat = np.array(A)
+    B_mat = np.array(B)
     if (A_mat.shape[0] != A_mat.shape[1] or
         A_mat.shape[0] != B_mat.shape[0]):
         raise ControlDimension("matrix dimensions are incorrect")
 
     # Compute the system eigenvalues and convert poles to numpy array
     system_eigs = np.linalg.eig(A_mat)[0]
-    placed_eigs = np.array(p);
+    placed_eigs = np.array(p)
+
+    # Need a character parameter for SB01BD
+    if dtime:
+        DICO = 'D'
+    else:
+        DICO = 'C'
+
+    if alpha is None:
+        # SB01BD ignores eigenvalues with real part less than alpha
+        # (if DICO='C') or with modulus less than alpha
+        # (if DICO = 'D').
+        if dtime:
+            # For discrete time, slycot only cares about modulus, so just make
+            # alpha the smallest it can be.
+            alpha = 0.0
+        else:
+            # Choosing alpha=min_eig is insufficient and can lead to an
+            # error or not having all the eigenvalues placed that we wanted.
+            # Evidently, what python thinks are the eigs is not precisely
+            # the same as what slicot thinks are the eigs. So we need some
+            # numerical breathing room. The following is pretty heuristic,
+            # but does the trick
+            alpha = -2*abs(min(system_eigs.real))
+    elif dtime and alpha < 0.0:
+        raise ValueError("Need alpha > 0 when DICO='D'")
 
-    # SB01BD sets eigenvalues with real part less than alpha
-    # We want to place all poles of the system => set alpha to minimum
-    alpha = min(system_eigs.real);
 
     # Call SLICOT routine to place the eigenvalues
     A_z,w,nfp,nap,nup,F,Z = \
         sb01bd(B_mat.shape[0], B_mat.shape[1], len(placed_eigs), alpha,
-               A_mat, B_mat, placed_eigs, 'C');
+               A_mat, B_mat, placed_eigs, DICO)
 
     # Return the gain matrix, with MATLAB gain convention
     return -F
diff --git a/control/statesp.py b/control/statesp.py
index 4a617abd2..6385642eb 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -385,7 +385,7 @@ def evalfr(self, omega):
         with input value s = i * omega.
 
         """
-        warn("StateSpace.evalfr(omega) will be depracted in a future "
+        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)
@@ -1248,7 +1248,7 @@ def tf2ss(*args):
 
 def rss(states=1, outputs=1, inputs=1):
     """
-    Create a stable **continuous** random state space object.
+    Create a stable *continuous* random state space object.
 
     Parameters
     ----------
@@ -1285,7 +1285,7 @@ def rss(states=1, outputs=1, inputs=1):
 
 def drss(states=1, outputs=1, inputs=1):
     """
-    Create a stable **discrete** random state space object.
+    Create a stable *discrete* random state space object.
 
     Parameters
     ----------
diff --git a/control/tests/canonical_test.py b/control/tests/canonical_test.py
index a450453c1..5c73c1f49 100644
--- a/control/tests/canonical_test.py
+++ b/control/tests/canonical_test.py
@@ -53,6 +53,36 @@ def test_unreachable_system(self):
         # Check if an exception is raised
         np.testing.assert_raises(ValueError, canonical_form, sys, "reachable")
 
+    def test_modal_form(self):
+        """Test the modal canonical form"""
+
+        # Create a system in the modal canonical form
+        A_true = np.diag([4.0, 3.0, 2.0, 1.0]) # order from the largest to the smallest
+        B_true = np.matrix("1.1 2.2 3.3 4.4").T
+        C_true = np.matrix("1.3 1.4 1.5 1.6")
+        D_true = 42.0
+
+        # Perform a coordinate transform with a random invertible matrix
+        T_true = np.matrix([[-0.27144004, -0.39933167,  0.75634684,  0.44135471],
+                            [-0.74855725, -0.39136285, -0.18142339, -0.50356997],
+                            [-0.40688007,  0.81416369,  0.38002113, -0.16483334],
+                            [-0.44769516,  0.15654653, -0.50060858,  0.72419146]])
+        A = np.linalg.solve(T_true, A_true)*T_true
+        B = np.linalg.solve(T_true, B_true)
+        C = C_true*T_true
+        D = D_true
+
+        # Create a state space system and convert it to the modal canonical form
+        sys_check, T_check = canonical_form(ss(A, B, C, D), "modal")
+
+        # Check against the true values
+        #TODO: Test in respect to ambiguous transformation (system characteristics?)
+        np.testing.assert_array_almost_equal(sys_check.A, A_true)
+        #np.testing.assert_array_almost_equal(sys_check.B, B_true)
+        #np.testing.assert_array_almost_equal(sys_check.C, C_true)
+        np.testing.assert_array_almost_equal(sys_check.D, D_true)
+        #np.testing.assert_array_almost_equal(T_check, T_true)
+        
     def test_observable_form(self):
         """Test the observable canonical form"""
 
diff --git a/control/tests/convert_test.py b/control/tests/convert_test.py
index 5d9012399..0340fa718 100644
--- a/control/tests/convert_test.py
+++ b/control/tests/convert_test.py
@@ -245,6 +245,29 @@ def testTf2SsDuplicatePoles(self):
         except ImportError:
             print("Slycot not present, skipping")
 
+    @unittest.skipIf(not slycot_check(), "slycot not installed")
+    def test_tf2ss_robustness(self):
+        """Unit test to make sure that tf2ss is working correctly.
+         Source: https://github.com/python-control/python-control/issues/240
+        """
+        import control
+        
+        num =  [ [[0], [1]],           [[1],   [0]] ]
+        den1 = [ [[1], [1,1]],         [[1,4], [1]] ]
+        sys1tf = control.tf(num, den1)
+        sys1ss = control.tf2ss(sys1tf)
+
+        # slight perturbation
+        den2 = [ [[1], [1e-10, 1, 1]], [[1,4], [1]] ]
+        sys2tf = control.tf(num, den2)
+        sys2ss = control.tf2ss(sys2tf)
+
+        # Make sure that the poles match for StateSpace and TransferFunction
+        np.testing.assert_array_almost_equal(np.sort(sys1tf.pole()),
+                                             np.sort(sys1ss.pole()))
+        np.testing.assert_array_almost_equal(np.sort(sys2tf.pole()),
+                                             np.sort(sys2ss.pole()))
+
 def suite():
    return unittest.TestLoader().loadTestsFromTestCase(TestConvert)
 
diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py
index a6fc49a7f..609a91e8f 100644
--- a/control/tests/freqresp_test.py
+++ b/control/tests/freqresp_test.py
@@ -141,6 +141,11 @@ def test_discrete(self):
          import warnings
          warnings.simplefilter('always', UserWarning)   # don't supress
          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")
+
+            # 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)
             print("len(w) =", len(w))
diff --git a/control/tests/matlab_test.py b/control/tests/matlab_test.py
index efde21c1d..d187f6125 100644
--- a/control/tests/matlab_test.py
+++ b/control/tests/matlab_test.py
@@ -414,6 +414,14 @@ def testAcker(self):
     @unittest.skipIf(not slycot_check(), "slycot not installed")
     def testLQR(self):
         (K, S, E) = lqr(self.siso_ss1.A, self.siso_ss1.B, np.eye(2), np.eye(1))
+
+        # Should work if [Q N;N' R] is positive semi-definite
+        (K, S, E) = lqr(self.siso_ss2.A, self.siso_ss2.B, 10*np.eye(3), \
+                            np.eye(1), [[1], [1], [2]])
+
+    @unittest.skip("check not yet implemented")
+    def testLQR_checks(self):
+        # Make sure we get a warning if [Q N;N' R] is not positive semi-definite
         (K, S, E) = lqr(self.siso_ss2.A, self.siso_ss2.B, np.eye(3), \
                             np.eye(1), [[1], [1], [2]])
 
diff --git a/control/tests/statefbk_test.py b/control/tests/statefbk_test.py
index 042bda701..35df769a2 100644
--- a/control/tests/statefbk_test.py
+++ b/control/tests/statefbk_test.py
@@ -6,7 +6,7 @@
 from __future__ import print_function
 import unittest
 import numpy as np
-from control.statefbk import ctrb, obsv, place, lqr, gram, acker
+from control.statefbk import ctrb, obsv, place, place_varga, lqr, gram, acker
 from control.matlab import *
 from control.exception import slycot_check, ControlDimension
 
@@ -186,7 +186,10 @@ def testPlace(self):
         np.testing.assert_raises(ValueError, place, A, B, P_repeated)
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
-    def testPlace_varga(self):
+    def testPlace_varga_continuous(self):
+        """
+        Check that we can place eigenvalues for dtime=False
+        """
         A = np.array([[1., -2.], [3., -4.]])
         B = np.array([[5.], [7.]])
 
@@ -202,6 +205,77 @@ def testPlace_varga(self):
         np.testing.assert_raises(ControlDimension, place, A[1:, :], B, P)
         np.testing.assert_raises(ControlDimension, place, A, B[1:, :], P)
 
+        # Regression test against bug #177
+        # https://github.com/python-control/python-control/issues/177
+        A = np.array([[0, 1], [100, 0]])
+        B = np.array([[0], [1]])
+        P = np.array([-20 + 10*1j, -20 - 10*1j])
+        K = place_varga(A, B, P)
+        P_placed = np.linalg.eigvals(A - B.dot(K))
+
+        # No guarantee of the ordering, so sort them
+        P.sort()
+        P_placed.sort()
+        np.testing.assert_array_almost_equal(P, P_placed)
+
+    @unittest.skipIf(not slycot_check(), "slycot not installed")
+    def testPlace_varga_continuous_partial_eigs(self):
+        """
+        Check that we are able to use the alpha parameter to only place
+        a subset of the eigenvalues, for the continous time case.
+        """
+        # A matrix has eigenvalues at s=-1, and s=-2. Choose alpha = -1.5
+        # and check that eigenvalue at s=-2 stays put.
+        A = np.array([[1., -2.], [3., -4.]])
+        B = np.array([[5.], [7.]])
+
+        P = np.array([-3.])
+        P_expected = np.array([-2.0, -3.0])
+        alpha = -1.5
+        K = place_varga(A, B, P, alpha=alpha)
+
+        P_placed = np.linalg.eigvals(A - B.dot(K))
+        # No guarantee of the ordering, so sort them
+        P_expected.sort()
+        P_placed.sort()
+        np.testing.assert_array_almost_equal(P_expected, P_placed)
+
+    @unittest.skipIf(not slycot_check(), "slycot not installed")
+    def testPlace_varga_discrete(self):
+        """
+        Check that we can place poles using dtime=True (discrete time)
+        """
+        A = np.array([[1., 0], [0, 0.5]])
+        B = np.array([[5.], [7.]])
+
+        P = np.array([0.5, 0.5])
+        K = place_varga(A, B, P, dtime=True)
+        P_placed = np.linalg.eigvals(A - B.dot(K))
+        # No guarantee of the ordering, so sort them
+        P.sort()
+        P_placed.sort()
+        np.testing.assert_array_almost_equal(P, P_placed)
+
+    @unittest.skipIf(not slycot_check(), "slycot not installed")
+    def testPlace_varga_discrete_partial_eigs(self):
+        """"
+        Check that we can only assign a single eigenvalue in the discrete
+        time case.
+        """
+        # A matrix has eigenvalues at 1.0 and 0.5. Set alpha = 0.51, and
+        # check that the eigenvalue at 0.5 is not moved.
+        A = np.array([[1., 0], [0, 0.5]])
+        B = np.array([[5.], [7.]])
+        P = np.array([0.2, 0.6])
+        P_expected = np.array([0.5, 0.6])
+        alpha = 0.51
+        K = place_varga(A, B, P, dtime=True, alpha=alpha)
+        P_placed = np.linalg.eigvals(A - B.dot(K))
+        P_expected.sort()
+        P_placed.sort()
+        np.testing.assert_array_almost_equal(P_expected, P_placed)
+
+
     def check_LQR(self, K, S, poles, Q, R):
         S_expected = np.array(np.sqrt(Q * R))
         K_expected = S_expected / R
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 35f511ab1..3e9adcca6 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -148,7 +148,11 @@ def testEvalFr(self):
         # Deprecated version of the call (should generate warning)
         import warnings
         with warnings.catch_warnings(record=True) as w:
-            warnings.simplefilter("always")
+            # 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)
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index 9a16e26a1..9e11a9d0c 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -15,6 +15,7 @@
 from control.statesp import *
 from control.xferfcn import TransferFunction, _convertToTransferFunction
 from control.dtime import c2d
+from control.exception import slycot_check
 
 class TestTimeresp(unittest.TestCase):
     def setUp(self):
@@ -233,6 +234,23 @@ def test_discrete_initial(self):
         t, yout = impulse_response(h1, np.arange(4))
         np.testing.assert_array_equal(yout[0], [0., 1., 0., 0.])
 
+    @unittest.skipIf(not slycot_check(), "slycot not installed")
+    def test_step_robustness(self):
+        "Unit test: https://github.com/python-control/python-control/issues/240"
+        # Create 2 input, 2 output system
+        num =  [ [[0], [1]],           [[1],   [0]] ]
+        
+        den1 = [ [[1], [1,1]],         [[1,4], [1]] ]
+        sys1 = TransferFunction(num, den1)
+
+        den2 = [ [[1], [1e-10, 1, 1]], [[1,4], [1]] ]   # slight perturbation
+        sys2 = TransferFunction(num, den2)
+
+        # Compute step response from input 1 to output 1, 2
+        t1, y1 = step_response(sys1, input=0)
+        t2, y2 = step_response(sys2, input=0)
+        np.testing.assert_array_almost_equal(y1, y2)
+
 def suite():
     return unittest.TestLoader().loadTestsFromTestCase(TestTimeresp)
 
diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py
index 204c6dfd8..44953e731 100644
--- a/control/tests/xferfcn_test.py
+++ b/control/tests/xferfcn_test.py
@@ -5,8 +5,8 @@
 
 import unittest
 import numpy as np
-from control.statesp import StateSpace, _convertToStateSpace
-from control.xferfcn import TransferFunction, _convertToTransferFunction
+from control.statesp import StateSpace, _convertToStateSpace, rss
+from control.xferfcn import TransferFunction, _convertToTransferFunction, ss2tf
 from control.lti import evalfr
 from control.exception import slycot_check
 # from control.lti import isdtime
@@ -536,6 +536,26 @@ def testMIMO(self):
         np.testing.assert_array_almost_equal(H.num[1][0], H2.num[1][0])
         np.testing.assert_array_almost_equal(H.den[1][0], H2.den[1][0])
 
+    @unittest.skipIf(not slycot_check(), "slycot not installed")
+    def testIndexing(self):
+        tm = ss2tf(rss(5, 3, 3))
+
+        # scalar indexing
+        sys01 = tm[0, 1]
+        np.testing.assert_array_almost_equal(sys01.num[0][0], tm.num[0][1])
+        np.testing.assert_array_almost_equal(sys01.den[0][0], tm.den[0][1])
+
+        # slice indexing
+        sys = tm[:2, 1:3]
+        np.testing.assert_array_almost_equal(sys.num[0][0], tm.num[0][1])
+        np.testing.assert_array_almost_equal(sys.den[0][0], tm.den[0][1])
+        np.testing.assert_array_almost_equal(sys.num[0][1], tm.num[0][2])
+        np.testing.assert_array_almost_equal(sys.den[0][1], tm.den[0][2])
+        np.testing.assert_array_almost_equal(sys.num[1][0], tm.num[1][1])
+        np.testing.assert_array_almost_equal(sys.den[1][0], tm.den[1][1])
+        np.testing.assert_array_almost_equal(sys.num[1][1], tm.num[1][2])
+        np.testing.assert_array_almost_equal(sys.den[1][1], tm.den[1][2])
+
     def testMatrixMult(self):
         """MIMO transfer functions should be multiplyable by constant
         matrices"""
diff --git a/control/timeresp.py b/control/timeresp.py
index 4c137c933..a0ecf9d93 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -592,7 +592,6 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None,
     else:
         new_X0 = X0
         U[0] = 1.
-        print("discrete", U)
 
     T, yout, _xout = forced_response(
         sys, T, U, new_X0,
diff --git a/control/xferfcn.py b/control/xferfcn.py
index 50b9877fa..9e272342c 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -491,6 +491,46 @@ def __pow__(self, other):
         if other < 0:
             return (TransferFunction([1], [1]) / self) * (self**(other+1))
 
+    def __getitem__(self, key):
+        key1, key2 = key
+
+        # pre-process
+        if isinstance(key1, int):
+            key1 = slice(key1, key1 + 1, 1)
+        if isinstance(key2, int):
+            key2 = slice(key2, key2 + 1, 1)
+        # dim1
+        start1, stop1, step1 = key1.start, key1.stop, key1.step
+        if step1 is None:
+            step1 = 1
+        if start1 is None:
+            start1 = 0
+        if stop1 is None:
+            stop1 = len(self.num)
+        # dim1
+        start2, stop2, step2 = key2.start, key2.stop, key2.step
+        if step2 is None:
+            step2 = 1
+        if start2 is None:
+            start2 = 0
+        if stop2 is None:
+            stop2 = len(self.num[0])
+
+        num = []
+        den = []
+        for i in range(start1, stop1, step1):
+            num_i = []
+            den_i = []
+            for j in range(start2, stop2, step2):
+                num_i.append(self.num[i][j])
+                den_i.append(self.den[i][j])
+            num.append(num_i)
+            den.append(den_i)
+        if self.isctime():
+            return TransferFunction(num, den)
+        else:
+            return TransferFunction(num, den, self.dt)
+
     def evalfr(self, omega):
         """Evaluate a transfer function at a single angular frequency.
 
@@ -510,7 +550,7 @@ def _evalfr(self, omega):
             # Convert the frequency to discrete time
             dt = timebase(self)
             s = exp(1.j * omega * dt)
-            if (omega * dt > pi):
+            if np.any(omega * dt > pi):
                 warn("_evalfr: frequency evaluation above Nyquist frequency")
         else:
             s = 1.j * omega
@@ -800,9 +840,13 @@ def _common_den(self, imag_tol=None):
                     num[i,j,0] = poleset[i][j][2]
             else:
                 # create the denominator matching this input
+                # polyfromroots gives coeffs in opposite order from what we use
+                # coefficients should be padded on right, ending at np
                 np = len(poles[j])
                 den[j,np::-1] = polyfromroots(poles[j]).real
                 denorder[j] = np
+
+                # now create the numerator, also padded on the right
                 for i in range(self.outputs):
                     # start with the current set of zeros for this output
                     nwzeros = list(poleset[i][j][0])
@@ -811,14 +855,15 @@ def _common_den(self, imag_tol=None):
                     for ip in chain(poleset[i][j][3],
                                     range(poleset[i][j][4],np)):
                         nwzeros.append(poles[j][ip])
-                    
+
                     numpoly = poleset[i][j][2] * polyfromroots(nwzeros).real 
-                    m = npmax - len(numpoly)
-                    #print(j,i,m,len(numpoly),len(poles[j]))
-                    if m < 0:
-                        num[i,j,::-1] = numpoly
-                    else:
-                        num[i,j,:m:-1] = numpoly   
+                    # print(numpoly, den[j])
+                    # polyfromroots gives coeffs in opposite order => invert
+                    # numerator polynomial should be padded on left and right
+                    #   ending at np to line up with what td04ad expects...
+                    num[i, j, np+1-len(numpoly):np+1] = numpoly[::-1]
+                    # print(num[i, j])
+
         if (abs(den.imag) > epsnm).any():
             print("Warning: The denominator has a nontrivial imaginary part: %f"
                       % abs(den.imag).max())
@@ -1017,7 +1062,7 @@ def _convertToTransferFunction(sys, **kw):
     If sys is an array-like type, then it is converted to a constant-gain
     transfer function.
 
-    >>> sys = _convertToTransferFunction([[1. 0.], [2. 3.]])
+    >>> sys = _convertToTransferFunction([[1., 0.], [2., 3.]])
 
     In this example, the numerator matrix will be
        [[[1.0], [0.0]], [[2.0], [3.0]]]
diff --git a/doc/README b/doc/README
index 9b829033e..f542badc5 100644
--- a/doc/README
+++ b/doc/README
@@ -1,7 +1,11 @@
 Sphinx Documentation
 --------------------
+This directory contains the user manual for the python-control
+toolbox.  The documentation is built using sphinx.  The master toctree
+document is index.rst.
 
-Note: Sphinx actually runs and imports python code, so broken code, or code not in conf.py sys.path, cannot be documented!
+Note: Sphinx actually runs and imports python code, so broken code, or
+code not in conf.py sys.path, cannot be documented!
 
 1.  Get Sphinx [http://sphinx.pocoo.org/]
     [python setup.py build/install]
@@ -14,7 +18,8 @@ Note: Sphinx actually runs and imports python code, so broken code, or code not
 4.  >> touch *.rst
     >> make html  [or make latex]
 
-Creating/updating manual on sourceforge:
-
-5.  >> rsync -rav _build/html/ user@shell.sourceforge.net:/home/project-web/python-control/htdocs/manual-N.mx/
+Creating/updating manual on readthedocs.org:
 
+5.  Log in to readthedocs.org and go to the 'Admin' menu for
+python-control.  Choose 'Versions' from the sidebar and mark the
+latest release as 'Active'.  Update the default version if needed.
diff --git a/doc/control.rst b/doc/control.rst
index 40cb407bb..1d0b14644 100644
--- a/doc/control.rst
+++ b/doc/control.rst
@@ -154,3 +154,5 @@ Utility functions and conversions
     timebase
     timebaseEqual
     unwrap
+    use_fbs_defaults
+    use_matlab_defaults
diff --git a/doc/conventions.rst b/doc/conventions.rst
index 048356436..7bdf3c628 100644
--- a/doc/conventions.rst
+++ b/doc/conventions.rst
@@ -15,7 +15,7 @@ LTI system representation
 Linear time invariant (LTI) systems are represented in python-control in
 state space, transfer function, or frequency response data (FRD) form.  Most
 functions in the toolbox will operate on any of these data types and
-functions for converting between between compatible types is provided.
+functions for converting between compatible types is provided.
 
 State space systems
 -------------------
@@ -45,8 +45,8 @@ transfer functions
 .. math::
 
   G(s) = \frac{\text{num}(s)}{\text{den}(s)}
-       = \frac{a_0 s^n + a_1 s^{n-1} + \cdots + a_n}
-              {b_0 s^m + b_1 s^{m-1} + \cdots + b_m},
+       = \frac{a_0 s^m + a_1 s^{m-1} + \cdots + a_m}
+              {b_0 s^n + b_1 s^{n-1} + \cdots + b_n},
 
 where n is generally greater than or equal to m (for a proper transfer
 function).
@@ -77,11 +77,10 @@ performed.
 
 Discrete time systems
 ---------------------
-By default, all systems are considered to be continuous time systems.  A
-discrete time system is created by specifying the 'time base' dt.  The time
-base argument can be given when a system is constructed:
+A discrete time system is created by specifying a nonzero 'timebase', dt.
+The timebase argument can be given when a system is constructed:
 
-* dt = None: no timebase specified
+* dt = None: no timebase specified (default)
 * dt = 0: continuous time system
 * dt > 0: discrete time system with sampling period 'dt'
 * dt = True: discrete time with unspecified sampling period
@@ -89,11 +88,16 @@ base argument can be given when a system is constructed:
 Only the :class:`StateSpace` and :class:`TransferFunction` classes allow
 explicit representation of discrete time systems.
 
-Systems must have the same time base in order to be combined.  For
-continuous time systems, the :func:`sample_system` function or the
-:meth:`StateSpace.sample` and :meth:`TransferFunction.sample` methods can be
-used to create a discrete time system from a continuous time system.  See
-:ref:`utility-and-conversions`.
+Systems must have compatible timebases in order to be combined.  A system
+with timebase `None` can be combined with a system having a specified
+timebase; the result will have the timebase of the latter system.
+Similarly, a discrete time system with unspecified sampling time (`dt =
+True`) can be combined with a system having a specified sampling time;
+the result will be a discrete time system with the sample time of the latter
+system.  For continuous time systems, the :func:`sample_system` function or
+the :meth:`StateSpace.sample` and :meth:`TransferFunction.sample` methods
+can be used to create a discrete time system from a continuous time system.
+See :ref:`utility-and-conversions`.
 
 Conversion between representations
 ----------------------------------
@@ -106,18 +110,23 @@ argument or using the explicit conversion functions :func:`ss2tf` and
 
 Time series data
 ================
-
-This is a convention for function arguments and return values that
-represent time series: sequences of values that change over time. It
-is used throughout the library, for example in the functions
+A variety of functions in the library return time series data: sequences of
+values that change over time.  A common set of conventions is used for
+returning such data: columns represent different points in time, rows are
+different components (e.g., inputs, outputs or states).  For return
+arguments, an array of times is given as the first returned argument,
+followed by one or more arrays of variable values.  This convention is used
+throughout the library, for example in the functions
 :func:`forced_response`, :func:`step_response`, :func:`impulse_response`,
 and :func:`initial_response`.
 
 .. note::
-    This convention is different from the convention used in the library
-    :mod:`scipy.signal`. In Scipy's convention the meaning of rows and columns
-    is interchanged.  Thus, all 2D values must be transposed when they are
-    used with functions from :mod:`scipy.signal`.
+    The convention used by python-control is different from the convention
+    used in the `scipy.signal
+    <https://docs.scipy.org/doc/scipy/reference/signal.html>`_ library. In
+    Scipy's convention the meaning of rows and columns is interchanged.
+    Thus, all 2D values must be transposed when they are used with functions
+    from `scipy.signal`_.
 
 Types:
 
@@ -177,7 +186,7 @@ conventions.  The currently configurable options allow the units for Bode
 plots to be set as dB for gain, degrees for phase and Hertz for frequency
 (MATLAB conventions) or the gain can be given in magnitude units (powers of
 10), corresponding to the conventions used in `Feedback Systems
-<http://www.cds.caltech.edu/~murray/FBSwiki>`_.
+<http://fbsbook.org>`_ (FBS).
 
 Variables that can be configured, along with their default values:
   * bode_dB (False): Bode plot magnitude plotted in dB (otherwise powers of 10)
@@ -190,6 +199,7 @@ Variables that can be configured, along with their default values:
 Functions that can be used to set standard configurations:
 
 .. autosummary::
-
+    :toctree: generated/
+    
     use_fbs_defaults
     use_matlab_defaults
diff --git a/doc/index.rst b/doc/index.rst
index 407e89817..7ea8fe1dd 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -12,7 +12,8 @@ implements basic operations for analysis and design of feedback control systems.
 - Time response: initial, step, impulse
 - Frequency response: Bode and Nyquist plots
 - Control analysis: stability, reachability, observability, stability margins
-- Control design: eigenvalue placement, linear quadratic regulator
+- Control design: eigenvalue placement, LQR, H2, Hinf
+- Model reduction: balanced realizations, Hankel singular values
 - Estimator design: linear quadratic estimator (Kalman filter)
 
 .. rubric:: Documentation
diff --git a/doc/intro.rst b/doc/intro.rst
index 7402576ba..9677135c1 100644
--- a/doc/intro.rst
+++ b/doc/intro.rst
@@ -7,36 +7,41 @@ Manual.  This manual contains information on using the python-control
 package, including documentation for all functions in the package and
 examples illustrating their use.
 
-Overview of the Toolbox
+Overview of the toolbox
 =======================
 
 The python-control package is a set of python classes and functions that
 implement common operations for the analysis and design of feedback control
 systems.  The initial goal is to implement all of the functionality required
 to work through the examples in the textbook `Feedback Systems
-<http://www.cds.caltech.edu/~murray/FBSwiki>`_ by Astrom and Murray. A
-MATLAB compatibility package (control.matlab) is available that provides
-many of the common functions corresponding to commands available in the
-MATLAB Control Systems Toolbox.
+<http://fbsbook.org>`_ by Astrom and Murray. A :ref:`matlab-module` is
+available that provides many of the common functions corresponding to
+commands available in the MATLAB Control Systems Toolbox.
 
-Some Differences from MATLAB
+Some differences from MATLAB
 ============================
-The python-control package makes use of NumPy and SciPy.  A list of general
-differences between NumPy and MATLAB can be found `here
+The python-control package makes use of `NumPy <http://www.numpy.org>`_ and
+`SciPy <https://www.scipy.org>`_.  A list of general differences between
+NumPy and MATLAB can be found `here
 <http://www.scipy.org/NumPy_for_Matlab_Users>`_.
 
 In terms of the python-control package more specifically, here are
 some thing to keep in mind:
 
 * You must include commas in vectors.  So [1 2 3] must be [1, 2, 3].
-* Functions that return multiple arguments use tuples
-* You cannot use braces for collections; use tuples instead
+* Functions that return multiple arguments use tuples.  
+* You cannot use braces for collections; use tuples instead.
 
 Installation
 ============
 
-The `python-control` package may be installed using pip, conda or the
-standard distutils/setuptools mechanisms.
+The `python-control` package can be installed using pip, conda or the
+standard distutils/setuptools mechanisms.  The package requires `numpy`_ and
+`scipy`_, and the plotting routines require `matplotlib
+<https://matplotlib.org>`_.  In addition, some routines require the `slycot
+<https://github.com/python-control/Slycot>`_ library in order to implement
+more advanced features (including some MIMO functionality).
+
 
 To install using pip::
 
@@ -54,22 +59,23 @@ correctly by running the command::
   python -c "import slycot"
 
 and verifying that no error message appears.  It may be necessary to install
-`slycot` from source, which requires a working FORTRAN compiler and the
-`lapack` library.  More information on the slycot package can be obtained
-from the `slycot project page <https://github.com/python-control/Slycot>`_.
+`slycot` from source, which requires a working FORTRAN compiler and either
+the `lapack` or `openplas` library.  More information on the slycot package
+can be obtained from the `slycot project page
+<https://github.com/python-control/Slycot>`_.
 
 For users with the Anaconda distribution of Python, the following
 commands can be used::
 
   conda install numpy scipy matplotlib    # if not yet installed
-  conda install -c python-control -c cyclus slycot control
+  conda install -c conda-forge control
 
-This installs `slycot` and `python-control` from the `python-control`
-channel and uses the `cyclus` channel to obtain the required `lapack`
-package. 
+This installs `slycot` and `python-control` from conda-forge, including the
+`openblas` package.
 
-Alternatively, to use setuptools, first `download the source <https://github.com/python-control/python-control/releases>`_ and unpack
-it.  To install in your home directory, use::
+Alternatively, to use setuptools, first `download the source
+<https://github.com/python-control/python-control/releases>`_ and unpack it.
+To install in your home directory, use::
 
   python setup.py install --user
 
@@ -78,11 +84,7 @@ or to install for all users (on Linux or Mac OS)::
   python setup.py build
   sudo python setup.py install
 
-The package requires `numpy` and `scipy`, and the plotting routines require
-`matplotlib`.  In addition, some routines require the `slycot` module,
-described above.
-
-Getting Started
+Getting started
 ===============
 
 There are two different ways to use the package.  For the default interface
diff --git a/setup.py b/setup.py
index 6751f64cd..cd4bcbf9f 100644
--- a/setup.py
+++ b/setup.py
@@ -19,9 +19,8 @@
 Programming Language :: Python :: 2
 Programming Language :: Python :: 2.7
 Programming Language :: Python :: 3
-Programming Language :: Python :: 3.3
-Programming Language :: Python :: 3.4
 Programming Language :: Python :: 3.5
+Programming Language :: Python :: 3.6
 Topic :: Software Development
 Topic :: Scientific/Engineering
 Operating System :: Microsoft :: Windows