/* -*- MaTX -*-
 *
 *  NAME
 *      lsim() - Simulation of continuous-time linear systems
 *
 *  SYNOPSIS
 *     	{Y,X} = lsim(A,B,C,D,U,T)
 *            Matrix Y,X;
 *            Matrix A,B,C,D;
 *            Matrix U,T;
 *
 *     	{Y,X} = lsim(A,B,C,D,U,T,x0)
 *            Matrix Y,X;
 *            Matrix A,B,C,D;
 *            Matrix U,T;
 *            Matrix x0;
 *
 *  DESCRIPTION
 *      lsim() calculates the time response of the continuous-time system:
 *         .
 *         x = Ax + Bu
 *         y = Cx + Du
 *
 *      to the input sequence U, where Rows(U) = Rows(u).  Each column of U 
 *      corresponds to a time point. 
 *
 *      lsim() returns the Rows(y)-by-Cols(U) output-history matrix Y
 *      and the stat-history matrix X.
 *
 *      lsim(..., x0) uses x0 to specify the initial conditions. 
 *
 *  SEE ALSO
 *      dlsim
 */

Func List lsim(A,B,C,D,U_,T,x0_, ...)
	Matrix A,B,C,D,x0_;
	Matrix U_,T;
{
	Integer m,n,N;
	Real dt;
	String msg;
	Matrix Aa,Ba,Ca,Da,Aad,Bad,X,Y,Xa,U,x0;
	
	error(nargchk(6, 7, nargs, "lsim"));

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

	{n,m} = size(B);
	N = Cols(U_);

	U = U_;

	if (nargs == 6) {
		x0 = Z(n,1);
	} else {
		x0 = x0_;
	}

	dt = T(2) - T(1);

	{Aa,Ba,Ca,Da} = series(Z(m),I(m),I(m),Z(m),A,B,C,D); // augmented system
	{Aad,Bad} = c2d(Aa,Ba,dt);

	// Initial states for the augmented system
	x0 = [[Z(m,1)]
		  [  x0  ]] + Ba*U(:,1);

	U(:,1:N-1) = (U(:,2:N) - U(:,1:N-1))/dt;
	U(:,N)     = U(:,N-1);

	Xa = ltitr(Aad,Bad,U,x0);
	Y = Ca*Xa + Da*U;
	
	X = Xa(m+1:m+n,:); // State of the original system

	return {Y,X};
}