function [I,J,P,N,rho,delta,n1,n2,n3,n4] = ...
   elimZW(METHOD,b2,c1,c2,d12,d21,d22,tau,Ts,fod)
%ELIMZW  Express W in terms of X,U when D22 is structurally nilpotent.
%
%   When D22 is nilpotent, w(t) has the explicit expression
%
%      w(t) = inv(I-Delta_tau*D22) * Delta_tau * v(t)
%             
%           = sum_m  delta(m)  e_i(m)  v_j(m) (t - theta(m))
%
%   where v(t) = C2 x(t) + D21 u(t).
%
%   ELIMZW collects all delays THETA(m) in this expansion as well as 
%     1) The fractional output delays (positive entries in FOD)
%     2) The combined delays THETA(m,p) = THETA(m)+FOD(p) arising in y(t)
%        for delayed output channels (FOD(p)>0).
%   Each delay is decomposed into its integral and fractional parts N and
%   RHO. Non-contributing delays are omitted.
%
%   Each delay has integration milestones t=k+rho (for terms changing the 
%   input through B2*w(t)) and/or t=k+1-rho (for terms impacting the outputs 
%   y or z). The following information is stored to facilitate these updates:
%      I,J     Vectors of i(m),j(m) values. I=J=0 for FOD delays.
%      P       Output channel. P=0 for THETA(m) delays.
%      DELTA   Vector of delta(m) values. DELTA=0 for FOD delays.

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

% Initialize with nonzero terms in Delta_tau * v(t)
nfd = length(tau);
Jnz = find(any([c2 d21],2));  % col
I = Jnz;
J = Jnz;
delta = ones(size(Jnz));
theta = tau(Jnz);

% Iteratively compute expansion of (Delta_tau*D22)^m * (Delta_tau*v)
I_m = I;  J_m = J;  delta_m = delta;  theta_m = theta;  
for m=1:nfd-1
   % Construct numel(theta)-by-nfd grid of new terms
   theta_m = theta_m + tau';            theta_m = theta_m(:);
   delta_m = delta_m .* d22(:,I_m)';    delta_m = delta_m(:);
   J_m = repmat(J_m,[nfd 1]);
   I_m = repmat(1:nfd,[numel(I_m) 1]);  I_m = I_m(:);
   inz = find(delta_m~=0);
   if isempty(inz)
      break
   else
      I_m = I_m(inz,:);
      J_m = J_m(inz,:);
      theta_m = theta_m(inz,:);
      delta_m = delta_m(inz,:);
      I = [I ; I_m]; %#ok<*AGROW>
      J = [J ; J_m];
      theta = [theta ; theta_m];
      delta = [delta ; delta_m];
   end
end

% Remove terms that do not contribute to [b2;d12]*w(t)
nzI = any([b2;d12],1);
inz = find(nzI(I));
I = I(inz,:);
J = J(inz,:);
theta = theta(inz,:);
delta = delta(inz,:); % all nonzero

% Decompose THETA delays           
[N,rho] = ltipack.util.splitDelay(theta,Ts);

% XDOT flags terms that contribute to the state equation.
% Y flags terms that contribute to the output equation.
nzB2 = any(b2,1)';
nzC2 = any(c2,2);
nzD21 = any(d21,2);
XDOT = nzB2(I) & nzD21(J);  % B2(:,i(m)) * D21(j(m),:)) nonzero
zfod = (fod==0);
if METHOD(1)=='i'
   Y = any(d12(zfod,I),1)' & nzC2(J);  % D12(zfod,i(m)) * C2(j(m),:) nonzero
else
   Y = any(d12(zfod,I),1)';  % D12(zfod,i(m)) * [C2 D21](j(m),:) nonzero
end

% Add contribution of fractional output delays (FOD) 
nt = numel(theta);
p = find(~zfod);  % fod>0
if isempty(p)
   P = zeros(nt,1);
else
   fod = fod(p,:);
   theta_mp = fod + (N+rho)';      theta_mp = theta_mp(:);
   P = repmat(p,[1 nt]);           P = P(:);
   M = repmat(1:nt,[numel(p) 1]);  M = M(:);

   % Decompose THETA(m,p) and add surviving FOD delays
   [Ny,rhoy] = ltipack.util.splitDelay([theta_mp;fod],1);

   % Merge with THETA set and update contribution flags
   XDOT = [XDOT ; false(numel(Ny),1)];
   D12pI = (d12(p,I)~=0);
   if METHOD(1)=='i'
      Y = [Y ; D12pI(:) & nzC2(J(M)) ; any(c1(p,:),2)];
   else
      Y = [Y ; D12pI(:) ; any(c1(p,:),2)];
   end
   pzero = zeros(size(p));
   I = [I ; I(M) ; pzero];
   J = [J ; J(M) ; pzero];
   P = [zeros(nt,1) ; P ; p];
   N = [N ; Ny];  rho = [rho ; rhoy];
   delta = [delta ; delta(M) ; pzero];
end

% Remove non-contributing terms
inz = find(XDOT | Y);
I = I(inz,:);  J = J(inz,:);  P = P(inz,:);  
N = N(inz,:);  rho = rho(inz,:);  delta = delta(inz,:);
XDOT = XDOT(inz,:);  Y = Y(inz,:);

% Sort terms according to METHOD requirements.
switch METHOD
   case {'zoh','foh'}
      %  1) First N1 entries are THETA delays contributing to B2*w(t) 
      %     (with constant or changing value/slope). The delays with 
      %     rho>0 appear first.
      %  2) Next N2 entries are THETA delays that don't contribute to B2*w(t)
      %  3) Next N3 entries are the delays THETA(m,p)>1
      %  4) Last N4 entries are the remaining THETA(m,p) delays and the FOD
      %     delays (all <=1). The FOD delays appear last.
      F1 = (P==0 & XDOT);
      F2 = (N+rho>1);
      is1 = [find(F1 & rho>0);find(F1 & rho==0)];
      is2 = find(P==0 & ~XDOT);
      is3 = find(P>0 & F2);   is4 = find(P>0 & ~F2);
   case 'imp'
      %  1) First N1 entries are THETA delays contributing a pulse in B2*w(t)
      %  2) Next N2 entries are THETA delays contributing to Y
      %  3) Next N3 entries are the delays THETA(m,p)>1
      %  4) Last N4 entries are the remaining THETA(m,p) delays and the FOD
      %     delays (all <=1). The FOD delays appear last.
      F = (N+rho>1);
      is1 = find(P==0 & XDOT);  is2 = find(P==0 & Y);
      is3 = find(P>0 & F);      is4 = find(P>0 & ~F);
end

is = [is1;is2;is3;is4];
I = I(is);  J = J(is);  P = P(is);
N = N(is);  rho = rho(is);   delta = delta(is);
n1 = numel(is1); n2 = numel(is2); n3 = numel(is3); n4 = numel(is4);
