scilab/modules/cacsd/macros/nicholschart.sci

   1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 // Copyright (C) 2010 - INRIA - Serge STEER
   3 // This file must be used under the terms of the CeCILL.
   4 // This source file is licensed as described in the file COPYING, which
   5 // you should have received as part of this distribution.  The terms
   6 // are also available at
   7 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
   8
   9 function nicholschart(modules,args,colors)
  10
  11   l10=log(10);
  12   ratio=%pi/180;
  13
  14   drawlater()
  15   ax=gca();
  16   nc=size(ax.children,"*")
  17   if nc==0 then
  18     ax.data_bounds=[-360,-40;0,40];
  19     ax.axes_visible="on";
  20     ax.box="on";
  21     ax.tight_limits="on"
  22     ax.title.text=_("Amplitude and phase contours of y/(1+y)")
  23     ax.x_label.text=_("phase(y) (degree)");
  24     ax.y_label.text=_("magnitude(y) (dB)");
  25   else
  26     ax.data_bounds(2,2)=max( ax.data_bounds(2,2),40)
  27   end
  28   ax.clip_state="clipgrf"
  29
  30   phi_min=ax.data_bounds(1,1)
  31   phi_max=ax.data_bounds(2,1)
  32   mod_min=ax.data_bounds(1,2)
  33   mod_max=ax.data_bounds(2,2)
  34
  35   defaultArgs = [1 5 10 20 30 50 90 120 150 180]
  36   defaultModules=[mod_min:20:-40 -12 -6  -3   -1 0 0.25 0.5  1  2.3 4  6  12];
  37
  38
  39   if exists("modules","local")==0 then
  40     modules=defaultModules
  41   else
  42     if type(modules)<>1|~isreal(modules) then
  43       error(msprintf("%s: Wrong type for imput argument ""%s"": real floating point array expected\n"),"nicholschart","modules");
  44     end
  45     modules=matrix(modules,1,-1)
  46   end
  47   if exists("args","local")==0 then
  48     args=defaultArgs
  49   else
  50     if type(args)<>1|~isreal(args) then
  51       error(msprintf("%s: Wrong type for imput argument ""%s"": real floating point array expected\n"),"nicholschart","args");
  52     end
  53     args=matrix(args,1,-1)
  54   end
  55   //
  56   if exists("colors","local")==0 then
  57     colors=[4 12];
  58   else
  59     if type(colors)<>1|~isreal(colors) then
  60       error(msprintf("%s: Wrong type for imput argument ""%s"": real floating point array expected\n"),"hallchart","colors");
  61     end
  62     if size(colors,"*")==1 then
  63       colors=colors*ones(1,2)
  64     end
  65   end
  66   // convert args to radian and insure negative
  67   args = -abs(args) * ratio;
  68
  69   //initialize handles array for chart entities
  70   chart_handles=[]
  71
  72   // Replication bounds
  73   k1=floor(phi_min/180)
  74   k2=ceil(phi_max/180)
  75
  76   //isogain curves: y as fixed gain and varying phase
  77   //-------------------------------------------------
  78   if modules<>[] then
  79     w=[linspace(-%pi,-0.1,100) linspace(-0.1,0,80) ]
  80     nw=size(w,"*")
  81     for i = 1:prod(size(modules)),
  82       att=modules(i);
  83       y=10^(att/20)*exp(%i*w);
  84       y(y==1)=[];//remove singular point if any
  85       rf=y./(1-y);
  86       [module, phi]=dbphi(rf)
  87       //use symetry and period to extend the curve on [k1*180 k2*180]
  88       p=[];m=[];
  89       S=[];cut=[]
  90       for k=k1:k2-1
  91         if pmodulo(k,2)==0 then
  92           p=[p cut k*180-phi($:-1:1)]
  93           m=[m cut module($:-1:1)]
  94           if att>0 then
  95             str=msprintf("%.2gdB",att)
  96             r=xstringl(0,0,str)
  97             xstring(k*180-phi($)-r(3)/2,module($),str,0,0),
  98             e=gce();
  99             e.font_foreground=colors(1)
 100             S=[e S]
 101           elseif att==0 then
 102             l=find(module>mod_max-r(4),1)
 103             if l<>[] then
 104               xstring(k*180-phi(l-1),module(l-1),"0dB",0,0),
 105               e=gce();
 106               e.font_foreground=colors(1)
 107               S=[e S]
 108             end
 109           end
 110         else
 111           p=[p cut ((k+1)*180)+phi]
 112           m=[m cut module]
 113           if att<0 then
 114             str=msprintf("%.2gdB",att)
 115             r=xstringl(0,0,str)
 116             xstring(p($)-r(3),m($),str,0,0),
 117             e=gce();
 118             e.font_foreground=colors(1)
 119             S=[e S]
 120           end
 121         end
 122         cut=%nan
 123       end
 124       xpoly(p,m)
 125       e=gce();
 126       e.foreground=colors(1),
 127       e.line_style=7;
 128       datatipInitStruct(e,"formatfunction","formatNicholsGainTip","gain",att)
 129
 130       if size(S,"*")>1 then S=glue(S),end
 131       chart_handles=[glue([S,e]),chart_handles];
 132     end;
 133   end
 134
 135   //isophase curves: y as fixed phase and varying gain
 136   //-------------------------------------------------
 137
 138   if args<>[] then
 139
 140     eps=10*%eps;
 141     for teta=args,
 142       //w = teta produce a 0 gain and consequently a singularity in module
 143       if teta < -%pi/2 then
 144         last=teta-eps,
 145       else
 146         last=teta+eps,
 147       end;
 148       //use logarithmic discretization to have more mesh points near low modules
 149       w=real(logspace(log10(-last),log10(170*ratio),150))
 150       w=-w($:-1:1)
 151
 152       n=prod(size(w));
 153       module=real(20*log((sin(w)*cos(teta)/sin(teta)-cos(w)))/l10)
 154       w=w/ratio
 155       //use symetry and period to extend the curve on [k1*180 k2*180]
 156       p=[];m=[];
 157       for k=k1:k2-1
 158         if pmodulo(k,2)==0 then
 159           p=[p %nan k*180-w($:-1:1)]
 160           m=[m %nan module($:-1:1)]
 161         else
 162           p=[p %nan ((k+1)*180)+w]
 163           m=[m %nan module]
 164         end
 165       end
 166       xpoly(p,m)
 167       e=gce();
 168       e.foreground=colors(2);
 169       e.line_style=7;
 170       datatipInitStruct(e,"formatfunction", ...
 171                         "formatNicholsPhaseTip","phase",teta*180/%pi)
 172       chart_handles=[e chart_handles]
 173     end;
 174   end
 175   chart_handles=glue(chart_handles)
 176   //reorder axes children to make chart drawn before the black curves if any
 177   for k=1:nc
 178     swap_handles(ax.children(k),ax.children(k+1))
 179   end
 180
 181   drawnow() ;
 182 endfunction
 183 function str=formatNicholsGainTip(curve,pt,index)
 184 //This function is called by the datatip mechanism to format the tip
 185 //string for the Nichols chart iso gain curves.
 186   ud=datatipGetStruct(curve);
 187   str=msprintf("%.2g"+_("dB"),ud.gain);
 188 endfunction
 189 function str=formatNicholsPhaseTip(curve,pt,index)
 190 //This function is called by the datatip mechanism to format the tip
 191 //string for the Nichols chart iso phase curves.
 192   ud=datatipGetStruct(curve);
 193   str=msprintf("%.2g"+_("°"),ud.phase)
 194 endfunction