/*  -*- MaTX -*-
 *
 *  NAME
 *      hankel() - Hankel matrix
 *
 *  SYNOPSIS
 *      H = hankel(x)
 *         Matrix H;
 *         Matrix x;
 *
 *      H = hankel(x, y)
 *         Matrix H;
 *         Matrix x;
 *         Matrix y;
 *
 *  DESCRIPTION
 *      hankel(x) is a square Hankel matrix whose first column is x and
 *      whose elements are zero below the first anti-diagonal.
 *
 *      hankel(x,y) is a Hankel matrix whose first column is x and whose
 *      last row is y.
 *
 *  SEE ALSO
 *     toeplitz and vander
 */


Func Matrix hankel(x_, y_, ...)
	Matrix x_, y_;
{
	Integer j, nx, ny;
	Matrix x,y,H;

	error(nargchk(1, 2, nargs, "hankel"));
	
	if (nargs == 1) {
		x = makecolv(x_);
	} else {
		x = makecolv(x_);
		y = makecolv(y_);
	}

	nx = length(x);
	if (nargs == 1) {
		H = Z(nx);
		for (j = 1; j <= nx; j++) {
			H(1:nx-j+1,j) = x(j:nx);
		}
	} else {
		ny = length(y);
		H = Z(nx,ny);
		if (x(nx) != y(1)) {
			warning("hankel(): Anti-diagonal confliction.\n");
		}
		for (j = 1; j <= ny; j++) {
			if (j <= nx) {
				H(:,j) = [[x(j:nx)]
						  [y(2:j) ]](1:nx);
			} else {
				H(:,j) = y(j-nx+1:j);
			}
		}
	}
	return H;
}