function [A,B,C,D,E,dx0,x0,u0,y0,Delay] = seriesDF(sys1,sys2,ScalarFlags,t,p1,p2)
% Constructs data function for series connection SYS1*SYS2 of LTVSS/LPVSS models.

%   Copyright 2022 The MathWorks, Inc.
if nargin<5
   [A1,B1,C1,D1,E1,dx01,x01,u01,y01,Delay1] = sys1.DataFunction(t);
   [A2,B2,C2,D2,E2,dx02,x02,u02,y02,Delay2] = sys2.DataFunction(t);
else
   [A1,B1,C1,D1,E1,dx01,x01,u01,y01,Delay1] = sys1.DataFunction(t,p1);
   [A2,B2,C2,D2,E2,dx02,x02,u02,y02,Delay2] = sys2.DataFunction(t,p2);
end

% Handle scalar flags
if ScalarFlags(1)
   [A1,B1,C1,D1,E1,dx01,x01,u01,y01,Delay1] = localDiagRep(size(sys2,1),...
      A1,B1,C1,D1,E1,dx01,x01,u01,y01,Delay1);
elseif ScalarFlags(2)
   [A2,B2,C2,D2,E2,dx02,x02,u02,y02,Delay2] = localDiagRep(size(sys1,2),...
      A2,B2,C2,D2,E2,dx02,x02,u02,y02,Delay2);
end
nx1 = size(A1,1);  nx2 = size(A2,1);
E = ltipack.util.blkdiagE(E1,E2,nx1,nx2);

% Combine state and state derivative offsets
if isempty(x01) && isempty(x02)
   x0 = [];
else
   [x01,x02] = ltvpack.util.matchOffset(x01,nx1,x02,nx2);
   x0 = [x01;x02];
end
if isempty(dx01) && isempty(dx02)
   dx0 = [];
else
   [dx01,dx02] = ltvpack.util.matchOffset(dx01,nx1,dx02,nx2);
   dx0 = [dx01;dx02];
end

if isempty(Delay1) && isempty(Delay2)
   % Delay-free case
   A = [A1 , B1*C2 ; zeros(nx2,nx1) , A2];
   B = [B1*D2 ; B2];
   C = [C1 , D1*C2];
   D = D1 * D2;
   u0 = u02;
   y0 = y01;
   if ~(isempty(u01) && isempty(y02))
      n = size(B1,2);
      [u01,y02] = ltvpack.util.matchOffset(u01,n,y02,n);
      aux = y02-u01;
      if isempty(dx0)
         dx0 = [B1*aux;zeros(nx2,1)];
      else
         dx0(1:nx1) = dx0(1:nx1)+B1*aux;
      end
      if isempty(y0)
         y0 = D1*aux;
      else
         y0 = y0+D1*aux;
      end
   end
   Delay = [];
else
   [Delay1,Delay2] = ltvpack.util.matchDelay(Delay1,size(D1),Delay2,size(D2));
   ny = numel(Delay1.Output);   nu = numel(Delay2.Input);
   % Combine and fold delays squeezed between SYS2 and SYS1
   fDelay = ltvpack.util.addDelay(Delay1.Input,Delay2.Output);
   n = numel(fDelay);
   % Fold into D2 for consistency with MTIMES
   [B2,C2,D2,u02,y02,Delay2] = ltipack.foldDelay('ltv',...
      B2,C2,D2,u02,y02,Delay2,NaN(nu,1),fDelay);
   Delay = struct('Input',Delay2.Input,'Output',Delay1.Output,...
      'Internal',[Delay1.Internal;Delay2.Internal]);
   nfd1 = numel(Delay1.Internal);   nfd2 = numel(Delay2.Internal);
   % Matrices
   A = [A1 , B1(:,1:n)*C2(1:n,:) ; zeros(nx2,nx1) , A2];
   aux = [B1(:,1:n)*D2(1:n,:) ; B2];
   B = [aux(:,1:nu) [B1(:,n+1:n+nfd1);zeros(nx2,nfd1)] aux(:,nu+1:nu+nfd2)];
   C = [C1 D1(:,1:n)*C2(1:n,:) ; zeros(nfd2,nx1) , C2(n+1:n+nfd2,:)];
   aux = [D1(:,1:n)*D2(1:n,:) ; D2(n+1:n+nfd2,:)];
   D = [aux(:,1:nu) [D1(:,n+1:n+nfd1);zeros(nfd2,nfd1)] aux(:,nu+1:nu+nfd2)];
   % Offsets
   if isempty(u01) && isempty(u02)
      u0 = [];
   else
      [u01,u02] = ltvpack.util.matchOffset(u01,n+nfd1,u02,nu+nfd2);
      u0 = [u02(1:nu,:) ; u01(n+1:n+nfd1,:) ; u02(nu+1:nu+nfd2)];
   end
   if isempty(y01) && isempty(y02)
      y0 = [];
   else
      [y01,y02] = ltvpack.util.matchOffset(y01,ny+nfd1,y02,n+nfd2);
      y0 = [y01 ; y02(n+1:n+nfd2)];
   end
   if ~(isempty(u01) && isempty(y02))
      [u01,y02] = ltvpack.util.matchOffset(u01,n+nfd1,y02,n+nfd2);
      aux = y02(1:n,:)-u01(1:n,:);
      if isempty(dx0)
         dx0 = [B1(:,1:n)*aux;zeros(nx2,1)];
      else
         dx0(1:nx1) = dx0(1:nx1)+B1(:,1:n)*aux;
      end
      if isempty(y0)
         y0 = [D1(:,1:n)*aux;zeros(nfd2,1)];
      else
         y0(1:ny+nfd1,:) = y0(1:ny+nfd1,:)+D1(:,1:n)*aux;
      end
   end
end
      
%-------------------
function [A,B,C,D,E,dx0,x0,u0,y0,Delay] = localDiagRep(N,A,B,C,D,E,dx0,x0,u0,y0,Delay)
% Constructs model of SYS*eye(N) for a SISO LTV/LPV model SYS
I = eye(N);
A = kron(I,A);
if ~isempty(E)
   E = kron(I,E);
end
if ~isempty(dx0)
   dx0 = repmat(dx0,[N 1]);
end
if ~isempty(x0)
   x0 = repmat(x0,[N 1]);
end

if isempty(Delay)
   % No internal delays
   B = kron(I,B);
   C = kron(I,C);
   D = kron(I,D);
   if ~isempty(u0)
      u0 = u0(ones(N,1),:);
   end
   if ~isempty(y0)
      y0 = y0(ones(N,1),:);
   end
else
   if isfield(Delay,'Internal') && ~isempty(Delay.Internal)
      [nyz,nuw] = size(D);
      B = [kron(I,B(:,1)) , kron(I,B(:,2:nuw))];
      C = [kron(I,C(1,:)) ; kron(I,C(2:nyz,:))];
      D = [kron(I,D(1,1)) kron(I,D(1,2:nuw)) ; ...
         kron(I,D(2:nyz,1)) kron(I,D(2:nyz,2:nuw))];
      Delay.Internal = repmat(Delay.Internal,[N 1]);
   else
      B = kron(I,B);
      C = kron(I,C);
      D = kron(I,D);
   end
   Delay.Input = repmat(Delay.Input,[N 1]);
   Delay.Output = repmat(Delay.Output,[N 1]);
   if ~isempty(u0)
      u0 = [u0(ones(N,1),:) ; repmat(u0(2:nuw,:),[N 1])];
   end
   if ~isempty(y0)
      y0 = [y0(ones(N,1),:) ; repmat(y0(2:nyz,:),[N 1])];
   end
end
