function [b,c,d,u0,y0,Delay] = foldDelay(CONTEXT,b,c,d,u0,y0,Delay,uDelay,yDelay)
% Turns subsets of the input and output delays into internal delays. 
% This function does the following:
%  1) For LTI models, all nonzero values in uDelay and yDelay are folded.
%  2) For LTV models, all non-NaN values in uDelay and yDelay are folded., whether zero or nonzero. This is necessary to 
%     ensure a constant number of internal delays.

%   Author(s): P. Gahinet
%   Copyright 1986-2023 The MathWorks, Inc.

% Decide which delays to fold
switch CONTEXT
   case 'lti'
      ifu = find(uDelay>0);
      ify = find(yDelay>0);
   case 'ltv'
      ifu = find(~isnan(uDelay));
      ify = find(~isnan(yDelay));
end
uDelay = uDelay(ifu,:);
yDelay = yDelay(ify,:);
nfu = numel(ifu);
nfy = numel(ify);

% Fold delays
if nfu>0 || nfy>0
   nx = size(b,1);
   [nyz,nuw] = size(d);
   % Augment w,z to [w;wu;wy] where wu = delayed u(ifu) and y(ify) = wy
   b = [b b(:,ifu) zeros(nx,nfy,'like',b)];  b(:,ifu) = 0;
   c = [c;zeros(nfu,nx,'like',c);c(ify,:)];  c(ify,:) = 0;
   d = [d d(:,ifu);zeros(nfu,nuw+nfu,'like',d);d(ify,:) d(ify,ifu)];
   d(:,ifu) = 0;  d(nyz+1:nyz+nfu,ifu) = eye(nfu,'like',d);  % zu = u(ifu)
   d(ify,:) = 0;  d(ify,end+1:end+nfy) = eye(nfy,'like',d);  % y(ify) = wy
   % Delays
   Delay.Input(ifu) = max(0,Delay.Input(ifu) - uDelay);
   Delay.Output(ify) = max(0,Delay.Output(ify) - yDelay);
   Delay.Internal = [Delay.Internal ; uDelay ; yDelay];
   % Offsets
   if ~isempty(u0)
      u0 = [u0 ; u0(ifu,:) ; zeros(nfy,1)];  u0(ifu) = 0;
   end
   if ~isempty(y0)
      y0 = [y0 ; zeros(nfu,1); y0(ify,:)];  y0(ify) = 0;
   end
end
