function Dd = c2dXC(Dc,Ts,options)
% Discretization with state and delay consistency across array.
% Given an SS array with constant number of states and internal delays,
% C2DXC returns a discretized SS array with the same property.

%   Copyright 2023 The MathWorks, Inc.
nsys = numel(Dc);
Dd = ltipack.ssdata.array(size(Dc));
METHOD = options.Method(1);

if METHOD=='t'
   % Tustin
   for ct=1:nsys
      Dd(ct) = utDiscretizeTustin(Dc(ct),Ts,options);
   end
else
   % ZOH, FOH, IMP
   % Consistently eliminate E
   for ct=1:nsys
      Dc(ct) = des2exp(Dc(ct),true);
      if ~isfinite(Dc(ct))
         error(message('Control:transformation:dss2ss1'))
      end
   end

   MINFD = 1e-3;  % min fractional delay
   if any(METHOD=='zf')
      % Fractional input delays gives rise to extra internal delays. Make
      % sure each input delay is consistently integral or fractional.
      FID = false;
      for ct=1:nsys
         [~,fid] = ltipack.util.splitDelay(Dc(ct).Delay.Input,Ts);
         FID = FID | (fid>0);
      end
      for ct=1:nsys
         [N,fid] = ltipack.util.splitDelay(Dc(ct).Delay.Input,Ts);
         fid(FID) = max(fid(FID),MINFD);
         Dc(ct).Delay.Input = (N+fid)*Ts;
      end
   end

   if METHOD=='f'
      % (1) Fractional output delays can give rise to extra internal delays.
      % Same for fractional internal delays which get mapped to output 
      % delays of H(s). Make sure each delay is consistently integral or 
      % fractional. 
      % (2) Internal delay tau<Ts gives rise to a zero discrete delay 
      % that must be upped to one for simulation, resulting in a total
      % delay of tau+Ts. Enforcing tau>=Ts is therefore more accurate
      FOD = false;  FXD = false;
      for ct=1:nsys
         [~,fod] = ltipack.util.splitDelay(Dc(ct).Delay.Output,Ts);
         [~,fxd] = ltipack.util.splitDelay(max(Dc(ct).Delay.Internal,Ts),Ts);
         FOD = FOD | (fod>0);
         FXD = FXD | (fxd>0);
      end
      for ct=1:nsys
         [Nod,fod] = ltipack.util.splitDelay(Dc(ct).Delay.Output,Ts);
         [Nxd,fxd] = ltipack.util.splitDelay(max(Dc(ct).Delay.Internal,Ts),Ts);
         fod(FOD) = max(fod(FOD),MINFD);
         fxd(FXD) = max(fxd(FXD),MINFD);
         Dc(ct).Delay.Output = (Nod+fod)*Ts;
         Dc(ct).Delay.Internal = (Nxd+fxd)*Ts;
      end
   end

   % Discretize
   switch METHOD
      case 'z'
         for ct=1:nsys
            Dd(ct) = utDiscretizeZOH(Dc(ct),Ts,options);
         end
      case 'f'
         for ct=1:nsys
            Dd(ct) = utDiscretizeFOH(Dc(ct),Ts,options);
         end
      case 'i'
         for ct=1:nsys
            Dd(ct) = utDiscretizeIMP(Dc(ct),Ts,options);
         end
   end

end


