1 #!/usr/bin/python

     2 # 13.6

     3 # Viscosity control problem of MacGregor

     4 

     5 import scipy as sp

     6 from pylab import plot, show

     7 from lqg import lqg

     8 from cl import cl

     9 

    10 A, dA = sp.array([1, -0.44]), 1

    11 B, dB = sp.array([0.51, 1.21]), 1

    12 F, dF = sp.array([1, -1]), 1

    13 C, dC = sp.array([1, -0.44]),1

    14 k = 1

    15 V, W = 1, 1

    16 dV, dW = 0, 0

    17 u_lqg = []

    18 y_lqg =[]

    19 uy_lqg = []

    20 int1 = 1

    21 

    22 for rho in sp.arange(0.001, 3, 0.1):

    23     R1,dR1,Sc,dSc = lqg(A,dA,B,dB,C,dC,k,rho,V,dV,W,dW,F,dF)

    24     Nu,dNu,Du,dDu,Ny,dNy,Dy,dDy,yvar,uvar = cl(A,dA,B,dB,C,dC,k,Sc,dSc,R1,dR1,int1)

    25     uvar = sp.atleast_1d(uvar.squeeze())

    26     yvar = sp.atleast_1d(yvar.squeeze())

    27     u_lqg = sp.concatenate((u_lqg, uvar))

    28     y_lqg = sp.concatenate((y_lqg, yvar))

    29     if sp.size(uy_lqg) == 0:

    30         uy_lqg = sp.array([rho, uvar[0], yvar[0]])

    31     else:

    32         uy_lqg = sp.vstack((uy_lqg, sp.array([rho, uvar[0], yvar[0]])))

    33 

    34 plot(u_lqg, y_lqg)

    35 show()

    36