function [y,t,focus,x,ysd] = timeresp(D,RespType,t,Config)
%TIMERESP  Time response of single SPARSS model with real coefficients.
%
%   [Y,T,FOCUS,X] = TIMERESP(D,RESPTYPE,T,X0) computes step, impulse, or initial
%   response based on the value of RESPTYPE ('step', 'impulse', or 'initial').
%   The input T is either 
%      * [] for infinite-horizon simulations
%      * a scalar specifying the final time Tf 
%      * a 1x2 vector T=[T0 Tf] specifying the start time and final time
%      * a vector of time samples.
%   TIMERESP assumes that ISREAL(D) is true.
%   
%   When T=[], FOCUS is set the preferred simulation horizon. Otherwise FOCUS is 
%   equal to T(END) unless the final time has to be reduced to prevent overflow 
%   (for unstable systems only).

%	 Author: P. Gahinet
%   Copyright 2022 The MathWorks, Inc.
INITIAL = strcmp(RespType,'initial');
[ny,nu] = iosize(D);
x = [];   % never compute X for sparse
ysd = []; % default handling for model containing no covariance info
Ts = D.Ts;
OPTS = getSolverConfig(D);
if Ts==0
   DAEfun = str2func(OPTS.DAESolver);
else
   DAEfun = @DSIM;
end

% Compatibility checks
if INITIAL
   D = getsubsys(D,':',[]); % eliminate input channels
   nu = 1;
else
   Config = validate(Config,nu);
   if ischar(Config.InputOffset)
      Config.InputOffset = D.u0;
   end
end
[D,xinit] = checkCompatibility(D,RespType,Config.InitialState);

% Interpret T (start, stop, step size)
[t,t0,tf,dt] = ltipack.util.getTimeInfo(t,Ts);
if isempty(t)
   if Ts==0
      t = linspace(t0,tf,1001).';
      dt = (tf-t0)/1000;
   else
      t = (t0:dt:tf).';
   end
end
Ns = numel(t);
focus = [t0,tf];

% Quick exit
if nu==0
   % Quick exit for systems w/o inputs (step/impulse only)
   y = zeros(Ns,ny,0);
   return
end

% Resolve initial condition including internal delay signals
if INITIAL
   % Past y,z values are all zero (truly zero, unknown, or unsteady)
   winit = zeros(numel(D.Delay.Internal),1);
   ypast = zeros(ny,1);
else
   [xinit,winit,~,ypast] = resolveIC(D,xinit,Config.InputOffset);
end

% Internal delays
[D,xinit,winit] = flushZeroDelay(D,xinit,winit);

% 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

% Configure parallelism
UseParallel = (nu>1 && ltipack.util.checkParallel(OPTS.UseParallel));
if UseParallel
   ParPool = gcp;
else
   ParPool = matlab.internal.threadPool;
end

% Simulate
a = D.a;  b = D.b;  c = D.c;  d = D.d;  e = D.e;  
dx0 = D.dx0;  x0 = D.x0;  u0 = D.u0;  y0 = D.y0;
tau = D.Delay.Internal;
switch RespType
   case 'step'
      % Step response
      u = repmat(Config.InputOffset,[1 Ns]);
      du = Config.Amplitude;
      % Total step delay in sampling periods
      % Note: Can't account for fractional delays (step occuring between
      % time steps) so round input delays
      if Ts==0
         nDelay = round((D.Delay.Input+Config.Delay)/dt);
      else
         nDelay = D.Delay.Input+round(Config.Delay/dt);
      end
      y = zeros(ny,Ns,nu);
      parfor (j=1:nu, ParPool)
         uj = u;
         uj(j,nDelay(j)+1:Ns) = uj(j,nDelay(j)+1:Ns) + du(j);
         y(:,:,j) = feval(DAEfun,a,b,c,d,e,dx0,x0,u0,y0,tau,uj,t,xinit,winit);
      end
      
   case 'impulse'
      % Impulse response
      % Note: In continuous time, state jumps (see PULSE below) are the 
      % sure way to accurately capture the response to a Dirac impulse.
      % Setting u=1/dt at impulse time works only for Delay>0 (half of 
      % spike area is missing when Delay=0) and can cause a huge spike at 
      % impulse time. 
      u = repmat(Config.InputOffset,[1 Ns]);
      du = Config.Amplitude;
      % Total impulse delay in sampling periods
      % Note: Need to round input delays to keep impulses in sync with
      % integration steps.
      if Ts==0
         nDelay = round((D.Delay.Input+Config.Delay)/dt);
      else
         nDelay = D.Delay.Input+round(Config.Delay/dt);
      end
      y = zeros(ny,Ns,nu);
      parfor (j=1:nu, ParPool)
         uj = u;
         if Ts==0
            uP = zeros(nu,1);   uP(j) = du(j);
            PULSE = struct('k',nDelay(j)+1,'uk',uP);
            y(:,:,j) = feval(DAEfun,a,b,c,d,e,dx0,x0,u0,y0,tau,uj,t,xinit,winit,PULSE);
         else
            uj(j,nDelay(j)+1) = uj(j,nDelay(j)+1) + du(j)/dt;
            y(:,:,j) = feval(DAEfun,a,b,c,d,e,dx0,x0,u0,y0,tau,uj,t,xinit,winit);
         end
      end
      
   case 'initial'
      % Initial response
      u = zeros(0,Ns);
      y = DAEfun(D.a,D.b,D.c,D.d,D.e,dx0,x0,u0,y0,tau,u,t,xinit,winit);
end

% Account for output delays
y = permute(y,[2 1 3]);
nsim = size(y,3);
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);
         aux = repmat(ypast(ct),[n 1]);
         for j=1:nsim
            y(:,ct,j) = [aux ; (1-f)*y(1:Ns-n,ct,j)+f*y(2:Ns-n+1,ct,j)];
         end
      end
   end
else
   for ct=1:ny
      n = D.Delay.Output(ct);
      if n>0
         % Shift by n samples
         aux = repmat(ypast(ct),[n 1]);
         for j=1:nsim
            y(:,ct,j) = [aux ; y(1:Ns-n,ct,j)];
         end
      end
   end
end
