code improved: roots(p), replaced by roots(p,'e') + unuseful test on continuous time...
Serge Steer [Mon, 29 Mar 2010 13:41:04 +0000 (15:41 +0200)]
scilab/modules/cacsd/macros/krac2.sci

index 0247c94..e4bc115 100644 (file)
@@ -7,32 +7,33 @@
 // are also available at    
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
-function kp=krac2(n)
-
-  select typeof(n)
+function [kp,s]=krac2(sys)
+//The denominator of the closed loop system is den(s)+K*num(s). So the
+//  the closed loops poles verify K(s)=-den(s)/num(s)
+//The real axis breakaway points occurs at the extrema of the K(s)
+// so at the point where K'=dK/ds = 0
+// K'=-(den'*num-den*num')/num^2
+// K'= 0 --> den'*num-den*num'=0
+//  http://www.scribd.com/doc/21374148/An-Introduction-to-Control-Systems
+  select typeof(sys)
   case 'rational' then
-    [n,d,dom]=n(['num','den','dt'])
   case 'state-space' then
-    n=ss2tf(n);[n,d,dom]=n(['num','den','dt'])
+    sys=ss2tf(sys);
   else
      error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n"),"krac2",1))
   end
-  if dom<>'c' then 
-    error(msprintf(gettext("%s: Wrong value for input argument #%d: Continuous time system expected.\n"),"krac2",1))
-  end
-  if size(n,'*')<>1 then
+  if size(sys,'*')<>1 then
     error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"krac2",1))
   end
-
-  x=[];
-  q1=derivat(n/d);s=roots(q1(2));
-  //
-  for a=s',
-    if abs(imag(a))<=10*%eps then 
-      x=[x;a],
-    end,
-  end
-  //x(x==0)=[]
-  if x==[] then;return,end
-  kp=gsort(-real(freq(d,n,real(x))))
+  num=sys.num
+  den=sys.den
+  s=roots(derivat(num)*den-derivat(den)*num,'e')
+  //collect the real roots only
+  i=find(abs(imag(s))<=10*%eps)
+  if i==[] then kp=[],s=[];return,end
+  s=s(i)'
+  kp=-real(freq(den,num,real(s)))
+  i=find(kp>=0)
+  kp=kp(i)
+  s=s(i)
 endfunction