diff options
Diffstat (limited to 'modules/cacsd/macros/trfmod.sci')
| -rwxr-xr-x | modules/cacsd/macros/trfmod.sci | 179 |
1 files changed, 179 insertions, 0 deletions
diff --git a/modules/cacsd/macros/trfmod.sci b/modules/cacsd/macros/trfmod.sci new file mode 100755 index 000000000..9c8da5dea --- /dev/null +++ b/modules/cacsd/macros/trfmod.sci @@ -0,0 +1,179 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) INRIA - +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + +function h=trfmod(h,job) + // hm=trfmod(h [,job]) + // To visualize the pole-zero structure of a SISO transfer function h + // job='p' : visualization of polynomials (default) + // job='f' : visualization of natural frequencies and damping + // + //! + select typeof(h) + case "rational" then + if size(h("num"))<>[1 1] then + error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"trfmod",1)) + end + flag="r" + case "state-space" then + if size(h("D"))<>[1 1] then + error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"trfmod",1)) + end + flag="lss" + den=real(poly(h("A"),"s")) + na=degree(den) + c=h(4) + [m,i]=max(abs(c)) + ci=c(i) + t=eye(h(2))*ci;t(i,:)=[-c(1:i-1), 1, -c(i+1:na)] + al=h(2)*t; + t=eye(h(2))/ci;t(i,:)=[c(1:i-1)/ci, 1, c(i+1:na)/ci] + al=t*al;ai=al(:,i), + b=t*h(3) + al(:,i)=ai+b + num=-(real(poly(al,"s"))-den)*ci + h=syslin(h(7),num+h(5)*den,den); + else + error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n"),"trfmod",1)) + end + + + // + ft = format(); + format("v", 15); + + [lhs,rhs]=argn(0) + if rhs==1 then job="p",end + // + if type(h("num"))==1 then h("num")=poly(h("num"),varn(h("den")),"c"),end + if type(h("den"))==1 then h("den")=poly(h("den"),varn(h("num")),"c"),end + + var=varn(h("num")),nv=length(var); + while part(var,nv)==" " then nv=nv-1,end;var=part(var,1:nv); + + fnum=polfact(h("num")) + fden=polfact(h("den")) + g=coeff(fnum(1))/coeff(fden(1)) + nn=prod(size(fnum)) + nd=prod(size(fden)) + // + num=[] + for in=2:nn + p=fnum(in) + if job=="p" then + num=[num;pol2str(p)] + else + if degree(p)==2 then + p=coeff(p) + omeg=sqrt(p(1)) + xsi=p(2)/(2*omeg) + num=[num;string(omeg)+" "+string(xsi)] + else + num=[num;string(-coeff(p,0))] + end + end + end + // + den=[]; + for id=2:nd + p=fden(id) + if job=="p" then + den=[den;pol2str(p)] + else + if degree(p)==2 then + p=coeff(p) + omeg=sqrt(p(1)) + xsi=p(2)/(2*omeg) + den=[den;string(omeg)+" "+string(xsi)] + else + den=[den;string(-coeff(p,0))] + end + end + end + + txt=[_("Gain :");string(g);_("Numerator :");num;_("Denominator :");den] + + id=[] + if job=="p" then + tit=[gettext("Irreducible Factors of transfer function (click below)")] + else + tit=[gettext("Irreducible Factors of transfer function natural frequency and damping factor (click below)")] + end + while id==[] then + t=x_dialog(tit,txt) + id=find(t==_("Denominator :")) + end + txt=t; + + tgain=txt(2) + tnum=txt(4:id-1) + tden=txt(id+1:prod(size(txt))) + execstr(var+"=poly(0,''"+var+"'')") + num=1 + for in=1:prod(size(tnum)) + txt=tnum(in) + if length(txt)==0 then txt=" ",end + if job=="p" then + t=" "; + for k=1:length(txt), + tk=part(txt,k), + if tk<>" " then t=t+tk,end + end + f=1;if t<>" " then f=evstr(t),end + else + if txt==part(" ",1:length(txt)) then + f=1 + else + f=evstr(txt) + select prod(size(f)) + case 1 then + f=poly(f,var) + case 2 then + f=poly([f(1)*f(1), 2*f(1)*f(2),1],var,"c") + else + error(msprintf(gettext("%s: Incorrect answer.\n"),"trfmod")) + end + end + end + num=num*f + end + // + den=1 + for id=1:prod(size(tden)) + txt=tden(id); + if length(txt)==0 then txt=" ",end + if job=="p" then + t=" "; + for k=1:length(txt), + tk=part(txt,k), + if tk<>" " then t=t+tk,end + end + f=1;if t<>" " then f=evstr(t),end + else + if txt==part(" ",1:length(txt)) then + f=1 + else + f=evstr(txt) + select prod(size(f)) + case 1 then + f=poly(f,var) + case 2 then + f=poly([f(1)*f(1), 2*f(1)*f(2),1],var,"c") + else + error(msprintf(gettext("%s: Incorrect answer.\n"),"trfmod")) + end + end + end + den=den*f + end + x=evstr(tgain)/coeff(den,degree(den)) + h("num")=num*x + h("den")=den/coeff(den,degree(den)) + format(ft(2),ft(1)); + if flag=="lss" then h=tf2ss(h),end +endfunction |