Revert "* Bug #12174 fixed - The function "routh_t" gives incorrect output for all...
[scilab.git] / scilab / modules / cacsd / macros / routh_t.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA - Serge STEER
3 // Copyright (C) 1999 - Lucien.Povy@eudil.fr (to get a good table)
4 // 
5 // This file must be used under the terms of the CeCILL.
6 // This source file is licensed as described in the file COPYING, which
7 // you should have received as part of this distribution.  The terms
8 // are also available at    
9 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10
11 function r=routh_t(h,k)
12 //r=routh_t(h,k) computes routh table of denominator of the
13 //system described by transfer matrix SISO continue h with the
14 //feedback by the gain k
15 //If  k=poly(0,'k') we will have a polynomial matrix with dummy variable
16 //k, formal expression of the Routh table.
17 //r=routh_t(d) computes Routh table of h :attention ! d=denom of system
18
19 //If a zero row appears, it means that there exist a pair of pure imaginary
20 //roots (oscillating system) or symmetric real roots. In this case, the pure imaginary roots are the
21 //imaginary roots of the bisquare polynomial given by the previous row. The
22 //routh table can be continued replacing this row by the derivative of this
23 //polynomial.
24 //see http://www.jdotec.net/s3i/TD_Info/Routh/Routh.pdf for degenerated
25 //cases
26   
27   //see also 
28 //Comments on the Routh-Hurwitz criterion, Shamash, Y.,Automatic Control, IEEE T.A.C
29 //Volume 25, Issue 1, Feb 1980 Page(s): 132 - 133
30   
31   //http://controls.engin.umich.edu/wiki/index.php/RouthStability
32   [lhs,rhs]=argn(0);
33   h1=h(1);
34   if rhs==2 then
35     if typeof(h)<>'rational' then
36       error(msprintf(gettext("%s: Wrong type for input argument #%d: rational fraction array expected.\n"),"routh_t",1));
37     end
38     [n,d]=h(2:3)
39     if size(n,'*')<>1 then
40       error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"routh_t",1))
41     end
42     nd=max([degree(d) degree(n)])+1;
43     cod=coeff(d,0:nd-1);//coeff du denominateur
44     con=coeff(n,0:nd-1);//coeff du numerateur
45     cobf=cod+k*con //coeff de la boucle fermee
46   else
47     if type(h)>2 then 
48       error(msprintf(gettext("%s: Wrong type for input argument #%d: Polynomial array expected.\n"),"routh_t",1));
49     end
50     if size(h,'*')<>1 then
51       error(msprintf(gettext("%s: Wrong size for input argument #%d: A polynomial expected.\n"),"routh_t",1))
52     end
53
54     nd=degree(h)+1;
55     cobf=coeff(h,0:nd-1)
56   end;
57   //
58   r1=cobf(nd:-2:1);
59   r2=cobf(nd-1:-2:1);
60   ncol=size(r1,'*');
61   if size(r2,'*')<>ncol then r2=[r2,0],end
62   r=[r1;r2]
63   if ncol<2 then r=[],return,end;
64   if rhs==2 then
65
66     for i=3:nd,
67       r(i,1:ncol-1)=[r(i-1,1),-r(i-2,1)]*[r(i-2,2:ncol);r(i-1,2:ncol)]
68     end;
69   else
70     for i=3:nd,
71      if and(r(i-1,:)==0) then
72         naux=nd-i+2 //order of previous polynomial
73         exponents=naux:-2:0
74         ncoeff=size(exponents,'*')
75         r(i-1,1:ncoeff)=r(i-2,1:ncoeff).*exponents //derivative of previous polynomial
76       end
77       if r(i-1,1)==0 then 
78         if rhs==1 then
79           if typeof(r)=='rational' then 
80             //scilab is not able to handle multivariable polynomials
81             r=horner(r,%eps^2); 
82           end
83           r(i-1,1)=poly(0,'eps')
84         else
85           r(i-1,1)=%eps^2,
86         end
87       end
88       r(i,1:ncol-1)=[1.,-r(i-2,1)/r(i-1,1)]*[r(i-2,2:ncol);r(i-1,2:ncol)]
89       
90     end;
91   end;
92 endfunction