scilab/modules/cacsd/macros/chart.sci

   1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 // Copyright (C) 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 chart(attenu,angl,flags)
  10
  11   titre=_("Amplitude and phase contours of y/(1+y)")
  12   l10=log(10);
  13   ratio=%pi/180;
  14   defaultAngles = [1:10,20:10:160];
  15   //
  16   [lhs,rhs]=argn(0)
  17
  18   select rhs
  19
  20   case 3 then
  21
  22   case 2 then,
  23     if type(angl)==15 then
  24           // angl actually stands for flags
  25       flags=angl
  26       angl = defaultAngles;
  27     else
  28       // angl is the phase
  29       flags=[]
  30     end
  31
  32   case 1 then
  33     if type(attenu)==15 then
  34       flags=attenu
  35       attenu=[-12 -8 -6 -5 -4 -3 -2 -1.4 -1 -.5 ,..
  36               0.25 0.5 0.7 1 1.4 2 2.3 3 4 5 6 8 12];
  37     else
  38       flags=list()
  39     end
  40     angl = defaultAngles;
  41   else
  42     flags=list()
  43     attenu=[-12 -8 -6 -5 -4 -3 -2 -1.4 -1 -.5 ,..
  44             0.25 0.5 0.7 1 1.4 2 2.3 3 4 5 6 8 12];
  45     angl = defaultAngles;
  46   end
  47
  48   // convert angles to radian
  49   angl = -angl * ratio;
  50
  51   c1=color('lightgrey');c2=c1
  52   select size(flags)
  53   case 0 then
  54     flags=list(0,-1,c1,c2)
  55   case 1 then
  56     flags=list(flags(1),-1,c1,c2)
  57   case 2 then
  58     flags=list(flags(1),flags(2),c1,c2)
  59   case 3 then
  60     flags(4)=c2
  61   end
  62   //
  63
  64
  65   if ~flags(1) then
  66     clf()
  67     ax=gca();
  68     ax.data_bounds=[-360,-50;0,40];
  69     ax.axes_visible='on';
  70     ax.box='on';
  71     ax.tight_limits='on'
  72   else
  73     ax=gca();
  74   end
  75   ax.clip_state="clipgrf"
  76 //  ax.data_bounds=[0,-50;360,40];
  77
  78   phi_min=ax.data_bounds(1,1)
  79   phi_max=ax.data_bounds(2,1)
  80   k1=floor(phi_min/180)
  81   k2=ceil(phi_max/180)
  82   //
  83
  84   fig=gcf();
  85   immediate_drawing=fig.immediate_drawing;
  86   fig.immediate_drawing="off";
  87
  88   if flags(2) then xtitle(titre,_("phase(y) (degree)"),_("magnitude(y) (Db)")),end
  89   llrect=xstringl(0,0,'1')
  90
  91   //isogain curves
  92   lambda=exp(l10*attenu/20)
  93   rayon=lambda./(lambda.*lambda-ones(lambda))
  94   centre=-lambda.*rayon
  95    //
  96    if attenu<>[]
  97      for i = 1:prod(size(attenu)),
  98        att=attenu(i);
  99        if att<0 then
 100          w=%eps:0.03:%pi;
 101        else
 102          w=-%pi:0.03:0;
 103        end;
 104        n=prod(size(w))
 105        rf=centre(i)*ones(w)+rayon(i)*exp(%i*w);
 106        phi=atan(imag(rf),real(rf))/ratio; //phi is in [-180 0]
 107        module=20*log(abs(rf))/l10;
 108
 109        //use symetry and period to extend the curve on [k1*180 k2*180]
 110        p=[];m=[];
 111        S=[]
 112        for k=k1:k2-1
 113          if pmodulo(k,2)==0 then
 114            p=[p %nan k*180-phi($:-1:1)]
 115            m=[m %nan module($:-1:1)]
 116            if att>0 then
 117              xstring(p($),m($),string(att),0,0),
 118              e=gce();
 119              if ~flags(1) then e.clip_state='off';end
 120              S=[e S]
 121            end
 122          else
 123            p=[p %nan ((k+1)*180)+phi]
 124            m=[m %nan module]
 125            if att<0 then
 126              xstring(p($),m($),string(att),0,0),
 127              e=gce();
 128              if ~flags(1) then e.clip_state='off';end
 129              S=[e S]
 130            end
 131          end
 132
 133        end
 134        xpoly(p,m)
 135        e=gce();e.foreground=flags(3);//e.line_style=3,
 136        if size(S,'*')>1 then S=glue(S),end
 137        S=glue([S,e]);S.user_data=att;
 138      end;
 139      glue(ax.children(size(attenu,'*'):-1:1))
 140    end
 141   //isophase curves
 142   if angl<>[] then
 143   eps=100*%eps;
 144   for teta=angl,
 145     if teta < -%pi/2 then
 146       last=teta-eps,
 147     else
 148       last=teta+eps,
 149     end;
 150     w=[-170*ratio:0.03:last last]
 151     n=prod(size(w));
 152     module=real(20*log((sin(w)*cos(teta)/sin(teta)-cos(w)))/l10)
 153     w=w/ratio
 154     p=[];m=[];
 155     for k=k1:k2-1
 156       if pmodulo(k,2)==0 then
 157         p=[p %nan k*180-w($:-1:1)]
 158         m=[m %nan module($:-1:1)]
 159       else
 160         p=[p %nan ((k+1)*180)+w]
 161         m=[m %nan module]
 162       end
 163     end
 164     xpoly(p,m)
 165     e=gce();e.foreground=flags(3);//e.line_style=3,
 166     e.user_data=teta
 167   end;
 168   glue(ax.children(size(angl,'*'):-1:1))
 169   end
 170   fig.immediate_drawing=immediate_drawing;
 171 endfunction