suggestion about evans
Iai Masafumi ax
iai at axelspace.com
Sun Oct 23 02:00:26 CEST 2011
I suggest to review the code of "evans" function to see if setting
"clip_state" off is useful to users.
What I wanted is to generate a figure like fig_clip_state_is_clipgrf.png
as attached. However with the current "evans" function, I only got a
figure like fig_clip_state_is_off.png, in which some lines are drawn
outside the plot area. This is because "clip_state" is set off. Changing
the value of "clip_state" of the current axes before or after calling
"evans" function did not affect the result. I was successful when I
commented out a line:
axes.clip_state = "off";
in evans.sci. This should be good unless other side effects occur.
So I think it needs re-consideration which value "clip_state" should be
set to.
The code below should generate the same figures as attachments. My
"evans_clip" function is identical to Scilab's evans function except
that a line [axes.clip_state = "off"] is commented out.
Iai
////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
function evans_clip(n,d,kmax)
// Seuil maxi et mini (relatifs) de discretisation en espace
// Copyright INRIA
smax=0.002;smin=smax/3;
nptmax=2000 //nbre maxi de pt de discretisation en k
//Check calling sequence
[lhs,rhs]=argn(0)
if rhs <= 0 then // demonstration
n=real(poly([0.1-%i 0.1+%i,-10],"s"));
d=real(poly([-1 -2 -%i %i],"s"));
evans(n,d,80);
return
end
select typeof(n)
case "polynomial" then
if rhs==2 then kmax=0,end
case "rational" then
if rhs==2 then kmax=d,else kmax=0,end
[n,d]=n(2:3)
case "state-space" then
if rhs==2 then kmax=d,else kmax=0,end
n=ss2tf(n);
[n,d]=n(2:3);n=clean(n);d=clean(d);
else
error(msprintf(_("%s: Wrong type for input argument #%d: A linear
dynamical system or a polynomial expected.\n"),"evans",1));
end
if prod(size(n))<>1 then
error(msprintf(_("%s: Wrong value for input argument #%d: Single
input, single output system expected.\n"),"evans",1));
end
if degree(n)==0°ree(d)==0 then
error(msprintf(_("%s: The given system has no poles and no
zeros.\n"),"evans"));
end
if kmax<=0 then
nm=min([degree(n),degree(d)])
fact=norm(coeff(d),2)/norm(coeff(n),2)
kmax=round(500*fact),
end
//
//Compute the discretization for "k" and the associated roots
nroots=roots(n);racines=roots(d);
if nroots==[] then
nrm=max([norm(racines,1),norm(roots(d+kmax*n),1)])
else
nrm=max([norm(racines,1),norm(nroots,1),norm(roots(d+kmax*n),1)])
end
md=degree(d)
//
ord=1:md;kk=0;nr=1;k=0;pas=0.99;fin="no";
klim=gsort(krac2(rlist(n,d,"c")),"g","i")
ilim=1
while fin=="no" then
k=k+pas
r=roots(d+k*n);r=r(ord)
dist=max(abs(racines(:,nr)-r))/nrm
//
point=%f
if dist <smax then //pas correct
if k-pas<klim(ilim)& k>klim(ilim) then,
k=klim(ilim);
r=roots(d+k*n);r=r(ord)
end
if k>klim(ilim) then ilim=min(ilim+1,size(klim,"*"));end
point=%t
else //Too big step or incorrect root order
// look for a root order that minimize the distance
ix=1:md
ord1=[]
for ky=1:md
yy=r(ky)
mn=10*dist*nrm
for kx=1:md
if ix(kx)>0 then
if abs(yy-racines(kx,nr)) < mn then
mn=abs(yy-racines(kx,nr))
kmn=kx
end
end
end
ix(kmn)=0
ord1=[ord1 kmn]
end
r(ord1)=r
dist=max(abs(racines(:,nr)-r))/nrm
if dist <smax then
point=%t,
ord(ord1)=ord
else
k=k-pas,pas=pas/2.5
end
end
if dist<smin then
//KToo small step
pas=2*pas;
end
if point then
racines=[racines,r];kk=[kk,k];nr=nr+1
if k>kmax then fin="kmax",end
if nr>nptmax then fin="nptmax",end
end
end
//draw the axis
x1 =[nroots;matrix(racines,md*nr,1)];
xmin=min(real(x1));xmax=max(real(x1))
ymin=min(imag(x1));ymax=max(imag(x1))
dx=abs(xmax-xmin)*0.05
dy=abs(ymax-ymin)*0.05
if dx<1d-10, dx=0.01,end
if dy<1d-10, dy=0.01,end
legs=[],lstyle=[];lhandle=[]
rect=[xmin-dx;ymin-dy;xmax+dx;ymax+dy];
f=gcf();
immediate_drawing= f.immediate_drawing;
f.immediate_drawing = "off";
a=gca();
if a.children==[]
a.data_bounds=[rect(1) rect(2);rect(3) rect(4)];
a.axes_visible="on";
a.title.text=_("Evans root locus");
a.x_label.text=_("Real axis");
a.y_label.text=_("Imaginary axis");
else //enlarge the boundaries
a.data_bounds=[min(a.data_bounds(1,:),[rect(1) rect(2)]);
max(a.data_bounds(2,:),[rect(3) rect(4)])];
end
if nroots<>[] then
xpoly(real(nroots),imag(nroots))
e=gce();e.line_mode="off";e.mark_mode="on";
e.mark_size_unit="point";e.mark_size=7;e.mark_style=5;
legs=[legs; _("open loop zeroes")]
lhandle=[lhandle; e];
end
if racines<>[] then
xpoly(real(racines(:,1)),imag(racines(:,1)))
e=gce();e.line_mode="off";e.mark_mode="on";
e.mark_size_unit="point";e.mark_size=7;e.mark_style=2;
legs=[legs;_("open loop poles")]
lhandle=[lhandle; e];
end
dx=max(abs(xmax-xmin),abs(ymax-ymin));
//plot the zeros locations
//computes and draw the asymptotic lines
m=degree(n);q=md-m
if q>0 then
la=0:q-1;
so=(sum(racines(:,1))-sum(nroots))/q
i1=real(so);i2=imag(so);
if prod(size(la))<>1 then
ang1=%pi/q*(ones(la)+2*la)
x1=dx*cos(ang1),y1=dx*sin(ang1)
else
x1=0,y1=0,
end
if md==2,
if coeff(d,md)<0 then
x1=0*ones(2),y1=0*ones(2)
end,
end;
if max(k)>0 then
xpoly(i1,i2);
e=gce();
legs=[legs;_("asymptotic directions")]
lhandle=[lhandle; e];
axes = gca();
axes.clip_state = "clipgrf";
for i=1:q,xsegs([i1,x1(i)+i1],[i2,y1(i)+i2]),end,
//axes.clip_state = "off"; // !! OMITTED !! //
end
end;
[n1,n2]=size(racines);
// assign the colors for each root locus
cmap=f.color_map;cols=1:size(cmap,1);
if a.background==-2 then
cols(and(cmap==1,2))=[]; //remove white
elseif a.background==-1 then
cols(and(cmap==0,2))=[]; //remove black
else
cols(a.background)=[];
end
cols=cols(modulo(0:n1-1,size(cols,"*"))+1);
//draw the root locus
xpolys(real(racines)',imag(racines)',cols)
//set info for datatips
E=gce();
for k=1:size(E.children,"*")
datatipInitStruct(E.children(k),"formatfunction","formatEvansTip","K",kk)
end
c=captions(lhandle,legs($:-1:1),"in_upper_right")
c.background=a.background;
f.immediate_drawing = immediate_drawing;
if fin=="nptmax" then
warning(msprintf(gettext("%s: Curve truncated to the first %d
discretization points.\n"),"evans",nptmax))
end
endfunction
function str=formatEvansTip(curve,pt,index)
//this function is called by the datatip mechanism to format the tip
//string for the evans root loci curves
ud=datatipGetStruct(curve);
if index<>[] then
K=ud.K(index)
else //interpolated
[d,ptp,i,c]=orth_proj(curve.data,pt)
K=ud.K(i)+(ud.K(i+1)-ud.K(i))*c
end
str=msprintf("r: %.4g %+.4g i\nK: %.4g", pt,K);
endfunction
figure(1);
clf;
subplot(2,2,1);
Gc=poly([-6+2*%i -6-2*%i],"s","roots")/poly(0,"s","roots")
G=1/poly([-2 -5],"s","roots")
H=G*Gc
evans_clip(H,100)
a=gca();
a.data_bounds = [[-15; 10] [-8; 8]];
subplot(2,2,2);
Gc=poly([-3.5+2*%i -3.5-2*%i],"s","roots")/poly(0,"s","roots")
G=1/poly([-2 -5],"s","roots")
H=G*Gc
evans_clip(H,100)
a=gca();
a.data_bounds = [[-15; 10] [-8; 8]];
subplot(2,2,3);
Gc=poly([-6+2*%i -6-2*%i],"s","roots")/poly(0,"s","roots")
G=1/poly([-2 -5],"s","roots")
H=G*Gc
evans_clip(H,100)
a=gca();
a.data_bounds = [[-15; 10] [-8; 8]];
subplot(2,2,4);
xs2png(gcf(), "fig1");
figure(2);
clf;
subplot(2,2,1);
Gc=poly([-6+2*%i -6-2*%i],"s","roots")/poly(0,"s","roots")
G=1/poly([-2 -5],"s","roots")
H=G*Gc
evans(H,100)
a=gca();
a.data_bounds = [[-15; 10] [-8; 8]];
subplot(2,2,2);
Gc=poly([-3.5+2*%i -3.5-2*%i],"s","roots")/poly(0,"s","roots")
G=1/poly([-2 -5],"s","roots")
H=G*Gc
evans(H,100)
a=gca();
a.data_bounds = [[-15; 10] [-8; 8]];
subplot(2,2,3);
Gc=poly([-6+2*%i -6-2*%i],"s","roots")/poly(0,"s","roots")
G=1/poly([-2 -5],"s","roots")
H=G*Gc
evans(H,100)
a=gca();
a.data_bounds = [[-15; 10] [-8; 8]];
subplot(2,2,4);
xs2png(gcf(), "fig2");
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fig_clip_state_is_clipgrf.png
Type: image/png
Size: 13612 bytes
Desc: not available
URL: <https://lists.scilab.org/pipermail/users/attachments/20111023/974353f6/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fig_clip_state_is_off.png
Type: image/png
Size: 13558 bytes
Desc: not available
URL: <https://lists.scilab.org/pipermail/users/attachments/20111023/974353f6/attachment-0001.png>
More information about the users
mailing list