1 #!/usr/bin/python

     2 # 7.1

     3 import os, sys

     4 sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"]

     5 

     6 import pylab as pl

     7 from scipy import signal

     8 

     9 def rlocus(H, k_vec=None):

    10     """Plots the Root Locus given 

    11     an lti system and 

    12     a vector of gain values"""

    13     if k_vec is None:

    14         k_vec = linspace(0, 100, 10000)

    15 

    16     num = pl.poly1d(H.num)

    17     den = pl.poly1d(H.den)

    18     

    19     roots = []

    20     for k in k_vec:

    21         p = den + k*num

    22         roots.append(pl.roots(p))

    23     roots = pl.array(roots)

    24 

    25     # Tried to do it with array slicing. Looks BAD

    26     # n = H.num

    27     # d = H.den

    28     # l = d.shape[0]-n.shape[0]

    29     # n_vec = n[pl.newaxis,:]*k_vec[:,pl.newaxis]

    30     # A = pl.zeros((k_vec.shape[0], d.shape[0]))

    31     # A = A+d

    32     # A[:,l:] = n_vec 

    33     # r = map(pl.roots, A)

    34     # r = pl.array(r)

    35 

    36     x = pl.linspace(0, 2*pl.pi, 100)

    37     pl.plot(pl.cos(x), pl.sin(x), 'b.', markersize=2)

    38 

    39     for i, k_root in enumerate(roots[0]):

    40         pl.plot(pl.real(roots[:,i]), pl.imag(roots[:,i]))

    41   

    42     xmin, xmax = pl.xlim()

    43     x = pl.linspace(int(xmin), int(xmax)) 

    44     ymin, ymax = pl.ylim()

    45     y = pl.linspace(int(ymin), int(ymax))

    46 

    47     pl.plot(x, pl.zeros_like(x), 'b.', pl.zeros_like(y), y, 'b.', markersize=2)

    48 

    49     pl.xlim(xmin, xmax)

    50     pl.ylim(ymin, ymax)

    51 

    52     pl.xlabel('Real part')

    53     pl.ylabel('Imaginary part')

    54 

    55     pl.show()

    56 

    57 n = [1]

    58 d = [1, -1, 0]

    59 k_vec = pl.linspace(0, 20, 10000)

    60 H = signal.lti(n, d)

    61 

    62 rlocus(H, k_vec)