function D1 = mtimes(D1,D2,ScalarFlags)
% Multiplies two LFT models: D = D1 * D2

%   Copyright 1986-2010 The MathWorks, Inc.
if isfinite(D1) && isfinite(D2)
   if nargin==3
      if ScalarFlags(1)
         [nz2,~] = iosize(D2);
         D1 = diagrep(D1,nz2);
      elseif ScalarFlags(2)
         [~,nw1] = iosize(D1);
         D2 = diagrep(D2,nw1);
      end
   end
   [nw1,~] = iosize(D1.Blocks);
   [nw2,nz2] = iosize(D2.Blocks);
   [ny2,nu2] = iosize(D2);
   nyz2 = ny2+nz2;  nuw2 = nu2+nw2;
   yperm = [1:ny2 nyz2+1:nyz2+nw1 ny2+1:nyz2];
   uperm = [1:nu2 nuw2+1:nuw2+nw1 nu2+1:nuw2];
   
   % Build IC
   if isnumeric(D1.IC)
      IC1 = blkdiag(D1.IC,eye(nz2));
      IC2 = blkdiag(D2.IC,eye(nw1));
      IC = IC1 * IC2(yperm,uperm);
   else
      IC1 = appendGain(D1.IC,eye(nz2));
      IC2 = appendGain(D2.IC,eye(nw1));
      IC = mtimes(IC1,ioperm(IC2,yperm,uperm),false(1,2));
   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
   % One model is NaN
   if nargin<3 || ~any(ScalarFlags)
      [rs,~] = iosize(D1);
      [~,cs] = iosize(D2);
      D1 = createGain(D1,NaN(rs,cs));
   elseif ScalarFlags(1)
      D1 = createGain(D1,NaN(iosize(D2)));
   else
      D1 = createGain(D1,NaN(iosize(D1)));
   end
end