56da3ce3d0025f11788eec49a5cfc690707b6ca8
[scilab.git] / scilab / modules / cacsd / macros / minreal.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 [a,b,c]=minreal(a,b,c,domaine,tol)
14     //
15
16     [lhs,rhs]=argn(0)
17     select typeof(a)
18     case "state-space" then
19         if lhs<>1 then
20             error(msprintf(gettext("%s: Wrong number of output arguments: %d expected.\n"),"minreal",1)),
21         end;
22         select rhs
23         case 1 then
24             istol = %f;
25         case 2 then
26             istol = %t,
27             tol = b,
28         else
29             error(msprintf(gettext("%s: Wrong number of input arguments: %d or %d expected.\n"),"minreal",1,2)),
30         end;
31         [a,b,c,d,x0,dom] = a(2:7);
32         if dom == [] then
33             error(96,1);
34         end
35         domaine="c";
36         if dom<>"c" then
37             domaine="d";
38         end
39     case "constant" then
40         if lhs<>3 then
41             error(msprintf(gettext("%s: Wrong number of output arguments: %d expected.\n"),"minreal",3)),
42         end;
43         select rhs
44         case 4 then istol = %f
45         case 5 then istol = %t,
46         else
47             error(msprintf(gettext("%s: Wrong number of input arguments: %d or %d expected.\n"),"minreal",4,5));
48         end;
49     else
50         error(91,1);
51     end;
52     //
53     wc = lyap(a', -b*b', domaine);
54     wo = lyap(a, -c'*c, domaine);
55     if ~istol then
56         [r,n] = equil1(wc,wo);
57     else
58         [r,n] = equil1(wc,wo,tol);
59     end;
60     n1 = n(1);
61     ri = inv(r);
62     r = r(1:n1,:);
63     ri = ri(:,1:n1);
64     a = r*a*ri;
65     b = r*b;
66     c = c*ri;
67     if lhs == 1 then
68         a = syslin(dom,a,b,c,d,r*x0);
69     end
70 endfunction