function [D,xkeep] = getsubsys(D,rowIndex,colIndex,varargin)
% Extracts subsystem.
%
%    Dsub = getsubsys(D,rowIndex,colIndex) extracts the subsystem 
%    sys(rowIndex,colIndex). Dsub and D have the same number of 
%    states and delays.
%
%    [Dsub,xkeep] = getsubsys(D,rowIndex,colIndex,'smin') further
%    eliminates states and delays that are structurally decoupled 
%    from the subsystem. The boolean vector XKEEP indicates which 
%    states were kept and which states were discarded.

%   Author(s): P. Gahinet
%   Copyright 1986-2020 The MathWorks, Inc.

% NOTE:
% 1) Don't systematically eliminate non-sminimal states. This messes
%    up applications like LQG design where gain and estimators need to
%    be designed with same number of states, see MILLDEMO.
% 2) Deleting inputs is interpreted as setting their value to zero.

% Sizes
[nyz,nuw] = size(D.d);
nfd = numel(D.Delay.Internal);
ny = nyz-nfd;
nu = nuw-nfd;

% Delays
D.Delay.Input = D.Delay.Input(colIndex,:);  % colon needed if scalar and colIndex=[1 1]
D.Delay.Output = D.Delay.Output(rowIndex,:);

% Absorb input offsets for discarded inputs (consistent with setting these
% input signals to zero)
jDel = 1:nu; jDel(:,colIndex) = [];
u0Del = D.u0(jDel,:);
if any(u0Del)
   D.dx0 = D.dx0-D.b(:,jDel)*u0Del;
   D.y0 = D.y0-D.d(:,jDel)*u0Del;
end

% Account for internal delays
if nfd>0
   % Map ROWINDEX and COLINDEX to integer indices into D.d
   if ischar(rowIndex)
      rowIndex = 1:ny;
   elseif isa(rowIndex,'logical')
      rowIndex = find(rowIndex);
   end
   if ischar(colIndex)
      colIndex = 1:nu;
   elseif isa(colIndex,'logical')
      colIndex = find(colIndex);
   end
   rowIndex = [reshape(rowIndex,[1 numel(rowIndex)]) ny+1:nyz];
   colIndex = [reshape(colIndex,[1 numel(colIndex)]) nu+1:nuw];
end

% Select I/O channels
D.d = D.d(rowIndex,colIndex);
D.b = D.b(:,colIndex);
D.u0 = D.u0(colIndex,:);
D.c = D.c(rowIndex,:);
D.y0 = D.y0(rowIndex,:);

% Optionally discard disconnected states and delays
if nargin>3
   [D,xkeep] = sminreal(D);
else
   xkeep = true(size(D.x0));
end
