def _compute_tfinal_and_dt(sys, is_step=True):
    '\n    Helper function to estimate a final time and a sampling period for\n    time domain simulations. It is essentially geared towards impulse response\n    but is also used for step responses.\n\n    For discrete-time models, obviously dt is inherent and only tfinal is\n    computed.\n\n    Parameters\n    ----------\n    sys : {State, Transfer}\n        The system to be investigated\n    is_step : bool\n        Scales the dc value by the magnitude of the nonzero mode since\n        integrating the impulse response gives ∫exp(-λt) = -exp(-λt)/λ.\n        Default is True.\n\n    Returns\n    -------\n    tfinal : float\n        The final time instance for which the simulation will be performed.\n    dt : float\n        The estimated sampling period for the simulation.\n\n    Notes\n    -----\n    Just by evaluating the fastest mode for dt and slowest for tfinal often\n    leads to unnecessary, bloated sampling (e.g., Transfer(1,[1,1001,1000]))\n    since dt will be very small and tfinal will be too large though the fast\n    mode hardly ever contributes. Similarly, change the numerator to [1, 2, 0]\n    and the simulation would be unnecessarily long and the plot is virtually\n    an L shape since the decay is so fast.\n\n    Instead, a modal decomposition in time domain hence a truncated ZIR and ZSR\n    can be used such that only the modes that have significant effect on the\n    time response are taken. But the sensitivity of the eigenvalues complicate\n    the matter since dλ = <w, dA*v> with <w,v> = 1. Hence we can only work\n    with simple poles with this formulation. See Golub, Van Loan Section 7.2.2\n    for simple eigenvalue sensitivity about the nonunity of <w,v>. The size of\n    the response is dependent on the size of the eigenshapes rather than the\n    eigenvalues themselves.\n\n    '
    sqrt_eps = np.sqrt(np.spacing(1.0))
    min_points = 100
    min_points_z = 20
    max_points = 10000
    max_points_z = 75000
    default_tfinal = 5
    total_cycles = 5
    pts_per_cycle = 25
    log_decay_percent = np.log(100)
    if sys._isgain:
        if sys._isdiscrete:
            return ((sys._dt * min_points_z), sys._dt)
        else:
            return (default_tfinal, (default_tfinal / min_points))
    if sys._isdiscrete:
        dt = sys._dt
        tfinal = default_tfinal
        p = eigvals(sys.a)
        m_u = (np.abs(p) >= (1 + sqrt_eps))
        (p_u, p) = (p[m_u], p[(~ m_u)])
        if (p_u.size > 0):
            m_u = ((p_u.real < 0) & (np.abs(p_u.imag) < sqrt_eps))
            t_emp = np.max((log_decay_percent / np.abs((np.log(p_u[(~ m_u)]) / dt))))
            tfinal = max(tfinal, t_emp)
        m_z = (np.abs(p) < sqrt_eps)
        p = p[(~ m_z)]
        m_nr = ((p.real < 0) & (np.abs(p.imag) < sqrt_eps))
        (p_nr, p) = (p[m_nr], p[(~ m_nr)])
        if (p_nr.size > 0):
            t_emp = np.max((log_decay_percent / np.abs((np.log(p_nr) / dt).real)))
            tfinal = max(tfinal, t_emp)
        m_int = (((p.real - 1) < sqrt_eps) & (np.abs(p.imag) < sqrt_eps))
        (p_int, p) = (p[m_int], p[(~ m_int)])
        m_w = (np.abs((np.abs(p) - 1)) < sqrt_eps)
        (p_w, p) = (p[m_w], p[(~ m_w)])
        if (p_w.size > 0):
            t_emp = (((total_cycles * 2) * np.pi) / np.abs((np.log(p_w) / dt)).min())
            tfinal = max(tfinal, t_emp)
        if (p.size > 0):
            t_emp = (log_decay_percent / np.abs((np.log(p) / dt).real).min())
            tfinal = max(tfinal, t_emp)
        if (p_int.size > 0):
            tfinal = (tfinal * 5)
        num_samples = (tfinal // dt)
        if (num_samples > max_points_z):
            tfinal = (dt * max_points_z)
        else:
            tfinal = (dt * num_samples)
        return (tfinal, dt)
    (b, (sca, perm)) = matrix_balance(sys.a, separate=True)
    (p, l, r) = eig(b, left=True, right=True)
    eig_sens = reciprocal(maximum(1e-12, einsum('ij,ij->j', l, r).real))
    eig_sens = minimum(1000000000000.0, eig_sens)
    p[(np.abs(p) < np.spacing((eig_sens * norm(b, 1))))] = 0.0
    l[perm, :] *= reciprocal(sca)[:, None]
    r[perm, :] *= sca[:, None]
    (w, v) = ((sys.c @ r), (l.T.conj() @ sys.b))
    origin = False
    wn = np.abs(p)
    if np.any((wn == 0.0)):
        origin = True
    dc = zeros_like(p, dtype=float)
    ok = (np.abs(eig_sens) <= (1 / sqrt_eps))
    dc[ok] = ((norm(v[ok, :], axis=1) * norm(w[:, ok], axis=0)) * eig_sens[ok])
    dc[(wn != 0.0)] /= (wn[(wn != 0)] if is_step else 1.0)
    dc[(wn == 0.0)] = 0.0
    dc[(p.imag != 0.0)] *= 2
    relevance = ((dc > (0.1 * dc.max())) | (~ ok))
    psub = p[relevance]
    wnsub = wn[relevance]
    (tfinal, dt) = ([], [])
    ints = (wnsub == 0.0)
    iw = ((psub.imag != 0.0) & (np.abs(psub.real) <= sqrt_eps))
    if np.any(iw):
        tfinal += (((total_cycles * 2) * np.pi) / wnsub[iw]).tolist()
        dt += (((2 * np.pi) / pts_per_cycle) / wnsub[iw]).tolist()
    texp_mode = (log_decay_percent / np.abs(psub[((~ iw) & (~ ints))].real))
    tfinal += texp_mode.tolist()
    dt += minimum((texp_mode / 50), (((2 * np.pi) / pts_per_cycle) / wnsub[((~ iw) & (~ ints))])).tolist()
    if (len(tfinal) == 0):
        return ((default_tfinal * 5), ((default_tfinal * 5) / min_points))
    tfinal = (np.max(tfinal) * (5 if origin else 1))
    dt = np.min(dt)
    dt = ((tfinal / max_points) if ((tfinal // dt) > max_points) else dt)
    tfinal = ((dt * min_points) if ((tfinal // dt) < min_points) else tfinal)
    return (tfinal, dt)
Helper function to estimate a final time and a sampling period for
time domain simulations. It is essentially geared towards impulse response
but is also used for step responses.

For discrete-time models, obviously dt is inherent and only tfinal is
computed.

Parameters
----------
sys : {State, Transfer}
    The system to be investigated
is_step : bool
    Scales the dc value by the magnitude of the nonzero mode since
    integrating the impulse response gives ∫exp(-λt) = -exp(-λt)/λ.
    Default is True.

Returns
-------
tfinal : float
    The final time instance for which the simulation will be performed.
dt : float
    The estimated sampling period for the simulation.

Notes
-----
Just by evaluating the fastest mode for dt and slowest for tfinal often
leads to unnecessary, bloated sampling (e.g., Transfer(1,[1,1001,1000]))
since dt will be very small and tfinal will be too large though the fast
mode hardly ever contributes. Similarly, change the numerator to [1, 2, 0]
and the simulation would be unnecessarily long and the plot is virtually
an L shape since the decay is so fast.

Instead, a modal decomposition in time domain hence a truncated ZIR and ZSR
can be used such that only the modes that have significant effect on the
time response are taken. But the sensitivity of the eigenvalues complicate
the matter since dλ = <w, dA*v> with <w,v> = 1. Hence we can only work
with simple poles with this formulation. See Golub, Van Loan Section 7.2.2
for simple eigenvalue sensitivity about the nonunity of <w,v>. The size of
the response is dependent on the size of the eigenshapes rather than the
eigenvalues themselves.
def _compute_tfinal_and_dt(sys, is_step=True):
    '\n    Helper function to estimate a final time and a sampling period for\n    time domain simulations. It is essentially geared towards impulse response\n    but is also used for step responses.\n\n    For discrete-time models, obviously dt is inherent and only tfinal is\n    computed.\n\n    Parameters\n    ----------\n    sys : {State, Transfer}\n        The system to be investigated\n    is_step : bool\n        Scales the dc value by the magnitude of the nonzero mode since\n        integrating the impulse response gives ∫exp(-λt) = -exp(-λt)/λ.\n        Default is True.\n\n    Returns\n    -------\n    tfinal : float\n        The final time instance for which the simulation will be performed.\n    dt : float\n        The estimated sampling period for the simulation.\n\n    Notes\n    -----\n    Just by evaluating the fastest mode for dt and slowest for tfinal often\n    leads to unnecessary, bloated sampling (e.g., Transfer(1,[1,1001,1000]))\n    since dt will be very small and tfinal will be too large though the fast\n    mode hardly ever contributes. Similarly, change the numerator to [1, 2, 0]\n    and the simulation would be unnecessarily long and the plot is virtually\n    an L shape since the decay is so fast.\n\n    Instead, a modal decomposition in time domain hence a truncated ZIR and ZSR\n    can be used such that only the modes that have significant effect on the\n    time response are taken. But the sensitivity of the eigenvalues complicate\n    the matter since dλ = <w, dA*v> with <w,v> = 1. Hence we can only work\n    with simple poles with this formulation. See Golub, Van Loan Section 7.2.2\n    for simple eigenvalue sensitivity about the nonunity of <w,v>. The size of\n    the response is dependent on the size of the eigenshapes rather than the\n    eigenvalues themselves.\n\n    '
    sqrt_eps = np.sqrt(np.spacing(1.0))
    min_points = 100
    min_points_z = 20
    max_points = 10000
    max_points_z = 75000
    default_tfinal = 5
    total_cycles = 5
    pts_per_cycle = 25
    log_decay_percent = np.log(100)
    if sys._isgain:
        if sys._isdiscrete:
            return ((sys._dt * min_points_z), sys._dt)
        else:
            return (default_tfinal, (default_tfinal / min_points))
    if sys._isdiscrete:
        dt = sys._dt
        tfinal = default_tfinal
        p = eigvals(sys.a)
        m_u = (np.abs(p) >= (1 + sqrt_eps))
        (p_u, p) = (p[m_u], p[(~ m_u)])
        if (p_u.size > 0):
            m_u = ((p_u.real < 0) & (np.abs(p_u.imag) < sqrt_eps))
            t_emp = np.max((log_decay_percent / np.abs((np.log(p_u[(~ m_u)]) / dt))))
            tfinal = max(tfinal, t_emp)
        m_z = (np.abs(p) < sqrt_eps)
        p = p[(~ m_z)]
        m_nr = ((p.real < 0) & (np.abs(p.imag) < sqrt_eps))
        (p_nr, p) = (p[m_nr], p[(~ m_nr)])
        if (p_nr.size > 0):
            t_emp = np.max((log_decay_percent / np.abs((np.log(p_nr) / dt).real)))
            tfinal = max(tfinal, t_emp)
        m_int = (((p.real - 1) < sqrt_eps) & (np.abs(p.imag) < sqrt_eps))
        (p_int, p) = (p[m_int], p[(~ m_int)])
        m_w = (np.abs((np.abs(p) - 1)) < sqrt_eps)
        (p_w, p) = (p[m_w], p[(~ m_w)])
        if (p_w.size > 0):
            t_emp = (((total_cycles * 2) * np.pi) / np.abs((np.log(p_w) / dt)).min())
            tfinal = max(tfinal, t_emp)
        if (p.size > 0):
            t_emp = (log_decay_percent / np.abs((np.log(p) / dt).real).min())
            tfinal = max(tfinal, t_emp)
        if (p_int.size > 0):
            tfinal = (tfinal * 5)
        num_samples = (tfinal // dt)
        if (num_samples > max_points_z):
            tfinal = (dt * max_points_z)
        else:
            tfinal = (dt * num_samples)
        return (tfinal, dt)
    (b, (sca, perm)) = matrix_balance(sys.a, separate=True)
    (p, l, r) = eig(b, left=True, right=True)
    eig_sens = reciprocal(maximum(1e-12, einsum('ij,ij->j', l, r).real))
    eig_sens = minimum(1000000000000.0, eig_sens)
    p[(np.abs(p) < np.spacing((eig_sens * norm(b, 1))))] = 0.0
    l[perm, :] *= reciprocal(sca)[:, None]
    r[perm, :] *= sca[:, None]
    (w, v) = ((sys.c @ r), (l.T.conj() @ sys.b))
    origin = False
    wn = np.abs(p)
    if np.any((wn == 0.0)):
        origin = True
    dc = zeros_like(p, dtype=float)
    ok = (np.abs(eig_sens) <= (1 / sqrt_eps))
    dc[ok] = ((norm(v[ok, :], axis=1) * norm(w[:, ok], axis=0)) * eig_sens[ok])
    dc[(wn != 0.0)] /= (wn[(wn != 0)] if is_step else 1.0)
    dc[(wn == 0.0)] = 0.0
    dc[(p.imag != 0.0)] *= 2
    relevance = ((dc > (0.1 * dc.max())) | (~ ok))
    psub = p[relevance]
    wnsub = wn[relevance]
    (tfinal, dt) = ([], [])
    ints = (wnsub == 0.0)
    iw = ((psub.imag != 0.0) & (np.abs(psub.real) <= sqrt_eps))
    if np.any(iw):
        tfinal += (((total_cycles * 2) * np.pi) / wnsub[iw]).tolist()
        dt += (((2 * np.pi) / pts_per_cycle) / wnsub[iw]).tolist()
    texp_mode = (log_decay_percent / np.abs(psub[((~ iw) & (~ ints))].real))
    tfinal += texp_mode.tolist()
    dt += minimum((texp_mode / 50), (((2 * np.pi) / pts_per_cycle) / wnsub[((~ iw) & (~ ints))])).tolist()
    if (len(tfinal) == 0):
        return ((default_tfinal * 5), ((default_tfinal * 5) / min_points))
    tfinal = (np.max(tfinal) * (5 if origin else 1))
    dt = np.min(dt)
    dt = ((tfinal / max_points) if ((tfinal // dt) > max_points) else dt)
    tfinal = ((dt * min_points) if ((tfinal // dt) < min_points) else tfinal)
    return (tfinal, dt)