function [A,B,C,D,E,dx0,x0,u0,y0,Delay] = parallelDF(sys1,sys2,t,p1,p2)
% Constructs data function for parallel IC of LTVSS/LPVSS models.

%   Copyright 2022-2024 The MathWorks, Inc.
if nargin<4
   [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
nx1 = size(A1,1);  nx2 = size(A2,1);
A = ltipack.util.blkdiag(A1,A2);
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
   [ny,nu] = size(D1);
   B = [B1 ; B2];
   C = [C1 C2];
   D = D1+D2;
   Delay = [];
   % Output offset
   if isempty(y01) && isempty(y02)
      y0 = [];
   else
      [y01,y02] = ltvpack.util.matchOffset(y01,ny,y02,ny);
      y0 = y01+y02;
   end
   % Reconcile input offsets
   if isempty(u01) && isempty(u02)
      u0 = [];
   else
      [u01,u02] = ltvpack.util.matchOffset(u01,nu,u02,nu);
      u0 = (u01+u02)/2;  du = (u02-u01)/2;
      % Note: Don't skip this when DU=0 (we can't guarantee u01=u02 for all t,p)
      [dx0,y0] = ltvpack.util.matchOffset(dx0,nx1+nx2,y0,ny);
      dx0 = dx0 +[B1*du;-B2*du];
      y0 = y0 + (D1-D2)*du;
   end
else
   % Note: All input and output delays not explicitly marked as absent (NaN)
   % are assumed varying and turned into internal delays.
   [Delay1,Delay2] = ltvpack.util.matchDelay(Delay1,size(D1),Delay2,size(D2));
   [B1,C1,D1,u01,y01,Delay1] = ltipack.foldDelay('ltv',...
      B1,C1,D1,u01,y01,Delay1,Delay1.Input,Delay1.Output);
   [B2,C2,D2,u02,y02,Delay2] = ltipack.foldDelay('ltv',...
      B2,C2,D2,u02,y02,Delay2,Delay2.Input,Delay2.Output);
   nu = numel(Delay1.Input);  nfd1 = numel(Delay1.Internal);
   ny = numel(Delay1.Output);  nfd2 = numel(Delay2.Internal);
   Delay = struct('Input',NaN(nu,1),'Output',NaN(ny,1),...
      'Internal',[Delay1.Internal;Delay2.Internal]);
   B = [[B1(:,1:nu) ; B2(:,1:nu)] , ...
      ltipack.util.blkdiag(B1(:,nu+1:nu+nfd1),B2(:,nu+1:nu+nfd2))];
   C = [[C1(1:ny,:) , C2(1:ny,:)] ;...
      ltipack.util.blkdiag(C1(ny+1:ny+nfd1,:),C2(ny+1:ny+nfd2,:))];
   D = [D1(1:ny,1:nu)+D2(1:ny,1:nu) , ...
      [D1(1:ny,nu+1:nu+nfd1),D2(1:ny,nu+1:nu+nfd2)] ; ...
      [D1(ny+1:ny+nfd1,1:nu);D2(ny+1:ny+nfd2,1:nu)] , ...
      ltipack.util.blkdiag(D1(ny+1:ny+nfd1,nu+1:nu+nfd1),D2(ny+1:ny+nfd2,nu+1:nu+nfd2))];
   % Output offset
   if isempty(y01) && isempty(y02)
      y0 = [];
   else
      [y01,y02] = ltvpack.util.matchOffset(y01,ny+nfd1,y02,ny+nfd2);
      y0 = [y01(1:ny,:)+y02(1:ny,:);y01(ny+1:ny+nfd1,:);y02(ny+1:ny+nfd2,:)];
   end
   % Reconcile input offsets
   if isempty(u01) && isempty(u02)
      u0 = [];
   else
      [uw01,uw02] = ltvpack.util.matchOffset(u01,nu+nfd1,u02,nu+nfd2);
      u01 = uw01(1:nu,:);  u02 = uw02(1:nu,:);
      u0 = [(u01+u02)/2 ; uw01(nu+1:nu+nfd1,:) ; uw02(nu+1:nu+nfd2,:)];
      du = (u02-u01)/2;
      % Note: Don't skip this when DU=0 (we can't guarantee u01=u02 for all t,p)
      [dx0,y0] = ltvpack.util.matchOffset(dx0,nx1+nx2,y0,ny+nfd1+nfd2);
      dx0 = dx0 +[B1(:,1:nu)*du;-B2(:,1:nu)*du];
      y0 = y0 + [D1(1:ny,1:nu)-D2(1:ny,1:nu) ; ...
         D1(ny+1:ny+nfd1,1:nu) ; -D2(ny+1:ny+nfd2,1:nu)] *du;
   end   
end