damp function added in CACSD module
[scilab.git] / scilab / modules / cacsd / macros / damp.sci
1 function [wn,z,p] = damp(R,dt1)
2 //Natural frequency and damping factor for continuous systems.
3 //   [Wn,Z,P] = damp(R) returns vectors Wn and Z containing the
4 //   natural frequencies and damping factors of R.
5 //   The variable R  must be a real or complex array of roots, a
6 //   polynomial array or  a linear dynamical system
7
8   if argn(2)<1 then
9     error(msprintf(_("%s: Wrong number of input arguments: %d or %d expected.\n"),"damp",1,2))
10   end
11   //handling optionnal argument dt1
12   if argn(2)==1 then
13     dt1=[];
14   else
15     if type(dt1)==1 then
16       if size(dt1,'*')<>1|dt1<0 then
17         error(msprintf(_("%s: Wrong type for input argument #%d: Real non negative scalar expected.\n"),..
18                        "damp",2))
19       end
20     elseif type(dt1)==10 then
21       if dt1=="c" then
22         dt1=0
23       elseif dt1=="d" then
24         dt1=1
25       else
26         error(msprintf(_("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),..
27                        "damp",2,"""c"", ""d"""))
28       end
29     else
30       error(msprintf(_("%s: Wrong type for input argument #%d: Scalar or string expected.\n"),..
31                      "damp",2))
32     end
33   end
34
35   toBeOrdered=%t;dt=[];
36   select typeof(R)
37   case "polynomial" then //polynomial array
38     p=[];
39     for k=1:size(R,'*')
40       p=[p;roots(R(k),"e")];
41     end
42   case "rational" then
43     dt=R.dt
44     if dt=="c" then
45       dt=0
46     elseif dt=="d" then
47       dt=1
48     end
49     p=roots(lcm(R.den))
50   case "state-space" then
51     dt=R.dt
52     if dt=="c" then
53       dt=0
54     elseif dt=="d" then
55       dt=1
56     end
57     p=spec(R.A)
58   case "constant" then
59     p=R;
60     toBeOrdered=%f
61   else
62     error(msprintf(_("%s: Wrong type for input argument #%d: Array of floats or Polynomial expected.\n"),..
63                    "damp",1))
64   end
65   if dt==[] then
66     //R does not furnish time domain
67     if dt1==[] then
68       //no user time domain specified, continuuous time assumed
69       dt=0
70     else
71       //user time domain specified
72       dt=dt1
73     end
74   elseif dt1<>[] then
75     warning(msprintf(_("%s: Input argument #%d ignored.\n"),"damp",2))
76   end
77   // Initialize
78   wn=zeros(p);
79   z=-ones(p);
80   im=ieee();ieee(2);//to allow inf and nan's
81   if dt>0 then // Discrete  time case
82     ind=find(abs(p-1)>10*%eps)
83     s=p(ind);
84     s=log(s)/dt;
85   else //continuous time case
86     ind=find(p<>0)
87     s=p(ind);
88   end
89   wn(ind)=abs(s)
90   z(ind)=-real(s)./abs(s)
91   ieee(im)
92   if toBeOrdered then
93     [wn,k]=gsort(wn,'g','i');
94     z=z(k);
95     p=p(k)
96   end
97 endfunction