code improved: roots(p), replaced by roots(p,'e') + unuseful test on continuous time...
[scilab.git] / scilab / modules / cacsd / macros / krac2.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 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 [kp,s]=krac2(sys)
11 //The denominator of the closed loop system is den(s)+K*num(s). So the
12 //  the closed loops poles verify K(s)=-den(s)/num(s)
13 //The real axis breakaway points occurs at the extrema of the K(s)
14 // so at the point where K'=dK/ds = 0
15 // K'=-(den'*num-den*num')/num^2
16 // K'= 0 --> den'*num-den*num'=0
17 //  http://www.scribd.com/doc/21374148/An-Introduction-to-Control-Systems
18   select typeof(sys)
19   case 'rational' then
20   case 'state-space' then
21     sys=ss2tf(sys);
22   else
23      error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n"),"krac2",1))
24   end
25   if size(sys,'*')<>1 then
26     error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"krac2",1))
27   end
28   num=sys.num
29   den=sys.den
30   s=roots(derivat(num)*den-derivat(den)*num,'e')
31   //collect the real roots only
32   i=find(abs(imag(s))<=10*%eps)
33   if i==[] then kp=[],s=[];return,end
34   s=s(i)'
35   kp=-real(freq(den,num,real(s)))
36   i=find(kp>=0)
37   kp=kp(i)
38   s=s(i)
39 endfunction