function [A,B,C,D,E,dx0,x0,u0,y0,Delay] = d2cDF(sys,Ts,w,t,varargin)
% Tustin-based D2C fo LTV models.

%   Copyright 2022 The MathWorks, Inc.

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

% To obtain D2C model with continuous y response, interpolate discrete-time 
% data between time samples.
% Note: Assumes p constant on [k*Ts,(k+1)*Ts]
DF = sys.DataFunction;
f = t/Ts;
k = floor(f);
[Ad,Bd,Cd,Dd,Ed,dx0d,x0d,u0d,y0d,Delay] = DF(k,varargin{:});
nx = size(Ad,1);   
nfd = nfdelay(sys);
[nyz,nuw] = size(Dd);
nu = nuw-nfd;
ny = nyz-nfd;
alpha = f-k;
if alpha>sqrt(eps)
   [Ad2,Bd2,Cd2,Dd2,Ed2,dx0d2,x0d2,u0d2,y0d2,Delay2] = DF(k+1,varargin{:});
   Ad = (1-alpha)*Ad+alpha*Ad2;
   Bd = (1-alpha)*Bd+alpha*Bd2;
   Cd = (1-alpha)*Cd+alpha*Cd2;
   Dd = (1-alpha)*Dd+alpha*Dd2;
   Ed = (1-alpha)*Ed+alpha*Ed2;
   dx0d = (1-alpha)*dx0d+alpha*dx0d2;
   x0d = (1-alpha)*x0d+alpha*x0d2;
   u0d = (1-alpha)*u0d+alpha*u0d2;
   y0d = (1-alpha)*y0d+alpha*y0d2;
   if ~isempty(Delay)
      Delay.Input = (1-alpha)*Delay.Input+alpha*Delay2.Input;
      Delay.Output = (1-alpha)*Delay.Output+alpha*Delay2.Output;
      if isfield(Delay,'Internal')
         Delay.Internal = (1-alpha)*Delay.Internal+alpha*Delay2.Internal;
      end
   end
end

% Prepare offsets
x0 = x0d;
u0 = u0d;
if isempty(x0d) && isempty(dx0d)
   dx0 = [];
   y0 = y0d;
else
   if isempty(dx0d)
      dx0d = zeros(nx,1);
   end
   if isempty(x0d)
      x0d = zeros(nx,1);
   end
   if isempty(y0d)
      y0d = zeros(ny,1);
   end
   % At this stage dx0d and y0d are both nonempty
end

% Compute continuous matrices
if isempty(Ed)
   % Explicit: Return explicit model
   M = eye(nx)+Ad;
   [aux,Y] = ltipack.util.safeMdivide(M, (2*r)*[Ad-eye(nx),Bd,dx0d-x0d], Cd);
   if ~(allfinite(aux) && allfinite(Y))
      error(message('Control:transformation:d2c06'))
   end
   A = aux(:,1:nx);        % 2/T*(Ad+I)\(Ad-I)
   B = aux(:,nx+1:nx+nuw);  % 2/T*(Ad+I)\B
   C = 2*Y;                % 2*C/(Ad+I)
   D = Dd - Y*Bd;          % D-C/(Ad+I)*B
   E = [];
   if ~isempty(dx0d)
      dx0 = aux(:,nx+nuw+1);
      y0 = y0d-(Cd*dx0)/(2*r);
   end
else
   % Descriptor: Return descriptor model
   E = (Ad+Ed)/2;
   A = r*(Ad-Ed);
   B = r*Bd;
   Y = matlab.internal.math.nowarn.mrdivide(Cd/2,E);  % Cd/(Ad+Ed)
   if ~allfinite(Y)
      error(message('Control:transformation:d2c06'))
   end
   C = 2*(Y*Ed);
   D = Dd - Y*Bd;
   if ~isempty(dx0d)
      dx0d = dx0d-Ed*x0d;
      dx0 = r*dx0d;
      y0 = y0d-Y*dx0d;
   end
end

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

