function D = ctranspose(D)
% Pertransposition of state-space models.
%
%   For a continuous-time model SYS with data (A,B,C,D), 
%   CTRANSPOSE produces the state-space model TSYS with
%   data (-A',-C',B',D').  If H(s) is the transfer function 
%   of SYS, then H(-s).' is the transfer function of TSYS.
%
%   For a discrete-time model SYS with data (A,B,C,D), TSYS
%   is the state-space model with data 
%       (AA, AA*C', -B'*AA, D'-B'*AA*C')  with AA=inv(A').
%   Equivalently, H(z^-1).' is the transfer function of TSYS
%   if H(z) is the transfer function of SYS.

%   Copyright 1986-2023 The MathWorks, Inc.
Ts = D.Ts;

% Compute pertranspose
% Note: This operation is applied to the incremental model with state
% x-x0, input u-u0, output y-y0
if Ts==0
   % Continuous-time case
   if any(D.dx0)
      error(message('Control:transformation:ctranspose2'))
   end
   D.a = -D.a';
   aux = D.b';
   D.b = -D.c';
   D.c = aux;
   D.d = D.d';
   D.e = D.e';
else
   % Discrete time
   if any(D.dx0-D.x0)
      error(message('Control:transformation:ctranspose2'))
   end
   [a,b,c,d,e] = getABCDE(D);
   nx = size(a,1);
   [ny,nu] = size(d);
   if nu<=ny
      % First realize H(1/z) as
      %   z [A B;0 0] [x;u] = [E 0;0 I] [x;u] + [0;-I] u
      %                   y = [C 0] [x;u] + D u
      % then transpose.
      aux = eye(nu,'like',a);
      D.e = ltipack.util.blkdiag([a';b'],zeros(0,nu));
      D.a = ltipack.util.blkdiag(e',aux);
      D.b = ltipack.util.blkdiag(c',zeros(nu,0));
      D.c = ltipack.util.blkdiag(zeros(0,nx),-aux);
      D.x0 = [D.x0;D.u0(1:nu,:)];
      D.dx0 = D.x0;
   else
      % First transpose H(z) to
      %   E'x = A'x + C'y,   u = B'x+D'y
      % then replace z->1/z to obtain
      %   z [A' C';0 0] [x;y] = [E' 0;0 I] [x;y] + [0;-I] y
      %                     u = [B' 0] [x;y] + D' y
      aux = eye(ny,'like',a);
      D.e = ltipack.util.blkdiag([a' c'],zeros(ny,0));
      D.a = ltipack.util.blkdiag(e',aux);
      D.b = ltipack.util.blkdiag(zeros(nx,0),-aux);
      D.c = ltipack.util.blkdiag(b',zeros(0,ny));
      D.x0 = [D.x0;D.y0(1:ny,:)];
      D.dx0 = D.x0;
   end
   D.d = D.d';
   D.Scaled = false;
   % Note: Do not reduce here. This requires scaling and users would have
   % no way to address scaling issues
end
id = D.Delay.Input;
D.Delay.Input = D.Delay.Output;
D.Delay.Output = id;
aux = D.u0;  D.u0 = D.y0;  D.y0 = aux;