Typo fixes
[scilab.git] / scilab / modules / cacsd / macros / nyquist.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 1984-2011 - INRIA - Serge STEER
3 // Copyright (C) 2012 - 2016 - Scilab Enterprises
4 //
5 // This file is hereby licensed under the terms of the GNU GPL v2.0,
6 // pursuant to article 5.3.4 of the CeCILL v.2.1.
7 // This file was originally licensed under the terms of the CeCILL v2.1,
8 // and continues to be available under such terms.
9 // For more information, see the COPYING file which you should have received
10 // along with this program.
11
12 function nyquist(varargin)
13     // Nyquist plot
14     //!
15     rhs=size(varargin);
16
17     if rhs == 0 then
18         //Hall chart as a grid for nyquist
19         s=poly(0,"s");
20         Plant=syslin("c",16000/((s+1)*(s+10)*(s+100)));
21         //two degree of freedom PID
22         tau=0.2;xsi=1.2;
23         PID=syslin("c",(1/(2*xsi*tau*s))*(1+2*xsi*tau*s+tau^2*s^2));
24         nyquist([Plant;Plant*PID],0.5,100,["Plant";"Plant and PID corrector"]);
25         hallchart(colors=color("light gray")*[1 1])
26         //move the caption in the lower right corner
27         ax=gca();Leg=ax.children(1);
28         Leg.legend_location="in_upper_left";
29         return;
30     end
31
32     symmetry=%t
33     if type(varargin(rhs))==4 then //symmetrization flag
34         symmetry=varargin(rhs)
35         rhs=rhs-1
36     end
37     if type(varargin(rhs))==10 then
38         comments=varargin(rhs);
39         rhs=rhs-1;
40     else
41         comments=[];
42     end
43     fname="nyquist";//for error messages
44     fmax=[];
45     if or(typeof(varargin(1))==["state-space" "rational"]) then
46         //sys,fmin,fmax [,pas] or sys,frq
47         refdim=1; //for error message
48         sltyp=varargin(1).dt;
49         if rhs==1 then
50             [frq,repf,splitf]=repfreq(varargin(1),1d-3,1d3);
51         elseif rhs==2 then //sys,frq
52             if size(varargin(2),2)<2 then
53                 error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"),fname,2,1))
54             end
55             [frq,repf]=repfreq(varargin(1:rhs));
56         elseif or(rhs==(3:4)) then //sys,fmin,fmax [,pas]
57             [frq,repf,splitf]=repfreq(varargin(1:rhs));
58         else
59             error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,1,5))
60         end
61     elseif  type(varargin(1))==1 then
62         //frq,db,phi [,comments] or frq, repf [,comments]
63         refdim=2;
64         sltyp="x";
65         splitf=[];
66         splitf=1;
67         select rhs
68         case 2 then //frq,repf
69             frq=varargin(1);
70             repf=varargin(2);
71             if size(frq,2)<2 then
72                 error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"),fname,1,1))
73             end
74             if size(frq,2)<>size(varargin(2),2) then
75                 error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),fname,1,2))
76             end
77
78         case 3 then  //frq,db,phi
79             frq=varargin(1);
80             if size(frq,2)<>size(varargin(2),2) then
81                 error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),fname,1,2));
82             end
83             if size(frq,2)<>size(varargin(3),2) then
84                 error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),fname,1,3));
85             end
86             repf=exp(log(10)*varargin(2)/20 + %pi*%i/180*varargin(3));
87
88         else
89             error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,2,4))
90         end
91     else
92         error(msprintf(_("%s: Wrong type for input argument #%d: Linear dynamical system or row vector of floats expected.\n"),fname,1));
93     end;
94     if size(frq,1)==1 then
95         ilf=0;
96     else
97         ilf=1;
98     end
99
100     [mn,n]=size(repf);
101     if and(size(comments,"*")<>[0 mn]) then
102         error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same number of elements expected.\n"),fname,refdim,rhs+1));
103     end
104     //
105
106     repi=imag(repf);
107     repf=real(repf);
108
109     // computing bounds of graphic window
110     mnx=min(-1,min(repf));// to make the critical point visible
111     mxx=max(-1,max(repf));
112
113     if symmetry then
114         mxy=max(0,max(abs(repi)));
115         mny=min(0,-mxy);
116     else
117         mxy=max(0,max(repi));
118         mny=min(0,min(repi));
119     end
120     dx=(mxx-mnx)/30;
121     dy=(mxy-mny)/30;
122     rect=[mnx-dx,mny-dy;mxx+dx,mxy+dy];
123
124     fig=gcf();
125     immediate_drawing=fig.immediate_drawing;
126     fig.immediate_drawing="off";
127
128     ax=gca();
129     if ax.children==[] then
130         ax.data_bounds=rect;
131         ax.axes_visible="on";
132         ax.grid=color("lightgrey")*ones(1,3)
133         ax.title.text=_("Nyquist plot");
134         if sltyp=="c" then
135             ax.x_label.text=_("Re(h(2iπf))");
136             ax.y_label.text=_("Im(h(2iπf))");
137         elseif sltyp=="x" then
138             ax.x_label.text=_("Re");
139             ax.y_label.text=_("Im");
140         else
141             ax.x_label.text=_("Re(h(exp(2iπf*dt)))");
142             ax.y_label.text=_("Im(h(exp(2iπf*dt)))");
143         end
144     else
145         rect= ax.data_bounds
146         mnx=rect(1,1);
147         mxx=rect(2,1)
148         mny=rect(1,2)
149         mxy=rect(2,2)
150     end
151
152     // drawing the curves
153     splitf($+1)=n+1;
154
155     ksplit=1;sel=splitf(ksplit):splitf(ksplit+1)-1;
156     R=[repf(:,sel)];  I=[repi(:,sel)];
157     F=frq(:,sel);
158     for ksplit=2:size(splitf,"*")-1
159         sel=splitf(ksplit):splitf(ksplit+1)-1;
160         R=[R %nan(ones(mn,1)) repf(:,sel)];
161         I=[I %nan(ones(mn,1)) repi(:,sel)];
162         F=[F %nan(ones(size(frq,1),1)) frq(:,sel)];
163     end
164     Curves=[]
165
166     kf=1
167     if symmetry then
168         for k=1:mn
169             xpoly([R(k,:) R(k,$:-1:1)],[I(k,:) -I(k,$:-1:1)]);
170             e=gce();e.foreground=k;
171             e.display_function = "formatNyquistTip";
172             e.display_function_data = [F(kf,:) -1*F(kf,$:-1:1)];
173             Curves=[Curves,e];
174             kf=kf+ilf;
175         end
176     else
177         for k=1:mn
178             xpoly(R(k,:),I(k,:));
179             e=gce();e.foreground=k;
180             e.display_function = "formatNyquistTip";
181             e.display_function_data = F(kf,:);
182             Curves=[Curves,e];
183             kf=kf+ilf;
184         end
185     end
186     clear R I
187
188     kk=1;p0=[repf(:,kk) repi(:,kk)];ks=1;d=0;
189     dx=rect(2,1)-rect(1,1);
190     dy=rect(2,2)-rect(1,2);
191     dx2=dx^2;
192     dy2=dy^2;
193
194     // collect significant frequencies along the curve
195     //-------------------------------------------------------
196     Ic=min(cumsum(sqrt((diff(repf,1,"c").^2)/dx2+ (diff(repi,1,"c").^2)/dy2),2),"r");
197     kk=1;
198     L=0;
199     DIc=0.2;
200     while %t
201         if isempty(Ic) then
202             break
203         end
204         ksup=find(Ic-L>DIc);
205         if ksup==[] then break,end
206         kk1=min(ksup);
207         L=Ic(kk1);
208         Ic(1:kk1)=[];
209         kk=kk+kk1;
210
211         if min(abs(frq(:,ks($))-frq(:,kk))./abs(frq(:,kk)))>0.001 then
212             if min(sqrt(((repf(:,ks)-repf(:,kk)*ones(ks)).^2)/dx2+..
213                 ((repi(:,ks)-repi(:,kk)*ones(ks)).^2)/dy2)) >DIc then
214                 ks=[ks kk];
215                 d=0;
216             end
217         end
218     end
219     if ks($)~=n then
220         if min(((repf(:,ks(1))-repf(:,n)).^2)/dx2+((repi(:,ks(1))-repi(:,n)).^2)/dy2)>0.01 then
221             ks=[ks n];
222         end
223     end
224     // display of parametrization (frequencies along the curve)
225     //-------------------------------------------------------
226     kf=1
227     if ks($)<size(repf,2) then last=$;else last=$-1;end
228     for k=1:mn,
229         L=[];
230         for kks=ks
231             xstring(repf(k,kks),repi(k,kks),msprintf("%-0.3g",frq(kf,kks)),0);
232             e=gce();e.font_foreground=k;
233             L=[e L];
234             if symmetry&(abs(repi(k,kks))>mxy/20) then //not to overlap labels
235                 xstring(repf(k,kks),-repi(k,kks),msprintf("%-0.3g",-frq(kf,kks)),0);
236                 e=gce();e.font_foreground=k;
237                 L=[e L];
238             end
239         end
240         L=glue(L);
241         A=[];
242
243         if size(ks,"*")>1 then
244             dr=repf(k,ks(1:last)+1)-repf(k,ks(1:last));
245             di=repi(k,ks(1:last)+1)-repi(k,ks(1:last));
246             dd=1500*sqrt((dr/dx).^2+(di/dy).^2);
247             dr=dr./dd;
248             di=di./dd;
249             // we should use xarrows or xsegs here.
250             // However their displayed arrow size depends
251             // on the data bounds and we want to avoid this
252             if symmetry then
253                 xx=[repf(k,ks(1:last))         repf(k,ks(last:-1:1))+dr($:-1:1) ;
254                 repf(k,ks(1:last))+dr      repf(k,ks(last:-1:1))]
255                 yy=[repi(k,ks(1:last))        -repi(k,ks(last:-1:1))-di($:-1:1) ;
256                 repi(k,ks(1:last))+di     -repi(k,ks(last:-1:1))]
257             else
258                 xx=[repf(k,ks(1:last))    ;
259                 repf(k,ks(1:last))+dr]
260                 yy=[repi(k,ks(1:last));
261                 repi(k,ks(1:last))+di]
262             end
263             xpolys(xx,yy)
264             //xarrows([repf(k,ks(1:last));repf(k,ks(1:last))+dr],..
265             //    [repi(k,ks(1:last));repi(k,ks(1:last))+di],1.5)
266             A=gce();
267             A.children.arrow_size_factor = 1.5;
268             A.children.polyline_style = 4;
269             A.children.foreground=k;
270         end
271
272         kf=kf+ilf;
273         glue([Curves(k) glue([L A])]);
274
275     end;
276
277     if comments<>[] then
278         legend(Curves, comments);
279     end
280     fig.immediate_drawing=immediate_drawing;
281 endfunction