bug 7095 fix, the margins unit_test fails, g_margin, p_margin and unit_test improved
[scilab.git] / scilab / modules / cacsd / tests / unit_tests / margins.tst
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
8 // <-- ENGLISH IMPOSED -->
9
10 //gain margin
11 //-----------
12 function ok=check_gmargin(h,fref,gref)
13   eps=1e-10
14   [g,f]=g_margin(h)
15   if fref==0 then
16      ok=abs(f - fref)<eps
17   else
18     ok=abs(f - fref)/abs(fref)<eps
19   end
20   if gref==0 then
21      ok=ok&abs(g - gref)<eps
22   else
23     ok=ok&abs(g - gref)/abs(gref)<eps
24   end
25 endfunction
26
27 //discrete time case
28 z=poly(0,'z');
29 //the references solutions are computed using the following Maple instructions
30 //Digits:=50:
31 //z:=exp(I*w):
32 //assume(w, 'real', w > 0)
33 //h:= rational fraction in z
34 //M := solve(Im(h) = 0, w)
35 //evalf(M/(2*Pi))
36 //evalf(eval(20*log10(1/abs(h)), w = M[k]))
37
38 h = syslin(0.1,((29/625)+(2399/50000)*z)/((1131/1250)+(-181/100)*z+z^2));
39 f_ref=atan(sqrt(9003479136639)/4963519)/(0.1*2*%pi);
40 g_ref=6.242379356591869534116;
41 if ~check_gmargin(h,f_ref,g_ref) then pause,end
42
43
44
45 h=syslin(0.1,((43/500)+(-161/1000)*z+(39/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3));
46 f_ref=atan((1/215)*sqrt(1136805-8330*sqrt(5970))/(833/43+(1/215)*sqrt(5970)))/(0.1*2*%pi);
47 g_ref=12.3470513198103944037870;
48 if ~check_gmargin(h,f_ref,g_ref) then pause,end
49
50 h=syslin(0.1,((129/500)+(-483/1000)*z+(117/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3));  
51 f_ref=atan((1/215)*sqrt(1136805-8330*sqrt(5970))/(833/43+(1/215)*sqrt(5970)))/(0.1*2*%pi);
52 g_ref=2.8046262254171456578864840;
53 if ~check_gmargin(h,f_ref,g_ref) then pause,end
54
55
56 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));  
57 f_ref=1/2;
58 g_ref=6.48227782514616029706190;
59 if ~check_gmargin(h,f_ref,g_ref) then pause,end
60
61 //continuuous time case
62 s=poly(0,'s');
63 //the reference solutions are computed using  the following Maple instructions
64 //s:=I*w:
65 //assume(w, 'real', w > 0)
66 //h:= rational fraction in z
67 //M := solve(Im(h) = 0, w)
68 //evalf(M/(2*Pi))
69 //evalf(eval(20*log10(1/abs(h)), w = M[k]))
70
71 h=syslin('c',(-1+s)/(3+2*s+s^2));  
72 f_ref=0;
73 g_ref=9.5424250943932487459005580;
74 if ~check_gmargin(h,f_ref,g_ref) then pause,end
75
76
77 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)); 
78 f_ref=0.0714555258202006740373134;
79 g_ref=-4.91687406933815400242335;
80 if ~check_gmargin(h,f_ref,g_ref) then pause,end
81
82
83 h=syslin('c',485000/(10000*s+200*s^2+s^3));  
84 f_ref=100/(2*%pi);
85 g_ref=12.305765141234350772862319;
86 if ~check_gmargin(h,f_ref,g_ref) then pause,end
87
88
89 h = syslin('c',1/(s+2*s^2+s^3));
90 f_ref=1/(2*%pi);
91 g_ref=6.0205999132796239042747779;
92 if ~check_gmargin(h,f_ref,g_ref) then pause,end
93
94 //phase margin
95 //-----------
96 function ok=check_pmargin(h,f_ref,p_ref)
97   eps=1e-9
98   [p,f]=p_margin(h)
99   if f_ref==0 then
100      ok=abs(f - f_ref)<eps
101   else
102     ok=abs(f - f_ref)/abs(f_ref)<eps
103   end
104   if p_ref==0 then
105      ok=ok&abs(p - p_ref)<eps
106   else
107     ok=ok&abs(p - p_ref)/abs(p_ref)<eps
108   end
109 endfunction
110
111 //discrete time case
112 //the reference solutions are computed using  the following Maple instructions
113 //z:=exp(I*w):
114 //assume(w, 'real', w > 0)
115 //P := solve(abs(h) = 1, w)
116 //evalf(-(eval(180-180*argument(h)/Pi, w = Re(P[k]))))+360
117
118 h = syslin(0.1,((29/625)+(2399/50000)*z)/((1131/1250)+(-181/100)*z+z^2)); 
119 f_ref=0.693016600315284442350578876;
120 p_ref=13.57115563612946355428439468;
121 if  ~check_pmargin(h,f_ref,p_ref) then pause,end
122
123
124 h=syslin(0.1,((43/500)+(-161/1000)*z+(39/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3)); 
125 [p ,f]=p_margin(h);
126 if p<>%inf  then pause,end
127 if f<>[]  then pause,end
128
129 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);
130 h=syslin(0.1,((129/500)+(-483/1000)*z+(117/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3)); 
131 f_ref=0.212336488950669705771059018;
132 p_ref=52.94967415965772478856630911;
133 if  ~check_pmargin(h,f_ref,p_ref) then pause,end
134
135 //continuous case
136 //the reference solutions are computed using  the following Maple instructions
137 //z:=I*w:
138 //assume(w, 'real', w > 0)
139 //P := solve(abs(h) = 1, w)
140 //evalf(-(eval(180-180*argument(h)/Pi, w = Re(P[k]))))+360
141
142 h=syslin('c',((11/10)+(12/5)*s+(7/10)*s^2)/(3+2*s+s^2));
143 f_ref=(1/51)*sqrt(15861-204*sqrt(3562))/(2*%pi);
144 p_ref=-148.547076202317410601324666;
145 if  ~check_pmargin(h,f_ref,p_ref) then pause,end
146
147
148 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));  
149 f_ref=0.09144216563554157543991;
150 p_ref=-13.1128497150069802772313;
151 if  ~check_pmargin(h,f_ref,p_ref) then pause,end