function [Dd,ICMap] = c2d(Dc,Ts,options)
%C2D  Continuous-to-discrete conversion of MECHSS models.

%   Copyright 2023 The MathWorks, Inc.
if ~strcmp(options.Method,'tustin')
   error(message('Control:transformation:c2d28'))
end

% Handle prewarp
w = options.PrewarpFrequency;
if w == 0
    h = Ts/2;
else
    % Handle prewarping
    h = tan(w*Ts/2)/w;
end
h2 = h^2;

% Get dimensions
m = Dc.m;  k = Dc.k;
if isempty(m)
   m = speye(size(k));
elseif nnz(m)==0
   % Use first-order Tustin when M=0
   [Dd,ICMap] = utDiscretizeTustin(sparss(Dc),Ts,options);
   Dd = mechss(Dd);
   return
end
c = Dc.c;  d = Dc.d;  f = Dc.f;  g = Dc.g;
nfd  = numel(Dc.Delay.Internal);
ny = size(d,1) - nfd;
nu = size(d,2) - nfd;

% Get Fractional Delays options
nf = options.ThiranOrder;
roundOK = (nf == 0);
delayOK = strcmpi(options.DelayModeling,'delay');

% Extract discrete external delays
[Delay,fid,fod,ffd] = ltipack.AbstractSS.discretizeDelay(Dc.Delay,Ts);

% Compute discrete matrices. The discretized DoFs are
%   qd[k] = q(k*Ts)-h*q'(k*Ts) - h (q'(k*Ts)-h*q''(k*Ts))
aux = m+h2*k;
md = aux+h*c;
kd = aux-h*c;
cd = 2*(h2*k-m);
bd = 4*h2*Dc.b;
[~,Y] = ltipack.util.safeMdivide(md,[],g+h*f);
fd = (Y*c)/2-g/(2*h);
gd = f/2-(Y*cd)/(4*h);
dd = d + h * Y * Dc.b;
Dd = ltipack.mechdata(md,cd,kd,bd,fd,gd,dd,Ts);
Dd.StateInfo = Dc.StateInfo;

% Handle fractional delays
if any(ffd) && ~roundOK
    % Approximate using Thiran filters
    appxInternalDelays = min(Delay.Internal,nf-1) .* (ffd>0);
    G = ltipack.ssdata.getThiranAppx(appxInternalDelays+ffd,Ts,delayOK);
    G.Delay.Output = Delay.Internal - appxInternalDelays;
    Dd = lft(Dd,mechss(G),nu+(1:nfd),ny+(1:nfd),1:nfd,1:nfd); % nfd size may decrease
else
    % Round delays to nearest multiple of Ts
    Dd.Delay = ltipack.utDelayStruct(ny,nu,true);
    Dd.Delay.Internal = Delay.Internal + round(ffd);    
end
% Note: At this point Dd is well formed with correct dimensions for
% Dd.Delay for the MTIMES operations below.

if any(fid) && ~roundOK
    % Approximate using Thiran filters
    appxInputDelays = min(Delay.Input,nf-1) .* (fid>0);
    G = ltipack.ssdata.getThiranAppx(appxInputDelays+fid,Ts,delayOK);
    Dd = mtimes(Dd,mechss(G)); % Dd*G
    Dd.Delay.Input = Delay.Input-appxInputDelays;
else
    % Round delays to nearest multiple of Ts
    Dd.Delay.Input = Delay.Input + round(fid);
end

if any(fod) && ~roundOK
   % Approximate using Thiran filters
   appxOutputDelays = min(Delay.Output,nf-1) .* (fod>0);
   G = ltipack.ssdata.getThiranAppx(appxOutputDelays+fod,Ts,delayOK);
   % Multiply by the Output Delays. This results in the first states being those of G
   Dd = mtimes(mechss(G),Dd); % G*Dd
   Dd.Delay.Output = Delay.Output-appxOutputDelays;
else
   Dd.Delay.Output = Delay.Output + round(fod);
end

if roundOK && (any(fid) || any(fod) || any(ffd))
   warning(message('Control:transformation:RoundedDelay'))
end

% Eliminate zero internal delays (may be introduced by rounding)
Dd = elimZeroDelay(Dd);

% IC Mapping q(0),q'(0),u(0) -> nu[0],nu[1] involves large dense matrices.
ICMap = [];
