Bug 3783 fixed + non reg test
[scilab.git] / scilab / modules / cacsd / macros / trfmod.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA - 
3 // 
4 // This file must be used under the terms of the CeCILL.
5 // This source file is licensed as described in the file COPYING, which
6 // you should have received as part of this distribution.  The terms
7 // are also available at    
8 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
9
10 function h=trfmod(h,job)
11 // hm=trfmod(h [,job])
12 // To visualize the pole-zero structure of a SISO transfer function h 
13 //     job='p' : visualization of polynomials (default)
14 //     job='f' : visualization of natural frequencies and damping
15 //
16 //!
17   select typeof(h)
18   case 'rational' then
19     if size(h('num'))<>[1 1] then
20       error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"trfmod",1))
21     end
22     flag='r'
23   case 'state-space' then
24     if size(h('D'))<>[1 1] then
25       error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"trfmod",1))
26     end
27     flag='lss'
28     den=real(poly(h('A'),'s'))
29     na=degree(den)
30     c=h(4)
31     [m,i]=maxi(abs(c))
32     ci=c(i)
33     t=eye(h(2))*ci;t(i,:)=[-c(1:i-1), 1, -c(i+1:na)]
34     al=h(2)*t;
35     t=eye(h(2))/ci;t(i,:)=[c(1:i-1)/ci, 1, c(i+1:na)/ci]
36     al=t*al;ai=al(:,i),
37     b=t*h(3)
38     al(:,i)=ai+b
39     num=-(real(poly(al,'s'))-den)*ci
40     h=syslin(h(7),num+h(5)*den,den);
41   else
42     error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n"),"trfmod",1))
43   end
44   
45   
46   //
47   format('v',15)
48   [lhs,rhs]=argn(0)
49   if rhs==1 then job='p',end
50   //
51   if type(h('num'))==1 then h('num')=poly(h('num'),varn(h('den')),'c'),end
52   if type(h('den'))==1 then h('den')=poly(h('den'),varn(h('num')),'c'),end
53   
54   var=varn(h('num')),nv=length(var);
55   while part(var,nv)==' ' then nv=nv-1,end;var=part(var,1:nv);
56   
57   fnum=polfact(h('num'))
58   fden=polfact(h('den'))
59   g=coeff(fnum(1))/coeff(fden(1))
60   nn=prod(size(fnum))
61   nd=prod(size(fden))
62   //
63   num=[]
64   for in=2:nn
65     p=fnum(in)
66     if job=='p' then
67       num=[num;pol2str(p)]
68     else
69       if degree(p)==2 then
70         p=coeff(p)
71         omeg=sqrt(p(1))
72         xsi=p(2)/(2*omeg)
73         num=[num;string(omeg)+'    '+string(xsi)]
74       else
75         num=[num;string(-coeff(p,0))]
76       end
77     end
78   end
79   //
80   den=[];
81   for id=2:nd
82     p=fden(id)
83     if job=='p' then
84       den=[den;pol2str(p)]
85     else
86       if degree(p)==2 then
87         p=coeff(p)
88         omeg=sqrt(p(1))
89         xsi=p(2)/(2*omeg)
90         den=[den;string(omeg)+'    '+string(xsi)]
91       else
92         den=[den;string(-coeff(p,0))]
93       end
94     end
95   end
96   
97   txt=[_("Gain :");string(g);_("Numerator :");num;_("Denominator :");den]
98   
99   id=[]
100   if job=='p' then
101     tit=[gettext("Irreducible Factors of transfer function (click below)")]
102   else
103     tit=[gettext("Irreducible Factors of transfer function natural frequency and damping factor (click below)")]
104   end
105   while id==[] then
106     t=x_dialog(tit,txt)
107     id=find(t==_("Denominator :"))
108   end
109   txt=t;
110   
111   tgain=txt(2)
112   tnum=txt(4:id-1)
113   tden=txt(id+1:prod(size(txt)))
114   execstr(var+'=poly(0,'''+var+''')')
115   num=1
116   for in=1:prod(size(tnum))
117     txt=tnum(in)
118     if length(txt)==0 then txt=' ',end
119     if job=='p' then
120       t=' ';
121       for k=1:length(txt),
122         tk=part(txt,k),
123         if tk<>' ' then t=t+tk,end
124       end
125       f=1;if t<>' ' then f=evstr(t),end
126     else
127       if txt==part(' ',1:length(txt)) then
128         f=1
129       else
130         f=evstr(txt)
131         select prod(size(f))
132         case 1 then
133           f=poly(f,var)
134         case 2 then
135           f=poly([f(1)*f(1), 2*f(1)*f(2),1],var,'c')
136         else 
137           error(msprintf(gettext("%s: Incorrect answer.\n"),"trfmod"))
138         end
139       end
140     end
141     num=num*f
142   end
143   //
144   den=1
145   for id=1:prod(size(tden))
146     txt=tden(id);
147     if length(txt)==0 then txt=' ',end
148     if job=='p' then
149       t=' ';
150       for k=1:length(txt),
151         tk=part(txt,k),
152         if tk<>' ' then t=t+tk,end
153       end
154       f=1;if t<>' ' then f=evstr(t),end
155     else
156       if txt==part(' ',1:length(txt)) then
157         f=1
158       else
159         f=evstr(txt)
160         select prod(size(f))
161         case 1 then
162           f=poly(f,var)
163         case 2 then
164           f=poly([f(1)*f(1), 2*f(1)*f(2),1],var,'c')
165         else 
166           error(msprintf(gettext("%s: Incorrect answer.\n"),"trfmod"))
167         end
168       end
169     end
170     den=den*f
171   end
172   x=evstr(tgain)/coeff(den,degree(den))
173   h('num')=num*x
174   h('den')=den/coeff(den,degree(den))
175   format(10)
176   if flag=='lss' then h=tf2ss(h),end
177 endfunction