summaryrefslogtreecommitdiff
path: root/macros/cplxreal.sci
diff options
context:
space:
mode:
Diffstat (limited to 'macros/cplxreal.sci')
-rw-r--r--macros/cplxreal.sci97
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