/*  -*- MaTX -*-
 *
 *  NAME
 *      unwrap()     - Correct phase angles for all elements
 *      unwrap_col() - Correct phase angles for each columun
 *      unwrap_row() - Correct phase angles for each row
 *
 *  SYNOPSIS
 *      ph2 = unwrap(ph)
 *         Array ph2;
 *         Array ph;
 *
 *      ph2 = unwrap(ph, tol)
 *         Array ph2;
 *         Array ph;
 *         Real tol;
 *
 *      ph2 = unwrap_col(ph)
 *         Array ph2;
 *         Array ph;
 *
 *      ph2 = unwrap_col(ph, tol)
 *         Array ph2;
 *         Array ph;
 *         Real tol;
 *
 *      ph2 = unwrap_row(ph)
 *         Array ph2;
 *         Array ph;
 *
 *      ph2 = unwrap_row(ph, tol)
 *         Array ph2;
 *         Array ph;
 *         Real tol;
 *
 *  DESCRIPTION
 *      unwrap(ph) unwraps phase angles by adding +-2*PI to ph so that 
 *      they will have a smooth transition. The phase must be in radians. 
 *
 *      unwrap(ph,tol) uses tol to check the gap between each phase.
 *
 *      When ph is a matrix,
 *
 *        unwrap_col(ph) corrects the phase angles down each column. 
 *        unwrap_row(ph) corrects the phase angles down each row. 
 *
 *      unwrap(ph) works for all elements. The result matrix is 
 *      the same size as the given matrix.
 *
 *  SEE ALSO
 *      angle and abs.
 */

Func Array unwrap(ph, tol, ...)
	Array ph;
	Real tol;
{
	Array x;
	Integer m, n;

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

	m = Rows(ph);
	n = Cols(ph);
	
	if (nargs == 1) {
		if (m == 1 || n == 1) {
			return unwrap_col(ph);
		} else {
			x = unwrap_col(reshape(ph, m*n, 1));
			return reshape(x, m, n);
		}
	} else {
		if (m == 1 || n == 1) {
			return unwrap_col(ph, tol);
		} else {
			x = unwrap_col(reshape(ph, m*n, 1), tol);
			return reshape(x, m, n);
		}
	}
}

Func Array unwrap_row(ph, tol, ...)
	Array ph;
	Real tol;
{
	error(nargchk(1, 2, nargs, "unwrap_row"));

	if (nargs == 1) {
		return trans(unwrap_col(trans(ph)));
	} else {
		return trans(unwrap_col(trans(ph), tol));

	}
}

Func Array unwrap_col(ph_, tol, ...)
	Array ph_;
	Real tol;
{
	Integer i,j,k,m,mo,no,n;
	Array ph,ph2,p,pd;
	Index idx;

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

	if (nargs == 1) {
		tol = PI*170.0/180.0;
	}
	
	ph = ph_;

	{m,n} = size(ph);
	mo = m;
	no = n;
	if (mo < no) {
		ph = trans(ph);
		{m,n} = size(ph);
	}

	ph2 = ph;

	for (j = 1; j <= n; j++) {
		p = ph(:,j);
		pd = diff(p);
		idx = find(abs(pd) > tol);
		for (i = 1; i <= length(idx); i++) {
			k = idx(i);
			p(k+1:m) = p(k+1:m) .- 2*PI*sgn(pd(k));
		}
		ph2(:,j) = p;
	}

	if (mo < no) {
		ph2 = trans(ph2);
	}
	return ph2;
}