%function y = dtrsp(sys,input,T,tfinal,x0)
%
%  Time response of a discrete system:
%	SYS:	packed form A,B,C,D system or constant.
%	INPUT:	varying format input vector or constant.
%	T:	sample time.
%	TFINAL:	final time value
%		(optional, default = input final time).
%	X0:	initial state
%		(optional, default = 0).
%
%   See also: TRSP, SAMHLD, TUSTIN, VINTERP, and VDCMATE.

function y = dtrsp(sys,input,T,tfinal,x0)
if nargin == 0,
    disp('usage: y = dtrsp(sys,input,T,tfinal,x0)')
    return
    end

if nargin < 3,
    error('insufficient number of input arguments')
    return
    end

[type,ny,nu,nx] = minfo(sys);
if type == 'syst',
    if nargin < 5,
        x0 = zeros(nx,1);
    else
        [m,n] = size(x0);
        if m ~= nx | n ~= 1,
            error('incorrect initial state vector dimension')
            return
       	    end
        end
    [a,b,c,d] = unpck(sys);
elseif type ~= 'cons',
    error('the system is not in packed or constant form')
    return
    end

%       check that the input vector is compatible.
 
[typeu,nru,ncu,nupts] = minfo(input);
if typeu == 'vary' & ncu == 1,
    [udat,uptr,ut] = vunpck(input);
elseif typeu == 'cons' & ncu == 1,
    nupts = 1;
    udat = input;
    ut = 0;
else
    error('invalid input time vector')
    return
    end

if nru ~= nu, 
    error('incompatible input vector dimension')
    return
    end

% 	an epsilon value for time interval calculation is created.

teps = T*1e-8;

%       The time vector is created.  The final time is derived from
%       this in case the user did not specify an integer number of
%       time steps.

if nargin < 4,
    yt = [ut(1):T:max(ut)].';
elseif tfinal < min(ut),
    error('final time less than initial time')
    return
else
    yt = [ut(1):T:tfinal].';
    end
tfinal = max(yt);
nypts = length(yt);

%       Check if the final time is less than the maximum time in the input.
%       If it is truncate the input.  Doing this now can save a
%       lot of time if the interpolation is required.

if max(ut) - tfinal > teps,
	nupts = max(find(ut <= tfinal+teps));
	ut = ut(1:nupts);
	udat = udat(1:uptr(nupts)+nru-1,:);
	end

%        If the input and output time vectors do not
%        match interpolate the input.

interflg = 0;
if nypts ~= nupts,
    interflg = 1;
elseif max(abs(yt - ut)) > teps,
    interflg = 1;
    end

if interflg == 1,
    uint = zeros(nru*nypts,1);
    disp('interpolating input vector (zero order hold)')
    if nypts <= nupts,
        for i = 1:nypts,
            idx = max(find(ut <= yt(i)+teps));
            uint((i-1)*nru+1:i*nru) = udat((idx-1)*nru+1:idx*nru);
            end
    elseif nupts == 1,
        uint = reshape(udat*ones(1,nypts),nru*nypts,1);
    else
        for i = 1:nupts-1,
            idx = find(yt >= ut(i) & yt < ut(i+1));
            uint((min(idx)-1)*nru+1:max(idx)*nru) = ...
              reshape(udat((i-1)*nru+1:i*nru)*ones(1,length(idx)), ...
              nru*length(idx),1);
            end
        idx = find(yt >= ut(nupts));
        uint((min(idx)-1)*nru+1:max(idx)*nru) = ...
          reshape(udat((i-1)*nru+1:i*nru)*ones(1,length(idx)), ...
          nru*length(idx),1);
        end
    udat = uint;
    end

%       reshape the input so that the matrix multiplications
%       line up.  At this point the ut should correspond to
%       yt.

udat = reshape(udat,nru,nypts);
	
%       the input vector has now been expanded to a matrix: u
%       constant systems are treated by a multiplication call

%  Initialize the x & y matrices to speed up the code

y = zeros(nypts*ny+1,2);
if type == 'syst',
    x = [x0 zeros(length(x0),nypts-1)];
    end

%       First, the single time point case is dealt with.  Here
%       only the d term of a system and x0 affect the output.

if nypts == 1,
    if type == 'syst',
        y(1:ny,1) = c*x + d*udat;
    else
        y(1:ny,1) = sys*udat;
        end

%       Now the multiple timepoint case is treated.  Here 
%       a check is made to determine if the system has dynamics.
%       The no dynamics case can be handled by a simple multiplication.

elseif type == 'syst',
    for i=2:nypts,
        x(:,i) = a*x(:,i-1) + b*udat(:,i-1);
        end
    y(1:ny*nypts,1) = reshape(c*x + d*udat,ny*nypts,1);
else
    y(1:ny*nypts,1) = reshape(sys*udat,ny*nypts,1);
    end

%       In the case of only a single
%       time point, a varying matrix is still created.  This is
%       after all, assumed to be a time response.

y(1:nypts,2) = yt;
y((ny*nypts)+1,2) = inf;
y((ny*nypts)+1,1) = nypts;

%-----------------------------------------------------------------------------


%
% Copyright MUSYN INC 1991,  All Rights Reserved
