bug 7566 fix: The cacsd module graphic functions (bode, black, nyquist,...)
[scilab.git] / 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 5 10 20 30 50 90 120 150 180]
39   defaultModules=[mod_min:20:-40 -12 -6  -3   -1 0 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