error(number): converting occurrences remaining in all .sce .sci files
[scilab.git] / scilab / modules / cacsd / macros / arhnk.sci
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
13 function [slm]=arhnk(a,ordre,tol)
14
15     [lhs,rhs]=argn(0),
16     if lhs<>1 then
17         msg = _("%s: Wrong number of output arguments: %d to %d expected.\n");
18         error(msprintf(msg, "arhnk", 0, 1));
19     end
20     if typeof(a)<>"state-space" then
21         msg = _("%s: Argument #%d: Linear state space expected.\n");
22         error(msprintf(msg, "arhnk", 1));
23     end;
24     if a.dt<>"c" then
25         msg = _("%s: Wrong type for input argument #%d: In continuous time expected.\n");
26         error(msprintf(msg, "arhnk",1));
27     end
28     select rhs
29     case 2 then istol=0;
30     case 3 then istol=1;
31     end;
32
33     [a,b,c,d,x0,dom]=a(2:7);
34     if(max(real(spec(a)))) > 0 then
35         msg = _("%s: Wrong values for input argument #%d: Stable system expected.\n");
36         error(msprintf(msg, "arhnk", 1));
37     end
38     domaine="c"
39     wc=lyap(a',-b*b',domaine)
40     wo=lyap(a,-c'*c,domaine)
41     if istol==0 then [t,nn]=equil1(wc,wo);
42     else [t,nn]=equil1(wc,wo,tol);
43     end;
44     n1=nn(1);
45     ti=inv(t);a=t*a*ti;b=t*b;c=c*ti
46     wc=t*wc*t';wo=ti'*wo*ti;
47     if ordre>n1 then
48         ordre=n1
49     end;
50     if ordre==n1 then
51         a=a(1:n1,1:n1);b=b(1:n1,:);c=c(:,1:n1);
52         if lhs==1 then a=syslin("c",a,b,c,d,0*ones(n1,1)),end
53         return,
54     end;
55     sigma=wc(ordre+1,ordre+1)
56
57     r=max(n1-ordre-1,1)
58
59     n=n1
60     sel=[1:ordre ordre+r+1:n];seleq=ordre+1:ordre+r
61     b2=b(seleq,:);c2=c(:,seleq);
62     u=-c(:,seleq)*pinv(b(seleq,:)')
63     a=a(sel,sel);b=b(sel,:);c=c(:,sel);
64     wo=wo(sel,sel);wc=wc(sel,sel);
65
66     Gamma=wc*wo-sigma*sigma*eye()
67     a=Gamma\(sigma*sigma*a'+wo*a*wc-sigma*c'*u*b')
68     b1=Gamma\(wo*b+sigma*c'*u)
69     c=c*wc+sigma*u*b';b=b1;
70     d=-sigma*u+d
71     //
72     [n,n]=size(a)
73     [u,m]=schur(a,"c")
74     a=u'*a*u;b=u'*b;c=c*u;
75     if m<n then
76         t=sylv(a(1:m,1:m),-a(m+1:n,m+1:n),-a(1:m,m+1:n),"c")
77         a=a(1:m,1:m)
78         b=b(1:m,:)-t*b(m+1:n,:)
79         c=c(:,1:m)
80     end;
81     //
82     slm=syslin("c",a,b,c,d,0*ones(m,1));
83 endfunction