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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
|
# GNUPLOT v3.6 beta multiplot script file
#
# Second Order System Characteristics
#
# D**2 + 2*zeta*wn*D + (wn**2)y = (wn**2)*x
#
# x input variable
# y output variable
# w frequency ratio (w/wn)
# wn natural frequency
# wd damped natural frequency
# zeta damping ratio
# mag(w) amplitude response
# phi(w) phase response
# wdwn damped natural frequency ratio
# wnt normalized time
#
# Plots:
# Frequency domain magnitude response
# phase response
#
# Time domain unit step response
# unit impulse response
#
#
# Created by: W. D. Kirby email: wdkirby@ix.netcom.com
# Date: 1/18/96
# Released to the public domain with no warranty of any kind
#
reset
set style function lines
set size 1.0, 1.0
set origin 0.0, 0.0
set multiplot
set size 0.5,0.5
set origin 0.0,0.5
set grid
unset key
set angles radians
set samples 250
# Plot Magnitude Response
set title "Second Order System Transfer Function - Magnitude"
mag(w) = -10*log10( (1-w**2)**2 + 4*(zeta*w)**2)
set dummy w
set logscale x
set xlabel "Frequency (w/wn)"
set ylabel "Magnitude (dB)" offset 1,0
set label 1 "Damping =.1,.2,.3,.4,.5,.707,1.0,2.0" at .14,17
set xrange [.1:10]
set yrange [-40:20]
plot \
zeta=.1,mag(w), \
zeta=.2,mag(w), \
zeta=.3,mag(w), \
zeta=.4,mag(w), \
zeta=.5,mag(w), \
zeta=.707,mag(w), \
zeta=1.0,mag(w), \
zeta=2.0,mag(w),-6
# Plot Phase Response
set size 0.5,0.5
set origin 0.0,0.0
set title "Second Order System Transfer Function - Phase"
set label 1 ""
set ylabel "Phase (deg)" offset 1,0
set ytics -180, 30, 0
set yrange [-180:0]
tmp(w) = (-180/pi)*atan( 2*zeta*w/(1-w**2) )
# Fix for atan function wrap problem
tmp1(w)= w<1?tmp(w):(tmp(w)-180)
phi(w)=zeta==1?(-2*(180/pi)*atan(w)):tmp1(w)
plot \
zeta=.1,phi(w), \
zeta=.2,phi(w), \
zeta=.3,phi(w), \
zeta=.4,phi(w), \
zeta=.5,phi(w), \
zeta=.707,phi(w), \
zeta=1,phi(w), \
zeta=2.0,phi(w), \
-90
# Plot Step Response
set size 0.5,0.5
set origin 0.5,0.5
set dummy wnt
unset logscale x
set title "Second Order System - Unit Step Response"
set ylabel "Amplitude y(wnt)" offset 1,0
set xlabel "Normalized Time (wnt)"
set xrange [0:20]
set xtics 0,5,20
set yrange [0:2.0]
set ytics 0, .5, 2.0
set mytics 5
set mxtics 10
wdwn(zeta)=sqrt(1-zeta**2)
shift(zeta) = atan(wdwn(zeta)/zeta)
alpha(zeta)=zeta>1?sqrt(zeta**2-1.0):0
tau1(zeta)=1/(zeta-alpha(zeta))
tau2(zeta)=1/(zeta+alpha(zeta))
c1(zeta)=(zeta + alpha(zeta))/(2*alpha(zeta))
c2(zeta)=c1(zeta)-1
y1(wnt)=zeta==1?1 - exp(-wnt)*(wnt + 1):0
y2(wnt)=zeta<1?(1 - (exp(-zeta*wnt)/wdwn(zeta))*sin(wdwn(zeta)*wnt + shift(zeta))):y1(wnt)
y(wnt)=zeta>1?1-c1(zeta)*exp(-wnt/tau1(zeta))+c2(zeta)*exp(-wnt/tau2(zeta)):y2(wnt)
plot \
zeta=.1,y(wnt), \
zeta=.2,y(wnt), \
zeta=.3,y(wnt), \
zeta=.4,y(wnt), \
zeta=.5,y(wnt), \
zeta=.707,y(wnt), \
zeta=1,y(wnt), \
zeta=2,y(wnt)
#
# Plot Impulse Response
set origin .5,0.
set title "Second Order System - Unit Impulse Response"
y(wnt)=exp(-zeta*wnt) * sin(wdwn(zeta)*wnt) / wdwn(zeta)
set yrange [-1. :1.]
set ytics -1,.5,1.
plot \
zeta=.1,y(wnt), \
zeta=.2,y(wnt), \
zeta=.3,y(wnt), \
zeta=.4,y(wnt), \
zeta=.5,y(wnt), \
zeta=.707,y(wnt), \
zeta=1,y(wnt), \
zeta=2,y(wnt)
unset multiplot
pause -1 "Hit return to continue"
#
# Clean up: reset parameter defaults
#
reset
|
