function [Dd,ICMap] = utDiscretizeIMP(Dc,Ts,options)
% Impulse-invariant discretization of affine state-space model:
%
%     dx/dt = dx0 + A (x-x0) + B [u-u0;w-w0] 
%     [y-y0;z-z0] = C (x-x0) + D [u-u0;w-w0] 
%     w(t) = z(t-tau)
% 
% Computes a discrete model whose response to a unit-area pulse of length 
% Ts and height 1/Ts matches the continuous-time response to the Dirac
% impulse u(t) = delta(t) * ej.
%
% ICMAP    Matrix mapping continuous-time initial conditions (xci,uci) to 
%          discrete-time initial condition xdi.
%
% Note: Assumes zero internal delays have already been eliminated.

%   Copyright 1986-2023 The MathWorks, Inc.

% Note: The state x[k] is x(k*Ts-) (so that x[0]=x0) but the output y[k] is
%       the post-impulse value y(k*Ts+).

% TOL: tolerance for comparing normalized delays (valued in [0,1])
tol = 1e4*eps;
CONSISTENT = (nargin>2 && strcmp(options.Consistency,'on'));

% Decompose external delays into entire + (normalized) fractional delays
[Delay,fid,fod,ffd] = ltipack.ssdata.discretizeDelay(Dc.Delay,Ts);
nu = length(fid);
ny = length(fod);
nfd = length(ffd);
nx = size(Dc.a,1);  % always preserved

% Discretization algorithm
if nfd>0
   % Model has internal delays
   if hasDelayDynamics(Dc,'pole')
      % General case: not supported
      error(message('Control:transformation:c2d04'))
   end
   % Exact discretization of models without state delays in state equation.
   % DDAE is then equivalent to ODE
   %     dx/dt = a x + sum bj u(t-theta_j) + f
   %     y = sum cj x(t-thetaj) + dj u(t-theta_j) + g
   if any(fid)
      % Turn residual fractional input delays into internal delays
      Dc = utFoldDelay(Dc,fid*Ts,[]);
   end
   [a,b1,b2,c1,c2,d11,d12,d21,d22] = getBlockData(Dc);
   tau = Dc.Delay.Internal;
   nfd = numel(tau);
   x0 = Dc.x0;  u0 = Dc.u0(1:nu,:); % absorb u0
   gw = (eye(nfd)-d22)\(Dc.y0(ny+1:end,:)-d21*u0-Dc.u0(nu+1:end,:)); % (I-D22)\(z0-w0)
   dx0 = Dc.dx0 - b1*u0 + b2*gw;
   y0 = Dc.y0(1:ny,:) - d11*u0 + d12*gw;

   if CONSISTENT || ~any([fid;fod;ffd])
      % Shortcut when all delays are multiple of Ts (there is a minimal
      % realization with NFD internal delays)
      [~,~,~,~,~,~,INFO] = smreal(d22,c2,b2,[]);
      i1 = find(INFO.xc);   n1 = numel(i1);
      i2 = find(~INFO.xc);  n2 = numel(i2);
      M = utScaledExpm(Ts*[[a dx0];zeros(1,nx+1)]);
      Ah = M(1:nx,1:nx);
      b1Ts = b1*Ts;  b2i2 = b2(:,i2);  c2i1 = c2(i1,:);
      % Note: Keep original order of internal delays to ensure
      % consistency in gridded LTVs.
      perm = [i2;i1]';  rp = [1:ny ny+perm];  cp = [1:nu nu+perm];
      Bd(:,cp) = [Ah*b1Ts Ah*b2i2 zeros(nx,n1)];
      Cd(rp,:) = [c1 ; zeros(n2,nx) ; c2i1];
      Dd(rp,cp) = [c1*b1Ts c1*b2i2 d12(:,i1);...
         d21(i2,:)*Ts d22(i2,i2) zeros(n2,n1);...
         c2i1*b1Ts c2i1*b2i2 d22(i1,i1)];
      Dd = ltipack.ssdata(Ah,Bd,Cd,Dd,[],Ts);
      Dd.dx0 = x0 + M(1:nx,nx+1);
      Dd.x0 = x0;
      Dd.u0 = zeros(nu+nfd,1);
      Dd.y0 = [y0;zeros(nfd,1)];
      if CONSISTENT
         Delay.Internal = max(1,Delay.Internal);
      end
      Dd.Delay = Delay;
   else
      % General case with fractional delays
      [I,J,P,N,rho,delta,n1,n2,n3,n4] = elimZW('imp',b2,c1,c2,d12,d21,d22,tau,Ts,fod);
      n = n1+n2+n3;
      ndelta = nnz(delta);  % excludes FOD terms
      zfod = (fod==0);

      % z[k] consists of the values
      %    D21(j(m),:) u[k]
      % for the first N1 delays (those contributing to B2*w(t)) and
      % the values
      %    C2(j(m),:) (x(k+)-x0) for rho(m)=0
      %    C2(j(m),:) (x(k+1-rho(m)+)-x0) for rho(m)>0
      % for the next N2+N3 delays.
      Nd = N(1:n);
      ix = n1+1:n1+n2;  Nd(ix) = Nd(ix) + (rho(ix)>0);
      ix = n1+n2+1:n;   Nd(ix) = Nd(ix) - (rho(ix)==0);

      % Initialize the incremental updates
      %    x(k+Ej+) = Aj*x[k]+B1j*u[k]+B2j*w[k]+fj
      % where x[k]=x(k-). The initial Aj,Bj,fj are for x(k+) and the
      % final values are for x[k+1]=x(k+1-).
      b2_I = b2(:,I(1:ndelta)) .* delta(1:ndelta).';
      Aj = eye(nx);  B1j = b1;  fj = zeros(nx,1);
      B2j = zeros(nx,n);
      m = find(rho(1:n1)==0);  B2j(:,m) = b2_I(:,m); % x(k-)->x(k+)

      % Initialize the known part of
      %   z[k] = z0 + C2j (x(k+)-x0) + D21j*u[k] + D22j w[k]
      C2j = zeros(n,nx);  D22j = zeros(n);
      D21j = [d21(J(1:n1),:) ; zeros(n2+n3,nu)];
      m = n1+find(rho(n1+1:n)==0);  c2jm = c2(J(m),:);
      C2j(m,:) = c2jm*Aj;  D21j(m,:) = c2jm*B1j;  D22j(m,:) = c2jm*B2j;
      z0 = zeros(n,1);

      % Initialize the known part of
      %   y[k/k+1] = y0 + C1j*(x(k+)-x0) + D11j*u[k] + D12j*w[k]
      d12_I = d12(:,I(1:ndelta)) .* delta(1:ndelta).';
      C1j = zeros(ny,nx);  D11j = zeros(ny,nu);  D12j = zeros(ny,n);
      fact = c1(zfod,:);
      C1j(zfod,:) = fact*Aj;  D11j(zfod,:) = fact*B1j;
      D12j(zfod,1:n1+n2) = [fact*B2j(:,1:n1) , d12_I(zfod,n1+1:n1+n2)];
      for m=n1+n2+1:n  % theta(m,p)>1
         D12j(P(m),m) = d12_I(P(m),m);
      end
      for m=n+1:ndelta
         if rho(m)==0 % theta(m,p)=1
            p = P(m);  aux = d12_I(p,m) * c2(J(m),:);
            C1j(p,:) = C1j(p,:) + aux*Aj;
            D11j(p,:) = D11j(p,:) + aux*B1j;
            D12j(p,:) = D12j(p,:) + aux*B2j;
         end
      end

      % Sort the integration events
      % Note: Skip 1-rho(m) event when rho(m)=0 or C2(j(m),:)=0
      E1 = rho(1:n1,:);  % input changes at E1
      E2 = 1-rho;
      zC2 = all(c2==0,2);  % Note: C1(p,:)~=0 by FOD selection
      SKIP2 = [zC2(J(1:ndelta)) ; false(n+n4-ndelta,1)] | rho==0;
      Events = unique([0;E1;E2(~SKIP2);1]);

      % Integrate state equation
      t0 = 0;
      T = [[a dx0];zeros(1,nx+1)];  % transition matrix
      for j=2:length(Events)
         t1 = Events(j);
         % Integrate state equation on [T0,T1]
         if t1>t0+tol
            h = (t1-t0)*Ts;
            M = utScaledExpm(h*T);
            Ah = M(1:nx,1:nx);
            % Update Aj,B1j,B2j,fj
            Aj = Ah * Aj;
            B1j = Ah * B1j;
            B2j = Ah * B2j;
            fj = Ah * fj + M(1:nx,nx+1);
            t0 = t1;
         end

         % Find impulses at t=T1 and record state jump
         ix = find(abs(E1-t1)<tol);  % in range [1,n1]
         B2j(:,ix) = b2_I(:,ix);  % these columns are zero until pulse triggers
         E1(ix) = -1;

         % Update rows of y(k+1+) and z[k]
         ix = find(abs(E2-t1)<tol & ~SKIP2);  % in range [1,n+n4]
         for ct=1:numel(ix)
            m = ix(ct);  p = P(m);
            if m>n
               % y[k+1] update
               if J(m)==0
                  % C1(p,:) x(k+1-FOD(p)+) term
                  fact = c1(p,:);
               else
                  % delta(m) D12(p,i(m)) C2(j(m),:) x(k+1-theta(m,p)+) term
                  % for theta(m,p)<1
                  fact = d12_I(p,m) * c2(J(m),:);
               end
               C1j(p,:) = C1j(p,:) + fact * Aj;
               D11j(p,:) = D11j(p,:) + fact * B1j;
               D12j(p,:) = D12j(p,:) + fact * B2j;
               y0(p) = y0(p) + fact * fj;
            else
               % z[k] update for C2(j(m),:)*x(k+1-rho(m)+) term
               c2jm = c2(J(m),:);
               C2j(m,:) = c2jm * Aj;
               D21j(m,:) = D21j(m,:) + c2jm * B1j;
               D22j(m,:) = c2jm * B2j;
               z0(m) = c2jm * fj;
            end
         end
         E2(ix) = -1;
      end

      Dd = ltipack.ssdata(Aj,[B1j*Ts B2j],[C1j;C2j],[D11j*Ts D12j;D21j*Ts D22j],[],Ts);
      Delay.Internal = Nd;
      Delay.Output(~zfod) = Delay.Output(~zfod)+1;
      Dd.dx0 = x0 + fj;
      Dd.x0 = x0;
      Dd.u0 = zeros(nu+n,1);
      Dd.y0 = [y0 ; z0];
      Dd.Delay = Delay;
      % Eliminate zero discrete internal delays
      Dd = elimZeroDelay(Dd);
   end

else
   % Model has only input and output delays (after eliminating zero internal delays)
   [a,b,c,d] = getABCD(Dc);
   % Input offset must be absorbed into dx0,y0 to correctly discretize
   % constant (step-like) terms.
   x0 = Dc.x0;  u0 = Dc.u0;
   dx0 = Dc.dx0-b*u0;
   y0 = Dc.y0-d*u0;

   % Perform IMP discretization
   fid = reshape(fid,[1 nu]);
   fod = reshape(fod,[1 ny]);
   zfid = (fid==0);
   zfod = (fod==0);
   if all(zfid) && all(zfod)
      % Delay-free case
      % Note: x(k+) = x[k] + B * u[k]
      %       y(k+) = y0 + C*x(k+) = y0 + C*(x[k] + B*u[k])
      M = utScaledExpm([[a dx0]*Ts; zeros(1,nx+1)]);
      Ah = M(1:nx,1:nx);
      bTs = b*Ts;
      Dd = ltipack.ssdata(Ah,Ah*bTs,c,c*bTs,[],Ts);
      Dd.dx0 = x0 + M(1:nx,nx+1);
      Dd.y0 = y0;
   else
      % Discretization with fractional input and output delays
      % Initialize the incremental formulas
      %    x(k+Ej+) = Aj*x[k]+Bj*u[k]+fj
      %    [Y_i(k+) or Y_i(k+1+)] = Cj*x[k]+Dj*u[k]+gj
      % where x[k]=x(k-). The initial Aj,Bj,fj are for x(k+) and the final
      % values are for x[k+1]=x(k+1-).
      Aj = eye(nx);  fj = zeros(nx,1);
      Bj = zeros(nx,nu);  Bj(:,zfid) = b(:,zfid);
      Cj = zeros(ny,nx);  Cj(zfod,:) = c(zfod,:);
      Dj = zeros(ny,nu);  Dj(zfod,:) = c(zfod,:)*Bj;
      gj = y0;

      % Sort the integration events
      E1 = fid;   E2 = 1-fod;
      Events = unique([0 E1 E2 1]);
      E2(zfod) = -1;  % to prevent output update at t=1

      % Piecewise integration over intervals [Events(j-1),Events(j)]
      % RE: To guard against events differing by o(eps), integrate
      % only when the step size is > TOL and treat nearby events
      % as a single event to synchronize updates
      t0 = 0;
      T = [[a dx0];zeros(1,nx+1)];  % transition matrix
      for j=2:length(Events)
         t1 = Events(j);
         % Integrate state equation on [T0,T1]
         if t1>t0+tol
            h = (t1-t0)*Ts;
            M = utScaledExpm(h*T);
            Ah = M(1:nx,1:nx);
            % Update Aj,Bj,fj
            Aj = Ah * Aj;
            Bj = Ah * Bj;
            fj = Ah * fj + M(1:nx,nx+1);
            t0 = t1;
         end

         % Find impulses at t=T1 and record state jump
         iu = find(abs(E1-t1)<tol);
         Bj(:,iu) = b(:,iu);  % these columns are zero until pulse triggers
         E1(iu) = -1;

         % Update rows of y(k+1+)
         iy = find(abs(E2-t1)<tol);
         ciy = c(iy,:);
         Cj(iy,:) = ciy*Aj;
         Dj(iy,:) = ciy*Bj;
         gj(iy) = gj(iy)+ciy*fj;
         E2(iy) = -1;
      end

      Dd = ltipack.ssdata(Aj,Bj*Ts,Cj,Dj*Ts,[],Ts);
      Dd.dx0 = x0 + fj;
      Dd.y0 = gj;
      % Apply z^-1 shift to all output channels with fractional delays
      Delay.Output(~zfod) = Delay.Output(~zfod)+1;
   end
   Dd.x0 = x0;
   Dd.Delay = Delay;
end

Dd.StateName = Dc.StateName;
Dd.StatePath = Dc.StatePath;
Dd.StateUnit = Dc.StateUnit;

% Mapping ICMAP from continuous initial conditions (xc0,uc0) to xd0
ICMap = [eye(nx) zeros(nx,nu)];
