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