function D1 = mtimes(D1,D2,ScalarFlags)
% Multiplies two state-space models D = D1 * D2

%   Copyright 1986-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
   % Note: 
   % (1) Do not consolidate delays at the inputs of D2 or the outputs of
   %     D1 when D2 or D1 is "scalar". This does not work with offsets and
   %     does not work with IC since it changes (delays) the state. For a 
   %     SISO model, the response with initial condition XINIT is not
   %     the same with InputDelay=2 and InputDelay=OutputDelay=1.
   % (2) An exception is made when the scalar model has no states (gain or
   %     pure delay).
   % (3) The InputDelay of D1 is not updated (not looked up again)
   if utIsIOScaling(D1)
      % D1=diag(exp(-tau*s)): Consolidate delays at outputs of D1
      D1.Delay.Output = D1.Delay.Output + D1.Delay.Input + D2.Delay.Output;
   elseif utIsIOScaling(D2)
      % D2=diag(exp(-tau*s)): Consolidate delays at inputs of D2
      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);
   nx1 = size(D1.a,1);
   nx2 = size(D2.a,1);
   
   % Compute realization for product D1 * D2
   D1.Delay.Input = D2.Delay.Input;
   D1.a = [D1.a , D1.b(:,1:n)*D2.c(1:n,:) ; zeros(nx2,nx1) , D2.a];
   D1.e = ltipack.util.blkdiagE(D1.e,D2.e,nx1,nx2);
   D1.x0 = [D1.x0 ; D2.x0];
   if nfd1==0 && nfd2==0
      % No internal delays
      %     [ a1  b1*c2 ]       [ b1*d2 ]
      % A = [  0    a2  ]   B = [   b2  ]
      %
      % C = [ c1  d1*c2 ]   D =  d1*d2
      aux = D2.y0-D1.u0;  % y02-u01
      D1.dx0 = [D1.dx0+D1.b*aux ; D2.dx0];
      D1.y0 = D1.y0+D1.d*aux;
      D1.u0 = D2.u0;
      D1.b = [D1.b*D2.d ; D2.b];
      D1.c = [D1.c , D1.d*D2.c];
      D1.d = D1.d * D2.d;
   else
      aux = D2.y0(1:n,:)-D1.u0(1:n,:);  % y02-u01
      D1.dx0 = [D1.dx0+D1.b(:,1:n)*aux ; D2.dx0];
      D1.y0 = [D1.y0+D1.d(:,1:n)*aux ; D2.y0(n+1:n+nfd2)];
      D1.u0 = [D2.u0(1:nu2,:) ; D1.u0(n+1:n+nfd1,:) ; D2.u0(nu2+1:nu2+nfd2)];
      b = [D1.b(:,1:n)*D2.d(1:n,:) ; D2.b];
      D1.b = [b(:,1:nu2) [D1.b(:,n+1:n+nfd1);zeros(nx2,nfd1)] b(:,nu2+1:nu2+nfd2)];
      D1.c = [D1.c , D1.d(:,1:n)*D2.c(1:n,:) ; ...
         zeros(nfd2,nx1) , D2.c(n+1:n+nfd2,:)];
      d = [D1.d(:,1:n)*D2.d(1:n,:) ; D2.d(n+1:n+nfd2,:)];
      D1.d = [d(:,1:nu2) [D1.d(:,n+1:n+nfd1);zeros(nfd2,nfd1)] d(:,nu2+1:nu2+nfd2)];
      D1.Delay.Internal = [D1.Delay.Internal ; D2.Delay.Internal];
   end
   D1.Scaled = false;
   if allfinite(D1.d)
      D1 = appendStateInfo(D1,D2,nx1,nx2);
   else
      % NaNs may arise in sys1/sys2
      D1 = createGain(D1,nan(iosize(D1)));
   end
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
