function [Dsim,dt,tf,SimInfo] = trespSetUp(D,RespType,dt,tf,SimInfo)
% Builds discrete-time model for STEP, IMPULSE, or INITIAL simulation.
%
% Inputs:
%   * D: continuous or discrete, SISO or MIMO state-space model with no 
%        E matrix.
%   * dt: user-defined sampling time ([] if none)
%   * tf: final time ([] if unspecified)
%   * SimInfo: Simulation parameters (preliminary)
%
% Outputs:
%   * Dsim: discrete-time model
%   * dt: sampling time of Dsim
%   * tf: adjusted final time
%   * SimInfo: Simulation parameters (finalized)

%	 Author: P. Gahinet
%   Copyright 1986-2023 The MathWorks, Inc.
isCT = (D.Ts==0);
NoTs = (isCT && isempty(dt));
NoTf = isempty(tf);

% When D is continuous, determine adequate step size for simulation
% Note: No attempt to assess stability, response may still converge when
% unstable, e.g., impulse(tf(1,[1 2 0])) or g339096
if NoTs
   % Compute largest delay
   MaxDelay = getMaxDelay(D);
   % Determine adequate step size dt
   [a,b,c] = getABCD(D);
   Nsmax = SimInfo.MaxSample;
   if NoTf
      Nsmax = Nsmax/2;  % leave room for settling detection
   end
   switch RespType
      case 'step'
         dt = ltipack.ssdata.timegrid(a,b,c,[],MaxDelay,Nsmax,tf);
      case 'impulse'
         dt = ltipack.ssdata.timegrid(a,[],c,b,MaxDelay,Nsmax,tf);
      case 'initial'
         dt = ltipack.ssdata.timegrid(a,[],c,SimInfo.xinit,MaxDelay,Nsmax,tf);
   end
   % Make sure to hit step or impulse time
   if SimInfo.Delay>0
      N = ceil(SimInfo.Delay/dt);
      dt = SimInfo.Delay/N;
   end
end  

% NOTE:
%   1) Past output value yinit are needed when output delays are present
%   2) SSRESP and SSIM simulate
%         x[k+1] =  A x[k] +  B1 u[k] +  B2 w[k] + f
%           y[k] = C1 x[k] + D11 u[k] + D12 w[k] + g1
%           z[k] = C2 x[k] + D21 u[k] + D22 w[k] + g2
%           w[k] = z[k-Delays]
%      with initial state xinit and zero u[k],w[k],z[k] prior to start.
%      The affine terms f,g account for input and internal signal offsets
%      in STEP and IMPULSE. They are needed for correct IMP discretization
%      when u(t) = uinit + du * delta(t).
[ny,nu] = iosize(D);
INITIAL = strcmp(RespType,'initial');
if INITIAL
   % Past y,z values are all zero (truly zero, unknown, or unsteady)
   SimInfo.yinit = zeros(ny,1);
else
   % Absorb input offset UINIT into F,G so that baseline input is U=0.
   % This is expected by the IMP discretization DSIM.
   uinit = SimInfo.uinit;
   D.dx0 = D.dx0 + D.b(:,1:nu) * uinit;
   D.y0 = D.y0 + D.d(:,1:nu) * uinit;
end

% Discretize
if ~isCT
   Dsim = D;
elseif strcmp(RespType,'impulse')
   Dsim = utDiscretizeIMP(D,dt);
else
   Dsim = utDiscretizeZOH(D,dt);
end

% Resolve initial condition and absorb WINIT into F,G
% Note: Must be done for DSIM to get correct WINIT
if ~INITIAL
   [SimInfo.xinit,winit,~,SimInfo.yinit] = resolveIC(Dsim,SimInfo.xinit,zeros(nu,1));
   Dsim.u0 = Dsim.u0-[zeros(nu,1);winit];
   Dsim.y0 = Dsim.y0-[zeros(ny,1);winit];
end
