Imlemented by Abinash Singh During FOSSEE Semester Long Fellowship 2024

1 files changed, 50 insertions, 47 deletions

diff --git a/macros/cplxreal.sci b/macros/cplxreal.sci
index a552998..b06fac3 100644
--- a/macros/cplxreal.sci
+++ b/macros/cplxreal.sci
@@ -1,48 +1,51 @@
-function [zc, zr] = cplxreal (z, thresh)
-//Function to divide vector z into complex and real elements, removing the one of each complex conjugate pair.
-//Calling Sequence
-//[zc, zr] = cplxreal (z, thresh)
-//[zc, zr] = cplxreal (z)
-//zc = cplxreal (z, thresh)
-//zc = cplxreal (z)
-//Parameters 
-//z: vector of complex numbers.
-//thresh: tolerance for comparisons.
-//zc: vector containing the elements of z that have positive imaginary parts.
-//zr: vector containing the elements of z that are real.
-//Description
-//This is an Octave function.
-//Every complex element of z is expected to have a complex-conjugate elsewhere in z. From the pair of complex-conjugates, the one with the negative imaginary part is removed.
-//If the magnitude of the imaginary part of an element is less than the thresh, it is declared as real.  
-//Examples
-//[zc, zr] = cplxreal([1 2 3+i 4 3-i 5])
-//zc =  3 + 1i
-//zr =
-//   1   2   4   5
-funcprot(0);
-lhs = argn(1)
-rhs = argn(2)
-if (rhs < 1 | rhs > 2)
-error("Wrong number of input arguments.")
-end
-
-select(rhs)
-	case 1 then
-		if(lhs==1)
-		zc = callOctave("cplxreal",z)
-		elseif (lhs==2)
-		[zc, zr] = callOctave("cplxreal",z)
-		else
-		error("Wrong number of output argments.")
-		end
-
-	case 2 then
-		if(lhs==1)
-		zc = callOctave("cplxreal",z, thresh)
-		elseif (lhs==2)
-		[zc, zr] = callOctave("cplxreal",z, thresh)
-		else
-		error("Wrong number of output argments.")
-		end
-	end
+/*Description
+        Sort the numbers z into complex-conjugate-valued and real-valued elements.
+        The positive imaginary complex numbers of each complex conjugate pair are returned in zc and the real numbers are returned in zr.
+        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:
+        [zc, zr] = cplxreal (z)
+        [zc, zr] = cplxreal (z, tol)
+        [zc, zr] = cplxreal (z, tol, dim)
+Parameters
+        Inputs
+             z  -  A vector of numbers or Matrix
+            tol - tol is a weighting factor in the range [0, 1) which determines the tolerance of the matching.
+            The default value is 100 * eps and the resulting tolerance for a given complex pair is tol * abs (z(i))).
+            dim - 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.
+        Outputs
+            zc - complex conjugate pair
+            zr -  real numbers
+Example:
+        with 2 real zeros, one of them equal to 0
+        [zc, zr] = cplxreal (roots ([1, 0, 0, 1, 0])) */
+function [zc, zr] = cplxreal (z, tol, dim)
+  if (nargin < 1 || nargin > 3)
+    error("invalid inputs");
+  end
+  if (isempty (z))
+    zc = zeros (size (z,1),size(z,2));
+    zr = zeros (size (z,1),size(z,2));
+    return;
+  end
+  if (nargin < 2 || isempty (tol))
+    tol = 100 * %eps ;
+  end
+  if (nargin >= 3)
+    zcp = cplxpair(z,tol,dim);
+  else
+    zcp = cplxpair (z , tol);
+  end
+  nz = max(size (z) );
+  idx = nz;
+  while ((idx > 0) && (zcp(idx) == 0 || (abs (imag (zcp(idx))) ./ abs (zcp(idx))) <= tol))
+    zcp(idx) = real (zcp(idx));
+    idx = idx - 1;
+  end
+  if (pmodulo (idx, 2) ~= 0)
+    error ("cplxreal: odd number of complex values was returned from cplxpair");
+  end
+  zc = zcp(2:2:idx);
+  zr = zcp(idx+1:nz);
 endfunction