/*****************************************************************************
*                                                                            *
*   CPLXPAIR.SPL Copyright (C) 2003-2015 DSP Development Corporation         *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Orders complex values by complex conjugates                 *
*                                                                            *
*   Revisions:   15 Oct 2003  RRR  Creation                                  *
*                12 Jun 2015  RRR  conjugate sorting, negative first         *
*                                                                            *
*****************************************************************************/

#if @HELP_CPLXPAIR

    CPLXPAIR

    Purpose: Groups complex values by conjugate pairs.

    Syntax:  CPLXPAIR(z, tol)

                 z - A series, the complex values.

               tol - Optional. A real, the tolerance to determine a purely 
                     real value. Defaults to 100 * eps.
           
    Returns: A series, the ordered paired complex conjugate values.

    Example:

             W1: 1..6
             W2: roots(w1)
             W3: cplxpair(w2)

             W2 == {0.551685 + 1.253349i,
                    0.551685 - 1.253349i,
                   -1.491798 + 0.000000i,
                   -0.805786 + 1.222905i,
                   -0.805786 - 1.222905i}

             W3 == {-0.805786 - 1.222905i,
                    -0.805786 + 1.222905i,
                     0.551685 - 1.253349i,
                     0.551685 + 1.253349i,
                    -1.491798 + 0.000000i}


             W1 represents the polynomial:

             y(x) = x^5 + 2x^4 + 3x^3 + 4x^2 + 5x + 6

             W2 contains the roots of y(x).

             W3 sorts the roots in ascending order with the purely real
             root at the end of the list.

    Remarks:
             CPLXPAIR groups complex conjugate pairs of complex numbers
             in ascending order. The value with a negative imaginary
             part precedes the paired value with a positive imaginary part. 
             Purely real values are placed at the end of the list.

             The tolerance tol specifies the threshold of the absolute
             value of the imaginary part below which the value is
             considered purely real.

             An error occurs if the values cannot be paired.

    See Also:
             Cartesian
             Imaginary
             Real
             Sort
#endif


/* return complex pairs */
ITERATE cplxpair(z, tol)
{
        local y, q, idx, n, i;

        if (argc < 1 || argc > 2) 
        {
                error("(z, n) = cplxpair(z [, tol]);"); 
        }

        if (argc < 2) tol = 100*eps;

        y = zeros(size(z));

        if (length(z) == 0) return({});

        // sort in ascending order
        idx = grade(real(z), 1);
        z   = z[idx];

        // purely real values at end of list
        idx = find(mag(imag(z)) / (mag(z)+realmin) < tol);
        n   = length(z) - length(idx);

        if (length(idx) > 0) 
        {
                y[n+1..length(z)] = z[idx];
                z[idx] = {};
        }

        // put remaining values and conjugates at the start of list

        if (n > 0) 
        {
                loop (i = 1..2..n)
                {
                        if (i+1 > n) 
                        {
                                error("cplxpair - could not pair all complex numbers");
                        }

                        (v, idx) = min(mag(z[i+1..n] - conj(z[i])));

                        if (v > tol) 
                        {
                                error ("cplxpair - could not pair all complex numbers");
                        }

                        if (imag(z[i]) < 0) 
                        {
                                y[{i, i+1}] = z[{i, idx+i}];
                        }
                        else 
                        {
                                y[{i, i+1}] = z[{idx+i, i}];
                        }

                        z[idx+i] = z[i+1];
                }
        }

        return(y, n);
}