a9125dcd3681580c49dc0e9d39704f58ea3b2ab2
[scilab.git] / scilab / modules / cacsd / macros / bode.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 1985 - 2016 - INRIA - Serge Steer
3 // Copyright (C) 2012 - 2016 - Scilab Enterprises
4 // Copyright (C) 2018 - Samuel GOUGEON
5 //
6 // This file is hereby licensed under the terms of the GNU GPL v2.0,
7 // pursuant to article 5.3.4 of the CeCILL v.2.1.
8 // This file was originally licensed under the terms of the CeCILL v2.1,
9 // and continues to be available under such terms.
10 // For more information, see the COPYING file which you should have received
11 // along with this program.
12
13 function bode(varargin)
14     rhs = size(varargin)
15
16     if rhs == 0 then
17         s = poly(0, "s");
18         h1 = syslin("c", (s^2+2*0.9*10*s+100)/(s^2+2*0.3*10.1*s+102.01));
19         num = 22801+4406.18*s+382.37*s^2+21.02*s^3+s^4;
20         den = 22952.25+4117.77*s+490.63*s^2+33.06*s^3+s^4;
21         h2 = syslin("c", num/den);
22
23         bode([h1; h2], 0.01, 100, ["h1"; "h2"]);
24         return;
25     end
26
27     rad = %f;
28     if type(varargin($)) == 10 then
29         if varargin($) == "rad" then
30             rad = %t;
31             rhs = rhs-1;
32             if type(varargin($-1)) == 10 then
33                 comments = varargin($-1);
34                 rhs = rhs-1;
35             else
36                 comments = [];
37             end
38         else
39             comments = varargin($);
40             rhs = rhs-1;
41         end
42     else
43         comments = [];
44     end
45     fname = "bode"; // For error messages
46     fmax = [];
47     discr = %f; // For Shannon limit
48     if or(typeof(varargin(1)) == ["state-space" "rational" "zpk"]) then
49         // sys, fmin, fmax [,pas] or sys, frq
50         refdim = 1; // for error message
51         discr = varargin(1).dt<>"c";
52         if rhs == 1 then // sys
53             [frq, repf] = repfreq(varargin(1), 1d-3, 1d3);
54         elseif rhs == 2 then // sys, frq
55             if size(varargin(2), 2) < 2 then
56                 error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"), ..
57                 fname, 2, 1))
58             end
59             [frq, repf] = repfreq(varargin(1:rhs));
60         elseif or(rhs == (3:4)) then // sys, fmin, fmax [,pas]
61             [frq, repf] = repfreq(varargin(1:rhs));
62         else
63             error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"), fname, 1, 5))
64         end
65         [phi, d] = phasemag(repf);
66         if rhs >= 3 then fmax = varargin(3); end
67     elseif type(varargin(1)) == 1 then
68         // frq, db, phi [,comments] or frq, repf [,comments]
69         refdim = 2;
70         select rhs
71         case 2 then // frq,repf
72             frq = varargin(1);
73             if size(frq, 2) < 2 then
74                 error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"), ..
75                 fname, 1, 1))
76             end
77             if size(frq, 2) <> size(varargin(2), 2) then
78                 error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"), ..
79                 fname, 1, 2))
80             end
81             [phi, d] = phasemag(varargin(2));
82         case 3 then  // frq, db, phi
83             [frq, d, phi] = varargin(1:rhs);
84             if size(frq, 2) <> size(d, 2) then
85                 error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"), ..
86                 fname, 1, 2))
87             end
88             if size(frq, 2) <> size(phi, 2) then
89                 error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"), ..
90                 fname, 1, 3))
91             end
92         else
93             error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"), fname, 2, 4))
94         end
95     else
96         ierr=execstr("%"+overloadname(varargin(1))+"_bode(varargin(:))","errcatch")
97         if ierr<>0 then
98             error(msprintf(_("%s: Wrong type for input argument #%d: Linear dynamical system or row vector of floats expected.\n"),fname,1))
99         end
100         return
101     end
102     frq = frq';
103     d = d';
104     phi = phi';
105     [n, mn] = size(d);
106
107     if comments <> [] & size(comments, "*") <> mn then
108         msg = _("%s: Incompatible input arguments #%d and #%d: Same number of elements expected.\n")
109         error(mprintf(msg, fname, refdim, rhs+1))
110     end
111
112     fig = gcf();
113     immediate_drawing = fig.immediate_drawing;
114     fig.immediate_drawing = "off";
115
116     sciCurAxes = gca();
117     axes = sciCurAxes;
118     wrect = axes.axes_bounds;
119
120     // Magnitude plot
121     // --------------
122     axes.axes_bounds = [wrect(1), wrect(2), wrect(3), wrect(4)*0.5];
123     axes.data_bounds = [min(frq), min(d); max(frq), max(d)];
124     axes.log_flags = "lnn";
125     axes.grid = color("lightgrey")*ones(1, 3);
126     axes.grid_style = [8 8];
127     axes.sub_ticks(1) = 8;
128     axes.axes_visible = "on";
129     axes.clip_state = "clipgrf";
130     if size(d, 2) > 1 & size(frq, 2) == 1 then
131         xpolys(frq(:, ones(1, mn)), d, 1:mn);
132     else
133         xpolys(frq, d, 1:mn);
134     end
135     // Set datatips info
136     e = gce();
137
138     for i=1:size(e.children, "*")
139         e.children(i).display_function = "formatBodeMagTip"
140     end
141
142     if discr & fmax <> [] & max(frq) < fmax then
143         xpoly(max(frq)*[1; 1], axes.y_ticks.locations([1 $]));
144         e = gce();
145         e.foreground = 5;
146     end
147     xtitle("", _("Frequency (Hz)"), _("Magnitude (dB)"));
148     axesM = axes;
149
150     // Phase plot
151     // ----------
152     axes = newaxes();
153     axes.axes_bounds = [wrect(1), wrect(2)+wrect(4)*0.5, wrect(3), wrect(4)*0.5];
154     axes.data_bounds = [min(frq), min(phi); max(frq), max(phi)];
155     axes.log_flags = "lnn";
156     axes.grid = color("lightgrey")*ones(1, 3);
157     axes.grid_style = [8 8];
158     axes.sub_ticks(1) = 8;
159     axes.axes_visible = "on";
160     axes.clip_state = "clipgrf";
161     if size(phi, 2) > 1 & size(frq, 2) == 1 then
162         xpolys(frq(:, ones(1, mn)), phi, 1:mn);
163     else
164         xpolys(frq, phi, 1:mn);
165     end
166     ephi = gce();
167     // Set datatips info
168     for i=1:size(ephi.children, "*")
169         ephi.children(i).display_function = "formatBodePhaseTip";
170     end
171
172     if discr & fmax <> [] & max(frq) < fmax then
173         xpoly(max(frq)*[1; 1], axes.y_ticks.locations([1 $]));
174         e = gce();
175         e.foreground = 5;
176     end
177     xtitle("", _("Frequency (Hz)"), _("Phase (degree)"));
178
179     // Create legend
180     // -------------
181     if comments <> [] then
182         c = captions(ephi.children, comments, "lower_caption");
183         c.background = get(gcf(), "background");
184     end
185
186     // Rendering
187     fig.immediate_drawing = immediate_drawing;
188
189     // Post-tuning the vertical distribution of Bode plots
190     // ---------------------------------------------------
191     // (only after rendering, otherwise the legend position is [0 0] (bug))
192     if comments <> [] then
193         LgH = (1-c.position(2))*wrect(4)*0.5; // Legend height
194         // Ah = Magnitude axes height.
195         // Ah + Ah+LgH*Ah/0.5 = wrect(4)  => Ah (1 + 1 + 2*LgH) = wrect(4)
196         Ah = wrect(4)/(2+2*LgH);
197         axes.axes_bounds([2 4]) = [wrect(2)+Ah, wrect(4)-Ah];
198         axesM.axes_bounds(4) = Ah*1.05;  // Magnitude plot
199     end
200
201     // Returns to the initial axes
202     sca(sciCurAxes);
203
204     // rad/s mode
205     // ----------
206     if rad == %t then
207         // This function modifies the Bode diagrams for a rad/s display instead of Hz.
208         // h is a hanlde of a figure containing Bode diagrams.
209         // Ref: http://forge.scilab.org/index.php/p/cpge/source/tree/HEAD/macros/bode_Hz2rad_2.sci
210         labels = [_("Phase (degree)"); _("Magnitude (dB)")];
211         pos_h = [9, 5];
212         for k=1:2
213             for i=1:size(fig.children(k).children, 1)
214                 if fig.children(k).children(i).type == "Compound"
215                     for j=1:size(fig.children(k).children(i).children, 1)
216                         fig.children(k).children(i).children(j).data(:, 1) = fig.children(k).children(i).children(j).data(:, 1)*2*%pi;
217                     end
218
219                     // h.children(k).title.text = h.children(k).y_label.text;
220                     xmin1 = fig.children(k).data_bounds(1)*2*%pi;
221                     xmax1 = fig.children(k).data_bounds(2)*2*%pi;
222                     ymin1 = fig.children(k).data_bounds(3);
223                     ymax1 = fig.children(k).data_bounds(4);
224
225                     rect = [xmin1, ymin1, xmax1, ymax1];
226                     nb_dec = log(xmax1/xmin1)/log(10);
227                     fig.children(k).x_label.text = _("Frequency (rad/s)");
228                     fig.children(k).x_location = "bottom";
229                     fig.children(k).y_label.text = labels(k);
230                     replot(rect, fig.children(k));
231                 end
232             end
233         end
234     end
235 endfunction