/* -*- MaTX -*-
 *
 *  NAME
 *      dimpulse() - Impulse response of discrete-time linear systems
 *
 *  SYNOPSIS
 *     {Y,X} = dimpulse(A,B,C,D,iu,N)
 *        Matrix Y, X;
 *        Matrix A,B,C,D;
 *        Integer iu,N;
 *
 *  DESCRIPTION
 *       dimpulse(A,B,C,D,iu,N) calculates the response of the system:
 *
 *         x[n+1] = Ax[n] + Bu[n]
 *         y[n]   = Cx[n] + Du[n]
 *
 *       to an unit impulse applied to the iu'th input.
 *
 *       dimpulse() returns the Rows(y)-by-N output-history matrix Y, and
 *       the state-history matrix X.
 *
 *       If iu = 0, dimpuse() returns matrices
 *
 *       [[Y for 1st input]          [[X for 1st input]
 *        [Y for 2nd input]           [X for 2nd input]
 *        [...............]           [...............]
 *        [Y for m'th input]]  and    [Y for m'th input]].
 *
 *  SEE ALSO
 *      impulse and dstep
 */

Func List dimpulse(A,B,C,D,iu,N)
	Matrix A,B,C,D;
	Integer iu,N;
{
	Integer i,m;
	Matrix X,Y,XX,YY;
	String msg;

	msg = abcdchk(A, B, C, D);
	if (length(msg) > 0) {
		error("dimpulse(): " + msg);
	}
	
	if (iu == 0) {
		m = Cols(B);
		YY = Z(Rows(C)*m,N);
		XX = Z(Rows(A)*m,N);
		if (version() == 4) {
			for (i = 1; i <= m; i++) {
				{Y,X} = dlsim(A, B(:,i), C, D(:,i), [ONE(1), Z(1,N-1)]);
				YY(i-1,0,Y) = Y;
				XX(i-1,0,X) = X;
			}
		} else {
			for (i = 1; i <= m; i++) {
				{Y,X} = dlsim(A, B(:,i), C, D(:,i), [ONE(1), Z(1,N-1)]);
				YY(i,1,Y) = Y;
				XX(i,1,X) = X;
			}
		}
	} else {
		{YY,XX} = dlsim(A, B(:,iu), C, D(:,iu), [ONE(1), Z(1,N-1)]);
	}

	return {YY,XX};
}