function D1 = plus(D1,D2)
% Adds two FINITE state-space models.

%   Copyright 1986-2020 The MathWorks, Inc.
[ny,nu] = iosize(D1);

% Equalize input and output delays
if hasdelay(D1) || hasdelay(D2)
   [D1,D2] = plusDelay(D1,D2,...
      struct('Input',~any([D1.b(:,1:nu);D1.d(:,1:nu)],1),...
      'Output',~any([D1.c(1:ny,:),D1.d(1:ny,:)],2)),...
      struct('Input',~any([D2.b(:,1:nu);D2.d(:,1:nu)],1),...
      'Output',~any([D2.c(1:ny,:),D2.d(1:ny,:)],2)));
end

% Sizes
nfd1 = length(D1.Delay.Internal);
nfd2 = length(D2.Delay.Internal);
nx1 = size(D1.a,1);
nx2 = size(D2.a,1);

% Equalize input offsets by absorbing deviation from average into dx0,y0
u10 = D1.u0(1:nu,:);
u20 = D2.u0(1:nu,:);
u0 = (u10+u20)/2;
du = (u20-u10)/2;

% Compute realization for sum
D1.a = ltipack.util.blkdiag(D1.a,D2.a);
D1.e = ltipack.util.blkdiagE(D1.e,D2.e,nx1,nx2);
D1.dx0 = [D1.dx0;D2.dx0];
D1.x0 = [D1.x0;D2.x0];
if nfd1==0 && nfd2==0
   % No internal delays
   D1.u0 = u0;
   D1.y0 = D1.y0+D2.y0;
   if any(du)
      D1.dx0 = D1.dx0 + [D1.b*du ; -D2.b*du];
      D1.y0 = D1.y0 + (D1.d-D2.d)*du;
   end
   D1.b = [D1.b ; D2.b];
   D1.c = [D1.c , D2.c];
   D1.d = D1.d + D2.d;
else
   D1.u0 = [u0;D1.u0(nu+1:nu+nfd1,:);D2.u0(nu+1:nu+nfd2,:)];
   D1.y0 = [D1.y0(1:ny,:)+D2.y0(1:ny,:);D1.y0(ny+1:ny+nfd1,:);D2.y0(ny+1:ny+nfd2,:)];
   if any(du)
      D1.dx0 = D1.dx0 + [D1.b(:,1:nu)*du ; -D2.b(:,1:nu)*du];
      D1.y0 = D1.y0 + [D1.d(1:ny,1:nu)-D2.d(1:ny,1:nu) ; ...
         D1.d(ny+1:ny+nfd1,1:nu) ; -D2.d(ny+1:ny+nfd2,1:nu)] * du;
   end
   D1.b = [[D1.b(:,1:nu) ; D2.b(:,1:nu)] , ...
      ltipack.util.blkdiag(D1.b(:,nu+1:nu+nfd1),D2.b(:,nu+1:nu+nfd2))];
   D1.c = [[D1.c(1:ny,:) , D2.c(1:ny,:)] ;...
      ltipack.util.blkdiag(D1.c(ny+1:ny+nfd1,:),D2.c(ny+1:ny+nfd2,:))];
   D1.d = [D1.d(1:ny,1:nu)+D2.d(1:ny,1:nu) , ...
      [D1.d(1:ny,nu+1:nu+nfd1),D2.d(1:ny,nu+1:nu+nfd2)] ; ...
      [D1.d(ny+1:ny+nfd1,1:nu);D2.d(ny+1:ny+nfd2,1:nu)] , ...
      ltipack.util.blkdiag(D1.d(ny+1:ny+nfd1,nu+1:nu+nfd1),D2.d(ny+1:ny+nfd2,nu+1:nu+nfd2))];
   D1.Delay.Internal = [D1.Delay.Internal ; D2.Delay.Internal];
end
D1.Scaled = false;
