bug 3796 fix: pb with tf2ss
[scilab.git] / scilab / modules / cacsd / macros / tf2ss.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2010 - INRIA - Serge Steer
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 [sl]=tf2ss(h,tol)
11 // Transfer function to state-space.
12 //Syntax : sl=tf2ss(h [,tol])
13 // h = transfer matrix
14 // sl = linear system in state-space representation (syslin list)
15 //!
16
17   [lhs,rhs]=argn(0)
18
19   if type(h)<=2 then
20     sl=syslin([],[],[],[],h);
21     return;
22   end
23   [num,den]=h(['num','den']);
24   //
25   [nd,md]=size(den)
26   a=[];b=[];c=[];d=[];n1=0; // s=[]
27   for k=1:md
28     for l=1:nd
29       [r,q]=pdiv(num(l,k),den(l,k));
30       dk(l)=q;
31       num(l,k)=r;
32     end
33     if nd<>1 then
34       [nk,pp]=cmndred(num(:,k),den(:,k));
35     else
36       pp=den(k);nk=num(k);
37     end;
38
39     slk=cont_frm(nk,pp);
40     [ak,bk,ck,dk1]=slk(2:5);
41     // [s sk]
42     [n2,m2]=size(bk);
43     if n2<>0 then
44       a(n1+1:n1+n2,n1+1:n1+n2)=ak;
45       b(n1+1:n1+n2,k)=bk;
46       c=[c ck];
47       n1=n1+n2;
48     else
49       if n1<>0 then b(n1,k)=0;end
50     end;
51     d=[d dk];
52   end;
53
54   if degree(d)==0 then d=coeff(d),end
55
56   if n1<>0 then
57     nrmb=norm(b,1);
58     nrmc=norm(c,1);
59     fact=sqrt(nrmc*nrmb);
60     b=b*fact/nrmb;
61     c=c*fact/nrmc;
62     atmp=a;
63     [a,u]=balanc(a);
64     //next lines commented out to fix bug 3796
65     //   if rcond(u)< %eps*n1*n1*100
66     //     nn=size(a,1);u=eye(nn,nn);a=atmp;
67     //   end
68     
69     //apply transformation u without matrix inversion
70     [k,l]=find(u<>0) //get the permutation
71     u=u(k,l);c=c(:,k)*u; b=diag(1 ./diag(u))*b(k,:);
72
73     if rhs<2 then
74       [no,u]=contr(a',c',%eps);
75     else
76       [no,u]=contr(a',c',tol);
77     end
78     u=u(:,1:no);
79     a=u'*a*u;b=u'*b;c=c*u;
80     sl=syslin(h('dt'),a,b,c,d);
81   else
82     sl=syslin(h('dt'),[],[],[],d)
83   end
84 endfunction