CMSCOPE changed

1 files changed, 49 insertions, 0 deletions

diff --git a/modules/cacsd/macros/kpure.sci b/modules/cacsd/macros/kpure.sci
new file mode 100755
index 000000000..428b6493d
--- /dev/null
+++ b/modules/cacsd/macros/kpure.sci
@@ -0,0 +1,49 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - Serge Steer
+//
+// 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.1-en.txt
+
+function [K,R]=kpure(sl,eps)
+    //if sl is a transfert function N(S)/D(s) kpure looks for K producing
+    //pure  imaginary roots for
+    //  P(s)=D(s)+K*N(s)
+    //There is a pair of pure imaginary poles if and only if
+    //  P(%i*q)=0  (1)
+    // and
+    //  P(-%i*q)=0 (2)
+    // because N and D are polynomials with real coefficients.
+
+    //Author: Serge Steer, INRIA
+    y=[];R=[];
+    msg=_("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n")
+    if argn(2)==1 then eps=1e-6,end
+    if size(eps,"*")==2 then  eps=eps(2),end //compatibility
+    select typeof(sl)
+    case "rational" then
+        if size(sl.num,"*") <> 1 then
+            error(msprintf(msg,"kpure",1))
+        end
+    case "state-space" then
+        if size(sl.D,"*") <> 1 then
+            error(msprintf(msg,"kpure",1))
+        end
+        sl=ss2tf(sl)
+    else
+        error(msprintf(msg,"kpure",1))
+    end
+
+    //(1) give K(s)=-D(s)/N(s)
+    s=poly(0,varn(sl.den))
+    K=-sl.den/sl.num;
+    // replace K by the previous value in (2) and find the roots
+    s=roots(numer(horner(sl.den,-s)+K*horner(sl.num,-s)),"e");
+    //retain pure imaginary roots
+    s=imag(s(abs(real(s))<eps));
+    R=(s(s>0).'*%i);
+    //find the K(s) values K(s)=-D(s)/N(s)
+    K=-real(freq(sl.den,sl.num,R))
+endfunction