function [y,x] = lsim(D,u,t,x0,~)
% Linear response simulation of state-space model with real coefficients.
% RE: Assume U is Ns x Nu and T has Ns samples.

%	 Author: P. Gahinet
%   Copyright 2020 The MathWorks, Inc.
[ny,nu] = iosize(D);
x = [];   % never compute X for sparse
Ts = D.Ts;
Ns = numel(t);
dt = t(2)-t(1);
OPTS = getSolverConfig(D);
if Ts==0
   DAEfun = str2func(OPTS.DAESolver);
else
   DAEfun = @DSIM;
end

% Computability and consistency checks
D = checkCompatibility(D,'lsim',x0,[],false);

% Apply input delays, interpolating for fractional delays
if Ts==0
   for ct=1:nu
      tau = D.Delay.Input(ct)/dt;
      if tau>0
         n = floor(tau);   f = tau-n;
         u(:,ct) = [zeros(n,1) ; (1-f)*u(1,ct) ; f*u(1:Ns-n-1,ct)+(1-f)*u(2:Ns-n,ct)];
      end
   end
else
   for ct=1:nu
      n = D.Delay.Input(ct);
      if n>0
         % Shift by n samples
         u(:,ct) = [zeros(n,1) ; u(1:Ns-n,ct)];
      end
   end
end
u = u.'; % now nu-by-Ns

% Initial condition
nq = size(D.k,1);
if isempty(x0)
   x0 = zeros(2*nq,1);
elseif numel(x0)==nq
   x0 = [x0 ; zeros(nq,1)];  % assume q'=0
end

% Internal delays
[D,x0] = flushZeroDelay(D,x0);
   
% Seek arrow structure
[D,perm] = xsort(D);
if ~isempty(perm)
   x0 = x0([perm;nq+perm],:);
end

% Scaling
if Ts==0
   if isempty(D.m)
      [D,x0] = scale(D,x0,speye(size(D.k))+dt*abs(D.c)+dt^2*abs(D.k));
   else
      [D,x0] = scale(D,x0,abs(D.m)+dt*abs(D.c)+dt^2*abs(D.k));
   end
end

% Simulate
tau = D.Delay.Internal;
w0 = zeros(size(tau));
y = DAEfun(D.m,D.c,D.k,D.b,D.f,D.g,D.d,tau,u,t,x0,w0);
y = y.';

% Account for output delays
if Ts==0
   for ct=1:ny
      tau = D.Delay.Output(ct)/dt;
      if tau>0
         % Note: Output is zero up to the delay.
         n = ceil(0.999999*tau);   f = max(0,n-tau);
         y(:,ct) = [zeros(n,1) ; (1-f)*y(1:Ns-n,ct)+f*y(2:Ns-n+1,ct)];
      end
   end
else
   for ct=1:ny
      n = D.Delay.Output(ct);
      if n>0
         % Shift by n samples
         y(:,ct) = [zeros(n,1) ; y(1:Ns-n,ct)];
      end
   end
end
