← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to macros/auto/ss2tf.sci

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
macros/auto/ss2tf.sci
 
1
function [h,num,den]=ss2tf(sl,rmax)
 
2
// State-space to transfer function.
 
3
// Syntax:
 
4
//   h=ss2tf(sl) 
 
5
//   h=ss2tf(sl,'b')
 
6
//   h=ss2tf(sl,rmax)
 
7
//
 
8
//   sl   : linear system (syslin list)
 
9
//   h    : transfer matrix 
 
10
//   rmax : optional parameter controlling the conditioning in
 
11
//          block diagonalization method is used (If 'b' is entered
 
12
//          a default value is used)
 
13
//   Method: By default, one uses characteristic polynomial and 
 
14
//   det(A+Eij)=det(A)+C(i,j) where C is the adjugate matrix of A
 
15
//   Other method used : block-diagonalization (generally 
 
16
//   this gives a less accurate result).
 
17
//
 
18
//!
 
19
// Copyright INRIA
 
20
if type(sl)==1|type(sl)==2 then
 
21
  h=sl
 
22
  return
 
23
end
 
24
flag=sl(1);
 
25
if (type(sl)<>16)|flag(1)<>'lss' then
 
26
  error('First argument must be in state-space form')
 
27
end
 
28
if sl(3)==[] then h=sl(5);num=sl(5);den=eye(sl(5));return;end
 
29
if sl(4)==[] then h=sl(5);num=sl(5);den=eye(sl(5));return;end
 
30
if size(sl(2),'*')==0 then
 
31
  h=sl(5)
 
32
  return
 
33
end 
 
34
//1
 
35
domaine=sl(7)
 
36
if type(domaine)==1 then var='z';end
 
37
if domaine=='c' then var='s';end;
 
38
if domaine=='d' then var='z';end;
 
39
if domaine==[] then
 
40
   var='s'; 
 
41
   if type(sl(5))==2 then var=varn(sl(5));end
 
42
end
 
43
//
 
44
[lhs,rhs]=argn(0);
 
45
meth='p'
 
46
if rhs==2 then
 
47
  if type(rmax)==10 then  
 
48
    meth=part(rmax,1),
 
49
  else 
 
50
    meth='b'
 
51
  end
 
52
end
 
53
//
 
54
select meth
 
55
case 'b' // 
 
56
a=sl(2)
 
57
[n1,n1]=size(a)
 
58
z=poly(0,var);
 
59
if rhs==1 then
 
60
  [a,x,bs]=bdiag(a)
 
61
else
 
62
  [a,x,bs]=bdiag(a,rmax)
 
63
end
 
64
k=1;m=[];v=ones(1,n1);den=1;
 
65
for n=bs';k1=k:k-1+n;
 
66
// leverrier
 
67
             h=z*eye(n,n)-a(k1,k1)
 
68
             f=eye(n,n)
 
69
             for kl=1:n-1,
 
70
                b=h*f,
 
71
                d=-sum(diag(b))/kl,
 
72
                f=b+eye()*d,
 
73
             end;
 
74
             d=sum(diag(h*f))/n
 
75
//
 
76
          den=den*d;
 
77
          l=[1:k-1,k+n:n1] ,
 
78
          if l<>[] then v(l)=v(l)*d,end
 
79
          m=[m,x(:,k1)*f];
 
80
          k=k+n;
 
81
end;
 
82
if lhs==3 then h=sl(5),num=sl(4)*m*diag(v)*(x\sl(3));return;end
 
83
m=sl(4)*m*diag(v)*(x \ sl(3))+sl(5)*den;
 
84
[m1,n1]=size(m);[m,den]=simp(m,den*ones(m1,n1))
 
85
h=syslin(domaine,m,den)
 
86
case 'p' then
 
87
  Den=real(poly(sl(2),var))
 
88
  na=degree(Den);den=[];
 
89
  [m,n]=size(sl(5))
 
90
  c=sl(4)
 
91
  for l=1:m
 
92
    [m,i]=maxi(abs(c(l,:)));
 
93
    if m<>0 then
 
94
      ci=c(l,i)
 
95
      t=eye(na,na)*ci;t(i,:)=[-c(l,1:i-1), 1, -c(l,i+1:na)]
 
96
      al=sl(2)*t;
 
97
      t=eye(na,na)/ci;t(i,:)=[c(l,1:i-1)/ci, 1, c(l,i+1:na)/ci]
 
98
      al=t*al;ai=al(:,i),
 
99
      b=t*sl(3)
 
100
      for k=1:n
 
101
        al(:,i)=ai+b(:,k);
 
102
        [nlk,dlk]=simp(real(poly(al,var)),Den)
 
103
        den(l,k)=dlk;
 
104
        num(l,k)=-(nlk-dlk)*ci
 
105
      end
 
106
    else
 
107
      num(l,1:n)=0*ones(1,n);
 
108
      den(l,1:n)=ones(1,n);
 
109
    end
 
110
  end
 
111
  if lhs==3 then h=sl(5);return;end
 
112
 w=num./den+sl(5);
 
113
 if type(w)==1 then h=w;return;end   //degenerate case
 
114
 h=syslin(domaine,w);
 
115
end