function [Dd,ICMap] = utDiscretizeFOH(Dc,Ts,options)
%
%     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)
%
%   ICMAP    Matrix mapping continuous-time initial conditions (xci,uci) to 
%            discrete-time initial condition xdi.

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

% TOLINT: tolerance for comparing delay times
tol = 1e4*eps;
CONSISTENT = (nargin>2 && strcmp(options.Consistency,'on'));
USESTATE = (nargin>2 && strcmp(options.DelayModeling,'state'));

% Decompose external delays into entire + 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
   DELAYDYN = hasDelayDynamics(Dc,'pole');
   if DELAYDYN || CONSISTENT || ~any([fid;fod;ffd])
      % General case: Discretize rational part H(s) of delay model.
      % Note: This is exact for explicit DDAEs when all delays are multiple of Ts
      if DELAYDYN && ~CONSISTENT
         warning(message('Control:transformation:C2dApproximate'))
      end

      % Absorb fractional internal delays into H(s)
      Dc.Delay.Internal = zeros(0,1);
      Dc.Delay.Output = [fod*Ts ; ffd*Ts];
      Dc.Delay.Input = [fid*Ts ; zeros(nfd,1)];
      
      % Discretize H(s)
      [Dd,ICMap] = utDiscretizeFOH(Dc,Ts);
      Dd.Delay.Internal = [...
         Delay.Internal + Dd.Delay.Input(nu+1:end,:) + Dd.Delay.Output(ny+1:end,:) ; ...
         Dd.Delay.Internal];
      Dd.Delay.Output = Delay.Output + Dd.Delay.Output(1:ny,:);
      Dd.Delay.Input = Delay.Input + Dd.Delay.Input(1:nu,:);
      ICMap = ICMap(:,1:nx+nu);
   else
      % 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);

      % Compute expansion 
      %    w(t)-w0 = gw + sum_m delta(m) e_i(m) v_j(m)(t-theta(m))
      %       v(t) = C2*(x-x0) + D21*(u-u0)
      % and collect all delays involved in w(t), FOD, and w(t-FOD(p)).
      [I,J,P,N,rho,delta,n1,n2,n3,n4] = elimZW('foh',b2,c1,c2,d12,d21,d22,tau,Ts,fod);
      n12 = n1+n2;  n = n12+n3;   ndelta = nnz(delta);  % excludes FOD terms
      gw = (eye(nfd)-d22)\(Dc.y0(ny+1:end,:)-Dc.u0(nu+1:end,:)); % (I-D22)\(z0-w0)
      dx0 = Dc.dx0 + b2*gw;
      x0 = Dc.x0;
      u0 = Dc.u0(1:nu,:);
      y0 = Dc.y0(1:ny,:) + d12*gw;
      zfod = (fod==0);

      % w[k] consists of
      %   C2(j(m),:) x(k-theta(m)),      m=1:n12
      %   C2(j(m),:) x(k+1-theta(m,p)),  m=n12+1:n
      %   D21(j(m),:) u[k-N(m)+1],       m=1:n
      %   D21(j(m),:) u[k-N(m)],         m=1:n
      %   D21(j(m),:) u[k-N(m)-1],       m=1:n
      % The first N entries of z[k] consist of
      %   C2(j(m),:) (x[k]-x0) for rho(m)=0
      %   C2(j(m),:) (x(k+1-rho(m))-x0) for rho(m)>0
      % and are delayed by Ndx. The remaining 3*N entries D21(j(m),:)*u[k]
      % are delayed by Ndu. Initialize the known part of
      %   z[k] = z0 + C2j (x[k]-x0) + D21j*(u[k]-u0) + D22j w[k] + D23j (u[k+1]-u0)
      Ndx = N(1:n) + (rho(1:n)>0);  % delays for first N entries of z[k]
      Ndx(n12+1:n) = Ndx(n12+1:n) - 1;
      Ndu = [N(1:n)-1 ; N(1:n) ; N(1:n)+1];
      C2j = zeros(4*n,nx);  D22j = zeros(4*n);  D23j = zeros(4*n,nu);
      m = find(rho(1:n)==0);  C2j(m,:) = c2(J(m),:);
      D21j = [zeros(n,nu) ; repmat(d21(J(1:n),:),[3 1])];
      z0 = zeros(4*n,1);

      % Initialize the known part of
      %   y[k/k+1] = y0 + C1j*(x[k]-x0) + D11j*(u[k]-u0) + D12j*w[k] + D13j*(u[k+1]-u0)
      d12_I = d12(:,I(1:ndelta)) .* delta(1:ndelta).';
      C1j = zeros(ny,nx);  D11j = zeros(ny,nu);  D12j = zeros(ny,4*n);  D13j = zeros(ny,nu);
      % No output delay
      C1j(zfod,:) = c1(zfod,:);  D11j(zfod,:) = d11(zfod,:);
      for m=1:n12
         % w(m)[k], u[k-N(m)], u[k-N(m)-1]
         D12j(zfod,[m 2*n+m 3*n+m]) = [1 1-rho(m) rho(m)] .* d12_I(zfod,m);
      end
      % With output delay
      ip = find(~zfod);
      D11j(ip,:) = fod(ip,:).*d11(ip,:);     % rho(m) D11(j(m),:) (u[k]-u0)
      D13j(ip,:) = (1-fod(ip,:)).*d11(ip,:); % (1-rho(m)) D11(j(m),:) (u[k+1]-u0)
      for m=n12+1:n % theta(m,p)>1
         D12j(P(m),[m n+m 2*n+m]) = [1 1-rho(m) rho(m)] * d12_I(P(m),m);
      end
      for m=n+1:ndelta % theta(m,p)<=1
         p = P(m);  j = J(m);
         if rho(m)==0 % theta(m,p)=1
            C1j(p,:) = C1j(p,:) + d12_I(p,m) * c2(j,:);
            D11j(p,:) = D11j(p,:) + d12_I(p,m) * d21(j,:);
         else
            D11j(p,:) = D11j(p,:) + rho(m) * d12_I(p,m) * d21(j,:);            
            D13j(p,:) = D13j(p,:) + (1-rho(m)) * d12_I(p,m) * d21(j,:);            
         end
      end

      % Update the relations
      %    x(k+Ej)-x0   = Aj (x[k]-x0) + B1j (u[k]-u0) + B2j w[k] + B3j (u[k+1]-u0)
      %    [u;sigma](k+Ej)  = [u;sigma](k) + R1j (u[k]-u0) + R2j w[k] + R3j (u[k+1]-u0)
      %    Ts d[u;sigma]/dt = H1j (u[k]-u0) + H2j w[k] + H3j (u[k+1]-u0)
      % while integrating the state equation over one sampling interval.
      % Note: These updates do not depend on the first n entries of w[k].
      nusig = nu+nfd;
      fj = zeros(nx,1);
      Aj = eye(nx); B1j = zeros(nx,nu); B2j = zeros(nx,4*n); B3j = zeros(nx,nu);
      R1j = [eye(nu);zeros(nfd,nu)];  R2j = zeros(nusig,4*n);  R3j = zeros(nusig,nu);
      H1j = [-eye(nu);zeros(nfd,nu)];  H2j = zeros(nusig,4*n);  H3j = [eye(nu);zeros(nfd,nu)];
      for m=1:n1
         ir = nu+I(m);
         if rho(m)>0
            R2j(ir,[2*n+m 3*n+m]) = [1-rho(m) rho(m)] * delta(m);
            H2j(ir,[2*n+m 3*n+m]) = [delta(m) -delta(m)]; % u[k-N(m)],u[k-N(m)-1]
         else
            R2j(ir,2*n+m) = delta(m);
            H2j(ir,[n+m 2*n+m]) = [delta(m) -delta(m)]; % u[k+1-N(m)],u[k-N(m)]
         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,:);  % slope changes at E1>0
      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]);
      t0 = 0;
      T = [a b1 b2 zeros(nx,nusig) dx0 ;...
         zeros(nusig,nx+nusig) eye(nusig)/Ts zeros(nusig,1);....
         zeros(nusig+1,nx+2*nusig+1)];  % transition matrix
      for j=2:numel(Events)
         t1 = Events(j);  % > 0
         % Integrate state equation on [T0,T1]
         if t1>t0+tol
            h = (t1-t0)*Ts;
            ehT = utScaledExpm(h*T);
            Ah = ehT(1:nx,1:nx);  Bh = ehT(1:nx,nx+1:nx+nusig);  
            Sh = ehT(1:nx,nx+nusig+1:nx+2*nusig);  fh = ehT(1:nx,end);   
            Rh = ehT(nx+1:nx+nusig,nx+nusig+1:nx+2*nusig);
            % Update Aj,Bj,fj
            Aj = Ah * Aj;
            B1j = Ah * B1j + Bh * R1j + Sh * H1j;
            B2j = Ah * B2j + Bh * R2j + Sh * H2j;
            B3j = Ah * B3j + Bh * R3j + Sh * H3j;
            fj = Ah * fj + fh;
            R1j = R1j + Rh * H1j;
            R2j = R2j + Rh * H2j;
            R3j = R3j + Rh * H3j;
            t0 = t1;
         end

         % Find terms causing B2*(w(t)-w0) to change slope at t=T1 and update Hj
         ix = find(abs(E1-t1)<tol);  % in range [1,n1]
         for ct=1:numel(ix)
            % u[k-Nm]-u[k-Nm-1] -> u[k-Nm+1]-u[k-Nm] for rho(m)>0
            m = ix(ct);
            ir = nu+I(m);
            H2j(ir,[n+m 2*n+m 3*n+m]) = [delta(m) -delta(m) 0];
         end
         E1(ix) = -1;  % make sure each delay is processed only once

         % Find y,z entries updated at t=T1
         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;
               D13j(p,:) = D13j(p,:) + fact * B3j;
               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,:) = D22j(m,:) + c2jm * B2j;
               D23j(m,:) = D23j(m,:) + c2jm * B3j;
               z0(m) = c2jm * fj;
            end
         end
         E2(ix) = -1;
      end

      % Move all u[k],u[k+1] terms out of w[k] (ensures Ndu>0)
      m0 = find(N(1:n)==0);
      m1 = find(N(1:n)==1);
      B3j = B3j + B2j(:,n+m0) * d21(J(m0),:);   
      B1j = B1j + B2j(:,[n+m1;2*n+m0]) * d21(J([m1;m0]),:);  
      B2j(:,[n+m0;n+m1;2*n+m0]) = 0;
      D13j = D13j + D12j(:,n+m0) * d21(J(m0),:);
      D11j = D11j + D12j(:,[n+m1;2*n+m0]) * d21(J([m1;m0]),:);
      D12j(:,[n+m0;n+m1;2*n+m0]) = 0;
      D23j = D23j + D22j(:,n+m0) * d21(J(m0),:);
      D21j = D21j + D22j(:,[n+m1;2*n+m0]) * d21(J([m1;m0]),:);
      D22j(:,[n+m0;n+m1;2*n+m0]) = 0;

      % Trim and optimize delay model
      Ju = repmat(J(1:n),[3 1]);
      [B2j,D12j,C2j,D21j,D22j,D23j,z0,Nd] = squeezeZW(...
         B2j,D12j,C2j,D21j,D22j,D23j,z0,Ndx,Ndu,Ju,n,d21);
      nfdd = numel(Nd);
      u0 = [u0 ; zeros(nfdd,1)];
      y0 = [y0 ; z0];

      % Eliminate u[k+1] from the state equation by change of state variable
      %   zeta[k] = x[k] - B3j * u[k]
      ICMap = [eye(nx) -B3j];
      B1j = B1j + Aj*B3j;
      D11j = D11j + C1j*B3j;
      D21j = D21j + C2j*B3j;
      x0 = x0-B3j*u0(1:nu,:);  % zeta0
      dx0 = x0 + fj;

      % Eliminate u[k+1] from the z equation.
      % 1) Replace occurrences delayed by N=1 by u[k] term in state/output equations
      nzd23 = any(D23j,2);
      ix = find(nzd23 & Nd==1);
      B1j = B1j + B2j(:,ix) * D23j(ix,:);
      D11j = D11j + D12j(:,ix) * D23j(ix,:);
      % 2) For occurrences delayed by N>1, obtain wm[k] as 
      %       wix[k] + D23j(m,:)*u[k] 
      %    delayed by N-1 samples, where wix[k] is zm[k] stripped of u[k+1]
      %    term and delayed by 1 sample.
      ix = find(nzd23 & Nd>1);  nix = numel(ix);
      if nix>0
         B2j = [B2j zeros(nx,nix)];   D12j = [D12j zeros(ny,nix)];
         C2j = [C2j ; C2j(ix,:)];    C2j(ix,:) = 0;
         D21j = [D21j ; D21j(ix,:)];   D21j(ix,:) = D23j(ix,:);
         D22j = [[D22j ; D22j(ix,:)] zeros(nfdd+nix,nix)];
         D22j(ix,1:nfdd+nix) = [zeros(nix,nfdd) eye(nix)];  % wix[k]+D23j(m,:)*u[k]
         Nd(ix) = Nd(ix)-1; % decrement delay
         Nd = [Nd ; ones(nix,1)];
         u0 = [u0 ; zeros(nix,1)];
         y0 = [y0 ; zeros(nix,1)];
         nfdd = nfdd+nix;
      end

      % Eliminate u[k+1] terms from output equation
      % * When D13j(p,:)=0, obtain yp[k] from yp[k+1] with unit output delay
      % * When D13j(p,:) nonzero, use internal delay
      aux = any(D13j,2);
      ix = find(~(zfod | aux));
      Delay.Output(ix) = Delay.Output(ix)+1;
      ix = find(aux);  nix = numel(ix);
      if nix>0
         % y(ix)[k] = gj(ix) + D13j(ix,:) (u[k]-u0) + wix[k]
         % Note: Nonzero rows of D13j are a subset of delayed outputs
         B2j = [B2j zeros(nx,nix)];
         C2j = [C2j ; C1j(ix,:)];   C1j(ix,:) = 0;
         D21j = [D21j ; D11j(ix,:)];   D11j(ix,:) = D13j(ix,:);
         D22j = [[D22j ; D12j(ix,:)] zeros(nfdd+nix,nix)];
         D12j(ix,1:nfdd+nix) = [zeros(nix,nfdd) eye(nix)];  % wix[k]
         Nd = [Nd ; ones(nix,1)];
         u0 = [u0 ; zeros(nix,1)];
         y0 = [y0 ; zeros(nix,1)];
      end

      Dd = ltipack.ssdata(Aj,[B1j B2j],[C1j;C2j],[D11j D12j;D21j D22j],[],Ts);
      Dd.dx0 = dx0;   Dd.x0 = x0;   Dd.u0 = u0;   Dd.y0 = y0;
      Delay.Internal = Nd;
      Dd.Delay = Delay;
   end

   % Eliminate zero discrete internal delays
   % Note: Internal delays < Ts give rise to zero discrete delays.
   if ~CONSISTENT
      Dd = elimZeroDelay(Dd);
      if ~isempty(Dd.e)
         error(message('Control:transformation:c2d12')) % non causal
      end
   end

else
   % Only input and output delays (after eliminating zero internal delays)
   [a,b,c,d] = getABCD(Dc);
   nxu = nx+nu;
   dx0 = Dc.dx0;  x0 = Dc.x0;  u0 = Dc.u0;  y0 = Dc.y0;

   % Perform ZOH 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.
      M = utScaledExpm(...
         [a*Ts b*Ts zeros(nx,nu) dx0*Ts ; ...
         zeros(nu,nx+nu)  eye(nu) zeros(nu,1) ; ...
         zeros(nu+1,nx+2*nu+1)]);
      Ah = M(1:nx,1:nx);     Bh = M(1:nx,nx+1:nxu);
      Sh = M(1:nx,nxu+1:nxu+nu);  fh = M(1:nx,end);
      Dd = ltipack.ssdata(Ah,Bh+Ah*Sh-Sh,c,d+c*Sh,[],Ts);
      % New state is zeta[k] = x[k]-Sh*u[k]
      x0 = x0 - Sh*u0;
      dx0 = x0 + fh;
      ICMap = [eye(nx) , -Sh];
      nw = 0;
   else
      % Discretization with fractional delays. Propagate the equations
      %    x(k+Ej) = Aj xk + Bj [uk;wk;sk] + fj
      %    v(k+Ej) = Rj [uk;wk;sk]
      %    Ts dv/dt = Hj [uk;wk;sk]
      %    y[k+(fod>0)] = Cj xk + Dj [uk;wk;sk] + gj
      % where wk = {uj[k-1],betaj>0} and sk = u[k+1]-u[k]
      jdelay = find(~zfid);  % delayed input channels
      jzero = find(zfid);    % input channels without delays
      nw = numel(jdelay);    % number of nonzero fractional input delays
      Aj = eye(nx);  Bj = zeros(nx,2*nu+nw);  fj = zeros(nx,1);
      Rj = [diag(1-fid) zeros(nu,nw+nu)];   Rj(~zfid,nu+1:nu+nw) = diag(fid(~zfid));
      Hj = zeros(nu,2*nu+nw);
      Hj(jdelay,jdelay) = eye(numel(jdelay));
      Hj(jdelay,nu+1:nu+nw) = -eye(numel(jdelay));
      Hj(jzero,nu+nw+jzero) = eye(numel(jzero));
      Cj = c;  Dj = d *Rj;   gj = Dc.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)]
      t0 = 0;
      T = [a b zeros(nx,nu) Dc.dx0 ; ...
         zeros(nu,nx+nu)  eye(nu)/Ts zeros(nu,1) ; ...
         zeros(nu+1,nx+2*nu+1)];  % transition matrix
      for j=2:numel(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);  Bh = M(1:nx,nx+1:nxu);  Sh = M(1:nx,nxu+1:nxu+nu);  
            fh = M(1:nx,end);   Rh = M(nx+1:nxu,nxu+1:nxu+nu);
            % Update Aj,B1j,...
            Aj = Ah * Aj;
            Bj = Ah * Bj + Bh * Rj + Sh * Hj;
            Rj = Rj + Rh * Hj;
            fj = Ah * fj + fh;
            t0 = t1;
         end

         % Find slopes updated at t=T1 (within TOL)
         iu = find(abs(E1-t1)<tol);
         if ~isempty(iu)
            Hj(iu,1:nu+nw) = 0;
            Hj(iu,nu+nw+iu) = eye(numel(iu));
            E1(iu) = -1;  % make sure each delay is processed only once
         end

         % Find delayed outputs updated at t=T1 and update Cj,D1j,.. accordingly
         %   y(T1) = y0+c*x(T1)+d*u(T1)
         % Note: Gives value of y[k+1]
         iy = find(abs(E2-t1)<tol);
         if ~isempty(iy)
            ciy = c(iy,:);
            Cj(iy,:) = ciy * Aj;
            Dj(iy,:) = ciy * Bj + d(iy,:) * Rj;
            % Note: Output delay modifies output offset!
            gj(iy,:) = gj(iy,:) + ciy * fj;
            E2(iy) = -1;
         end
      end

      % Set z[k]=u[k] for fid>0 so that w[k]=z[k-1]=u[k-1]
      Cj = [Cj;zeros(nw,nx)];
      aux = zeros(nw,2*nu+nw);  aux(:,jdelay) = eye(nw);
      Dj = [Dj;aux];
      w0 = u0(jdelay);
      y0 = [gj ; w0];
      u0 = [u0 ; w0];

      % Eliminate u[k+1] from the state equation by the change of 
      % state variable
      %    zeta[k] = x[k] - Bj3 * u[k]
      Bj3 = Bj(:,nu+nw+1:end);  Dj3 = Dj(:,nu+nw+1:end);
      ICMap = [eye(nx) -Bj3];
      Bj = Bj(:,1:nu+nw) + [Aj*Bj3-Bj3 , zeros(nx,nw)];  
      Dj = Dj(:,1:nu+nw) + [Cj*Bj3-Dj3 , zeros(ny+nw,nw)];
      x0 = x0-Bj3*u0(1:nu,:);  % zeta0
      dx0 = x0 + fj;

      % Handle remaining u[k+1] terms in output equation
      %  * When Dj3(p,:)=0, obtain yp[k] from yp[k+1] with unit output delay
      %  * When Dj3(p,:) nonzero, use internal delay
      aux = any(Dj3(1:ny,:),2);
      if any(aux) && CONSISTENT
         % Consistently use internal delay
         aux = true(ny,1);
      end
      ix = find(~(zfod | aux'));
      Delay.Output(ix) = Delay.Output(ix)+1;
      ix = find(aux);  nix = numel(ix);
      if nix>0
         % Note: nonzero rows of Dj3 are a subset of delayed outputs
         Bj = [Bj zeros(nx,nix)];
         Cj = [Cj ; Cj(ix,:)];   Cj(ix,:) = 0;
         % y(ix)[k] = gj(ix) + Dj3(ix,:) (u[k]-u0) + wix[k]
         Dj = [Dj ; Dj(ix,:)];
         Dj(ix,1:nu+nw+nix) = [Dj3(ix,:) zeros(nix,nw) eye(nix)];
         nw = nw+nix;
         u0 = [u0 ; zeros(nix,1)];
         y0 = [y0 ; zeros(nix,1)];
      end

      Dd = ltipack.ssdata(Aj,Bj,Cj,Dj,[],Ts);
      Delay.Internal = ones(nw,1);  % w[k]=z[k-1]
   end
   Dd.dx0 = dx0;
   Dd.x0 = x0;
   Dd.u0 = u0;
   Dd.y0 = y0;
   Dd.Delay = Delay;
end

% When original model had no internal delays, map new internal delays 
% to states if requested
if nfd==0 && nw>0 && USESTATE
   [Dd,xmap] = elimDelay(Dd,0,0,Dd.Delay.Internal);
   ICMap0 = ICMap;
   ICMap = zeros(numel(xmap),nx+nu);
   ICMap(xmap,:) = ICMap0;
end
