/* -*- MaTX -*-
 *
 *  NAME
 *      impulse() - Impulse response of continuous-time linear system
 *
 *  SYNOPSIS
 *      {Y,X} = impulse(A,B,C,D,iu,T)
 *          Matrix X,Y;
 *          Matrix A,B,C,D;
 *          Integer iu;
 *          Matrix T;
 *
 *  DESCRIPTION
 *       impulse(A,B,C,D,iu,T) calculates the response of the system:
 *      	.
 *      	x = Ax + Bu
 *      	y = Cx + Du
 *
 *       to an unit impulse applied to the iu'th input.  T is a regularly 
 *       spaced time vector which specifies the time axis.
 * 
 *       impulse() returns the Rows(y)-by-length(T) output-history matrix Y 
 *       and the state-history matrix X.
 *
 *       If iu == 0, impulse() 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
 *      step and dimpulse
 */

Func List impulse(A,B,C,D,iu,T)
	Matrix A,B,C,D;
	Integer iu;
	Matrix T;
{
	Integer i,m,N;
	Real dt;
	Matrix Ad,Bd,b,x0,X,Y,XX,YY,U;
	String msg;

	if (length((msg = abcdchk(A, B, C, D))) > 0) {
		error("impulse(): " + msg);
	}

	N = length(T);
	U = Z(1,N);
	dt = T(2) - T(1);

	if (iu == 0) {
		m = Cols(B);
		YY = Z(Rows(C)*m,N);
		XX = Z(Rows(A)*m,N);

		for (i = 1; i <= m; i++) {
			x0 = b = B(:,i);
			{Ad, Bd} = c2d(A, b, dt);
			X = ltitr(Ad, Bd, U, x0);
			Y = C*X;
			
			if (version() == 4) {
				YY(i-1,0,Y) = Y;
				XX(i-1,0,X) = X;
			} else {
				YY(i,1,Y) = Y;
				XX(i,1,X) = X;
			}
		}
	} else {
		x0 = b = B(:,iu);
		{Ad, Bd} = c2d(A, b, dt);
		XX = ltitr(Ad, Bd, U, x0);
		YY = C*XX;
	}

	return {YY,XX};
}