/*  -*- MaTX -*-
 *
 *  NAME
 *      ccpair() -  Make complex conjugate pairs
 *
 *  SYNOPSIS
 *      y = ccpair(x)
 *         CoMatrix y;
 *         CoMatrix x;
 *
 *      y = ccpair(x, tol)
 *         CoMatrix y;
 *         CoMatrix x;
 *         Real tol;
 *
 *  DESCRIPTION
 *      ccpair(x) sorts the elements of x so that complex numbers 
 *      are collected into matched pairs of complex conjugates.  The pairs 
 *      are ordered by increasing real part.  All real elements are placed 
 *      after the complex pairs.
 *
 *      ccpair(X,tol) uses tol for the tolerance of comparison.
 *
 *  SEE ALSO
 *      sort
 */

Func CoMatrix ccpair(x_, tol, ...)
	CoMatrix x_;
	Real tol;
{
	Integer k,nc,n,nr;
	Complex j;
	Matrix x2,c,cr,r,xi,tmp,R;
	CoMatrix x,y;
	Index idx_c,xind,idx_re,i1,i2;

	error(nargchk(1, 2, nargs, "ccpair"));

	j = (0,1);

	if (nargs == 1) {
		tol = 100*EPS;
	}

	y = Z(x_);
	if (Rows(x_) < Cols(x_)) {
		x = trans(x_);
	} else {
		x = x_;
	}

	// Find real numbers
	idx_re = find(abs(Im(x)) .<= tol*abs(x));
	nr = length(idx_re);

	if (nr != 0) {
		{R} = sort(Re(x(idx_re)));
		x(idx_re) = [];
	}
	nc = length(x);

	if (nc == 0) {
		y = R;
		return y;
	}
	if (rem(nc,2)) {
		error("ccpair(): Complex numbers can't be paired\n");
	}

	// Upper half plane
	c = Re(x) + j*abs(Im(x));

	// Sort by real part
	{cr, idx_c} = sort(Re(c));
	c = cr + j*Im(c(idx_c));

	// Real parts are at least in pair ?
	if (any(abs(cr(1:2:nc)-cr(2:2:nc)) .> tol*abs(c(1:2:nc)))) {
		error("ccpair(): Complex numbers can't be paired\n");
	}
	x = x(idx_c);

	k = 1;
	while (k <= nc/2) {
		i1 = find(abs(cr(1:2:nc) .- cr(2*k)) .<= tol*abs(c(2*k)));

		if (length(i1) == 0) {
			error("ccpair(): Complex numbers can't be paired\n");
		} else if (length(i1) == 1) {
			// One pair
			if (Re(x(2*k) - conj(x(2*k-1))) > tol*abs(x(2*k))) {
				error("ccpair(): Complex numbers can't be paired\n");
			} else {
				tmp = [Re(x(2*k)) - j*abs(Im(x(2*k)))];
				x(2*k-1:2*k) = [[tmp][conj(tmp)]];
			}
		} else {
			// Multi-pairs with same real part
			i2 = find(abs(cr(1:nc) .- cr(2*k)) .<= tol*abs(c(2*k)));
			
			n = max(1,length(i2));
			x2 = x(i2);
			{xi, xind} = sort(Im(x2));
			x2 = x2(xind);

			if (any(abs(x2-conj(x2(n:-1:1))) .> tol*abs(x2))) {
				error("ccpair(): Complex numbers can't be paired\n");
			} else {
				x(i2(1):2:i2(n-1)) = x2(1:n/2);
				x(i2(2):2:i2(n)) = conj(x2(1:n/2));
			}
		}

		k = k + length(i1);
	}

	if (nc) {
		y(1:nc,1) = x(1:nc);
	}
	if (nr) {
		y(nc+1:nc+nr,1) = R;
	}
	return y;
}