diff --git a/control/config.py b/control/config.py
index 27a5712a3..8703af417 100644
--- a/control/config.py
+++ b/control/config.py
@@ -15,7 +15,7 @@
 
 # Package level default values
 _control_defaults = {
-    # No package level defaults (yet)
+    'control.default_dt':0
 }
 defaults = dict(_control_defaults)
 
@@ -59,6 +59,9 @@ def reset_defaults():
     from .statesp import _statesp_defaults
     defaults.update(_statesp_defaults)
 
+    from .iosys import _iosys_defaults
+    defaults.update(_iosys_defaults)
+
 
 def _get_param(module, param, argval=None, defval=None, pop=False):
     """Return the default value for a configuration option.
@@ -170,5 +173,6 @@ def use_legacy_defaults(version):
     """
     if version == '0.8.3': 
         use_numpy_matrix(True) # alternatively: set_defaults('statesp', use_numpy_matrix=True)
+        set_defaults('control', default_dt=None)
     else:
         raise ValueError('''version number not recognized. Possible values are: ['0.8.3']''') 
\ No newline at end of file
diff --git a/control/dtime.py b/control/dtime.py
index 89f17c4af..725bcde47 100644
--- a/control/dtime.py
+++ b/control/dtime.py
@@ -53,21 +53,19 @@
 
 # Sample a continuous time system
 def sample_system(sysc, Ts, method='zoh', alpha=None, prewarp_frequency=None):
-    """Convert a continuous time system to discrete time
-
-    Creates a discrete time system from a continuous time system by
-    sampling.  Multiple methods of conversion are supported.
+    """
+    Convert a continuous time system to discrete time by sampling
 
     Parameters
     ----------
-    sysc : linsys
+    sysc : LTI (StateSpace or TransferFunction)
         Continuous time system to be converted
-    Ts : real
+    Ts : real > 0
         Sampling period
     method : string
-        Method to use for conversion: 'matched', 'tustin', 'zoh' (default)
+        Method to use for conversion, e.g. 'bilinear', 'zoh' (default)
 
-    prewarp_frequency : float within [0, infinity)
+    prewarp_frequency : real within [0, infinity)
         The frequency [rad/s] at which to match with the input continuous-
         time system's magnitude and phase
 
@@ -78,13 +76,13 @@ def sample_system(sysc, Ts, method='zoh', alpha=None, prewarp_frequency=None):
 
     Notes
     -----
-    See `TransferFunction.sample` and `StateSpace.sample` for
+    See :meth:`StateSpace.sample` or :meth:`TransferFunction.sample`` for
     further details.
 
     Examples
     --------
     >>> sysc = TransferFunction([1], [1, 2, 1])
-    >>> sysd = sample_system(sysc, 1, method='matched')
+    >>> sysd = sample_system(sysc, 1, method='bilinear')
     """
 
     # Make sure we have a continuous time system
@@ -95,35 +93,39 @@ def sample_system(sysc, Ts, method='zoh', alpha=None, prewarp_frequency=None):
 
 
 def c2d(sysc, Ts, method='zoh', prewarp_frequency=None):
-    '''
-    Return a discrete-time system
+    """
+    Convert a continuous time system to discrete time by sampling
 
     Parameters
     ----------
-    sysc: LTI (StateSpace or TransferFunction), continuous
-        System to be converted
+    sysc : LTI (StateSpace or TransferFunction)
+        Continuous time system to be converted
+    Ts : real > 0
+        Sampling period
+    method : string
+        Method to use for conversion, e.g. 'bilinear', 'zoh' (default)
 
-    Ts: number
-        Sample time for the conversion
+    prewarp_frequency : real within [0, infinity)
+        The frequency [rad/s] at which to match with the input continuous-
+        time system's magnitude and phase
 
-    method: string, optional
-        Method to be applied,
-        'zoh'        Zero-order hold on the inputs (default)
-        'foh'        First-order hold, currently not implemented
-        'impulse'    Impulse-invariant discretization, currently not implemented
-        'tustin'     Bilinear (Tustin) approximation, only SISO
-        'matched'    Matched pole-zero method, only SISO
+    Returns
+    -------
+    sysd : linsys
+        Discrete time system, with sampling rate Ts
+
+    Notes
+    -----
+    See :meth:`StateSpace.sample` or :meth:`TransferFunction.sample`` for
+    further details.
 
-    prewarp_frequency : float within [0, infinity)
-            The frequency [rad/s] at which to match with the input continuous-
-            time system's magnitude and phase
+    Examples
+    --------
+    >>> sysc = TransferFunction([1], [1, 2, 1])
+    >>> sysd = sample_system(sysc, 1, method='bilinear')
+    """
 
-    '''
     #  Call the sample_system() function to do the work
     sysd = sample_system(sysc, Ts, method, prewarp_frequency)
 
-    # TODO: is this check needed?  If sysc is  StateSpace, sysd is too?
-    if isinstance(sysc, StateSpace) and not isinstance(sysd, StateSpace):
-        return _convertToStateSpace(sysd)       # pragma: no cover
-
     return sysd
diff --git a/control/iosys.py b/control/iosys.py
index 1fe10346f..528f1a347 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -38,12 +38,15 @@
 
 from .statesp import StateSpace, tf2ss
 from .timeresp import _check_convert_array
-from .lti import isctime, isdtime, _find_timebase
+from .lti import isctime, isdtime, common_timebase
+from . import config
 
 __all__ = ['InputOutputSystem', 'LinearIOSystem', 'NonlinearIOSystem',
            'InterconnectedSystem', 'input_output_response', 'find_eqpt',
            'linearize', 'ss2io', 'tf2io']
 
+# Define module default parameter values
+_iosys_defaults = {}
 
 class InputOutputSystem(object):
     """A class for representing input/output systems.
@@ -69,9 +72,11 @@ class for a set of subclasses that are used to implement specific
     states : int, list of str, or None
         Description of the system states.  Same format as `inputs`.
     dt : None, True or float, optional
-        System timebase.  None (default) indicates continuous time, True
-        indicates discrete time with undefined sampling time, positive number
-        is discrete time with specified sampling time.
+        System timebase. 0 (default) indicates continuous
+        time, True indicates discrete time with unspecified sampling
+        time, positive number is discrete time with specified
+        sampling time, None indicates unspecified timebase (either  
+        continuous or discrete time).
     params : dict, optional
         Parameter values for the systems.  Passed to the evaluation functions
         for the system as default values, overriding internal defaults.
@@ -86,9 +91,11 @@ class for a set of subclasses that are used to implement specific
         Dictionary of signal names for the inputs, outputs and states and the
         index of the corresponding array
     dt : None, True or float
-        System timebase.  None (default) indicates continuous time, True
-        indicates discrete time with undefined sampling time, positive number
-        is discrete time with specified sampling time.
+        System timebase. 0 (default) indicates continuous
+        time, True indicates discrete time with unspecified sampling
+        time, positive number is discrete time with specified
+        sampling time, None indicates unspecified timebase (either  
+        continuous or discrete time).
     params : dict, optional
         Parameter values for the systems.  Passed to the evaluation functions
         for the system as default values, overriding internal defaults.
@@ -109,7 +116,7 @@ class for a set of subclasses that are used to implement specific
 
     """
     def __init__(self, inputs=None, outputs=None, states=None, params={},
-                 dt=None, name=None):
+                 name=None, **kwargs):
         """Create an input/output system.
 
         The InputOutputSystem contructor is used to create an input/output
@@ -134,10 +141,11 @@ def __init__(self, inputs=None, outputs=None, states=None, params={},
         states : int, list of str, or None
             Description of the system states.  Same format as `inputs`.
         dt : None, True or float, optional
-            System timebase.  None (default) indicates continuous
-            time, True indicates discrete time with undefined sampling
+            System timebase. 0 (default) indicates continuous
+            time, True indicates discrete time with unspecified sampling
             time, positive number is discrete time with specified
-            sampling time.
+            sampling time, None indicates unspecified timebase (either  
+            continuous or discrete time).
         params : dict, optional
             Parameter values for the systems.  Passed to the evaluation
             functions for the system as default values, overriding internal
@@ -153,7 +161,7 @@ def __init__(self, inputs=None, outputs=None, states=None, params={},
         """
         # Store the input arguments
         self.params = params.copy()     # default parameters
-        self.dt = dt                    # timebase
+        self.dt = kwargs.get('dt', config.defaults['control.default_dt'])                    # timebase
         self.name = name                # system name
 
         # Parse and store the number of inputs, outputs, and states
@@ -200,10 +208,8 @@ def __mul__(sys2, sys1):
                              "inputs and outputs")
 
         # Make sure timebase are compatible
-        dt = _find_timebase(sys1, sys2)
-        if dt is False:
-            raise ValueError("System timebases are not compabile")
-
+        dt = common_timebase(sys1.dt, sys2.dt)
+        
         # Return the series interconnection between the systems
         newsys = InterconnectedSystem((sys1, sys2))
 
@@ -478,10 +484,8 @@ def feedback(self, other=1, sign=-1, params={}):
                              "inputs and outputs")
 
         # Make sure timebases are compatible
-        dt = _find_timebase(self, other)
-        if dt is False:
-            raise ValueError("System timebases are not compabile")
-
+        dt = common_timebase(self.dt, other.dt)
+        
         # Return the series interconnection between the systems
         newsys = InterconnectedSystem((self, other), params=params, dt=dt)
 
@@ -601,10 +605,11 @@ def __init__(self, linsys, inputs=None, outputs=None, states=None,
         states : int, list of str, or None, optional
             Description of the system states.  Same format as `inputs`.
         dt : None, True or float, optional
-            System timebase.  None (default) indicates continuous
-            time, True indicates discrete time with undefined sampling
+            System timebase. 0 (default) indicates continuous
+            time, True indicates discrete time with unspecified sampling
             time, positive number is discrete time with specified
-            sampling time.
+            sampling time, None indicates unspecified timebase (either  
+            continuous or discrete time).
         params : dict, optional
             Parameter values for the systems.  Passed to the evaluation
             functions for the system as default values, overriding internal
@@ -670,7 +675,8 @@ class NonlinearIOSystem(InputOutputSystem):
 
     """
     def __init__(self, updfcn, outfcn=None, inputs=None, outputs=None,
-                 states=None, params={}, dt=None, name=None):
+                 states=None, params={},  
+                 name=None, **kwargs):
         """Create a nonlinear I/O system given update and output functions.
 
         Creates an `InputOutputSystem` for a nonlinear system by specifying a
@@ -722,10 +728,10 @@ def __init__(self, updfcn, outfcn=None, inputs=None, outputs=None,
             operating in continuous or discrete time.  It can have the
             following values:
 
-            * dt = None       No timebase specified
-            * dt = 0          Continuous time system
-            * dt > 0          Discrete time system with sampling time dt
-            * dt = True       Discrete time with unspecified sampling time
+            * dt = 0: continuous time system (default)
+            * dt > 0: discrete time system with sampling period 'dt'
+            * dt = True: discrete time with unspecified sampling period
+            * dt = None: no timebase specified 
 
         name : string, optional
             System name (used for specifying signals).
@@ -741,6 +747,7 @@ def __init__(self, updfcn, outfcn=None, inputs=None, outputs=None,
         self.outfcn = outfcn
 
         # Initialize the rest of the structure
+        dt = kwargs.get('dt', config.defaults['control.default_dt'])
         super(NonlinearIOSystem, self).__init__(
             inputs=inputs, outputs=outputs, states=states,
             params=params, dt=dt, name=name
@@ -865,10 +872,10 @@ def __init__(self, syslist, connections=[], inplist=[], outlist=[],
             operating in continuous or discrete time.  It can have the
             following values:
 
-            * dt = None       No timebase specified
-            * dt = 0          Continuous time system
-            * dt > 0          Discrete time system with sampling time dt
-            * dt = True       Discrete time with unspecified sampling time
+            * dt = 0: continuous time system (default)
+            * dt > 0: discrete time system with sampling period 'dt'
+            * dt = True: discrete time with unspecified sampling period
+            * dt = None: no timebase specified 
 
         name : string, optional
             System name (used for specifying signals).
@@ -881,19 +888,14 @@ def __init__(self, syslist, connections=[], inplist=[], outlist=[],
         # Check to make sure all systems are consistent
         self.syslist = syslist
         self.syslist_index = {}
-        dt = None
         nstates = 0; self.state_offset = []
         ninputs = 0; self.input_offset = []
         noutputs = 0; self.output_offset = []
         system_count = 0
+
         for sys in syslist:
             # Make sure time bases are consistent
-            # TODO: Use lti._find_timebase() instead?
-            if dt is None and sys.dt is not None:
-                # Timebase was not specified; set to match this system
-                dt = sys.dt
-            elif dt != sys.dt:
-                raise TypeError("System timebases are not compatible")
+            dt = common_timebase(dt, sys.dt)
 
             # Make sure number of inputs, outputs, states is given
             if sys.ninputs is None or sys.noutputs is None or \
diff --git a/control/lti.py b/control/lti.py
index 8db14794b..e41fe416b 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -9,14 +9,16 @@
 isdtime()
 isctime()
 timebase()
-timebaseEqual()
+common_timebase()
 """
 
 import numpy as np
 from numpy import absolute, real
+from warnings import warn
 
-__all__ = ['issiso', 'timebase', 'timebaseEqual', 'isdtime', 'isctime',
-           'pole', 'zero', 'damp', 'evalfr', 'freqresp', 'dcgain']
+__all__ = ['issiso', 'timebase', 'common_timebase', 'timebaseEqual', 
+           'isdtime', 'isctime', 'pole', 'zero', 'damp', 'evalfr', 
+           'freqresp', 'dcgain']
 
 class LTI:
     """LTI is a parent class to linear time-invariant (LTI) system objects.
@@ -157,9 +159,59 @@ def timebase(sys, strict=True):
 
     return sys.dt
 
+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 
+        or StateSpace system)
+
+    Returns
+    -------
+    dt: number
+        The common timebase of dt1 and dt2, as specified in 
+        :ref:`conventions-ref`. 
+        
+    Raises
+    ------
+    ValueError
+        when no compatible time base can be found
+    """
+    # explanation: 
+    # if either dt is None, they are compatible with anything
+    # 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 
+    if hasattr(dt1, 'dt'):
+        dt1 = dt1.dt
+    if hasattr(dt2, 'dt'):
+        dt2 = dt2.dt
+
+    if dt1 is None: 
+        return dt2
+    elif dt2 is None: 
+        return dt1
+    elif dt1 is True: 
+        if dt2 > 0:
+            return dt2
+        else: 
+            raise ValueError("Systems have incompatible timebases")
+    elif dt2 is True: 
+        if dt1 > 0: 
+            return dt1
+        else: 
+            raise ValueError("Systems have incompatible timebases")
+    elif np.isclose(dt1, dt2):
+        return dt1
+    else: 
+        raise ValueError("Systems have incompatible timebases")
+
 # Check to see if two timebases are equal
 def timebaseEqual(sys1, sys2):
-    """Check to see if two systems have the same timebase
+    """
+    Check to see if two systems have the same timebase
 
     timebaseEqual(sys1, sys2)
 
@@ -168,7 +220,10 @@ def timebaseEqual(sys1, sys2):
     discrete or continuous timebase systems.  If two systems have a discrete
     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", 
+         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
@@ -178,27 +233,6 @@ def timebaseEqual(sys1, sys2):
     else:
         return sys1.dt == sys2.dt
 
-# Find a common timebase between two or more systems
-def _find_timebase(sys1, *sysn):
-    """Find the common timebase between systems, otherwise return False"""
-
-    # Create a list of systems to check
-    syslist = [sys1]
-    syslist.append(*sysn)
-
-    # Look for a common timebase
-    dt = None
-
-    for sys in syslist:
-        # Make sure time bases are consistent
-        if (dt is None and sys.dt is not None) or \
-           (dt is True and isdiscrete(sys)):
-            # Timebase was not specified; set to match this system
-            dt = sys.dt
-        elif dt != sys.dt:
-            return False
-    return dt
-
 
 # Check to see if a system is a discrete time system
 def isdtime(sys, strict=False):
diff --git a/control/statesp.py b/control/statesp.py
index 522d187a9..b71248af7 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -61,7 +61,7 @@
 import scipy as sp
 from scipy.signal import lti, cont2discrete
 from warnings import warn
-from .lti import LTI, timebase, timebaseEqual, isdtime
+from .lti import LTI, common_timebase, isdtime
 from . import config
 from copy import deepcopy
 
@@ -71,7 +71,6 @@
 # Define module default parameter values
 _statesp_defaults = {
     'statesp.use_numpy_matrix': True,
-    'statesp.default_dt': None,
     'statesp.remove_useless_states': True, 
     }
 
@@ -143,22 +142,29 @@ class StateSpace(LTI):
     `numpy.ndarray` objects.  The :func:`~control.use_numpy_matrix` function
     can be used to set the storage type.
 
-    Discrete-time state space system are implemented by using the 'dt'
-    instance variable and setting it to the sampling period.  If 'dt' is not
-    None, then it must match whenever two state space systems are combined.
-    Setting dt = 0 specifies a continuous system, while leaving dt = None
-    means the system timebase is not specified.  If 'dt' is set to True, the
-    system will be treated as a discrete time system with unspecified sampling
-    time. The default value of 'dt' is None and can be changed by changing the 
-    value of ``control.config.defaults['statesp.default_dt']``.
-
+    A discrete time system is created by specifying a nonzero 'timebase', dt 
+    when the system is constructed:
+
+    * dt = 0: continuous time system (default)
+    * dt > 0: discrete time system with sampling period 'dt'
+    * dt = True: discrete time with unspecified sampling period
+    * dt = None: no timebase specified 
+
+    Systems must have compatible timebases in order to be combined. 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. Similarly, 
+    a system with timebase `None` can be combined with a system having any 
+    timebase; the result will have the timebase of the latter system. 
+    The default value of dt can be changed by changing the value of 
+    ``control.config.defaults['control.default_dt']``.
     """
 
     # Allow ndarray * StateSpace to give StateSpace._rmul_() priority
     __array_priority__ = 11     # override ndarray and matrix types
 
 
-    def __init__(self, *args, **kw):
+    def __init__(self, *args, **kwargs):
         """
         StateSpace(A, B, C, D[, dt])
 
@@ -171,13 +177,13 @@ def __init__(self, *args, **kw):
         call StateSpace(sys), where sys is a StateSpace object.
 
         """
+        # first get A, B, C, D matrices
         if len(args) == 4:
             # The user provided A, B, C, and D matrices.
             (A, B, C, D) = args
-            dt = config.defaults['statesp.default_dt']
         elif len(args) == 5:
             # Discrete time system
-            (A, B, C, D, dt) = args
+            (A, B, C, D, _) = args
         elif len(args) == 1:
             # Use the copy constructor.
             if not isinstance(args[0], StateSpace):
@@ -187,15 +193,11 @@ def __init__(self, *args, **kw):
             B = args[0].B
             C = args[0].C
             D = args[0].D
-            try:
-                dt = args[0].dt
-            except NameError:
-                dt = config.defaults['statesp.default_dt']
         else:
-            raise ValueError("Needs 1 or 4 arguments; received %i." % len(args))
+            raise ValueError("Expected 1, 4, or 5 arguments; received %i." % len(args))
 
         # Process keyword arguments
-        remove_useless = kw.get('remove_useless', config.defaults['statesp.remove_useless_states'])
+        remove_useless = kwargs.get('remove_useless', config.defaults['statesp.remove_useless_states'])
 
         # Convert all matrices to standard form
         A = _ssmatrix(A)
@@ -207,12 +209,33 @@ def __init__(self, *args, **kw):
         D = _ssmatrix(D)
 
         # TODO: use super here?
-        LTI.__init__(self, inputs=D.shape[1], outputs=D.shape[0], dt=dt)
+        LTI.__init__(self, inputs=D.shape[1], outputs=D.shape[0])
         self.A = A
         self.B = B
         self.C = C
         self.D = D
 
+        # now set dt
+        if len(args) == 4: 
+            if 'dt' in kwargs:
+                dt = kwargs['dt']
+            elif self.is_static_gain(): 
+                dt = None
+            else:
+                dt = config.defaults['control.default_dt']
+        elif len(args) == 5:
+            dt = args[4]
+            if 'dt' in kwargs:
+                warn('received multiple dt arguments, using positional arg dt=%s'%dt)
+        elif len(args) == 1:
+            try:
+                dt = args[0].dt
+            except NameError:
+                if self.is_static_gain(): 
+                    dt = None
+                else:
+                    dt = config.defaults['control.default_dt']
+        self.dt = dt
         self.states = A.shape[1]
 
         if 0 == self.states:
@@ -316,14 +339,7 @@ def __add__(self, other):
                     (self.outputs != other.outputs)):
                 raise ValueError("Systems have different shapes.")
 
-            # Figure out the sampling time to use
-            if self.dt is None and other.dt is not None:
-                dt = other.dt       # use dt from second argument
-            elif (other.dt is None and self.dt is not None) or \
-                    (timebaseEqual(self, other)):
-                dt = self.dt        # use dt from first argument
-            else:
-                raise ValueError("Systems have different sampling times")
+            dt = common_timebase(self.dt, other.dt)
 
             # Concatenate the various arrays
             A = concatenate((
@@ -372,16 +388,8 @@ def __mul__(self, other):
             # Check to make sure the dimensions are OK
             if self.inputs != other.outputs:
                 raise ValueError("C = A * B: A has %i column(s) (input(s)), \
-but B has %i row(s)\n(output(s))." % (self.inputs, other.outputs))
-
-            # Figure out the sampling time to use
-            if (self.dt == None and other.dt != None):
-                dt = other.dt       # use dt from second argument
-            elif (other.dt == None and self.dt != None) or \
-                    (timebaseEqual(self, other)):
-                dt = self.dt        # use dt from first argument
-            else:
-                raise ValueError("Systems have different sampling times")
+                    but B has %i row(s)\n(output(s))." % (self.inputs, other.outputs))
+            dt = common_timebase(self.dt, other.dt)
 
             # Concatenate the various arrays
             A = concatenate(
@@ -453,9 +461,8 @@ 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:
+            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
@@ -512,9 +519,8 @@ def freqresp(self, omega):
         # 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:
+            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
@@ -617,14 +623,7 @@ def feedback(self, other=1, sign=-1):
         if (self.inputs != other.outputs) or (self.outputs != other.inputs):
                 raise ValueError("State space systems don't have compatible inputs/outputs for "
                                  "feedback.")
-
-        # Figure out the sampling time to use
-        if self.dt is None and other.dt is not None:
-            dt = other.dt       # use dt from second argument
-        elif other.dt is None and self.dt is not None or timebaseEqual(self, other):
-            dt = self.dt        # use dt from first argument
-        else:
-            raise ValueError("Systems have different sampling times")
+        dt = common_timebase(self.dt, other.dt)
 
         A1 = self.A
         B1 = self.B
@@ -694,14 +693,7 @@ def lft(self, other, nu=-1, ny=-1):
         # dimension check
         # TODO
 
-        # Figure out the sampling time to use
-        if (self.dt == None and other.dt != None):
-            dt = other.dt       # use dt from second argument
-        elif (other.dt == None and self.dt != None) or \
-                timebaseEqual(self, other):
-            dt = self.dt        # use dt from first argument
-        else:
-            raise ValueError("Systems have different time bases")
+        dt = common_timebase(self.dt, other.dt)
 
         # submatrices
         A = self.A
@@ -815,8 +807,7 @@ def append(self, other):
         if not isinstance(other, StateSpace):
             other = _convertToStateSpace(other)
 
-        if self.dt != other.dt:
-            raise ValueError("Systems must have the same time step")
+        self.dt = common_timebase(self.dt, other.dt)
 
         n = self.states + other.states
         m = self.inputs + other.inputs
@@ -929,6 +920,12 @@ def dcgain(self):
             # eigenvalue at DC
             gain = np.tile(np.nan, (self.outputs, self.inputs))
         return np.squeeze(gain)
+    
+    def is_static_gain(self):
+        """True if and only if the system has no dynamics, that is, 
+        if A and B are zero. """
+        return not np.any(self.A) and not np.any(self.B)
+
 
 
 # TODO: add discrete time check
@@ -1246,8 +1243,7 @@ def _mimo2simo(sys, input, warn_conversion=False):
 
     return sys
 
-
-def ss(*args):
+def ss(*args, **kwargs):
     """ss(A, B, C, D[, dt])
 
     Create a state space system.
@@ -1292,8 +1288,7 @@ def ss(*args):
         Output matrix
     D: array_like or string
         Feed forward matrix
-    dt: If present, specifies the sampling period and a discrete time
-        system is created
+    dt: If present, specifies the timebase of the system
 
     Returns
     -------
@@ -1324,7 +1319,7 @@ def ss(*args):
     """
 
     if len(args) == 4 or len(args) == 5:
-        return StateSpace(*args)
+        return StateSpace(*args, **kwargs)
     elif len(args) == 1:
         from .xferfcn import TransferFunction
         sys = args[0]
@@ -1336,7 +1331,7 @@ def ss(*args):
             raise TypeError("ss(sys): sys must be a StateSpace or \
 TransferFunction object.  It is %s." % type(sys))
     else:
-        raise ValueError("Needs 1 or 4 arguments; received %i." % len(args))
+        raise ValueError("Needs 1, 4, or 5 arguments; received %i." % len(args))
 
 
 def tf2ss(*args):
diff --git a/control/tests/config_test.py b/control/tests/config_test.py
index 2fdae22e4..749fcc51b 100644
--- a/control/tests/config_test.py
+++ b/control/tests/config_test.py
@@ -218,16 +218,28 @@ def test_legacy_defaults(self):
         assert(isinstance(ct.ss(0,0,0,1).D, np.ndarray))
     
     def test_change_default_dt(self):
-        ct.set_defaults('statesp', default_dt=0)
-        self.assertEqual(ct.ss(0,0,0,1).dt, 0)
-        ct.set_defaults('statesp', default_dt=None)
-        self.assertEqual(ct.ss(0,0,0,1).dt, None)
-        ct.set_defaults('xferfcn', default_dt=0)
-        self.assertEqual(ct.tf(1, 1).dt, 0)
-        ct.set_defaults('xferfcn', default_dt=None)
-        self.assertEqual(ct.tf(1, 1).dt, None)
+        # test that system with dynamics uses correct default dt
+        ct.set_defaults('control', default_dt=0)
+        self.assertEqual(ct.tf(1, [1,1]).dt, 0)
+        self.assertEqual(ct.ss(1,0,0,1).dt, 0)
+        self.assertEqual(ct.iosys.NonlinearIOSystem(
+            lambda t, x, u: u*x*x, 
+            lambda t, x, u: x, inputs=1, outputs=1).dt, 0)
+        ct.set_defaults('control', default_dt=None)
+        self.assertEqual(ct.tf(1, [1,1]).dt, None)
+        self.assertEqual(ct.ss(1,0,0,1).dt, None)
+        self.assertEqual(ct.iosys.NonlinearIOSystem(
+            lambda t, x, u: u*x*x, 
+            lambda t, x, u: x, inputs=1, outputs=1).dt, None)
         
+        # test that static gain systems always have dt=None
+        ct.set_defaults('control', default_dt=0)
+        self.assertEqual(ct.tf(1, 1).dt, None)
+        self.assertEqual(ct.ss(0,0,0,1).dt, None)        
+        # TODO: add in test for static gain iosys
 
+        ct.reset_defaults()
+        
     def tearDown(self):
         # Get rid of any figures that we created
         plt.close('all')
diff --git a/control/tests/discrete_test.py b/control/tests/discrete_test.py
index 6598e3a81..bd0050fbc 100644
--- a/control/tests/discrete_test.py
+++ b/control/tests/discrete_test.py
@@ -7,7 +7,7 @@
 import numpy as np
 from control import StateSpace, TransferFunction, feedback, step_response, \
     isdtime, timebase, isctime, sample_system, bode, impulse_response, \
-    timebaseEqual, forced_response
+    common_timebase, forced_response
 from control import matlab
 
 class TestDiscrete(unittest.TestCase):
@@ -16,9 +16,10 @@ class TestDiscrete(unittest.TestCase):
     def setUp(self):
         """Set up a SISO and MIMO system to test operations on."""
 
-        # Single input, single output continuous and discrete time systems
+        # Single input, single output systems with unspecified and 
+        # discrete timebases
         sys = matlab.rss(3, 1, 1)
-        self.siso_ss1 = StateSpace(sys.A, sys.B, sys.C, sys.D)
+        self.siso_ss1 = StateSpace(sys.A, sys.B, sys.C, sys.D, None)
         self.siso_ss1c = StateSpace(sys.A, sys.B, sys.C, sys.D, 0.0)
         self.siso_ss1d = StateSpace(sys.A, sys.B, sys.C, sys.D, 0.1)
         self.siso_ss2d = StateSpace(sys.A, sys.B, sys.C, sys.D, 0.2)
@@ -29,7 +30,7 @@ def setUp(self):
         B = [[1., 4.], [-3., -3.], [-2., 1.]]
         C = [[4., 2., -3.], [1., 4., 3.]]
         D = [[-2., 4.], [0., 1.]]
-        self.mimo_ss1 = StateSpace(A, B, C, D)
+        self.mimo_ss1 = StateSpace(A, B, C, D, None)
         self.mimo_ss1c = StateSpace(A, B, C, D, 0)
 
         # Two input, two output discrete time system
@@ -39,19 +40,24 @@ def setUp(self):
         self.mimo_ss2d = StateSpace(A, B, C, D, 0.2)
 
         # Single input, single output continuus and discrete transfer function
-        self.siso_tf1 = TransferFunction([1, 1], [1, 2, 1])
+        self.siso_tf1 = TransferFunction([1, 1], [1, 2, 1],  None)
         self.siso_tf1c = TransferFunction([1, 1], [1, 2, 1], 0)
         self.siso_tf1d = TransferFunction([1, 1], [1, 2, 1], 0.1)
         self.siso_tf2d = TransferFunction([1, 1], [1, 2, 1], 0.2)
         self.siso_tf3d = TransferFunction([1, 1], [1, 2, 1], True)
 
-    def testTimebaseEqual(self):
-        self.assertEqual(timebaseEqual(self.siso_ss1, self.siso_tf1), True)
-        self.assertEqual(timebaseEqual(self.siso_ss1, self.siso_ss1c), True)
-        self.assertEqual(timebaseEqual(self.siso_ss1, self.siso_ss1d), True)
-        self.assertEqual(timebaseEqual(self.siso_ss1d, self.siso_ss1c), False)
-        self.assertEqual(timebaseEqual(self.siso_ss1d, self.siso_ss2d), False)
-        self.assertEqual(timebaseEqual(self.siso_ss1d, self.siso_ss3d), False)
+    def testCompatibleTimebases(self):
+        common_timebase(self.siso_ss1.dt, self.siso_tf1.dt)  # raises error 
+        common_timebase(self.siso_ss1.dt, self.siso_ss1c.dt) # if incompat
+        common_timebase(self.siso_ss1d.dt, self.siso_ss1.dt)
+        common_timebase(self.siso_ss1.dt, self.siso_ss1d.dt)
+        common_timebase(self.siso_ss1.dt, self.siso_ss1d.dt)
+        common_timebase(self.siso_ss1d.dt, self.siso_ss3d.dt)
+        common_timebase(self.siso_ss3d.dt, self.siso_ss1d.dt)
+        self.assertRaises(ValueError, common_timebase, 
+            self.siso_ss1d.dt, self.siso_ss1c.dt) # cont + discrete
+        self.assertRaises(ValueError, common_timebase, 
+            self.siso_ss1d.dt, self.siso_ss2d.dt) # incompatible discrete
 
     def testSystemInitialization(self):
         # Check to make sure systems are discrete time with proper variables
@@ -69,12 +75,24 @@ def testSystemInitialization(self):
         self.assertEqual(self.siso_tf2d.dt, 0.2)
         self.assertEqual(self.siso_tf3d.dt, True)
 
+        # keyword argument check
+        # dynamic systems
+        self.assertEqual(TransferFunction(1, [1, 1], dt=0.1).dt, 0.1)
+        self.assertEqual(TransferFunction(1, [1, 1], 0.1).dt, 0.1)
+        self.assertEqual(StateSpace(1,1,1,1, dt=0.1).dt, 0.1)
+        self.assertEqual(StateSpace(1,1,1,1, 0.1).dt, 0.1)
+        # static gain system, dt argument should still override default dt
+        self.assertEqual(TransferFunction(1, [1,], dt=0.1).dt, 0.1)
+        self.assertEqual(TransferFunction(1, [1,], 0.1).dt, 0.1)
+        self.assertEqual(StateSpace(0,0,1,1, dt=0.1).dt, 0.1)
+        self.assertEqual(StateSpace(0,0,1,1, 0.1).dt, 0.1)
+
     def testCopyConstructor(self):
         for sys in (self.siso_ss1, self.siso_ss1c, self.siso_ss1d):
-            newsys = StateSpace(sys);
+            newsys = StateSpace(sys)
             self.assertEqual(sys.dt, newsys.dt)
         for sys in (self.siso_tf1, self.siso_tf1c, self.siso_tf1d):
-            newsys = TransferFunction(sys);
+            newsys = TransferFunction(sys)
             self.assertEqual(sys.dt, newsys.dt)
 
     def test_timebase(self):
@@ -97,16 +115,17 @@ def test_timebase(self):
 
     def test_timebase_conversions(self):
         '''Check to make sure timebases transfer properly'''
-        tf1 = TransferFunction([1,1],[1,2,3])       # unspecified
+        tf1 = TransferFunction([1,1],[1,2,3], None) # unspecified
         tf2 = TransferFunction([1,1],[1,2,3], 0)    # cont time
         tf3 = TransferFunction([1,1],[1,2,3], True) # dtime, unspec 
-        tf4 = TransferFunction([1,1],[1,2,3], 1)    # dtime, dt=1
+        tf4 = TransferFunction([1,1],[1,2,3], 0.1)    # dtime, dt=1
 
         # Make sure unspecified timebase is converted correctly
         self.assertEqual(timebase(tf1*tf1), timebase(tf1))
         self.assertEqual(timebase(tf1*tf2), timebase(tf2))
         self.assertEqual(timebase(tf1*tf3), timebase(tf3))
         self.assertEqual(timebase(tf1*tf4), timebase(tf4))
+        self.assertEqual(timebase(tf3*tf4), timebase(tf4))
         self.assertEqual(timebase(tf2*tf1), timebase(tf2))
         self.assertEqual(timebase(tf3*tf1), timebase(tf3))
         self.assertEqual(timebase(tf4*tf1), timebase(tf4))
@@ -123,41 +142,23 @@ def test_timebase_conversions(self):
         self.assertEqual(timebase(tf3*tf3), timebase(tf3))
         self.assertEqual(timebase(tf3*tf4), timebase(tf4))
         self.assertEqual(timebase(tf3+tf3), timebase(tf3))
-        self.assertEqual(timebase(tf3+tf3), timebase(tf4))
+        self.assertEqual(timebase(tf3+tf4), timebase(tf4))
         self.assertEqual(timebase(feedback(tf3, tf3)), timebase(tf3))
         self.assertEqual(timebase(feedback(tf3, tf4)), timebase(tf4))
 
         # Make sure all other combinations are errors
-        try:
-            tf2*tf3             # Error; incompatible timebases
-            raise ValueError("incompatible operation allowed")
-        except ValueError:
-            pass
-        try:
-            tf2*tf4             # Error; incompatible timebases
-            raise ValueError("incompatible operation allowed")
-        except ValueError:
-            pass
-        try:
-            tf2+tf3             # Error; incompatible timebases
-            raise ValueError("incompatible operation allowed")
-        except ValueError:
-            pass
-        try:
-            tf2+tf4             # Error; incompatible timebases
-            raise ValueError("incompatible operation allowed")
-        except ValueError:
-            pass
-        try:
-            feedback(tf2, tf3)  # Error; incompatible timebases
-            raise ValueError("incompatible operation allowed")
-        except ValueError:
-            pass
-        try:
-            feedback(tf2, tf4)   # Error; incompatible timebases
-            raise ValueError("incompatible operation allowed")
-        except ValueError:
-            pass
+        with self.assertRaises(ValueError): tf2*tf3
+        with self.assertRaises(ValueError): tf3*tf2
+        with self.assertRaises(ValueError): tf2*tf4
+        with self.assertRaises(ValueError): tf4*tf2
+        with self.assertRaises(ValueError): tf2+tf3
+        with self.assertRaises(ValueError): tf3+tf2
+        with self.assertRaises(ValueError): tf2+tf4
+        with self.assertRaises(ValueError): tf4+tf2
+        with self.assertRaises(ValueError): feedback(tf2, tf3)
+        with self.assertRaises(ValueError): feedback(tf3, tf2)
+        with self.assertRaises(ValueError): feedback(tf2, tf4)
+        with self.assertRaises(ValueError): feedback(tf4, tf2)        
         
     def testisdtime(self):
         # Constant
@@ -214,13 +215,12 @@ def testAddition(self):
         sys = self.siso_ss1c + self.siso_ss1c
         sys = self.siso_ss1d + self.siso_ss1d
         sys = self.siso_ss3d + self.siso_ss3d
+        sys = self.siso_ss1d + self.siso_ss3d
         self.assertRaises(ValueError, StateSpace.__add__, self.mimo_ss1c,
                           self.mimo_ss1d)
         self.assertRaises(ValueError, StateSpace.__add__, self.mimo_ss1d,
                           self.mimo_ss2d)
-        self.assertRaises(ValueError, StateSpace.__add__, self.siso_ss1d,
-                          self.siso_ss3d)
-
+        
         # Transfer function addition
         sys = self.siso_tf1 + self.siso_tf1d
         sys = self.siso_tf1 + self.siso_tf1c
@@ -229,12 +229,11 @@ def testAddition(self):
         sys = self.siso_tf1c + self.siso_tf1c
         sys = self.siso_tf1d + self.siso_tf1d
         sys = self.siso_tf2d + self.siso_tf2d
+        sys = self.siso_tf1d + self.siso_tf3d
         self.assertRaises(ValueError, TransferFunction.__add__, self.siso_tf1c,
                           self.siso_tf1d)
         self.assertRaises(ValueError, TransferFunction.__add__, self.siso_tf1d,
                           self.siso_tf2d)
-        self.assertRaises(ValueError, TransferFunction.__add__, self.siso_tf1d,
-                          self.siso_tf3d)
 
         # State space + transfer function
         sys = self.siso_ss1c + self.siso_tf1c
@@ -252,12 +251,11 @@ def testMultiplication(self):
         sys = self.siso_ss1d * self.siso_ss1
         sys = self.siso_ss1c * self.siso_ss1c
         sys = self.siso_ss1d * self.siso_ss1d
+        sys = self.siso_ss1d * self.siso_ss3d
         self.assertRaises(ValueError, StateSpace.__mul__, self.mimo_ss1c,
                           self.mimo_ss1d)
         self.assertRaises(ValueError, StateSpace.__mul__, self.mimo_ss1d,
                           self.mimo_ss2d)
-        self.assertRaises(ValueError, StateSpace.__mul__, self.siso_ss1d,
-                          self.siso_ss3d)
 
         # Transfer function addition
         sys = self.siso_tf1 * self.siso_tf1d
@@ -266,12 +264,11 @@ def testMultiplication(self):
         sys = self.siso_tf1d * self.siso_tf1
         sys = self.siso_tf1c * self.siso_tf1c
         sys = self.siso_tf1d * self.siso_tf1d
+        sys = self.siso_tf1d * self.siso_tf3d
         self.assertRaises(ValueError, TransferFunction.__mul__, self.siso_tf1c,
                           self.siso_tf1d)
         self.assertRaises(ValueError, TransferFunction.__mul__, self.siso_tf1d,
                           self.siso_tf2d)
-        self.assertRaises(ValueError, TransferFunction.__mul__, self.siso_tf1d,
-                          self.siso_tf3d)
 
         # State space * transfer function
         sys = self.siso_ss1c * self.siso_tf1c
@@ -290,10 +287,10 @@ def testFeedback(self):
         sys = feedback(self.siso_ss1d, self.siso_ss1)
         sys = feedback(self.siso_ss1c, self.siso_ss1c)
         sys = feedback(self.siso_ss1d, self.siso_ss1d)
+        sys = feedback(self.siso_ss1d, self.siso_ss3d)
         self.assertRaises(ValueError, feedback, self.mimo_ss1c, self.mimo_ss1d)
         self.assertRaises(ValueError, feedback, self.mimo_ss1d, self.mimo_ss2d)
-        self.assertRaises(ValueError, feedback, self.siso_ss1d, self.siso_ss3d)
-
+        
         # Transfer function addition
         sys = feedback(self.siso_tf1, self.siso_tf1d)
         sys = feedback(self.siso_tf1, self.siso_tf1c)
@@ -301,10 +298,10 @@ def testFeedback(self):
         sys = feedback(self.siso_tf1d, self.siso_tf1)
         sys = feedback(self.siso_tf1c, self.siso_tf1c)
         sys = feedback(self.siso_tf1d, self.siso_tf1d)
+        sys = feedback(self.siso_tf1d, self.siso_tf3d)
         self.assertRaises(ValueError, feedback, self.siso_tf1c, self.siso_tf1d)
         self.assertRaises(ValueError, feedback, self.siso_tf1d, self.siso_tf2d)
-        self.assertRaises(ValueError, feedback, self.siso_tf1d, self.siso_tf3d)
-
+        
         # State space, transfer function
         sys = feedback(self.siso_ss1c, self.siso_tf1c)
         sys = feedback(self.siso_tf1c, self.siso_ss1c)
@@ -340,11 +337,28 @@ def test_sample_system(self):
             sysd = sample_system(sysc, 1, method="matched")
             self.assertEqual(sysd.dt, 1)
 
+        # bilinear approximation with prewarping test
+        wwarp = 50
+        Ts = 0.025
+        # test state space version
+        plant = self.siso_ss1c
+        plant_d_warped = plant.sample(Ts, 'bilinear', prewarp_frequency=wwarp)
+        np.testing.assert_array_almost_equal(
+            plant.evalfr(wwarp), 
+            plant_d_warped.evalfr(wwarp))
+        # test transfer function version
+        plant = self.siso_tf1c
+        plant_d_warped = plant.sample(Ts, 'bilinear', prewarp_frequency=wwarp)
+        np.testing.assert_array_almost_equal(
+            plant.evalfr(wwarp), 
+            plant_d_warped.evalfr(wwarp))
+
         # Check errors
         self.assertRaises(ValueError, sample_system, self.siso_ss1d, 1)
         self.assertRaises(ValueError, sample_system, self.siso_tf1d, 1)
         self.assertRaises(ValueError, sample_system, self.siso_ss1, 1, 'unknown')
 
+
     def test_sample_ss(self):
         # double integrators, two different ways
         sys1 = StateSpace([[0.,1.],[0.,0.]], [[0.],[1.]], [[1.,0.]], 0.)
diff --git a/control/tests/matlab_test.py b/control/tests/matlab_test.py
index f8b481248..e12167b70 100644
--- a/control/tests/matlab_test.py
+++ b/control/tests/matlab_test.py
@@ -613,20 +613,20 @@ def testCombi01(self):
         # margin command should remove the solution for w = nearly zero
 
         # Example is a concocted two-body satellite with flexible link
-        Jb = 400;
-        Jp = 1000;
-        k = 10;
-        b = 5;
+        Jb = 400
+        Jp = 1000
+        k = 10
+        b = 5
 
         # can now define an "s" variable, to make TF's
-        s = tf([1, 0], [1]);
-        hb1 = 1/(Jb*s);
-        hb2 = 1/s;
-        hp1 = 1/(Jp*s);
-        hp2 = 1/s;
+        s = tf([1, 0], [1])
+        hb1 = 1/(Jb*s)
+        hb2 = 1/s
+        hp1 = 1/(Jp*s)
+        hp2 = 1/s
 
         # convert to ss and append
-        sat0 = append(ss(hb1), ss(hb2), k, b, ss(hp1), ss(hp2));
+        sat0 = append(ss(hb1), ss(hb2), k, b, ss(hp1), ss(hp2))
 
         # connection of the elements with connect call
         Q = [[1, -3, -4],  # link moment (spring, damper), feedback to body
@@ -635,9 +635,9 @@ def testCombi01(self):
              [4,  1, -5],  # damper input
              [5,  3,  4],  # link moment, acting on payload
              [6,  5,  0]]
-        inputs = [1];
-        outputs = [1, 2, 5, 6];
-        sat1 = connect(sat0, Q, inputs, outputs);
+        inputs = [1]
+        outputs = [1, 2, 5, 6]
+        sat1 = connect(sat0, Q, inputs, outputs)
 
         # matched notch filter
         wno = 0.19
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
index 96404d79f..d7a3d948f 100644
--- a/control/tests/statesp_test.py
+++ b/control/tests/statesp_test.py
@@ -245,6 +245,26 @@ def test_freq_resp(self):
         np.testing.assert_almost_equal(phase, true_phase)
         np.testing.assert_equal(omega, true_omega)
 
+    def test_is_static_gain(self):
+        A0 = np.zeros((2,2))
+        A1 = A0.copy()
+        A1[0,1] = 1.1
+        B0 = np.zeros((2,1))
+        B1 = B0.copy()
+        B1[0,0] = 1.3
+        C0 = A0
+        C1 = np.eye(2)
+        D0 = 0
+        D1 = np.ones((2,1))
+        self.assertTrue(StateSpace(A0, B0, C1, D1).is_static_gain()) # True
+        # fix this once remove_useless_states is false by default
+        #print(StateSpace(A1, B0, C1, D1).is_static_gain()) # should be False when remove_useless is false
+        self.assertFalse(StateSpace(A0, B1, C1, D1).is_static_gain()) # False
+        self.assertFalse(StateSpace(A1, B1, C1, D1).is_static_gain()) # False
+        self.assertTrue(StateSpace(A0, B0, C0, D0).is_static_gain()) # True
+        self.assertTrue(StateSpace(A0, B0, C0, D1).is_static_gain()) # True
+        self.assertTrue(StateSpace(A0, B0, C1, D0).is_static_gain()) # True
+        
     @unittest.skipIf(not slycot_check(), "slycot not installed")
     def test_minreal(self):
         """Test a minreal model reduction."""
@@ -653,25 +673,5 @@ def test_copy_constructor(self):
         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..57e6205ef 100644
--- a/control/tests/xferfcn_test.py
+++ b/control/tests/xferfcn_test.py
@@ -349,7 +349,7 @@ def test_divide_siso(self):
 
     def test_div(self):
         # Make sure that sampling times work correctly
-        sys1 = TransferFunction([1., 3., 5], [1., 6., 2., -1])
+        sys1 = TransferFunction([1., 3., 5], [1., 6., 2., -1], None)
         sys2 = TransferFunction([[[-1., 3.]]], [[[1., 0., -1.]]], True)
         sys3 = sys1 / sys2
         self.assertEqual(sys3.dt, True)
@@ -409,6 +409,26 @@ def test_evalfr_siso(self):
             sys._evalfr(32.),
             np.array([[0.00281959302585077 - 0.030628473607392j]]))
 
+    def test_is_static_gain(self):
+        numstatic = 1.1
+        denstatic = 1.2
+        numdynamic = [1,1]
+        dendynamic = [2,1]
+        numstaticmimo = [[[1.1,], [1.2,]], [[1.2,], [0.8,]]]        
+        denstaticmimo = [[[1.9,], [1.2,]], [[1.2,], [0.8,]]]
+        numdynamicmimo = [[[1.1, 0.9], [1.2]], [[1.2], [0.8]]]        
+        dendynamicmimo = [[[1.1, 0.7], [0.2]], [[1.2], [0.8]]]
+        self.assertTrue(TransferFunction(numstatic,denstatic).is_static_gain()) # True
+        self.assertTrue(TransferFunction(numstaticmimo,denstaticmimo).is_static_gain()) # True
+        
+        self.assertFalse(TransferFunction(numstatic,dendynamic).is_static_gain()) # False
+        self.assertFalse(TransferFunction(numdynamic,dendynamic).is_static_gain()) # False
+        self.assertFalse(TransferFunction(numdynamic,denstatic).is_static_gain()) # False
+        self.assertFalse(TransferFunction(numstatic,dendynamic).is_static_gain()) # False
+        
+        self.assertFalse(TransferFunction(numstaticmimo,dendynamicmimo).is_static_gain()) # False
+        self.assertFalse(TransferFunction(numdynamicmimo,denstaticmimo).is_static_gain()) # False
+        
     # 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
@@ -659,7 +679,6 @@ def test_minreal_3(self):
         g = TransferFunction([1,1],[1,1]).minreal()
         np.testing.assert_array_almost_equal(1.0, g.num[0][0])
         np.testing.assert_array_almost_equal(1.0, g.den[0][0])
-        np.testing.assert_equal(None, g.dt)
 
     def test_minreal_4(self):
         """Check minreal on discrete TFs."""
@@ -913,26 +932,5 @@ def test_repr(self):
                         H.den[p][m], H2.den[p][m])
             self.assertEqual(H.dt, H2.dt)
     
-    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 = ss2tf(plant)
-        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/xferfcn.py b/control/xferfcn.py
index f50d5141d..2d3078fdc 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -62,18 +62,16 @@
 from warnings import warn
 from itertools import chain
 from re import sub
-from .lti import LTI, timebaseEqual, timebase, isdtime
+from .lti import LTI, common_timebase, isdtime
 from . import config
 
 __all__ = ['TransferFunction', 'tf', 'ss2tf', 'tfdata']
 
 
 # Define module default parameter values
-_xferfcn_defaults = {
-    'xferfcn.default_dt': None}
+_xferfcn_defaults = {}
 
 class TransferFunction(LTI):
-
     """TransferFunction(num, den[, dt])
 
     A class for representing transfer functions
@@ -89,13 +87,22 @@ class TransferFunction(LTI):
     means that the numerator of the transfer function from the 6th input to the
     3rd output is set to s^2 + 4s + 8.
 
-    Discrete-time transfer functions are implemented by using the 'dt'
-    instance variable and setting it to something other than 'None'.  If 'dt'
-    has a non-zero value, then it must match whenever two transfer functions
-    are combined.  If 'dt' is set to True, the system will be treated as a
-    discrete time system with unspecified sampling time. The default value of 
-    'dt' is None and can be changed by changing the value of 
-    ``control.config.defaults['xferfcn.default_dt']``.
+    A discrete time transfer function is created by specifying a nonzero 
+    'timebase' dt when the system is constructed:
+
+    * dt = 0: continuous time system (default)
+    * dt > 0: discrete time system with sampling period 'dt'
+    * dt = True: discrete time with unspecified sampling period
+    * dt = None: no timebase specified 
+
+    Systems must have compatible timebases in order to be combined. 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. Similarly, 
+    a system with timebase `None` can be combined with a system having any 
+    timebase; the result will have the timebase of the latter system. 
+    The default value of dt can be changed by changing the value of 
+    ``control.config.defaults['control.default_dt']``.
 
     The TransferFunction class defines two constants ``s`` and ``z`` that
     represent the differentiation and delay operators in continuous and
@@ -104,9 +111,9 @@ class TransferFunction(LTI):
 
     >>> s = TransferFunction.s
     >>> G  = (s + 1)/(s**2 + 2*s + 1)
-
     """
-    def __init__(self, *args):
+    
+    def __init__(self, *args, **kwargs):
         """TransferFunction(num, den[, dt])
 
         Construct a transfer function.
@@ -124,7 +131,6 @@ def __init__(self, *args):
         if len(args) == 2:
             # The user provided a numerator and a denominator.
             (num, den) = args
-            dt = config.defaults['xferfcn.default_dt']
         elif len(args) == 3:
             # Discrete time transfer function
             (num, den, dt) = args
@@ -136,11 +142,6 @@ def __init__(self, *args):
                                 % type(args[0]))
             num = args[0].num
             den = args[0].den
-            # TODO: not sure this can ever happen since dt is always present
-            try:
-                dt = args[0].dt
-            except NameError:   # pragma: no coverage
-                dt = config.defaults['xferfcn.default_dt']
         else:
             raise ValueError("Needs 1, 2 or 3 arguments; received %i."
                              % len(args))
@@ -198,12 +199,36 @@ def __init__(self, *args):
                 if zeronum:
                     den[i][j] = ones(1)
 
-        LTI.__init__(self, inputs, outputs, dt)
+        LTI.__init__(self, inputs, outputs)
         self.num = num
         self.den = den
 
         self._truncatecoeff()
 
+        # get dt
+        if len(args) == 2:
+            # no dt given in positional arguments
+            if 'dt' in kwargs:
+                dt = kwargs['dt']
+            elif self.is_static_gain():
+                dt = None
+            else: 
+                dt = config.defaults['control.default_dt']
+        elif len(args) == 3:
+            # Discrete time transfer function
+            if 'dt' in kwargs:
+                warn('received multiple dt arguments, using positional arg dt=%s'%dt)
+        elif len(args) == 1:
+            # TODO: not sure this can ever happen since dt is always present
+            try:
+                dt = args[0].dt
+            except NameError:   # pragma: no coverage
+                if self.is_static_gain():
+                    dt = None
+                else: 
+                    dt = config.defaults['control.default_dt']
+        self.dt = dt
+
     def __call__(self, s):
         """Evaluate the system's transfer function for a complex variable
 
@@ -369,14 +394,7 @@ def __add__(self, other):
                 "The first summand has %i output(s), but the second has %i."
                 % (self.outputs, other.outputs))
 
-        # Figure out the sampling time to use
-        if self.dt is None and other.dt is not None:
-            dt = other.dt       # use dt from second argument
-        elif (other.dt is None and self.dt is not None) or \
-             (timebaseEqual(self, other)):
-            dt = self.dt        # use dt from first argument
-        else:
-            raise ValueError("Systems have different sampling times")
+        dt = common_timebase(self.dt, other.dt)
 
         # Preallocate the numerator and denominator of the sum.
         num = [[[] for j in range(self.inputs)] for i in range(self.outputs)]
@@ -420,15 +438,8 @@ def __mul__(self, other):
         inputs = other.inputs
         outputs = self.outputs
 
-        # Figure out the sampling time to use
-        if self.dt is None and other.dt is not None:
-            dt = other.dt       # use dt from second argument
-        elif (other.dt is None and self.dt is not None) or \
-             (self.dt == other.dt):
-            dt = self.dt        # use dt from first argument
-        else:
-            raise ValueError("Systems have different sampling times")
-
+        dt = common_timebase(self.dt, other.dt)
+        
         # Preallocate the numerator and denominator of the sum.
         num = [[[0] for j in range(inputs)] for i in range(outputs)]
         den = [[[1] for j in range(inputs)] for i in range(outputs)]
@@ -471,14 +482,7 @@ def __rmul__(self, other):
         inputs = self.inputs
         outputs = other.outputs
 
-        # Figure out the sampling time to use
-        if self.dt is None and other.dt is not None:
-            dt = other.dt       # use dt from second argument
-        elif (other.dt is None and self.dt is not None) \
-                or (self.dt == other.dt):
-            dt = self.dt        # use dt from first argument
-        else:
-            raise ValueError("Systems have different sampling times")
+        dt = common_timebase(self.dt, other.dt)
 
         # Preallocate the numerator and denominator of the sum.
         num = [[[0] for j in range(inputs)] for i in range(outputs)]
@@ -518,14 +522,7 @@ def __truediv__(self, other):
                 "TransferFunction.__truediv__ is currently \
                 implemented only for SISO systems.")
 
-        # Figure out the sampling time to use
-        if self.dt is None and other.dt is not None:
-            dt = other.dt       # use dt from second argument
-        elif (other.dt is None and self.dt is not None) or \
-             (self.dt == other.dt):
-            dt = self.dt        # use dt from first argument
-        else:
-            raise ValueError("Systems have different sampling times")
+        dt = common_timebase(self.dt, other.dt)
 
         num = polymul(self.num[0][0], other.den[0][0])
         den = polymul(self.den[0][0], other.num[0][0])
@@ -625,9 +622,8 @@ def _evalfr(self, omega):
         # 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):
+            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
@@ -691,9 +687,8 @@ def freqresp(self, omega):
         # 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:
+            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])
@@ -736,15 +731,7 @@ def feedback(self, other=1, sign=-1):
             raise NotImplementedError(
                 "TransferFunction.feedback is currently only implemented "
                 "for SISO functions.")
-
-        # Figure out the sampling time to use
-        if self.dt is None and other.dt is not None:
-            dt = other.dt       # use dt from second argument
-        elif (other.dt is None and self.dt is not None) or \
-             (self.dt == other.dt):
-            dt = self.dt        # use dt from first argument
-        else:
-            raise ValueError("Systems have different sampling times")
+        dt = common_timebase(self.dt, other.dt)
 
         num1 = self.num[0][0]
         den1 = self.den[0][0]
@@ -1025,7 +1012,7 @@ def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None):
 
         Returns
         -------
-        sysd : StateSpace system
+        sysd : TransferFunction system
             Discrete time system, with sampling rate Ts
 
         Notes
@@ -1091,6 +1078,18 @@ def _dcgain_cont(self):
                         gain[i][j] = np.nan
         return np.squeeze(gain)
 
+    def is_static_gain(self):
+        """returns True if and only if all of the numerator and denominator 
+        polynomials of the (possibly MIMO) transfer funnction are zeroth order, 
+        that is, if the system has no dynamics. """
+        for list_of_polys in self.num, self.den: 
+            for row in list_of_polys:
+                for poly in row:
+                    if len(poly) > 1: 
+                        return False
+        return True
+
+
 
 # c2d function contributed by Benjamin White, Oct 2012
 def _c2d_matched(sysC, Ts):
@@ -1302,7 +1301,7 @@ def _convert_to_transfer_function(sys, **kw):
     raise TypeError("Can't convert given type to TransferFunction system.")
 
 
-def tf(*args):
+def tf(*args, **kwargs):
     """tf(num, den[, dt])
 
     Create a transfer function system. Can create MIMO systems.
@@ -1392,7 +1391,7 @@ def tf(*args):
     """
 
     if len(args) == 2 or len(args) == 3:
-        return TransferFunction(*args)
+        return TransferFunction(*args, **kwargs)
     elif len(args) == 1:
         # Look for special cases defining differential/delay operator
         if args[0] == 's':
@@ -1413,7 +1412,7 @@ def tf(*args):
         raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
 
 
-def ss2tf(*args):
+def ss2tf(*args, **kwargs):
     """ss2tf(sys)
 
     Transform a state space system to a transfer function.
@@ -1478,7 +1477,7 @@ def ss2tf(*args):
     from .statesp import StateSpace
     if len(args) == 4 or len(args) == 5:
         # Assume we were given the A, B, C, D matrix and (optional) dt
-        return _convert_to_transfer_function(StateSpace(*args))
+        return _convert_to_transfer_function(StateSpace(*args, **kwargs))
 
     elif len(args) == 1:
         sys = args[0]
@@ -1564,7 +1563,6 @@ def _clean_part(data):
 
     return data
 
-
 # Define constants to represent differentiation, unit delay
 TransferFunction.s = TransferFunction([1, 0], [1], 0)
 TransferFunction.z = TransferFunction([1, 0], [1], True)
diff --git a/doc/conventions.rst b/doc/conventions.rst
index f07b51238..4cccf351e 100644
--- a/doc/conventions.rst
+++ b/doc/conventions.rst
@@ -80,27 +80,24 @@ Discrete time systems
 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 (default)
-* dt = 0: continuous time system
+* dt = 0: continuous time system (default)
 * dt > 0: discrete time system with sampling period 'dt'
 * dt = True: discrete time with unspecified sampling period
+* dt = None: no timebase specified 
 
 Only the :class:`StateSpace`, :class:`TransferFunction`, and
 :class:`InputOutputSystem` classes allow explicit representation of
 discrete time systems.
 
-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
+Systems must have compatible timebases in order to be combined. 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.  Similarly, 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. 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`. The default value of 'dt' can be changed by
-changing the values of ``control.config.defaults['statesp.default_dt']`` and 
-``control.config.defaults['xferfcn.default_dt']``.
+changing the value of ``control.config.defaults['control.default_dt']``.
 
 Conversion between representations
 ----------------------------------