function [A,B,C,D,E,dx0,x0,u0,y0,Delay] = ...
   feedbackDF(sys1,sys2,indu,indy,sign,t,p1,p2)
% Data function for feedback interconnection of LTVSS/LPVSS models.
% Note: Expects INDU/INDY to be row vectors.

%   Copyright 2022 The MathWorks, Inc.
if nargin<7
   [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);  nx = nx1+nx2;
nu = numel(indu);
ny = numel(indy);

% Combine 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

% Handle delays
if isempty(Delay1) && isempty(Delay2)
   Delay = [];   nfd2 = 0;
else
   [Delay1,Delay2] = ltvpack.util.matchDelay(Delay1,size(D1),Delay2,size(D2));
   % Fold I/O delays of SYS1 in feedback path
   uDelay = NaN(size(Delay1.Input));
   uDelay(indu,:) = Delay1.Input(indu,:);
   yDelay = NaN(size(Delay1.Output));
   yDelay(indy,:) = Delay1.Output(indy,:);
   [B1,C1,D1,u01,y01,Delay1] = ltipack.foldDelay('ltv',...
      B1,C1,D1,u01,y01,Delay1,uDelay,yDelay);
   Delay1.Input(indu,:) = NaN;   % now always zero
   Delay1.Output(indy,:) = NaN;
   % Fold all I/O delays of SYS2
   [B2,C2,D2,u02,y02,Delay2] = ltipack.foldDelay('ltv',...
      B2,C2,D2,u02,y02,Delay2,Delay2.Input,Delay2.Output);
   % Resulting delays
   Delay = struct('Input',Delay1.Input,'Output',Delay1.Output,...
      'Internal',[Delay1.Internal;Delay2.Internal]);
   nfd2 = numel(Delay2.Internal);
end

% Input and output offsets
[nyz1,nuw1] = size(D1);
[nyz2,nuw2] = size(D2);
if isempty(u01) && isempty(u02)
   u0 = [];
else
   [u01,u02] = ltvpack.util.matchOffset(u01,nuw1,u02,nuw2);
   u0 = [u01;u02(ny+1:nuw2,:)];
end
if isempty(y01) && isempty(y02)
   y0 = [];
else
   [y01,y02] = ltvpack.util.matchOffset(y01,nyz1,y02,nyz2);
   y0 = [y01;y02(nu+1:nyz2,:)];
end

% Compute realization for feedback interconnection (see notes)
E = ltipack.util.blkdiagE(E1,E2,nx1,nx2);
M = [eye(ny) -D1(indy,indu) ; -D2(1:nu,1:ny) sign*eye(nu)];
if nfd2>0
   R = matlab.internal.math.nowarn.mrdivide(...
      [zeros(nx1,ny) B1(:,indu);B2(:,1:ny) zeros(nx2,nu);...
      zeros(nyz1,ny) D1(:,indu);D2(nu+1:nyz2,1:ny) zeros(nfd2,nu)],M);
   aux = R * ltipack.util.blkdiag(C1(indy,:),C2(1:nu,:));
   A = ltipack.util.blkdiag(A1,A2) + aux(1:nx,:);
   C = ltipack.util.blkdiag(C1,C2(nu+1:nyz2,:)) + aux(nx+1:end,:);
   aux = R * ltipack.util.blkdiag(D1(indy,:),D2(1:nu,ny+1:nuw2));
   B = ltipack.util.blkdiag(B1,B2(:,ny+1:nuw2)) + aux(1:nx,:);
   D = ltipack.util.blkdiag(D1,D2(nu+1:nyz2,ny+1:nuw2)) + aux(nx+1:end,:);
else
   R = matlab.internal.math.nowarn.mrdivide(...
      [zeros(nx1,ny) B1(:,indu);B2 zeros(nx2,nu);zeros(nyz1,ny) D1(:,indu)],M);
   aux = R * ltipack.util.blkdiag(C1(indy,:),C2);
   A = ltipack.util.blkdiag(A1,A2) + aux(1:nx,:);
   C = [C1 , zeros(nyz1,nx2)] + aux(nx+1:end,:);
   aux = R * [D1(indy,:) ; zeros(nu,nuw1)];
   B = [B1 ; zeros(nx2,nuw1)] + aux(1:nx,:);
   D = D1 + aux(nx+1:end,:);
end

% Finalize offsets
% fs0 = [y01(indy)-u02(1:ny) ; y02(1:nu)]
fs0 = zeros(ny+nu,1);
if ~isempty(y0)
   fs0 = fs0 + [y01(indy,:) ; y02(1:nu,:)];
end
if ~isempty(u0)
   fs0 = fs0 - [u02(1:ny,:) ; zeros(nu,1)];
end
if any(fs0)
   aux = R*fs0;
   [dx0,y0] = ltvpack.util.matchOffset(dx0,nx,y0,nyz1+nfd2);
   dx0 = dx0 + aux(1:nx,:);
   y0 = y0 + aux(nx+1:end,:);
end
