function D = setsubsys(D,rowIndex,colIndex,rhs)
% Modifies portion of model via M(i,j) = RHS.

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

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

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 eliminate blocks and delays that are disconnected from 
   % the model I/Os and do not contribute to the state matrix (full state
   % vector is preserved)
   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
      if NumFlag
         rhs = repmat(rhs,iosize(D));
      else
         rhs = iorep(rhs,iosize(D));
      end
   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
      if isnumeric(rhs)
         rhs = repmat(rhs,[numel(rowIndex),numel(colIndex)]);
      else
         rhs = iorep(rhs,[numel(rowIndex),numel(colIndex)]);
      end
   elseif ~isempty(iRepeat)
      % Deal with repeated indices
      rhs = getsubsys(rhs,iRepeat{:});
   end
   
   % Reassign (2,2) block of SYS(RPERM,CPERM):
   %    [sys11  sys12 ;         [sys11  sys12;
   %     sys21  sys22 ]   -->    sys21   rhs ]
   % Try row-wise and column-wise approaches and pick result with fewer blocks
   if isempty(frows) && isempty(fcols),
      % Full reassignment
      D = rhs;
   elseif isempty(frows)
      % Assemble result as [sys11 , rhs]
      sys11 = getsubsys(D,rperm,fcols,'smin');
      D = iocat(2,sys11,rhs);
   elseif isempty(fcols)
      % 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]]
      D1_ = getsubsys(D,frows,cperm,'smin');
      D21 = getsubsys(D,rowIndex,fcols,'smin');
      Nblk1 = numel(D1_.Blocks) + numel(D21.Blocks);
      
      % Compute structural order of [[sys11 ; sys21] , [sys12 ; rhs]]
      D_1 = getsubsys(D,rperm,fcols,'smin');
      D12 = getsubsys(D,frows,colIndex,'smin');
      Nblk2 = numel(D_1.Blocks) + numel(D12.Blocks);
      
      % Construct minimum-order realization for modified lhs
      if Nblk1<=Nblk2,
         % Assemble result as [[sys11 , sys12] ; [sys21 , rhs]]
         D = iocat(1,D1_,iocat(2,D21,rhs));
      else
         % Assemble result as [[sys11 ; sys21] , [sys12 ; rhs]]
         D = iocat(2,D_1,iocat(1,D12,rhs));
      end
   end

   % Undo I/O permutation
   [nu,ny] = iosize(D.Blocks);
   [nz,nw] = iosize(D);
   rp = 1:ny+nz;  rp(rperm) = 1:nz;
   cp = 1:nu+nw;  cp(cperm) = 1:nw;   
   if NumFlag
      D.IC = D.IC(rp,cp);
   else
      D.IC = ioperm(D.IC,rp,cp);
   end
   
   % Try simplifying uncertain blocks based on their AutoSimplify settings
   if any(isUncertain(D.Blocks))
      D = simplify(D,-1,Inf);
   end
end