CMSCOPE changed

1 files changed, 84 insertions, 0 deletions

diff --git a/modules/cacsd/macros/tf2ss.sci b/modules/cacsd/macros/tf2ss.sci
new file mode 100755
index 000000000..c953185c7
--- /dev/null
+++ b/modules/cacsd/macros/tf2ss.sci
@@ -0,0 +1,84 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - INRIA - Serge Steer
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [sl]=tf2ss(h,tol)
+    // Transfer function to state-space.
+    //Syntax : sl=tf2ss(h [,tol])
+    // h = transfer matrix
+    // sl = linear system in state-space representation (syslin list)
+    //!
+
+    [lhs,rhs]=argn(0)
+
+    if type(h)<=2 then
+        sl=syslin([],[],[],[],h);
+        return;
+    end
+    [num,den]=h(["num","den"]);
+    //
+    [nd,md]=size(den)
+    a=[];b=[];c=[];d=[];n1=0; // s=[]
+    for k=1:md
+        for l=1:nd
+            [r,q]=pdiv(num(l,k),den(l,k));
+            dk(l)=q;
+            num(l,k)=r;
+        end
+        if nd<>1 then
+            [nk,pp]=cmndred(num(:,k),den(:,k));
+        else
+            pp=den(k);nk=num(k);
+        end;
+
+        slk=cont_frm(nk,pp);
+        [ak,bk,ck,dk1]=slk(2:5);
+        // [s sk]
+        [n2,m2]=size(bk);
+        if n2<>0 then
+            a(n1+1:n1+n2,n1+1:n1+n2)=ak;
+            b(n1+1:n1+n2,k)=bk;
+            c=[c ck];
+            n1=n1+n2;
+        else
+            if n1<>0 then b(n1,k)=0;end
+        end;
+        d=[d dk];
+    end;
+
+    if degree(d)==0 then d=coeff(d),end
+
+    if n1<>0 then
+        nrmb=norm(b,1);
+        nrmc=norm(c,1);
+        fact=sqrt(nrmc*nrmb);
+        b=b*fact/nrmb;
+        c=c*fact/nrmc;
+        atmp=a;
+        [a,u]=balanc(a);
+        //next lines commented out to fix bug 3796
+        //   if rcond(u)< %eps*n1*n1*100
+        //     nn=size(a,1);u=eye(nn,nn);a=atmp;
+        //   end
+
+        //apply transformation u without matrix inversion
+        [k,l]=find(u<>0) //get the permutation
+        u=u(k,l);c=c(:,k)*u; b=diag(1 ./diag(u))*b(k,:);
+
+        if rhs<2 then
+            [no,u]=contr(a',c',%eps);
+        else
+            [no,u]=contr(a',c',tol);
+        end
+        u=u(:,1:no);
+        a=u'*a*u;b=u'*b;c=c*u;
+        sl=syslin(h("dt"),a,b,c,d);
+    else
+        sl=syslin(h("dt"),[],[],[],d)
+    end
+endfunction