function [y,x] = lsim(D,u,t,xinit,~)
% 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',xinit,[],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
nx = size(D.a,1);
if isempty(xinit)
   xinit = zeros(nx,1);
end

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

% Scaling
if Ts==0
   if isempty(D.e)
      [D,xinit] = scale(D,xinit,speye(size(D.a))+dt*abs(D.a));
   else
      [D,xinit] = scale(D,xinit,abs(D.e)+dt*abs(D.a)); % |E|+dt*|A|
   end
end

% Simulate
tau = D.Delay.Internal;
winit = zeros(size(tau));
y = DAEfun(D.a,D.b,D.c,D.d,D.e,D.dx0,D.x0,D.u0,D.y0,tau,u,t,xinit,winit);
y = y.';

% Account for output delays
if Ts==0
   for ct=1:ny
      tau = D.Delay.Output(ct)/dt;
      if tau>0
         % Output is zero up to the delay. Watch for ceil(1+o(eps))=2
         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
