bug 7566 fix: The cacsd module graphic functions (bode, black, nyquist,...)
[scilab.git] / scilab / modules / cacsd / macros / evans.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA
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
10 function evans(n,d,kmax)
11 // Seuil maxi et mini (relatifs) de discretisation en espace
12 // Copyright INRIA
13   
14   smax=0.002;smin=smax/3;
15   nptmax=2000 //nbre maxi de pt de discretisation en k
16   
17   //analyse de la liste d'appel
18   
19   [lhs,rhs]=argn(0)
20   
21   if rhs <= 0 then   // demo
22     n=real(poly([0.1-%i 0.1+%i,-10],'s'));
23     d=real(poly([-1 -2 -%i %i],'s'));
24     evans(n,d,80);
25     return
26   end
27
28   select typeof(n)
29   case 'constant'  then
30     if rhs==2 then kmax=0,end
31   case 'polynomial'  then
32     if rhs==2 then kmax=0,end
33     //-compat next case retained for list/tlist compatibility
34   case 'rational' then
35     if rhs==2 then kmax=d,else kmax=0,end
36     [n,d]=n(2:3)
37   case 'state-space' then
38     if rhs==2 then kmax=d,else kmax=0,end
39     n=ss2tf(n);
40     [n,d]=n(2:3);n=clean(n);d=clean(d);
41   else 
42      error(msprintf(gettext("%s: Wrong type for input argument #%d: A linear dynamical system or a polynomial expected.\n"),"evans",1));
43   end
44   if prod(size(n))<>1 then
45     error(msprintf(gettext("%s: Wrong value for input argument #%d: Single input, single output system expected.\n"),"evans",1));
46   end
47   if kmax<=0 then
48     nm=min([degree(n),degree(d)])
49     fact=norm(coeff(d),2)/norm(coeff(n),2)
50     kmax=round(500*fact),
51   end
52   //
53   //calcul de la discretisation en k et racines associees
54   nroots=roots(n);racines=roots(d);
55   if nroots==[] then
56     nrm=max([norm(racines,1),norm(roots(d+kmax*n),1)])
57   else
58     nrm=max([norm(racines,1),norm(nroots,1),norm(roots(d+kmax*n),1)])
59   end
60   md=degree(d)
61   //
62   ord=1:md;kk=0;nr=1;k=0;pas=0.99;fin='no';
63
64
65   while fin=='no' then
66     k=k+pas
67     r=roots(d+k*n);r=r(ord)
68     dist=max(abs(racines(:,nr)-r))/nrm
69     //
70     point='nok'
71     if dist <smax then //pas correct
72       point='ok'
73     else //pas trop grand ou ordre incorrect
74          // on cherche l'ordre qui minimise la distance
75          ix=1:md
76          ord1=[]
77          for ky=1:md
78            yy=r(ky)
79            mn=10*dist*nrm
80            for kx=1:md
81              if ix(kx)>0 then
82                if  abs(yy-racines(kx,nr)) < mn then
83                  mn=abs(yy-racines(kx,nr))
84                  kmn=kx
85                end
86              end
87            end
88            ix(kmn)=0
89            ord1=[ord1 kmn]
90          end
91          r(ord1)=r
92          dist=max(abs(racines(:,nr)-r))/nrm
93          if dist <smax then
94            point='ok',
95            ord(ord1)=ord
96          else
97            k=k-pas,pas=pas/2.5
98          end
99     end
100     if dist<smin then
101       //pas trop petit
102       pas=2*pas;
103     end
104     if point=='ok' then
105       racines=[racines,r];kk=[kk,k];nr=nr+1
106       if k>kmax then fin='kmax',end
107       if nr>nptmax then fin='nptmax',end
108     end
109   end
110   //draw the axis
111   x1 =[nroots;matrix(racines,md*nr,1)];
112   xmin=min(real(x1));xmax=max(real(x1))
113   ymin=min(imag(x1));ymax=max(imag(x1))
114   dx=abs(xmax-xmin)*0.05
115   dy=abs(ymax-ymin)*0.05
116   if dx<1d-10, dx=0.01,end
117   if dy<1d-10, dy=0.01,end
118   legs=[],lstyle=[];lhandle=[]
119   rect=[xmin-dx;ymin-dy;xmax+dx;ymax+dy];
120   fig=gcf();
121   immediate_drawing = fig.immediate_drawing;
122   fig.immediate_drawing = 'off';
123   a=gca()
124   a.data_bounds=[rect(1) rect(2);rect(3) rect(4)]
125   if nroots<>[] then 
126     xpoly(real(nroots),imag(nroots))
127     e=gce();e.line_mode='off';e.mark_mode='on';
128     e.mark_size_unit="point";e.mark_size=7;e.mark_style=5;
129     legs=[legs; _("open loop zeroes")]
130     lhandle=[lhandle; e];
131   end
132   if racines<>[] then 
133     xpoly(real(racines(:,1)),imag(racines(:,1)))
134     e=gce();e.line_mode='off';e.mark_mode='on';
135     e.mark_size_unit="point";e.mark_size=7;e.mark_style=2;
136     legs=[legs;_("open loop poles")]
137     lhandle=[lhandle; e];
138   end
139
140   dx=max(abs(xmax-xmin),abs(ymax-ymin));
141   //plot the zeros locations
142
143
144   //computes and draw the asymptotic lines
145   m=degree(n);q=md-m
146   if q>0 then
147     la=0:q-1;
148     so=(sum(racines(:,1))-sum(nroots))/q
149     i1=real(so);i2=imag(so);
150     if prod(size(la))<>1 then
151       ang1=%pi/q*(ones(la)+2*la)
152       x1=dx*cos(ang1),y1=dx*sin(ang1)
153     else
154       x1=0,y1=0,
155     end
156     if md==2,
157       if coeff(d,md)<0 then
158         x1=0*ones(2),y1=0*ones(2)
159       end,
160     end;
161     if max(k)>0 then
162       xpoly(i1,i2);
163       e=gce();
164       legs=[legs;_("asymptotic directions")]
165       lhandle=[lhandle; e];
166       axes = gca();
167       axes.clip_state = "clipgrf";
168       for i=1:q,xsegs([i1,x1(i)+i1],[i2,y1(i)+i2]),end,
169       
170       axes.clip_state = "off";
171     end
172   end;
173
174   //lieu de evans
175   [n1,n2]=size(racines);
176
177   plot2d(real(racines)',imag(racines)',style=2+(1:n2));
178   legend(lhandle,legs);
179   xtitle(_("Evans root locus"),_("Real axis"),_("Imaginary axis"));
180   fig.immediate_drawing = immediate_drawing;
181
182   if fin=='nptmax' then
183     warning(msprintf(gettext("%s: Curve truncated to the first %d discretization points.\n"),"evans",nptmax))
184   end
185 endfunction