function D = setsubsys(D,rowIndex,colIndex,rhs)
% Modifies subsystem via SYS(i,j) = RHS.
% D and RHS are @mechdata objects with the same sampling time

%   Copyright 2020 The MathWorks, Inc.

%RE: Assumes SYS(I,J) is a subsystem of SYS (no size growing)
[ny,nu] = iosize(D);

if isequal(rhs,[])
   % SYS(i,:) = [] or SYS(:,j) = []
   iscolon = strcmp({rowIndex,colIndex},':');
   rowSel = 1:ny;
   colSel = 1:nu;
   if all(iscolon)
      % sys(:,:) = [] produces a 0-by-nu system as for matrices
      rowSel = [];
   elseif iscolon(1)
      colSel(colIndex) = [];
   else
      rowSel(rowIndex) = [];
   end
   % Use GETSUBSYS to get structurally minimal number of internal delays
   % RE: Do not eliminate any state to be consistent with SUBSREF
   D = getsubsys(D,rowSel,colSel);

elseif strcmp(rowIndex,':') && strcmp(colIndex,':')
   % Optimized code for sys(:,:,ind) = rhs
   if all(iosize(rhs)==1) && (nu~=1 || ny~=1)
      % Scalar expansion
      rhs = iorep(rhs,iosize(D));
   end
   D = rhs;
   
elseif ~(isempty(rowIndex) || isempty(colIndex))
   % SYS(ROWINDEX,COLINDEX) = RHS
   % Use GETSUBSYS and IOCAT to correctly handle delays
   [rowIndex,colIndex,rperm,cperm,frows,fcols,ScalarFlag,iRepeat] = ...
      ltipack.checkRowColIndices(rowIndex,colIndex,iosize(D),iosize(rhs));
   if ScalarFlag
      % Scalar expansion
      rhs = iorep(rhs,[numel(rowIndex),numel(colIndex)]);
   elseif ~isempty(iRepeat)
      % Deal with repeated indices
      rhs = LocalGetSubSys(rhs,iRepeat{:});
   end
   % Reassign (2,2) block of SYS(RPERM,CPERM):
   %    [sys11  sys12 ;         [sys11  sys12;
   %     sys21  sys22 ]   -->    sys21   rhs ]
   if isempty(frows) && isempty(fcols)
      % Full reassignment
      D = rhs;
   elseif isempty(frows)
      % 2-Block case (e.g. Case for ny = 0)
      % Assemble result as [sys11 , rhs]
      sys11 = getsubsys(D,rperm,fcols,'smin');
      D = iocat(2,sys11,rhs);
   elseif isempty(fcols)
      % 2-Block Case (e.g. Case for nu = 0)
      % Assemble result as [sys11 ; rhs]
      sys11 = getsubsys(D,frows,cperm,'smin');
      D = iocat(1,sys11,rhs);
   else
      % Compute structural order of [[sys11 sys12] ; [sys21 rhs]]
      sysfr = LocalGetSubSys(D,frows,cperm);
      sys21 = LocalGetSubSys(D,rowIndex,fcols);
      Naug1 = order(sysfr) + order(sys21);
      
      % Compute structural order of [[sys11 ; sys21] , [sys12 ; rhs]]
      sysfc = LocalGetSubSys(D,rperm,fcols);
      sys12 = LocalGetSubSys(D,frows,colIndex);
      Naug2 = order(sysfc) + order(sys12);
      
      % Construct minimum-order realization for modified lhs
      if Naug1<=Naug2
         % Assemble result as [[sys11 , sys12] ; [sys21 , rhs]]
         D = iocat(1,sysfr,iocat(2,sys21,rhs));
      else
         % Assemble result as [[sys11 ; sys21] , [sys12 ; rhs]]
         D = iocat(2,sysfc,iocat(1,sys12,rhs));
      end
   end
   
   % Undo I/O permutation
   D.b(:,cperm) = D.b(:,1:nu);
   D.d(:,cperm) = D.d(:,1:nu);
   D.Delay.Input(cperm,:) = D.Delay.Input;
   D.f(rperm,:) = D.f(1:ny,:);
   D.g(rperm,:) = D.g(1:ny,:);
   D.d(rperm,:) = D.d(1:ny,:);
   D.Delay.Output(rperm,:) = D.Delay.Output;
end

%-------- Local Functions ------------------------

function D = LocalGetSubSys(D,rowIndex,colIndex)
% Wrapper around GETSUBSYS to discard all states when rowIndex
% or colIndex is empty (see TSUBS:lvlTwo2)
nr = length(rowIndex);
nc = length(colIndex);
if nr==0 || nc==0
   % Blast all states and delays
   D = createGain(D,zeros(nr,nc));
else
   D = getsubsys(D,rowIndex,colIndex,'smin');
end
