Merge remote-tracking branch 'origin/6.1'
[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 // Copyright (C) 2012 - 2016 - Scilab Enterprises
4 // Copyright (C) 2018 - Samuel GOUGEON
5 //
6 // This file is hereby licensed under the terms of the GNU GPL v2.0,
7 // pursuant to article 5.3.4 of the CeCILL v.2.1.
8 // This file was originally licensed under the terms of the CeCILL v2.1,
9 // and continues to be available under such terms.
10 // For more information, see the COPYING file which you should have received
11 // along with this program.
12
13 function handles = nicholschart(modules, args, colors)
14
15     fname = "nicholschart"
16     [lhs,rhs]=argn(0);
17
18     l10=log(10);
19     d2rad=%pi/180;
20
21     fig=gcf();
22     immediate_drawing=fig.immediate_drawing;
23     fig.immediate_drawing="off";
24
25     ax=gca();
26     blankAxes = length(ax.children)==0;
27     old_data_bounds = ax.data_bounds;
28     nc=size(ax.children,"*")
29     if nc==0 then
30         ax.data_bounds=[-360,-40;0,40];
31         ax.axes_visible="on";
32         ax.box="on";
33         ax.tight_limits="on";
34         ax.title.text=_("Amplitude and phase contours of y/(1+y)");
35         ax.x_label.text=_("phase(y) (degree)");
36         ax.y_label.text=_("magnitude(y) (dB)");
37     else
38         ax.data_bounds(2,2)=max( ax.data_bounds(2,2),40);
39     end
40     ax.clip_state="clipgrf"
41
42     phi_min=ax.data_bounds(1,1)
43     phi_max=ax.data_bounds(2,1)
44     mod_min=ax.data_bounds(1,2)
45     mod_max=ax.data_bounds(2,2)
46
47     defaultArgs = [1 2 5 10 20 30 50 70 90 120 140 160 180];
48     defaultModules=[-30:-10:mod_min -25 -20 -15 -12 -9 -6 -3 -2 -1 -0.5 -0.2 -0.1 0 0.1 0.2 0.5  1  2 3 6 12];
49
50
51     if exists("modules","local") == 0 | modules == [] then
52         modules=defaultModules
53     else
54         if type(modules)<>1|~isreal(modules) then
55             msg = _("%s: Wrong type for input argument ""%s"": real floating point array expected\n")
56             error(msprintf(msg, fname, "modules"));
57         end
58         modules=matrix(modules,1,-1)
59     end
60     if exists("args","local")==0 | args == [] then
61         args=defaultArgs
62     else
63         if type(args)<>1|~isreal(args) then
64             msg = _("%s: Wrong type for input argument ""%s"": real floating point array expected\n")
65             error(msprintf(msg, fname, "args"));
66         end
67         args=matrix(args,1,-1)
68     end
69
70     // colors
71     if exists("colors","local")==0 | colors == [] then
72         colors = "grey85";
73     end
74     c = iscolor(colors);
75     if or(isnan(c))
76         msg = _("%s: Argument #%d: Wrong color specification.\n")
77         error(msprintf(msg, fname, 3));
78     end
79     if size(c,1)==1 then
80         c = [c ; c];    // Same color for both subframes
81     else
82         c = c(1:2, :);  // selects only the 2 first inputs (no warning)
83     end
84     if size(c,2)==3
85         colors = addcolor(c);
86     else
87         colors = c
88     end
89     if length(colors)<2 then
90         colors = colors * [1 1]
91     end
92
93     // convert args to radian and insure negative
94     args = -abs(args) * d2rad;
95
96     //initialize handles array for chart entities
97     chart_handles=[]
98
99     // Replication bounds
100     k1=floor(phi_min/180)
101     k2=ceil(phi_max/180)
102
103     //isogain curves: y as fixed gain and varying phase
104     //-------------------------------------------------
105     c = addcolor(min(1, gcf().color_map(colors(1),:)*0.8));  // labels are a bit darker
106     gainCurves = [];
107     gainLabels = [];
108     phaseCurves = [];
109
110     if modules<>[] then
111         w=[linspace(-%pi,-0.1,100) linspace(-0.1,0,80) ]
112         nw=size(w,"*")
113
114         // Calibrating labels dimensions
115         str = msprintf("%.2gdB",0.25)
116         r = xstringl(0,0,str)
117
118         // Main loop on gains
119         for i = 1:prod(size(modules)),
120             att=modules(i);
121             y=10^(att/20)*exp(%i*w);
122             y(y==1)=[];//remove singular point if any
123             rf=y./(1-y);
124             [module, phi]=dbphi(rf)
125             //use symetry and period to extend the curve on [k1*180 k2*180]
126             p=[];m=[];
127             S=[];cut=[]
128             for k = k1:k2-1
129                 if pmodulo(k,2)==0 then
130                     p = [p cut k*180-phi($:-1:1)]
131                     m = [m cut module($:-1:1)]
132                     if att>0 then
133                         str = msprintf("%.2gdB",att)
134                         xstring(k*180-phi($)-r(3)/6, module($)+r(4)/6, str, 0, 0),
135                         e=gce();
136                         e.font_foreground = c;
137                         e.text_box_mode = "centered";
138                         S=[e S]
139                     elseif att==0 then
140                         l = find(module>mod_max-r(4),1)
141                         if l<>[] then
142                             xstring(k*180-phi(l-1),module(l-1)+r(4)/5,"0dB",0,0),
143                             e=gce();
144                             e.font_foreground = c;
145                             e.text_box_mode = "centered";
146                             S=[e S]
147                         end
148                     end
149                 else
150                     p=[p cut ((k+1)*180)+phi]
151                     m=[m cut module]
152                     if att<0 then
153                         str = msprintf("%.2gdB",att)
154                         xstring(p($)-r(3)/2,m($)+r(4)/5,str,0,0),
155                         e = gce();
156                         e.font_foreground = c;
157                         e.text_box_mode = "centered";
158                         S=[e S]
159                     end
160                 end
161                 cut=%nan
162             end
163             gainLabels = [S gainLabels];
164
165             xpoly(p,m)
166             e=gce();
167             e.foreground=colors(1),
168             e.line_style = 1;
169             e.display_function = "formatNicholsGainTip";
170             e.display_function_data = att;
171             gainCurves = [e gainCurves];
172
173             if size(S,"*")>1 then S=glue(S),end
174             chart_handles=[glue([S,e]),chart_handles];
175         end
176     end
177
178     //isophase curves: y as fixed phase and varying gain
179     //-------------------------------------------------
180     if args<>[] then
181
182         eps=10*%eps;
183         for teta=args,
184             //w = teta produce a 0 gain and consequently a singularity in module
185             if teta < -%pi/2 then
186                 last=teta-eps,
187             else
188                 last=teta+eps,
189             end;
190             //use logarithmic discretization to have more mesh points near low modules
191             w=real(logspace(log10(-last),log10(170*d2rad),150))
192             w=-w($:-1:1)
193
194             n=prod(size(w));
195             module=real(20*log((sin(w)*cos(teta)/sin(teta)-cos(w)))/l10)
196             w=w/d2rad
197             //use symetry and period to extend the curve on [k1*180 k2*180]
198             p=[];m=[];
199             for k=k1:k2-1
200                 if pmodulo(k,2)==0 then
201                     p=[p %nan k*180-w($:-1:1)]
202                     m=[m %nan module($:-1:1)]
203                 else
204                     p=[p %nan ((k+1)*180)+w]
205                     m=[m %nan module]
206                 end
207             end
208             xpoly(p,m)
209             e=gce();
210             e.foreground=colors(2);
211             e.line_style = 1;
212             e.display_function = "formatNicholsPhaseTip";
213             e.display_function_data = teta * 180 / %pi;
214             phaseCurves = [e phaseCurves];
215             chart_handles=[e chart_handles]
216         end;
217     end
218     chart_handles=glue(chart_handles)
219     //reorder axes children to make chart drawn before the black curves if any
220     for k=1:nc
221         swap_handles(ax.children(k),ax.children(k+1))
222     end
223     fig.immediate_drawing=immediate_drawing;
224
225     // reset data_bounds
226     if rhs == 0 | blankAxes then
227         ax.data_bounds = [-360,-40;0,40];
228     else
229         ax.data_bounds = old_data_bounds;
230     end
231
232     // output
233     handles = struct("phaseFrame", phaseCurves($:-1:1), "gainFrame", gainCurves($:-1:1), "gainLabels", gainLabels($:-1:1));
234 endfunction