function [y,x] = lsim(D,u,t,x0,InterpRule)
% 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 1986-2012 The MathWorks, Inc.
Ns = size(u,1);
dt = t(2)-t(1);
Ts = D.Ts;

% Computability and consistency checks
ComputeX = (nargout>1);
D = checkCompatibility(D,'lsim',x0,[],ComputeX);

% Select interpolation rule if unspecified
if Ts==0 && strcmp(InterpRule,'auto')
   InterpRule = utSelectInterp(u);
end

% Check for complex data
[CplxFlag,D,u,x0] = localComplex2Real(D,u,x0);
[ny,nu] = iosize(D);

% Branch based on the integration technique
if Ts==0 && hasDelayDynamics(D,'pole')
   % Simulate with DDAE solver when delay dynamics
   if ComputeX
      [y,x] = ddaesim(D,u,t,InterpRule);
   else
      y = ddaesim(D,u,t,InterpRule);
   end
else
   % Explicit integration: D is discrete or can be discretized
   nx = size(D.a,1);
   if isempty(x0)
      x0 = zeros(nx,1);
   end
   u = u.'; % now nu-by-Ns

   % Discretize continuous models
   if Ts==0
      % Discretize
      switch InterpRule
         case 'zoh'
            [D,ICMap] = utDiscretizeZOH(D,dt);
         case 'foh'
            % FOH discretization produces the model
            %     z[k+1] = exp(A*dt) * z[k] + Bd * u[k]
            %       y[k] =     Cd    * z[k] + Dd * u[k]
            % where z[k] = x[k] + G * u[k].  For simulation, the initial condition
            % must be set to z[0]=x0+G*u[0], and the state trajectory is obtained
            % as z[k] - G * u[k].
            [D,ICMap] = utDiscretizeFOH(D,dt);
      end
      % Update initial condition
      u0 = u(:,1);
      u0(D.Delay.Input>0,:) = 0;
      x0 = ICMap * [x0;u0];
   end

   % Simulate with SSSIM
   % RE: * D is discrete at this point
   %     * Limit discrete delays to the number of samples (g172142)
   [f,g] = getAffineOffsets(D);
   x0 = linsimstate('ss',x0,min(D.Delay.Input,Ns),...
      min(D.Delay.Output,Ns),min(D.Delay.Internal,Ns));
   if ComputeX
      [y,~,z] = sssim(D.a,D.b,D.c,D.d,[],f,g,u,x0);
      % x[k] = z[k] - G * u[k]
      x = z(1:nx,:);
      if Ts==0
         for j=1:nu
            Din = min(D.Delay.Input(j),Ns);
            u(j,:) = [zeros(1,Din) , u(j,1:Ns-Din)];
         end
         x = x - ICMap(1:nx,nx+1:nx+nu) * u;
      end
   else
      y = sssim(D.a,D.b,D.c,D.d,[],f,g,u,x0);
   end

   % Format output arrays
   y = y.';
   if ComputeX
      x = x.';
   end
end

% Return complex data for complex inputs
if CplxFlag
   ny = ny/2;
   y = complex(y(:,1:ny),y(:,ny+1:2*ny));
   if ComputeX
      nx = size(x,2)/2;
      x = complex(x(:,1:nx),x(:,nx+1:2*nx));
   end
end

end

%------------ Local Functions ------------------------------

function [CplxFlag,D,u,x0] = localComplex2Real(D,u,x0)
% Convert complex data into real data with twice as many states
% and outputs.
CplxFlag = ~(isreal(u) && isreal(x0) && isreal(D));
if CplxFlag
   D = cplx2real(D);
   if ~isempty(x0)
      x0 = [real(x0);imag(x0)];
   end
   u = [real(u) , imag(u)];
end
end
