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