תורת הבקרה/עקום בודה/תכנות

from numpy import *
import pylab as p
 
def G(s):
	return 2/((s+1)*(s+2))		# Transfer function
G_latex='$\\frac{2}{(s+1)(s+2)}$'	# To be used for plotting
 
f_i = 0.01				# There's singularuty in f=0
f_f = 100				# because log(0) is -inf.
f_steps = 100
freq=logspace(log10(f_i),log10(f_f),f_steps)
 
log10 = log(10)
i=1j
 
magnitude = array([20*log(abs(G(n*i)))/log10 for n in freq])	# Magnitude
phase = array([angle(G(n*i))*180/pi for n in freq])		# Phase
 
freq = freq.real
magnitude = magnitude.real
phase = phase.real
 
# Magnitude sub-plot
p.subplot(2,1,1)
p.semilogx(freq,magnitude)
p.grid(True)
p.ylabel('Magnitude [dB]')
p.title('Bode Diagram of ' + G_latex)
 
# Phase sub-plot
p.subplot(2,1,2)
p.semilogx(freq,phase)
p.grid(True)
p.xlabel('Frequency [rad/s]')
p.ylabel('Phase [deg]')
 
p.show()