function [D1,SingularFlag] = feedback(D1,D2,indu2,indy2,sign)
% Feedback interconnection of LFT models.
%
%                      +--------+
%         u1 --------->|        |--------> y1
%                      |   D1   |
%         u2 --->O---->|        |----+---> y2
%                |     +--------+    |
%                |                   |
%                +-----[   D2   ]<---+

%   Author(s): P. Gahinet
%   Copyright 1986-2011 The MathWorks, Inc.
if nargin<3
   [nu2,ny2] = iosize(D2); 
   indy2 = 1:ny2;
   indu2 = 1:nu2;
   sign = -1;
else
   ny2 = length(indy2);
   nu2 = length(indu2);
end
[ny12,nu12,nz1,nw1] = getSignalSize(D1);
nyz1 = ny12+nz1;  nuw1 = nu12+nw1;  % sizes of D1.IC

if isfinite(D1) && isfinite(D2)
   % Center blocks to prevent artificial singularities
   if nw1>0 && nz1>0
      D1 = centerBlocks(D1);
   end
   if ~isempty(D2.Blocks)
      D2 = centerBlocks(D2);
   end
   
   % Build IC
   if isnumeric(D1.IC)
      iu = [1:nuw1 indu2]; % replicate u2
      iy = [1:nyz1 indy2]; % replicate y2
      IC2 = D2.IC;
      if sign<0
         % Note: Do not multiply block channels by SIGN
         IC2(1:nu2,:) = -IC2(1:nu2,:);
      end
      D1.IC = lft(D1.IC(iy,iu),IC2,nuw1+1:nuw1+nu2,nyz1+1:nyz1+ny2,1:ny2,1:nu2);
      SingularFlag = ~allfinite(D1.IC);
   else
      % Note: Do not replicate u2,y2 and use LFT (leads to duplicating the input
      % and output delays of D1 in the feedback path)
      [D,SingularFlag] = feedbackG(D1.IC,D2.IC,indu2,indy2,1:ny2,1:nu2,sign,+1);
      % Eliminate I/Os of D2 except block channels
      [rs2,cs2] = iosize(D2.IC);
      D1.IC = getsubsys(D,[1:nyz1 nyz1+nu2+1:nyz1+rs2],[1:nuw1 nuw1+ny2+1:nuw1+cs2]);
   end
   D1.Blocks = append(D1.Blocks,D2.Blocks);
   
   % Try simplifying uncertain blocks based on their AutoSimplify settings
   if any(isUncertain(D1.Blocks))
      D1 = simplify(D1,-1,Inf);
   end
else
   % NaN systems
   D1 = createGain(D1,NaN(nyz1,nuw1));
   SingularFlag = false;
end
