function [D1,SingularFlag] = lft(D1,D2,indr1,inds1,indr2,inds2)
% LFT interconnection of LFT models.
%
%                        +-------+
%            u1 -------->|       |-------> y1
%                        |   M   |
%            r1    +---->|       |-----+   s1
%                  |     +-------+     |
%                  |                   |
%                  |     +-------+     |
%            s2    +-----|       |<----+   r2
%                        |   N   |
%           y2 <---------|       |-------- u2
%                        +-------+

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

% Channel sizes
ns1 = length(inds1);
nr1 = length(indr1);
[ny1,nu1,nz1,nw1] = getSignalSize(D1);
[ny2,nu2,nz2,nw2] = getSignalSize(D2);
ny1 = ny1-ns1;  nu1 = nu1-nr1;
ny2 = ny2-nr1;  nu2 = nu2-ns1;
nuw1 = nu1+nw1;  nyz1 = ny1+nz1;

if isfinite(D1) && isfinite(D2)
   % Center blocks to prevent artificial singularities
   if nw1>0 && nz1>0
      D1 = centerBlocks(D1);
   end
   if nw2>0 && nz2>0
      D2 = centerBlocks(D2);
   end
   
   % Build IC
   uperm = [1:nu1 nuw1+1:nuw1+nu2 nu1+1:nuw1 nuw1+nu2+1:nuw1+nu2+nw2];
   yperm = [1:ny1 nyz1+1:nyz1+ny2 ny1+1:nyz1 nyz1+ny2+1:nyz1+ny2+nz2];
   if isnumeric(D1.IC)
      IC = lft(D1.IC,D2.IC,indr1,inds1,indr2,inds2);
      SingularFlag = ~allfinite(IC);
      IC = IC(yperm,uperm);
   else
      [IC,SingularFlag] = lft(D1.IC,D2.IC,indr1,inds1,indr2,inds2);
      IC = ioperm(IC,yperm,uperm);
   end
   D1.IC = IC;
   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(ny1+ny2,nu1+nu2));
   SingularFlag = false;
end