7100001c8872d8e755d2b5d9d71947bc3b2ebaf6
[scilab.git] / 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   drawlater()
84   if flags(2) then xtitle(titre,_("phase(y) (degree)"),_("magnitude(y) (Db)")),end
85   llrect=xstringl(0,0,'1')
86   
87   //isogain curves
88   lambda=exp(l10*attenu/20)
89   rayon=lambda./(lambda.*lambda-ones(lambda))
90   centre=-lambda.*rayon
91    //
92    if attenu<>[]
93      for i = 1:prod(size(attenu)),
94        att=attenu(i);
95        if att<0 then 
96          w=%eps:0.03:%pi;
97        else 
98          w=-%pi:0.03:0;
99        end;
100        n=prod(size(w))
101        rf=centre(i)*ones(w)+rayon(i)*exp(%i*w);
102        phi=atan(imag(rf),real(rf))/ratio; //phi is in [-180 0]
103        module=20*log(abs(rf))/l10;
104        
105        //use symetry and period to extend the curve on [k1*180 k2*180]
106        p=[];m=[];
107        S=[]
108        for k=k1:k2-1
109          if pmodulo(k,2)==0 then
110            p=[p %nan k*180-phi($:-1:1)]
111            m=[m %nan module($:-1:1)]
112            if att>0 then 
113              xstring(p($),m($),string(att),0,0),
114              e=gce();
115              if ~flags(1) then e.clip_state='off';end
116              S=[e S]
117            end
118          else
119            p=[p %nan ((k+1)*180)+phi]
120            m=[m %nan module]
121            if att<0 then 
122              xstring(p($),m($),string(att),0,0),
123              e=gce();
124              if ~flags(1) then e.clip_state='off';end
125              S=[e S]
126            end
127          end
128          
129        end
130        xpoly(p,m)
131        e=gce();e.foreground=flags(3);//e.line_style=3,
132        if size(S,'*')>1 then S=glue(S),end
133        S=glue([S,e]);S.user_data=att;
134      end;
135      glue(ax.children(size(attenu,'*'):-1:1))
136    end
137   //isophase curves
138   if angl<>[] then
139   eps=100*%eps;
140   for teta=angl,
141     if teta < -%pi/2 then
142       last=teta-eps,
143     else
144       last=teta+eps,
145     end;
146     w=[-170*ratio:0.03:last last]
147     n=prod(size(w));
148     module=real(20*log((sin(w)*cos(teta)/sin(teta)-cos(w)))/l10)
149     w=w/ratio
150     p=[];m=[];
151     for k=k1:k2-1
152       if pmodulo(k,2)==0 then
153         p=[p %nan k*180-w($:-1:1)]
154         m=[m %nan module($:-1:1)]
155       else
156         p=[p %nan ((k+1)*180)+w]
157         m=[m %nan module]
158       end
159     end
160     xpoly(p,m)
161     e=gce();e.foreground=flags(3);//e.line_style=3,
162     e.user_data=teta
163   end;
164   glue(ax.children(size(angl,'*'):-1:1))
165   end
166   drawnow() ;
167 endfunction