da7d7f02e2e233ff38b4d01def1686882b4ecd4d
[scilab.git] / scilab / modules / cacsd / macros / pfss.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA - 
3 // 
4 // This file must be used under the terms of the CeCILL.
5 // This source file is licensed as described in the file COPYING, which
6 // you should have received as part of this distribution.  The terms
7 // are also available at    
8 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
9
10 function elts=pfss(S,rmax,cord)
11 //Syntax : elts=pfss(S)
12 //Partial fraction decomposition of the linear system S (in state-space form):
13 // elts is the list of linear systems which add up to S
14 // i.e. elts=list(S1,S2,S3,...,Sn) with S1 + S2 +... +Sn = S
15 // Each Si contains some poles of S according to the block-diagonalization
16 // of the A matrix of S.
17 // If S is given in transfer form, it is first converted into state-space
18 // and each subsystem is then converted in transfer form.
19 //!
20   select argn(2)
21   case 1 then
22     rmax=[];cord=[]
23   case 2 the
24     if type(rmax)==10 then cord=rmax;end
25     if type(rmax)==1 then cord=[];end
26   end
27   
28   if and(typeof(S)<>['rational','state-space']) then
29     error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n"),"pfss",1))
30   end
31   if typeof(S)=='rational' then 
32     S=tf2ss(S),flag=%f
33   else
34     flag=%t
35   end
36
37   [f,g,h,dd,dom]=S([2:5,7]);
38    [n,n]=size(f);
39   if rmax==[] then
40     [f,x,bs]=bdiag(f);
41   else
42     [f,x,bs]=bdiag(f,rmax);
43   end
44   h=h*x;g=x\g;
45   k=1;ll=0;
46   elts=list();
47   for l=bs',
48     ind=k:k+l-1;
49     f1l=f(ind,ind);
50     gl=g(ind,:);
51     hl=h(:,ind);
52     elts(ll+1)=syslin('c',f1l,gl,hl)
53     ll=ll+1;k=k+l;
54   end;
55   if argn(2)==2  then
56     select cord
57     case 'c'
58       nb=size(bs,'*');
59       class=[];
60       for k=1:nb
61         oneortwo=bs(k);ss=elts(k);A=ss(2);
62         if oneortwo==1 then 
63           class=[class,real(spec(A))];
64         end
65         if oneortwo>1 then 
66           class=[class,mini(real(spec(A)))];
67         end
68       end;
69       [cl,indi]=sort(-class);
70       elts1=elts;
71       for k=1:size(elts);
72         elts(k)=elts1(indi(k));
73       end
74     case 'd'
75       nb=size(bs,'*');
76       class=[];
77       for k=1:nb
78         oneortwo=bs(k);ss=elts(k);A=ss(2);
79         if oneortwo==1 then 
80           class=[class,abs(spec(A))];
81         end
82         if oneortwo>1 then 
83           class=[class,maxi(abs(spec(A)))];
84         end
85       end;
86       [cl,indi]=sort(-class);
87       elts1=elts;
88       for k=1:size(elts);
89         elts(k)=elts1(indi(k));
90       end
91     end
92   end
93   if type(dd)==1 then
94     if norm(dd,'fro')<>0 then elts(ll+1)=dd,end
95   end
96   if type(dd)==2 then 
97     if norm(coeff(dd),1) > %eps then elts(ll+1)=dd,end
98   end
99   if flag then
100     k=size(elts);
101     for kk=1:k,elts(kk)=ss2tf(elts(kk));end
102   end
103 endfunction