* Bug #12793 fixed (proposal 2) - Cacsd: improve the bode() plots
[scilab.git] / scilab / modules / cacsd / macros / bode_asymp.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
3 // Copyright (C) 09/2013 - A. Khorshidi
4 //
5 // This file must be used under the terms of the CeCILL.
6 // This source file is licensed as described in the file COPYING, which
7 // you should have received as part of this distribution.  The terms
8 // are also available at
9 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10
11 function [] = bode_asymp(sl, w_min, w_max)
12     // Calling sequence:
13     //     syntax: bode_asymp(sl, [w_min, w_max])
14     // Arguments:
15     //     sl: a linear system given by its transfer function representation or state-space form and defined by syslin.
16     //     w_min, w_max: lower and upper frequency bounds (in rad/s)
17     // Description:
18     //     This function plots asymptotes for the magnitude and phase graphs.
19     // Examples:
20     //     s = %s; h = syslin("c", 1/(s+1)); scf(10001); bode(h); bode_asymp(h)
21     //     s = %s; h = syslin("c", 1/(s+1)); scf(10002); gainplot(h); bode_asymp(h)
22     // Ref:
23     //     http://forge.scilab.org/index.php/p/cpge/source/tree/HEAD/macros/REP_FREQ_pre_simulate.sci
24
25     rhs = argn(2);
26
27     if and(typeof(sl) <> ["state-space" "rational"]) then
28         msg = _("Wrong type for argument #%d: Rational or State-space matrix expected.\n");
29         error(msprintf(msg, 1))
30         return;
31     end
32
33     typeSL = typeof(sl);
34
35     for elem=1:size(sl, "r") // Individually draw each asymptote of "sl" elements (row vector).
36         if typeSL == "rational" then
37             h = sl(elem, 1);
38         else
39             h = clean(ss2tf(sl(elem, 1)), 1e-8);
40             // Also removing all the coefficients smaller than < 1e-8
41         end
42
43         root_num = roots(h.num);
44         root_den = roots(h.den);
45
46         if (find(root_num==0))
47             disp("Problem class of system is negative")
48         end
49         rac_nul = find(root_den==0);
50         alpha = length(rac_nul);
51         var = "s";
52         domain = h.dt;
53         if type(domain) == 1 then
54             var = "z";
55         elseif domain == "c" then
56             var = "s";
57         elseif domain == "d" then
58             var = "z";
59         end
60         if domain == [] then
61             var = "s";
62             if type(h.D) == 2 then
63                 var = varn(h.D);
64             end
65         end
66         s = poly(0, var);
67         msg = _("%s: Problem evaluating the first argument.\nn")
68         try K = horner(h*s^alpha, 0); catch error(msprintf(msg, "bode_asymp")); end
69
70         root_den(rac_nul) = []; // Removing the zeros
71
72         if length([find(abs(root_num)>1e8) find(abs(root_num)<1e-8)]) > 0 then
73             disp("Extreme root removed: "+string(root_num([find(abs(root_num)>1e8) find(abs(root_num)<1e-8)])'))
74             root_num([find(abs(root_num)>1e8) find(abs(root_num)<1e-8)]) = [];
75         end
76         if length([find(abs(root_den)>1e8) find(abs(root_den)<1e-8)]) > 0 then
77             disp("Extreme root removed: "+string(root_den([find(abs(root_den)>1e8) find(abs(root_den)<1e-8)])'))
78             root_den([find(abs(root_den)>1e8) find(abs(root_den)<1e-8)]) = [];
79         end
80
81         i = 1;
82         puls = [];
83         order = [];
84         while i<=length(root_num) // Real root
85             if isreal(root_num(i), 0) then
86                 puls = [puls -root_num(i)];
87                 order = [order 1];
88                 i = i+1;
89             else // Complex root
90                 xi = 1/sqrt(1+(imag(root_num(i))/real(root_num(i)))^2);
91                 puls = [puls -real(root_num(i))/xi];
92                 i = i+2;
93                 order = [order 2];
94             end
95         end
96         i = 1;
97         while i<=length(root_den) // Real root
98             if isreal(root_den(i), 0) then
99                 puls = [puls -root_den(i)];
100                 order = [order -1];
101                 i = i+1;
102             else // Complex root
103                 xi = 1/sqrt(1+(imag(root_den(i))/real(root_den(i)))^2);
104                 puls = [puls -real(root_den(i))/xi];
105                 order = [order -2];
106                 i = i+2;
107             end
108         end
109
110         [puls, ind] = gsort(puls, "g", "i");
111         order = order(ind);
112
113         asym = [-20*alpha];
114         phas = [-90*alpha];
115         i = 1;
116         while i<=length(puls)
117             new_dir = asym($)+order(i)*20;
118             asym = [asym new_dir];
119             new_arg = phas($)+order(i)*90;
120             phas = [phas new_arg];
121             i = i+1;
122         end
123
124         // bode(h, w_min, w_max)
125         fig = gcf();
126         try sca(fig.children(1)); catch
127             msg = _("ERROR: %s: bode() or gainplot() must be called before bode_asymp() and graphic window must still be running.\n");
128             error(msprintf(msg, "bode_asymp"))
129         end
130
131         if length(fig.children) == 1 then
132             flag = 1;
133         else
134             flag = 0;
135         end
136
137         if rhs == 1 then
138             wmin = fig.children(1).data_bounds(1, 1); // Minimal frequence, w_min
139             wmax = fig.children(1).data_bounds(2, 1); // Maximal frequence, w_max
140         else
141             wmin = w_min;
142             wmax = w_max;
143         end
144         puls_to_plot = [];
145
146         for p=real(puls)
147             if p >= wmin & p <= wmax then
148                 puls_to_plot($+1) = p;
149             end
150         end
151         puls = [wmin puls_to_plot' wmax];
152         // End change DV
153
154         eq_asymp = [20*log10(K/wmin^alpha)];
155         puls_p = [];
156         phas($+1) = phas($);
157         eq_phas = [phas(1)];
158         i = 2;
159         while i<=length(puls)
160             eq_asymp = [eq_asymp eq_asymp($)+asym(i-1)*log10(puls(i)/puls(i-1))];
161             puls_p = [puls_p puls(i-1) puls(i)];
162             eq_phas = [eq_phas phas(i-1) phas(i)];
163             i = i+1;
164         end
165
166         // Adapt abscissae to current unit (Hertz or rad/s)
167         if fig.children(1).x_label.text == "Frequency (Hz)" then
168             puls   = puls./(2*%pi);
169             puls_p = puls_p./(2*%pi);
170         end
171
172         // Change the color of asymptote:
173         color_asymp  = ["b--", "g--", "c--", "m--", "r--"];
174         rep = modulo(elem, 5); // This line will be useful if the number of systems is more than five.
175         if rep == 0 then
176             rep = 1;
177         end
178
179         if flag == 0 then
180             plot(puls_p, eq_phas(1:$-1), color_asymp($-rep));
181             sca(fig.children(2));
182             plot(puls, eq_asymp, color_asymp(rep));
183         elseif flag == 1 then
184             plot(puls, eq_asymp, color_asymp(rep));
185             fig.children.title.text = "Gainplot";
186         end
187     end
188
189 endfunction