1// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2// Copyright (C) INRIA -
3//
4// Copyright (C) 2012 - 2016 - Scilab Enterprises
5//
6// This file is hereby licensed under the terms of the GNU GPL v2.0,
7// pursuant to article 5.3.4 of the CeCILL v.2.1.
8// This file was originally licensed under the terms of the CeCILL v2.1,
9// and continues to be available under such terms.
10// For more information, see the COPYING file which you should have received
11// along with this program.
12
13function [n]=linf(g,eps,tol)
14    //linf(g,[eps],[tol])  L_infinity norm
15    //    n=sup [sigmax(g(jw)] (sigmax largest singular value).
16    //       w
17    //-- g is a syslin system.
18    //-- eps is error tolerance on n.
19    //-- tol threshold for imaginary axis poles.
20    // See also: h_norm
21    //!
22    if type(g)==1,if norm(g)==0,n=0,return,end,end,
23
24    if and(typeof(g)<>["rational","state-space"]) then
25        error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n"),"linf",1))
26    end
27    if g.dt<>"c"&g.dt<>[] then
28        error(msprintf(gettext("%s: Wrong value for input argument #%d: Continuous time system expected.\n"),"linf",1))
29    end
30    g.dt="c"
31    if typeof(g)=="rational" then g=tf2ss(g),end
32
33    [lhs,rhs]=argn(0),
34    select rhs,
35    case 1 then eps=1e-7,tol=1000*%eps,
36    case 2 then tol=1000*%eps,
37    end,
38    [a,b,c,d]=g(2:5),[t,t]=size(a),
39    p=ctr_gram(g),q=obs_gram(g);
40    //Algorithm:
41    //----------
42    //1. min , max.
43    //----------------------------------
44    [slp,slm]=dtsi(g),
45    if slp==0 then pp=0,qq=0,tp=1,
46        pm=p,qm=q,tm=t,
47    else
48        if slm==0 then pm=0,qm=0,tm=1,
49            pp=p,qq=q,tp=t,
50        else
51            [tp,tp]=size(slp(2)),[tm,tm]=size(slm(2)),
52            pp=ctr_gram(slp),qq=obs_gram(slp),
53            pm=ctr_gram(slm),qm=obs_gram(slm),
54        end,
55    end,
56    hsvp=sqrt(spec(pp*qq)),hsvp=gsort(real(hsvp)),
57    hsvm=sqrt(spec(pm*qm)),hsvm=gsort(real(hsvm)),
58    gl=max([norm(d),hsvp(tp),hsvm(tm)]),
59    gu=norm(d)+2*(sum(hsvp)+sum(hsvm)),
60    //2. binary search
61    //----------------------
62    while gu-gl>2*eps*gl,
63        x=(gl+gu)/2,
64        r=d'*d-(x*x)*eye(),s=d*d'-(x*x)*eye(),
65        mx=[a-b/r*d'*c, -x*b/r*b'; ..
66        x*c'/s*c,   -a'+c'*d/r*b'],
67        mp=abs(real(spec(mx))),mp=min(mp),
68        if mp>tol then gu=x, else gl=x, end,
69    end;
70    n=(gu+gl)/2
71endfunction
72