# Hs.p ---  Gnuplot script to plot the transfer function of an LTI
#           dynamic system in the s-plane 
#           The system is two coupled oscillators

set dummy g,w	# independent variables ... real and imaginary parts of s-plane

i = sqrt(-1);				# imaginary number

s(g,w) = g + i*w;			# Laplace variable

# -----------------------------------   system properties

m1 = 1.0;  c1 = 5.00;   k1 =  40.0;
m2 = 3.0;  c2 = 0.50;   k2 = 200.0;

#    Natural                     Damped  
#   Frequency                   Frequency      Eigenvalue 
#   (cyc/sec)      Damping      (cyc/sec)     real        imag 
#   -----------------------------------------------------
#      0.47957     0.16822       0.47273     -0.50690     2.97027  
#      0.47957     0.16822       0.47273     -0.50690    -2.97027  
#      2.72756     0.13575       2.70231     -2.32644    16.97914  
#      2.72756     0.13575       2.70231     -2.32644   -16.97914  

# ----------------------------------   plot lay out

set xrange [-3.09:0]
set yrange [-20.0:20.0]
set noborder
set hidden
set view 70,45
set isosamples 50,201
set surface
set nokey
set noxtic
set noytic
set noztic

set title "Transfer Function Magnitude in the {/Times-Italic s}-plane \n \
 {/Times-Italic s = }{/Symbol s + }{/Times-Italic i}{/Symbol w}" 

# ------------------------------------- arrows

set noarrow
 set arrow from 0,0,0.001 to  0,0,10
 set arrow from 0,0,0.001 to  0,25,0.001
 set arrow from 0,0,0.001 to  0,-25,0.001
 set arrow from 0,0,0.001 to  1,0,0.001
 set arrow from 0,0,0.001 to -0.6,0,0.001
 
 set arrow from -0.50690,  2.97027,1.0 to -0.50690, 2.97027,5 
 set arrow from -0.50690, -2.97027,1.0 to -0.50690,-2.97027,5 
 set arrow from -2.32644, 16.97914,0.2 to -2.32644, 16.97914,1
 set arrow from -2.32644,-16.97914,0.2 to -2.32644,-16.97914,1
 
# ------------------------------------- labels

set nolabel
 set label "log(|{/Times-Italic H}(s)|)" at 0.2,0,8
 set label "{/Times-Italic i}{/Symbol w}"  at 0.1, 20,0.001
 set label "{/Symbol -}{/Times-Italic i}{/Symbol w}" at 0.1,-22,0.001
 set label "{/Symbol s}"  at  1.100,0,0.001
 set label "{/Symbol -s}" at -0.800,0,0.001
 set label "{/Times-Italic p}_1"   at -0.50690,  2.97027, 7
 set label "{/Times-Italic p}_1^*" at -0.50690, -2.97027, 7
 set label "{/Times-Italic p}_2"   at -2.32644, 16.97914, 1.5
 set label "{/Times-Italic p}_2^*" at -2.32644,-16.97914, 1.5

 set label "x"   at -0.50690,  2.97027, 0.001
 set label "x"   at -0.50690, -2.97027, 0.001
 set label "x"   at -2.32644, 16.97914, 0.001
 set label "x"   at -2.32644,-16.97914, 0.001

# ---------------------------------------- functions to plot

D(g,w) = ( m1*s(g,w)*s(g,w) + (c1+c2)*s(g,w) + (k1+k2) ) * \
	 ( m2*s(g,w)*s(g,w) + c2*s(g,w) + k2 );

# Hs(g,w) = s(g,w)*s(g,w)*( c2*s(g,w) + k2 ) / ( D(g,w) - c2*s(g,w) - k2 ) ;

Hs(g,w)   = ( c2*s(g,w) + k2 ) / ( D(g,w) - (c2*s(g,w) + k2)**2.0 ) ;

# ----------------------------------------- make the plots

set nolog;  splot real(Hs(g,w)); pause -1;	# surface plot of real part

set nolog;  splot imag(Hs(g,w)); pause -1;	# surface plot of imag part

set nolog;  splot arg(Hs(g,w));  pause -1;	# surface plot of phase

set log z;  splot abs(Hs(g,w)) ;		# surface plot of magnitude

# --------------------------------------- save the plots as post script

  set out 'Hs.ps'
  set terminal postscript landscape enhanced "Helvetica" 14 
  replot
  set size 1,1
  set terminal x11
 
  !gs Hs.ps