bug 7566 fix: The cacsd module graphic functions (bode, black, nyquist,...)
[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   
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