1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) 2010 - INRIA - Serge Steer
4 //
5 //  This file is distributed under the same license as the Scilab package.
6 // =============================================================================
7 // <-- ENGLISH IMPOSED -->
8 //gain margin
9 //-----------
10 function ok=check_gmargin(h,fref,gref)
11   eps=1e-10
12   [g,f]=g_margin(h)
13   if fref==0 then
14      ok=abs(f - fref)<eps
15   else
16     ok=abs(f - fref)/abs(fref)<eps
17   end
18   if gref==0 then
19      ok=ok&abs(g - gref)<eps
20   else
21     ok=ok&abs(g - gref)/abs(gref)<eps
22   end
23 endfunction
24 //discrete time case
25 z=poly(0,'z');
26 //the references solutions are computed using the following Maple instructions
27 //Digits:=50:
28 //z:=exp(I*w):
29 //assume(w, 'real', w > 0)
30 //h:= rational fraction in z
31 //M := solve(Im(h) = 0, w)
32 //evalf(M/(2*Pi))
33 //evalf(eval(20*log10(1/abs(h)), w = M[k]))
34 h = syslin(0.1,((29/625)+(2399/50000)*z)/((1131/1250)+(-181/100)*z+z^2));
35 f_ref=atan(sqrt(9003479136639)/4963519)/(0.1*2*%pi);
36 g_ref=6.242379356591869534116;
37 if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
38 h=syslin(0.1,((43/500)+(-161/1000)*z+(39/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3));
39 f_ref=atan((1/215)*sqrt(1136805-8330*sqrt(5970))/(833/43+(1/215)*sqrt(5970)))/(0.1*2*%pi);
40 g_ref=12.3470513198103944037870;
41 if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
42 h=syslin(0.1,((129/500)+(-483/1000)*z+(117/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3));
43 f_ref=atan((1/215)*sqrt(1136805-8330*sqrt(5970))/(833/43+(1/215)*sqrt(5970)))/(0.1*2*%pi);
44 g_ref=2.8046262254171456578864840;
45 if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
46 h=syslin(1,((21/500)+(-3933/100000)*z+(-15407/100000)*z^2+(9259/50000)*z^3+(6939/(10^12))*z^4)/((-21/500)+(11/125)*z+(57/500)*z^2+(-9/25)*z^3+(1/5)*z^4));
47 f_ref=1/2;
48 g_ref=6.48227782514616029706190;
49 if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
50 //continuuous time case
51 s=poly(0,'s');
52 //the reference solutions are computed using  the following Maple instructions
53 //s:=I*w:
54 //assume(w, 'real', w > 0)
55 //h:= rational fraction in z
56 //M := solve(Im(h) = 0, w)
57 //evalf(M/(2*Pi))
58 //evalf(eval(20*log10(1/abs(h)), w = M[k]))
59 h=syslin('c',(-1+s)/(3+2*s+s^2));
60 f_ref=0;
61 g_ref=9.5424250943932487459005580;
62 if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
63 h = syslin('c',((1/5)+(4/5)*s+(0/1)*s^2+(3/10)*s^3)/((409/10000)+(1827/10000)*s+(5129/4000)*s^2+(31909/10000)*s^3+(64/25)*s^4 +s^5));
64 f_ref=0.0714555258202006740373134;
65 g_ref=-4.91687406933815400242335;
66 if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
67 h=syslin('c',485000/(10000*s+200*s^2+s^3));
68 f_ref=100/(2*%pi);
69 g_ref=12.305765141234350772862319;
70 if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
71 h = syslin('c',1/(s+2*s^2+s^3));
72 f_ref=1/(2*%pi);
73 g_ref=6.0205999132796239042747779;
74 if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
75 //phase margin
76 //-----------
77 function ok=check_pmargin(h,f_ref,p_ref)
78   eps=1e-9
79   [p,f]=p_margin(h)
80   if f_ref==0 then
81      ok=abs(f - f_ref)<eps
82   else
83     ok=abs(f - f_ref)/abs(f_ref)<eps
84   end
85   if p_ref==0 then
86      ok=ok&abs(p - p_ref)<eps
87   else
88     ok=ok&abs(p - p_ref)/abs(p_ref)<eps
89   end
90 endfunction
91 //discrete time case
92 //the reference solutions are computed using  the following Maple instructions
93 //z:=exp(I*w):
94 //assume(w, 'real', w > 0)
95 //P := solve(abs(h) = 1, w)
96 //evalf(-(eval(180-180*argument(h)/Pi, w = Re(P[k]))))+360
97 h = syslin(0.1,((29/625)+(2399/50000)*z)/((1131/1250)+(-181/100)*z+z^2));
98 f_ref=0.693016600315284442350578876;
99 p_ref=13.57115563612946355428439468;
100 if  ~check_pmargin(h,f_ref,p_ref) then bugmes();quit;end
101 h=syslin(0.1,((43/500)+(-161/1000)*z+(39/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3));
102 [p ,f]=p_margin(h);
103 if p<>%inf  then bugmes();quit;end
104 if f<>[]  then bugmes();quit;end
105 h=syslin(0.1,3*(0.086-0.161*%z+0.078*%z^2),-0.67+2.286*%z-2.61*%z^2+%z^3);
106 h=syslin(0.1,((129/500)+(-483/1000)*z+(117/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3));
107 f_ref=0.212336488950669705771059018;
108 p_ref=52.94967415965772478856630911;
109 if  ~check_pmargin(h,f_ref,p_ref) then bugmes();quit;end
110 //continuous case
111 //the reference solutions are computed using  the following Maple instructions
112 //z:=I*w:
113 //assume(w, 'real', w > 0)
114 //P := solve(abs(h) = 1, w)
115 //evalf(-(eval(180-180*argument(h)/Pi, w = Re(P[k]))))+360
116 h=syslin('c',((11/10)+(12/5)*s+(7/10)*s^2)/(3+2*s+s^2));
117 f_ref=(1/51)*sqrt(15861-204*sqrt(3562))/(2*%pi);
118 p_ref=-148.547076202317410601324666;
119 if  ~check_pmargin(h,f_ref,p_ref) then bugmes();quit;end
120 h = syslin('c',((1/5)+(4/5)*s+(3/10)*s^3)/((409/10000)+(1827/10000)*s+(5129/4000)*s^2+(31909/10000)*s^3+(64/25)*s^4+s^5));
121 f_ref=0.09144216563554157543991;
122 p_ref=-13.1128497150069802772313;
123 if  ~check_pmargin(h,f_ref,p_ref) then bugmes();quit;end