c6e79fadee08680296b0a1c192560f9de23bf80f
[scilab.git] / scilab / modules / cacsd / macros / tf2ss.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 [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);return;
21 end
22 [num,den]=h(['num','den']);
23 //
24 [nd,md]=size(den)
25 a=[];b=[];c=[];d=[];n1=0; // s=[]
26 for k=1:md
27   for l=1:nd
28     [r,q]=pdiv(num(l,k),den(l,k));
29     dk(l)=q;
30     num(l,k)=r;
31   end
32   if nd<>1 then 
33     [nk,pp]=cmndred(num(:,k),den(:,k));
34   else
35     pp=den(k);nk=num(k);
36   end;
37   
38   slk=cont_frm(nk,pp);//
39   [ak,bk,ck,dk1]=slk(2:5);
40   // [s sk]
41   [n2,m2]=size(bk);
42   if n2<>0 then
43     a(n1+1:n1+n2,n1+1:n1+n2)=ak;
44     b(n1+1:n1+n2,k)=bk;
45     c=[c ck];
46     n1=n1+n2;
47   else
48     if n1<>0 then b(n1,k)=0;end
49   end;
50   d=[d dk];
51 end;
52
53 if degree(d)==0 then d=coeff(d),end
54
55 if n1<>0 then
56   nrmb=norm(b,1);nrmc=norm(c,1);fact=sqrt(nrmc*nrmb);
57   b=b*fact/nrmb;c=c*fact/nrmc;atmp=a;
58   [a,u]=balanc(a);
59   if rcond(u)< %eps*n1*n1*100
60   nn=size(a,1);u=eye(nn,nn);a=atmp;
61   end
62   //apply transformation u without matrix inversion
63   [k,l]=find(u<>0) //get the permutation
64   u=u(k,l);c=c(:,k)*u; b=diag(1 ./diag(u))*b(k,:);
65   //c=c*u;b=u^(-1)*b; 
66
67   if rhs<2 then 
68     [no,u]=contr(a',c',%eps);
69   else
70     [no,u]=contr(a',c',tol);
71   end
72   u=u(:,1:no);
73   a=u'*a*u;b=u'*b;c=c*u;
74   //[a,u]=balanc(a);
75   
76   //apply transformation u without matrix inversion
77   //[k,l]=find(u<>0) //get the permutation
78   //u=u(k,l);c=c(:,k)*u; b=diag(1 ./diag(u))*b(k,:);
79   //c=c*u;b=u^(-1)*b;
80
81   sl=syslin(h('dt'),a,b,c,d);
82 else
83   sl=syslin(h('dt'),[],[],[],d)
84 end
85 endfunction