/*  -*- MaTX -*-
 *
 *  NAME
 *      tzero() - Transmission zeros
 *
 *  SYNOPSIS
 *      z = tzero(A,B,C,D)
 *         CoMatrix z;
 *         Matrix A,B,C,D;
 *
 *  DESCRIPTION
 *      tzero() returns the transmission zeros of the system
 *   		.
 *   		x = Ax + Bu
 *   		y = Cx + Du
 *
 *  SEE ALSO
 *      poles and zeros
 */

Func CoMatrix tzero(A,B,C,D)
	Matrix A,B,C,D;
{
	Integer i,m,n,p;
	Real mu;
	String msg;
	Matrix AA,A1,A2,BB;
	CoMatrix tz,g1,g2;

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

	n = Rows(A);
	m = Cols(B);
	p = Rows(C);

	AA = [[A B]
		  [C D]];

	if (m == p) {
		BB = diag(I(n), Z(p,m));
		tz = eigval(AA,BB);
		tz = tz( ! isnan(tz) && isfinite(tz));
	} else {
		mu = norm(AA,1);
		if (m > p) {
			A1 = [[           AA            ]
				  [mu*(rand(m-p,n+m) .- 0.5)]];
			A2 = [[           AA            ]
				  [mu*(rand(m-p,n+m) .- 0.5)]];
		} else {
			A1 = [AA mu*(rand(n+p,p-m) .- 0.5)];
			A2 = [AA mu*(rand(n+p,p-m) .- 0.5)];
		}
		BB = diag(I(n), Z(Rows(A1)-n,m));
		g1 = eigval(A1,BB);
		g2 = eigval(A2,BB);
		g1 = g1( ! isnan(g1) && isfinite(g1));
		g2 = g2( ! isnan(g2) && isfinite(g2));

		tz = [];
		for (i = 1; i <= length(g1); i++) {
			if (any(abs(Array(g1(i) .- g2)) .< mu*sqrt(EPS))) {
				tz = [[tz][g1(i)]];
			}
		}
	}
	return tz;
}