function [sysT, Q] = output_state_transform(sys, output_name, opts)
    arguments
        sys ss;
        output_name (1,1) string;
        opts.Unitary (1,1) logical = false; % Ensures state is identical to output, setting to false may leave scaling in C, but may be numerically better behaved.
        opts.Orthonormal (1,1) logical = false;
        opts.PreconditionQ (1,1) logical = true;
    end

    yidx = startsWith(sys.OutputName, output_name); ny = sum(yidx);
    if ny == 0, warning("No output starts with '%s'", output_name); end

    if opts.PreconditionQ
        Cy = sys.C(yidx, :);
        T = diag(fillmissing(sum(abs(Cy)) ./ sum(Cy ~= 0), "constant", 1));
        sys = ss2ss(sys, T);
    end

    Cy = sys.C(yidx, :);
    Ct = Cy(1:rank(Cy), :);

    nullC = null(Cy);
    if opts.Orthonormal
        colnorm = zeros(1, size(nullC, 2));
        for i = 1:size(nullC, 2)
           colnorm(i) = norm(nullC(:,i));
        end
        nullC = nullC ./ colnorm;
    end

    P = Ct*Ct';
    if ~opts.Unitary
        P = P / norm(P);
    end
    Q = [nullC, Ct'/P]; % It is crucial that this matrix is well conditioned.

    % essentially do "sysT = ss2ss(sys, inv(Q));" but in numerically more stable:
    sysT = sys;
    sysT.A = double(Q \ sys.A * Q);
    sysT.B = double(Q \ sys.B);
    sysT.C = double(sys.C * Q);
end
