49487e21c9ad4f77a772ba865bb73f7d1dedd489
[scilab.git] / scilab / modules / cacsd / macros / nyquist.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 nyquist(varargin)
10 // Nyquist plot
11 //!
12   rhs=size(varargin);
13
14   if type(varargin(rhs))==10 then
15     comments=varargin(rhs);
16     rhs=rhs-1;
17   else
18     comments=[];
19   end
20   fname="nyquist";//for error messages
21   fmax=[];
22   if or(typeof(varargin(1))==["state-space" "rational"]) then
23     //sys,fmin,fmax [,pas] or sys,frq
24     refdim=1; //for error message
25     sltyp=varargin(1).dt;
26     if rhs==1 then
27       [frq,repf,splitf]=repfreq(varargin(1),1d-3,1d3);
28     elseif rhs==2 then //sys,frq
29       if size(varargin(2),2)<2 then
30         error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"),fname,2,1))
31       end
32       [frq,repf]=repfreq(varargin(1:rhs));
33     elseif or(rhs==(3:4)) then //sys,fmin,fmax [,pas]
34       [frq,repf,splitf]=repfreq(varargin(1:rhs));
35     else
36       error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,1,5))
37     end
38   elseif  type(varargin(1))==1 then
39     //frq,db,phi [,comments] or frq, repf [,comments]
40     refdim=2;
41     sltyp="x";
42     splitf=[];
43     splitf=1;
44     select rhs
45     case 2 then //frq,repf
46       frq=varargin(1);
47       repf=varargin(2);
48       if size(frq,2)<2 then
49         error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"),fname,1,1))
50       end
51       if size(frq,2)<>size(varargin(2),2) then
52         error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),fname,1,2))
53       end
54
55     case 3 then  //frq,db,phi
56       frq=varargin(1);
57       if size(frq,2)<>size(varargin(2),2) then
58         error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),fname,1,2));
59       end
60       if size(frq,2)<>size(varargin(3),2) then
61         error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),fname,1,3));
62       end
63       repf=exp(log(10)*varargin(2)/20 + %pi*%i/180*varargin(3));
64
65     else
66       error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,2,4))
67     end
68   else
69     error(msprintf(_("%s: Wrong type for input argument #%d: Linear dynamical system or row vector of floats expected.\n"),fname,1));
70   end;
71   if size(frq,1)==1 then
72     ilf=0;
73   else
74     ilf=1;
75   end
76
77   [mn,n]=size(repf);
78   if and(size(comments,"*")<>[0 mn]) then
79     error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same number of elements expected.\n"),fname,refdim,rhs+1));
80   end
81   //
82
83   repi=imag(repf);
84   repf=real(repf);
85
86   // computing bounds of graphic window
87   mnx=min(-1,min(repf));// to make the critical point visible
88   mxx=max(-1,max(repf));
89
90   mxy=max(0,max(abs(repi)));
91   mny=min(0,-mxy);
92
93   dx=(mxx-mnx)/30;
94   dy=(mxy-mny)/30;
95   rect=[mnx-dx,mny-dy;mxx+dx,mxy+dy];
96
97   drawlater()
98   ax=gca();
99   if ax.children==[] then
100     ax.data_bounds=rect;
101     ax.axes_visible="on";
102     ax.grid=color("lightgrey")*ones(1,3)
103     ax.title.text=_("Nyquist plot");
104     if sltyp=="c" then
105       ax.x_label.text=_("Re(h(2iπf))");
106       ax.y_label.text=_("Im(h(2iπf))");
107     elseif sltyp=="x" then
108       ax.x_label.text=_("Re");
109       ax.y_label.text=_("Im");
110     else
111       ax.x_label.text=_("Re(h(exp(2iπf*dt)))");
112       ax.y_label.text=_("Im(h(exp(2iπf*dt)))");
113     end
114   else
115     ax.data_bounds=[min(ax.data_bounds(1,:),rect(1,:));max(ax.data_bounds(2,:),rect(2,:))];
116   end
117   // drawing the curves
118   splitf($+1)=n+1;
119
120   ksplit=1;sel=splitf(ksplit):splitf(ksplit+1)-1;
121   R=[repf(:,sel)];  I=[repi(:,sel)];
122   F=frq(:,sel);
123   for ksplit=2:size(splitf,"*")-1
124     sel=splitf(ksplit):splitf(ksplit+1)-1;
125     R=[R %nan(ones(mn,1)) repf(:,sel)];
126     I=[I %nan(ones(mn,1)) repi(:,sel)];
127     F=[F %nan(ones(size(frq,1),1)) frq(:,sel)];
128   end
129   Curves=[]
130
131   kf=1
132   for k=1:mn
133     xpoly([R(k,:) R(k,$:-1:1)],[I(k,:) -I(k,$:-1:1)]);
134     e=gce();e.foreground=k;
135     datatipInitStruct(e,"formatfunction","formatNyquistTip","freq",[F(kf,:) F(kf,$:-1:1)])
136     Curves=[Curves,e];
137     kf=kf+ilf;
138   end
139   clear R I
140
141   kk=1;p0=[repf(:,kk) repi(:,kk)];ks=1;d=0;
142   dx=rect(2,1)-rect(1,1);
143   dy=rect(2,2)-rect(1,2);
144   dx2=dx^2;
145   dy2=dy^2;
146
147   // collect significant frequencies along the curve
148   //-------------------------------------------------------
149   Ic=min(cumsum(sqrt((diff(repf,1,"c").^2)/dx2+ (diff(repi,1,"c").^2)/dy2),2),"r");
150   kk=1;
151   L=0;
152   DIc=0.2;
153   while %t
154     ksup=find(Ic-L>DIc);
155     if ksup==[] then break,end
156     kk1=min(ksup);
157     L=Ic(kk1);
158     Ic(1:kk1)=[];
159     kk=kk+kk1;
160
161     if min(abs(frq(:,ks($))-frq(:,kk))./abs(frq(:,kk)))>0.001 then
162       if min(sqrt(((repf(:,ks)-repf(:,kk)*ones(ks)).^2)/dx2+..
163                    ((repi(:,ks)-repi(:,kk)*ones(ks)).^2)/dy2)) >DIc then
164         ks=[ks kk];
165         d=0;
166       end
167     end
168   end
169   if ks($)~=n then
170     if min(((repf(:,ks(1))-repf(:,n))^2)/dx2+((repi(:,ks(1))-repi(:,n))^2)/dy2)>0.01 then
171       ks=[ks n];
172     end
173   end
174   // display of parametrization (frequencies along the curve)
175   //-------------------------------------------------------
176   kf=1
177   if ks($)<size(repf,2) then last=$;else last=$-1;end
178   for k=1:mn,
179     L=[];
180     for kks=ks
181       xstring(repf(k,kks),repi(k,kks),msprintf("%-0.3g",frq(kf,kks)),0);
182       e=gce();e.font_foreground=k;
183       L=[e L];
184       if abs(repi(k,kks))>mxy/20 then //not to overlap labels
185         xstring(repf(k,kks),-repi(k,kks),msprintf("%-0.3g",-frq(kf,kks)),0);
186         e=gce();e.font_foreground=k;
187         L=[e L];
188       end
189     end
190     L=glue(L);
191     A=[];
192
193     if size(ks,"*")>1 then
194       dr=repf(k,ks(1:last)+1)-repf(k,ks(1:last));
195       di=repi(k,ks(1:last)+1)-repi(k,ks(1:last));
196       dd=1500*sqrt((dr/dx).^2+(di/dy).^2);
197       dr=dr./dd;
198       di=di./dd;
199       // we should use xarrows or xsegs here.
200       // However their displayed arrow size depends
201       // on the data bounds and we want to avoid this
202       xx=[repf(k,ks(1:last))         repf(k,ks(last:-1:1))+dr($:-1:1) ;
203           repf(k,ks(1:last))+dr      repf(k,ks(last:-1:1))]
204       yy=[repi(k,ks(1:last))        -repi(k,ks(last:-1:1))-di($:-1:1) ;
205           repi(k,ks(1:last))+di     -repi(k,ks(last:-1:1))]
206       xpolys(xx,yy)
207       //xarrows([repf(k,ks(1:last));repf(k,ks(1:last))+dr],..
208       //    [repi(k,ks(1:last));repi(k,ks(1:last))+di],1.5)
209       A=gce();
210       A.children.arrow_size_factor = 1.5;
211       A.children.polyline_style = 4;
212       A.children.foreground=k;
213     end
214
215     kf=kf+ilf;
216     glue([Curves(k) glue([L A])]);
217
218   end;
219
220   if comments<>[] then
221     legend(Curves($:-1:1),comments);
222   end
223   drawnow()
224 endfunction
225 function str=formatNyquistTip(curve,pt,index)
226 //This function is called by the datatip mechanism to format the tip
227 //string for the nyquist curves.
228   ud=datatipGetStruct(curve);
229   if index<>[] then
230     f=ud.freq(index);
231   else //interpolated
232     [d,ptp,i,c]=orthProj(curve.data,pt);
233     f=ud.freq(i)+(ud.freq(i+1)-ud.freq(i))*c;
234   end
235   str=msprintf("%.4g%+.4gi\n%.4g"+_("Hz"), pt,f);
236 endfunction