function [Ad,Bd,Cd,Dd,Ed,dx0d,x0d,u0d,y0d,Delay] = c2dDF(sys,Ts,w,k,varargin)
% Tustin discretization fo LTV models.

%   Copyright 2022-2023 The MathWorks, Inc.

% Handle prewarp
if w==0
   T = Ts;
else
   % Handle prewarping
   T = 2*tan(w*Ts/2)/w;
end
tau = T/2;

% Get continuous-time data at t=k*Ts
[A,B,C,D,E,dx0,x0,u0,y0,Delay] = sys.DataFunction(k*Ts,varargin{:});
nfd = nfdelay(sys);
nx = size(A,1);
[nyz,nuw] = size(D);
nu = nuw-nfd;
ny = nyz-nfd;

% Prepare offsets
x0d = x0;
u0d = u0;
if isempty(dx0) && isempty(x0)
   dx0d = [];
   y0d = y0;
else
   if isempty(dx0)
      dx0 = zeros(nx,1);
   end
   if isempty(x0)
      x0 = zeros(nx,1);
   end
   if isempty(y0)
      y0 = zeros(ny,1);
   end
end

% Compute discrete matrices
% Note: Discretized state is w[k] = x(k*Ts) - (T/2) x'(k*Ts)
if isempty(E)
   % Explicit: Return explicit model
   M = eye(nx)-tau*A;
   [aux,Cd] = ltipack.util.safeMdivide(M, [eye(nx)+tau*A,T*B,T*dx0], C);
   if ~(allfinite(Cd) && allfinite(aux))
      error(message('Control:transformation:c2d10'))
   end
   Ad = aux(:,1:nx);        % (I - a*T/2)\(I + a*T/2)
   Bd = aux(:,nx+1:nx+nuw);  % (I - a*T/2)\b*T
   Dd = (C*Bd)/2 + D;       % (T/2)*c*(I - a*T/2)\b + d
   Ed = [];
   if ~isempty(dx0)
      dx0d = aux(:,nx+nuw+1);
      y0d = y0 + (C*dx0d)/2;
      dx0d = x0 + dx0d;
   end
else
   % Descriptor: Return descriptor model
   Ed = E-tau*A;
   Ad = E+tau*A;
   Bd = T*B;
   aux = matlab.internal.math.nowarn.mrdivide(C,Ed); % c/(e - a*T/2)
   if ~allfinite(aux)
      error(message('Control:transformation:c2d10'))
   end
   Cd = aux*E;
   Dd = D + (aux*Bd)/2;
   if ~isempty(dx0)
      dx0d = T*dx0 + Ed*x0;
      y0d = y0+(T/2)*(aux*dx0);
   end
end

if ~isempty(Delay)
   Delay.Input = round(Delay.Input/Ts);
   Delay.Output = round(Delay.Output/Ts);
   if isfield(Delay,'Internal')
      Delay.Internal = round(Delay.Internal/Ts);
   end
end
