function D1 = append(D1,D2)
% Appends inputs and outputs of two models.

%   Copyright 2020 The MathWorks, Inc.
if isfinite(D1) && isfinite(D2)
   % Sizes
   nfd1 = length(D1.Delay.Internal);
   nfd2 = length(D2.Delay.Internal);
   nq1 = size(D1.k,1);
   nq2 = size(D2.k,1);
   
   % Compute realization for concatenation
   D1.k = ltipack.util.blkdiag(D1.k,D2.k);
   D1.c = ltipack.util.blkdiag(D1.c,D2.c);
   D1.m = ltipack.util.blkdiagE(D1.m,D2.m,nq1,nq2);
   if nfd1==0 && nfd2==0
      % No internal delays
      D1.b = ltipack.util.blkdiag(D1.b,D2.b);
      D1.f = ltipack.util.blkdiag(D1.f,D2.f);
      D1.g = ltipack.util.blkdiag(D1.g,D2.g);
      D1.d = ltipack.util.blkdiag(D1.d,D2.d);
   else
      % Account for internal delays
      [rs1,cs1] = size(D1.d);
      [rs2,cs2] = size(D2.d);
      ny1 = rs1-nfd1;  nu1 = cs1-nfd1;
      ny2 = rs2-nfd2;  nu2 = cs2-nfd2;
      D1.b = [ltipack.util.blkdiag(D1.b(:,1:nu1),D2.b(:,1:nu2)) ,...
         ltipack.util.blkdiag(D1.b(:,nu1+1:cs1),D2.b(:,nu2+1:cs2))];
      D1.f = [ltipack.util.blkdiag(D1.f(1:ny1,:),D2.f(1:ny2,:)) ;...
         ltipack.util.blkdiag(D1.f(ny1+1:rs1,:),D2.f(ny2+1:rs2,:))];
      D1.g = [ltipack.util.blkdiag(D1.g(1:ny1,:),D2.g(1:ny2,:)) ;...
         ltipack.util.blkdiag(D1.g(ny1+1:rs1,:),D2.g(ny2+1:rs2,:))];
      D1.d = [ltipack.util.blkdiag(D1.d(1:ny1,1:nu1),D2.d(1:ny2,1:nu2)) , ...
         ltipack.util.blkdiag(D1.d(1:ny1,nu1+1:cs1),D2.d(1:ny2,nu2+1:cs2)) ; ...
         ltipack.util.blkdiag(D1.d(ny1+1:rs1,1:nu1),D2.d(ny2+1:rs2,1:nu2)) , ...
         ltipack.util.blkdiag(D1.d(ny1+1:rs1,nu1+1:cs1),D2.d(ny2+1:rs2,nu2+1:cs2))];
      D1.Delay.Internal = [D1.Delay.Internal ; D2.Delay.Internal];
   end
   D1.Delay.Input = [D1.Delay.Input ; D2.Delay.Input];
   D1.Delay.Output = [D1.Delay.Output ; D2.Delay.Output];
   D1.StateInfo = [D1.StateInfo;D2.StateInfo];
   D1.SolverConfig = ltioptions.sparss.reconcile(D1.SolverConfig,D2.SolverConfig);
else
   % One system is NaN
   D1 = createGain(D1,NaN(iosize(D1)+iosize(D2)));
end