function [D,T] = compreal(D,type)
% Computes companion realization.

%   Copyright 2023 The MathWorks, Inc.
if hasInternalDelay(D)
   throw(ltipack.utNoDelaySupport('compreal',D.Ts,'internal'))
end
EXPLICIT = isempty(D.e);

% Check properness and eliminate E
[isProper,D] = isproper(D,true);
if ~isProper
   error(message('Control:general:NotSupportedImproperSys','compreal'))
end

if type=='c'
   a = D.a;  b = D.b;  c = D.c;
else
   a = D.a';  b = D.c';  c = D.b';
end

[nx,nu] = size(b);
if nu==0 || nx==0
   T = eye(nx);
else
   T = ctrb(a,b(:,1));
   % Note: For some reason, computing the last column of A as T\(A^n^b) 
   % gives a better match than using poly(A)
   aux = matlab.internal.math.nowarn.mldivide(T,[a*T(:,nx) b(:,2:nu)]);
   if ~allfinite(aux)
      error(message('Control:transformation:compreal1'))
   end
   a = [[zeros(1,nx-1);eye(nx-1)] , aux(:,1)];
   b = [[1;zeros(nx-1,1)] , aux(:,2:nu)];
   c = c*T;
end

if type=='c'
   aux = T\[D.dx0 D.x0];
   D.dx0 = aux(:,1);  D.x0 = aux(:,2);
else
   a = a';  tmp = b;  b = c';  c = tmp';  T = T';
   D.dx0 = T*D.dx0;  D.x0 = T*D.x0;  
end

D.a = a; D.b = b; D.c = c; D.e = [];
D.StateName = strings(0,1);
D.StatePath = strings(0,1);
D.StateUnit = strings(0,1);
D.Scaled = false;
if ~EXPLICIT
   T = [];
end
