1// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2// Copyright (C) 1996 - 2016 - INRIA - Serge Steer
3// Copyright (C) 2012 - 2016 - Scilab Enterprises
4//
5// This file is hereby licensed under the terms of the GNU GPL v2.0,
6// pursuant to article 5.3.4 of the CeCILL v.2.1.
7// This file was originally licensed under the terms of the CeCILL v2.1,
8// and continues to be available under such terms.
9// For more information, see the COPYING file which you should have received
10// along with this program.
11
12function [kp,s]=krac2(sys)
13    //The denominator of the closed loop system is den(s)+K*num(s). So the
14    //  the closed loops poles verify K(s)=-den(s)/num(s)
15    //The real axis breakaway points occurs at the extrema of the K(s)
16    // so at the point where K'=dK/ds = 0
17    // K'=-(den'*num-den*num')/num^2
18    // K'= 0 --> den'*num-den*num'=0
19    //  http://www.scribd.com/doc/21374148/An-Introduction-to-Control-Systems
20    if and(typeof(sys)<>["state-space" "rational" "zpk"]) then
21        ierr=execstr("[kp,s]=%"+typeof(sys,"overload")+"_krac2(sys)","errcatch")
22        if ierr<>0 then
23            error(msprintf(_("%s: Wrong type for input argument #%d: Linear dynamical system expected.\n"),"krac2",1))
24        end
25        return
26    end
27
28    if size(sys,"*")<>1 then
29        error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"krac2",1))
30    end
31
32    if typeof(sys)=="state-space" then
33        sys=ss2tf(sys)
34    elseif typeof(sys)=="zpk" then
35        sys=zpk2tf(sys)
36    end
37    num=sys.num;
38    den=sys.den;
39
40    s=roots(derivat(num)*den-derivat(den)*num,"e")
41    //collect the real roots only
42    i=find(abs(imag(s))<=10*%eps)
43    if i==[] then kp=[],s=[];return,end
44    s=real(s(i))'
45    s=s(horner(num,s)<>0);
46    kp=-real(freq(den,num,real(s)));
47    i=find(kp>=0);
48    kp=kp(i)
49    s=s(i)
50endfunction
51
52