--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/specs/rlocus_ex1.py	Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,62 @@
+#!/usr/bin/python
+# 7.1
+import os, sys
+sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"]
+
+import pylab as pl
+from scipy import signal
+
+def rlocus(H, k_vec=None):
+    """Plots the Root Locus given 
+    an lti system and 
+    a vector of gain values"""
+    if k_vec is None:
+        k_vec = linspace(0, 100, 10000)
+
+    num = pl.poly1d(H.num)
+    den = pl.poly1d(H.den)
+    
+    roots = []
+    for k in k_vec:
+        p = den + k*num
+        roots.append(pl.roots(p))
+    roots = pl.array(roots)
+
+    # Tried to do it with array slicing. Looks BAD
+    # n = H.num
+    # d = H.den
+    # l = d.shape[0]-n.shape[0]
+    # n_vec = n[pl.newaxis,:]*k_vec[:,pl.newaxis]
+    # A = pl.zeros((k_vec.shape[0], d.shape[0]))
+    # A = A+d
+    # A[:,l:] = n_vec 
+    # r = map(pl.roots, A)
+    # r = pl.array(r)
+
+    x = pl.linspace(0, 2*pl.pi, 100)
+    pl.plot(pl.cos(x), pl.sin(x), 'b.', markersize=2)
+
+    for i, k_root in enumerate(roots[0]):
+        pl.plot(pl.real(roots[:,i]), pl.imag(roots[:,i]))
+  
+    xmin, xmax = pl.xlim()
+    x = pl.linspace(int(xmin), int(xmax)) 
+    ymin, ymax = pl.ylim()
+    y = pl.linspace(int(ymin), int(ymax))
+
+    pl.plot(x, pl.zeros_like(x), 'b.', pl.zeros_like(y), y, 'b.', markersize=2)
+
+    pl.xlim(xmin, xmax)
+    pl.ylim(ymin, ymax)
+
+    pl.xlabel('Real part')
+    pl.ylabel('Imaginary part')
+
+    pl.show()
+
+n = [1]
+d = [1, -1, 0]
+k_vec = pl.linspace(0, 20, 10000)
+H = signal.lti(n, d)
+
+rlocus(H, k_vec)