function [y,Nfocus,x] = tresp(Dsim,RespType,Ns,SimInfo)
% Step or impulse response of discrete MIMO transfer functions.
% 
% Inputs:
%   * Dsim: discrete MIMO @tfdata model
%   * Resptype: 'step' or 'impulse'
%   * Ns: number of time steps (set to [] for automatic selection)
%   * SimInfo: structure containing DC info for settling time detection 
%     (used only for Ns=[])
%
% Outputs:
%   * y: output history (ns-by-ny-by-nu)
%   * nFocus: Number of time steps to be shown by default in plots
%     (-1 if divergent, 0 if the response has neither settled nor 
%      diverged)

%	 Author: P. Gahinet
%   Copyright 1986-2022 The MathWorks, Inc.
[ny,nu] = iosize(Dsim);
x = [];
Ts = abs(Dsim.Ts);

% Note: TF simulation always starts from steady state (SimInfo.xinit=[]).
%       Under this assumption, offsets can be ignored for simulation and
%       the problem decouples into NU independent simulations of each 
%       input channel.
if isempty(Ns)
   % Infinite-horizon simulation (terminates when response settles
   % or diverges)
   [y,Nfocus] = tfresp(Dsim,RespType,SimInfo);
   y = permute(y,[3 1 2]);
   
else
   % Finite-horizon simulation
   Nfocus = Ns;
   % Limit delays to simulation horizon to avoid "out of memory" errors in TFSIM
   iod = min(getIODelay(Dsim,'total'),Ns);
   ut = zeros(Ns,1);
   y = zeros(Ns,ny,nu);
   for j=1:nu
      nj = Dsim.num(:,j);
      dj = Dsim.den(:,j);
      iodj = iod(:,j);
      % Simulate
      kChg = 1+round(SimInfo.Delay/Ts);
      switch RespType
         case 'step'
            ut(kChg:Ns,:) = SimInfo.du(j);
         case 'impulse'
            ut(kChg) = SimInfo.du(j)/Ts;
      end
      InitState = ltipack.tfsimIC(nj,dj,iodj);
      y(:,:,j) = tfsim(nj,dj,iodj,ut,InitState);
   end
end

% Apply output offsets corresponding to input offsets
y = y + reshape(SimInfo.y0,[1 ny]);
