function G = modsum(H)
% Sums models in H.

%   Copyright 2023 The MathWorks, Inc.
nH = numel(H);
H1 = H(1);
if hasInternalDelay(H1)
   throw(ltipack.utNoDelaySupport('modsum',H1.Ts,'internal'))
end
[ny,nu] = size(H1.d);
Delay = H1.Delay;
u0 = H1.u0;  y0 = H1.y0;
EXPLICIT = true;

% Determine block sizes
bs = zeros(nH,1);
for ct=1:nH
   bs(ct) = size(H(ct).a,1);
end

% Allocate
nx = sum(bs);
a = zeros(nx);  e = eye(nx);
b = zeros(nx,nu);  c = zeros(ny,nx);  d = zeros(ny,nu); 
x0 = zeros(nx,1);  dx0 = zeros(nx,1);

% Populate
i = 0;
for ct=1:nH
   if ~isequal(H(ct).Delay,Delay)
      error(message('Control:transformation:modalsum1'))
   elseif ~(isequal(H(ct).u0,u0) && isequal(H(ct).y0,y0))
      error(message('Control:transformation:modalsum3'))      
   end
   idx = i+1:i+bs(ct);
   a(idx,idx) = H(ct).a;
   b(idx,:) = H(ct).b;
   c(:,idx) = H(ct).c;
   d = d + H(ct).d;
   if ~isempty(H(ct).e)
      EXPLICIT = false;
      e(idx,idx) = H(ct).e;
   end
   dx0(idx,:) = H(ct).dx0;
   x0(idx,:) = H(ct).x0;
   i = i+bs(ct);
end

% Build result
if EXPLICIT
   G = ltipack.ssdata(a,b,c,d,[],H1.Ts);
else
   G = ltipack.ssdata(a,b,c,d,e,H1.Ts);
end
G.u0 = u0;  G.y0 = y0;  G.dx0 = dx0;  G.x0 = x0;
G.Delay = Delay;
