def modal_form(xsys, pmax=None):
    """Convert a system into modal canonical form

    Parameters
    ----------
    xsys : StateSpace object
        System to be transformed, with state `x`
    pmax: float
        An upper bound on individual transformations.  Default eps ** -0.25.

    Returns
    -------
    zsys : StateSpace object
        System in modal canonical form, with state `z`
    T : matrix
        Coordinate transformation: z = T * x

    Notes
    -----
    pmax is passed to slycot.mb03rd.  The pmax value allows for a
    trade-off: larger values give smaller modal blocks, but
    worse-conditioned T.
    """

    if pmax is None:
        pmax = finfo(float64).eps ** -0.5

    aschur, tschur = schur(xsys.A)
    _, tmodal, _, _ = mb03rd(aschur.shape[0], aschur, tschur, pmax=pmax)
    return similarity_transform(xsys, tmodal, inverse=True), tmodal


def similarity_transform(xsys, T, timescale=1, inverse=False):
    """Perform a similarity transformation, with option time rescaling.

    Transform a linear state space system to a new state space representation
    z = T x, where T is an invertible matrix.

    Parameters
    ----------
    xsys : StateSpace object
           System to transform
    T : 2D invertible array
        The matrix `T` defines the new set of coordinates z = T x.
    timescale : float
        If present, also rescale the time unit to tau = timescale * t
    inverse: boolean
        If True (default), transform so z = T x.  If False, transform
        so x = T z.

    Returns
    -------
    zsys : StateSpace object
        System in transformed coordinates, with state 'z'

    """
    # Create a new system, starting with a copy of the old one
    zsys = StateSpace(xsys)

    # Define a function to compute the right inverse (solve x M = y)
    def rsolve(M, y):
        return transpose(solve(transpose(M), transpose(y)))

    # Update the system matrices
    if not inverse:
        zsys.A = rsolve(T, dot(T, zsys.A)) / timescale
        zsys.B = dot(T, zsys.B) / timescale
        zsys.C = rsolve(T, zsys.C)
    else:
        zsys.A = solve(T, zsys.A).dot(T) / timescale
        zsys.B = solve(T, zsys.B) / timescale
        zsys.C = zsys.C.dot(T)

    return zsys

from control import canonical_form, ss
import numpy as np

from scipy.linalg import schur
from slycot import mb03rd
from control import similarity_transform

import control as ct

s = ct.tf([1,0],[1])
g_tf = 1 / (s+1)**3 * 1 / (s**2 + s + 100)**2

g = ct.ss(g_tf)
gmod, tmod = ct.modal_form(g)

with np.printoptions(precision=3, suppress=True):
    print(gmod.A)

print(f'{np.linalg.cond(tmod) = }')