diff options
Diffstat (limited to 'macros/mpoles.sci')
| -rw-r--r-- | macros/mpoles.sci | 194 |
1 files changed, 99 insertions, 95 deletions
diff --git a/macros/mpoles.sci b/macros/mpoles.sci index bba9997..4556120 100644 --- a/macros/mpoles.sci +++ b/macros/mpoles.sci @@ -1,124 +1,128 @@ // Copyright (C) 2018 - IIT Bombay - FOSSEE -// // This file must be used under the terms of the CeCILL. // This source file is licensed as described in the file COPYING, which // you should have received as part of this distribution. The terms // are also available at // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt -// Author:[insert name] +// Original Source : https://octave.sourceforge.io/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 // Organization: FOSSEE, IIT Bombay // Email: toolbox@scilab.in +function [multp, idxp] = mpoles (p, tol, reorder) + if (nargin < 1) + error("mpoles: Invalid number of input arguments"); + end -function [multp, indx] = mpoles (p, tol, reorder) - -//Calling sequence: -// [multp, idxp] = mpoles (p) -// [multp, idxp] = mpoles (p, tol) -// [multp, idxp] = mpoles (p, tol, reorder) -// Identify unique poles in p and their associated multiplicity. -// -// The output is ordered from largest pole to smallest pole. -// -// If the relative difference of two poles is less than tol then they are -// considered to be multiples. The default value for tol is 0.001. -// -// If the optional parameter reorder is zero, poles are not sorted. -// -// The output multp is a vector specifying the multiplicity of the poles. -// multp(n) refers to the multiplicity of the Nth pole -// p(idxp(n)). -// -// -// - -// test case: - -// p = [2 3 1 1 2]; -// [m, n] = mpoles (p) -// n = -// 2. -// 5. -// 1. -// 4. -// 3. -// m = -// 1. -// 1. -// 2. -// 1. -// 2. -// - - -[nargout,nargin]=argn(); - - if (nargin < 1 | nargin > 3) - error("wrong number of input arguments"); + if ~( type(p)== 1) + error ("mpoles: P must be a single or double floating point vector"); end - if (nargin < 2 | isempty (tol)) - tol = 0.001; - end + if (nargin < 2 || isempty (tol)) + tol = 0.001; + elseif (~(isscalar (tol) && isreal (tol) && tol > 0)) + error ("mpoles: TOL must be a real scalar greater than 0"); + end - if (nargin < 3 | isempty (reorder)) - reorder = %t; - end + if (nargin < 3 || isempty (reorder)) + reorder = %t; + elseif (~(isscalar (reorder) && isreal (reorder))) + error ("mpoles: REORDER must be a numeric or logical scalar"); + end Np = length (p); - - // Force the poles to be a column vector. - - p = p(:); - - // Sort the poles according to their magnitidues, largest first. + p = p(:); // force poles to be a column vector if (reorder) - // Sort with smallest magnitude first. - [p, ordr] = gsort (p,"r","i"); - // Reverse order, largest maginitude first. - n = Np:-1:1; - p = p(n); - ordr = ordr(n); + //// sort with largest magnitude first + [_, order] = gsort (abs(p)); + p = p(order); else - ordr = 1:Np; + order = (1:Np).'; end - // Find pole multiplicty by comparing the relative differnce in the - // poles. + //// Create vector of tolerances for use in algorithm. + vtol = zeros (Np, 1, typeof(p)); + p_nz = (p ~= 0); // non-zero poles + vtol(~p_nz) = tol; // use absolute tolerance for zero poles - multp = zeros (Np, 1); - indx = []; + //// Find pole multiplicity by comparing relative difference of poles. + multp = zeros (Np, 1, typeof(p)); + idxp = []; n = find (multp == 0, 1); - - - - - - while (n) - dp = abs (p-p(n)); - if (p(n) == 0.0) - if (or (abs (p) > 0 & (abs(p)<%inf))) - p0 = mean (abs (p(abs (p) > 0 & abs(p)<%inf))); - else - p0 = 1; - end - else - p0 = abs (p(n)); + dp = abs (p - p(n)); + vtol(p_nz) = tol * abs (p(n)); + k = find (dp < vtol); + //// Poles can only be members of one multiplicity group. + if (length (idxp)) + k = k(~ismember (k, idxp)); end - k = find (dp < tol * p0)'; - // Poles can only be members of one multiplicity group. -// if (length(indx)) -// mk=members(k,indx); -// k = (~ bool2s(mk~=0)); -// end m = 1:length (k); - multp(k) = m'; - indx = [indx; k]; + multp(k) = m; + // disp("k") + // disp(k) + // disp("idxp") + // disp(idxp) + idxp = [idxp; k(:)]; n = find (multp == 0, 1); end - multp = multp(indx); - indx = ordr(indx); + multp = multp(idxp); + idxp = order(idxp); +endfunction +function y = ismember(a, b) // FIXME : do this in a more efficient . + y = zeros(size(a,1),size(a,2)); + for i = 1:length(a) + for j = 1:length(b) + if a(i) == b(j) + y(i) = 1; + break; + end + end + end endfunction + + +/* + + // test + [mp, ip] = mpoles ([0 0], 0.01); + assert_checkequal (mp, [1; 2]); + + // test + [mp, ip] = mpoles ([-1e4, -0.1, 0]); + assert_checkequal (mp, [1; 1; 1]); + assert_checkequal (ip, [1; 2; 3]); + +// Test single inputs +// test + [mp, ip] = mpoles ([-1e4, -0.1, 0]); + assert_checkequal (mp,[1; 1; 1]); + assert_checkequal (ip, [1; 2; 3]); + +// Test relative tolerance criteria +// test + [mp, ip] = mpoles ([1, 1.1, 1.3], .1/1.1); + assert_checkequal (mp, [1; 1; 1]); + [mp, ip] = mpoles ([1, 1.1, 1.3], .1/1.1 + %eps); + assert_checkequal (mp, [1; 1; 2]); + +// Test absolute tolerance criteria with a zero pole +// test + [mp, ip] = mpoles ([0, -0.1, 0.3], .1); + assert_checkequal (mp, [1; 1; 1]); + [mp, ip] = mpoles ([0, -0.1, 0.3], .1 + %eps); + assert_checkequal (mp, [1; 1; 2]); + +//// Test input validation +error <Invalid call> mpoles () +error <P must be a single or double floating point vector> mpoles (uint8 (1)) +error <TOL must be a real scalar greater than 0> mpoles (1, [1, 2]) +error <TOL must be a real scalar greater than 0> mpoles (1, 1i) +error <TOL must be a real scalar greater than 0> mpoles (1, 0) +error <REORDER must be a numeric or logical scalar> mpoles (1, 1, [1, 2]) +error <REORDER must be a numeric or logical scalar> mpoles (1, 1, {1}) + +*/
\ No newline at end of file |