1 files changed, 68 insertions, 0 deletions
diff --git a/modules/cacsd/macros/linf.sci b/modules/cacsd/macros/linf.sci
new file mode 100755
index 000000000..7cb6ef315
--- /dev/null
+++ b/modules/cacsd/macros/linf.sci
@@ -0,0 +1,68 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA -
+//
+// 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 [n]=linf(g,eps,tol)
+    //linf(g,[eps],[tol])  L_infinity norm
+    //    n=sup [sigmax(g(jw)] (sigmax largest singular value).
+    //       w
+    //-- g is a syslin system.
+    //-- eps is error tolerance on n.
+    //-- tol threshold for imaginary axis poles.
+    // See also: h_norm
+    //!
+    if type(g)==1,if norm(g)==0,n=0,return,end,end,
+
+    if and(typeof(g)<>["rational","state-space"]) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n"),"linf",1))
+    end
+    if g.dt<>"c"&g.dt<>[] then
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Continuous time system expected.\n"),"linf",1))
+    end
+    g.dt="c"
+    if typeof(g)=="rational" then g=tf2ss(g),end
+
+    [lhs,rhs]=argn(0),
+    select rhs,
+    case 1 then eps=1e-7,tol=1000*%eps,
+    case 2 then tol=1000*%eps,
+    end,
+    [a,b,c,d]=g(2:5),[t,t]=size(a),
+    p=ctr_gram(g),q=obs_gram(g);
+    //Algorithm:
+    //----------
+    //1. min , max.
+    //----------------------------------
+    [slp,slm]=dtsi(g),
+    if slp==0 then pp=0,qq=0,tp=1,
+        pm=p,qm=q,tm=t,
+    else
+        if slm==0 then pm=0,qm=0,tm=1,
+            pp=p,qq=q,tp=t,
+        else
+            [tp,tp]=size(slp(2)),[tm,tm]=size(slm(2)),
+            pp=ctr_gram(slp),qq=obs_gram(slp),
+            pm=ctr_gram(slm),qm=obs_gram(slm),
+        end,
+    end,
+    hsvp=sqrt(spec(pp*qq)),hsvp=gsort(real(hsvp)),
+    hsvm=sqrt(spec(pm*qm)),hsvm=gsort(real(hsvm)),
+    gl=max([norm(d),hsvp(tp),hsvm(tm)]),
+    gu=norm(d)+2*(sum(hsvp)+sum(hsvm)),
+    //2. binary search
+    //----------------------
+    while gu-gl>2*eps*gl,
+        x=(gl+gu)/2,
+        r=d'*d-(x*x)*eye(),s=d*d'-(x*x)*eye(),
+        mx=[a-b/r*d'*c, -x*b/r*b'; ..
+        x*c'/s*c,   -a'+c'*d/r*b'],
+        mp=abs(real(spec(mx))),mp=min(mp),
+        if mp>tol then gu=x, else gl=x, end,
+    end;
+    n=(gu+gl)/2
+endfunction