/* -*- MaTX -*-
 *
 *  NAME
 *      obsf() - Observable-unobservable decomposition
 *	
 *  SYNOPSIS
 *      {Ab,Bb,Cb,T,K} = obsf(A,B,C)
 *         Matrix Ab,Bb,Cb,T,K;
 *         Matrix A,B,C;
 *
 *      {Ab,Bb,Cb,T,K} = obsf(A,B,C,tol)
 *         Matrix Ab,Bb,Cb,T,K;
 *         Matrix A,B,C;
 *         Real tol;
 *
 *  DESCRIPTION
 *      obsf() returns a decomposition into the observable/unobservable
 *      subspaces.
 *
 *      If rank([[A][C]]) < Rows(A), then there is a similarity
 *      transformation  z = T x  such that
 *
 *         Ab = T*A*T# ,  Bb = T*B  ,  Cb = C*T#.
 *
 *      and the transformed system has the form
 *
 *              [Ano|A12]        [Bno]
 *        Ab =  [---+---],  Bb = [---],  Cb = [0|Co].
 *              [ 0 |Ao ]        [Bo ]
 *	                                             -1          -1
 *      where (Co,Ao) is observable, and Co(sI-Ao)Bo = C(sI-A)B.
 *
 *      obsf(..., tol) uses tolerance tol to determine the order of
 *      the subspaces.
 *
 *  SEE ALSO
 *      ctrf and obsm
 */

Func List obsf(A, B, C, tol, ...)
	Matrix A, B, C;
	Real tol;
{
	Matrix Ab, Bb, Cb, T, K;
	Matrix Ac, Bc, Cc;
	String msg;

	error(nargchk(3, 4, nargs, "obsf"));

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

	if (nargs == 3) {
		{Ac, Bc, Cc, T, K} = ctrf(A#, C#, B#);
	} else {
		{Ac, Bc, Cc, T, K} = ctrf(A#, C#, B#, tol);
	}

	Ab = Ac#;
	Bb = Cc#;
	Cb = Bc#;

	return {Ab, Bb, Cb, T, K};
}