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