# set terminal pngcairo  transparent enhanced font "arial,8" fontscale 1.0 size 610, 480 
# set output 'multiplt.1.png'
set bar 1.000000 front
set style circle radius graph 0.02, first 0.00000, 0.00000 
set style ellipse size graph 0.05, 0.03, first 0.00000 angle 0 units xy
set dummy wnt, y
set grid nopolar
set grid xtics nomxtics ytics nomytics noztics nomztics \
 nox2tics nomx2tics noy2tics nomy2tics nocbtics nomcbtics
set grid layerdefault   lt 0 linewidth 0.500,  lt 0 linewidth 0.500
unset key
set label 1 "" at 0.140000, 17.0000, 0.00000 left norotate back nopoint
set style textbox transparent margins  1.0,  1.0 border
unset logscale
set samples 250, 250
set size ratio 0 0.5,0.5
set origin 0.5,0
set mxtics 10.000000
set mytics 5.000000
set xtics  norangelimit 0.00000,5,20.0000
set ytics  norangelimit -1.00000,0.5,1.00000
unset paxis 1 tics
unset paxis 2 tics
unset paxis 3 tics
unset paxis 4 tics
unset paxis 5 tics
unset paxis 6 tics
unset paxis 7 tics
set title "Second Order System - Unit Impulse Response" 
set xlabel "Normalized Time (wnt)" 
set xrange [ 0.00000 : 20.0000 ] noreverse nowriteback
set ylabel "Amplitude y(wnt)" 
set ylabel  offset character 1, 0, 0 font "" textcolor lt -1 rotate by -270
set yrange [ -1.00000 : 1.00000 ] noreverse nowriteback
set paxis 1 range [ * : * ] noreverse nowriteback
set paxis 2 range [ * : * ] noreverse nowriteback
set paxis 3 range [ * : * ] noreverse nowriteback
set paxis 4 range [ * : * ] noreverse nowriteback
set paxis 5 range [ * : * ] noreverse nowriteback
set paxis 6 range [ * : * ] noreverse nowriteback
set paxis 7 range [ * : * ] noreverse nowriteback
set colorbox vertical origin screen 0.9, 0.2, 0 size screen 0.05, 0.6, 0 front  noinvert bdefault
mag(w) = -10*log10( (1-w**2)**2 + 4*(zeta*w)**2)
tmp(w) = (-180/pi)*atan( 2*zeta*w/(1-w**2) )
tmp1(w)= w<1?tmp(w):(tmp(w)-180)
phi(w)=zeta==1?(-2*(180/pi)*atan(w)):tmp1(w)
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)=exp(-zeta*wnt) * sin(wdwn(zeta)*wnt) / wdwn(zeta)
zeta = 2
GPFUN_mag = "mag(w) = -10*log10( (1-w**2)**2 + 4*(zeta*w)**2)"
GPFUN_tmp = "tmp(w) = (-180/pi)*atan( 2*zeta*w/(1-w**2) )"
GPFUN_tmp1 = "tmp1(w)= w<1?tmp(w):(tmp(w)-180)"
GPFUN_phi = "phi(w)=zeta==1?(-2*(180/pi)*atan(w)):tmp1(w)"
GPFUN_wdwn = "wdwn(zeta)=sqrt(1-zeta**2)"
GPFUN_shift = "shift(zeta) = atan(wdwn(zeta)/zeta)"
GPFUN_alpha = "alpha(zeta)=zeta>1?sqrt(zeta**2-1.0):0"
GPFUN_tau1 = "tau1(zeta)=1/(zeta-alpha(zeta))"
GPFUN_tau2 = "tau2(zeta)=1/(zeta+alpha(zeta))"
GPFUN_c1 = "c1(zeta)=(zeta + alpha(zeta))/(2*alpha(zeta))"
GPFUN_c2 = "c2(zeta)=c1(zeta)-1"
GPFUN_y1 = "y1(wnt)=zeta==1?1 - exp(-wnt)*(wnt + 1):0"
GPFUN_y2 = "y2(wnt)=zeta<1?(1 - (exp(-zeta*wnt)/wdwn(zeta))*sin(wdwn(zeta)*wnt + shift(zeta))):y1(wnt)"
GPFUN_y = "y(wnt)=exp(-zeta*wnt) * sin(wdwn(zeta)*wnt) / wdwn(zeta)"
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)