bug 7566 fix: The cacsd module graphic functions (bode, black, nyquist,...)
[scilab.git] / scilab / modules / cacsd / macros / black.sci
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
9
10
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
14 //Syntax:
15 //
16 //           black( sl,fmin,fmax [,pas] [,comments] )
17 //           black(frq,db,phi [,comments])
18 //           black(frq, repf  [,comments])
19 //
20 //  sl       : SIMO linear system (see syslin). In case of multi-output
21 //             system the outputs are plotted with differents symbols.
22 //
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.
28 //
29 //  frq      : (row)-vector of frequencies (in Hz) or (SIMO case) matrix
30 //             of frequencies.
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.
34
35 //To plot the grid of iso-gain and iso-phase of y/(1+y) use abaque()
36 //%Example
37 //  s=poly(0,'s")
38 //  h=syslin("c",(s**2+2*0.9*10*s+100)/(s**2+2*0.3*10.1*s+102.01))
39 //  abaque();
40 //  black(h,0.01,100,"(s**2+2*0.9*10*s+100)/(s**2+2*0.3*10.1*s+102.01)")
41 //  //
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"])
44 //See also:
45 //  bode nyquist  abaque freq repfreq
46 //!
47   rhs=size(varargin)
48   if type(varargin($))==10 then
49     comments=varargin($),rhs=rhs-1;
50   else
51     comments=[];
52   end
53   fname="black";//for error messages
54   fmax=[]
55   if or(typeof(varargin(1))==["state-space" "rational"]) then
56     //sys,fmin,fmax [,pas] or sys,frq
57     refdim=1 //for error message
58     if rhs==1 then
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"),..
63                        fname,2,1))
64       end
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))
68     else
69       error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,1,5))
70     end
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]
75     refdim=2
76     select rhs
77     case 2 then //frq,repf
78       frq=varargin(1);
79       if size(frq,2)<2 then
80         error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"),..
81                        fname,1,1))
82       end
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"),..
85                        fname,1,2))
86       end
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"),..
92                        fname,1,2))
93       end
94       if size(frq,2)<>size(phi,2) then
95         error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
96                        fname,1,3))
97       end
98     else
99       error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,2,4))
100     end
101   else
102     error(msprintf(_("%s: Wrong type for input argument #%d: Linear dynamical system or row vector of floats expected.\n"),fname,1))
103   end;
104
105   if size(frq,1)==1 then
106     ilf=0
107   else
108     ilf=1
109   end
110
111   [mn,n]=size(phi);
112   if and(size(comments,"*")<>[0 mn]) then
113     error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same number of elements expected.\n"),...
114                    fname,refdim,rhs+1))
115   end
116
117   //
118   xmn=floor(min(phi)/90)*90
119   xmx=ceil(max(phi)/90)*90
120   ymn=min(d)
121   ymx=max(d)
122
123
124   kf=1
125   phi1=phi+5*ones(phi);
126
127   kk=1;p0=[phi(:,kk) d(:,kk)];ks=1;dst=0;
128   dx=max(%eps,xmx-xmn);
129   dy=max(%eps,ymx-ymn);
130   dx2=dx^2;dy2=dy^2
131
132   while kk<n
133     kk=kk+1
134     dst=dst+min(sqrt(((phi(:,kk-1)-phi(:,kk))^2)/dx2+((d(:,kk-1)-d(:,kk))^2)/dy2))
135     if dst>0.2 then
136       if min(abs(frq(:,ks(prod(size(ks))))-frq(:,kk))./frq(:,kk))>0.2 then
137         ks=[ks kk]
138         dst=0
139       end
140     end
141   end
142   kf=1
143   fig=gcf();
144   immediate_drawing=fig.immediate_drawing;
145   fig.immediate_drawing="off";
146
147   ax=gca();
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)")
152   else
153     ax.data_bounds=[min([xmn ymn],ax.data_bounds(1,:));
154                     max([xmx ymx],ax.data_bounds(2,:))];
155   end
156   ax.axes_visible="on";
157   ax.clip_state="clipgrf";
158   r=xstringl(0,0,"m");r=r(3)
159   E=[]
160   if ks($)<size(phi,2) then last=$;else    last=$-1;end
161   for k=1:mn
162     e2=[]
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);
167       if dd>0 then
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]);
173         ea=gce();
174         ea.children.foreground=k;
175         ea.children.polyline_style = 4;
176         ea.children.arrow_size_factor = 1.5;
177
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
182         el=[];
183
184         for l=ks
185           xstring(phi(k,l)+r,d(k,l),msprintf("%-5.2g",frq(kf,l)))
186           e=gce()
187           e.font_foreground=k;
188           el=[e,el]
189         end
190         e2=glue([el ea])
191       end
192     end
193
194     xpoly(phi(k,:),d(k,:));e1=gce()
195     e1.foreground=k;
196     datatipInitStruct(e1,"formatfunction","formatBlackTip","freq",frq(kf,:))
197
198     // glue entities relative to a single black curve
199     E=[E glue([e2 e1])]
200     kf=kf+ilf
201   end
202
203
204   //  2.3 db curve
205   mbf=2.3;
206   lmda=exp(log(10)/20*mbf);
207   r=lmda/(lmda**2-1);
208   npts=100;
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;
213   if comments<>[] then
214     c=[];for k=1:mn,c=[E(k).children(1),c];end
215     legend([c e]',["2.3"+_("dB");comments(:)])
216   end
217   fig.immediate_drawing=immediate_drawing;
218 endfunction
219
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);
224   if index<>[] then
225     f=ud.freq(index)
226   else //interpolated
227     [d,ptp,i,c]=orthProj(curve.data,pt)
228     f=ud.freq(i)+(ud.freq(i+1)-ud.freq(i))*c
229   end
230   str=msprintf("%.4g°\n%.4g"+_("dB")+"\n%.4g"+_("Hz"), pt,f);
231 endfunction