bug 7566 fix: The cacsd module graphic functions (bode, black, nyquist,...)
[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   fig=gcf();
98   immediate_drawing=fig.immediate_drawing;
99   fig.immediate_drawing="off";
100
101   ax=gca();
102   if ax.children==[] then
103     ax.data_bounds=rect;
104     ax.axes_visible="on";
105     ax.grid=color("lightgrey")*ones(1,3)
106     ax.title.text=_("Nyquist plot");
107     if sltyp=="c" then
108       ax.x_label.text=_("Re(h(2iπf))");
109       ax.y_label.text=_("Im(h(2iπf))");
110     elseif sltyp=="x" then
111       ax.x_label.text=_("Re");
112       ax.y_label.text=_("Im");
113     else
114       ax.x_label.text=_("Re(h(exp(2iπf*dt)))");
115       ax.y_label.text=_("Im(h(exp(2iπf*dt)))");
116     end
117   else
118     ax.data_bounds=[min(ax.data_bounds(1,:),rect(1,:));max(ax.data_bounds(2,:),rect(2,:))];
119   end
120   // drawing the curves
121   splitf($+1)=n+1;
122
123   ksplit=1;sel=splitf(ksplit):splitf(ksplit+1)-1;
124   R=[repf(:,sel)];  I=[repi(:,sel)];
125   F=frq(:,sel);
126   for ksplit=2:size(splitf,"*")-1
127     sel=splitf(ksplit):splitf(ksplit+1)-1;
128     R=[R %nan(ones(mn,1)) repf(:,sel)];
129     I=[I %nan(ones(mn,1)) repi(:,sel)];
130     F=[F %nan(ones(size(frq,1),1)) frq(:,sel)];
131   end
132   Curves=[]
133
134   kf=1
135   for k=1:mn
136     xpoly([R(k,:) R(k,$:-1:1)],[I(k,:) -I(k,$:-1:1)]);
137     e=gce();e.foreground=k;
138     datatipInitStruct(e,"formatfunction","formatNyquistTip","freq",[F(kf,:) F(kf,$:-1:1)])
139     Curves=[Curves,e];
140     kf=kf+ilf;
141   end
142   clear R I
143
144   kk=1;p0=[repf(:,kk) repi(:,kk)];ks=1;d=0;
145   dx=rect(2,1)-rect(1,1);
146   dy=rect(2,2)-rect(1,2);
147   dx2=dx^2;
148   dy2=dy^2;
149
150   // collect significant frequencies along the curve
151   //-------------------------------------------------------
152   Ic=min(cumsum(sqrt((diff(repf,1,"c").^2)/dx2+ (diff(repi,1,"c").^2)/dy2),2),"r");
153   kk=1;
154   L=0;
155   DIc=0.2;
156   while %t
157     ksup=find(Ic-L>DIc);
158     if ksup==[] then break,end
159     kk1=min(ksup);
160     L=Ic(kk1);
161     Ic(1:kk1)=[];
162     kk=kk+kk1;
163
164     if min(abs(frq(:,ks($))-frq(:,kk))./abs(frq(:,kk)))>0.001 then
165       if min(sqrt(((repf(:,ks)-repf(:,kk)*ones(ks)).^2)/dx2+..
166                    ((repi(:,ks)-repi(:,kk)*ones(ks)).^2)/dy2)) >DIc then
167         ks=[ks kk];
168         d=0;
169       end
170     end
171   end
172   if ks($)~=n then
173     if min(((repf(:,ks(1))-repf(:,n))^2)/dx2+((repi(:,ks(1))-repi(:,n))^2)/dy2)>0.01 then
174       ks=[ks n];
175     end
176   end
177   // display of parametrization (frequencies along the curve)
178   //-------------------------------------------------------
179   kf=1
180   if ks($)<size(repf,2) then last=$;else last=$-1;end
181   for k=1:mn,
182     L=[];
183     for kks=ks
184       xstring(repf(k,kks),repi(k,kks),msprintf("%-0.3g",frq(kf,kks)),0);
185       e=gce();e.font_foreground=k;
186       L=[e L];
187       if abs(repi(k,kks))>mxy/20 then //not to overlap labels
188         xstring(repf(k,kks),-repi(k,kks),msprintf("%-0.3g",-frq(kf,kks)),0);
189         e=gce();e.font_foreground=k;
190         L=[e L];
191       end
192     end
193     L=glue(L);
194     A=[];
195
196     if size(ks,"*")>1 then
197       dr=repf(k,ks(1:last)+1)-repf(k,ks(1:last));
198       di=repi(k,ks(1:last)+1)-repi(k,ks(1:last));
199       dd=1500*sqrt((dr/dx).^2+(di/dy).^2);
200       dr=dr./dd;
201       di=di./dd;
202       // we should use xarrows or xsegs here.
203       // However their displayed arrow size depends
204       // on the data bounds and we want to avoid this
205       xx=[repf(k,ks(1:last))         repf(k,ks(last:-1:1))+dr($:-1:1) ;
206           repf(k,ks(1:last))+dr      repf(k,ks(last:-1:1))]
207       yy=[repi(k,ks(1:last))        -repi(k,ks(last:-1:1))-di($:-1:1) ;
208           repi(k,ks(1:last))+di     -repi(k,ks(last:-1:1))]
209       xpolys(xx,yy)
210       //xarrows([repf(k,ks(1:last));repf(k,ks(1:last))+dr],..
211       //    [repi(k,ks(1:last));repi(k,ks(1:last))+di],1.5)
212       A=gce();
213       A.children.arrow_size_factor = 1.5;
214       A.children.polyline_style = 4;
215       A.children.foreground=k;
216     end
217
218     kf=kf+ilf;
219     glue([Curves(k) glue([L A])]);
220
221   end;
222
223   if comments<>[] then
224     legend(Curves($:-1:1),comments);
225   end
226   fig.immediate_drawing=immediate_drawing;
227 endfunction
228
229 function str=formatNyquistTip(curve,pt,index)
230 //This function is called by the datatip mechanism to format the tip
231 //string for the nyquist curves.
232   ud=datatipGetStruct(curve);
233   if index<>[] then
234     f=ud.freq(index);
235   else //interpolated
236     [d,ptp,i,c]=orthProj(curve.data,pt);
237     f=ud.freq(i)+(ud.freq(i+1)-ud.freq(i))*c;
238   end
239   str=msprintf("%.4g%+.4gi\n%.4g"+_("Hz"), pt,f);
240 endfunction