bug 6744 fix 03/303/2
Serge Steer [Wed, 28 Apr 2010 13:21:23 +0000 (15:21 +0200)]
Change-Id: I403539682d39f1ad88efee9962802a8d77141790

scilab/CHANGES_5.3.X
scilab/modules/cacsd/macros/p_margin.sci
scilab/modules/cacsd/tests/nonreg_tests/bug_6744.dia.ref [new file with mode: 0644]
scilab/modules/cacsd/tests/nonreg_tests/bug_6744.tst [new file with mode: 0644]

index 86fa87f..d10368f 100644 (file)
@@ -243,6 +243,8 @@ Bug fixes:
 * bug 6740 fixed - It was not possible to launch Scilab as Minimized 
                    or Maximized Window.
 
+* bug 6744 fixed - p_margin() returned an erroneous result.
+
 * bug 6784 fixed - It was not possible to move a Scilab installation without
                    breaking the (previously installed) ATOMS packages load.
 
index 3b706fb..2238abf 100644 (file)
@@ -1,52 +1,67 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) INRIA -  Author: Serge Steer
-// 
+// 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    
+// are also available at
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
 function [phm,fr]=p_margin(h)
   //compute the phase margin of a SISO transfer function
   select typeof(h)
-  case 'rational' then 
-  case 'state-space' then 
-    h=ss2tf(h)
-  else 
+  case 'rational' then
+  case 'state-space' then
+    h=ss2tf(h);
+  else
     error(97,1)
   end
-  if or(size(h)<>[1 1]) then 
+  if or(size(h)<>[1 1]) then
    error(msprintf(_("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"p_margin",1))
   end
-  eps=1.e-7// threshold used for testing if complex numbers are real or pure imaginary
-  
+  eps=1.e-7;// threshold used for testing if complex numbers are real or pure imaginary
+
   if h.dt=='c' then  //continuous time case
     w=poly(0,'w');
-    niw=horner(h.num,%i*w);diw=horner(h.den,%i*w);
+    niw=horner(h.num,%i*w);
+    diw=horner(h.den,%i*w);
     // |n(iw)/d(iw)|=1 <-- (n(iw)*n(-iw))/(d(iw)*d(-iw))=1 <--  (n(iw)*n(-iw)) - (d(iw)*d(-iw))=0
-    w=roots(real(niw*conj(niw)-diw*conj(diw)))
+    w=roots(real(niw*conj(niw)-diw*conj(diw)));
     //select positive real roots
-    ws=w(find((abs(imag(w))<eps)&(real(w)>0))) //frequency points with unitary modulus
-    if ws==[] then phm=[],fr=[],return,end
+    ws=w(find((abs(imag(w))<eps)&(real(w)>0))); //frequency points with unitary modulus
+    if ws==[] then 
+      phm=[];
+      fr=[];
+      return
+    end
+    ws=min(real(ws)); //keep only the first point where the open-loop phase first crosses 180°
     f=horner(h,%i*ws);
   else  //discrete time case
-    if h.dt=='d' then dt=1,else dt=h.dt,end
+    if h.dt=='d' then 
+      dt=1;
+    else 
+      dt=h.dt;
+    end
     // |h(e^(i*w*dt))|=1 <-- h(e^(i*w*dt))*h(e^(-i*w*dt))
     z=poly(0,varn(h.den));
-    sm=simp_mode();simp_mode(%f);hh=h*horner(h,1/z)-1;simp_mode(sm)
+    sm=simp_mode();
+    simp_mode(%f);
+    hh=h*horner(h,1/z)-1;
+    simp_mode(sm);
     //find the numerator roots
     z=roots(hh.num);
     z(abs(abs(z)-1)>eps)=[];// retain only roots with modulus equal to 1
     w=log(z)/(%i*dt);
     ws=real(w(abs(imag(w))<eps&real(w)>0)); //frequency points with unitary modulus
-    if ws==[] then phm=%inf,fr=[],return,end
+    if ws==[] then 
+      phm=%inf;
+      fr=[];
+      return
+    end
+    ws=min(real(ws)); //keep only the first point where the open-loop phase first crosses 180°
     f=horner(h,exp(%i*ws*dt));
   end
-  phm=atan(imag(f),real(f))// phase of the frequency response in [-%pi %pi]
-  phm=pmodulo(phm,2*%pi)-%pi
-  phm=phm*180/%pi; //transform in degree
-  [phm,k]=min(phm);ws=ws(k) //select the minimum
-  //disp([phm,ws])
-  fr=real(ws)/(2*%pi) 
+  phm=atand(imag(f),real(f));// phase of the frequency response in [-180 180]
+  phm=pmodulo(phm,360)-180;
+  fr=real(ws)/(2*%pi);
 endfunction
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_6744.dia.ref b/scilab/modules/cacsd/tests/nonreg_tests/bug_6744.dia.ref
new file mode 100644 (file)
index 0000000..a07c4f4
--- /dev/null
@@ -0,0 +1,23 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - INRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- Non-regression test for bug 6744 -->
+// <-- JVM NOT MANDATORY --> 
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=6744
+//
+// <-- Short Description -->
+//  p_margin() returns an erroneous result
+z=poly(0,'z');
+num=(0.0373327-0.0405775*z-0.1413155*z^2+0.1783403*z^3);
+den=(-0.0373327+0.0405775*z+0.1413155*z^2-0.1783403*z^3+0.04867*z^4-0.04007*z^5-0.17482*z^6+0.2*z^7);
+H=syslin(1.0,num/den);
+[m,fr]=p_margin(H);
+r=repfreq(H,fr);
+//check if fr give a response with modulus=1
+if abs(abs(r)-1)>1d-10 then bugmes();quit;end
+if abs(m-58.03)>0.01 then bugmes();quit;end
+if abs(m-pmodulo(atand(imag(r),real(r))-180,360))>1d-10 then bugmes();quit;end
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_6744.tst b/scilab/modules/cacsd/tests/nonreg_tests/bug_6744.tst
new file mode 100644 (file)
index 0000000..6e9a160
--- /dev/null
@@ -0,0 +1,26 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - INRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- Non-regression test for bug 6744 -->
+// <-- JVM NOT MANDATORY -->
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=6744
+//
+// <-- Short Description -->
+//  p_margin() returns an erroneous result
+z=poly(0,'z');
+num=(0.0373327-0.0405775*z-0.1413155*z^2+0.1783403*z^3);
+den=(-0.0373327+0.0405775*z+0.1413155*z^2-0.1783403*z^3+0.04867*z^4-0.04007*z^5-0.17482*z^6+0.2*z^7);
+H=syslin(1.0,num/den);
+[m,fr]=p_margin(H);
+
+
+r=repfreq(H,fr);
+//check if fr give a response with modulus=1
+if abs(abs(r)-1)>1d-10 then pause,end
+if abs(m-58.03)>0.01 then pause,end
+if abs(m-pmodulo(atand(imag(r),real(r))-180,360))>1d-10 then pause,end