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