/* -*- MaTX -*-
 *
 *  NAME
 *      step_ss()  - Step response of continuous-time linear systems
 *      step_tf()  - Step response of continuous-time linear systems
 *      step_tfn() - Step response of continuous-time linear systems
 *      step_tfm() - Step response of continuous-time linear systems
 *
 *  SYNOPSIS
 *      {Y,X} = step_ss(A,B,C,D,iu,T)
 *         Matrix X,Y;
 *         Matrix A,B,C,D;
 *         Integer iu;
 *         Matrix T;
 *
 *      {Y,X} = step_tf(num,den,T)
 *         Matrix X,Y;
 *         Matrix num,den;
 *         Matrix T;
 *
 *      {Y,X} = step_tfn(g,T)
 *         Matrix X,Y;
 *         Rational g;
 *         Matrix T;
 *
 *      {Y,X} = step_tfm(G,iu,T)
 *         Matrix X,Y;
 *         RaMatrix G;
 *         Integer iu;
 *         Matrix T;
 *
 *  DESCRIPTION
 *       step_ss(A,B,C,D,iu,T) returns the response of the system:
 *      	.
 *      	x = Ax + Bu
 *      	y = Cx + Du
 *
 *       to an unit step applied to the iu-th input.  T is a regularly 
 *       spaced time vector which specifies the time axis.
 *
 *       step_ss() returns the Rows(y)-by-length(T) output-history matrix Y
 *       and the state-history matrix X.
 *
 *       If iu == 0, step_ss() 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
 *      impluse and dstep
 */

Func List step_ss(A,B,C,D,iu,T)
	Matrix A,B,C,D;
	Integer iu;
	Matrix T;
{
	Integer i,m,N;
	Matrix X,Y,XX,YY;
	Matrix A2,B2,C2,D2;
	String msg;

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

	if (iu == 0) {
		m = Cols(B);
		N = length(T);
		YY = Z(Rows(C)*m,N);
		XX = Z(Rows(A)*m,N);
		for (i = 1; i <= m; i++) {
			{A2,B2,C2,D2} = series([0],[1],[1],[0],A,B(:,i),C,D(:,i));
			{Y,X} = impulse(A2,B2,C2,D2,1,T);
			X = X(2:,:); // remove the state of integrator

			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 {
		{A2,B2,C2,D2} = series([0],[1],[1],[0],A,B(:,iu),C,D(:,iu));
		{YY,XX} = impulse(A2,B2,C2,D2,1,T);
		XX = XX(2:,:); // remove the state of integrator
	}

	return {YY,XX};
}

Func List step_tfn(g, T)
	Rational g;
	Matrix T;
{
	Matrix A,B,C,D;
	Matrix X, Y;
	
	{A, B, C, D} = tfn2ss(g);
	{Y, X} = step_ss(A, B, C, D, 1, T); 
	return {Y,X};
}

Func List step_tfm(G, iu, T)
	RaMatrix G;
	Integer iu;
	Matrix T;
{
	Matrix A,B,C,D;
	Matrix X, Y;
	
	if (iu == 0) {
		{A, B, C, D} = tfm2ss(G);
		{Y, X} = step_ss(A, B, C, D, iu, T); 
	} else {
		{A, B, C, D} = tfm2ss(G, iu);
		{Y, X} = step_ss(A, B, C, D, 1, T); 
	}
	return {Y, X};
}

Func List step_tf(num, den, T)
	Matrix num, den;
	Matrix T;
{
	Matrix A,B,C,D;
	Matrix X, Y;
	
	{A, B, C, D} = tf2ss(num,den);
	{Y, X} = step_ss(A, B, C, D, 1, T); 
	return {Y, X};
}

Func List step(A,B,C,D,iu,T)
	Matrix A,B,C,D;
	Integer iu;
	Matrix T;
{
	Matrix XX,YY;

	{YY,XX} = step_ss(A,B,C,D,iu,T);
	return {YY,XX};
}

Func List Step_tf(g, T)
	Rational g;
	Matrix T;
{
	Matrix XX,YY;

	pause "Step_tf() is bosolete.  Use step_tfn()", 1;
	print "\n";

	{YY,XX} = step_tfn(g,1,T);
	return {YY,XX};
}

Func List Step_tfm(G, iu, T)
	RaMatrix G;
	Integer iu;
	Matrix T;
{
	Matrix XX, YY;

	pause "Step_tfm() is bosolete.  Use step_tfm()", 1;
	print "\n";

	{YY,XX} = step_tfm(G,iu,T);
	return {YY,XX};
}