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

%   Copyright 2020 The MathWorks, Inc.
if isfinite(D1) && isfinite(D2)
   if nargin>2
      if ScalarFlags(1)
         D1 = diagrep(D1,numel(D2.Delay.Output));
      elseif ScalarFlags(2)
         D2 = diagrep(D2,numel(D1.Delay.Input));
      end
   end

   % Fold in delays squeezed between D1 and D2
   if utIsIOScaling(D1)
      % D1=diag(exp(-tau*s))
      D1.Delay.Output = D1.Delay.Output + D1.Delay.Input + D2.Delay.Output;
   elseif utIsIOScaling(D2)
      % D2=diag(exp(-tau*s))
      D2.Delay.Input = D2.Delay.Input + D2.Delay.Output + D1.Delay.Input;
   else
      D2 = utFoldDelay(D2,[],D1.Delay.Input + D2.Delay.Output);
   end
   
   % Internal delays and states
   [n,nu2] = iosize(D2);
   nfd1 = length(D1.Delay.Internal);
   nfd2 = length(D2.Delay.Internal);
   nq1 = size(D1.k,1);
   nq2 = size(D2.k,1);
   nq = nq1+nq2;
   
   % Compute realization 
   %   (M s^2 + C s + K) X(s) = B U(s)
   %   Y(s) = (F + G s) X(s) + D U(s)
   % for product D1 * D2:
   %   [ m1 s^2 + c1 s + k1        0           -b1(:,1:n) | b1(:,n+1:e)    0    ]
   %   [         0         m2 s^2 + c2 s + k2      0      |     0         b2    ]
   %   [         0         (f2 + g2 s)(1:n,:)     -I      |     0    -d2(1:n,:) ]
   %   [------------------------------------------------------------------------]
   %   [     f1 + g1 s             0            d1(:,1:n) | d1(:,n+1:e)    0    ]
   %   [         0        (f2 + g2 s)(n+1:e,:)     0      |      0  d2(n+1:e,:) ]
   % The internal signals (outputs of D2 = inputs of D1) are added as 
   % extra states to preserve sparsity.
   m = ltipack.util.blkdiagE(D1.m,D2.m,nq1,nq2);
   c = ltipack.util.blkdiag(D1.c,D2.c);
   k = ltipack.util.blkdiag(D1.k,D2.k);
   bF = [D1.b(:,1:n);sparse(nq2,n)];
   fF = [sparse(n,nq1) D2.f(1:n,:)];
   gF = [sparse(n,nq1) D2.g(1:n,:)];
   b = ltipack.util.blkdiag(D1.b(:,n+1:n+nfd1),D2.b);
   f = ltipack.util.blkdiag(D1.f,D2.f(n+1:n+nfd2,:));
   g = ltipack.util.blkdiag(D1.g,D2.g(n+1:n+nfd2,:));
   xF = [D1.d(:,1:n) ; sparse(nfd2,n)];
   yF = [sparse(n,nfd1) , D2.d(1:n,:)];
   d = ltipack.util.blkdiag(D1.d(:,n+1:n+nfd1),D2.d(n+1:n+nfd2,:));
   nALG = n;
   
   % Use Schur complement to eliminate algebraic states that can be 
   % eliminated without altering M,C,K
   idel = find(~(any(bF,1) & (any(fF,2) | any(gF,2))'));
   if ~isempty(idel)
%       fprintf('mtimes: eliminating %d AV, dA = %.3g\n',...
%          numel(idel),nnz(bF(:,idel) * fF(idel,:))+nnz(bF(:,idel) * gF(idel,:)))
      ikeep = 1:n;  ikeep(:,idel) = [];
      xF2 = xF(:,idel);  yF2 = yF(idel,:);
      b = b + bF(:,idel) * yF2;
      f = f + xF2 * fF(idel,:);
      g = g + xF2 * gF(idel,:);
      d = d + xF2 * yF2;
      bF = bF(:,ikeep);  fF = fF(ikeep,:);  gF = gF(ikeep,:);
      xF = xF(:,ikeep);  yF = yF(ikeep,:);
      nALG = numel(ikeep);
   end
   
   % Finalize realization
   D1.c = [[c ; gF] , sparse(nq+nALG,nALG)];
   D1.k = [k -bF ; fF -speye(nALG)];
   D1.b = [b ; -yF];
   D1.f = [f , xF];
   D1.g = [g , sparse(size(g,1),nALG)];
   D1.d = d;
   if nALG>0
      if isequal(m,[])
         D1.m = sparse(1:nq,1:nq,1,nq+nALG,nq+nALG);
      else
         D1.m = ltipack.util.blkdiag(m,sparse(nALG,nALG));
      end
      D1.StateInfo = [D1.StateInfo;D2.StateInfo;struct('Type',3,'Name',"",'Size',nALG)];
   else
      D1.m = m;
      D1.StateInfo = [D1.StateInfo;D2.StateInfo];
   end
   
   % Delays
   D1.Delay.Input = D2.Delay.Input;
   D1.Delay.Internal = [D1.Delay.Internal ; D2.Delay.Internal];
   if nfd1>0
      % If there are internal delays, reorder u and y to put them last
      uperm = [nfd1+1:nfd1+nu2 1:nfd1 nfd1+nu2+1:nfd1+nfd2+nu2];
      D1.b = D1.b(:,uperm);
      D1.d = D1.d(:,uperm);
   end
   
   D1.SolverConfig = ltioptions.sparss.reconcile(D1.SolverConfig,D2.SolverConfig);
else
   % One model is NaN
   if nargin<3 || ~any(ScalarFlags)
      D1 = createGain(D1,NaN(size(D1.d,1),size(D2.d,2)));
   elseif ScalarFlags(1)
      D1 = createGain(D1,NaN(iosize(D2)));
   else
      D1 = createGain(D1,NaN(iosize(D1)));
   end
end
