Merge pull request #17 from avinashlalotra/masterHEADmaster

Abinash's Work

1 files changed, 107 insertions, 0 deletions

diff --git a/macros/cplxpair.sci b/macros/cplxpair.sci
new file mode 100644
index 0000000..16fceac
--- /dev/null
+++ b/macros/cplxpair.sci
@@ -0,0 +1,107 @@
+/*Description
+        Sort the numbers z into complex conjugate pairs ordered by increasing real part.
+        The negative imaginary complex numbers are placed first within each pair. All real numbers (those with abs (imag (z) / z) < tol) are placed after the complex pairs.
+        and the resulting tolerance for a given complex pair is 100 * eps (abs (z(i))).
+        Signal an error if some complex numbers could not be paired. Signal an error if all complex numbers are not exact conjugates (to within tol).
+        Note that there is no defined order for pairs with identical real parts but differing imaginary parts
+Calling Sequence
+         cplxpair (z)
+         cplxpair (z, tol)
+         cplxpair (z, tol, dim)
+Parameters:
+         z is a matrix or vector.
+         tol is a weighting factor which determines the tolerance of matching. The default value is 100.
+         By default the complex pairs are sorted along the first non-singleton dimension of z. If dim is specified, then the complex pairs are sorted along this dimension.
+Dependencies: ipermute
+Example:
+        The following code
+        [ cplxpair(exp(2*%i*%pi*[0:4]'/5)), exp(2*%i*%pi*[3; 2; 4; 1; 0]/5) ]
+        Produces the following output
+                ans =
+             -0.80902 - 0.58779i  -0.80902 - 0.58779i
+             -0.80902 + 0.58779i  -0.80902 + 0.58779i
+              0.30902 - 0.95106i   0.30902 - 0.95106i
+              0.30902 + 0.95106i   0.30902 + 0.95106i
+              1.00000 + 0.00000i   1.00000 + 0.00000i
+*/
+function zsort = cplxpair (z, tol, dim)
+  if (nargin < 1)
+    error("Invalid inputs");
+  end
+  // default  double
+  realmin = 2.2251e-308
+  if (isempty (z))
+    zsort = zeros (size (z,1) , size (z,2));
+    return;
+  end
+  if (nargin < 2 || isempty (tol))
+    tol = 100* %eps;
+  elseif (~ isscalar (tol) || tol < 0 || tol >= 1)
+    error ("cplxpair: TOL must be a scalar number in the range 0 <= TOL < 1");
+  end
+  nd = ndims (z);
+  if (nargin < 3)
+    // Find the first singleton dimension.
+    sz = size (z);
+    dim = find (sz > 1, 1);
+    if isempty(dim) 
+        dim  = 1;
+    end
+  else
+    dim = floor (dim);
+    if (dim < 1 || dim > nd)
+      error ("cplxpair: invalid dimension DIM");
+    end
+  end
+  // Move dimension to analyze to first position, and convert to a 2-D matrix.
+  perm = [dim:nd, 1:dim-1];
+  z = permute (z, perm);
+  sz = size (z);
+  n = sz(1);
+  m = prod (sz) / n;
+  z = matrix (z, n, m);
+  // Sort the sequence in terms of increasing real values.
+  [temp, idx] = gsort (real (z), 1 , "i");
+  z = z(idx + n * ones (n, 1) * [0:m-1]);
+  // Put the purely real values at the end of the returned list.
+  [idxi, idxj] = find (abs (imag (z)) ./ (abs (z) + realmin) <= tol);
+  // Force values detected to be real within tolerance to actually be real.
+  z(idxi + n*(idxj-1)) = real (z(idxi + n*(idxj-1)));
+  //if idxi and idxj are not scalers
+  if ~isscalar(idxi) then
+      v = ones(size(idxi,1),size(idxi,2));
+  else
+      v = 1 ;
+  end
+  q = sparse ([idxi' idxj'], v, [n m]);
+  nr = sum (q, 1);
+  [temp, idx] = gsort (q, 'r','i');
+  midx = idx + size (idx,1) * ones (size (idx,1), 1) * [0:size(idx,2)-1];
+  z = z(midx);
+  zsort = z;
+  // For each remaining z, place the value and its conjugate at the start of
+  // the returned list, and remove them from further consideration.
+  for j = 1:m
+    p = n - nr(j);
+    for i = 1:2:p
+      if (i+1 > p)
+        error ("cplxpair: could not pair all complex numbers");
+      end
+      [v, idx] = min (abs (z(i+1:p,j) - conj (z(i,j))));
+      if (v >= tol * abs (z(i,j)))
+        error ("cplxpair: could not pair all complex numbers");
+      end
+      // For pairs, select the one with positive imaginary part and use it and
+      // it's conjugate, but list the negative imaginary pair first.
+      if (imag (z(i,j)) > 0)
+        zsort([i, i+1],j) = [conj(z(i,j)), z(i,j)];
+      else
+        zsort([i, i+1],j) = [conj(z(idx+i,j)), z(idx+i,j)];
+      end
+      z(idx+i,j) = z(i+1,j);
+    end
+  end
+  // Reshape the output matrix.
+  zsort = ipermute (matrix (zsort, sz), perm);
+endfunction
+