function D = xelim(D,elim,method)
% Eliminates states

%   Author(s): J.N. Little, P. Gahinet
%   Copyright 1986-2023 The MathWorks, Inc.
nx = size(D.a,1);
if nx==0
   return
elseif hasInternalDelay(D)
   throw(ltipack.utNoDelaySupport('xelim',D.Ts,'internal'))
end
keep = 1:nx;
keep(elim) = [];
nx1 = length(keep);
Ts = D.Ts;

switch method
   case 'MatchDC'
      % Matched DC gains: partition into x1, to be kept, and x2, to be eliminated:
      a = D.a;  b = D.b;  c = D.c;  d = D.d;  e = D.e;
      if ~D.Scaled
         % Scale data (required by ELIMAV and orthogonal transformation below)
         [a,b,c,e] = ltipack.xscale(a,b,c,d,e,Ts,'Warn',false);
      end
      % Put eliminated states last
      xperm = [keep(:);elim(:)];
      ie = nx1+1:nx;
      a = a(xperm,xperm);   
      b = b(xperm,:);
      c = c(:,xperm);
      if isempty(e)
         if Ts~=0
            % Discrete-time system: DC corresponds to z=1
            a(ie,ie) = a(ie,ie) - eye(nx-nx1);
         end
      else
         e = e(xperm,xperm);
         % Transform to E=[E11 E12;0 E22] so that zeroing out dx2 is equivalent
         % to zeroing out E12 and E22
         if any(e(ie,1:nx1),'all')
            [Q,~] = qr(e(:,1:nx1));
            e = Q'*e;  a = Q'*a;  b = Q'*b;
         end
         if Ts~=0
            % Discrete-time system: DC corresponds to setting x2[k+1]-x2[k]=0
            a(:,ie) = a(:,ie) - e(:,ie);
         end
         % Keep only E11 for elimAV
         e(ie,:) = 0;   e(:,ie) = 0;
      end
      % Eliminate algebraic variables
      [D.a,D.b,D.c,D.d,D.e,SingularFlag] = elimAV(a,b,c,d,e,nx1,eps^(1/3));
      if SingularFlag
         % Error when MatchDC leads to singular (A,E) pencil
         error(message('Control:transformation:xelim1'))
      end
      nxr = size(D.a,1);
      D.dx0 = zeros(nxr,1);
      D.x0 = zeros(nxr,1);
      
   case 'Truncate'
      % Simply delete specified states
      D.a = D.a(keep,keep);
      D.b = D.b(keep,:);
      D.c = D.c(:,keep);
      if ~isempty(D.e)
         D.e = D.e(keep,keep);
      end
      D.x0 = D.x0(keep,:);
      D.dx0 = D.dx0(keep,:);
      nxr = nx1;
end

% Build output
if ~isempty(D.StateName)
   D.StateName = [D.StateName(keep) ; strings(nxr-nx1,1)];
end
if ~isempty(D.StatePath)
   D.StatePath = [D.StatePath(keep) ; strings(nxr-nx1,1)];
end
if ~isempty(D.StateUnit)
   D.StateUnit = [D.StateUnit(keep) ; strings(nxr-nx1,1)];
end
D.Scaled = false;
