bug 6394 fix + datatips custimization + revision of help pages
[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 mini(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   drawlater()
144   ax=gca();
145   if size(ax.children,"*")==0 then
146     ax.data_bounds=[xmn ymn;xmx ymx];
147     ax.x_label.text=_("Phase (deg)");
148     ax.y_label.text=_("Magnitude (dB)")
149   else
150     ax.data_bounds=[min([xmn ymn],ax.data_bounds(1,:));
151                     max([xmx ymx],ax.data_bounds(2,:))];
152   end
153   ax.axes_visible="on";
154   ax.clip_state="clipgrf";
155   r=xstringl(0,0,"m");r=r(3)
156   E=[]
157   if ks($)<size(phi,2) then last=$;else    last=$-1;end
158   for k=1:mn
159     e2=[]
160     if size(ks,"*") >1 then
161       d_phi=phi(k,ks(1:last)+1)-phi(k,ks(1:last));
162       d_d=d(k,ks(1:last)+1)-d(k,ks(1:last));
163       dd=400*sqrt((d_phi/dx).^2+(d_d/dy).^2);
164       if dd>0 then
165         // we should use xarrows or xsegs here.
166         // However their displayed arrow size depends
167         // on the data bounds and we want to avoid this.
168         xpolys([phi(k,ks(1:last));phi(k,ks(1:last))+d_phi./dd],..
169                [d(k,ks(1:last));d(k,ks(1:last))+d_d./dd]);
170         ea=gce();
171         ea.children.foreground=k;
172         ea.children.polyline_style = 4;
173         ea.children.arrow_size_factor = 1.5;
174
175         //xarrows([phi(k,ks(1:last));phi(k,ks(1:last))+d_phi./dd],..
176         //    [d(k,ks(1:last));d(k,ks(1:last))+d_d./dd],60)
177         //ea=gce();ea.segs_color=k*ones(dd);
178         //add a frequency label near each arrow
179         el=[];
180
181         for l=ks
182           xstring(phi(k,l)+r,d(k,l),msprintf("%-5.2g",frq(kf,l)))
183           e=gce()
184           e.font_foreground=k;
185           el=[e,el]
186         end
187         e2=glue([el ea])
188       end
189     end
190
191     xpoly(phi(k,:),d(k,:));e1=gce()
192     e1.foreground=k;
193     datatipInitStruct(e1,"formatfunction","formatBlackTip","freq",frq(kf,:))
194
195     // glue entities relative to a single black curve
196     E=[E glue([e2 e1])]
197     kf=kf+ilf
198   end
199
200
201   //  2.3 db curve
202   mbf=2.3;
203   lmda=exp(log(10)/20*mbf);
204   r=lmda/(lmda**2-1);
205   npts=100;
206   crcl=exp(%i*(-%pi:(2*%pi/npts):%pi));
207   lgmt=log(-r*crcl+r*lmda*ones(crcl));
208   xpoly([180*(imag(lgmt)/%pi-ones(lgmt))],[(20/log(10)*real(lgmt))])
209   e=gce();e.foreground=2;e.line_style=3;
210   if comments<>[] then
211     c=[];for k=1:mn,c=[E(k).children(1),c];end
212     legend([c e]',["2.3"+_("dB");comments(:)])
213   end
214   drawnow()
215 endfunction
216 function str=formatBlackTip(curve,pt,index)
217 //This function is called by the datatip mechanism to format the tip
218 //string for black curves.
219   ud=datatipGetStruct(curve);
220   if index<>[] then
221     f=ud.freq(index)
222   else //interpolated
223     [d,ptp,i,c]=orthProj(curve.data,pt)
224     f=ud.freq(i)+(ud.freq(i+1)-ud.freq(i))*c
225   end
226   str=msprintf("%.4g"+_("°")+"\n%.4g"+_("dB")+"\n%.4g"+_("Hz"), pt,f);
227 endfunction