1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 1998-2010 - INRIA - Serge Steer
3 // Copyright (C) 2010 - DIGITEO - Yann COLLETTE
4 // This file must be used under the terms of the CeCILL.
5 // This source file is licensed as described in the file COPYING, which
6 // you should have received as part of this distribution. The terms
7 // are also available at
8 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11 function black(varargin)
12 //Black's diagram (Nichols chart) for a linear system sl.
13 //sl can be a continuous-time, discrete-time or sampled SIMO system
16 // black( sl,fmin,fmax [,pas] [,comments] )
17 // black(frq,db,phi [,comments])
18 // black(frq, repf [,comments])
20 // sl : SIMO linear system (see syslin). In case of multi-output
21 // system the outputs are plotted with differents symbols.
23 // fmin : minimal frequency (in Hz).
24 // fmax : maximal frequency (in Hz).
25 // pas : logarithmic discretization step. (see calfrq for the
26 // choice of default value).
27 // comments : character strings to comment the curves.
29 // frq : (row)-vector of frequencies (in Hz) or (SIMO case) matrix
31 // db : matrix of modulus (in Db). One row for each response.
32 // phi : matrix of phases (in degrees). One row for each response.
33 // repf : matrix of complex numbers. One row for each response.
35 //To plot the grid of iso-gain and iso-phase of y/(1+y) use abaque()
38 // h=syslin("c",(s**2+2*0.9*10*s+100)/(s**2+2*0.3*10.1*s+102.01))
40 // black(h,0.01,100,"(s**2+2*0.9*10*s+100)/(s**2+2*0.3*10.1*s+102.01)")
42 // h1=h*syslin("c",(s**2+2*0.1*15.1*s+228.01)/(s**2+2*0.9*15*s+225))
43 // black([h1;h],0.01,100,["h1";"h"])
45 // bode nyquist abaque freq repfreq
48 if type(varargin($))==10 then
49 comments=varargin($),rhs=rhs-1;
53 fname="black";//for error messages
55 if or(typeof(varargin(1))==["state-space" "rational"]) then
56 //sys,fmin,fmax [,pas] or sys,frq
57 refdim=1 //for error message
59 [frq,repf]=repfreq(varargin(1),1d-3,1d3)
60 elseif rhs==2 then //sys,frq
61 if size(varargin(2),2)<2 then
62 error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"),..
65 [frq,repf]=repfreq(varargin(1:rhs))
66 elseif or(rhs==(3:4)) then //sys,fmin,fmax [,pas]
67 [frq,repf]=repfreq(varargin(1:rhs))
69 error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,1,5))
71 [phi,d]=phasemag(repf)
72 if rhs>=3 then fmax=varargin(3),end
73 elseif type(varargin(1))==1 then
74 //frq,db,phi [,comments] or frq, repf [,comments]
77 case 2 then //frq,repf
80 error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"),..
83 if size(frq,2)<>size(varargin(2),2) then
84 error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
87 [phi,d]=phasemag(varargin(2))
88 case 3 then //frq,db,phi
89 [frq,d,phi]=varargin(1:rhs)
90 if size(frq,2)<>size(d,2) then
91 error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
94 if size(frq,2)<>size(phi,2) then
95 error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
99 error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,2,4))
102 error(msprintf(_("%s: Wrong type for input argument #%d: Linear dynamical system or row vector of floats expected.\n"),fname,1))
105 if size(frq,1)==1 then
112 if and(size(comments,"*")<>[0 mn]) then
113 error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same number of elements expected.\n"),...
118 xmn=floor(min(phi)/90)*90
119 xmx=ceil(max(phi)/90)*90
125 phi1=phi+5*ones(phi);
127 kk=1;p0=[phi(:,kk) d(:,kk)];ks=1;dst=0;
128 dx=max(%eps,xmx-xmn);
129 dy=max(%eps,ymx-ymn);
134 dst=dst+min(sqrt(((phi(:,kk-1)-phi(:,kk))^2)/dx2+((d(:,kk-1)-d(:,kk))^2)/dy2))
136 if min(abs(frq(:,ks(prod(size(ks))))-frq(:,kk))./frq(:,kk))>0.2 then
144 immediate_drawing=fig.immediate_drawing;
145 fig.immediate_drawing="off";
148 if size(ax.children,"*")==0 then
149 ax.data_bounds=[xmn ymn;xmx ymx];
150 ax.x_label.text=_("Phase (deg)");
151 ax.y_label.text=_("Magnitude (dB)")
153 ax.data_bounds=[min([xmn ymn],ax.data_bounds(1,:));
154 max([xmx ymx],ax.data_bounds(2,:))];
156 ax.axes_visible="on";
157 ax.clip_state="clipgrf";
158 r=xstringl(0,0,"m");r=r(3)
160 if ks($)<size(phi,2) then last=$;else last=$-1;end
163 if size(ks,"*") >1 then
164 d_phi=phi(k,ks(1:last)+1)-phi(k,ks(1:last));
165 d_d=d(k,ks(1:last)+1)-d(k,ks(1:last));
166 dd=400*sqrt((d_phi/dx).^2+(d_d/dy).^2);
168 // we should use xarrows or xsegs here.
169 // However their displayed arrow size depends
170 // on the data bounds and we want to avoid this.
171 xpolys([phi(k,ks(1:last));phi(k,ks(1:last))+d_phi./dd],..
172 [d(k,ks(1:last));d(k,ks(1:last))+d_d./dd]);
174 ea.children.foreground=k;
175 ea.children.polyline_style = 4;
176 ea.children.arrow_size_factor = 1.5;
178 //xarrows([phi(k,ks(1:last));phi(k,ks(1:last))+d_phi./dd],..
179 // [d(k,ks(1:last));d(k,ks(1:last))+d_d./dd],60)
180 //ea=gce();ea.segs_color=k*ones(dd);
181 //add a frequency label near each arrow
185 xstring(phi(k,l)+r,d(k,l),msprintf("%-5.2g",frq(kf,l)))
194 xpoly(phi(k,:),d(k,:));e1=gce()
196 datatipInitStruct(e1,"formatfunction","formatBlackTip","freq",frq(kf,:))
198 // glue entities relative to a single black curve
206 lmda=exp(log(10)/20*mbf);
209 crcl=exp(%i*(-%pi:(2*%pi/npts):%pi));
210 lgmt=log(-r*crcl+r*lmda*ones(crcl));
211 xpoly([180*(imag(lgmt)/%pi-ones(lgmt))],[(20/log(10)*real(lgmt))])
212 e=gce();e.foreground=2;e.line_style=3;
214 c=[];for k=1:mn,c=[E(k).children(1),c];end
215 legend([c e]',["2.3"+_("dB");comments(:)])
217 fig.immediate_drawing=immediate_drawing;
220 function str=formatBlackTip(curve,pt,index)
221 //This function is called by the datatip mechanism to format the tip
222 //string for black curves.
223 ud=datatipGetStruct(curve);
227 [d,ptp,i,c]=orthProj(curve.data,pt)
228 f=ud.freq(i)+(ud.freq(i+1)-ud.freq(i))*c
230 str=msprintf("%.4g°\n%.4g"+_("dB")+"\n%.4g"+_("Hz"), pt,f);